Chemical reactions in multidimensional driven systems are typically described by a time-dependent rank-1 saddle associated with one reaction and several orthogonal coordinates (including the solvent bath). To investigate reactions in such systems, we develop a fast and robust method—viz., local manifold analysis (LMA)—for computing the instantaneous decay rate of reactants. Specifically, it computes the instantaneous decay rates along saddle-bound trajectories near the activated complex by exploiting local properties of the stable and unstable manifold associated with the normally hyperbolic invariant manifold (NHIM). The LMA method offers substantial reduction in numerical effort and increased reliability in comparison with direct ensemble integration. It provides an instantaneous flux that can be assigned to every point on the NHIM and which is associated with a trajectory—regardless of whether it is periodic, quasiperiodic, or chaotic—that is bound on the NHIM. The time average of these fluxes in the driven system corresponds to the average rate through a given local section containing the corresponding point on the NHIM. We find good agreement between the results of the LMA and direct ensemble integration obtained using numerically constructed, recrossing-free dividing surfaces.

## I. INTRODUCTION

The dynamics of a chemical system can often be described by the classical equations of motion for selected coordinates driven by a Born-Oppenheimer potential. The transition from reactants to products in such a chemical reaction is typically marked by a barrier region with a rank-1 saddle that has exactly one unstable direction called the reaction coordinate, while the remaining degrees of freedom are locally stable and are associated with other bound internal motions and external bath coordinates. The reaction can be described both qualitatively and quantitatively within the framework of Transition State Theory (TST). This theory rests on the identification of a recrossing-free dividing surface (DS) in the barrier region, which separates reactants from products,^{1–6} and has been applied in a wide range of fields, as we have noted in our previous work (see, e.g., in Ref. 7), and in recent advances in chemical reaction dynamics.^{8–11}

The construction of recrossing-free DSs is a nontrivial task that can be achieved by using perturbation theory to construct good action-angle variables associated with the DS.^{12–14} The theory was generalized for chemical reactions under time-dependent conditions, arising from driving, noise, or both,^{15,16} and has been applied to chemical reactions with few degrees of freedom including, e.g., H + H_{2},^{17,18} LiCN,^{19} and ketene isomerization.^{20}

Typical chemical reactions, however, require representations with higher dimensionality coupled to complex environments. In a multidimensional autonomous Hamiltonian system, a recrossing-free DS attached to the normally hyperbolic invariant manifold (NHIM) contains all trajectories that are trapped in the saddle region when propagated both forward and backward in time.^{7,21} This manifold can be constructed approximately using normal form expansions^{12,14,17,22–28} or exactly by numerical application of Lagrangian descriptors,^{7,21} the binary contraction method,^{29} and machine learning algorithms.^{7,30}

In a driven system with one degree of freedom, the rate constant can be obtained by propagating a large ensemble of trajectories in the vicinity of the Transition State (TS) trajectory—which coincides with the time-dependent NHIM in a one-dimensional system—and observing the decay in the number of reactant trajectories, which have not yet crossed the DS. As a numerically less expensive alternative, the rate constant can be computed using the Floquet exponents of the TS trajectory.^{31} However, the situation is less clear for systems with *d* ≥ 2 degrees of freedom because the dimensionality of the time-dependent NHIM is (2*d* − 2) > 0, and hence there is not a unique trajectory from which to obtain a single set of stability exponents. Indeed, the additional degrees of freedom enable a nontrivial dynamics on the NHIM^{32} and substantially extend the possibilities for trajectories to cross a DS close to the NHIM.

