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

Stereochemistry Aldegunde et al. (2006, 2011); Cernuto et al. (2018); Ascenzi et al. (2019); Falcinelli et al. (2019, 2021) along with coherent and optimal control of molecular processes Tannor, Kosloff, and Rice (1986); Rice and Zhao (2000); Levis, Menkir, and Rabitz (2001); Shapiro and Brumer (2012) 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 Cattaneo et al. (2018)– to an ultra–slow sub–Kelvin nuclear dynamics, in the context of either cold or ultracold molecular collisions Krems (2008); Balakrishnan (2016).

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₂ and HD + H₂) in the same beam stream enable cold collisions with relative kinetic energies of about 1 Kelvin Perreault, Mukherjee, and Zare (2017, 2018). The technique, combined with Stark–induced adiabatic Raman passage (SARP) Mukherjee and Zare (2011) 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 $\mathbf{j}$ combined with selection of specific $\mathbf{j}$–projections, say $m$, along a given quantization axis of choice Dong, Mukherjee, and Zare (2013); Mukherjee, Dong, and Zare (2014). As a consequence, the lower collision energies provided by the co–expansion can selectively discard higher and undesired Perreault, Mukherjee, and Zare (2017) external collisional orbital angular momentum states, $\ell$, 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₂ collisions with the latter revealing strong signatures of an $\ell=2$ partial wave resonance Perreault, Mukherjee, and Zare (2019); Zhou et al. (2021). 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. Theoretical prospects on that matter have been available for a while Aldegunde et al. (2005, 2006, 2011), and revisited recently in the context of coherent control Devolder, Tscherbul, and Brumer (2020).

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 de Miranda et al. (2011).

