Transition state theory for activated systems with driven anharmonic barriers

The external control of reactions^1–7 is increasingly a focus in the chemical sciences because novel products, synthesized at improved yields and selectivity, are needed for advanced industrial applications. Catalysts are often used to accelerate reaction rates and increase the efficiency of industrial processes. Unfortunately, most catalysts can only be used to increase the rate of chosen reactions; they do not allow a finer control of the reaction pathway or outcome.

Several different strategies have been developed to control chemical reactions in a more flexible and accurate manner. For example, mechanical forces have been used to control the height of the activation barrier.⁸ In ultracold systems, reaction rates have been seen to be strongly influenced by the internal state of the reactants.⁹ However, a much more general approach towards chemical control is available through lasers. The amazing recent developments of this technology have opened the route for a new, and much more accurate, control of chemical reactivity that was simply inconceivable a short time ago.¹⁰ Such mode-selective chemistry permits carrying a system to a desired product state^3,7,11,12 under the influence of controlled laser pulses. This huge variety of laser-driven systems all have in common that the energetic barrier dividing the phase space of the system into distinct regions, reactants and products, is time-dependent. Thus, the phase space bottleneck formed at the barrier top, which is the limiting step for the reaction to take place, becomes also time-dependent and, as a consequence, traditional transition state theory (TST)^13–20 fails to give accurate rate predictions.

For fixed barriers, TST describes the reaction mechanism approximately by specifying the rate-limiting bottleneck in phase space. It allows the identification of reactive trajectories and leads to increased accuracy in the computation of reaction rates. The fundamental challenge to the implementation of TST is to remove the approximation by constructing an optimal dividing surface (DS) that is crossed once and only once by all reactive trajectories.^21–23 A calculation of the flux through this DS gives the reaction rates. However, if the DS is not strictly recrossing-free, the reactive flux is overestimated and TST gives an upper bound of the true reaction rate. The discoveries of periodic orbit dividing surfaces (PODS)^24–27 in systems with two degrees of freedom and normally hyperbolic invariant manifolds (NHIMs) in higher dimension^19,28–40 have allowed the computation of DSs that are free of recrossings. Unfortunately, these constructions only apply for autonomous Hamiltonian systems that describe gas phase reactions. Variational TST (VTST) optimizes the choice of DS by minimizing the flux.⁴¹ Both TST and VTST can nevertheless give large overestimates to the true reaction rate when the system is strongly coupled to an environment such as a solvent.

A strictly recrossing-free DS can be constructed in driven reactions with harmonic barriers using the so-called transition state (TS) trajectory.⁴² The TS trajectory is the only solution of the equations of motion (EoM) that remains close to the energetic saddle point for all times. This optimal DS can be obtained for anharmonic barriers using a perturbative scheme.^43–45 Time-dependent manifolds attached to the hyperbolic TS trajectory separate the phase space into regions of selected reactivity. The determination of these manifolds in driven activated systems^11,46–53 would allow one to enhance or suppress the formation of specific products through tailored control pulses, and that is the goal of this work.

In reactions with time-dependent energy barriers, a fixed DS placed close to the barrier top leads to a large number of recrossings and is therefore ill-suited for describing the reaction dynamics. Recently, Craven et al.^54–56 have shown how a time-dependent DS can be constructed for dissipative systems over barriers that are moving periodically in time. In an athermal periodically driven system, the TS trajectory is a periodic orbit that moves close to the barrier top region, following the motion determined by the external driving field. Furthermore, the asymptotic reaction rate is given by the difference of the Floquet exponents of this distinguished trajectory.

In this paper, we extend the previous work reported in Refs. 54–56 to include the treatment of a reaction occurring in a Markovian solvent by constructing the exact TS trajectory over an anharmonic energy barrier that is moving in time. Using this trajectory, we can identify a strictly recrossing-free DS for an anharmonic system in the presence of noise. Likewise, we have succeeded in classifying trajectories as reactive or nonreactive without the necessity of any numerical simulation by comparing their initial velocity to the critical value that they must exceed in order to react. We also demonstrate that the rate of a chemical reaction in the presence of these external forces is given by the difference of the characteristic Lyapunov exponents of the TS trajectory. This is directly analogous to, but extends, our earlier finding⁵⁶ that the rates are given by the differences of the characteristic Floquet exponents of the TS trajectory in the periodically driven case.

