State-selective cross sections from ring polymer molecular dynamics

Understanding the influence of different forms of energy (eg, translational, vibrational, rotational) on chemical reactions is a key goal and great challenge in physical chemistry. Very recently, we proposed a new approach to obtain state-selective cross sections that approximately include quantum effects such as zero-point energy and tunneling. The method is a combination of the widely used quasiclassical trajectory approach (QCT) and the ring polymer molecular dynamics method and thus is numerically very efficient and easily employed. Here, we present a detailed description of the method and exhaustive tests of its accuracy and applicability. The robustness of the approach is tested, as well as the convergence with the number of beads. The approach is then applied to several prototypical X + H 2 ( ν = 0, 1), X = Mu, H, D, F, Cl reactions over a wide range of collision energies. Good agreement with rigorous quantum dynamics simulations is found for most cases. Encouraging improvement over QCT results is found for particular cases, while only a small increase in numerical cost is required. H and results for reactivity to the of effects contributing to the reactivity, the vibrationally adiabatic barrier. for the the results overestimates the reactivity for E tot 0.95 eV. RPMD vibrational excitation approach underestimates the QM scattering results for higher collision energies but approximates tunneling contributions to some extent the exact quantum ICSs better than the different QCT variants. that our model with RPMD approximates the tunneling contribution through the = 1 vibrationally adiabatic barrier well as the finite extension of the ring polymer stretches over the are found at higher collision energies, where the results coincide with cross sections obtained from QCT. QCT-RPMD for the F + H 2 reaction slightly improves over QCT simulations but cannot fully reproduce exact QM simulations. For low energies, this difference can partially be due to reactive resonance effects in this reaction, but overall, there are also some shortcomings in our vibrational excitation scheme at low collision energies, which needs improvement. Ideas to improve the current approach in this regard are under investigation. The method can be generalized to polyatomic systems, is numerically efficient and easily implemented, and shows encouraging improvement over QCT simulations. Several extensions to include electronically nonadiabatic effects and the resolution of final states are planned for future work.

the approach is tested, as well as the convergence with the number of beads. The approach is then applied to several prototypical X + H 2 (ν = 0, 1), X = Mu, H, D, F, Cl reactions over a wide range of collision energies. Good agreement with rigorous quantum dynamics simulations is found for most cases. Encouraging improvement over QCT results is found for particular cases, while only a small increase in numerical cost is required.
integral cross sections, nuclear quantum effects, ring polymer molecular dynamics, stateresolved reaction dynamics