The purpose of this paper is to address this challenge in obtaining rate constants in multidimensional driven systems with a rank-1 saddle coupled to a formally arbitrary number of orthogonal modes. We find that we can assign instantaneous reaction rates to every point on the time-dependent NHIM by observing the change in the number of trajectories crossing the numerically constructed time-dependent DSs. We also address the question of how appropriate initial conditions of trajectory ensembles must be chosen to compute the instantaneous rate through a local section of the NHIM. Trajectories having just enough energy to overcome the time dependent barrier will cross the DS close to the NHIM. They correspond to the slowest flux and are therefore critical in determining instantaneous rates, which are approximately given by the decay rates of the TS, at least in the local vicinity of the NHIM. The time averages of these instantaneous rates yield rate constants associated with periodic, quasiperiodic, or even chaotic trajectories on the NHIM. Furthermore, we present a very efficient numerical method for the computation of instantaneous reaction rates associated with arbitrary points on the NHIM, which does not require the propagation of trajectory ensembles. This local manifold analysis (LMA) method uses local properties of the stable and unstable manifolds at points on the NHIM. Reaction rates calculated using the LMA method agree with those obtained through the more numerically expensive sampling of a trajectory ensemble.

The paper is organized as follows. In Sec. II, we introduce our methods for the computation of reaction rates based on the propagation of trajectory ensembles in Sec. II A, on the Floquet exponents of trajectories on the NHIM in Sec. II B, and on local properties of the stable and unstable manifolds at points on the NHIM in Sec. II C. Numerical results and comparisons are presented in Sec. III. As this work is focused on advancing methods so as to ultimately calculate rate constants of multidimensional chemical reactions, numerical results are restricted to a periodically driven two-dimensional model system useful for verifying the theory for periodic and general nonperiodic trajectories in Secs. III A and III B, respectively.

## II. THEORY AND METHODS

In this section, we investigate rate constants for a driven chemical reaction in a system with *d* ≥ 2 degrees of freedom. More precisely, the system is described by a moving rank-1 saddle with one unstable reaction coordinate *x* and *d* − 1 stable bath coordinates ** y**. The calculation of rate constants in Transition State Theory (TST) usually implies that the pathway of reactive trajectories is located close to the Transition State (TS), as illustrated in Fig. 1. Exceptions are roaming reactions

^{33–37}with pathways far away from the TS, which we do not consider in this paper. We also do not investigate the Kramers rate

^{38}for an ensemble of trajectories prepared in a reactant well far away from the saddle. While the Kramers rate depends on the thermally distributed kinetic energy of the trajectories in the well, in this paper, the term

*rate*refers to the decay rate of individual trajectories in the activated complex near the TS, which is directly related to geometrical and dynamical properties of the barrier itself, but not necessarily to properties of the reacting ensemble directly.

The effect of the barrier and its driving on reaction rates is high for trajectories that react very closely to the TS and therefore stay in the saddle region for a long time. In one-dimensional static systems with an energy barrier, the pointlike TS defines the minimum energy required for a reaction. It also yields the precise separation between reactants and products, and therefore a recrossing-free dividing surface (DS) can be attached to this point. If the system is subject to time-dependent driving, this point becomes the one-dimensional TS trajectory, and a time-dependent recrossing-free DS can be attached to the TS trajectory to separate reactants from products.^{39}

In multidimensional systems with a time-dependent rank-1 saddle, the situation is more complex. The normally hyperbolic invariant manifold (NHIM) becomes a higher-dimensional manifold and is no longer a single TS trajectory. Consequently, reacting trajectories have additional degrees of freedom associated with the crossing of the DS close to the NHIM. In Ref. 32, the reaction rate for a specific periodic trajectory of a higher-dimensional system has been obtained using the Floquet exponents of that trajectory. Here, we present three methods to calculate reaction rates either using or inspired by this trajectory-based approach focused on the flux through arbitrary points on the NHIM. The first method discussed in Sec. II A is based on the propagation of trajectory ensembles with appropriately chosen initial conditions. The second, numerically less expensive method presented in Sec. II B is an extension of the Floquet method of Ref. 31 to multidimensional systems. The third method derived in Sec. II C uses local properties of the stable and unstable manifolds for the computation of the instantaneous reaction rates.

### A. Ensemble method

#### 1. Preparation of trajectory ensembles