The outline of this paper is as follows: In Sec. II, we describe the model that will serve as the paradigm for the present study: the Langevin equation on a time-dependent anharmonic potential. The method for the calculation of the TS trajectory and the associated recrossing-free DS for this system is developed in Sec. III. Expressions for rate constants are determined in Sec. IV in terms of the Lyapunov exponents of the TS trajectory. The numerical results of these calculations are reported for various driving forms and thermal strengths in Sec. V. Both of the averaging procedures for the calculation of reaction rates described earlier are shown to agree in periodically driven systems, in systems driven by two incommensurable frequencies, and when the system strongly interacts with its environment through a stochastic force. The perturbative calculation of reaction rates also provides good agreement over a range of parameter values. A key finding is that all of these methods also support the conjecture that the rates are determined by the stability of the TS trajectory as quantified through the difference of the Lyapunov exponents.

The activated dynamics of a collective reaction coordinate on a time-varying energy landscape can be described by a Langevin equation^57,58

where x(t) is the mass-scaled position coordinate, $γ$ is the friction, $ξ α ( t )$ is the fluctuating force exerted by the environment, and U(x, t) is a time-dependent potential barrier. As in previous work, the driving is induced by a horizontal movement of the barrier at the instantaneous mass-scaled position X(t). For simplicity, we retain the symmetry of the barrier and introduce anharmonicity through a quartic term. The potential is taken to be

where $ω b$ is the inverse barrier frequency and $𝜖$ represents the anharmonicity. Except where noted, we consider sinusoidal driving through

such that the barrier motion is periodic in time.

The strength of the thermal noise in the Langevin equation (1) is varied through the parameter $σ$⁠. In equilibrium at temperature T, it is related to the friction by the fluctuation-dissipation theorem,

In some cases, detailed below, it is convenient to relax this restriction and regard $σ$ as an independent parameter. We consider only Markovian solvent-solute interactions, although history-dependent effects can be included.^45,59,60 The noise is therefore white (uncorrelated in time) such that

and

where $… α$ denotes an average over the realizations $α$ of the noise.

The Langevin equation (1) is a second-order differential equation that can be rewritten as a system of two first-order differential equations. It takes the form

where we have specified the potential through Eq. (2). Example trajectories evolved through Eq. (7) are shown in Fig. 1 for a dissipative environment with (⁠$σ=3$⁠) and without (⁠$σ=0$⁠) noise. Dimensionless units, obtained by setting a = 1 and $ω b =1$⁠, are employed without loss of generality. Throughout the paper, the temperature of the thermal bath is taken as $k B T=1$⁠, and the initial phase of the periodic driving is $φ=0$⁠. In the noiseless case, the TS trajectory, starting at $x ‡ ( 0 ) =0$⁠, is a periodic orbit with the period of the driving as seen in the top panel of Fig. 1. In the noisy case, the position at time 0 and evolution of the TS trajectory depend on the past and future driving of the system. For the particular noise sequence randomly selected in creating the bottom panel of Fig. 1, the position $x ‡ ( 0 )$ of the TS trajectory at time 0 is found to be near 1.3.

We first consider the harmonic limit of Eq. (7) for a stationary barrier, a dissipative environment, and no noise. The last of these conditions effectively eliminates the temperature, and hence the system can be called athermal. The solution of the resulting autonomous EoM can be obtained by changing to the diagonal coordinates defined by

with eigenvalues

The inverse transformation is given by

In the coordinate system defined by u and s, Eq. (7) transforms into

where

Equation (11) has a unique solution $u ‡$ and $s ‡$ that remains bounded as $t→±∞$⁠. This solution is called the TS trajectory.

