Lucien Dupuy

David Lauvergnat

Yohann Scribano
email: yohann.scribano@umontpellier.fr

Smolyak representations with absorbing boundary conditions for reaction path Hamiltonian model of reactive scattering

come

Smolyak representations with absorbing boundary conditions for reaction path Hamiltonian model of reactive scattering

I. INTRODUCTION

Cumulative Reaction Probability (CRP), here designed as N (E), is the key ingredient for the description of elementary chemical reactive scattering and represents one of the major research field of theoretical chemistry and chemical physics [START_REF] Levine | Molecular Reaction Dynamics[END_REF]. In terms of the S matrix theory [START_REF] Tannor | Introduction to Quantum Mechanics: A Time-Dependent Perspective[END_REF], the CRP is defined as the sum of all state-to-state reaction probabilities at a given energy by:

N (E) = nr np S nrnp (E) 2 (1)
The sum is running over all open reactants and products channels (respectively asymptotic states n r and n p ).

The most rigorous way to compute CRP of Eq.( 1), is to perform Time Independent Quantum Mechanical (TIQM) simulations for each collisional energy using the Close-Coupling method [START_REF] Truhlar | Quantum Mechanics of the H + H 2 Reaction: Exact Scattering Probabilities for Collinear Collisions[END_REF][START_REF] Launay | Hyperspherical close-coupling calculation of integral cross sections for the reaction H + H 2 → H 2 + H[END_REF][START_REF] Schatz | Quantum mechanical reactive scattering for threedimensional atom plus diatom systems. II. Accurate cross sections for H + H 2[END_REF]. This formalism is very well suited at low collisional energy for which a small number of reactive and non-reactive channels are involved and leads to the determination of the full state-to-state scattering matrix elements S nrnp (E). This approach was successfully applied to both abstraction and insertions chemical reactions [START_REF] Honvault | Theory of Chemical Reaction Dynamics[END_REF][START_REF] Bañares | Dynamics of the S( 1 D) + H 2 Insertion Reaction: A Combined Quantum Mechanical and Quasiclassical Trajectory Study[END_REF][START_REF] Lique | Ortho-para-H 2 conversion by hydrogen exchange: Comparison of theory and experiment[END_REF][START_REF] Honvault | State-to-State Quantum Mechanical Calculations of Rate Coefficients for the D + + H 2 → HD + H + reaction at Low Temperature[END_REF][START_REF] González-Lezana | The D + + H 2 Reaction: Differential and Integral Cross Sections at Low Energy and Rate Constants at Low Temperature[END_REF].

Many years ago, Miller and co-workers have proposed a direct calculation of the thermal rate constant from the formalism of the flux correlation functions [START_REF] Yamamoto | Quantum Statistical Mechanical Theory of the Rate of Exchange Chemical Reactions in the Gas Phase[END_REF][START_REF] Miller | Quantum mechanical transition state theory and a new semiclassical model for reaction rate constants[END_REF]. The main advantage of the Miller re-formulation is its efficiency for directly computing N (E) (and also the thermal rate constant) which does not require the knowledge of individual S-matrix elements unlike the S-matrix Kohn technique [START_REF] Colbert | A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix Kohn method[END_REF]. Based on this efficient reformulation of the CRP, Seideman and Miller [START_REF] Seideman | Calculation of the cumulative reaction probability via a discrete variable representation with absorbing boundary conditions[END_REF][START_REF] Seideman | Quantum mechanical reaction probabilities via a discrete variable representation-absorbing boundary condition Green's function[END_REF] have then introduced a Discrete Variable Representation (DVR) for both the Hamiltonian and the flux operator. By imposing absorbing boundary conditions (ABC) they were able to get a well-behaved representation of the outgoing wave Green's function operator (leading to the DVR-ABC method). Moreover, Manthe and Miller [START_REF] Manthe | The cumulative reaction probability as eigenvalue problem[END_REF] have done an important progress in the reformulation of the CRP calculation by proposing to compute it as the trace of an operator, the reaction probability operator P (E):

