Effect of nuclear spin symmetry in cold and ultracold reactions: D + para/ortho-H$_2$

We report results for reaction and vibrational quenching of the collision D with para-H$_2$($v,j=0$) and ortho-H$_2$($v,j=1$) at cold and ultracold temperatures. We investigate the effect of nuclear spin symmetry for barrier dominated processes ($0\le v\le 4$) and for one barrierless case ($v=5$). We find resonant structures for energies in the range corresponding to 0.01--10 K, which depend on the nuclear spin of H$_2$, arising from contributions of specific partial waves. We discuss the implications on the results in this benchmark system for ultracold chemistry.


Introduction
Research with cold and ultracold molecules has witnessed an explosive growth since the first predictions [1,2] and observations [3,4] of ultracold molecules in the late 1990s. Already, several books [5,6,7] have been published on cold and ultracold molecules, and numerous reviews have been written on more specialized topics [8,9,10,11,12], illustrating the continuously broadening scope of this field, as exemplified by this special issue.
In addition to new proposals to produce ultracold molecules (as in Ref. [13]), this rapidly expanding field has motivated theoretical studies of the chemical structure and properties of ultracold molecules, e.g., for the interaction of two KRb molecules [14]. The extraordinary level of control of all degrees of freedom, from internal to translational degrees, reached in ultracold systems has allowed the investigation of the effect of the nuclear spin on atom-diatom scattering, as demonstrated in Ref. [15] with KRb+K.
In this article, we explore such questions using a benchmark system exhibiting both quenching and reaction, namely H 2 +D, which can react to produce HD+H. Because of its small mass, this system offers the opportunity to study the contributions of low partial waves even at cold temperatures, as opposed to heavier systems where ultracold temperatures are required. This particular system is fundamental in quantum chemistry, and is being actively investigated in experiments capable of reaching the energy range explored in this work [16]. The detection of the resonance features we predict would be a valuable test of the ab initio potential energy surface calculation, and their dependence on the nuclear symmetry provide an even stricter comparison between experiment and theory, the need of which is discussed in Ref. [17]. In addition, H 2 +D is also relevant to astrophysics, especially in the astrochemistry of cold interstellar clouds [18], and even in the early universe [19]. Although many studies involving the scattering of H 2 with various atoms have been performed in the ultracold regime [20,21,22,23,24], the effect of the nuclear symmetry has not been studied for this reactive system at low temperatures. However, nuclear symmetry has been considered in studies of paraand ortho-H 2 , e.g., the rovibrational relaxation and energy transfer in ultracold H 2 +H 2 collisions [25,26,27], and in studies of vibrational quenching of O 2 when colliding with He [28], or spin-changing scattering between two ground state O 2 molecules [29].
In Sec. 2 we give a brief description of the theoretical and numerical tools used, as well as the properties of this benchmark system. We present and analyze the results in Sec. 3, and conclude in Sec. 4.

Theoretical and computational details
We are interested in the effect of the nuclear spin symmetry of homonuclear molecules in reactive scattering processes at low temperatures. To this end, we consider the benchmark reaction H 2 +D which involves the simplest and most fundamental diatomic molecule, H 2 . The nuclear part of the wave function of H 2 must be anti-symmetric with respect to the permutation of the two protons, and since the nuclear spin wave function of two spin i = 1/2 nuclei with total nuclear spin I acquires a factor (−1) 2i−I under exchange, the rotational states of molecular hydrogen are restricted: para-Hydrogen corresponds to the singlet spin state with (I = 0, M I = 0) which allows only even-j rotational states, while the nuclear triplet spin state of ortho-Hydrogen with (I = 1, M I = 0, ±1) allows only odd-j rotational states.