In the original coordinates (x and $v$⁠), the TS trajectory satisfies the integral equations

where f is given by Eq. (12), and S denotes the S-functional⁶¹ defined as

The bounded TS trajectory can be obtained directly from these functionals as they suppress the exponential contribution of the solution of the EoM, returning only the part that remains in the vicinity of the barrier top for all times.

In the harmonic limit (⁠$𝜖=0$⁠), the TS trajectory of Eq. (13) can be separated into two parts,

where $x α ‡ ( 0 )$ is a stochastic term that depends on the noise,

and $x T ‡ ( 0 ) ( t )$ is the deterministic and periodic part^54–56 that depends on the external periodic driving as

with

In the anharmonic setting (⁠$𝜖>0$⁠), the TS trajectory $x ‡ ( t )$ appears in the argument of the S-functionals on the right hand side of Eqs. (13) and (14). As a consequence, these relations do not give an explicit expression for the TS trajectory. Two different methods have been developed in order to obtain the TS trajectory in this case. The first method is based on perturbation theory (PT). It is described in Sec. III A. The second, presented in Sec. III B, is based on a numerical iteration of Eqs. (13) and (14) and yields the numerically exact TS trajectory.

The perturbative calculation of the TS trajectory starts by expanding it in powers of $𝜖$⁠, the perturbative parameter that accounts for the anharmonic strength in the potential [cf. Eq. (2)], as

where $x ‡ ( 0 ) ( t )$ is given by Eq. (16). The TS trajectory is calculated order-by-order in $𝜖$⁠, starting with the lowest order term and continuing up to the desired order. As usual, this procedure is accomplished by matching the terms in the EoM that have the same order in $𝜖$⁠. In the case of a noiseless system (⁠$σ=0$⁠), the lowest order correction to the harmonic TS trajectory (16) is given by

with the coefficients

and $b= γ 2 + Ω 2$⁠. The explicit expression for the second order correction term $x ‡ ( 2 ) ( t )$ is contained in the Mathematica notebook that is provided as supplementary material.

The top panel of Fig. 2 compares the time evolution of the TS trajectory for a system without noise (⁠$σ=0$⁠) computed using PT and the method reported in Sec. III B. The agreement between the PT at order 3 and the numerical results is quite remarkable in view of the large value of the expansion parameter (⁠$𝜖=10$⁠). Unfortunately, in the presence of noise (⁠$σ>0$⁠), the accuracy of the PT is restricted to much smaller values of the perturbative parameter. An example of a noisy TS trajectory is shown in the bottom panel of Fig. 2 for $𝜖=0.1$⁠, a value where the PT still holds. The reference (exact) TS trajectory has been computed using the method reported below in Sec. III B. The breakdown of PT occurs for smaller values of the perturbative parameter in the presence of noise. The noise causes the TS trajectory to deviate much farther from the barrier top than what is seen in deterministic driving (at the given parameter values) and leads it to explore more of the anharmonic regions of the barrier. As a consequence, the higher order perturbative expansion coefficients $x ‡ ( n )$⁠, which contribute to the TS trajectory with an amplitude of order $𝜖 n$⁠, are larger in the noisy system. For the noiseless system, $x ‡ ( n + 1 )$ is, in general, one order of magnitude smaller than $x ‡ ( n )$⁠, while for the noisy case it can even be larger than $x ‡ ( n )$⁠. Of course the magnitude of $x ‡ ( n )$ will shrink for smaller values of $σ$⁠, and the range of validity of PT will correspondingly get wider.

In order to be able to calculate the noisy TS trajectory for large values of $𝜖$⁠, we have also developed another, purely numerical procedure. For this purpose, one must first compute the harmonic TS trajectory $x 0 ‡ ≡ x ‡ ( 0 )$ given by the sum of Eqs. (17) and (18). This TS trajectory can be substituted into the right-hand side of Eq. (13) to obtain an improved approximation to $x ‡$ for a non-zero value of $𝜖$⁠. For sufficiently small anharmonicity $𝜖$⁠, iterating through the order of the approximation converges to the true TS trajectory.