| INTRODUCTION
The detailed study and understanding of reaction processes is a central challenge of chemical physics and theoretical chemistry. Chemisorption of a molecule onto a surface and bimolecular reactions in the gas phase are two important classes of reactions studied. Their detailed understanding is of relevance in many areas, for example, atmospheric and interstellar chemistry, combustion, and catalysis. [1][2][3][4][5][6][7][8][9][10] The most detailed results, initial state-selected and fully quantum state-resolved reaction probabilities and cross sections, can be computed from full-dimensional quantum dynamics simulations for reactions involving only few atoms. The biggest system treated today is the H + CH 4 !H 2 +CH 3 reaction. [11][12][13][14][15][16] To this end, the multiconfigurational time-dependent Hartree approach [17][18][19][20][21][22] was employed for the wavepacket propagation, and reaction probabilities were obtained using the quantum transition state concept. [23][24][25][26][27][28][29][30] However, these simulations are numerically very expensive and take many months of simulation time for each reaction. They also require fitted or interpolated potential energy surfaces (PESs) that can be evaluated efficiently and dedicated curvilinear coordinate systems with complicated kinetic energy operators. [31][32][33] To avoid the high numerical cost of rigorous, full-dimensional quantum scattering simulations, many reactive scattering calculations employ the quasiclassical trajectory (QCT) approach. [34][35][36][37][38][39] QCT treats all nuclei classically evolving on accurate PES or using gradients from on-the-fly ab initio calculations. [40] Initial conditions are chosen to mimic the quantum state under investigation. This is done by adding the correct amount of zero-point energy (ZPE) and energy for ro-vibrational excitations to the system. The QCT approach has been employed successfully to study many bimolecular and gas-surface reactions in detail. Reactions with many degrees of freedom, even including surface motion, as well as reactions with several heavy atoms, can be investigated and understood by employing QCT simulations, [41][42][43][44][45] where effects such as the depression of reactivity by collision energy were studied in detail. [46] The most accurate results are obtained for energies well above the threshold as tunneling is often not important for these cases. Due to the purely classical dynamics involved, the QCT approach is numerically very efficient, and large systems can be treated. Yet, due to the classical nature of the QCT simulations, quantum effects such as ZPE effects and tunneling are not systematically taken into account. Please note that there exist many advanced procedures for QCT, which can remedy those effects. [47][48][49] However, QCT simulations can be inaccurate for particular cases where quantum effects become important. For these cases, the inaccuracies can often be traced back to the following three points. First, the classical nature of the simulations allows any ZPE put into a certain degree of freedom initially to artificially leak out into any other degrees of freedom. If the energy leaks into the reaction coordinate, it can increase the reactivity in the QCT calculations. In certain cases, this may lead to an increase of the reactivity compared to rigorous quantum calculations. [50,51] Second, classical simulations do not respect ZPE constraints for the products, which can allow forbidden product outcomes and can lead to an overestimation of the reactivity. Third, the QCT simulations do not include tunneling effects. This can lead to an underestimation of the cross sections, particularly when the reaction includes a proton transfer and in the threshold region. Different problems can be manifest in a single reaction when looking at the reactivity of different initial vibrational states. [52] The approximate inclusion of nuclear quantum effects (NQEs) such as ZPE and tunneling into classical-like dynamics simulations has been a very active field of research. Various efficient methods along this way have emerged. One of the widely used approaches is the ring polymer molecular dynamics (RPMD) approach. [53][54][55][56][57][58][59] In RPMD, an initial quantum Boltzmann distribution is represented by an ensemble of discrete representations of the imaginary-time path integral, the so-called ring polymer. Each ring polymer consists of classical replicas (termed beads) of the system joined by harmonic springs. Real-time correlation functions are then obtained from classical dynamics of these ring polymers, that is, classical dynamics in an extended ring polymer phase space. RPMD has several appealing features. It is exact in the short-time, classical, and hightemperature limits. [60,53] Exact correlation functions can also be obtained for harmonic potentials. ZPE effects are incorporated through the harmonic springs connecting the beads, and RPMD therefore respects ZPE constraints and does not suffer from ZPE leakage. [61][62][63] Tunneling effects are also incorporated due to the finite extension of the ring polymer, which can stretch over a barrier and therefore lower its effective height. [54,55,[64][65][66] Originally, RPMD was formulated for the calculation of real-time correlation functions associated with equilibrium initial conditions. This allowed the efficient calculation of thermal rate constants, [54,55,[64][65][66][67][68][69][70][71][72]59] diffusion coefficients, [56,57,73,74] and-to some extentvibrational spectra. [75][76][77][78] Recently, an approach to obtain quantities associated with different nonequilibrium initial conditions, for example, a momentum kick or a vertical excitation, within the RPMD framework was proposed. [79] It was shown that RPMD can be obtained from the more general Matsubara dynamics framework. [80,81] Despite its success, RPMD also suffers from several drawbacks. Due to the classical nature of the dynamics simulations, it cannot describe real-time coherences. Furthermore, when calculating vibrational spectra, it suffers from the spuriousmode effect. [75][76][77] Very recently, the nonequilibrium initial conditions have been used to obtain microcanonical rate constants from RPMD simulations; however, initial state selectivity could not be addressed in these works. [82,62] We have proposed an efficient approach that combines the RPMD approach with QCT simulations to obtain initial state-selective cross sections for bimolecular reactions. [83] Initial state-selective cross sections for the Mu/H/D + H 2 (ν = 0,1) reactions were calculated close to the threshold energy. It was shown that the approach approximately incorporates ZPE constraints and includes some tunneling contributions. In this manuscript, we explore the approach in more detail. In particular, we test the robustness of our approach with respect to the choice of the ring polymer parameters (n, β) and the vibrational excitation scheme. Integral cross sections (ICSs) for the five X + H 2 , X = Mu, H, D, F, Cl reactions are studied for H 2 in its ground state, as well as vibrationally excited state, and for a wide range of collision energies.
2 | THEORY 2.1 | Quasiclassical trajectory method QCT has been widely reviewed, and thus, we will only give a brief overview of the approach here. [9,10,35] We focus on general triatomic systems A + BC(ν, j) and mostly follow Karplus et al. [34] The molecule BC is initially in a ro-vibrational state with vibrational and rotational quantum numbers ν and j, respectively. Quasiclassical initial positions and momenta mimicking state (ν, j) are obtained as follows: The center of mass of BC is set at the origin of the coordinate system. The intramolecular distance r = j rj, where r is the vector between B and C, is set to either of the classical turning-point values r + or r − , corresponding to the energy of the ro-vibrational state, E ν, j . The total internal molecular momentum P = μ BC (v C − v B ), where μ BC is the reduced mass of BC, is set so that P 2 = j j+ 1 ð Þ=r 2 AE and P Á r = 0. BC is arbitrarily oriented along r × κ, where κ is the unit vector along the x-axis. A vibrational phase is added by propagating BC for a random fraction of its vibrational period, and finally, BC is randomly oriented. The procedure described above is restricted to diatomic molecules. For polyatomic molecules, one can initialize the vibrational state via normal-mode sampling. [84,85] Initial position and momenta of each mass-weighted normal mode are chosen as where ν i is the quantum vibrational number corresponding to the ith mode, and ϕ i [0, 2π] is a random phase. One can then adiabatically switch the system to the full potential for better vibrational state accuracy. [84] Cartesian coordinates are then obtained by inverse normal mode transformation.
The atom A is initially placed at a fixed distance from the center of mass of BC. The distance is chosen such that there is no interaction between A and BC. The impact parameter b is chosen so that b 2 =b 2 max is uniformly sampled in [0, 1], where b max is the maximum impact parameter so that no reactions occur for any b > b max . To initiate the reaction, the relative momentum between A and BC is set as P rel = ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2μE col p κ, where E col and μ are the collision energy and the reduced mass of the system, respectively. The system is then propagated by solving Newton's equations of motion.
ICSs for a given total energy E tot are obtained as where N and N R are the number of total and reactive trajectories, respectively. [34] The number of reactive trajectories is obtained by evaluating suitable distances criteria that distinguish the products.

| Ring polymer molecular dynamics
RPMD is based on the ring polymer extended phase space used to represent imaginary-time path integrals. For a quantum mechanical system, the canonical partition function Z at inverse temperature β = 1/k B T can be computed as with n the number of replicas in the ring polymer, β n = β/n, and the ring polymer Hamiltonian H RP n reading as where N is the number of atoms, q k i and p k i the position and momentum vectors of the kth bead for the ith atom, and V the molecular potential interaction. [86] The Boltzmann distribution in the ring polymer phase space can be effectively sampled using molecular dynamics.
The quantum mechanical Kubo-transformed equilibrium real-time correlation function withÂ andB arbitrary operators, is approximated within RPMD as C with (q t , p t ) propagated with H RP n and A n q, p ð Þ= 1 n where A(q i , p i ) is the classical function corresponding toÂ evaluated at the ith bead coordinates. RPMD is exact in the short-time, high-temperature, and harmonic oscillator limits. In addition, it has been shown that RPMD can be derived from Matsubara dynamics. [80,81] Recently, it has been shown that RPMD can be used to approximate nonequilibrium correlation functions and expectation values, [79] for example, for cases with an initial momentum impulse or "kick". The nonequilibrium expectation valuê is approximated as with (q t , p t ) propagated with H RP 1,n .

| Initialization of BC
The initialization of BC mimicking the initial quantum state is performed in two steps. Please note that, for the initialization, we will employ a harmonic approximation for the PES of isolated BC where ω and Q are the harmonic frequency and corresponding mass-weighted normal-mode coordinate of BC. We will discuss the approach for a diatomic molecule here. The generalization to polyatomic molecules is discussed further in the following sections.
In the first step, initial ring polymer configurations for the mass-weighted normal mode are obtained by sampling the well-known thermal distribution for a harmonic potential. Cartesian coordinates for the beads are obtained by inverse normal-mode transformation. Sampling the normal mode guarantees that the molecule will not undergo spurious rotations or center-of-mass motion and that it will have the correct vibrational ZPE.
Yet, due to sampling the initial Boltzmann distribution, BC will have some additional thermal energy contributions. The thermal contributions will be minimized by choosing a low initial temperature, that is, high β value, for the RPMD simulations. The detailed choice of β will be discussed later in the subsection "Choice of β." By employing the path integral approach, we obtain correct mass-weighted initial quantum fluctuations for the position and momentum. For each bead k [1, n], we obtain In the second step, to mimic a vibrationally excited state, we modify the initial position and momentum of the ring polymers' centroids Q, P to add the correct amount of vibrational energy. This step is inspired by the QCT approach described in Section 2.1. Please note that the initial sampling already contains the correct amount of ZPE, and we only need to add additional vibrational energy. In this step, we employ mass-weighted normal modes (of the system) and ring polymer normal modes:Q k ,P k , k [0, n − 1]. To the position and momentum centroids (k = 0), we add the following quantities: with ϕ [0, 2π] a random phase. The resulting effect can be analyzed by employing the primitive path integral energy estimator ℰ 0 , which reads where N = 2 is the number of atoms in BC. As we employ a harmonic potential for the initialization, we can diagonalize the estimator using massweighted normal modes (of the system) and ring polymer normal modes. The estimator then reads Adding the two quantities given in Equation (13) to the path integral centroids, the path integral energy estimation reads Separating the centroids' contributions, we obtain and when averaging over the initial distribution, we obtain where the second term on the right-hand side vanishes asQ The mass-weighted position and momentum fluctuations for each bead k [1, n] are The initial ring polymer positions and velocities of BC are then used in the scattering calculations. To this end, one can directly switch to the full PES, which can depend on the anharmonicity of the full PES or one can employ an adiabatic switching procedure to obtain more accurate initial configurations. [84] However, for the direct switching, in some cases, the intramolecular distance of BC will be too small at the time of the switch due to the anharmonicity of the full PES, thus providing nonnegligible extravibrational energy, which can influence the reactivity at very low collision energies. We discard those trajectories where the centroid distance for BC falls below d ν = 1 min , corresponding to the lower initial length of excited BC in QCT at the time of the switch. This ad hoc discarding procedure is performed to limit the inaccuracies from the direct switching approximation but is not an intrinsic limitation of the overall approach as discussed later. The center-of-mass motion of BC is set following the QCT approach via a "kick" to the velocity centroids which induces a center-of-mass motion without perturbing the internal dynamics of BC. The centroid of the center of mass of BC is initially set at the origin of the coordinate system. For the initialization of a polyatomic molecule, the generalization is straightforward, and we sample in the same way each normal mode independently and set its centroid center of mass at rest and its angular momentum to zero using a standard iterative modification to the velocities routinely employed in QCT simulations to each ring polymer centroid of the molecule to compute the angular velocity Ω with Q and V the centroid position and centroid velocity vector, respectively. Direct switching is then performed. Employing adiabatic switching for polyatomic systems will be the focus of a forthcoming study.

| Initialization of A
A is initialized far away from the position of BC, with its noncentroid ring polymer normal modes sampled from the free ring polymer distribution along the direction of the initial QCT momentum kick κ. Then, the initial centroid position of A is set similar to QCT with an extra increment to the initial distance between the reactants due to the finite spatial extension of the ring polymers. The initial distance between the centroid of A and the centroid of the center of mass of BC d(A, BC) is set to where d(A, BC) QCT is the distance choice for QCT simulations, and R g A ,R g B ,R g C are the radii of gyration for the ring polymer of A, B, and C, respectively. Finally, the momentum centroid is then "kicked" following the QCT approach. The initial centroid velocity is set to

| Reactions and cross sections
After initializing both reaction partners, we propagate the system by solving the classical equations of motions in the extended ring polymer phase space. Considering the distances between the centroids of the ring polymers representing the different atoms, the trajectories are distinguished between reactive and nonreactive. Cross sections are then obtained following the QCT approach (see Equation 2) as where N and N R are the number of total and reactive trajectories, respectively, and E tot = E col + E reac and E reac is the total internal energy of the reactants. In this work, we focus on X + H 2 reactions and vibrational excitation and always set the angular momentum of H 2 , j to zero and employ the harmonic approximation. Therefore, we have E reac = E H2 = ω H2 ν + 0:5 ð Þ . This constitutes a reliable and consistent approximation in the case of a diatomic molecule, as we shall see later with the results for classical-like systems such as D + H 2 (ν = 0, 1).

| Choice of β
We touched on the importance of β for the current approach above. Here, we motivate the choice of β for the present approach. The decision is guided by four main considerations.
First, as the choice of β will influence the additional thermal energy in the initialization of the molecule, it should therefore be high enough to limit those thermal contributions. As the system is expected to behave classically when the collision energy becomes large, β should reach a lower limit, β − . Second, the dynamics of the ring polymer will be influenced by the choice of β through the spring constants, and therefore, β should be set such that the energy and momentum of each reactant reflect the temperature associated with β.
Third, β should be inversely proportional to the additional energy we give to the system in the initialization process by kicking the centroids, which comprise initiating the collision process and adding vibrational excitation.
Fourth, when ν = 0, we assume that β lies in the interval [β A , β BC ] for m BC > m A (and inversely otherwise) with where β X is obtained via equating j V X j with a typical thermal speed corresponding to different values of α (α = 4 π for the mean of the magnitude, α = 3 2 for the root mean square, and α = 1 for the most probable magnitude) at inverse temperature β X for a mass m X . To obtain a common β for all the ring polymers, we choose β as an average of both values, following a symmetry argument as, in the center-of-mass frame, we have jp A j = j p BC j.
So far, our choice of β for a given initial collision energy and vibrational state reads where α 4 π , 3 2 ,1 È É . This is the ansatz we used in our previous work (with α = 4 π and β − = 300). [83] Inducing a lower cut-off for β is necessary for high collision energies. In our previous work, β − was chosen to a fixed value so that the extra thermal energy contributions remained within 1% of the ZPE of H 2 . This can result in an abrupt cut-off. Nevertheless, we can avoid this abrupt cut-off by imposing the thermal contribution to the internal energy of the H 2 vibrational energy to be within a few percentages of the total energy instead of H 2 's vibrational energy. This way, β will smoothly decrease with the total energy instead of meeting a fixed cut-off value. Doing so slightly improves the results for the H/D+H 2 (ν = 1; j = 0) ICSs around the threshold energy values. We limit the thermal contributions below p% of the total energy, and thus: where E tot = E col + E H2 and E H2 is the total internal energy of H 2 . In this work, we focus on vibrational excitation and always set the angular momentum of H 2 to zero and therefore have E H2 = ω H2 ν + 0:5 ð Þ .
This choice can lead to high values of β at very low collision energies for ν = 0, which is only desired for cases where there is a large increase of ZPE from reactants to products, that is, for cases where there is a large decrease in "non-ZPE energy" for a given total energy, for example, for reactions like Mu + H 2 . In these cases, one requires a large value of β to describe the product side of the reaction well. However, for cases where the ZPE does not change much or where it decreases, the values of β employed are too large. We therefore correct our ansatz by only keeping the terms 2 + mBC mA from mtot μ for such systems when m A > m BC and ν = 0. This modification is a first attempt at obtaining a reliable and systematic choice of β for our approach and will require further study and, possibly, adjustments. For vibrationally excited systems, the change of ZPE is typically small compared to the vibrational energy added to the system, and thus, we do not employ this modification. All things considered, we have: The β profiles for the systems studies in this manuscript are shown in Figure 1. The range of collision energies employed in this work is between 0.08 and 1.55 eV, and the centroid kick (see Equation 13 with ν = 1) vibrational excitation energy of H 2 is 0.55 eV.

| System details
We will use our method to compute initial state-resolved ICSs for the following triatomic systems: D + H 2 (ν = 0, 1), H + H 2 (ν = 0, 1), Mu + H 2 (ν = 0, 1), Cl + H 2 (ν = 0, 1), and F + H 2 (ν = 0, 1), with the mass of Mu mass being 0.113 m H . In all cases, the angular momentum of H 2 is set to zero, and β is chosen according to Equation (28) (equivalently, Equation 26 for the systems Mu/H/D + H 2 ) with α = 4 π . The employed harmonic frequency of H 2 (ω H2 ) is 4407 cm −1 . For H 2 (BKMP2 PES), the difference between the exact energies for the vibrational ground and excited states and their harmonic counterparts are around 20 and 250 cm −1 , respectively. Thus, only a very small error is introduced for the ground state, while for the excited state, this error might play a small role around the threshold for the systems Cl/F + H 2 (ν = 1), which will be discussed in the Results section. Unless explicitly specified otherwise, the number of beads n employed throughout this work is expressed as an integer multiple of the dimensionless quantity βω H2 AE Ç , where dxe gives the least integer greater or equal to x. dβωe is often used in PIMD or RPMD simulations as an indicator of the minimum value for n such that the approximation in Equation (3) can be justified. We set p = 0.01 for β − so that the thermal contributions never exceed 1% of the total energy. For our simulations, the lowest value of β is reached for ν = 0 at E tot = 1.4 eV, where β = β − = 181. An upper limit for β, β max = 1500, is introduced for computational convenience. As for d ν = 1 min , the lower intermolecular distance corresponding to a vibrationally excited H 2 , it is set to 1.078 bohr.
To study the D/H/Mu+H 2 systems, we employ the BKMP2 PES. [87] The minimal initial distance between the reactants is 10 bohr. The maximal impact parameter in the QCT-RPMD simulations is typically about a factor of 1.3 for Mu+H 2 and 1.1 for D/H+H 2 larger than in the respective QCT simulations at the lowest collision energies when there is still reactivity in both cases. For each collision energy, the ICSs are computed by running 40 000 trajectories per collision energy using a modified Velocity-Verlet integrator with a time step of 0.02 fs. To study the systems Cl + H 2 and F + H 2 , we employ the Capecchi-Werner PES and the Li-Werner-Alexander-Lique PES, respectively. [88,89] As for the maximal impact parameters, in the QCT-RPMD simulations, it is noticeably higher than QCT for Cl + H 2 (ν = 1) by a factor of 1.5 at most and typically a factor of 1.15 for F + H 2 (ν = 0). The highest difference is found for the reaction F + H 2 (ν = 1), where QCT-RPMD trajectories at E col = 0.3 eV can react for impact parameters as high as 8 bohr, corresponding to a factor of 3 compared to QCT. The initial distance between the reactants is set to 10 and 12 bohr for Cl/F + H 2 , respectively. The same number of trajectories are run per collision energy with the same integrator but, this time, with a time step of 0.05 fs. In Figure 1, we plot the profiles of our choice of β for each system and vibrational state of H 2 over the total energy.

| RESULTS
In this section, we will discuss the results of the proposed approach for five different triatomic reactions: Mu/H/D + H 2 (ν = 0, 1) and Cl/F + H 2 (ν = 0, 1). Yet, before we go into the details of the specific reactions, we will investigate the robustness of our vibrational excitation scheme, the choice of β, and the convergence with the number of beads. To this end, we employ the Mu + H 2 (ν = 0, 1) reaction as a test case as this reaction shows the biggest difference between exact quantum dynamics and QCT.

| Convergence studies
First, let us consider the robustness of the proposed scheme to mimic the vibrational states ν = 0, 1 of H 2 . To this end, we investigate the time evolution of the primitive internal energy estimator averaged over the ring polymer configurations, initialized as described in the "Initialization of F I G U R E 1 β values for each reactive system treated BC" subsection and then propagated on the full PES. In this case, the second reactant is placed very far away from H 2 such that there is no interaction during the full simulation. The primitive internal energy estimators for H 2 reads where h…i is an average over the initial conditions. Here, we use 50 000 different initial conditions. We consider the primitive estimator regardless of the slow convergence as both virial variations are not suited for the ν = 1 initialization due to the absence of a suitable kinetic term to take into account the centroid stretch in Equation (15). The time evolution E RP ν q, p ð Þ t ð Þ is shown in Figure 2. For the ground state ν = 0 (blue and green lines), we observe almost stable average energy upon switching to the full PES at t = 0. As for the excited state ν = 1, yellow lines, we observe a small dependence on β. For β = 200, we find a stable average energy around the harmonic value. For β = 300, we find small jumps to lower energies at t = 0. This difference is due to the coupling of the beads with the now much more energetic centroids and the abrupt nature of the potential switch, for which the direct use of the same estimator before and after the instantaneous switch makes the evaluation of the energy approximate.
We also observe a very small energy damping from the centroids of around 3 cm −1 every 10 fs, which is too slow to interfere with the reaction process in the time frame of our simulations and is typically less than 250 fs. The damping also indicates that the leakage of energy from the centroid to the beads due to the anharmonicity of the potential is very small.
Second, we study the influence of β on the calculated cross sections. In Figure 3A, we display the cross sections for Mu + H 2 (ν = 0) in the vicinity the exact threshold energy for a fixed value of β. For each β, the ICS curve behaves quite differently, which indicates, as expected, the necessity to select β depending on the total energy of the system. The fact that the ICS calculated with β = 200 does not coincide with the one from QCT-histogram binning (HB) calculations is a confirmation that some nonclassical effects are still encompassed for this case, even at higher collision energies. As we use higher β, the ICS curves match the quantum threshold better for E tot < 0.9 eV and then tend to underestimate more and more the reactivity for higher total energy. This indicates that higher β values better capture the low energy reactivity and inversely so. A similar test for the system Mu + H 2 (ν = 1) is shown in Figure 3B. The same qualitative behavior of the ICS with respect to β is observed.
Third, we test the convergence of the scattering results with the number of beads employed. Figure 4A π ) for different number of beads. We observe graphical convergence with n ≥ 3n r resulting in fewer than 100 beads for all cases considered here (see Table 1). The biggest differences of about 23% can be observed around the threshold at E col ≈ 1.05 eV for ν = 0 and also around the threshold for the ν = 1 case, with a difference of around 20% at E col ≈ 0.2 eV. The presence of convergence with the number of beads shows that the method is consistent thorough its different steps.
F I G U R E 2 Energy estimators hE RP i for H 2 (ν = 0, 1) as a function of time. Reference values for the harmonic regime and exact energy levels of H 2 (ν = 0, 1) are indicated as solid and dashed lines, respectively