The experiments of Perreault et al. have been accompanied by the theoretical works of Croft et al. Croft et al. (2018); Croft and Balakrishnan (2019) on HD + H₂ and Morita et al. Morita and Balakrishnan (2020a, b) on He + HD collisions. Further, Jambrina et al. Jambrina et al. (2019) and Morita et al. Morita et al. (2020) have explored the control of isolated partial wave resonances through stereodynamic preparations of HD + H₂ and HCl + H₂ 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 Jambrina et al. (2021). In this work we extend these ideas to low–energy reactive collisions taking the prototypical chemical reactions of H + D${}_{2}\longleftrightarrow$ D + HD with polarization of the initial rotational state $\mathbf{j}$ of D₂ and HD. Thus, in what follows, we address questions such as (i) the role of bond axis orientation of D₂ (HD) relative to the initial relative velocity; (ii) whether it is possible to orient the reactants in a manner to enhance or suppress the production of HD (D₂); (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 Pack and Parker (1987); Kendrick et al. (1999); Kendrick (2003a) that has been used extensively to model atom–diatom reactive collisions, combined with the legacy BKMP2 potential energy surface (PES) for the H₃ system Boothroyd et al. (1996). 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. Kuppermann, Schatz, and Baer (1974) the H + H₂ 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. Goswami et al. (2020) who reported wavepacket calculations of the H + H₂ exchange reaction; Desrousseaux et al. Desrousseaux et al. (2018) 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. Zhou et al. (2020); and quasi–classical trajectories simulations by Bossion et al. Bossion et al. (2018, 2020). Kendrick and co–workers have reported extensive studies of geometric phase (GP) effects in cold and ultracold H + H₂, H + HD and D + HD collisions with vibrationally excited reactants that revealed a new interference mechanism that controls reactivity in the ultracold regime Hazra, Kendrick, and Balakrishnan (2016); Kendrick, Hazra, and Balakrishnan (2016a, b). Non–adiabatic effects in the H + H₂ and H + HD reactions with vibrationally excited H₂ and HD molecules for vibrational levels $\upsilon=$ 4–9 have also been reported by Kendrick Kendrick (2018a, b, 2019). Several theoretical and combined experiment–theory studies on H + D₂ are also available Wrede and Schnieder (1997); Banares et al. (1998); Kendrick et al. (1999); Fernández-Alonso, Bean, and Zare (1999); Kendrick (2000, 2001); Harich et al. (2002); Kendrick (2003b); Pomerantz et al. (2004).

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.

The methodology used herein is, in many aspects, very similar to that utilized in the previous works of Croft et al. Croft et al. (2018); Croft and Balakrishnan (2019) and Morita et al. Morita and Balakrishnan (2020a, b) except that reactive channels are also taken into consideration as a collision–induced outcome. Thus, upon a H + D${}_{2}\left(\upsilon,j\right)$ reactive collision (or likewise D + HD), the probability of detecting a scattered HD$\left(\upsilon^{\prime},j^{\prime}\right)$ product (or D₂) by a detector placed at $\Omega=(\theta,\phi)$, at a given collision energy $E_{\mathrm{coll}}$, is governed by the state–to–state scattering amplitude Kendrick (2003a) $f^{J}_{nn^{\prime}}\left(\theta,E_{\mathrm{coll}}\right)$, i.e.

where we have multiplied the original scattering amplitude expression of Kendrick, Eq. (22) Kendrick (2003a), by $\sqrt{k_{n\prime}/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 $\mathbf{J}$:

Also, the sums over $\ell$ and $\ell^{\prime}$ effectively average over the azimuthal angle $\phi$ so that Eq. (1) is independent of $\phi$. The initial and final states are denoted by $n=\left(\alpha,\upsilon,j,m\right)$ and $n^{\prime}=\left(\alpha^{\prime},\upsilon^{\prime},j^{\prime},m^{\prime}\right)$ with $\alpha$ denoting the different atom–diatom arrangements, i.e. $\alpha=1$ for H + D₂ and $\alpha=2,3$ for D + HD; $\upsilon 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 $\theta$ measures the deflection between the asymptotic velocity vectors of H (relative to D₂) and of D (relative to HD); $S_{n\ell n^{\prime}\ell^{\prime}}^{J}$ is the energy–dependent scattering matrix for a given total angular momentum, written in the $\ell$–representation $\left(\mathbf{J}=\mathbf{j}+\mathbf{\ell}\right)$, 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 $\mu$ 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_{\mathrm{exp}}$ of the molecular state $\left|\alpha\upsilon jm_{\mathrm{exp}}\right\rangle$, prepared experimentally, corresponds to a coherent superposition of all possible pure projections $m$, in $\left\{\left|\alpha\upsilon jm\right\rangle\right\}$, i.e.

where $\beta$ is the angle between the internuclear axis of D₂ (or HD) and the incident velocity vector. In the experiments of Perreault et al. $\beta$ denotes the orientation between the polarization lasers and the molecular beam axis; where choices of $\beta$ correspond to either a parallel orientation $\left(\beta=0\right)$, known as a H–SARP preparation, or a perpendicular orientation $\left(\beta=90^{\circ}\right)$, 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=\pm 1$, was also adopted in their recent work on He + HD and He + D₂ collisions Perreault, Mukherjee, and Zare (2019); Zhou et al. (2021). For exploring possible control scenarios, several values of $\beta$ within $\left[0,90^{\circ}\right]$ are considered below.

As discussed by Croft and Balakrishnan Croft and Balakrishnan (2019), since the experimental results correspond to an integral over the azimuthal angle $\phi$, one can express the differential cross sections (DCS) for any SARP–prepared state as a sum of contributions from each $m$–state weighted by $\left|d_{mm_{\mathrm{exp}}}^{j}\left(\beta\right)\right|^{2}$. For the state–to–state DCS, this yields

where, $\widetilde{n}=\left(\alpha,\upsilon,j,m_{\mathrm{exp}}\right)$. Likewise, for the integral cross section (ICS),

Figure (1) shows the $\beta$–dependence of the Wigner’s reduced rotation matrices for the hypothetical cases of $j=1,2,3$ along with their pair of $\left(m,m_{\mathrm{exp}}\right)$ values (negative values are not shown). Upon a simple inspection of Fig. (1), the following aspects are apparent: (i) For $\beta=0$, $d_{mm_{\mathrm{exp}}}^{j}=0$ whenever off–diagonal elements are concerned, i.e. $m_{\mathrm{exp}}\neq m$; (ii) similarly, $d_{mm_{\mathrm{exp}}}^{j}=1$, for diagonal terms at $\beta=0$ and, as a consequence, parallel orientations of the reactants are likely to be composed by terms driven by $m_{\mathrm{exp}}=m$, i.e. $\left|\alpha\upsilon jm_{\mathrm{exp}}\right\rangle=\left|\alpha\upsilon jm\right\rangle$; (iii) in the particular case of $m_{\mathrm{exp}}=0$ and $\beta=0$, due to (ii), the expansion of Eq. (4) collapses to a single and isotropic term in which $m_{\mathrm{exp}}=m=0$; (iv) $d_{00}^{0}\left(\beta\right)=1$ for any value of $\beta$ – not shown in Fig. (1) – and, thus, as expected, processes associated to rotationless reactants are invariant upon changes of $\beta$.

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 $\upsilon$ with a projection $m_{\mathrm{exp}}=0$, and subjected to two orientations of its internuclear axis relative to the incoming relative velocity vector: $\beta=0$ and $\beta=45^{\circ}$. The first case, $\beta=0$, has been addressed above and the resultant state is simply $\left|\alpha\upsilon j;m_{\mathrm{exp}}=0\right\rangle=\left|\alpha\upsilon j;m=0\right\rangle$. 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 $\beta=0$ for a $m_{\mathrm{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 $\beta=45^{\circ}$ alignment is made, the same $m_{\mathrm{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=\pm 1$ and $m=\pm 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_{mm_{\mathrm{exp}}}^{j}(\beta)$ matrix elements corresponding to the chosen $\beta$ values. If a $m_{\mathrm{exp}}=2$ preparation is considered at the same $\beta=45^{\circ}$ 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_{\mathrm{exp}}=0$ at the two typical orientations of $\beta=0$ and $\beta=90^{\circ}$ (V–SARP and H–SARP) has been discussed in great detail by Perreault et al. Perreault, Mukherjee, and Zare (2017, 2018).

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 Pack and Parker (1987); Kendrick et al. (1999); Kendrick (2000, 2001, 2003b). Details of the computational formalism have been discussed in great detail in prior works of Kendrick Kendrick et al. (1999); Kendrick (2000, 2001, 2003b) and in recent works of Kendrick et al. Kendrick, Hazra, and Balakrishnan (2016b) on the D + HD and H + D₂ reactions. We shall remark only on some key details here. The hyper–radius $\rho$ 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 Pack and Parker (1987). Smith–Johnson hyperspherical coordinates Smith (1962); Kendrick et al. (1999) are used in this region where collision–induced re–arrangements are more likely to occur, and Delves hyperspherical coordinates Delves (1958, 1960); Parker et al. (2002) 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 Kendrick et al. (1999); Kendrick (2000, 2001, 2003b). Briefly, the six–dimensional three–body problem is reduced to a set of one–dimensional coupled equations along the scattering coordinate $\rho$, 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 $\rho$. 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 Sorensen (1992); Maschhoff and Sorensen (1996). The corresponding eigenvalues within the Delves region are evaluated using a 1D Numerov propagator Johnson (1977). 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 Johnson (1973), first from $\rho_{\mathrm{min}}$ to $\rho_{\mathrm{match}}$. At $\rho_{\mathrm{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 $\rho_{\mathrm{match}}$ to $\rho_{\mathrm{max}}$, a sufficiently large value of $\rho$ where the interaction potential is negligible. At $\rho_{\mathrm{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_{n\ell n^{\prime}\ell^{\prime}}^{J}$ matrix Pack and Parker (1987). 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₂ 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 $\left(p,P,J\right)$ 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₂).

As both H + D₂ and D + HD chemical reactions (as well as other isotopic combinations) have been modeled with the APH3D code in the past Kendrick (2001), 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 Kendrick (2000, 2001, 2003b); Kendrick, Hazra, and Balakrishnan (2016a, b). Thus, the APH region ranges from $\rho_{\mathrm{min}}=1.6$ $a_{0}$ to $\rho_{\mathrm{match}}=7.8$ $a_{0}$ with 64 sectors varying logarithmically with a step size of $\triangle\rho_{\mathrm{aph}}=0.025$ $a_{0}$. The Delves region, starting at $\rho_{\mathrm{match}}$, extends to $\rho_{\mathrm{max}}=78$ $a_{0}$ with 350 sectors varying linearly with a step size $\triangle\rho_{\mathrm{delves}}=0.2$ $a_{0}$. Reactions such as H + D₂ are known to require a few dozens of partial waves whenever collision energies exceed 1 eV  Kendrick (2001, 2003b). 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_{\mathrm{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 $\rho_{\mathrm{min}}$ and $\rho_{\mathrm{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 Kendrick et al. (1999): (i) $l_{\mathrm{max}}=103$, $m_{\mathrm{max}}=190$ for $\rho=$ 1.6–2.9; (ii) $l_{\mathrm{max}}=115$, $m_{\mathrm{max}}=214$ for $\rho=$ 2.9–3.7; (iii) $l_{\mathrm{max}}=123$, $m_{\mathrm{max}}=232$ for $\rho=$ 3.7–4.6; (iv) $l_{\mathrm{max}}=135$, $m_{\mathrm{max}}=250$ for $\rho=$ 4.6–5.8; and, (v) $l_{\mathrm{max}}=143$, $m_{\mathrm{max}}=274$ for $\rho=$ 5.8–7.8. Finally, 350 $\left(J=0\right)$, 650 $\left(J=1\right)$, 950 $\left(J=2\right)$, 1250 $\left(J=3\right)$, 1550 $\left(J=4\right)$ and 1850 $\left(J=5\right)$ channels are propagated for 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 $\left(J=6,7\right)$, 2700 $\left(J=8,9\right)$ channels for even $P$; 1860 $\left(J=6\right)$, 2400 $\left(J=7,8\right)$ and 3000 $\left(J=9\right)$ channels for odd $P$.

In what follows the collision–induced exothermic reaction $\left(\alpha\neq\alpha^{\prime}\right)$ and inelastic de–excitation $\left(\alpha=\alpha^{\prime}\right)$ are studied for the reactants in the initial rovibrational state $\upsilon=4,j=$ 1, 2 $(J_{\mathrm{max}}=4)$, 3 $(J_{\mathrm{max}}=7)$ and 4 $(J_{\mathrm{max}}=9)$; whereas the internal states of the products are chosen to be either $\upsilon^{\prime}=\upsilon$ or $\upsilon^{\prime}=\upsilon-1$ as well as $j^{\prime}=$ 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. Hazra, Kendrick, and Balakrishnan (2016), the effective reaction pathway along the set of adiabatic coupled potentials, for hydrogen exchange reactions, becomes barrierless for the vibrational manifolds above $\upsilon=3$, which is the primary motivation underlying the choice of $\upsilon=4$ herein, since a natural boost in terms of reactivity is expected for low energy collisions Kendrick, Hazra, and Balakrishnan (2016b).

From an experimental perspective there is a growing effort towards addressing even higher vibrational states Perreault, Mukherjee, and Zare (2016); Mukherjee, Perreault, and Zare (2017) in order to explore the inherent delocalization of such weakly–bound systems Perreault et al. (2020). 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₂. 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 $\upsilon=4,j=0$ and below Kendrick (2000); Kendrick, Hazra, and Balakrishnan (2016a, b). 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. Croft et al. (2018); Croft and Balakrishnan (2019) and Morita et al.  Morita and Balakrishnan (2020a, b) 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 $\upsilon=4,j=1-4$ thresholds of D₂ 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 resonances. A detailed characterization in terms of the angular momentum partial waves, width, lifetime and bound state structure of the effective potential well will be addressed elsewhere. Here, we focus on the sensitivity of these resonances to stereodynamic preparation, in particular, to choices of both $m_{\mathrm{exp}}$ and $\beta$ in a manner to simulate physical conditions that resemble typical experimental conditions of Perreault et al. Despite strict choices of $m_{\mathrm{exp}}=0$, $\beta=0,\beta=90^{\circ}$ and $E_{\mathrm{coll}}\approx$ 1–5 K in the experiments reported so far, we take some (theoretical) freedom to explore other possibilities.

In this work we report extensive quantum dynamics calculations of the H + D₂ 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. Perreault, Mukherjee, and Zare (2017, 2018, 2019); Zhou et al. (2021) applied to rotational quenching of HD by H₂, D₂ and He and D₂ 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, $\upsilon=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 $\ell=$ 0–2 dominating at lower collision energies (below 10 K) and $\ell=$ 4–6 at higher energies. To our knowledge, resonances for the initial states considered here $\left(\upsilon=4,j>0,E_{\mathrm{coll}}=0-50\,\mathrm{K}\right)$ 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₂ 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 Perreault et al. (2020).

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\left(m^{\prime}\approx 0\right)$. 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 $\mathbf{j}$ of reactants, a particular $m$– contribution (denoted by $m_{\mathrm{exp}}$ in this work) can be filtered by rotating its internuclear axis.

(iii) The stereodynamically prepared state becomes a pure state as $\beta\rightarrow 0$, i.e., $\pm m_{\mathrm{exp}}\rightarrow\pm m$ as $\beta\rightarrow 0$, where $\beta$ 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_{\mathrm{exp}}$ becomes a broad mixture of all $m$–terms as $\beta\rightarrow 90^{\circ}$.

(v) Due to (iii), the $m=0$ contribution (largest for the present system) are accessible by a $m_{\mathrm{exp}}=0,\beta=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^{\prime}$) values, a preparation combining higher $m_{\mathrm{exp}}$ with $\beta\rightarrow 90^{\circ}$ can be used.

(vii) A lower value of $m_{\mathrm{exp}}$ combined with a higher value of $\beta$ (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 $\beta$ have a strong effect on the $m\rightarrow m^{\prime}$ branching ratios and the $m$–composition of a given $m_{\mathrm{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^{\prime}$ 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 $\left(E_{\mathrm{coll}}\ll 1\,\mathrm{K}\right)$, 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 Avdeenkov, Bortolotti, and Bohn (2006); Wang and Quéméner (2015); González-Martínez, Bohn, and Quéméner (2017). 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_{\mathrm{exp}}$ and $\beta$ – namely, lower (higher) $m_{\mathrm{exp}}$ and higher (lower) $\beta$ – is mostly effective in “turning off” the chemical reaction in a very general fashion.

Finally, in terms of actual experimentation using H + D₂ and D + HD, our calculations suggest that choices of either $j=2$ or $j=1\left(\upsilon=4\right)$ 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₂ 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\left(\upsilon=4\right)$ are better candidates as they are predicted to possess a series of overlapping resonant structures over a broader range of collision energies.

There are no conflicts to declare.

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