If $𝜖$ is too large to achieve convergence directly, we can add another layer of iteration over the value of $𝜖$⁠. One first computes the TS trajectory for a small value $𝜖=Δ𝜖≪1$⁠. One then uses $x Δ 𝜖 ‡$ as an initial approximation in the same manner in order to calculate $x 2 Δ 𝜖 ‡$⁠, and the iteration is continued up to the $𝜖$ value of interest. This algorithm can be used to compute the TS trajectory for any form of external driving. It is summarized in pseudocode as follows:

1. procedure TS-TRAJECTORY

2. $$Compute the harmonic TS trajectory $x 0 ‡$

3. $$Choose a small increment $Δ𝜖$ of $𝜖$

4. $$Set l = 1

5. $$while $l×Δ𝜖≤𝜖$ do

   • 5a

     Compute iteratively $x l Δ 𝜖 ‡$ using Eqs. (13) and (14) up to the required degree of precision

   • 5b

     l = l + 1

6. $$return $x 𝜖 ‡$

7. end procedure

If convergence is not achieved in step 6, then the step size $Δ𝜖$ must be reduced. All our computations have been performed using a step size equal to $Δ𝜖=0.1$ and a required degree of precision in the computation of the transition state trajectory equal to 10⁻¹⁰. For this set of parameters, the number of iterations needed for convergence was typically less than 100.

Figure 3 shows the time evolution of the TS trajectory computed for various friction values and noise strengths given by the fluctuation-dissipation theorem (4). The top left panel shows the results for $γ=0$⁠, which is associated with a Hamiltonian system. In this case, large values of the anharmonic parameter $𝜖$ generally correlate with large amplitudes. The value of $𝜖$ has no influence on the period of the TS trajectory, which is always the same as the external driving field $T=2π∕Ω$⁠. When the system is dissipative (⁠$γ>0$⁠), the TS trajectory is no longer periodic due to the stochastic interaction with the environment. Figure 3 also presents the TS trajectories computed for $γ∈ 0.1 , 1,3$⁠. Notice that the larger the dissipation, the more violently the TS trajectory oscillates due to stronger coupling with the environment. The time-averaged amplitude of the TS trajectory shown in Fig. 4 as a function of the anharmonic parameter $𝜖$ provides a quantitative and qualitative estimate of the magnitude of the oscillations. The amplitude always increases with the friction constant $γ$ and decreases with increasing frequency $Ω$⁠. The dependence of the amplitude on the anharmonic parameter is more difficult to analyze because it depends strongly on both parameters. For example, while for $γ=0.1$ and $Ω=1$⁠, the amplitude increases linearly with $𝜖$⁠, for $γ=3$ and $Ω=10$⁠, it always decreases in the considered $𝜖$ range. In general, however, the time-averaged amplitude first decreases with $𝜖$⁠, and then increases as can be observed, for example, in the curve for $γ=1$ and $Ω=3$⁠. It is notable that the onset of the irregularity of the anharmonic trajectories occurs when the friction $γ$ is on the order of the inverse period of the periodic driving, $1∕T=2π∕Ω≈0.5$⁠.

The method outlined here can be also applied to the computation of the TS trajectory of a system driven by several frequencies, even if they are incommensurable. For this purpose, one needs only to substitute the appropriate function X(t) into Eq. (12). The two top panels of Fig. 5 show the time evolution of the aperiodic pulse