Numerical approach
The expression for the state-to-state cross sections, integrated over all scattering directions, averaged over the initial rotational states of the reactant dimer, and summed over the final rotational states of the product, reads where the generic notation n = (avj) stands for the arrangement label "a" and quantum numbers (vj) of the molecular H 2 states, and k n = √ 2µ a E kin is the initial momentum (h = 1, atomic units are used); µ a the reduced mass of the binary system atom-diatom in the initial arrangement (a), the initial kinetic energy is E kin = E − ε n , and ε n are the rovibrational energies of the dimers. E is the (total) collision energy, J is the total angular momentum, and ℓ is angular momentum for the relative motion. The primes shown: open channels with energies smaller than the absolute scattering energy E for quenching and reactions, closed channels with energies larger than E affecting the results, and the neglected channels above a certain truncation energy E max . The inset shows a zoom of the levels near v = 1, with the scattering energy ε relative to the entrance channel (v = 1, j = 0).
indicate the corresponding quantities in the exit channel, with n ′ = (a ′ v ′ j ′ ) and angular momentum ℓ ′ . The T-matrix T J n ′ ℓ ′ nℓ = δ n ′ n δ ℓ ′ ℓ −S J n ′ ℓ ′ nℓ is given in terms of the S-matrix. We obtain an energy dependent rate constant by multiplying the cross sections by the relative velocity v rel. in the initial arrangement a (D + para/ortho-H 2 ) of the incoming scattering channel, namely where v rel. = 2E kin /µ a . The total rate constant is obtained by summing over the appropriate final states n ′ . If the arrangement remains the same (a ′ = a), the scattering process corresponds to quenching of H 2 by D to lower states, and if a ′ = a, a reaction occurred leading to HD + H. The corresponding rate constants are simply We note that if n ′ = n, we obtain an elastic collision: this is discussed in more detail in § 3.4.
To solve for the S-matrix and thus Eq.(1), we modified the ABC reactive scattering code developed by Manolopoulos and coworkers [30], which is based on Delves [31] hyperspherical coordinates, with a separate set of Delves hyperspherical coordinates for each arrangement sharing the same hyperradius coordinate. The full wave function is expanded in eigenbases corresponding to internal motion coordinates for each arrangement (after the combined multiple-arrangement basis is orthogonalized), and the resulting hyperradial coupled-channel equations are solved using the log derivative method [32]. The propagation starts at ρ = ρ min (where the potential is highly repulsive), and ends at a sufficiently large ρ = ρ max , where the S-matrix is extracted by imposing asymptotic boundary conditions. We note that the ABC code is well suited for system involving barriers; it was tested at high energies for a number of benchmark systems such as H + H 2 , F + H 2 , Cl + H 2 , and their isotopic counterparts [33,34,35,36,37], and modified versions to study certain benchmark chemical reactions in the ultracold regime [21,38,22,39,40].
To obtain fully converged results is expensive computationally, and a sufficient number of terms in the sum in Eq.(1) are required for higher energies. In Sec. 3 our results will show clearly that a maximum value of J max = 4 is required for energies in the kelvin regime, i.e., for E ≤ 10 K. For much higher energies, which are outside the scope of this work, the value of J max increases steadily. We remark that although the blocks for each J are separate (due to conservation of total angular momentum), the coupled problem for a given J becomes rapidly prohibitive as J increases. As described in [41], we modified the ABC code's structure and adapted it to the cold and ultracold regime; most important in our implementation is the ability to follow the progress of convergence in great detail with a save-and-restart feature, which allows us to go back and extend the radial propagation of a given run if needed. We typically monitor the convergence simultaneously for a large number of collision energies and for many initial states in a single run. Typically, we propagate to ρ max ≥ 40 a.u. to obtain converged results using small integration step size (∆ρ ≤ 0.002 a.u.) [41]. Finally, in addition to truncating the sum over J (and partial waves), we also restrict the number of channels by fixing a maximum energy E max . As illustrated in Fig. 1 for the case of the (v = 1, j = 0) initial channel of H 2 , closed channels above E max are neglected. This figure also depicts the difference in the number of molecular states between para-H 2 and HD; The density of states in HD is larger not only because it is more massive, but also because the restriction on rotational states j is lifted.

Interactions and Potential Energy Surface (PES)
The benchmark system we consider here, H 2 +D, has already been studied at ultracold [41] and higher temperatures [42]. Accurate ab initio potential energy surfaces (PES) exist [43,42] for this system. As in our previous studies [41,44], we adopted the electronic ground state PES of Ref. [43].
This system possesses an activation barrier, which is illustrated by a contour plot of the vibrational levels in Fig. 2(a,b); basically, the barrier is present for any geometry for v = 0, requiring tunneling for reactions to occur, while v = 1 is partially affected by the barrier (i.e. the barrier is present for a finite range of angles only). For higher v, the barrier does not block reactions (although it still affects the scattering rates). This figure also shows the reaction coordinates ξ. Following the line ξ, we plot the barrier for different angles θ and the lowest energy levels of H 2 and HD in Fig. 2(c). We can clearly see that v = 0 requires tunneling to produce HD, while v = 1 is higher than the barrier at angles starting at 112 o .
Hyperfine interactions are not considered in our treatment. Several studies have pointed out the effect of hyperfine coupling in cold and ultracold molecular dynamics, mainly in the inelastic processes of non-reactive systems such as YbF+He [45], NH+Mg [46], MnH+He [47], or in reactive systems such as Na 2 +Na in which only the atomic hyperfine state is modified while the molecular hyperfine state remains unchanged [48]. Other studies involving H 2 at higher temperatures for the interpretation of astrophysical spectra, such as CN or HCN interacting with para-H 2 [49], or HCl with para-or ortho-H 2 [50], give approximate methods to account for hyperfine interactions. However, like in our case, the hyperfine structure of H 2 is not considered. The H 2 hyperfine splittings in its electronic ground state X 1 Σ + g are negligible, so that the H 2 +D scattering will be dominated by the hyperfine structure of D (327 MHz or 15.7 mK). Since the focus of this work is the structure found in the range of 100 mK to roughly 10 K, the effect of hyperfine interaction can safely be neglected in our study of nuclear symmetry effects. However, at sub-millikelvin temperatures, the hyperfine effects will become relevant, in particular if a resonance is located in the regime of extremely low energies.

