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.

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 + H2,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 expansions12,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 (2d − 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 NHIM32 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.

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 reactions33–37 with pathways far away from the TS, which we do not consider in this paper. We also do not investigate the Kramers rate38 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.

FIG. 1.

Sketch of a reacting trajectory moving close to the TS trajectory (black solid line) in the barrier region. The trajectory changes from reactant (blue dashed-dotted line) to product (dashed red line) when crossing the DS.

FIG. 1.

Sketch of a reacting trajectory moving close to the TS trajectory (black solid line) in the barrier region. The trajectory changes from reactant (blue dashed-dotted line) to product (dashed red line) when crossing the DS.

Close modal

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.

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.

FIG. 2.

Sketch of the stable and unstable manifolds Ws,u in the (x, vx) plane at fixed parameters (y, vy, t) of bath coordinates and time of a multidimensional driven system. The intersection of the two manifolds is a point (x,vx)(NHIM)(y,vy,t) on the time-dependent NHIM. The DS attached to this point is marked by a vertical dashed line. The arrows depict the four reactive and nonreactive regions separated by the manifolds. The red points mark the initial conditions of a reactant trajectory ensemble at time t arranged equidistantly along a line segment parallel to the unstable manifold. The diamonds indicate the ensemble at a later time t + Δt, where this ensemble has partially crossed the DS (blue diamonds).

FIG. 2.

Sketch of the stable and unstable manifolds Ws,u in the (x, vx) plane at fixed parameters (y, vy, t) of bath coordinates and time of a multidimensional driven system. The intersection of the two manifolds is a point (x,vx)(NHIM)(y,vy,t) on the time-dependent NHIM. The DS attached to this point is marked by a vertical dashed line. The arrows depict the four reactive and nonreactive regions separated by the manifolds. The red points mark the initial conditions of a reactant trajectory ensemble at time t arranged equidistantly along a line segment parallel to the unstable manifold. The diamonds indicate the ensemble at a later time t + Δt, where this ensemble has partially crossed the DS (blue diamonds).

Close modal

As highlighted in Fig. 2 by the red dots, an ensemble of Nreact(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.

Further discussion on the NHIM and recrossing-free DSs is found in Ref. 7 and references therein. Algorithms to find the NHIM or to find Ws, as well as the construction of the DS, are presented in Refs. 7, 21, and 29.

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 Nreact(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 Nreact will decrease exponentially fast

Nreact(t)exp(kt),
(1)

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 Nreact(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 Nreact(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,

ddtNreact(t)=k(t)Nreact(t),
(2)

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

k(t)=1Nreact(t)ddtNreact(t)=ddtln(Nreact(t)).
(3)

Due to the way we measure Nreact(t), an exponential decay of reactants also translates into an exponential decay of measurement density for Nreact(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.

FIG. 3.

Rate constant for the TS of the static two-dimensional model system according to Eq. (12) with x^=0. An ensemble of 400 reactive trajectories was prepared with Δx = 10−3 according to Fig. 2 close to the TS and dynamically propagated. Each time one of the trajectories reacts, the number of reactants Nreact(t) is decreased by one. The red dashed line shows a fit according to Eq. (1) to the reaction events, yielding a rate constant of k = 2.7617. The step-function shown as a black solid line as a guide to the eye is to visualize that not only Nreact(t) but also the density of reaction events decreases exponentially fast.

FIG. 3.

Rate constant for the TS of the static two-dimensional model system according to Eq. (12) with x^=0. An ensemble of 400 reactive trajectories was prepared with Δx = 10−3 according to Fig. 2 close to the TS and dynamically propagated. Each time one of the trajectories reacts, the number of reactants Nreact(t) is decreased by one. The red dashed line shows a fit according to Eq. (1) to the reaction events, yielding a rate constant of k = 2.7617. The step-function shown as a black solid line as a guide to the eye is to visualize that not only Nreact(t) but also the density of reaction events decreases exponentially fast.

Close modal

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 Nreact(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.

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

k¯F=μuμs,
(4)

which depends on the Floquet exponents

μu,s=1Tln|mu,s(T)|
(5)

of that trajectory. Here, mu,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, mu,s in Eq. (5) are the two eigenvalues of the (2d × 2d)-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 γ̇=f(γ,t), where γR2d 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

Ṁ=JM,
(6)

with initial conditions M(0) = 12d.

For a Hamiltonian system with time-dependent potential V(q, t) and without friction, Eq. (6) takes the form

Ṁ=0d1d2Vq20dM.
(7)

In Hamiltonian systems, the monodromy matrix M is symplectic, and hence, if λ is an eigenvalue of M, then also 1/λ is an eigenvalue, as well as their complex conjugates λ¯ and 1/λ¯.

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

μu,s=limt1tln|mu,s(t)|.
(8)

In numerical simulations, the time t must be chosen sufficiently large to obtain converged estimates of the Floquet rate k¯F 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.

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 Δxs(t + Δt)/Δxs(t) pulls the line closer to the unstable manifold and a stretching in the direction of the unstable manifold by the factor Δxu(t + Δt)/Δxu(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 Nreact(t). By choosing Δxs(t) = Δxu(t) ≡ Δx as marked in Fig. 2, the number of reactants at time t + Δt is given as

Nreact(t+Δt)=Δxs(t+Δt)Δxu(t+Δt)Nreact(t).
(9)

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 Δxu(t) and Δxs(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

k(t)=ddtlnNreact(t)=ddtlnΔxu(t)lnΔxs(t)=Δẋu(t)Δxu(t)Δẋs(t)Δxs(t).
(10)

Using the linearized equations of motion (7) for the stability matrix M, it can be shown that Δẋs(t)=Δvxs(t) and Δẋu(t)=Δvxu(t). Inserting this into Eq. (10) and also explicitly writing the dependencies on the bath coordinates y and velocities vy, which have been dropped in the above derivations, we finally obtain

k(y,vy,t)=ΔvxuΔxu(y,vy,t)ΔvxsΔxs(y,vy,t).
(11)

This means that for any point on the NHIM parameterized by the bath coordinates y and velocities vy, the instantaneous rate k(y, vy, 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 J 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.

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 xDS(y, vy, t + Δt). Addressing such challenging cases is a subject for future work.

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

V(x,y,t)=Ebexpxx^sinωxt2+ωy22y2πarctan2x2,
(12)

where a time-dependent oscillating Gaussian barrier with height Eb 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(2x). For simplicity, we use dimensionless units in which the parameters are Eb = 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 nN. 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.

FIG. 4.

Floquet rates k¯F obtained for trajectories initialized at t = 0 on the corresponding NHIM of the time-periodically driven potential (12). The rate is visualized as color encoding on the (y, vy)-surface representing the coordinates on the NHIM. The white dots represent a PSOS, where selected trajectories are propagated for many oscillation periods of the potential. Their instantaneous positions y(t = n T) and velocities vy(t = n T) with nN at integer multiples of the potential’s oscillation period T = 2 are marked with these white dots, indicating the stable tori of the dynamics around the TS trajectory, that is located in the center of the tori.

FIG. 4.

Floquet rates k¯F obtained for trajectories initialized at t = 0 on the corresponding NHIM of the time-periodically driven potential (12). The rate is visualized as color encoding on the (y, vy)-surface representing the coordinates on the NHIM. The white dots represent a PSOS, where selected trajectories are propagated for many oscillation periods of the potential. Their instantaneous positions y(t = n T) and velocities vy(t = n T) with nN at integer multiples of the potential’s oscillation period T = 2 are marked with these white dots, indicating the stable tori of the dynamics around the TS trajectory, that is located in the center of the tori.

Close modal

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 Nreact(t), initialized with 200 reacting trajectories at t0 = 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 Nreact(t) are shown as straight lines, yielding the instantaneous rates k(t) at the respective time. As expected, the decay Nreact(t) is not purely exponential in the driven case.

FIG. 5.

(a) Reactant populations Nreact(t) when propagating an ensemble of 200 reactive trajectories. The trajectories are initialized at t0 = 0 with Δx = 10−3 according to the method explained in Fig. 2 close to the TS trajectory of the two-dimensional model system (12). The straight lines yield k(t) as instantaneous derivatives of Nreact(t), taken at several times t and visualized by crosses in the respective color. (b) Instantaneous ensemble rate k(t) (thin black line), obtained according to Sec. II A by numerical differentiation of Nreact(t) for ensembles initialized close to the TS trajectory of the two-dimensional model system (12) at different times t0. Here, T = 2 is the period of the TS trajectory. The rate kM(t), calculated via the LMA, as described in Sec. II C, is given as a thin, red dashed line that coincides perfectly with k(t), hence yielding the same mean value k¯M=k¯ when averaged over a full period of the oscillating barrier. The Floquet rate, k¯F=3.806 (thick, red dashed line) obtained by a linear regression of the red line in Fig. 7 is in perfect accordance to the mean values k¯ (thick black line) of both k(t) and kM(t).

FIG. 5.

(a) Reactant populations Nreact(t) when propagating an ensemble of 200 reactive trajectories. The trajectories are initialized at t0 = 0 with Δx = 10−3 according to the method explained in Fig. 2 close to the TS trajectory of the two-dimensional model system (12). The straight lines yield k(t) as instantaneous derivatives of Nreact(t), taken at several times t and visualized by crosses in the respective color. (b) Instantaneous ensemble rate k(t) (thin black line), obtained according to Sec. II A by numerical differentiation of Nreact(t) for ensembles initialized close to the TS trajectory of the two-dimensional model system (12) at different times t0. Here, T = 2 is the period of the TS trajectory. The rate kM(t), calculated via the LMA, as described in Sec. II C, is given as a thin, red dashed line that coincides perfectly with k(t), hence yielding the same mean value k¯M=k¯ when averaged over a full period of the oscillating barrier. The Floquet rate, k¯F=3.806 (thick, red dashed line) obtained by a linear regression of the red line in Fig. 7 is in perfect accordance to the mean values k¯ (thick black line) of both k(t) and kM(t).

Close modal

The instantaneous rate k(t), which is given as the most accurate result from the numerical differentiation of the reactant decays Nreact(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 kmin ≈ 3.4 and kmax ≈ 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 kM(t) is shown as a thin red dashed line in Fig. 5(b). The perfect agreement between k(t) and kM(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¯F can also be obtained. It is included in Fig. 5(b) as a thick red dashed horizontal line. Since the Floquet rate k¯F 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 kM(t) seem to oscillate around k¯F. For a more detailed comparison, the mean value of the instantaneous rates k¯=k(t)¯=kM(t)¯ is included as a thick black horizontal line in Fig. 5(b). The almost perfect agreement between k¯=3.8067 and k¯F=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.

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 t0 ∈ [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 Nreact(t) of these reactive ensembles—cf. the thin black line in Fig. 6(a).

FIG. 6.

Instantaneous rate k(t) (thin black line), obtained using a numerical differentiation of the reactant population Nreact(t) of 80 ensembles with 200 trajectories each, launched equidistantly in time in between a time interval of Δt = 40. All ensembles are prepared according to Sec. II A 1 with Δx = 10−3 close to a trajectory that is initialized for t = 0 at position y = 0.4 and vy = −0.75 (above) or y = 0.7 and vy = −0.75 (below) on the NHIM of the time-dependently driven system defined in Eq. (12). According to the procedure explained in Sec. II C, kM(t) is obtained and displayed as a thin, red dashed line that perfectly matches k(t). The Floquet rate, calculated as explained in Sec. II C or as displayed in Fig. 7, is given as a thick, red dashed line. When averaging over three quasioscillations of k(t), visualized by the two vertical, blue dashed lines, the Floquet rate perfectly matches the mean ensemble rate k¯ of both k(t) and kM(t).

FIG. 6.

Instantaneous rate k(t) (thin black line), obtained using a numerical differentiation of the reactant population Nreact(t) of 80 ensembles with 200 trajectories each, launched equidistantly in time in between a time interval of Δt = 40. All ensembles are prepared according to Sec. II A 1 with Δx = 10−3 close to a trajectory that is initialized for t = 0 at position y = 0.4 and vy = −0.75 (above) or y = 0.7 and vy = −0.75 (below) on the NHIM of the time-dependently driven system defined in Eq. (12). According to the procedure explained in Sec. II C, kM(t) is obtained and displayed as a thin, red dashed line that perfectly matches k(t). The Floquet rate, calculated as explained in Sec. II C or as displayed in Fig. 7, is given as a thick, red dashed line. When averaging over three quasioscillations of k(t), visualized by the two vertical, blue dashed lines, the Floquet rate perfectly matches the mean ensemble rate k¯ of both k(t) and kM(t).

Close modal

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 kM(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 kM(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¯F=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.

FIG. 7.

(a) Difference ln |mu| − ln |ms|, where mu,s are the eigenvalues of the time-dependent monodromy matrix along the unstable and stable directions of the driven barrier according to Eq. (12). The monodromy matrix is obtained time-dependently according to Eq. (6) for the periodic TS trajectory of Fig. 5 (red line), as well as for the two nonperiodic trajectories introduced in Fig. 6 (black dotted line and blue dashed-dotted line). The Floquet rates k¯F in the limit t → ∞ are obtained according to Eq. (8) by a linear regression to these curves. (b) Black line: calculated Floquet rates for trajectories initialized on a cross section at the velocity vy of the TS trajectory of the NHIM at time t = 0. This corresponds to a horizontal line through the center of the tori in Fig. 4. To verify these results, the mean ensemble rates (red crosses) are obtained according to Sec. II A for trajectories initialized at several values y on this cross section. The colored circles mark the corresponding Floquet calculations of Fig. 7(a).

FIG. 7.

(a) Difference ln |mu| − ln |ms|, where mu,s are the eigenvalues of the time-dependent monodromy matrix along the unstable and stable directions of the driven barrier according to Eq. (12). The monodromy matrix is obtained time-dependently according to Eq. (6) for the periodic TS trajectory of Fig. 5 (red line), as well as for the two nonperiodic trajectories introduced in Fig. 6 (black dotted line and blue dashed-dotted line). The Floquet rates k¯F in the limit t → ∞ are obtained according to Eq. (8) by a linear regression to these curves. (b) Black line: calculated Floquet rates for trajectories initialized on a cross section at the velocity vy of the TS trajectory of the NHIM at time t = 0. This corresponds to a horizontal line through the center of the tori in Fig. 4. To verify these results, the mean ensemble rates (red crosses) are obtained according to Sec. II A for trajectories initialized at several values y on this cross section. The colored circles mark the corresponding Floquet calculations of Fig. 7(a).

Close modal

Averaging over the three quasioscillations of k(t) or kM(t), indicated by the blue dashed vertical lines in Fig. 6(a), yields the average rate k¯=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¯F=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 kM(t) but also the accordance between their average rate k¯=4.255 and the Floquet rate k¯F=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.

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 KCN44 or Ketene20 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.

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.

1.
H.
Eyring
,
J. Chem. Phys.
3
,
107
(
1935
).
2.
E. P.
Wigner
,
J. Chem. Phys.
5
,
720
(
1937
).
3.
P.
Pechukas
,
Annu. Rev. Phys. Chem.
32
,
159
(
1981
).
4.
D. G.
Truhlar
,
B. C.
Garrett
, and
S. J.
Klippenstein
,
J. Phys. Chem.
100
,
12771
(
1996
).
5.
R. G.
Mullen
,
J.-E.
Shea
, and
B.
Peters
,
J. Chem. Phys.
140
,
041104
(
2014
).
6.
S.
Wiggins
,
Regular Chaotic Dyn.
21
,
621
(
2016
).
7.
M.
Feldmaier
,
P.
Schraft
,
R.
Bardakcioglu
,
J.
Reiff
,
M.
Lober
,
M.
Tschöpe
,
A.
Junginger
,
J.
Main
,
T.
Bartsch
, and
R.
Hernandez
,
J. Phys. Chem. B
123
,
2070
(
2019
).
8.
J. C.
Lorquet
,
J. Chem. Phys.
146
,
134310
(
2017
).
9.
S.
Patra
and
S.
Keshavamurthy
,
Phys. Chem. Chem. Phys.
20
,
4970
(
2018
).
10.
V.
Krajňák
and
H.
Waalkens
,
J. Math. Chem.
56
,
2341
2378
(
2018
).
11.
Y.
Tamiya
,
R.
Watanabe
,
H.
Noji
,
C.-B.
Li
, and
T.
Komatsuzaki
,
Phys. Chem. Chem. Phys.
20
,
1872
(
2018
).
12.
R.
Hernandez
and
W. H.
Miller
,
Chem. Phys. Lett.
214
,
129
(
1993
).
13.
S.
Wiggins
,
L.
Wiesenfeld
,
C.
Jaffe
, and
T.
Uzer
,
Phys. Rev. Lett.
86
,
5478
(
2001
).
14.
T.
Uzer
,
C.
Jaffé
,
J.
Palacián
,
P.
Yanguas
, and
S.
Wiggins
,
Nonlinearity
15
,
957
(
2002
).
15.
T.
Bartsch
,
R.
Hernandez
, and
T.
Uzer
,
Phys. Rev. Lett.
95
,
058301
(
2005
).
16.
T.
Bartsch
,
J. M.
Moix
,
R.
Hernandez
,
S.
Kawai
, and
T.
Uzer
,
Adv. Chem. Phys.
140
,
191
(
2008
).
17.
R.
Hernandez
,
J. Chem. Phys.
101
,
9534
(
1994
).
18.
A.
Allahem
and
T.
Bartsch
,
J. Chem. Phys.
137
,
214310
(
2012
).
19.
A.
Junginger
,
P. L.
Garcia-Muller
,
F.
Borondo
,
R. M.
Benito
, and
R.
Hernandez
,
J. Chem. Phys.
144
,
024104
(
2016
).
20.
G. T.
Craven
and
R.
Hernandez
,
Phys. Chem. Chem. Phys.
18
,
4008
(
2016
).
21.
M.
Feldmaier
,
A.
Junginger
,
J.
Main
,
G.
Wunner
, and
R.
Hernandez
,
Chem. Phys. Lett.
687
,
194
(
2017
).
22.
E.
Pollak
and
P.
Pechukas
,
J. Chem. Phys.
69
,
1218
(
1978
).
23.
P.
Pechukas
and
E.
Pollak
,
J. Chem. Phys.
71
,
2062
(
1979
).
24.
C.
Jaffé
,
S. D.
Ross
,
M. W.
Lo
,
J.
Marsden
,
D.
Farrelly
, and
T.
Uzer
,
Phys. Rev. Lett.
89
,
011101
(
2002
).
25.
H.
Teramoto
,
M.
Toda
, and
T.
Komatsuzaki
,
Phys. Rev. Lett.
106
,
054101
(
2011
).
26.
C.-B.
Li
,
A.
Shoujiguchi
,
M.
Toda
, and
T.
Komatsuzaki
,
Phys. Rev. Lett.
97
,
028302
(
2006
).
27.
H.
Waalkens
and
S.
Wiggins
,
J. Phys. A: Math. Gen.
37
,
L435
(
2004
).
28.
U.
Çiftçi
and
H.
Waalkens
,
Phys. Rev. Lett.
110
,
233201
(
2013
).
29.
R.
Bardakcioglu
,
A.
Junginger
,
M.
Feldmaier
,
J.
Main
, and
R.
Hernandez
,
Phys. Rev. E
98
,
032204
(
2018
).
30.
P.
Schraft
,
A.
Junginger
,
M.
Feldmaier
,
R.
Bardakcioglu
,
J.
Main
,
G.
Wunner
, and
R.
Hernandez
,
Phys. Rev. E
97
,
042309
(
2018
).
31.
G. T.
Craven
,
T.
Bartsch
, and
R.
Hernandez
,
J. Chem. Phys.
141
,
041106
(
2014
).
32.
M.
Tschöpe
,
M.
Feldmaier
,
J.
Main
, and
R.
Hernandez
, “
Neural network approach for the dynamics on the normally hyperbolic invariant manifold of periodically driven systems
,”
Phys. Rev. E
(submitted).
33.
D.
Townsend
,
S. A.
Lahankar
,
S. K.
Lee
,
S. D.
Chambreau
,
A. G.
Suits
,
X.
Zhang
,
J. L.
Rheinecker
,
L. B.
Harding
, and
J. M.
Bowman
,
Science
306
,
1158
(
2004
).
34.
I. S.
Ulusoy
,
J. F.
Stanton
, and
R.
Hernandez
,
J. Phys. Chem. A
117
,
10567
(
2013
).
35.
F. A. L.
Maugière
,
P.
Collins
,
G.
Ezra
,
S. C.
Farantos
, and
S.
Wiggins
,
Chem. Phys. Lett.
592
,
282
(
2014
).
36.
J. M.
Bowman
and
A. G.
Suits
,
Phys. Today
64
(
11
),
33
(
2011
).
37.
38.
P.
Hänggi
,
P.
Talkner
, and
M.
Borkovec
,
Rev. Mod. Phys.
62
,
251
(
1990
) and references therein.
39.
G. T.
Craven
,
T.
Bartsch
, and
R.
Hernandez
,
J. Chem. Phys.
142
,
074108
(
2015
).
40.
S. K.
Upadhyay
,
Chemical Kinetics and Reaction Dynamics
(
Springer Netherlands
,
2006
).
41.
K. A.
Connors
,
Chemical Kinetics: The Study of Reaction Rates in Solution
(
Wiley-VCH
,
1990
).
42.
F.
Revuelta
,
G. T.
Craven
,
T.
Bartsch
,
F.
Borondo
,
R. M.
Benito
, and
R.
Hernandez
,
J. Chem. Phys.
147
,
074104
(
2017
).
43.
A.
Junginger
,
L.
Duvenbeck
,
M.
Feldmaier
,
J.
Main
,
G.
Wunner
, and
R.
Hernandez
,
J. Chem. Phys.
147
,
064101
(
2017
).
44.
H.
Párraga
,
F. J.
Arranz
,
R. M.
Benito
, and
F.
Borondo
,
J. Phys. Chem. A
122
,
3433
(
2018
).