Real-time hydrogen molecular dynamics satisfying the nuclear spin statistics of a quantum rotor

Apparent presence of the nuclear-spin species of a hydrogen molecule, para-hydrogen and ortho-hydrogen, associated with the quantum rotation is a manifestation of the nuclear quantum nature of hydrogen, governing not only molecular structures but also physical and chemical properties of hydrogen molecules. It has been a great challenge to observe and calculate real-time dynamics of such molecularized fermions. Here, we developed the non-empirical quantum molecular dynamics method that enables real-time molecular dynamics simulations of hydrogen molecules satisfying the nuclear spin statistics of the quantum rotor. While reproducing the species-dependent quantum rotational energy, population ratio, specific heat, and H-H bond length and frequency, we found that their translational, orientational and vibrational dynamics becomes accelerated with the higher rotational excitation, concluding that the nuclear quantum rotation stemmed from the nuclear spin statistics can induce various kinds of dynamics and reactions intrinsic to each hydrogen species.


S1. Computational methods: the nuclear and electron wave packet molecular dynamics method with quantum rotation
We developed the non-empirical quantum molecular dynamics method satisfying the nuclear spin statistics of a quantum rotor called the nuclear and electron wave packet molecular dynamics method with quantum rotation (the NEWPMD-QR method) by extending the original NEWPMD method free from the nuclear spin statistics and will be called the Gaussian NEWPMD method (the G-NEWPMD method) in this Article. As will be explained below, the NEWPMD-QR method provides equations of motion (EOMs) to time-evolve a nuclear-electron wave function of a rotationally-excited hydrogen molecule.

S1A. Time-dependent nuclear-electron wave function satisfying the nuclear spin statistics of a quantum rotor
The NEWPMD-QR method describes two nuclei, A and B, by nuclear wave packets (NWPs) via the time-dependent Hartree approach, while it expresses two electrons, L and S, by two Gaussian electron wave packets (EWPs) through the perfect-pairing (PP) valence bond (VB) theory that appropriately treats the Pauli exclusion energy. The time-dependent nuclearelectron wave function for a hydrogen molecule composed of the two nuclei A and B and two electrons L and S with the rotational quantum number J is introduced as Ψ L,S AB,J (t) = Φ AB,J (Q 1 , Q 2 , t)ψ L,S (q 1 , q 2 , t), where Q 1 , Q 2 , q 1 and q 2 denote A-, B-, L-and S-coordinates, respectively. The normalized time-dependent PP VB electronic wave function for the two electrons is expressed as ψ L,S (q 1 , q 2 , t) = 1 Here, α and β are electron spin functions with the two electrons forming the singlet. The normalized Gaussian EWP φ k (q, t) is specified by the position of the EWP center R k e (t) = (X k e (t), Y k e (t), Z k e (t)) and the EWP width ρ k (t): with the normalization factor N k = (2πρ 2 k ) 3/4 for k =L and S. For simplicity,h = 1 and the electron charge and mass are set as unity for scaling. S LS (t) is the overlap integral between φ L and φ S : In order to obtain the time-dependent nuclear wave function for the two hydrogen nuclei A and B with the rotational quantum number J, Φ AB,J (Q 1 , Q 2 , t), we start with a normalized Gaussian NWP specified by the NWP center position R k (t) and width Ω k (t) accompanied with their conjugate momenta P k (t) and Π k (t), respectively, as for K =A or B. The time-dependent nuclear-electron wave function for a hydrogen molecule introduced in the G-NEWPMD method can be written by simply adopting the Gaussian NWPs and Gaussian EWPs as Figure 1(a) shows the nuclear distribution of the stable hydrogen molecule, calculated by the G-NEWPMD method. Its nuclear distribution has the dumbbell-shape nuclear delocalization composed of the two distinct Gaussian NWPs. Therefore, we will call this hydrogen molecule calculated by the G-NEWPMD method Gaussian throughout this Article.
The Hartree product of the two normalized Gaussian NWPs, can be transformed with the center-of-mass (COM) coordinates and the relative coordinates Here, we defined the COM position as R COM (t) = (R A (t)+R B (t))/2 and the relative position defined by the distance between the two NWP center positions as and their conjugate momenta as P COM (t) and P rel (t), respectively. The three dimensional absolute coordinates Q 1 and Q 2 were transformed to the COM and relative coordinates Q COM and Q rel to effectively express the nuclear quantum rotation.
We calculate hydrogen molecules in the ground (J = 0), first-excited (J = 1), secondexcited (J = 2) and third-excited (J = 3) quantum rotational states that will be called para, ortho, para-2, and ortho-2 throughout this Article, respectively. Although para on the ground rotational state should have a completely spherical nuclear wave function, eq.(12) itself is a simple deformation of the Hartree product of the paired Gaussian NWPs, and cannot describe such spherical nuclear wave function. Therefore, to obtain the timedependent nuclear wave function with the nuclear quantum rotation, we modify eq.(12) as where Q rel = |Q rel |, R rel (t) = |R A (t) − R B (t)| and P rel (t) = |P rel (t)|. Now eq.(13) can express the completely spherical shell-shape nuclear wave function for para owing to the spherical rotational wave function, exp[−(Q rel − R rel (t)) 2 /8Ω 2 A (t)]. Since the total nuclear wave function must be symmetric and antisymmetric with the anti-parallel and parallel nuclear spins, respectively, we further multiply the spherical harmonics, Y Jm (χ, ω), to eq.(13) as the angular wave function for para (J = 0), ortho (J = 1), para-2 (J = 2), and ortho-2 It should be noted that χ and ω in Y Jm (χ, ω) do not specify the actual H-H molecular axis but denote the angular coordinates of Q rel from the x-axis and from the y-axis on the yz-plane, respectively, that is, Q rel,x = Q rel cos χ, Q rel,y = Q rel sin χ cos ω, Q rel,z = Q rel sin χ sin ω.

