On the use of stereodynamical effects to control cold chemical reactions: the H + D$_{2}\longleftrightarrow$ D + HD case study

Quantum calculations are reported for the stereodynamic control of the H + D$_{2}\longleftrightarrow$ D + HD chemical reaction in the energy range of 1-50 K. Stereodynamic control is achieved by a formalism similar to that reported by Perreault et al. [Nature Chem. 10, 561 (2018)] in recent experimental works in which the alignment of the molecular bond axis relative to the incident relative velocity is controlled by selective preparations of the molecule in a specific or superposition of magnetic projection quantum numbers of the initial molecular rotational level. The approach presented here generalizes the experimental scheme of Perreault et al. and offers additional degree of control through various experimental preparation of the molecular alignment angle. Illustrative results presented for the H + D$_{2}$ and D + HD reactions show significant control with the possibility of turning the reaction completely on or off with appropriate stereodynamic preparation of the molecular state. Various scenarios for maximizing and minimizing the reaction outcomes are identified with selective preparation of molecular rotational states.


I. INTRODUCTION
Stereochemistry 1-6 along with coherent and optimal control of molecular processes 7-10 are a manifestation of a long standing goal of controlling chemical reactions. The problem has been addressed independently from the perspective of a variety of research fields; from exploiting ultra-fast light-induced electron dynamics -within an attosecond time scale 11 -to an ultra-slow sub-Kelvin nuclear dynamics, in the context of either cold or ultracold molecular collisions 12,13 .
In a recent series of seminal experiments, Perreault and co-workers have demonstrated a relatively simple setup in which a co-expansion of two colliding molecules (HD + D 2 and HD + H 2 ) in the same beam stream enable cold collisions with relative kinetic energies of about 1 Kelvin 14,15 . The technique, combined with Stark-induced adiabatic Raman passage (SARP) 16 allows selective preparation of the alignment of the diatomic internuclear axis relative to the incident velocity vector. This is made possible by a light-assisted spatial polarization of the diatomic angular momentum vector j combined with selection of specific j-projections, say m, along a given quantization axis of choice 17,18 . As a consequence, the lower collision energies provided by the co-expansion can selectively discard higher and undesired 14 external collisional orbital angular momentum states, , whereas SARP can selectively prepare a specific internal rovibrational state with a given m or superposition of m states. More recently the same methodology was also applied to cold He + HD and He + D 2 collisions with the latter revealing strong signatures of an = 2 partial wave resonance 19,20 . While these experiments have so far been used to study only inelastic processes, they open up new avenues for possible stereodynamic control of cold reactive collisions. Theoreti-cal prospects on that matter have been available for a while 1,2,21 , and revisited recently in the context of coherent control 22 .
It is worthwhile to highlight important aspects of these experiments, beyond requiring a simpler experimental setup compared to current molecule cooling and trapping techniques: (i) it allows controlled low-kinetic energy collisions of spatially-polarized non-polar molecules in a complete field-free environment; (ii) the production of colliding partners described by a single internal quantum state eliminates the need of averaging any observable over a large number of initial states; (iii) due to (i) and (ii), the actual interaction potential can be addressed experimentally as it is, within such physical conditions, mainly driven by the anisotropy caused by the target-projectile relative orientation; and (iv), unveils marked stereodynamical effects previously constrained to scattering experiments using polar molecules 23 .
The experiments of Perreault et al. have been accompanied by the theoretical works of Croft et al. 24,25 on HD + H 2 and Morita et al. 26,27 on He + HD collisions. Further, Jambrina et al. 28 and Morita et al. 29 have explored the control of isolated partial wave resonances through stereodynamic preparations of HD + H 2 and HCl + H 2 collisions, respectively, as in the experiments of Perreault et al. In a recent work, Jambrina et al. examined stereodynamic control of CO + HD collisions in which the HD molecule is prepared in a stereodynamic state but both molecules undergo a change in rotational angular momentum 30 . In this work we extend these ideas to low-energy reactive collisions taking the prototypical chemical reactions of H + D 2 ←→ D + HD with polarization of the initial rotational state j of D 2 and HD. Thus, in what follows, we address questions such as (i) the role of bond axis orientation of D 2 (HD) relative to the initial relative velocity; (ii) whether it is possible to orient the reactants arXiv:2111.08155v1 [physics.chem-ph] 16 Nov 2021 in a manner to enhance or suppress the production of HD (D 2 ); (iii) can we selectively induce the formation of products in a given internal state by controlling the orientation/alignment of the reactants during the collisions? (iv) identify any propensity rules that accompany these changes. These questions are addressed through a systematic analysis of typical scattering characteristics upon specific choices of so-called stereodynamical parameters -defined in more details below. Collisions are modeled using a time-independent quantum reactive scattering formalism as implemented in the APH3D code 31-33 that has been used extensively to model atom-diatom reactive collisions, combined with the legacy BKMP2 potential energy surface (PES) for the H 3 system 34 . As we shall see below, we have found strong evidence that supports a nearly-total stereodynamic control of a cold chemical reaction, from signifcant enhancement to near complete suppression.
Since the early works of Kuppermann et al. 35 the H + H 2 exchange reaction and its isotopic variants have been extensively studied in the literature both experimentally and theoretically and we shall focus mainly on aspects related to stereodynamic preparations and collisional outcomes. Interested readers are referred to recent works of Goswami et al. 36 who reported wavepacket calculations of the H + H 2 exchange reaction; Desrousseaux et al. 37 on rotational excitation in H + HD collisions at temperatures in the range of 10-1000 K; recent prediction of a Feshbach resonance on H + HD collisions by Zhou et al. 38 ; and quasi-classical trajectories simulations by Bossion et al. 39,40 . Kendrick and co-workers have reported extensive studies of geometric phase (GP) effects in cold and ultracold H + H 2 , H + HD and D + HD collisions with vibrationally excited reactants that revealed a new interference mechanism that controls reactivity in the ultracold regime [41][42][43] . Non-adiabatic effects in the H + H 2 and H + HD reactions with vibrationally excited H 2 and HD molecules for vibrational levels υ = 4-9 have also been reported by  . Several theoretical and combined experiment-theory studies on H + D 2 are also available 32,47-54 .
The paper is organized as follows: Section II provides a brief outline of the theoretical formalism for the stereodynamic preparation. This is followed by results and discussion in Section III and a summary and conclusions in Section IV.