Results and discussion
We give here an overview of the influence of the nuclear spin state on the scattering properties (reaction, quenching, and elastic processes) of H 2 with D. The nuclear singlet spin state of para-Hydrogen with (I = 0, M I = 0) allows only even-j rotational states, while the nuclear triplet spin state of ortho-Hydrogen with (I = 1, M I = 0, ±1) allows only odd-j rotational states. We consider here only the lowest rotational states j = 0 and j = 1 for the initial para-and ortho-H 2 , respectively. In what follows, we show results for scattering energies ranging from 1 µK to 100 K; as mentioned in the previous section, hyperfine interactions will become relevant for energies below roughly 15 mK, and thus our results below that energy are approximate.

Overview
We computed the energy-dependent rate constant K Q n (E) for quenching and K R n (E) for reaction given by Eqs.(3) and (4), respectively, for scattering of D on initially prepared para-H 2 (v, j = 0) and ortho-H 2 (v, j = 1). Here, n = (a = 1, v, j = 0) for para-H 2 and (a = 1, v, j = 1) for ortho-H 2 , with a = 1 labeling the arrangement H 2 +D (a = 2 stands for the arrangement HD+H). Figure 3 shows the overview of the results for the six lowest vibrational levels v of H 2 . Note that, for v = 0, only the reaction channel is opened for both para-and ortho-H 2 . Overall, we find that K v (E) of para-and ortho-H 2 share many features; for example, the position of resonances occur at roughly the same energy, since they are due to shape resonances in the entrance channel. However, we also notice differences, such as the values of K v (E) and the absence of certain resonant features in a few cases.
We find structures in the range of scattering energies corresponding roughly to 100 mK and 10 K. These features are due to low partial waves, with the more pronounced ones found at very low energies, near 15 mK, arising from p-wave resonances. At higher energies (above 10 K) no structure or resonance feature can be found, because the much higher centrifugal barrier (for higher partial waves, ℓ > 3) in the entrance channel lifts the shallow attractive well entirely (see Fig. 4 for D+para-H 2 in v = j = 0), and thus removes any possible bound states (i.e., no van der Waals complexes exist for large ℓ). For collision energies between 10 K and 100 K, the rate constant follows a simple behavior (nearly constant). For much higher energies (above 100 K, not shown here) an exponential increase is typical (especially for lower vibrational initial states). Finally, we also notice that both reaction and quenching rate constants for ortho-H 2 are slightly larger than the corresponding rates for para-H 2 . This is to be expected, since internal energy of ortho-H 2 initially in j = 1 is slightly higher than that of para-H 2 in j = 0, thus reducing the effect of the reaction barrier, and hence increasing the reaction and quenching processes.  In the next sections, we describe the results for both para-and ortho-H 2 in more details, but we first make a general remark about the magnitude of the rate coefficients in Fig. 3. As we discussed in our previous work [41], tunneling through the reaction barrier for v = 0 and v = 1 is responsible for the small values of K R (E) in the energy range considered here; in particular, for v = 1, we have K R < K Q . However, for v = 2 and higher, tunneling is no longer important; moreover, due to the restriction on the rotational levels for para-and ortho-H 2 , and also because of the increased reduced mass of HD, there are many more reaction channels (v ′ , j ′ ) for the product HD than for quenching. Consequently, we have K R > K Q for v ≥ 2.