for $Ω 1 =1$ (left) and $Ω 1 =7$ (right). The TS trajectories associated with these two pulses appear in the two bottom panels. As for the noiseless system shown in Fig. 3, the oscillation amplitude of the TS trajectory increases with the anharmonic parameter $𝜖$⁠. The TS trajectory is the only trajectory that remains close to the barrier top for all times and never falls down to the reactant or the product regions (cf. Figs. 3 and 5). The time evolution in Fig. 1 of a swarm of trajectories for a system in athermal (top) and thermal (bottom) environments exhibits this effect. The TS trajectory, shown in white, is the only solution of the EoM that remains close to the barrier top. All other trajectories fall on one of the sides of the DS after a certain time. Specifically, the trajectories that have an initial velocity $v$(0) larger than the critical velocity, $v ‡ ( 0 ) = V ‡$⁠, fall into the product well defined by $x>0$⁠, while the trajectories with a velocity $v ( 0 ) < V ‡$ fall into the reactant well defined by $x<0$⁠. The trajectories with an initial velocity $v$(0) close to $V ‡$ take longer to decay into the reactant of product states as they remain close to the TS trajectory for long times (Fig. 6).

The rate constant k of a chemical reaction is a measure of the speed at which reactants (R) are transformed into products (P). On complex energy landscapes, the calculation of a rate constant can be computationally very demanding as it usually requires the propagation of a very large number of trajectories to get statistically significant results. For this purpose, one usually takes initial conditions in the reactant well, with some distribution in the position and velocity, and solves the EoM of the system while following the survival probability of each trajectory. In equilibrated thermal systems, the velocity distribution p(⁠$v$⁠) takes the Boltzmann form

In the calculations presented below, we prepare the system with this distribution in order to model a system that is initially at thermal equilibrium. As a reference for modeling nonequilibrium states, we will also use a uniform distribution $U$ centered on the stable manifold of the TS trajectory.

A trajectory is reactive if it crosses the barrier and remains in the product region, and it is nonreactive otherwise. The population of the product state at time t is

where $x ‡$ is the position of the time-dependent DS that is attached to the TS trajectory (see Sec. III), p_x(x, t) is the distribution of x at time t, and $Θ ( x )$ is the Heaviside step function,

It has previously been demonstrated that given some initial position x = x₀, reactive trajectories at this position are those that have an initial velocity that exceeds a critical value.^{42–44,54–56} This critical velocity $V ‡$ is defined by the intersection of the stable manifold with the line x = x₀ of the trajectories’ initial positions. Reactive trajectories have an initial velocity $v ( 0 ) = v 0 ≥ V ‡$⁠, while nonreactive trajectories are those defined by $v 0 < V ‡$⁠. Consequently, in a rate calculation, one only has to account for those trajectories with $v 0 ≥ V ‡$⁠. Thus given a state prepared with initial distribution in position $δ ( x − x 0 )$⁠, the asymptotic (⁠$t→∞$⁠) population for any bounded distribution is

For a reaction with one reactant state and one product state, the corresponding probabilities are related by $P R ( ∞ ) =1− P P ( ∞ )$⁠. The rate of a reaction can be obtained from the decay of the population. In transient times, fluctuations due to the barrier motion and the thermal driving⁵⁶ manifest as nonstatistical effects in the decay rate.^62–64 However, if a rate exists, at long times this decay will settle to a constant value that is independent of the initial distribution.⁶² 

A computation of rate constants based on the numerical simulation of a large number of reactive trajectories is very time consuming. However, the correct identification of the geometrical structures that determine the dynamics of the reaction in the vicinity of the barrier top provides an alternative approach for this task that circumvents this problem: the stability of these guiding phase space objects^65–70 has been shown to determine the rate of state transitions.^55,56,71

The asymptotic decay rate is determined by those trajectories close to the stable manifold, whose dynamics are described by the linear equations of motion, and the corresponding stability matrix M(t) as

where $Δx=x− x ‡$ and $Δv=v= v ‡$ are the relative coordinates. In what follows, it is also useful to define the relative coordinates $Δu=u− u ‡$ and $Δs=s− s ‡$ in the diagonal space. The stability matrix accounts for the linearized motion about the TS trajectory and fulfills

where J is the Jacobian matrix

that satisfies

For the case of a periodically driven, noiseless system with period T, the stability matrix is a monodromy matrix. In this case, the eigenvalues of M(T) are the Floquet multipliers m_u,s. In systems driven only by periodic forces, the rate of reaction is given by the difference between the corresponding Floquet exponents, $μ u , s =log| m u , s |∕T$⁠, of the TS trajectory.^54–56 