N (E) = Tr P (E) = k p k (E) = 4 Tr ˆ r Ĝ(E) ˆ p Ĝ(E) † (2)
with the natural reaction probabilities p k (E) [START_REF] Manthe | S-matrix decomposition, natural reaction channels, and the quantum transition state approach to reactive scattering[END_REF][START_REF] Ellerbrock | Natural reaction channels in H + CHD 3 → H 2 + CD 3[END_REF] which all lie between 0 and 1. In the above equation, ˆ r,p are absorbing potentials at the boundaries of the interaction region, respectively towards the reactant and product regions of the potential energy surface and Ĝ(E) = (E + iˆ -Ĥ) -1 is the Green's function operator and Ĥ the Hamiltonian operator.

This approach was used by many other groups [START_REF] Leforestier | Quantum mechanical calculation of the rate constant for the reaction H + O 2 → OH + O[END_REF][START_REF] Viel | Quantum mechanical calculation of the rate constant for the reaction H + O 2 → OH + O[END_REF][START_REF] Poirier | Quantum reactive scattering for three-body systems via optimized preconditioning, as applied to the O + HCl reaction[END_REF][START_REF] Woittequand | Iterative time independent calculation of the cumulative reaction probability within a basis adapted preconditioner[END_REF] for (3D) reactive processes but was not applied to more challenging problems (with large number of Degrees of Freedom, DOF) due to its intrinsic limitation. Indeed, as in any other eigenvalue problems (bound states calculations), performing simulations for systems with many degrees of freedom require the use of huge sets of basis functions. This is the famous curse of dimensionality, a well-know problem met in vibrational bound state computations.

In this work, we propose to extend the ability to use the Miller and co-workers's method (the direct calculation of the CRP) to multi-D scattering. To achieve this goal, we propose to use an efficient representation basis based on the concept of sparse grids and Smolyak representations introduced many years ago by Smolyak.

[23] The Smolyak's schemes were more recently introduced by Avila and Carrington [START_REF] Avila | Nonproduct quadrature grids for solving the vibrational Schrödinger equation[END_REF][START_REF] Avila | Using nonproduct quadrature grids to solve the vibrational Schrödinger equation in 12D[END_REF][START_REF] Avila | Using a pruned basis, a non-product quadrature grid, and the exact Watson normal-coordinate kinetic energy operator to solve the vibrational Schrödinger equation for C 2 H 4[END_REF] for vibrational spectroscopy of semi-rigid molecules and also proposed by Lauvergnat and Nauts [START_REF] Lauvergnat | Torsional energy levels of nitric acid in reduced and full dimensionality with ElVibRot and Tnum[END_REF][START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF][START_REF] Nauts | Numerical on-the-fly implementation of the action of the kinetic energy operator on a vibrational wave function: application to methanol[END_REF] for molecular spectroscopy of systems with one DOF with large amplitude motion. More recently, it was extended to the eigenvalue calculation of confined diatomic molecule [START_REF] Felker | Intramolecular stretching vibrational states and frequency shifts of (H 2 ) 2 confined inside the large cage of clathrate hydrate from an eight-dimensional quantum treatment using small basis sets[END_REF][START_REF] Lauvergnat | and D 2 in the small cage of structure II clathrate hydrate: Vibrational frequency shifts from fully coupled quantum six-dimensional calculations of the vibration-translation-rotation eigenstates[END_REF][START_REF] Benoit | Does cage quantum delocalisation influence the translationrotational bound states of molecular hydrogen in clathrate hydrate?[END_REF][START_REF] Powers | The effect of the condensed-phase environment on the vibrational frequency shift of a hydrogen molecule inside clathrate hydrates[END_REF] and to van der Waals complexes. [START_REF] Avila | Toward breaking the curse of dimensionality in (ro)vibrational computations of molecular systems with multiple large-amplitude motions[END_REF][START_REF] Avila | Full-dimensional (12D) variational vibrational states of CH 4 •F : Interplay of anharmonicity and tunneling[END_REF][START_REF] Avila | Exact quantum dynamics background of dispersion interactions: case study for CH 4 •Ar in full (12) dimensions[END_REF] The present work involves the use of a Smolyak scheme for the representation of the cumulative reaction probability operator P (E) in combination with an iterative Lanczos algorithm in order to determine its eigenvalues.