Para-H 2
In the case of para-H 2 with j = 0, we have ℓ = J only, which reduces the summation over ℓ to a single term. The cross section (1) becomes The resulting rate constant K v (E) = v rel σ for both reaction and quenching processes are shown in Fig. 5. For scattering energies ranging from 1 µK to roughly 100 K, only J = 0 to 4 contribute significantly to the cross sections and rate constants: these individual contributions are also shown. We now discuss each individual initial vibrational level v for both reaction and quenching.
For v = 0, only reactions into HD(v ′ = 0, j ′ ≤ 2) are possible for this range of energies: quenching to lower states is not possible since H 2 is already in its lowest state, and reactions to higher v ′ or even higher j ′ would correspond to excitations and require larger scattering energies; e.g., about 500 K to excite H 2 (v ′ = 0, j ′ = 2) and 350 K to excite HD(v ′ = 0, j ′ = 3). Due to the presence of a barrier (see Fig. 2), reactions take place through tunneling, leading to a very small rate constant. We find a large resonant feature near 200 mK arising from the J = 1 contribution (or p-wave ℓ = 1 since ℓ = J for para-H 2 initially in j = 0). As shown in Fig. 4, the resonance occurs because the centrifugal interaction pushes the bound state of the van der Waals complex near the scattering threshold. Below 10 mK, the main contribution is from J = 0 (or s-wave ℓ = 0), while J = 2 (or d-wave ℓ = 2) impacts K(E) at roughly 3-4 K. At higher energies, the other partial waves add up (including possible excitations), but no structure is present.
For v > 0, both reaction and quenching processes have qualitatively the same behaviors, with similar contributions from the different J. Except for v = 1, the reaction rate constant is larger than the quenching rate constant, simply reflecting the larger number of exit channels to form HD as compared to the number of channels to remain as H 2 , roughly 2-4 in that range of levels v (see Fig. 1). Level v = 1 is slightly different because the presence of the reaction barrier is still relevant for a sizable range of relative orientations: as shown in Fig. 2, the reaction barrier disappears for angle larger than 112 o . The presence of the barrier for angle less than 112 o implies a reduced reaction rate, since it requires tunneling to occur.
The J = 0 contribution, which corresponds to s-wave (ℓ = J = 0) scattering for H 2 in initial states j = 0, dominates the rate constant at very low energies, as is well known. However, note that for energies above 100 mK, the s-wave contribution has a minimum; thus, higher J-contributions do become dominant as the scattering energy increases. It begins with J = 1 corresponding to p-wave scattering (ℓ = J = 1), followed by J = 2 (d-wave scattering with ℓ = J = 2), and so on. The exact energy at which these higher J contributions become important depends on the energy surface and the position of levels.
For both v = 0 (reaction only) and v = 1, the main feature is a scattering resonance due to J = 1 (p-wave scattering with ℓ = J = 1), arising from a quasi-bound state in the van der Waals H 2 · · ·D complex. From v = 0 to v = 1, the quasi-bound level gets closer to the threshold, leading to a shift of the resonance to a lower scattering energy (from 200 mK for v = 0 to 15 mK for v = 1). This quasi-bound level becomes bound for v = 2, removing the resonant feature from the rate constant, and becomes more deeply bound as v increases, moving the J = 1 contribution to higher scattering energies. The J = 2 contribution becomes important at higher energies in the range of a few K, and shifts slightly to lower scattering energies as v grows; it starts to be important for v = 2, and becomes the key feature at higher v. Similarly, J = 3 (f -wave scattering with ℓ = J = 3) starts contributing at roughly 10 K, with the exact position of its maximum contribution shifting slightly with v; it leads to a perceptible feature only for the highest v, its contribution being otherwise a simple addition to all higher J giving the total rate constant. This is illustrated by J = 4 (g-wave scattering with ℓ = J = 4) which simply adds up to the total rate constant without leading to features.

Ortho-H 2
In the case of ortho-H 2 with j = 1, since ℓ = |J − j|, · · · , J + j, three partial waves contribute to a given J, except for J = 0. In fact, for J = 0, we have only the ℓ = j = 1 p-wave scattering contribution, while J = 1 includes the ℓ = 0, 1 and 2 (s, p, and dwaves) contributions, J = 2 the ℓ = 1, 2, and 3 (p, d and f -waves) contributions, and so on. The cross section (1) becomes The resulting rate constant K v (E) = v rel σ for both reaction and quenching processes are shown in Fig. 6. Many of the features and behaviors described for para-H 2 apply for ortho-H 2 as well; only J = 0 to 4 (also individually shown) contribute significantly to the cross sections and rate constants for scattering energies ranging from 1 µK to roughly 100 K (recall that the hyperfine interactions omitted here will become relevant below roughly 15 mK). For the same reason as for para-H 2 , v = 0 leads to reaction only, since quenching to lower states of H 2 is not possible (para-ortho mixing due to nuclear spin-rotation coupling is extremely weak and is neglected here), and the value is very small because the reaction must occur via tunneling through a barrier. Similarly, except for v = 1 for which the reaction barrier still influences the rate, the reaction rate constant is larger than the quenching constant for v > 1. However, the exact features are different from those of para-H 2 . This can be understood in terms of the different couplings between the various states. We now describe in more detail the results for various initial v. Except for J = 0, all other J' include contributions from three initial partial waves ℓ. However, even if a given partial wave ℓ participates in various J-contributions, their effect is slightly different since each total angular moment J leads to different coupling strengths. For example, the p-wave ℓ = 1 initial partial wave occurs in the contributions of J = 0, 1, and 2. For v = 0 (reaction only), the p-wave resonance appears only in the J = 2 contribution, while the p-wave scattering corresponds to a bound state in J = 0, and has almost no influence over the s-wave dominated J = 1 contribution. We note that the higher J are usually favored because of the 2J + 1 factor appearing in the expression for the cross section (6). Still for v = 0, the s-wave scattering is present only in the J = 1 contribution, while the higher J = 3 and 4 contain ℓ = 2, 3, and 4, and 3, 4 and 5, respectively. As in para-H 2 , the effect of partial waves higher than d-wave (ℓ = 2) is hidden in the total rate; they do not produce noticeable structures.
For all other v, the J = 2 does not contain a p-wave resonance; traces of it appears only in the J = 1 contribution of v = 1 and 2. In more detail, for both v = 1 and 2, J = 0 includes only p-wave scattering without resonance. Also, J = 1 contains the s, p, and d-wave scattering; the position of the maximum for p-wave is shifted relative to the other J, illustrating the sensitivity of the position of the bound level in the van der Waals complex to the couplings. This is less so for the other partial waves since the bound states are not close to the threshold. The J = 2 contribution includes results from p, d, and f -scattering, with the ℓ = 1 and 3 components being dominant, while the l = 2 is obscured by them. This occurs for all v. As for v = 0, the higher contributions J = 3 and 4 are dominated by ℓ = 2 and ℓ = 3, respectively, with signals from ℓ = 3 and higher being hidden in the total rate constant. For v > 2, the p-wave feature of J = 1 is absent, preventing a double-peak structure in the total rates.

