When a chemical reaction is driven by an external field, the transition state that the system must pass through as it changes from reactant to product—for example, an energy barrier—becomes time-dependent. We show that for periodic forcing the rate of barrier crossing can be determined through stability analysis of the non-autonomous transition state. Specifically, strong agreement is observed between the difference in the Floquet exponents describing stability of the transition state trajectory, which defines a recrossing-free dividing surface [G. T. Craven, T. Bartsch, and R. Hernandez, “Persistence of transition state structure in chemical reactions driven by fields oscillating in time,” Phys. Rev. E **89**, 040801(R) (2014)], and the rates calculated by simulation of ensembles of trajectories. This result opens the possibility to extract rates directly from the intrinsic stability of the transition state, even when it is time-dependent, without requiring a numerically expensive simulation of the long-time dynamics of a large ensemble of trajectories.

Controlling the rate at which reactants transform to products, either to accelerate a chemical process or to bias a reaction toward a certain pathway, is fundamental to chemical physics. Such kinetic control can be achieved through forcing from an external field, leading to emergent behavior in molecular structure assembly,^{1–4} organic synthesis,^{5} ultracold chemical reactions,^{6} and single molecule spectroscopy.^{7} In these processes, reaction rates are typically obtained through transition state theory (TST).^{8–11} There are two major obstacles to the implementation of TST. First, reactive trajectories must be identified and, second, the flux of these reactive trajectories though a phase space dividing surface (DS) must be calculated. If this DS is recrossed by reactive trajectories, TST overestimates the rate. Only in cases where this DS is recrossing-free is TST formally exact.

In autonomous systems, the optimal DS is determined by a normally hyperbolic invariant manifold (NHIM).^{10,12–22} The study of NHIMs is the principle focus of modern reaction dynamics in so far as knowledge of their geometry inherently contains the determining characteristics of the reaction. However, even when a recrossing free DS can be found, a rate calculation can be intractable, especially for systems with many degrees of freedom, as large numbers of trajectories must be integrated to yield statistically relevant results.

When a reaction is subjected to a time varying external force, the geometric structures of TST are known to exist in several cases, though they become time dependent.^{23–29} For chemical reactions that are induced solely by an external field, the coupling of the field with the reacting molecule's dipole moment can accelerate the reaction rate,^{30} even for systems that dissipate energy through a spontaneous emission process.^{31}

An example of a molecular process where an external force influences the transition state (TS) geometry, and thus reaction rates, is the photoinduced isomerization between *cis* and *trans* stilbene (Ph-C=C-Ph).^{32–34} Its unimolecular reaction path can be parameterized through the torsion angle of the C=C double bond. Changing the energetics along this path through photoinduction alters the isomerization reaction rate.

We show here that when a chemical reaction is periodically forced by an external field (such as a laser), the reaction rates are determined directly by the stability of the transition state. We calculate the reaction rate of a model system by simulating large ensembles of trajectories and compare this result with the rate predicted by Floquet analysis of the transition state trajectory. Corresponding to the “chemical method” where the reactant concentration is followed as a function of time,^{35} we obtain reaction rates from the decay of a given initial distribution. These rates are well-defined because the decay is exponential when averaged over a period of the driving and independent of the choice of distribution. A major result of this work is that the rates can be obtained from a Floquet analysis of the transition state trajectory, an unstable periodic orbit (PO) close to the barrier top. This agreement suggests that chemical reaction rates can be extracted directly from the transition state without knowledge of the dynamics of the reactive population. This general result could have been anticipated from the known connection between the stability of periodic orbits of Hamiltonian systems and rates,^{36–38} but is here established even in the case of driven systems.

To model barrier crossings in chemical reactions driven by a time-dependent external field *E*(*t*) we consider a particle of unit mass with an initial position *x*_{0} on the reactant side of a moving energy barrier. The chosen barrier is a quartic potential of the form

which leads to the equations of motion