| D/H/Mu + H 2
These systems were studied in our previous work, which dealt with the ICSs for the corresponding systems around their thresholds. [83] Here, we recomputed the ICSs with a different β, that is, with the previous smooth energy-dependent cut-off instead of the abrupt β − = 300 cut-off, as well as a wider range of collision energies. First, we focus on the reactivity of H/D with H 2 in its vibrational ground state. The ICSs for the two reactions are shown in Figure 5. For those systems, β reaches its cut-off values β − (E tot ) around the thresholds, that is, for total energies of 0.62 eV F I G U R E 3 A, Mu + H 2 (ν = 0; j = 0) and B, Mu + H 2 (ν = 1; j = 0) QCT-RPMD ICSs for several fixed values of β with n = 4βω ωH 2 , where βω ωH 2 refers to the least integer i such that i ≥ βω ωH 2 . The corresponding numbers of beads used for β = 200, 400, 600, 700 are 16, 32, 48, and 56, respectively. ICSs, integral cross sections; QCT, quasiclassical trajectory approach; RPMD, ring polymer molecular dynamics method F I G U R E 4 A, Mu + H 2 (ν = 0; j = 0) and B, Mu + H 2 (ν = 1; j = 0) QCT-RPMD ICSs for several fixed values of β and β chosen according to Equation (26) with α = 4 π with varying number of beads given in Table 1, n r = βω H2 . ICSs, integral cross sections; QCT, quasiclassical trajectory approach; RPMD, ring polymer molecular dynamics method (threshold at 0.57 eV) and 0.55 eV (threshold at 0.58 eV) for H and D, respectively. As in our previous paper, the QCT-RPMD ICSs match closely with the quantum ICSs around the threshold, and we now observed the same match for high collision energies up to 1.4 eV. The almost classical behavior of the ICSs for all collision energies above the threshold is consistent with β reaching its lower cut-off value close to the threshold. However, we like to point out that the Gaussian binning (GB) approach, which typically increases the accuracy of QCT simulations (see Mu + H 2 below), decreases the accuracy here. QCT-RPMD, on the other hand, does not have a need for specifically designed binning schemes. The good agreement at higher collision energies justifies the proposed cut-off scheme for β at high collision energies.
For the calculation of ICSs for H/D with H 2 in its vibrational excited state, we have β = β − for all collision energies considered, which is consistent with QCT-HB ICSs and QM ICSs being very similar for those energies. The ICSs are shown in Figure 6. In our previous work, we studied this reaction up to a total energy of 1.1 eV and observed good agreements with QM results. We now match for total energies up to 1.4 eV. As in our previous work, we see that RPMD accurately reproduces the threshold energy. Only a slight underestimation for the H + H 2 (ν = 1) reaction is found. The ICSs at higher collision energies are very close to the QM results. This indicates that the proposed ring polymer vibrational excitation scheme for H 2 (ν = 1) works acceptably well from the threshold energies to reasonably high energies.
We switch our focus to the reactivity of Mu with H 2 (ν = 0, 1). For ν = 0, quantum effects govern the reactivity for the full energy range considered here. For ν = 1, quantum and QCT results start coming close at E tot = 2.2 eV. We studied these two reactions in our previous paper for T A B L E 1 The values of β (in Hartree −1 ) and the number of beads, n, employed in the Mu + H 2 (ν = 0,1) simulations Note: β is chosen according to Equation (26), and the number of beads is set relative to β and the vibrational frequency of H 2 , ω H2 to allow for systematic convergence tests (see Figure 4 and Section 3.2).
F I G U R E 5 A, D + H 2 (ν = 0; j = 0) and B, H + H 2 (ν = 0; j = 0) QCT-RPMD ICSs for n = 6βω H2 . ICSs, integral cross sections; QCT, quasiclassical trajectory approach; RPMD, ring polymer molecular dynamics method total energies close to the threshold. In this study, we extend the energy range to see how RPMD performs at approximating NQE for higher collision energies. The corresponding ICSs are displayed in Figure 7. For the reaction of H 2 in its ground state ( Figure 7A), we restate that the QCT-RPMD scattering results very well reproduce the exact QM ICSs calculations around the threshold. The QCT calculations using HB, that is, treating every trajectory equally, show a much lower reactivity threshold compared to the QM calculations due to the absence of ZPE constraints. [90][91][92] This problem can be partly avoided by making use of GB as can be seen in Figure 7. [93] GB lowers the contributions to the reactivity for trajectories, which violates the ZPE constraint when considering their products' classical energies, that is, more precisely, for which the products energies are far from the quantum vibrational energies. QCT with GB well predicts the threshold energy but underestimates the cross section away from the threshold, and for higher energies, it starts to overestimate the ICSs. For higher collision energies away from the threshold, we note that the RPMD results tend to underestimate the exact ICS. This is mostly noticeable above a total energy of 1.1 eV. In Figure 7B, we show the RPMD results for the same reaction but with H 2 vibrationally excited, ν = 1. Both QCT variants (HB and GB) cannot reproduce the exact QM results and display a higher threshold for reactivity due to the absence of quantum effects contributing to the reactivity, that is, tunneling through the vibrationally adiabatic barrier. [92] As mentioned in our previous paper, QCT-RPMD finds a similar threshold for the reaction as the QM results but slightly overestimates the reactivity for E tot < 0.95 eV. Our RPMD vibrational excitation approach slightly underestimates the QM scattering results for higher collision energies but approximates tunneling contributions to some extent and describes the exact quantum ICSs better than the different QCT variants. This shows that our model with RPMD approximates the tunneling contribution through the ν = 1 vibrationally adiabatic barrier well as the finite extension of the ring polymer stretches over the barrier and reduces the effective barrier height.
Only a slight increase in computational cost is required, and fewer than 60 beads are needed for convergence.
In order to study how the spatial extension of Mu accounts for the ZPE effects, we ran the same simulations using two different values of β for each reactant, that is, β = 10 au for Mu (shrinking the gyration radius of Mu to around 0.1 bohr) and β = 200,400 au for H 2 (ν = 0) and H 2 (ν = 1), respectively. A comparison of these results with QCT-HB is shown in Figure 8. For H 2 (ν = 0), the QCT-HB results are recovered with only a slight underestimation of the ICSs past E tot = 1.1 eV. This emphasizes the role of the spatial extension of the Mu ring polymer to describe the ZPE effects for the energy range considered. A comparison with the β = 200 QCT-RPMD ICS curve from Figure 4 shows that the sampling of ring polymer H 2 at high enough β to avoid thermal contributions to the energy is essential at low collision energies. For H 2 (ν = 1), we observe that a "point-like" Mu decreases the reactivity, but the RPMD results still show a lower threshold than the QCT results. The spatial extension of Mu only partly accounts for the NQE important for this reaction.