Elastic cross sections
Here, we describe very briefly the results for the elastic cross section, which reads For para-H 2 initially in j = 0 (with ℓ = J) we have and for ortho-H 2 initially in j = 1 ( with ℓ = |J − 1|, · · · , J + 1) The corresponding results are shown in Fig. 7 for v ≤ 5. As before, there is some structure in the range of roughly 10 mK to about 10 K. The cross sections for both paraand ortho-H 2 +D have the same order of magnitude for all v, except in cases for which the p-wave scattering leads to a sharp feature. This occurs for v = 0 and 1 of both para-and ortho-H 2 , and for v = 2 of ortho-H 2 . Overall, the elastic cross sections have basically the same value around 2000 a 2 0 , with very small variation about that value. This is expected, since the entrance channels are very similar for all those initial states (see Fig. 8), as opposed to the reaction and quenching processes where the position of the exit channels and the coupling strengths dictate the large range of constant rate values.

State-to-state results
In this section, we give an overview of the richness of the state-to-state processes. We first discuss the branching ratios of product formation from each initial vibrational state v of H 2 into individual vibrational levels v ′ (summed over j ′ ) for quenching and reaction separately, together with the branching ratio into all product channels (i.e., quenching and reaction combined). This is followed by an example of a typical case for rotationally-resolved rate constants. In all cases, we contrast the results of para-and ortho-H 2 .
3.5.1. Products branching ratios. We first give the branching ratios for para-H 2 initially in the (v, j = 0) state. They are depicted in Fig. 9; the left panels show ratios for quenching alone (dashed lines), the middle panels for reaction alone (solid lines), and   the right panels for the combined inelastic processes (quenching: dashed lines, reaction: solid lines). We omit the initial v = 0 panels since there is no quenching possible and all reactions end up in v ′ = 0 of HD, and the panel for v = 1 of quenching alone, since all quenching goes into v ′ = 0 of H 2 in this particular case. We find that the branching ratios are basically insensitive to the scattering energy within the 1 mK to 10 K range, even though there is structure in the corresponding rate constant shown in Fig. 5. This simply signifies that all the flux into the various product channels are changing in unison with the appearance of scattering features; these are due to shape resonances and higher partial waves in the entrance channel that raise all scattering amplitudes in the same proportion. In detail, starting from v = 5, the largest fraction (32%) of quenching is in the uppermost level v ′ = 4, followed by 28% in v ′ = 3, 19% in v ′ = 2, 12% in v ′ = 1, and 9% in v ′ = 0. For reaction alone, v ′ = 5 dominates the reaction products formed with 41%, followed by 25% in v ′ = 4, 14% in v ′ = 3, while each of v ′ = 2, 1, and 0 have between 6% and 8%. When all product channels are combined, it becomes clear that the most probable products stem from the reaction forming HD in v ′ = 5 (29%), and 4 (17%), with the rest of reaction and quenching being comparable at 10% or less. As discussed in Section 3.6 below, this is to be expected, since the number of possible exit channels is larger in the reaction arrangement than in the quenching arrangement [41].
The same general behavior takes place for v > 1: for v = 4, both v ′ = 3 and 2 have the same quenching ratios (∼32%), while v ′ = 4 leads the reaction ratios with 50%, and the most probable products are HD in v ′ = 4 (39%) and 3 (19%), the remaining being distributed over the other channels. For v ′ = 3, quenching is largest into v ′ = 2 (39%), reaction into v ′ = 3 (60%) and 2 (22%), and the most probable products are HD in v ′ = 3 (51%) and 2 (18%). Similarly for v = 2, 62% of the quenching ends up in v ′ = 1 and 38% in 0, while 64% of reaction goes into v ′ = 2, 23% in 1, and 13% in 0: most of the products are HD in v ′ = 2 (52%) and 1 (19%), with the rest distributed over the remaining channels. Now, the behavior for v = 1 reflects the considerable effect of the barrier on the scattering: 80% of the reaction into HD produces the lowest level v ′ = 0 as compared to 20% in the highest open level v ′ = 1. The reaction fractions are much smaller than the quenching fraction; specifically, 74% for quenching into H 2 (v ′ = 0), 20% and 6% into v ′ = 0 and 1 of HD, respectively. Finally, for v = 0 (not shown in the figure), only reaction into v ′ = 0 of HD is possible.
In general, we find that the products with the largest (although not dominant) branching ratios are into highest v ′ allowed: v ′ = v for reaction, and v ′ = v − 1 for quenching. Also, reaction products are more likely than quenching, which is expected since there are more possible exit channels available in the HD reaction arrangement as compared to the H 2 quenching arrangement. These conclusions apply to all initial states v, except when tunneling through the barrier is important, as discussed in Section 3.6.
The same overall conclusions apply to the case of ortho-H 2 initially in (v, j = 1) shown Fig. 10: in general, the branching ratios are rather insensitive to the scattering energy in the range of 1 mK to 10 K, the largest branching ratios are into products with the highest v ′ allowed (with reaction to form HD more probable than quenching of H 2 ), and tunneling through the barrier affects the lowest initial levels v = 1 and 0. However, we note that the branching ratios are changing for both v = 2 and 1. In v = 2, the quenching into v ′ = 1 increases from 60% to 85% near 50 mK, while quenching into 0 decreases accordingly from 40% to 15%. For the reaction products, the branching ratios into v ′ = 2 and 0 decrease slightly at the same energy, while it increases into v ′ = 1. The net effect in the combined branching ratios shows a sizable increase of the quenching into v ′ = 1 (reaching 25%) while reaction into v ′ = 2 decreases to 40%. For the initial v = 1 state of H 2 , the net result is an increase of the branching ratio into H 2 (v ′ = 0) from 60% to 75% near 200 mK, while the reaction into HD in v ′ = 0 and 1 are both decreasing accordingly. This is an example showing that para-and ortho-H 2 have different behaviors. To better understand the origin of those variations, we investigate rotationally-resolved rate constants for a specific case, namely v = 2, in the next section.