where γ is a dissipative emission parameter, ω_{b} is the barrier frequency, and ε is an anharmonic coefficient. The anharmonic coefficient is restricted to values ε ⩾ 0 such that there is a single maximum in the potential located at the barrier top (BT). The time dependent, instantaneous position of the BT is specified by *E*(*t*). Figure 1 shows the time evolution of *x*(*t*) for an ensemble of trajectories following Eq. (2). Each trajectory either crosses the energy barrier forming product or remains on the reactant side, never surmounting the barrier. The normalized flux of *reactive* trajectories through the phase-space bottleneck—the TS—is the reaction rate.^{8}

Every realization of the forcing *E*(*t*) has a special trajectory imbedded in the dynamics (2) that remains close to the BT for all time. This bounded trajectory, termed the TS trajectory,^{39–43} will never descend into the product or reactant regions.^{29} As illustrated in Fig. 1, the TS trajectory does not follow the time evolution of the energetic maximum given by the BT. It is instead a specific trajectory that responds to motion of the BT in such a way that it remains bounded for all time. When *E*(*t*) is a periodic function with period *T* such that *E*(*t*) = *E*(*t* + *T*) for all *t*, the resulting TS trajectory is a PO with the same period *T*.

Attached to the TS trajectory are stable and unstable manifolds. The stable manifold intersects a line of initial conditions *x* = *x*_{0} at a critical velocity *V*‡.^{42,43} A trajectory will surmount the energy barrier, moving from the reactant state to the product state, if *v*_{0} > *V*‡. If *v*_{0} < *V*‡, the trajectory is nonreactive. The extension of this point to all values of *x*_{0} creates a critical curve

To calculate rates, the TST methodology is concerned with creating a DS that is crossed once and only once by reactive trajectories and then evaluating the flux through that surface. For the case when

^{29}the configuration space projection of the TS trajectory is free of recrossings.

For the case of a harmonic barrier (ε = 0), Eq. (2) can be solved analytically with eigenvalues

in terms of the *S* functionals^{40,44}

that guarantee the appropriate boundary conditions for *t* → ±∞. The TS solution for any barrier motion is given by Eq. (3).

For anharmonic barriers (ε ≠ 0), the TS trajectory will be an unstable PO close to the barrier top, as in the harmonic case. Its period will typically coincide with the period *T* of the external driving. The anharmonic equations of motion (2) are not amenable to an exact analytical solution, although approximate analytical methods have previously been employed.^{42,43} Instead, we obtain the TS trajectory

The barrier crossing rates for Eq. (1) were calculated by simulating ensembles of trajectories driven by an external field of the form *E*(*t*) = *a*sin (Ω*t* + φ). For single mode sinusoidal driving, the TS trajectory is a PO with period 2π/Ω. Physical units were set by normalizing *a* and ω_{b} to unity, making all other parameters dimensionless. Each trajectory was given an initial position *x*_{0} = −0.1 to the left of the instantaneous barrier top and *v*_{0} was sampled from two separate distributions: (1) a Boltzmann distribution *q*_{B} with *k*_{B}*T* = 1, and (2) a uniform distribution *q*_{U} (bounded over the region [*V*‡ − 1/2, *V*‡ + 1/2]). For each parameter set {Ω, γ, ε}, 10^{8}-10^{9} trajectories were simulated. The normalized reactant population *P*_{R}(*t*) is obtained from a histogram of those trajectories that are on the reactant side of the TS trajectory at time *t*. Assuming first order kinetics, the scaled logarithm of the normalized population, −ln [*P*_{R}(*t*) − *P*_{R}(∞)], should be linear in time. As illustrated in Fig. 3, after transient trajectories have crossed, the decay of the logarithmic population is linear up to periodic modulation, and the first order assumption is confirmed. Periodic fluctuations are noticeable for small driving frequencies (Ω ⪅ 2) and large anharmonicities due to effect of higher order terms in the asymptotic decay of *P*_{R}(*t*). The slope of a least squares fit to the non-transient section of the data gives the reaction rates calculated from simulation *k*_{f}.