A conjecture, supported by the observed persistence of this rate theory over complex periodic driving forms,⁵⁶ is that the previous findings can be generalized to thermally driven (noisy) systems. In the case of aperiodic driving, the stability matrix M(t) is no longer periodic and it must be computed for very long times. The corresponding eigenvalues can then be used to calculate the Lyapunov exponents of the TS trajectory,

and their difference, we posit, once again gives the reaction rate. For noisy systems, we will assume that $μ s <0< μ u$ as in the harmonic limit of the noiseless system, where $0< μ u = λ u <1$ and $μ s = λ s <0$⁠. Notice that this condition could in principle break down if the anharmonicities become too strong. However, this case never occurred in our simulations; we always found that both Lyapunov exponents were real and of opposite sign. In order to reduce the numerical errors, we compute the stable Lyapunov exponent, $μ s$⁠, by making a backward time evolution changing t to −t in Eq. (31), as $μ s$ is then associated with an unstable direction in the phase space (see below).

The eigenvectors $v u ( t )$ and $v s ( t )$ of the M(t) matrix define two very particular directions to study the dynamics in the vicinity of the TS trajectory: A trajectory displaced in the direction of the unstable vector $v u$ will typically separate from the TS trajectory with an expansion factor $e μ u t$⁠. Thus, this eigenvector defines an unstable direction in the phase space. On the other hand, a trajectory with an initial condition placed at the stable vector, $v s$⁠, will typically approach the TS trajectory with a contraction factor $e μ s t$⁠. Consequently, it defines a stable direction in the phase space. The vectors $v u$ and $v s$ determine a new set of coordinates, $Δu$ and $Δs$⁠, that provide the natural framework to study the dynamics close to the TS trajectory.

Let us now consider the decay of a set of trajectories, each with initial position x(0) = x₀. In the long time limit, the reactive flux of this set, and hence the reaction rate, is determined by the trajectories with initial velocities close to the stable manifold. We further assume that the distance from this manifold is small enough that the linear approximation given by Eq. (29) satisfactorily describes the dynamics close to the TS trajectory. In analogy to the case of periodic driving,^54–56 we thus conjecture that the number of trajectories that cross the moving DS is proportional to $e ( μ u − μ s ) t$⁠. Thus, the reactive flux, which is the time derivative of the population decay, is proportional to the same factor. It follows from this that the reaction rate is given by the difference in the unstable and stable Lyapunov exponents of the TS trajectory,

In Sec. V, we explore the validity of Eq. (32) by comparing the rates given by stability analysis with those obtained directly using numerical simulation.

To measure the reaction rates $k f$ of the system modeled by Eq. (7), the survival probability of a large number (N = 10⁸ to 10⁹) of trajectories with initial positions x_i(0) in the reactant region was followed as a function of time. The reactant population is calculated as

where the Heaviside function $Θ$⁠, defined in Eq. (25), distinguishes between reactive and nonreactive trajectories using a DS attached to the TS trajectory calculated in Sec. III. The reaction rate is obtained from the decay of $P R ( t )$⁠. After the initial transient, the decay of the scaled logarithmic population $−ln P R ( t ) − P R ( ∞ )$ is approximately linear in time. The first-order reaction rate, as calculated from numerical simulation, is then given by the slope of this line. In all cases, we have found this DS free of recrossings.

The inclusion of noise gives rise to fluctuations in the population decay over time-intervals where our computational resources allowed sampling. These fluctuations depend strongly on the realization of the driving and noise, giving rise to a high degree of variation in the measured rates. Thus, an averaging method was implemented to obtain the overall reaction rate: the decay of a statistical sampling of populations is first obtained, rates are determined at different initial times, and all of the rates are averaged to obtain an overall rate of decay. In more detail, the decay of the population at t = t₀ is first followed to obtain the rate k(t₀), the decay of the population beginning at $t= t 1 > t 0$ is then followed to obtain k(t₁), and so forth at each t_i. The initial times t_i are evenly spaced with $Δt= t i − t j =5$ and chosen sufficiently far apart such that the transient section of the population decay leading to the rate k(t_i) is complete before t_i+1. Consequently, the different values k(t_i) represent the influence of independent noise sequences. The overall rate we report is the average of all these rates,