Rotationally-resolved rate constants for H
Because of the large number of exit channels for a given entrance channel, we restrict this analysis to a typical entrance channel showing structure, while still having a manageable amount of exit channels. We selected to show the state-to-state inelastic (reaction and quenching) rate constant for para-and ortho-H 2 initially in v = 2. Other entrance channels would exhibit the same type of behavior, though the exact details will vary. However, as will become obvious below, the rapidly growing number of channels make it difficult to show those results.
We begin with para-H 2 initially in (v = 2, j = 0), and show the state-to-state rate constant for quenching in Fig. 11, and for reaction in Fig. 12. For quenching, we show the results in v ′ = 1 (top left panel) and v ′ = 0 (bottom left panel): in each of exits levels, we show the individual contributions of rotational states j ′ to the total rate constant (with the quantum numbers labels following the same ordering as the partial rate constants). To the right of those two panels, we show the contribution of individual J total angular momenta for each channel j ′ . In all cases, all axes have the same range.
For v ′ = 1, the most important contribution comes from j ′ = 2, followed by 4, 0, 6, and 8: each have a very similar energy dependence. The individual contribution of each J are shown in the right panel: they all have the same energy dependence, with the main contribution arising from J = 1 and 2 (at higher energies). The results for v ′ = 0 are very similar, although the ordering of the j ′ contributions changes a bit and an additional j ′ = 10 appears.
For reaction (see Fig. 12), very similar behaviors are found, although the number of exit channels v ′ now include 2 as well, and both odd and even j ′ are allowed. The exact ordering of contributions differs from the quenching case, but the overall system exhibit the same dependence. The same information for the case of ortho-H 2 initially in (v = 2, j = 1) is shown in Fig. 13 for quenching and Fig. 14 for reaction. For quenching alone, the left panels in Fig. 13 depict the rate constant for v ′ = 1 and 0, each for the various odd j ′ contributions to the total value. Although these contributions show similar dependence on the scattering energy for both v ′ , one can notice that specific j ′ have larger increases than others near the resonant feature at 50 mK. For example, in v ′ = 1, j ′ = 7 becomes dominant at the resonance, even though it is the smallest contribution away from it. The same is true to a lesser extent for the other j ′ . This is more apparent in v ′ = 0, where not only j ′ = 7 (which becomes dominant) but also 9 and 11 have sharper increases at the resonant structure. The origin of those variations is due to the different couplings between the various states, as discussed in Section 3.3. Basically, this feature arises from a resonance in the p-wave (ℓ = 1) initial partial wave, which occurs in the contributions  Figure 12. Same as Fig. 11 for reaction. The quantum numbers v ′ and j ′ label the rovibrational channels of the product HD.
of J = 0, 1, and 2, and their effect is slightly different since each total angular moment J leads to different coupling strengths. For reaction alone, shown in Fig. 14, we have a very similar story, though there are many more possible product states (including v ′ = 2, and even and odd j ′ ). Again, the structure is due to the p-wave in the entrance channel, and the contributions of the various j ′ increase different near the resonant feature at 50 mK. For example, for v ′ = 2, j ′ = 3 and 4 have sharper increases, while the same is true for j ′ = 2, 6 and 8 of v ′ = 1. These "unequal" changes in the contributions of the total rate constant in each product v ′ explains the difference between the state-to-state scattering of para-and ortho-H 2 .
Rotational product distributions for a given scattering energy E, i.e., the probabilities p n a ′ ,v ′ (j ′ , E) to form a product in a rotational state j ′ (in vibrational state v ′ ) in arrangement a ′ can be obtained by simply dividing the rate constant for this specific channel by the total rate constant for the corresponding product channel v ′ , where n = (a, v, j) stands for the entrance channel (v, j) in the reactant arrangement a. By summing over all open exit channel within one arrangement a ′ , one can find the rovibrational distribution for quenching (Q: a ′ = a) or reaction (R: a ′ = a), or simply p n Q(R) (v ′ , j ′ , E) = K Q(R),v ′ ,j ′ ←n (E)/K Q(R) n (E). Finally, the distribution over all possible channels can be obtained by dividing K Q(R),v ′ ,j ′ ←n (E) by total rate constant K n (E) = K Q n (E) + K R n (E).