| Cl + H 2
We now apply our method to study the reactivity of Cl with H 2 in its vibrational ground and first excited states. This reaction involves vibrational nonadiabaticity, and the atom A is now much more massive. [94] This is the first time that the vibrational effects for Cl + H 2 are studied using RPMD. This makes it a good benchmark system to further test the applicability and the robustness of our method. Electronically, nonadiabatic F I G U R E 6 A, D + H 2 (ν = 1; j = 0) and B, H + H 2 (ν = 1; j = 0) QCT-RPMD ICSs for n = 6βω H2 . ICSs, integral cross sections; QCT, quasiclassical trajectory approach; RPMD, ring polymer molecular dynamics method effects and spin-orbit couplings are not considered. The results for the RPMD ICSs, compared with QCT-HB and the quantum dynamics results, are displayed in Figure 9. For Cl + H 2 (ν = 0), we observe a good match with the QM results for all collision energies considered, which is expected as QCT-HB only slightly underestimates the reactivity. Yet, we find a small improvement over the QCT results for all energies considered.
As for Cl + H 2 (ν = 1), which involves vibrational nonadiabaticity effects that help overcoming the barrier, we observe that RPMD correctly predicts the threshold, contrary to QCT, which-just as with Mu + H 2 (ν = 1)-overestimates the threshold. [94] RPMD reproduces the threshold correctly and predicts the exact quantum ICSs for E tot < 0.95 eV well, while QCT-HB underestimates the ICSs for all energies considered. In this regime, RPMD describes the contribution of vibrational excitation to the reactivity correctly, while QCT seems to be, to a certain extent, vibrationally adiabatic. [94] This correct prediction of the threshold behavior by our approach has to be partially tempered and put into perspective considering the presence of the direct switching approximation, which can potentially lead to a slight spurious shift of the threshold predicted by RPMD. Nevertheless, the RPMD results would remain an improvement over the QCT predictions. Between 0.95 and 1.10 eV for the total energy, we observe that RPMD results start to underestimate the QM ICSs and tends toward the QCT-HB results. For higher total energies (≥ 1.1 eV), QCT-RPMD ICSs coincide with QCT ICSs, which hints that, with the high masses involved here, RPMD becomes more classical at higher collision energies contrary to Mu + H 2 (ν = 0, 1), where RPMD approximated the quantum effects very well, even at similar high energies.