The outline of the paper is as follows. In Sec.II, we first describe the theoretical and computational scheme used within this work as well as the multidimensional chemical reaction model used as a benchmark on which we applied it. Sec.III presents our results using the proposed numerical scheme and finally, Sec.IV concludes and presents perspectives of this work.

The concept of Smolyak scheme was introduced [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF] many years ago and has seen renewed interest in recent years. In its translated English paper in 1963 [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF], Smolyak proposed two equivalent schemes to circumvent the exponential scaling of the basis set or grid sizes when the number of 1D-basis sets or the number of degrees of freedom, n, increases. With respect to others basis set selection or pruning approaches [START_REF] Dawes | How to choose one-dimensional basis functions so that a very efficient multidimensional basis may be extracted from a direct product of the one-dimensional functions: energy levels of coupled systems with as many as 16 coordinates[END_REF], the Smolyak schemes enable to use a consistent way to push away the exponential scaling for both the basis set and the grid.

In its second scheme, Smolyak has shown that a large direct-product can be substituted by a sum of small direct-products:

S rep = L-n+1≤| |≤L (-1) L-| | C L-| | n-1 • S 1 1 ⊗ S 2 2 ⊗ • • • ⊗ S n n with | | = n i=1 i (3)
where

C L-| |
n-1 are binomial coefficients. In this Eq. ( 3), one can select the small directproducts (S

1 1 ⊗ S 2 2 ⊗ • • • S n n ) kept
L -n + 1 ≤ | | ≤ L (4)
where, L is a parameter which controls the number of Smolyak terms in the sum and therefore the basis set (N b ) and grid (N q ) sizes. In Eq. (3), the S i i can be a finite basis set (B i i with nb i i basis functions) or finite quadrature grid points (G i i with nq i i grid points and weights).

Formally, nb i

i and nq i i can be defined by any increasing sequences in i . In the present study, their expressions are given by:

nb i i = A i + B i • i and nq i i = nb i i ( 5
)
Within this work, we keep the same value of the L parameter for both Smolyak's representations (for the grid representation and the finite basis representation). This choice was done regarding several benchmark calculations which have shown that no better convergence is obtained introducing two distincts parameter L B and L G as it was done previously in vibrational bound states calculations studies. [START_REF] Lauvergnat | Quantum dynamics with sparse grids: a combination of Smolyak scheme and cubature. Application to methanol in full dimensionality[END_REF][START_REF] Lauvergnat | and D 2 in the small cage of structure II clathrate hydrate: Vibrational frequency shifts from fully coupled quantum six-dimensional calculations of the vibration-translation-rotation eigenstates[END_REF] Then, any vector |u can be represented as a "weighted" sum of Smolyak contributions u 1 , 2 ••• n as:

|u = L-n+1≤| |≤L (-1) L-| | C L-| | n-1 • u 1 , 2 ••• n (6)
and each contribution is expanded on a small direct-product

S 1 1 ⊗ S 2 2 ⊗ • • • S n
n parametrized by the i of either restricted basis sets or restricted grids. In terms of basis functions, the contribution sizes, are n i=1 nb i i and, they can be different from one another. Furthermore, since a multidimensional basis function can present in several Smolyak contributions, we need a mapping table between this basis function in a Smolyak contribution and the one in |u . Therefore, the size of |u is smaller than than the sum of the n i=1 nb i i . This representation, as a sum of Smolyak contributions u 1 , 2 ••• n of any vector is fundamental so that the action of an operatorcan be split on the different Smolyak contributions.

B. Time-independent quantum mechanical method

In the present work, we compute the cumulative reaction probability as defined by Seideman and Miller [START_REF] Seideman | Calculation of the cumulative reaction probability via a discrete variable representation with absorbing boundary conditions[END_REF]:

N (E) = Tr Γr . Ĝ(E). Γp . Ĝ † (E) (7)
with Ĝ(E) = (E -Ĥ + iˆ ) -1 and ˆ = 1 2 Γr + Γp [START_REF] Lique | Ortho-para-H 2 conversion by hydrogen exchange: Comparison of theory and experiment[END_REF] and where Γα = 2ˆ α for α ≡ {r, p} (reactant and product). In the case of systems with many degrees of freedom, it was well established that the trace of the cumulative reaction operator as defined by Eq.( 2) can be computed without explicitly computing (and storing) the full matrix representation of cumulative reaction probability operator. Indeed, we only need to evaluate the action of the operator P (E) on a vector and its eigenvalues will thus be easily computed by an iterative eigenvalue solver. The iterative method used in this work is the one implemented in the ARPACK library [START_REF] Lehoucq | ARPACK Users' Guide[END_REF] and based on the Implicitly Restarted Lanczos algorithm (IRL). As a Krylov space method, it estimates a set of extremal eigenvalues by computing successive applications of P (E) to a randomized starting vector. Given that the basis set is well suited for its representation, a low number of Lanczos iterations in comparison to the initial size of the working basis is enough to get an accurate estimate of N (E).

The number of required iterations scales with the number of non negligible eigenvalues which needs to be accounted for.

The main task of the iterative Lanczos scheme to compute the CRP is thus linked to our ability on applying; at a given total energy E; the Green's function Ĝ(E) on a given vector expanded in the working basis set. Since we have an efficient compact Smolyak finite basis representation, we can build and store the full matrix representation of any operator involved in the definition of the cumulative reaction operator. This allows us to easily evaluate the action of this former operator by simple matrix-vector products. Each action of the cumulative reaction probability operator on a vector requires two operations of the Green's function operator onto a vector (acting sequentially all various matrices).

The action of the Green's function operator (using its matrix representation) onto a vector requires that we solve the set of linear equations :

G -1 (E).w = u (9)
where G -1 (E) = (E + i -H) and the vector u is known. We solved the linear equation system using a LU decomposition [START_REF] Press | Numerical Recipes in FORTRAN[END_REF] of the inverted Green's function matrix appearing in Eq.( 9). This decomposition is done only once for a given energy E. The action of the Green's function (leading to the vector w) being subsequently obtained by backward substitution.

C. Multidimensional chemical reaction models

To highlight the efficiency of our methodology, we apply it to a multidimensional model of a chemical reaction with many degrees of freedom. We consider the nD dimensional bottleneck Hamiltonian model defined as:

Ĥ = p2 x 2m + V 0 (x) + n i=2 p 2 y i 2m + V c (x, y) (10)
where V 0 (x) is the one-dimensional potential defined along the rectilinear reaction coordinate

x and V c (x, y) represents the multidimensional potential describing the coupling between the