Statistical model
Scattering problems with rearrangement are computationally demanding when a full quantum mechanical approach is used; thus, approximate models have been developed, often based on statistical formalism and semi-classical treatments [51,52]. Here, we compare the results of our full quantum computation with a statistical treatment based on the phase-space theory (PST) [51].
Because the entrance channel plays a special role at low scattering energies, it is advantageous to interchange the order of the J and ℓ summations in Eq.(1), and rewrite the cross section as Following Miller [51], we can write the T -matrix element T J n ′ ℓ ′ nℓ for an inelastic process (i.e., n = n ′ ) as where P J n,ℓ (E) is the energy-dependent capture probability in the channel n for a given partial wave ℓ of a particular total angular momentum J, with similar definitions for P J n ′ ,ℓ ′ (E) and P J n ′′ ,ℓ ′′ (E). Labeling the initial arrangement H 2 +D by a, the rate constant K n ′ ←n (E) = v rel. σ n ′ ←n (E) with v rel. = 2E kin /µ a =hk n /µ a is simply where ℓ ′ = |J − j ′ |, . . . , J + j ′ and ℓ ′′ = |J − j ′′ |, . . . , J + j ′′ , respectively. As in the case of Eqs. (3) and (4), we obtain the total quenching and reaction rate constant by summing over the appropriate channels and arrangements where with P J tot (E) = P J Q (E) + P J R (E) + ℓ ′′ P J n,ℓ ′′ (E). Due to the small mass of our system, we use the PST capture model to determine P J n ′ ,ℓ ′ (E) in the product channels only, where the kinetic energy is relatively large and the effect of near-threshold states below dynamical barriers may be neglected [52], i.e., The van der Waals coefficient C 6,n ′ and reduced mass µ a ′ appearing in the expression of the centrifugal barrier height V b,n ′ depends in the particular channel n ′ = (a ′ , v ′ , j ′ ). For P J n,ℓ (E) in the entrance channel, we use the low-energy expression P J n,ℓ (E) = A J n,ℓ (E)k 2ℓ+1 where the coefficient A J n,ℓ (E) is related to Jost functions [53].
At low-energy, ℓ ′′ P J n,ℓ ′′ (E) is small and we can omit it in evaluating P J tot (E) ≃ P J Q (E) + P J R (E). Furthermore, P J Q , P J R , and P J tot are essentially E-independent at lowenergy; figures 15 and 16 show their variation with J for the initial state n = (v, j = 0) for para-H 2 and n = (v, j = 1) for ortho-H 2 , respectively. Although P J Q , P J R , and P J tot all increase with J, the ratios P J Q /P J tot and P J R /P J tot are nearly constant (roughly 73% reaction and 27% quenching), allowing us to further simplify eq (13) and (14) as where P J Q(R) /P J tot ≃ P Q(R) /P tot does not depend on J. In particular, their values for J = 0 is simply given by the corresponding number of open channels N Q for quenching and N R for reaction, so that P Q(R) /P tot ≃ P J=0 Q(R) /P J=0 tot = N Q(R) /N tot with N tot = N Q + N R . Although the rate constants for quenching and reaction processes depends on the details of the entrance channel via the sums over ℓ and J, their ratio will not; the double-sum cancels out, and we find Figure 17 compares the results obtained using this statistical model and the full quantum treatment, We find that the agreement is roughly within 10-20% for v ≥ 2, but quickly becomes worse when the reaction barrier plays an important role for very low v. Overall, this statistical model based on PST approximately gives the right branching ratio between quenching and reaction processes. As discussed in previous sections, the magnitude of the rate constants at low scattering energies are mainly given by the s-wave contribution: figures 5 and 6 show that this is the case for most of initial v up to 10 mK and higher. For para-H 2 in j = 0, the s-wave ℓ = 0 implies J = 0 only, so that (2J + 1)/(2j + 1) = 1, while for ortho-H 2 in j = 1, ℓ = 0 implies J = 1 only, so that (2J + 1)/(2j + 1) = 1 as well. In both cases, both sums over ℓ and J in Eq. (18) reduce to a single term; For ℓ = 0 with J = j, the coefficient A J=j n,ℓ=0 (E → 0) is simply related to the imaginary component β v,j of the complex scattering length a v,j = α v,j − iβ v,j [54] in the entrance channel n = (a, v, j) of the initial arrangement a=H 2 +D and total angular momentum J = j [55], such that the rate constant takes the simple form  Figure 15. Top panels: P J Q(R) and P J tot for para-H 2 as a function of J for different entrance channels n = (j = 0, v) with v = 1, . . . , 5 (v = 0 is not shown since only reaction is allowed, so that P J R /P J tot = 1). Lower panels: corresponding ratios P J Q(R) /P J tot . Also shown is the number of molecular states opened to quenching and reaction, N Q , and N R respectively, and their ratios. These are J-independent and equal to the P J=0 Q(R) values. Table 1 contains the values of α v,j and β v,j , which were extracted from the full quantum calculations [56]. They can be compared to an approximate value based on quantum reflection [57,58,59], which is often used for barrierless reactions, Using the value C v=0,j=0 6 = 7.053 a.u. [44], together with µ a ≃ 1836 a.u., we obtain β ≃ 6.06 a.u., which is many order of magnitude larger than the values of β from the full quantum treatment. This illustrates well the fact that, although the statistical model gives approximately the right branching ratio of quenching to reaction processes, the approximation in Eq. (23) cannot be used to obtain the magnitude of those processes. Although such an approximation can be applied to barrierless chemical systems, the reaction barrier reduces the value of β considerably, as shown in Table 1.   Figure 17. Comparison of the branching ratios of quenching and reaction processes between the full quantum treatment K Q(R) n /K tot n (dash lines) and the statistical model N Q(R) /N tot (solid lines) for para-H 2 and ortho-H 2 . Results for v = 0 are not shown, since only reactive scattering is allowed. Ratios agree within 10-20%, except for v = 1 for which the effect of the reaction barrier is more important. Table 1. Real and imaginary contributions to the scattering length a v,j = α v,j − iβ v,j for para-H 2 (j = 0) and ortho-H 2 (j = 1) extracted from the quantum results. D + para-H 2 (v, j = 0) D + ortho-H 2 (v, j = 1)