II. THEORETICAL APPROACH
The methodology used herein is, in many aspects, very similar to that utilized in the previous works of Croft et al. 24,25 and Morita et al. 26,27 except that reactive channels are also taken into consideration as a collision-induced outcome. Thus, upon a H + D 2 (υ, j) reactive collision (or likewise D + HD), the probability of detecting a scattered HD(υ , j ) product (or D 2 ) by a detector placed at Ω = (θ, φ), at a given collision en-ergy E coll , is governed by the state-to-state scattering amplitude 33 where we have multiplied the original scattering amplitude expression of Kendrick, Eq. (22) 33 , by k n /k n and have explicitly expressed the spherical harmonic in terms of the Wigner d-function. We note that the total scattering amplitude is obtained by summing up all of the contributions from each value of the total angular momentum J: Also, the sums over and effectively average over the azimuthal angle φ so that Eq. (1) is independent of φ. The initial and final states are denoted by n = (α, υ, j, m) and n = (α , υ , j , m ) with α denoting the different atom-diatom arrangements, i.e. α = 1 for H + D 2 and α = 2, 3 for D + HD; υjm are the vibration, rotation and rotational projection quantum numbers, respectively. Moreover, the vector-coupling coefficients, C, are the well known Clebsch-Gordan coefficients and the scattering angle θ measures the deflection between the asymptotic velocity vectors of H (relative to D 2 ) and of D (relative to HD); S J n n is the energy-dependent scattering matrix for a given total angular momentum, written in the -representation (J = j + ), and k n is the magnitude of the wave vector associated to the entrance channel n of the reactants, where E n denotes the corresponding diatomic eigenvalue, E is the total energy, and µ is the reduced mass of the reactants.
In the space-fixed (SF) reference frame, used for evaluating the scattering matrix, the coordinate origin lies at the center-of-mass (COM) of the triatomic system. Often, a body-fixed (BF) reference frame is convenient for the scattering calculations in which the BF quantization axis lies parallel to the atom-diatom relative motion coordinate and coincides with the direction of the initial relative velocity vector. Thus, in this set of coordinates, the rotational projection quantum number m exp of the molecular state |αυjm exp , prepared experimentally, corresponds to a coherent superposition of all possible pure projections m, in {|αυjm }, i.e.
where β is the angle between the internuclear axis of D 2 (or HD) and the incident velocity vector. In the experiments of Perreault et al. β denotes the orientation between the polarization lasers and the molecular beam axis; where choices of β correspond to either a parallel orientation (β = 0), known as a H-SARP preparation, or a perpendicular orientation (β = 90 • ), a V-SARP preparation. Here, H-and V-denote horizontal and vertical polarizations, respectively. A cross polarization state, referred to as X-SARP, corresponding to a superposition of m = ±1, was also adopted in their recent work on He + HD and He + D 2 collisions 19,20 . For exploring possible control scenarios, several values of β within [0, 90 • ] are considered below. As discussed by Croft and Balakrishnan 25 , since the experimental results correspond to an integral over the azimuthal angle φ, one can express the differential cross sections (DCS) for any SARP-prepared state as a sum of contributions from each m-state weighted by For the state-to-state DCS, this yields where, n = (α, υ, j, m exp ). Likewise, for the integral cross section (ICS), Figure (1) shows the β-dependence of the Wigner's reduced rotation matrices for the hypothetical cases of j = 1, 2, 3 along with their pair of (m, m exp ) values (negative values are not shown). Upon a simple inspection of Fig. (1), the following aspects are apparent: (i) For β = 0, d j mmexp = 0 whenever off-diagonal elements are concerned, i.e. m exp = m; (ii) similarly, d j mmexp = 1, for diagonal terms at β = 0 and, as a consequence, parallel orientations of the reactants are likely to be composed by terms driven by m exp = m, i.e. |αυjm exp = |αυjm ; (iii) in the particular case of m exp = 0 and β = 0, due to (ii), the expansion of Eq. (4) collapses to a single and isotropic term in which m exp = m = 0; (iv) d 0 00 (β) = 1 for any value of β -not shown in Fig. (1) -and, thus, as expected, processes associated to rotationless reactants are invariant upon changes of β.
For the sake of clarity, it is worthwhile to illustrate at least a few possible linear combinations shown by Eq. (4). Consider a reactant prepared in the rotational state j = 2 of a given vibrational manifold υ with a projection m exp = 0, and subjected to two orientations of its internuclear axis relative to the incoming relative velocity vector: β = 0 and β = 45 • . The first case, β = 0, has been addressed above and the resultant state is simply |αυj; m exp = 0 = |αυj; m = 0 . Often, the m = 0 contribution in the ICS (DCS), for uncontrolled collisions, is dominant but suppressed by the typical average made up of all 2j + 1 initial m-terms. Thus, we may expect that a fixed orientation of β = 0 for a m exp = 0 preparation, in which only m = 0 contributes, to enhance substantially the collision-induced state-to-state branching ratio of interest. If, instead, a β = 45 • alignment is made, the same m exp = 0 preparation is expected to be composed by a broad mixture of all m-terms, namely: which is now mainly dominated by contributions from m = ±1 and m = ±2 with a minor contribution from m = 0. Note that in Eq.(7) (and in similar expressions below) we have used numerical values of d j mmexp (β) matrix elements corresponding to the chosen β values. If a m exp = 2 preparation is considered at the same β = 45 • alignment, the linear combination of m-terms yields In this case, the prepared state is mostly composed of m = +2 and m = +1, and the m = 0 term. The preparation for j = 2, m exp = 0 at the two typical orientations of β = 0 and β = 90 • (V-SARP and H-SARP) has been discussed in great detail by Perreault et al. 14,15 .
The scattering amplitudes in Eq. (1) are computed using a numerically exact time-independent quantum reactive scattering formalism based on hyperspherical coordinates as implemented in the APH3D code 31,32,50,51,53 . Details of the computational formalism have been discussed in great detail in prior works of Kendrick 32,50,51,53 and in recent works of Kendrick et al. 43 on the D + HD and H + D 2 reactions. We shall remark only on some key details here. The hyper-radius ρ describing the radial motion of the triatomic system is divided into an inner region where all three atoms are in close proximity and an outer region where the different atom-diatom arrangement channels are uncoupled. In the inner region the three-body Hamiltonian is expressed in the adiabatically adjusting principal axis hyperspherical (APH) method of Pack and Parker 31 . Smith-Johnson hyperspherical coordinates 32,55 are used in this region where collision-induced re-arrangements are more likely to occur, and Delves hyperspherical coordinates [56][57][58] are used in the outer region. An in-depth discussion of the hyperspherical coupled channel equations and numerical strategies to solve them are available elsewhere 32,50,51,53 . Briefly, the six-dimensional three-body problem is reduced to a set of one-dimensional coupled equations along the scattering coordinate ρ, that is discretized in a grid of N sectors. The eigenvalue problem associated to the remaining five (internal) degrees of freedom is solved adiabatically, at each sector, yielding an effective set of coupled potentials driving the relative motion along ρ. The five-dimensional (5D) eigenvalue problem is solved in the APH region by means of an implicitly restarted Lanczos method for sparse matrices as implemented in the PARPACK library 59,60 . The corresponding eigenvalues within the Delves region are evaluated using a 1D Numerov propagator 61 . Once the eigenvalue problem is solved in both regions for all sectors and the sector-to-sector overlap and transformation matrices evaluated, the resulting set of radial coupled equations are solved using Johnson's log-derivative method 62 , first from ρ min to ρ match . At ρ match the numerical solutions, from the outermost sector of the APH region, are projected onto solutions at the innermost sector of the Delves region. The propagation is continued from ρ match to ρ max , a sufficiently large value of ρ where the interaction potential is negligible. At ρ max all channels (from all arrangements) are numerically decoupled, and scattering boundary conditions are applied to the log-derivative matrix in order to evaluate the scattering S-matrix. This involves projecting the log-derivative matrix onto solutions associated with each asymptotic diatomic state, in appropriate Jacobi coordinates, yielding the S J n n matrix 31 . The procedure is repeated independently for each value of the total angular momentum quantum number J, J-parity P (good quantum numbers in the absence of external forces) and, each j-parity p for the D 2 molecule. It should be noted that the basis sets for both APH and Delves regions are independent of collision energy and, therefore, evaluated only once for each (p, P, J) sets. Also, both the APH and Delves solutions are properly symmetrized with respect to the identical particle permutation symmetry (e.g., the two D nuclei in HD 2 ).
As both H + D 2 and D + HD chemical reactions (as well as other isotopic combinations) have been modeled with the APH3D code in the past 51 , the optimal set of parameters required to extract fairly converged (< 5%) scattering characteristics are well established. In this work we follow closely the parameterization used previously by our group 42,43,50,51,53 . Thus, the APH region ranges from ρ min = 1.6 a 0 to ρ match = 7.8 a 0 with 64 sectors varying logarithmically with a step size of ρ aph = 0.025 a 0 . The Delves region, starting at ρ match , extends to ρ max = 78 a 0 with 350 sectors varying linearly with a step size ρ delves = 0.2 a 0 . Reactions such as H + D 2 are known to require a few dozens of partial waves whenever collision energies exceed 1 eV 51,53 . Here however, we shall simulate the reaction at low and moderate collision energies, between 1 K (about 10 −5 eV), a typical value in the experiments of Perreault and co-workers, and 50 K (about 10 −3 eV). Thus, to reduce computational efforts associated with nonzero initial rotational levels, we limit our calculations to sets of J max = 4, 7 and 9 depending on the entrance channel taken into consideration. Numerical convergence is estimated within 1% and, as we shall see below, most of the interesting physics is captured at resonant energies where the discrete contribution of specific (and often lower) values of J is more relevant.
In order to compute the APH basis set, the 64 sectors between ρ min and ρ match are further subdivided into five groups each of which with increasing maximum values of l and m (not to be confused with the orbital angular momentum and its projection quantum numbers), two parameters that define the primitive basis sets used in the APH region 32 : (i) l max = 103, m max = 190 for ρ = 1.6-2.9; (ii) l max = 115, m max = 214 for ρ = 2.9-3.7; (iii) l max = 123, m max = 232 for ρ = 3.7-4.6; (iv) l max = 135, m max = 250 for ρ = 4.6-5. batches of 50 collision energies at each of the entrance channels described below and every parity of J and j. For the J > 5 case a parity-adapted basis-set is utilized: 2170 (J = 6, 7), 2700 (J = 8, 9) channels for even P ; 1860 (J = 6), 2400 (J = 7, 8) and 3000 (J = 9) channels for odd P .