S1B. Kinetic energy and electrostatic interaction energy
The total molecular energy, E tot,J , the sum of the kinetic energy of the hydrogen nuclei and electrons as well as the three electrostatic interaction energy of electron-electron, nucleusnucleus and nucleus-electron, is indispensable for deriving the EOMs. Because the electronic wave function for the two electrons, eq.(2), is common in the NEWPMD-QR and G-NEWPMD methods, the kinetic energy of the two electrons with the one-electron kinetic energy and the electron-electron interaction energy with the two-electron integral where α = 1/4ρ 2 i , β = 1/4ρ 2 j , γ = 1/4ρ 2 k , and δ = 1/4ρ 2 l , and are the same both in the NEWPMD-QR and G-NEWPMD methods. Note that Ψ L,S AB,J (t)| · · · |Ψ L,S AB,J (t) means a quantum expectation with the nuclear-electron wave function Ψ L,S AB,J (t).
On the contrary, the kinetic energy of the hydrogen nuclei and the electrostatic interaction energy of nucleus-nucleus and nucleus-electron significantly differ between the NEWPMD-QR and G-NEWPMD methods. The kinetic energy of the hydrogen nuclei depending on the rotational quantum number J is derived in the NEWPMD-QR method as . M nuc is relative mass of the hydrogen nucleus to the electron.
Interestingly, the electrostatic interaction energy of nucleus-nucleus calculated by the NEWPMD-QR method was found to be independent of the rotational quantum number J, ] This is because the radial nucleus-nucleus interaction depends only on the internuclear distance Q rel , not being influenced by the angular coordinates χ and ω.
The electrostatic interaction energy of nucleus-electron is also independent of the rotational quantum number J due to their radial Coulomb interaction free from the angular coordinates, being derived as The explicit analytical form of the nucleus-electron integral V LS,AB in the NEWPMD-QR method was calculated as with and The remained integral by Q rel in eq.(41) can be performed by expanding the last line as Then, the nucleus-electron integral finally becomes This expression does not change even if the integrand in eq.(35), 1/|q 1 − Q COM − Q rel /2|, is replaced by 1/|q 1 − Q COM + Q rel /2| corresponding to the other nucleus.
For the shell-type species in the different nuclear rotational states, we freely optimized

S1C. Time-dependent quantum variational principle for the nuclearelectron wave function with the nuclear quantum rotation
The EOMs to time-evolve R COM (t), R rel (t), Ω A (t), P COM (t), P rel (t), and Π A (t) that specify dynamics of both the NWPs and EWPs can be derived through the time-dependent quantum variational principle. The time-dependent quantum variational principle minimizes the action integral for the time-dependent nuclear-electron wave function Ψ L,S AB,J (t) in eq.(1) defined with the quantum Lagrangian L. Since the Hamiltonian operatorĤ is composed of the kinetic energy operators of the electrons (eq.(18)) and hydrogen nuclei (eq.(26)) as well as the three electrostatic interaction energy operators of electron-electron (eq.(21)), nucleusnucleus (eq.(29)) and nucleus-electron (eq.(32)), Ψ L,S AB,J (t)|Ĥ|Ψ L,S AB,J (t) = E tot,J from the definition, eq.(17). The time-dependent variational principle, δΓ/δR COM (t) = 0, etc., yields the EOMs, where , Ω A (t), P COM (t), P rel (t), Π A (t)}. The explicit forms of these EOMs are written as and It should be mentioned that, since the difference in Ψ L,S AB,J (t) among the shell-type species is only the spherical harmonics Y Jm (χ, ω), the terms on the left-hand side of eqs.(48)-(53) are the same regardless of the shell-type species (J-independent), and that it is the total molecular energy E tot,J that distinguishes the shell-type species. The EOMs (48)-(53) can be drastically simplified if we adopt the reasonable approximations R rel (t) Ω A (t) and P rel (t) Π A (t), which actually always holds during all molecular dynamics simulations reported in this Article, resulting iṅ andΩ It is remarkable that the EOMs (54)-(59) are universal and can be systematically applicable not only to the current four shell-type species (para, ortho, para-2, and ortho-2) but also to another shell-type species possessing higher-energy nuclear quantum rotation. The NEWPMD-QR method is the first computational method to calculate real-time molecular dynamics on a nuclear excited state not on an electronically excited state.

S1D. Interaction potential energy between a rotationally-excited hydrogen molecule and carbon nano tube
In order to simulate real-time collision dynamics of the shell-type species with a single- The interaction potential energy between the shell-type species and CNT(15,0), E CNT,J (r, θ), was calculated by further averaging the total energy function E CNT (r, χ CNT ) with respect to the rotational degrees of freedom, Here, the J-dependent normalization factors are

S1E. Equations of motion
The anisotropic interaction with the CNT(15,0) surface breaks the spatial translational symmetry on the cross section of CNT(15,0), enabling the introduction of a laboratory coordinate system. Therefore, by transforming the EOMs (54)-(59), we derived the final EOMs aṡ andΩ Note that we set the origin as the center of CNT(15,0) in the above EOMs by replacing r with D CNT /2 − r as seen in E CNT,J (D CNT /2 − r, θ). of Gaussian and the shell-type species (para, ortho, para-2, and ortho-2) obtained from the peak frequencies of their H-H vibrational power spectra displayed in Fig.5(a-d). The numbers in the parentheses denote the deviation between the CNT-0 and experimental values.  Fig.2(a), respectively. of the EWP width on R rel is identical regardless of the shell-type species.