We now focus on analysis of the TS trajectory and the determination of reaction rates from its intrinsic stability. With the bounded TS trajectory now defined, a relative coordinate system can be introduced. In relative coordinates,

the equations of motion read

The last term represents a time-dependent driving for the relative dynamics that does not depend on the current trajectory. It ensures that the relative equations of motion have a fixed point

*x*= Δ

*v*= 0, i.e., on the TS trajectory.

The long-time decay rate of *P*_{R}(*t*) is determined by the behavior of trajectories close to the stable manifold. Once a trajectory is sufficiently close to the TS trajectory, it can be described by a linearization of the equations of motion (6),

where *a*(*t*) = *U*(*x*^{‡}(*t*)). In the phase space vector coordinate

where

is the Jacobian of Eq. (6) about

where the fundamental matrix solution ** σ**(

*t*) is a 2 × 2 matrix that satisfies

where ** I** is the identity matrix.

The fundamental matrix for one period of

*m*

_{u, s}are called Floquet multipliers. The Floquet exponents μ

_{u, s}= 1/

*T*ln |

*m*

_{u, s}| give the rates by which nearby trajectories approach or recede from

^{45}For a harmonic barrier, the multipliers are bounded according to

_{u}and a negative exponent μ

_{s}. We will assume that this qualitative condition is also satisfied for the anharmonic barriers; we neglect the possibility that for strong anharmonicities bifurcations of the TS trajectory might occur.

Let

**. By Floquet's theorem and the positivity of the Floquet multipliers, the vectors**

*M*are periodic in time with period *T*. In the coordinate system defined by these vectors,

the linearized equations of motion (8) read

with the solution

Therefore, the vectors

where α_{u, s} are the first components of the vectors

*T*. A trajectory with given initial conditions

*C*

_{u}and

*C*

_{s}will cross the moving dividing surface Δ

*x*= 0 at time

*t*determined by

If the initial condition *C*_{s} is fixed and a trajectory with a certain value of *C*_{u} crosses the moving DS at time *t*, Eq. (17) shows that a trajectory with initial value

*t*+

*T*. Iteration then leads to the existence of trajectories with initial values

*t*+

*nT*.

Now consider an arbitrary ensemble of initial conditions with fixed *x*(0) on the reactant side and with a fixed value *C*_{s} < 0 small enough to be in the region of phase space where the linear approximation (8) is accurate. In this region, the phase space density is constant up to linear corrections in the distance from the stable manifold and the number of trajectories that cross the DS in a given time interval is proportional to the width of the strip that contains these trajectories. From one period to the next this width decreases by a factor

*P*

_{R}(

*t*) is proportional to

*k*

_{f}= μ

_{u}− μ

_{s}, which states that the rate of barrier crossing is the difference in the Floquet exponents. Note that we have made no assumption for the energy distribution and thus this rate is independent of the ensemble of initial conditions.

A comparison between the rates calculated from numerical simulation *k*_{f}, for both the Boltzmann *q*_{B} and uniform *q*_{U} distributions, and rates predicted by the Floquet exponents μ_{u} − μ_{s} is shown in Fig. 4. For all values of the forcing frequency Ω, dissipative parameter γ, and anharmonic strength ε, the numerical rate is in agreement with rate predicted by stability analysis. This result opens the possibility that when chemical reactions are forced by periodic external fields the reaction rates can be extracted from knowledge of the stability of the TS trajectory. The extension of TS trajectory stability analysis to aperiodically forced or thermally activated reactions is a focus of our future research.

This work has been partially supported by the National Science Foundation (NSF) through Grant No. NSF-CHE-1112067. Travel between partners was partially supported through the People Programme (Marie Curie Actions) of the European Union's (EU) Seventh Framework Programme FP7/2007-2013/ under REA Grant Agreement No. 294974.