which is taken over the initial times. The total number N of trajectories is equally distributed over these time points. At each time point t_i, the initial position of each trajectory was taken the same distance from the TS trajectory, i.e., $x 0 = x ‡ ( t i ) −0.1$⁠, with a Boltzmann distribution in velocity, although this rate is independent of the chosen distribution⁶² as we previously confirmed.⁵⁵ 

Random fluctuations in the long-time decay of $P R ( t )$ were found to be significant in the presence of strong noise. When the population decay becomes noisy (reaching the asymptotic value) before multiple fluctuations can be averaged over, the rate calculation is taken over the time period after the transient decay but before the onset of the breakdown of the sampled logarithmic population. In this intermediate region, exponential behavior is observed as seen in the typical population decays of Fig. 7. The decay rates can be extracted from the intermediate time interval that excludes the initial transients and the long-time numerically unconverged behavior.

In the regime of strong noise and low frequency driving, the calculated rate is highly dependent on the choice of the time window used to sample the decay. Although not shown, we confirmed that starting each trajectory with the same potential energy, as opposed to the same distance from the TS trajectory, did not significantly alter the measured rate. While this change in initial conditions drastically altered the shape of the transient section, the asymptotic region still exhibited approximately linear behavior. Beginning all trajectories with the same fixed initial position would also yield the same asymptotic decay rate, but would require the integration of significantly more trajectories to achieve similar numerical accuracy because the number of reactive trajectories decreases as the initial position is moved further from the barrier top region.

As postulated earlier, the difference in the unstable and stable Lyapunov exponents of the TS trajectory gives the rate of barrier crossing. We test this conjecture here on a system in an athermal environment under quasiperiodic driving as it represents an important limit of the most general conditions. In a thermal environment, violent fluctuations can occur due to the solvent-solute interaction and can thereby drown the effects of the driving conditions. By removing the solvent fluctuations while keeping the dissipation, the driving terms still dominate the dynamics allowing us to confirm that they are properly accounted for in the rate theory developed in Sec. IV. The rates obtained for a quasiperiodically driven barrier in a dissipative environment are shown in Fig. 8. The reaction rate increases with the anharmonic strength $𝜖$⁠. This characteristic behavior is independent of the driving form and environment, arising from the geometry of the energy surface as shown in the following examples. Across all irrational frequencies $Ω 1$ studied here, the numerical results are in agreement with the rates obtained from the Lyapunov exponents of the TS trajectory. This finding supports our assertion that the stability of the guiding phase space objects directly dictates the reaction rate.

Decay rates in a thermal system with a stationary barrier obtained directly from numerical integration of population decays are compared in Fig. 9 to those obtained from the Lyapunov exponents of the TS trajectory. There is good agreement across all parameter values. However, we observe that the difference in the Lyapunov exponents systematically overestimates the rates, if only slightly. With the inclusion of thermal driving, the rates calculated at each time point t_i are highly dependent on the realization of the noise. To verify that the average (34) includes a representative sample of noise sequences, the number of time points used for sampling was increased from 10 to about 50 with 10⁷ trajectories sampled at each t_i with spacing between time points $Δt= t i − t j =1$⁠. For all noise sequences, the same time window was used to extract the decay rate. This increased sampling resulted in only marginally better agreement with the rates calculated from stability analysis for $γ=0.1$ and $γ=1$⁠, implying that the numerical rates calculated using 10 time points have converged to a statistical accuracy relevant for at least qualitative comparison; in fact, the rates obtained from stability analysis are quantitatively close to the numerical data. In the strong-noise (⁠$γ=3$⁠) system, increasing the number of sampling points resulted in decreased variance in the numerical rates. This agrees with the previous observation that convergence of the rates requires sampling over many noise sequences.