x coordinate and all the other bath coordinates y i (that we can represent by y in a vectorial

V 0 (x) = V 0 sech 2 (αx) and V c (x, y) = m n i=2 ω 2 0 2 1 + b sech 2 (αx) y 2 i ( 11
)
where along the x coordinate the potential shape is that of an Eckart barrier, while along the other y i coordinates the potential is shaped as an harmonic oscillator. Parameters involved in the potential function (m is the mass of the particle) are chosen as m = 1060 a.u., α = 1 a.u., V 0 = ω 0 = 0.425 eV. An illustration of the potential is reported on Fig. 1,

where contour lines of a 2D potential (defined by a reaction coordinate x coupled with a single harmonic mode y) are plotted.

D. Computational details

The calculation of the CRP is performed within the approximation that the Green's function representation is done with optical potentials located in the two asymptotic regions, the left (reactants) and the right (products). Considering the Implicitly Restarted Lanczos parameters, the number of requested eigenvalues was taken equal to the number of eigenreaction probabilities n oc at a given energy, and the length of the Arnoldi factorization was set to 2n oc as suggested in ARPACK documentation. The number of restarts needed to achieve convergence is between 2 and 16 depending on the system and energy considered.

In this work, we have used an optical potential proposed by González-Lezana et al. [START_REF] Gonzalez-Lezana | Quantum reactive scattering with a transmission-free absorbing potential[END_REF] and derived from a semi-classical approximation [START_REF] Manolopoulos | Derivation and reflection properties of a transmission-free absorbing potential[END_REF] given by:

r,p (x) = 0 4 (κ r,p (x) + c) 2 + 4 (κ r,p (x) -c) 2 - 8 c 2 (12)
with

κ p (x) =      x -x p W p , x > x p 0 , x < x p κ r (x) =      - x -x r W r , x < x r 0 , x > x r (13)
and where c = 2.62206 as proposed in Ref. [START_REF] Gonzalez-Lezana | Quantum reactive scattering with a transmission-free absorbing potential[END_REF].

In Eq.( 13), x designs the coordinate along which the optical potential is introduced and

x r,p are the onset positions respectively in the reactants and products asymptotes of the absorbing regions (x p > 0 > x r ). The width of this optical potential form is given by c × W r,p . We chose aforementioned parameters to ensure Green's function representation will be accurate down to a minimum collision energy E min = 0.0105 a.u., which is associated to a value of N (E min ) = 1.7028 × 10 -3 in the one-dimensional scattering through the Eckart barrier V 0 (x). This criterion sets 0 = E min and W p,r = 1.33 a.u.. Onset positions are

x p,r = ± 4 a.u. so that the absorbing bands are far enough from the interaction region.

In the following chemical reaction model, we use the Smolyak representation presented in section II A for all coordinates. The Krylov space associated to the probability reaction operator is built in the finite basis Smolyak representation, using specific basis for each kind of coordinates. For the reactive coordinate x we use particule-in-a-box (PIB) basis functions ψ j (x) = 2/a sin(jπ(x/a + 1/2)) where a = x px r + (W p + W r ).c is the width of the box so that the range of x is [-a/2, a/2]. For all the others remaining coordinates, we use harmonic oscillator (HO) basis functions with a coordinate rescaling parameter β = √ mω 0 suited for the potential. In the following calculations, the same constraint parameter L is used for both grid and finite basis representation. For both basis set types, we use quadrature grid points adapted to the basis: a set of equally space grid points for the PIB basis set and gauss Hermite grid points for the HO basis set.

III. RESULTS AND DISCUSSION

Multidimensional uncoupled model: comparison to analytical results

In order to benchmark our approach, we start with a un-coupled multidimensional model which can be solved exactly. In the uncoupled case (b = 0) an analytical result is obtainable by microcanonical convolution of the one-dimensional CRP through the Eckart Barrier, for which an analytical formula is known [START_REF] Seideman | Calculation of the cumulative reaction probability via a discrete variable representation with absorbing boundary conditions[END_REF][START_REF] Eckart | The Penetration of a Potential Barrier by Electrons[END_REF].

N nD (E) = v N 1D (E -E v ) with E v = n i=2 h ω 0 v i + 1 2 ( 14
)
where E is the total energy. The sum in Eq.( 14) carries on all energetically allowed channels (E v < E). The Smolyak parameters of the sequences seen in Eq.( 5) are such as (A x = 50, B x = 50) for the x reaction coordinate and (A i = 1, B i = 1) for all other coordinates (for i = 2, . . . , n). For the (1+2)D system model (n = 3), both the complete sum and the channel bundles associated to the same total vibrational number |v| = The lowest eigenvalues computed in each case are given in Tab.I along with N (E). Both cumulative reaction probability N (E) and individual eigenreaction probabilities p k are well converged.

Multidimensional coupled model

To benchmark our approach, for multidimensional coupled system, we first consider a coupled 3-dimensional system model with a coupling parameter b = 0.1. Now, the 1-D basis functions associated to y i coordinates are no longer partial eigenfunctions of V c (x, y). Hence, a larger constraint parameter L will be needed to construct an accurate representation of scattering states and obtain converged CRP values. Smolyak parameters are identical to the ones in the previous section, with decreased reaction coordinate parameters (A x = 40, B x = 40). On Fig. 3 a comparison of N (E) obtained with Smolyak constraint parameters from 1 to 3 and an ABC-DVR calculation with energy cutoff is drawn.

In the range considered, that is up to |v| = 2, an excellent agreement is observed between the results for L = 3 and the calculations using direct product DVR. By decreasing L, scattering states associated to higher vibrational excitations |v| = L cannot be accurately represented with the remaining Smolyak terms as they lack excited HO basis functions.

The corresponding eigenreaction probabilities are very similar to those of the uncoupled potential model, as one would expect. In this simple case, choosing L = |v max | + 1 is enough to maintain accuracy.

Secondly, we applied our numerical scheme for the CRP calculation of a 25D system model. The same convergence with respect to L was observed in those calculations keeping other Smolyak parameters unchanged. The CRP of this 25-dimensional coupled bottleneck system model obtained with L = 2 is drawn in Fig. 4. In this figure, we have also reported two approximated results obtained from reduced dimensional calculations using adiabatic and energy shifting models as proposed by Bowman. [START_REF] Bowman | Reduced Dimensionality Theory of Quantum Reactive Scattering[END_REF] The adiabatic model shows a global better agreement with our exact results than the energy shifting model. Our full dimensional calculation highlights the ability of our procedure to investigate high dimensional scattering processes. In Tab. II are reported the sizes of the Smolyak basis and grid for multidimensional bottleneck system models considered (from 3D to 25D), with L = 2, and the corresponding direct product grid size is also given as a comparison. Smolyak grid size benefits from a substantially better scaling law with the system's size. Ending up 10 2 times smaller than the direct product grid for the 9D system, its use allows a valuable reduction in numerical requirements. This efficiency gap between a direct product grid and an adapted Smolyak scheme continues to widen for bigger systems. The ratio is of the order of 10 8 for 25 DOF, making the computation unfeasible using the direct product grid. The Smolyak basis compactness makes it the most advantageous to carry the final steps of linear algebra described in section II B.

On Fig. 5 we report the size of the Smolyak finite basis representation used to compute the eigenvalues of the reaction probability operator. We show that the size grows reasonably with the dimensionality of the system model. Even for 25 DOF, the size of the basis (14 040 basis functions) is very far from the size that could reach a direct product basis. We also oberve on this figure that advantage of the Smolyak basis is even more confirmed if we increases the value of the Smolyak parameter L since the Smolyak basis do not increase dramatically increasing L by unit.

IV. CONCLUSIONS

In this paper, we have applied an efficient Smolyak scheme in a multidimensional chemical reactive scattering model (bottleneck potential). Our 25D simulations have highlighted that strong convergence of the CRP calculations can be reached using moderate number of Smolyak basis functions. The high reduction of the basis size allow us to use the full matrix representation of all operators used to define the probability reaction operator P (E). Within this work, we ensure a converged CRP calculations using of a basis truncation criteria, the Smolyak parameter L.

The efficient contracted Smolyak representation with absorbing boundary condition (SR-ABC method) used in this work will allow us to investigate high dimensional chemical reactions for which most of the previous theoretical studied are based on a the use of small reduced dimensional models. The feasibility of our method to increase the number of active coordinates in any reduced dimensional scheme should be as important than to increase the accuracy of the electronic structure calculation as reported by Farahani et al. in their study of a S N 2 mechanism. [START_REF] Farahani | Ab initio quantum mechanical calculation of the reaction probability for the Cl -+ PH 2 Cl → ClPH 2 + Clreaction[END_REF] The present results are very encouraging and the new computational scheme will be used to investigate much more complex reactive systems such as complex abstraction chemical reactions characterised by the formation of a long-lived intermediate complex. Indeed, this approach should be much relevant for more complex chemical reactions such as insertion reactions for which calculations are notoriously more difficult. Further extensions of the present work could be investigated for a more realistic representation of the chemical reactivity. It is for example possible to relax the actual rectilinear model using a curvilinear Reaction Path Hamiltion (RPH) formulation. We can indeed use an exact numerical derivation of the Kinetic Energy Operator for the RPH using the most adapted set of curvilinear coordinates as provided by the numerical code TNUM. [START_REF] Lauvergnat | ElVibRot-TnumTana, quantum dynamics code[END_REF] function application on a given Lanczos vector. Indeed, the LU matrix decomposition scales in CPU time as N 3 and requires storage of at least the N × N Hamiltonian matrix, will be a limiting factor for very high dimensional problems. The use of efficient preconditioner [START_REF] Poirier | Optimized preconditioners for Green function evaluation in quantum reactive scattering calculations[END_REF] in the Green's function calculation is also a possibility to increase again the efficiency of our numerical scheme of CRP's calculation. Such development is actually investigated in our laboratories. Moreover the present methodology based on the Green's function calculation with absorbing boundary conditions can be used for other processes than bimolecular reactions, as for example photodissociation processes (or radiative association, the inverse process), which can be investigated in a very similar way to what has been presented in this article.

V. ACKNOWLEDGMENTS

N 3D (E) E [hω 0 ] Analytical v = 0 v = 1 v = 2 L = 1 L = 2

  in the sum satisfying the constraint on | | related to the definition of the Smolyak scheme:

n i=2 v i are represented in Fig. 2 .

 2 Additionally, two calculations using different parameters L B are drawn. The calculation with L B = 1 (red circles) is accurate in the energy range between 1.5 and 3.5hω 0 . It does not account for the contribution of additional opened channels at higher energies, as this Smolyak scheme contains eigenfunctions of harmonic oscillators up to |v| = 1, which are sufficient to represent scattering states from the 2 first bundles only. Oscillations seen above E = 4hω 0 are due to a lack of basis functions along the reaction coordinate x associated to higher momentum. Of course, Smolyak parameters A x = 50 and B x = 50 could ( and should ) be decreased to match the energy domain of accuracy for the y modes Smolyak parameters (A i , B i ) in order to minimize computational cost. Still, it is illustrative to see the effects of insufficient basis sizes for reaction and bath coordinates separately. Increasing L B to 2 increases both the maximum |v| for which harmonic eigenfunctions are included in the Smolyak sum and the highest collisional energy along x accurately resolved. 3 more channels associated to |v| = 2 are accounted for (red dots), in agreement with the analytical result. In Tab. (I) are reported the number of converged eigenvalues p k requested from the Lanczos algorithm for different values of E. Each eigenvalue is associated to a different energetically allowed channel from equation (14) and satisfies 0 ≤ p k ≤ 1.

FIG. 2 .FIG. 4 .

 24 FIG.2. Cumulative reaction probability for the 3D uncoupled bottleneck. In black continuous curve is the analytical result, in dashed curves are the analytical contributions from each bundle of channels associated to a total vibrational number. In red circles is the numerical result obtained for L=1, in red dots for L=2.

FIG. 5 .

 5 FIG. 5. Smolyak FBR basis size as a function number of degrees of freedom (DOF) and for different values of L. Values of other Smolyak parameters correspond to the coupled (b = 0.1) bottleneck case.

TABLE I .

 I This work has been supported by the Agence Nationale de la Recherche (ANR-HYTRAJ), Contract No. ANR-19-CE30-0039-01. Y.S. and L.D. thank support from the High Performance Computing Platform MESO@LR at the University of Montpellier. Cumulative reaction probabilities computed for the 3D uncoupled bottleneck potential model with Smolyak constraint parameter L = 2 compared with the analytical result obtained by microcanonical convolution (see the text). For each energy is also reported the number of eigenreaction probabilities n oc as well as its smallest value p min k .

	Energy (hω 0 ) 0.1	n oc		p min k					CRP (3D)
								This work		Analytical
	1.667		1		1.50280×10 -3		1.50280×10 -3	1.51411×10 -3
	2.163	0.05	1		9.51417×10 -1		9.51417×10 -1	9.51155×10 -1
	2.863		3		8.06611×10 -2		1.16132		1.16245
	3.1545 3.738 4.117 y [a.u.]	0	3 6 6		9.44637×10 -1 6.98826×10 -3 8.98919×10 -1		2.88927 3.02095 5.69676		2.88849 3.02103 5.70027
		-0.05
		-0.1 -4	-3	-2	-1	0 x [a.u.]	1	2	3	4

FIG. 1. Contour lines illustrating the coupled bottleneck potential. The reaction coordinate is labelled by x and y is one of the equivalent coupled harmonic degrees of freedom.