Reaction rates are usually obtained by propagating a large number of trajectories and by measuring their flux from a reactant to a product state. When propagating ensembles, the accurate preparation in the reactant state is crucial to obtaining conclusive results for the instantaneous rate. As these rates are a property of the (possibly time-dependent) saddle, it is important that an ensemble is prepared in a way that many trajectories are affected by the saddle. If most of the trajectories would react at an energy much higher than the barrier, they will most likely cross the barrier region very fast and their dynamics will be mostly unaffected by the saddle. Because of this, it is important to initialize an ensemble close to the stable manifold, where trajectories closely approach the NHIM and are therefore strongly affected by the saddle-bound dynamics on the NHIM.

In Fig. 2, a schematic of the stable $Ws$ and unstable manifold $Wu$ in a (*x*, $vx$) section of a multidimensional system is shown. The limits of this section are small enough, such that the dynamics close to the NHIM are linear, and this requires Δ*x* to also be small. We attach a recrossing-free DS to the NHIM, which is shown as a vertical, dashed line through the intersection of $Ws$ and $Wu$ in Fig. 2.

As highlighted in Fig. 2 by the red dots, an ensemble of *N*_{react}(0) equidistant reactants is initialized close to the NHIM, on a line segment parallel to $Wu$, spanning from a point on $Ws$ to the DS. If a system is time-dependent, the ensemble preparation procedure can be repeated for any point on the NHIM at a time *t*. How often and at which times *t* such an ensemble needs to be prepared is determined by the rate calculation and will be explained Sec. II A 2.

#### 2. Instantaneous ensemble rates

Having prepared an ensemble of reactive trajectories according to Sec. II A 1, these trajectories need to be numerically propagated in time. Starting with a reactant population of *N*_{react}(0), each time a trajectory reacts or, more precisely, pierces the recrossing-free DS, the reactant population decreases by one. For a system that is not time-dependent, the number of reactants *N*_{react} will decrease exponentially fast

with *k* being the rate constant of this reaction. If a barrier is time-dependently driven, depending on the driving itself as well as on the initial time *t*, an ensemble is prepared (here, e.g., for periodically driven systems, the phase relative to the barrier’s oscillation is important), the decrease in *N*_{react}(*t*) will be more or less exponential. So fitting a rate constant *k* according to Eq. (1) may be approximately possible at best in most cases. If chosen wrong, the parameter Δ*x* introduced in Sec. II A 1 when preparing an ensemble will also have a huge influence on the decay of reactants.

The problem of obtaining a precise rate for a nonexponentially decaying reactant population *N*_{react}(*t*) can be solved by using a more fundamental definition of rates. In the framework of chemical kinetics,^{38,40,41} the change in a reactant’s population in first order, unimolecular reactions is proportional to the reactant’s population,

where we allow *k*(*t*) to be time-dependent. In the time-independent case, Eq. (2) would be solved by Eq. (1). However, the instantaneous rate *k*(*t*) yields

Due to the way we measure *N*_{react}(*t*), an exponential decay of reactants also translates into an exponential decay of measurement density for *N*_{react}(*t*) over time, as seen in the step function behavior in Fig. 3. If the time intervals between data points are too small or too large, numerical methods of differentiation may produce erroneous results, which has to be considered during evaluation.

Since the decay is exponential, the number of trajectories necessary to reduce the time intervals between individual reaction events, i.e., to increase the time range where the differentiation of *N*_{react}(*t*) produces accurate results, grows exponentially. Hence, increasing the number of trajectories for a single ensemble will produce diminishing results. For this reason, we start several ensembles at different initial times *t* and piece together the results of each ensemble to fully resolve *k*(*t*) of a specific trajectory on the NHIM with good accuracy everywhere.

Whereas in other references, thousands to millions of trajectories had to be propagated to somehow resolve a rate constant,^{7,21,30,42} with the methods introduced here, we are able to adequately resolve the time-dependent reaction rate *k*(*t*) propagating just a few hundreds of trajectories.

### B. Floquet method