| F + H 2
We now test our method with the important and notoriously intricate F + H 2 (ν = 0, 1) benchmark reactions. The F + H 2 system has a very low barrier, a bent transition state geometry, and is governed by quantum effects at low energies such as tunneling and resonances. [95,96] The ICSs increase very quickly with collision energy, and the considerable shift in the threshold between F + H 2 (ν = 0) and F + H 2 (ν = 1) indicates vibrational adiabaticity and a low vibrational enhancement. The resulting ICSs are displayed in Figure 10.
Overall agreement between QCT-RPMD and QM simulations is not as good for this case. For ν = 0 and for total energies below 0.65 eV, QCT-RPMD predicts slightly higher reactivity than QCT but underestimates the reactivity from QM simulations. It was partially expected as RPMD is not suited to describe quantum dynamical resonance effects, which are important in this system and are known to enhance the tunneling at very low collision energies for this reaction. [89,95] For the ν = 1 case, QCT-RPMD also underestimates the reactivity but does a bit better than QCT-HB around the threshold, partially capturing the enhancing quantum effects near the threshold. It was shown, contrary to the ν = 0 F I G U R E 9 A, Cl + H 2 (ν = 0; j = 0) and B, Cl + H 2 (ν = 1; j = 0) QCT-RPMD ICSs for n = 6βω H2 . ICSs, integral cross sections; QCT, quasiclassical trajectory approach; RPMD, ring polymer molecular dynamics method F I G U R E 1 0 A, F + H 2 (ν = 0; j = 0) and B, F + H 2 (ν = 1; j = 0) QCT-RPMD ICSs for n = 6βω H2 and β max = 1500. ICSs, integral cross sections; QCT, quasiclassical trajectory approach; RPMD, ring polymer molecular dynamics method case, that no reactive resonance occurs in the F + H 2 (ν = 1, j = 0) reaction, with this work instead indicating that the reaction mainly occurs via direct abstraction with a very strong dependence on the nature of the quantum H 2 -vibration. [97] Thus, it is a bit surprising that the QCT-RPMD approach underestimates the reactivity considerably. However, as the F + H 2 (ν = 1) reaction is nearly barrierless, it proceeds with minimal collision energy (eg, E col ≈ 0.003 eV at the threshold). This low collision energy leads to significantly increased propagation time, and therefore, the ring polymer H 2 vibrational energy undergoes a damping, which can lead to a loss of about 5% (0.04 eV) of its energy. This loss is nonnegligible for this reaction as the reactivity drastically increases low collision energies. The presence of the possible spurious shift of total energy due to the direct switching approximation also needs to be acknowledged when considering the threshold predictions of our approach. Those considerations constitute a limit of our current approach for very low collision energies and can contribute to QCT-RPMD underestimating the reactivity in this particular case. Extensions of the present approach to reduce those effects need to be investigated.

