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 x0 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
The time evolution of x(t) for a swarm of trajectories following Eq. (2) with E(t) = asin (Ωt + φ) are shown in black (below). The potential surface of Eq. (1) is shown above with a contour plot shown below. The BT and TS trajectories are shown in dashed and solid white, respectively. Time is shown in units of τ = Ωt/2π + 3/4. The initial velocities are sampled from qB. Parameters are ε = 1, Ω = 3, γ = 4, and φ = 0.
The time evolution of x(t) for a swarm of trajectories following Eq. (2) with E(t) = asin (Ωt + φ) are shown in black (below). The potential surface of Eq. (1) is shown above with a contour plot shown below. The BT and TS trajectories are shown in dashed and solid white, respectively. Time is shown in units of τ = Ωt/2π + 3/4. The initial velocities are sampled from qB. Parameters are ε = 1, Ω = 3, γ = 4, and φ = 0.
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 = x0 at a critical velocity V‡.42,43 A trajectory will surmount the energy barrier, moving from the reactant state to the product state, if v0 > V‡. If v0 < V‡, the trajectory is nonreactive. The extension of this point to all values of x0 creates a critical curve
![FIG. 2. Phase space plots for a swarm of trajectories following Eq. (2) with E(t) = asin (Ωt + φ). The initial position for every trajectory, x0, is shown as a dashed black line. Reactive trajectories are colored in blue and nonreactive trajectories are colored in orange with respective basins separated by \documentclass[12pt]{minimal}\begin{document}$V^\ddag _c$\end{document}Vc‡ (solid red). The TS trajectory \documentclass[12pt]{minimal}\begin{document}$\bm {\Gamma }^\ddag$\end{document}Γ‡ is shown in black. The critical velocity V‡ is indicated by a red circle at the intersection of the dashed red line and x0. The initial velocities are sampled from qB. Parameters are Ω = 5, γ = 2, and φ = 0.](/view-large/figure/62084991/041106_1_f2.jpg)
Phase space plots for a swarm of trajectories following Eq. (2) with E(t) = asin (Ωt + φ). The initial position for every trajectory, x0, is shown as a dashed black line. Reactive trajectories are colored in blue and nonreactive trajectories are colored in orange with respective basins separated by
Phase space plots for a swarm of trajectories following Eq. (2) with E(t) = asin (Ωt + φ). The initial position for every trajectory, x0, is shown as a dashed black line. Reactive trajectories are colored in blue and nonreactive trajectories are colored in orange with respective basins separated by
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
For the case of a harmonic barrier (ε = 0), Eq. (2) can be solved analytically with eigenvalues
in terms of the S functionals40,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) = asin (Ω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 x0 = −0.1 to the left of the instantaneous barrier top and v0 was sampled from two separate distributions: (1) a Boltzmann distribution qB with kBT = 1, and (2) a uniform distribution qU (bounded over the region [V‡ − 1/2, V‡ + 1/2]). For each parameter set {Ω, γ, ε}, 108-109 trajectories were simulated. The normalized reactant population PR(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 [PR(t) − PR(∞)], 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 PR(t). The slope of a least squares fit to the non-transient section of the data gives the reaction rates calculated from simulation kf.
Time dependence of the scaled logarithm of the reactant population for Ω = 2 and Ω = 7 with v0 sampled from qB. The slope of each dashed line is the barrier crossing rate kf, corresponding to a respective ε value. Parameters are γ = 1 and φ = 0.
Time dependence of the scaled logarithm of the reactant population for Ω = 2 and Ω = 7 with v0 sampled from qB. The slope of each dashed line is the barrier crossing rate kf, corresponding to a respective ε value. Parameters are γ = 1 and φ = 0.
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
The long-time decay rate of PR(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
Let
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
If the initial condition Cs is fixed and a trajectory with a certain value of Cu crosses the moving DS at time t, Eq. (17) shows that a trajectory with initial value
Now consider an arbitrary ensemble of initial conditions with fixed x(0) on the reactant side and with a fixed value Cs < 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
A comparison between the rates calculated from numerical simulation kf, for both the Boltzmann qB and uniform qU 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.
The barrier crossing rates for Hamiltonian (top) and dissipative (bottom) systems following Eq. (2). The numerically calculated rates kf for the distributions qB (circles) and qU (squares) are shown with units given on the left axes. The solid curves are the rates predicted by the difference in the Floquet exponents μu − μs of the TS trajectory with units given on the right axes.
The barrier crossing rates for Hamiltonian (top) and dissipative (bottom) systems following Eq. (2). The numerically calculated rates kf for the distributions qB (circles) and qU (squares) are shown with units given on the left axes. The solid curves are the rates predicted by the difference in the Floquet exponents μu − μs of the TS trajectory with units given on the right axes.
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.