In a thermal system where the motion of the barrier is periodic, the deterministic driving introduces nonequilibrium fluctuations that are not counteracted through dissipation, as is the case for the stochastic forces. However, the general trends for the reaction rates of this system remain the same (see Fig. 10). With increasing driving frequency or barrier anharmonicity, the reaction rate increases. Across all values of thermal strength studied, for high frequency driving the rates given by stability analysis are in excellent agreement with the numerically measured rates. This agreement at high frequencies has also been observed in a system driven only by periodic forces.^55,56 When the barrier is driven at a frequency close to the resonant frequency, the difference in the Lyapunov exponents underestimates the rate for $Ω=1$ and overestimates the rate for $Ω=3$⁠. This turnover corresponds to an increase in amplitude for the TS trajectory by a factor $∼2$⁠, as illustrated in Fig. 4. We conjecture that this amplitude effect results in more reactant population being pushed away from the barrier-top region over consecutive periods of oscillation. This gives rise to fluctuations in the decay, and correspondingly, the observed deviations from linearity in the logarithmic population. Thus, the interplay between friction $γ$⁠, thermal strength $σ$⁠, and driving frequency $Ω$ dictates the reaction rate through nonlinear relations. Based on the generality of our methodology, we expect the same trends to be observed in other reactions on energy surfaces with different anharmonic or more global geometries.

The PT rates computed for a noiseless and a noisy system subjected to periodic driving are compared in Fig. 11 to the difference between the Floquet exponents of the TS trajectory (an unstable periodic orbit) (top) and to the difference between the Lyapunov exponents of the TS trajectory (bottom), respectively. The rate constants extracted from the TS trajectory computed using PT of order n = 3 are in good agreement with the numerically exact result over a wide range of anharmonic couplings (shown in Fig. 11 for values upto $𝜖=10$⁠). In the case of a noisy system, the validity of PT is limited to much smaller anharmonicities. As should be expected, the larger the order in the PT expansion, the better the agreement of the rates to the numerical results. Although not shown, we have confirmed that this trend continues with the inclusion of higher-order terms.

In this paper, we have developed a new approach for computing the reaction rate of a thermally activated and externally driven chemical reaction across an anharmonic potential energy surface. It complements alternative approaches by other groups, cited above, which rely more directly on trajectories through the dividing surface, and extends our recent work^54–56 which relied on periodic driving. Thus, the current work is an advance in two major respects: First, an exact DS is defined that is never recrossed by any reactive trajectory. This recrossing-free DS allows one to identify all trajectories that cross the DS as reactive, as has been the goal of TST from its inception. It is interesting to note that this optimal DS is time-dependent because it is attached to the TS trajectory, which is the only trajectory that remains in the vicinity of the barrier top for all times. Second, we have also demonstrated that the reaction rate is given by the difference between the Lyapunov exponents of the TS trajectory. Third, we have shown that this methodology can be applied to systems driven by pulses with several frequencies, even if they are incommensurable. Finally, we are currently extending the methodology reported here to laser-driven systems in the high-dimensional phase space and also to systems that interact with the environment through colored noise. Extension and testing of our theory on chemical reactions in which the phase space density in the transition state region does not vanish due to the presence of metastable energy wells on the underlying energy surface is another possible direction for future work.

See supplementary material for the equations described in this work to obtain the TS trajectory which have been coded in Mathematica. A text file in the Mathematica format (.nb) and a printout in a portable document format (.pdf) have been made available.

This work has been partially supported by the the National Science Foundation (NSF) (USA) through Grant No. CHE-1700749, the Ministry of Economy and Competitiveness (MINECO, Spain) under Contract No. MTM2015-63914-P, and by ICMAT Severo Ochoa under Contract No. SEV-2015-0554. Travel between partners was partially supported by people mobility programs, and most recently by the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No. 734557.