Propagating an ensemble of trajectories to obtain reaction rates is numerically expensive as these ensembles, usually with a large number of reactive trajectories, have to be propagated in time. An alternative method to obtain reaction rates, which does not require the propagation of a trajectory ensemble, has been introduced in 2014 by Craven *et al.* for one-dimensional systems with periodic driving.^{31} The method is based on the stability analysis of the TS trajectory yielding the *Floquet rate*

which depends on the Floquet exponents

of that trajectory. Here, *m*_{u,s} are the eigenvalues of the (2 × 2)-dimensional stability (or monodromy) matrix along the unstable and stable directions, respectively, and *T* is the period of the external driving. The method has already been extended and applied by Tschöpe *et al.* to the TS trajectory of a higher-dimensional system with periodic driving.^{32} In a system with *d* degrees of freedom, *m*_{u,s} in Eq. (5) are the two eigenvalues of the (2*d* × 2*d*)-dimensional stability matrix of the TS trajectory, or more general, of any periodic orbit with any period *T*, with the highest and lowest absolute values.

For a first order, two-dimensional differential equation $\gamma \u0307=f(\gamma ,t)$, where ** γ** ∈

*R*

^{2d}is a vector in phase space, the stability matrix

**(**

*M**t*) is obtained through linearization of the dynamics, which is in turn characterized by the Jacobian

**(**

*J***(**

*γ**t*),

*t*) =

*∂*