| DISCUSSION
So far, our method describes well the exact quantum ICSs around the threshold and medium collision energies for the reactive systems Cl/D/H/ Mu + H 2 (ν = 0, 1). However, RPMD tends to underestimate the reactivity at higher collision energies. We saw above in the previous convergence studies that ICSs around the threshold are well reproduced by RPMD for low values of β and that RPMD ICSs tend to differ more with different β values at higher collision energies. Our ansatz for β is a first attempt at inferring the right physical choice for β but, as we mentioned for the ICSs, seems to overestimate its values at high collision energies. Taking α = 1, on the other hand, better renders the ICSs at higher collision energies but fails to reproduce the right shape of the threshold. We are currently working on a more systematic choice of β in the framework of our method.
To the best of our knowledge, the only attempt at tuning β in the function of the energy was done with a single ring polymer to compute the microcanonical rate for the Eckart barrier with RPMD. [82] Combining QCT and RPMD to describe initial state-selective quantities for triatomic reactive systems was the first natural step. The method presented here can be straightforwardly generalized to polyatomic molecule reactants and can also be directly applied to the study of gas-surface reactions. The initialization of well-defined vibrational initial states using an adiabatic switching procedure adapted for RPMD is currently under investigation. This approach is currently under investigation. In addition, the initialization of distinct rotational and ro-vibrational initial states on the basis of the previous step and the latest QCT techniques adapted to RPMD will improve the applicability of the approach. In future work, the resolution of final ro-vibrational states needs to be worked out. To this end, one can also adapt successful QCT procedures. [36,98] This will allow for the subsequent computation of the state-to-state differential cross sections.