III. RESULTS
In what follows the collision-induced exothermic reaction (α = α ) and inelastic de-excitation (α = α ) are studied for the reactants in the initial rovibrational state υ = 4, j = 1, 2 (J max = 4), 3 (J max = 7) and 4 (J max = 9); whereas the internal states of the products are chosen to be either υ = υ or υ = υ − 1 as well as j = 0-3. Choices of initial and final states are somewhat arbitrary and intended as representative cases upon which an actual experiment could be designed. However, it should be noted that, as pointed out by Hazra et al. 41 , the effective reaction pathway along the set of adiabatic coupled potentials, for hydrogen exchange reactions, becomes barrierless for the vibrational manifolds above υ = 3, which is the primary motivation underlying the choice of υ = 4 herein, since a natural boost in terms of reactivity is expected for low energy collisions 43 .
From an experimental perspective there is a growing effort towards addressing even higher vibrational states 63,64 in order to explore the inherent delocalization of such weakly-bound systems 65 . Indeed, Perreault et al. have recently reported complete transfer of population of H 2 (v = 0, j = 0) to H 2 (v = 7, j = 0) via the SARP method opening pathways for collisional experiments involving highly vibrationally excited H 2 . From a theoretical standpoint however, the description of highly excited vibrational levels during collisions would likely require a treatment which takes into consideration either (i) couplings to higher electronic states or (ii) geometric phase effects whenever the nuclear motion encircles the conical intersection between the two lowest electronic states. In particular, Kendrick and co-workers have discussed both the sources of geometrical phase effects and their consequences on the title reactions for reactants prepared at υ = 4, j = 0 and below 42,43,50 . The effects are found to play a negligible role for exchange reactions at collision energies above 1 K but to have a somewhat stronger influence on the inelastic outcome in the cold and ultracold regimes.
Previous studies of Croft et al. 24,25 and Morita et al. 26,27 have shown that collision-induced resonances are more sensitive to the relative orientation of the colliding partners. Therefore, we have first scanned collision energies for the various entrance channels mentioned above to survey for any resonances. To our knowledge resonances in the domain up to 50 K above the υ = 4, j = 1 − 4 thresholds of D 2 and HD have not yet been reported. We have identified several instances of these and, given the energy regimes involved, are likely to be shape reso- 32 K and 47 K are mostly dominated by the m = 0 component (panel a) of mid-range partial waves, J = 2-6, whereas the J = 6 term alone is responsible for the resonance at 47 K. Also, it is worthwhile to notice that J = 7-9 provide just a minor correction which suggests that numerical convergence is attained at J max = 6. Yet panels (a)-(c) offer a first glimpse on the stereodynamical effects underlying the reaction where hollow symbols denote the ICS for which the initial rotational state of D 2 , j = 4, is prepared in a specific projection m exp , hence a weighted sum of various 2j + 1 m-terms in −4 ≤ m ≤ 4 contribute, as shown in Eq. (5). Results are shown for internuclear axis orientation β = 0 (black circles) and β = 90 • (magenta circles) with respect to the incoming relative velocity vector. Regardless of the choice of m exp , the orientation of the internuclear axis alone is capable of modulating the ICS and, as we shall see below, the appropriate choice of a given m exp and β combined, can effectively enhance (suppress) the chemical reaction at the resonant energy. The m exp = 0 and β = 0 (panel a) represents the special case, mentioned in the previous section, in which all off-diagonal elements of d j mmexp are zero, whenever m = 0, i.e. d 4 m0 (0) = 0. Therefore, the ICS expansion is reduced to a sole term associated to m = 0, where d 4 00 (0) = 1. In order to better visualize why the same m exp = 0 preparation yields a much smaller cross section at β = 90 • , panels (a)-(b), we invoke Eq. (4) again to inspect the actual linear expansion for the prepared state:  For convenience, in the remainder of this paper, results presented will correspond to a sum over all J-values. However, whenever appropriate, we shall address the various partial wave contributions in terms of the orbital angular momentum quantum number with |J − j| ≤ ≤ (J + j).
Panels (a)-(e) of Fig (3) present the ICSs for each case of m exp ≤ j and β = 0-90 • , whereas panel (f) shows the so-called isotropic ICS, usually evaluated in ordinary studies of reactive scattering and represents a collision with no control imposed upon the rotational state j. The isotropic case corresponds, therefore, to an ICS summed over all final m (in this particular example happens to be only m = 0 due to j = 0) and averaged over all 2j + 1 initial m-values. The peak predicted at about 32 K is mostly composed of partial waves = 0-7 whereas the one at 47 K is mainly composed by even terms, = 2, 4, 6 and 8. All contributions from partial waves higher than = 9 are negligible.
Inspection of panels (a)-(e), unveils a propensity in which reactivity associated to lower projections of j, m exp ≈ 0-1, may be significantly enhanced by near parallel orientations of the reactants, say β ≤ 15 • , whereas higher projection values, m exp > 2, are boosted by higher orientation angles. Almost full suppression of the chemical reaction is attained by preparing the reactants in an appropriate superposition of rotational projection states corresponding to specific orientation of its internuclear axis (with respect the incident relative velocity vector). Typical cases are depicted in Fig (3 Results presented in panels (a)-(e) of Fig. (4) show a strong stereodynamical effect for different choices of m exp and β as for the j = 4 case. It is worthwhile to reiterate that this dynamical feature is entirely driven by the properties of the Wigner coefficients, d j mmexp (β), regardless of other scattering characteristics under consideration, such as collision energy. Thus, at the main resonant peak (≈ 41 K), reactivity is subjected to an enhancement of about a factor of 6 for a m exp = 0, β = 0 preparation and a complete suppression for m exp = 3, 9 K whose composition is dominated by = 1 and = 3 with a minor background contribution from = 2. At about 13 K a second (and more intense) peak is predicted to be formed from = 0, 2 and 4, followed by another occurrence at 20 K ( = 1, 3, 5) and 32 K ( = 4, 6). At energies as low as 1 K, most of the contribution is expected from = 0 and = 2 followed by a minor, yet significant, contribution from = 1. In terms of an actual experimental realization at E coll ≈ 5-10 K, therefore, j = 2 is likely to offer the best conditions when considering H + D 2 collisions.
Cross sections for the inelastic process -panels (a)-(d) of Fig. (6) -are presented to illustrate how the exact same machinery described so far applies in these cases as well. The inelastic cross sections are about a factor of 60 larger than their reactive counterparts, due in part, to the direct and completely barrierless pathway for the inelastic process. Panel (c) of both figures unveil another intriguing aspect, where the first resonant feature (≈ 9 K) is suppressed entirely for a m exp = 2, β ≈ 0 preparation. As before, from Eq. which is strongly dominated by the m = +2 component. However, at about 9 K, for = 1, 3 and j = 2, the most relevant contributions arise from m = 0 and m = ±1 for J = 1-4, whereas for J > 2, contributions from m = ±2 become more relevant above 10 K. Hence, all resonances above 10 K predicted for a D 2 state prepared with j = 2, m exp = 2, β ≈ 0 shall survive mostly due to their m = 2 character. Therefore, these calculations suggests that controlling the relative orientations of reactants can be also exploited in terms of addressing specific resonances, systematically, due to their m-,and J-composition, even if they happen to occur in a band composed of overlapping structures. Moreover, an inspection of the energy dependence of β-dependent ICS presented up to this point reveals an intricate competition among different β values. As the energy varies, the actual m → m branching ratio in the S-matrix is likely to vary whereas, at same time, β controls the weight of various m-terms in a given j, m exp preparation. Morita and Balakrishnan noted a similar observation 26 , in terms of an interference effect due to cross terms in the square of the scattering amplitude, in their analysis of stereodynamical effects in rotationally inelastic collisions of HD and He. We wrap up the discussion of ICSs by presenting results for the reverse reaction, D + HD → H + D 2 , Fig.  (7), in which HD is prepared in j = 3 (J max = 7) with m exp ≤ 3 for all β values discussed previously. We limit our analysis to j = 3 of D 2 to inspect all 16 possible (m exp , m )-cases, where m ≤ 3 (negative values are not shown for the sake of simplicity).
Many aspects of the results highlighted for H + D 2 are also observed here, and some are highlighted below: (i) In panel m exp = 0, m = 0, the maximum enhancement of reactivity is observed for the β = 0 preparation and cross sections are suppressed for higher β values.
(ii) In panel m exp = 3, m = 0: enhancement is found for higher β values.
(iii) In panel m exp = 0 − 3, m = 2: a unique resonant feature at about 15 K is predicted to emerge from the m = 2 case, whereas smaller traces of it are also expected for other m > 0 cases.
(iv) In panels m exp = 0, m = 1 − 3: selective suppression of a specific resonance occurs with increase in m at E coll ≈ 20-40 K as β varies from 0 to 90 • .
The richness, in terms of resonant and stereodynamical features, manifested in the ICS for a collision of D with HD, with HD prepared in j = 3, m exp , offers various handles for stereodynamic control of this reaction, in particular, for those experiments targeting collision energies near 1 K. Similarly, the j = 1 rotational levels from both HD and D 2 are also predicted to yield higher reactivity due to a dominant = 1 and = 2 contribution near E coll ≈ 1 K. This is illustrated in Fig. (8) for the D + HD(υ = 4, j = 1, m exp ) −→ H + D 2 (υ = 4, j = 1, m = 0) reaction (J max = 4). A second resonant feature, whose properties are modulated by m = ±1 and a = 3 partial wave, occurs at about 6 K. This resonance is also very sensitive to a m exp , β = 0 preparation (comparing black curves from panels (a) and (b)). Due to such a marked characteristic and presence of an isolated partial wave occurrence, a collision of D + HD(υ = 4, j = 1, m exp ) at about 6 K, with m exp = 0, β = 0 and m exp = 1, β = 0 preparations, appears to be a good candidate for an experimental study. Panel (c) of Fig. (8) depicts the isotropic cross section, where a smaller resonance is also expected to occur at about 14 K ( = 4). In Fig. (9) we present the partial wave composition for the isotropic ICS (summed over J = 0-4) for the D + HD(υ = 4, j = 1) → H + D 2 (υ , j ) and H + We would like to stress that the marked resonance profiles discussed in this work are a manifestation of a given entrance channel n. Thus, other choices of exit channels, n , shall present a similar energy-dependent behavior and no compromise has been made by choosing a specific υ or j product state.

B. Differential cross sections
Now we consider how the relative orientation of reactants influences the deflection of products towards a given scattering angle θ. To this end we shall focus on the m exp = 0 preparation as it has been shown to provide the largest enhancements, upon proper choices of β. As energy-dependent measurements of scattering processes are inherently more complicated, we shall also fix the collision energies to 1 K, 5 K and 10 K.  Fig. (11) present a similar result for the j = 2 and j = 0 case. As before, for comparison purposes, circles are used for the averaged isotropic DCS and squares are used for the pure m = 0, m = 0 component (summed over J = 0-4). A comparison between the isotropic case (circles) and the m-to-m component (squares) unveils once again a similar mechanism discussed above in terms of ICS: i.e., the collision-induced branching ratio associated to lower m-values, mainly m = 0, is likely to contribute the most but is often suppressed due to the typical average over all 2j + 1 m-terms. However, the dominant contribution from m = 0 is unlocked by choosing β = 0 (black curves). As the molecular alignment angle changes from β = 0 → 90 • , the relative contributions of m = 0 become smaller in |αυj; m exp = 0 yielding, therefore, an overall smaller cross section. The preference for (backward) forward scattering (θ → 0, θ → π) is seen to be mostly captured by the m = 0, m = 0 component and this preference is either magnified (β → 0) or suppressed (β → 90 • ) by a given choice of β. By comparing results depicted by circles and squares again, one notices that the small isotropic DCS value at θ = 90 • includes a negligible, if any, contribution from the m = 0, m = 0 component. Thus, the DCS is the smallest and almost independent of the scattering angle for β = 90 • .
For the sake of completeness, Fig. (12) illustrates a typical DCS for the H + D 2 (υ = 4, j = 2, m exp = 0) −→ D + HD(υ = 3, j = 0, m = 0) reaction. In this case, the DCS at 1 K is peaked near 90 • but fairly negligible overall. As the collision energy increase to 5 K and 10 K, the scattering becomes preferentially backward with maximum reactivity attained at 10 K, in which case, from Fig. (9), a = 3 character is predicted. At all three energies β ≤ 15 • favors the reaction.

IV. CONCLUSION AND PROSPECTS
In this work we report extensive quantum dynamics calculations of the H + D 2 and D + HD chemical reactions with an aim of controlling the reaction outcome through stereodynamic preparation of the molecular rotational state. Our approach extends the formalism presented in recent experimental studies of Perreault et al. 14,15,19,20 applied to rotational quenching of HD by H 2 , D 2 and He and D 2 by He to state-to-state chemistry for the first time. We report vast control through stereodynamic preparation of the molecular state, essentially allowing switching the reaction on or off through appropriate alignment of the molecular bond axis relative to the direction of the incident relative velocity vector.
The computations are performed using a state-of-the-art time-independent quantum reactive scattering formalism in hyperspherical coordinates at a moderate range of collision energies, 1-50 K. Both diatomic molecules were prepared in excited rovibrational states, υ = 4, j = 1-4, and several resonant features have been predicted to occur. Intensities of these resonances were found to be strongly sensitive to the stereodynamic preparation of the molecular rotational states. Several partial wave resonances were identified in this study with = 0-2 dominating at lower collision energies (below 10 K) and = 4-6 at higher energies. To our knowledge, resonances for the initial states considered here (υ = 4, j > 0, E coll = 0 − 50 K) are predicted for the first time and deserve additional investigations in order to characterize properties such as lifetime, quantum numbers etc. Our findings are significant in light of recent experimental progress in preparing vibrationally excited H 2 in the v = 7 vibrational level through the SARP techniques allowing the prospects of studying chemical reactions involving highly vibrationally excited molecules that are also prepared in a stereodynamic state 65 .
Our results reveal the following aspects which may or may not be universal to the specific system investigated here: (i) The isotropic reactive cross section, associated to freely rotating reactants, is composed of contributions from smaller m-values, in particular, m = 0 (m ≈ 0). However, as one averages over all 2j + 1 initial m states, these higher contributions are often washed out.
(ii) By imposing a certain spatial polarization on the rotational state j of reactants, a particular m-contribution (denoted by m exp in this work) can be filtered by rotating its internuclear axis.
(iii) The stereodynamically prepared state becomes a pure state as β → 0, i.e., ±m exp → ±m as β → 0, where β is the alignment angle of the molecular bond axis relative to the incident velocity vector, i.e. orientation of the SARP polarization lasers relative to the molecular beam in the experiments of Perreault et al.
(iv) The prepared state m exp becomes a broad mixture of all m-terms as β → 90 • .
(v) Due to (iii), the m = 0 contribution (largest for the present system) are accessible by a m exp = 0, β = 0 preparation, which generally enhances the collision-induced branching ratio of interest. This is equivalent to the H-SARP preparation in the SARP experiments.
(vi) Likewise, due to (iv), in order to enhance unfavored branching ratios, often associated to higher m (and m ) values, a preparation combining higher m exp with β → 90 • can be used.
(vii) A lower value of m exp combined with a higher value of β (and vice-versa) mostly suppresses the reaction.
Possibly, without loss of generality, points (i)-(vii) can be viewed as a set of propensity rules that might be exploited to effectively control any cold chemical reaction, as well as inelastic collisions.
At collision energies above a few Kelvin a diverse array of stereodynamical effects are predicted to occur due to the participation of various partial waves and different m-components of the rotationally excited states that allow additional control of the reaction outcome. Our results show that both the collision energy and the alignment angle β have a strong effect on the m → m branching ratios and the m-composition of a given m exp state, respectively. Likewise, we have found strong evidences suggesting that particular resonances can be selectively controlled, by means of orienting the reactant appropriately, due to their specific partial wave, m and m compositions.
Such a degree of control of a chemical reaction involving neutral systems, in a completely field-free environment, requires only ro-vibratational excitation, which is relatively easy to realize experimentally and it is also a universal property in contrast to, for instance, the need of working with electrically polar molecules. In particular, two-body losses in ultracold traps (E coll 1 K), due to chemical reactions and inelastic collisions, have been a long standing problem, mostly circumvent by means of an elegant use of external fields that prevents the colliding partners to approach one another [66][67][68] . Besides their effectiveness, most of these methods are not general and a given molecular system of interest must comply with specific electric properties. In contrast, we have demonstrated here how appropriate choices of m exp and βnamely, lower (higher) m exp and higher (lower) β -is mostly effective in "turning off" the chemical reaction in a very general fashion.
Finally, in terms of actual experimentation using H + D 2 and D + HD, our calculations suggest that choices of either j = 2 or j = 1 (υ = 4) are the best compromises in order to explore resonances at sufficiently low energies, where only few partial waves contribute. In particular, near 1 K, preparing either D + HD or H + D 2 with j = 1 is a promising choice. However, if the presence of higher partial waves are less of a problem experimentally, and higher collisional energies are of interest, then j = 2 and j = 3 (υ = 4) are better candidates as they are predicted to possess a series of overlapping resonant structures over a broader range of collision energies.

CONFLICTS OF INTEREST
There are no conflicts to declare.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available within the article.