Conclusion
The nuclear spin symmetry restricts the possible rotational states of H 2 , leading to para-H 2 with even j and ortho-H 2 with odd j. We investigated the effect of these restrictions on the scattering of H 2 with D, which can lead to elastic collisions, to quenching H 2 , or reaction to produce HD. We computed these various processes for the six lowest vibrational levels v of H 2 . We found structures in the energy range corresponding to 100 mK and 10 K, due to low partial waves. Except for cases involving a sharp p-wave resonance, the results for para-and ortho-H 2 are very similar, with structures due to d or f -waves located at roughly the same energies in both cases, and with rate constants of comparable magnitude. There are marked differences for the cases in which a p-wave resonance is present, namely the three lowest v levels. For example, in the reactive case, while the sharp peak is present for both para-and ortho-H 2 in v = 0, it is absent for ortho-H 2 in v = 1 and para-H 2 in v = 2. The same is true for the quenching rate constant (though there is no quenching for v = 0, since it is the ground state). We also find the same overall behavior for the elastic cross sections. We also note that the exact position of the resonance is slightly affected by the initial value of j.
We also investigated the state-to-state product formation, and found that the branching ratios are not sensitive to the scattering energy in the range considered here for para-H 2 . The same is true for ortho-H 2 except for the lower initial levels (v = 1, 2), where the p-wave resonance affects the branching ratios. By examining the rotationally-resolved rate constants of a specific initial case, namely v = 2, we can trace this dependence to sharper increases of the inelastic (either quenching or reaction) rate constant for certain rotational products near the resonance. These variations are due to the different coupling strengths of each total angular momentum J = 0, 1, and 2 contributing to p-wave (ℓ = 1) initial partial wave for ortho-H 2 , as opposed to the case of para-H 2 , when only J = 1 contributes.
The range of energies where the structure is present should be reachable experimentally, and thus could probe the effect of the nuclear symmetry on reaction (or quenching). In addition, since the occurrence of p-wave resonances is very sensitive to the exact shape of the potential energy surface, detecting any structure in this benchmark system will provide a unique tool to compare with ab initio quantum chemistry calculations. Finally, based on our results, we suspect that the effect of the nuclear spin symmetry on cold and ultracold chemistry will be significant in systems displaying resonant features, and less so otherwise. However, if the nuclear spin is involved in the dynamics of the system, e.g., through an interaction term in the Hamiltonian, as in the studies of Refs. [45,46,47,48,49,50], the nuclear spin symmetry could play an even more important role. Finally, we compared the results of our full quantum treatment to a statistical model, and found that although the model reproduces the proportion of quenching and reaction processes within 20% in most cases, it fails to give the magnitude of the rate constants.