| CONCLUSIONS
In this manuscript, we discussed the robustness, convergence, and accuracy of a very recently proposed method to approximately calculate quantum initial state-selected reactive cross sections, which combine QCT and RPMD. For testing, we focused on prototypical triatomic systems X + H 2 , X = Mu, H, D, Cl, F, with H 2 in it vibrational ground or first excited state, considering collision energies of up to 1.5 eV. We find that the proposed excitation scheme is robust and reasonably accurate and that, even for very light reaction partners, such as Muonium, fewer than 100 ring polymer beads are required for convergence.
Good agreement for the Mu/H/D + H 2 (ν = 0,1) reactions with exact quantum scattering calculations is found over the whole collision energy range. In particular, the case of Mu + H 2 (ν = 0) shows that the proposed approach correctly describes the ZPE constraints posed by the MuH product, and the case of Mu + H 2 (ν = 1) shows that the new approach can approximately describe tunneling through a vibrationally nonadiabatic barrier. The latter feature is also confirmed by the good results for the Cl + H 2 (ν = 1) reaction at lower collision energies. However, larger deviations are found at higher collision energies, where the results coincide with cross sections obtained from QCT.
QCT-RPMD for the F + H 2 reaction slightly improves over QCT simulations but cannot fully reproduce exact QM simulations. For low energies, this difference can partially be due to reactive resonance effects in this reaction, but overall, there are also some shortcomings in our vibrational excitation scheme at low collision energies, which needs improvement. Ideas to improve the current approach in this regard are under investigation. The method can be generalized to polyatomic systems, is numerically efficient and easily implemented, and shows encouraging improvement over QCT simulations. Several extensions to include electronically nonadiabatic effects and the resolution of final states are planned for future work.