**/**

*f**∂*

**of the system**

*γ*with initial conditions ** M**(0) = 1

_{2d}.

For a Hamiltonian system with time-dependent potential *V*(** q**,

*t*) and without friction, Eq. (6) takes the form

In Hamiltonian systems, the monodromy matrix ** M** is symplectic, and hence, if

*λ*is an eigenvalue of

**, then also 1/**

*M**λ*is an eigenvalue, as well as their complex conjugates $\lambda \xaf$ and $1/\lambda \xaf$.

The Floquet method can also be used to calculate the reaction rates of nonperiodic trajectories on the NHIM. In that case, Eq. (5) must be replaced with

In numerical simulations, the time *t* must be chosen sufficiently large to obtain converged estimates of the Floquet rate $k\xafF$ in Eq. (4).

The Floquet method described in Ref. 31 can only be used to obtain mean rate constants averaged over one period of a periodic trajectory or obtained in the limit *t* → ∞ but cannot provide the instantaneous rates *k*(*t*) discussed in Sec. II A 2. The reason is that the directions of the stable and unstable manifolds associated with points on the NHIM change with time, and thus, in general cannot be obtained as eigenvectors of the stability matrix.

### C. Local manifold analysis

Consider an initial ensemble of trajectories in the local vicinity of the stable and unstable manifolds of an arbitrary point on the NHIM in the same way as described in Sec. II A 1, sketched by the red circles in Fig. 2. The idea is not to propagate the individual trajectories of the ensemble by numerical integration of the equations of motion but to use the linearized dynamics in the local vicinity of the stable and unstable manifold to propagate the whole line segment. When going from time *t* to time *t* + Δ*t*, the linearized dynamics can be split into two parts. A compression toward the NHIM in the direction of the stable manifold by the factor Δ*x*^{s}(*t* + Δ*t*)/Δ*x*^{s}(*t*) pulls the line closer to the unstable manifold and a stretching in the direction of the unstable manifold by the factor Δ*x*^{u}(*t* + Δ*t*)/Δ*x*^{u}(*t*) lengthens the line segment. The result of this motion is shown by the red and blue diamonds in Fig. 2. With simple geometry, we can now determine the ratio between the length of the line segment left of the DS (red diamonds in Fig. 2) and the length of the entire line. Due to the properties of the linear mapping that occurred, we know that the density of reactive trajectories along the line is constant. Thus, the obtained ratio of line segments is proportional to the number of reactants *N*_{react}(*t*). By choosing Δ*x*^{s}(*t*) = Δ*x*^{u}(*t*) ≡ Δ*x* as marked in Fig. 2, the number of reactants at time *t* + Δ*t* is given as

It is important to note that the stable and unstable manifolds shown in Fig. 2 can independently rotate around their intersection point as functions of time. However, Eq. (9) stays valid in these cases. One can always separate the linear mapping that occurs into two mappings. First, one that transforms the parallelogram depicted in Fig. 2, such that these rotations are accounted for, all the while keeping Δ*x*^{u}(*t*) and Δ*x*^{s}(*t*) constant, which in turn ensures that the line segment does not cross further into the DS before the next step. Subsequently, we can perform the stretching and compression of the appropriate lengths, as discussed before, to complete the mapping of the linearized dynamics.

From Eq. (9), we obtain the instantaneous rate

Using the linearized equations of motion (7) for the stability matrix ** M**, it can be shown that $\Delta x\u0307s(t)=\Delta vxs(t)$ and $\Delta x\u0307u(t)=\Delta vxu(t)$. Inserting this into Eq. (10) and also explicitly writing the dependencies on the bath coordinates

**and velocities $vy$, which have been dropped in the above derivations, we finally obtain**

*y*This means that for any point on the NHIM parameterized by the bath coordinates ** y** and velocities $vy$, the instantaneous rate

*k*(

**, $vy$,**

*y**t*) is simply given by the difference of the slopes of the stable and unstable manifolds in the (

*x*, $vx$) diagram shown in Fig. 2. Note that for non-Hamiltonian systems, Δ

*ẋ*

^{u,s}depend on the contents of the Jacobian

**in Eq. (6), which, in general, emerges as a prefactor to the slopes in Eq. (11). To determine the slope, one can determine the positions of these manifolds in the immediate vicinity of the NHIM with methods discussed in Ref. 29, where stable and unstable manifolds are interpreted as boundaries between reactive and nonreactive regimes.**

*J*In the derivation of the local manifold analysis (LMA) method, we have assumed that the position of the DS is the same for all particles of the propagated ensemble at time *t* + Δ*t* (see the diamonds in Fig. 2). This assumption is satisfied in the system discussed in Sec. III but may not be so in arbitrary systems. For example, in some multidimensional systems, the orthogonal modes may be coupled to the velocity $vx$ of the reaction coordinate, and thus propagated positions will be associated with different values (** y**, $vy$) of the DS

*x*

_{DS}(

**, $vy$,**

*y**t*+ Δ

*t*). Addressing such challenging cases is a subject for future work.

## III. RESULTS AND DISCUSSION

To verify the methods described in Sec. II, we use the same model for a driven chemical reaction investigated earlier.^{7,21,29,30,32} The system is described by the two-dimensional potential

where a time-dependent oscillating Gaussian barrier with height *E*_{b} separates an open reactant from an open product basin. The barrier is moving along the *x*-coordinate with frequency *ω*_{x} and amplitude $x^$. In the *y*-direction, the dynamics is bound by a harmonic potential with frequency *ω*_{y}. The *x* and *y* coordinates are nonlinearly coupled so that the minimum energy path for a reaction over the saddle is given by *y* = (2/*π*) arctan(2*x*). For simplicity, we use dimensionless units in which the parameters are *E*_{b} = 2, *ω*_{x} = *π*, *ω*_{y} = 2, and $x^=0.4$ if not stated otherwise. The open reactant and product basins of the potential (12) allow us to study different methods on how to obtain reaction rates without having to worry about global recrossings that arise if one or both basins are closed.^{43} However, open reactant and product basins are not a prerequisite for the application of the methods presented here as these methods are based on the local dynamics near the normally hyperbolic invariant manifold (NHIM).

The NHIM of the model system (12) is a two-dimensional time-dependent manifold. The dynamics of trajectories on the NHIM has been studied in Ref. 32. For a periodically driven system with *d* = 2 degrees of freedom, it can be visualized by a Poincaré surface of section (PSOS), e.g., the points (*y*, $vy$) (*t* = *nT*) obtained by using a stroboscopic map in time, with *T* = 2 being the period of the driving and $n\u2208N$. The PSOS for system (12) is presented by the white points in Fig. 4 and exhibits regular toruslike structures. The fixed point with coordinates *y* = 0, $vy$ = −0.72 located at the center of the tori corresponds to a periodic trajectory with the same period as the external driving, and we call this orbit the Transition State (TS) trajectory. Trajectories on the surrounding tori typically show a quasiperiodic behavior. As discussed in Sec. II, rate constants can be assigned to trajectories on the NHIM. This is illustrated for periodic and arbitrary nonperiodic trajectories on the NHIM in Secs. III A and III B, respectively.

### A. Rate calculation for periodic trajectories on the NHIM

We first employ the direct ensemble calculations along the TS trajectory of the two-dimensional model system (12) at different times *t*, according to Sec. II A 1. A typical example of a reactant decay curve *N*_{react}(*t*), initialized with 200 reacting trajectories at *t*_{0} = 0 with a distance Δ*x* = 10^{−3}, is displayed in Fig. 5(a). For several times *t*, visualized by the different colored crosses in Fig. 5(a), the instantaneous derivatives of *N*_{react}(*t*) are shown as straight lines, yielding the instantaneous rates *k*(*t*) at the respective time. As expected, the decay *N*_{react}(*t*) is not purely exponential in the driven case.

The instantaneous rate *k*(*t*), which is given as the most accurate result from the numerical differentiation of the reactant decays *N*_{react}(*t*) of several ensembles (see Sec. II A 2), is therefore not constant in time, as seen in Fig. 5(b). For this particular result, 16 ensembles with 200 trajectories have been launched at equidistant times *t* ∈ [0, *T*], where *T* is the oscillation period of the saddle potential (12).

According to Fig. 5(b), the instantaneous rate *k*(*t*) of the TS trajectory shows an oscillation between *k*_{min} ≈ 3.4 and *k*_{max} ≈ 4.3. The oscillation of *k*(*t*) is twice as fast as the oscillation of the driving in Eq. (12) with *T* = 2 and, hence, of the TS trajectory,

Introduced in Sec. II C, the LMA allows us to obtain instantaneous rates using just the geometrical properties of phase space structures near a single trajectory on the NHIM, without time-consuming propagation of reactive trajectory ensembles. Based on Eq. (11), the instantaneous manifold rate *k*_{M}(*t*) is shown as a thin red dashed line in Fig. 5(b). The perfect agreement between *k*(*t*) and *k*_{M}(*t*) confirms the reliability of the LMA.

Based on a stability analysis of the TS trajectory summarized in Sec. II B, the Floquet rate $k\xafF$ can also be obtained. It is included in Fig. 5(b) as a thick red dashed horizontal line. Since the Floquet rate $k\xafF$ according to Eqs. (4) and (5) is an integrated quantity over the full (periodic) TS trajectory, it is only a single value, valid for the full trajectory. The instantaneous rates *k*(*t*) and *k*_{M}(*t*) seem to oscillate around $k\xafF$. For a more detailed comparison, the mean value of the instantaneous rates $k\xaf=k(t)\xaf=kM(t)\xaf$ is included as a thick black horizontal line in Fig. 5(b). The almost perfect agreement between $k\xaf=3.8067$ and $k\xafF=3.8064$ is clear. This result shows that the additional information encoded in the instantaneous reaction rate nevertheless has a mean value that coincides with the Floquet rate. The latter is an integrated quantity for the stability of the (periodic) trajectory which remains encoded in the integral of the rates obtained by the other two methods.

### B. Rate calculation for arbitrary trajectories on the NHIM

To explain the procedure for how to obtain reaction rates associated with nonperiodic trajectories on the NHIM, we investigate a quasiperiodic trajectory on a stable torus around the TS trajectory shown by the PSOS in Fig. 4. As in the periodic case reported in Sec. III A, an instantaneous rate can be obtained by propagating several ensembles of reactive trajectories that are initialized close to the respective trajectory on the stable torus.

In Fig. 6(a), the instantaneous rate *k*(*t*) is obtained by propagating 80 ensembles of 200 trajectories. Each ensemble is initialized according to the procedure explained in Sec. II A 1 with Δ*x* = 10^{−3} close to the trajectory with initial conditions *t* = 0 on the NHIM at position *y* = 0.4 and $vy$ = −0.75. The initial times *t*_{0} ∈ [0, 40] of any of these ensembles are chosen in equally spaced intervals of Δ*t* = 0.5. Four ensembles are launched per period *T* = 2 of the time-periodically driven saddle according to Eq. (12). As already explained in Secs. II A and III A, the instantaneous rate *k*(*t*) is obtained by numerically deriving the reactant population *N*_{react}(*t*) of these reactive ensembles—cf. the thin black line in Fig. 6(a).

Whereas the instantaneous rate of the periodic TS trajectory in Fig. 5(b) is itself periodic with the period of the external driving, the instantaneous rate *k*(*t*) of the nonperiodic trajectory in Fig. 6(a) clearly does not show such behavior on the time-scale of *T* = 2. Indeed, there seems to be an approximate periodicity on a much larger time-scale with three oscillations in between the two vertical, blue dashed lines. The reason is given by the dynamics itself as the corresponding trajectory on the NHIM is quasiperiodic, meaning it approximately orbits the stable torus three times in this time interval.

Using the LMA explained in Sec. II C, the manifold rate *k*_{M}(*t*) is obtained—cf. the thin red dashed line in Fig. 6(a). Just like in the periodic case of Sec. III A, both *k*(*t*) and *k*_{M}(*t*) coincide perfectly.

As introduced in Sec. II B, according to Eq. (8), a Floquet rate can be obtained from the stability of a nonperiodic trajectory. As infinite time integration cannot be achieved numerically, one has to pay special attention to the convergence of the Floquet rate. The time-dependent logarithmic differences of the eigenvalues of the monodromy matrix are shown in Fig. 7(a). They are obtained according to Eq. (8) and shown as the black dotted curve for an integration along the same trajectory as in Fig. 6(a). Here, a very long integration time (up to *t* = 200 in dimensionless units and corresponding to 100 quasiperiods) is chosen to ensure convergence. The convergence value is extrapolated using linear regression, and a Floquet rate of $k\xafF=3.725$ is obtained as the slope of the data shown in Fig. 7(a). The Floquet rate of the trajectory is also shown as a thick red dashed horizontal line in Fig. 6(a), and the agreement is remarkably good.

Averaging over the three quasioscillations of *k*(*t*) or *k*_{M}(*t*), indicated by the blue dashed vertical lines in Fig. 6(a), yields the average rate $k\xaf=3.720$, included as a thick black horizontal line in Fig. 6(a), that again is in quite good agreement with the Floquet rate $k\xafF=3.725$ obtained by the stability analysis of the underlying trajectory on the stable torus. Hence, even for nonperiodic trajectories, the mean value of the instantaneous rates corresponds to the Floquet rate obtained according to Eq. (8).

This provides a hint to another important result. The mean rate obtained for the periodic TS trajectory in Fig. 5 is different than the rate for the nonperiodic trajectory of Fig. 6(a) on a stable torus around the TS trajectory. This implies that the rate depends on the position of a trajectory on the NHIM. To test this, an alternative calculation not-quite-similar to the one of Fig. 6(a) was performed on a trajectory that starts for *t* = 0 at position *y* = 0.7 and velocity $vy$ = −0.75 further away from the TS trajectory on the NHIM. The results obtained with the same procedure as in Fig. 6(a) are shown in Fig. 6(b). Again, not only is the instantaneous rate *k*(*t*) obtained by propagating an ensemble of reactive trajectories in near perfect agreement with the manifold rate *k*_{M}(*t*) but also the accordance between their average rate $k\xaf=4.255$ and the Floquet rate $k\xafF=4.258$ obtained by linear regression of the blue dashed-dotted curve in Fig. 7(a) is pretty good. Again, the mean rates obtained here are different to the mean rates of the periodic TS trajectory in Fig. 5 or the nonperiodic torus trajectory of Fig. 6(a).

To investigate the trajectory dependence of the rates further, we employed the Floquet method to calculate the mean rates for many trajectories on the NHIM initialized in the intervals *y* ∈ [−0.8, 0.8] and $vy$ ∈ [−2, 2] at time *t* = 0. The results are given as color encoded as per Fig. 4. Note that although we show Floquet rates here, we could just as well use the mean rates obtained from one of the instantaneous rate methods because all of these methods lead to the same average rates.

For the quasiperiodic dynamics on the tori in Fig. 4, any torus is associated with its own mean rate. In an area around the TS trajectory, we observe that the mean rates first decrease, before increasing significantly. Taking a closer look, the rate of the TS trajectory of the model system (12) is a local maximum in the global minimum of reaction rates on the NHIM. This can be seen better by using a *y* cross section through Fig. 4 at the velocity $vy$ of the TS trajectory; see Fig. 7(b). Here, we support the Floquet calculations (black line) with the average rates obtained by the more time-consuming ensemble propagations (red crosses). We also highlight the TS trajectory regarded in Fig. 5 with a red circle and the two trajectories of Fig. 6 with a black and a blue circle.

## IV. CONCLUSION AND OUTLOOK

In this paper, we have shown how to obtain phase-space resolved rates in driven multidimensional chemical reactions. Three methods are compared. The first is based on measuring the flux of trajectories through a recrossing-free dividing surface (DS) attached to the time-dependently moving normally hyperbolic invariant manifold (NHIM) in the barrier region of a rank-1 saddle separating reactants from products. For the second local manifold analysis (LMA) method, rates are obtained directly using just the geometrical properties of phase space close to a single trajectory located within the NHIM of the saddle. No time-consuming trajectory propagation is needed here. The third method is a multidimensional generalization of the Floquet method that has already been applied to obtain the rates of a one-dimensional system.^{31} We provide evidence that the first two methods lead to the same result for the instantaneous reaction rate. Their mean (or long time) average corresponds perfectly to the generalized Floquet rates that are based on the stability analysis of a single trajectory.

Whereas the rates in this paper are the decay rates of trajectories contained in the activated complex near the Transition State (TS) in between reactants and products, the rates measured in real and thermally activated systems correspond to the Kramers rate. Finding a connection between these two different rates can be the subject of future work. So far, the instantaneous decay rates of trajectories in the TS can be determined exactly using the LMA, and these decay rates are also valid in a certain area close to the TS as has been shown via the ensemble method. In the nearby neighborhood of the NHIM, reactive trajectories correspond to the slowest flux and, hence, are critical in determining the rate of reactions through a time-dependent DS. Since this DS is attached to the NHIM it is, at least in that nearby neighborhood, free of recrossings. The question on how large this “nearby” region needs to be to yield conclusive results is also subject to further work and likely depends on the specific system.

It would also be interesting to determine how reactions can be controlled and what effects can be seen when changing the parameters of the external driving. A different challenge for future research is the application of the methods developed here to specific model chemical reactions. For example, the rank-1 barrier in between the two (meta-) stable conformations in the LiCN ↔ LiNC isomerization reaction may be an easy example to test the methods presented here. An application to KCN^{44} or Ketene^{20} is also conceivable but since these systems contain more than one barrier in between their reactant and product states, the dynamics is expected to be much more complex and chaotic.

## ACKNOWLEDGMENTS

The German portion of this collaborative work was supported by Deutsche Forschungsgemeinschaft (DFG) through Grant No. MA1639/14-1. R.H.’s contribution to this work was supported by the National Science Foundation (NSF) through Grant No. CHE-1700749. M.F. is grateful for support from the Landesgraduiertenförderung of the Land Baden-Württemberg. This collaboration has also benefited from support by the European Union’s Horizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie Grant Agreement No. 734557.