[AIP

The Journal of Chemical Physics

Chemical reactions induced by oscillating external fields in weak thermal environments

Galen T. Craven, Thomas Bartsch, and Rigoberto Hernandez

Citation: The Journal of Chemical Physics 142, 074108 (2015); doi: 10.1063/1.4907590 View online: http://dx.doi.org/10.1063/14907590

View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/7?ver=pdfcov Published by the AIP Publishing

Articles you may be interested in

Reaction efficiency effects on binary chemical reactions J. Chem. Phys. 141, 104103 (2014); 10.1063/1.4894791

Validity conditions for moment closure approximations in stochastic chemical kinetics J. Chem. Phys. 141, 084103 (2014); 10.1063/1.4892838

Communication: Transition state trajectory stability determines barrier crossing rates in chemical reactions

induced by time-dependent oscillating fields

J. Chem. Phys. 141, 041106 (2014); 10.1063/1.4891471

Effect of an external electric field on the diffusion-influenced geminate reversible reaction of a neutral particle and a charged particle in three dimensions. III. Ground-state ABCD reaction J. Chem. Phys. 139, 194107 (2013); 10.1063/1.4830401

Qualitative dynamics of generalized Langevin equations and the theory of chemical reaction rates J. Chem. Phys. 116, 2516 (2002); 10.1063/1.1436116

Launching in 2016!

The future of applied photonics research is here

Photonics

THE JOURNAL OF CHEMICAL PHYSICS 142, 074108 (2015)

CrossMark

£elick for updates

Chemical reactions induced by oscillating external fields in weak thermal environments

Galen T. Craven,1 Thomas Bartsch,2 and Rigoberto Hernandez1,3'

1 Center for Computational Molecular Science and Technology, School of Chemistry and Biochemistry, Georgia Institute of Technology, Atlanta, Georgia 30332-0400, USA

2Department of Mathematical Sciences, Loughborough University, Loughborough LE11 3TU, United Kingdom (Received 8 November 2014; accepted 24 January 2015; published online 18 February 2015)

Chemical reaction rates must increasingly be determined in systems that evolve under the control of external stimuli. In these systems, when a reactant population is induced to cross an energy barrier through forcing from a temporally varying external field, the transition state that the reaction must pass through during the transformation from reactant to product is no longer a fixed geometric structure, but is instead time-dependent. For a periodically forced model reaction, we develop a recrossing-free dividing surface that is attached to a transition state trajectory [T. Bartsch, R. Hernandez, and T. Uzer, Phys. Rev. Lett. 95, 058301 (2005)]. We have previously shown that for single-mode sinusoidal driving, the stability of the time-varying transition state directly determines the reaction rate [G. T. Craven, T. Bartsch, and R. Hernandez, J. Chem. Phys. 141, 041106 (2014)]. Here, we extend our previous work to the case of multi-mode driving waveforms. Excellent agreement is observed between the rates predicted by stability analysis and rates obtained through numerical calculation of the reactive flux. We also show that the optimal dividing surface and the resulting reaction rate for a reactive system driven by weak thermal noise can be approximated well using the transition state geometry of the underlying deterministic system. This agreement persists as long as the thermal driving strength is less than the order of that of the periodic driving. The power of this result is its simplicity. The surprising accuracy of the time-dependent noise-free geometry for obtaining transition state theory rates in chemical reactions driven by periodic fields reveals the dynamics without requiring the cost of brute-force calculations. © 2015 AIP Publishing LLC. [http://dx.doi.org/10.106371.4907590]

I. INTRODUCTION

Optimal control of reaction pathways in systems undergoing configurational changes can be achieved through forcing from tailored external fields. These fields can be tuned to induce specific deformations on a potential energy surface, providing control of state-to-state transitions.1-3 In these processes, a formally exact classical rate calculation can be obtained through modern-day transition state theory (TST).4-9 The principal assumption of microcanonical TST is that there exists a hypersurface between reactant and product confirmations that is crossed only once by reactive trajectories during the traversal of a free energy barrier separating these basins. The TST reaction rate is calculated from the flux through this dividing surface (DS). If the DS is recrossed by reactive trajectories, TST will give an overestimate to the classical reaction rate. Determination of a recrossing-free DS therefore leads to classically exact microcanonical TST rates. If there exists a Boltzmann distribution of reactant states during the reaction, then the Laplace transform of these rates leads to a canonical TST rate that is also exact.

A phase space DS that is free of recrossings can be constructed in conservative systems at energies close to the reaction threshold. In systems with two degrees of freedom,

a)Author to whom correspondence should be addressed. Electronic mail: hernandez@chemistry.gatech.edu.

the optimal DS is the configuration space projection of an unstable periodic orbit (PO).10-13 In systems with higher dimensionality, the generalization of this PO is a normally hyperbolic invariant manifold (NHIM).14-22 The NHIM bounds the TS, being one less in dimension.16 It defines a recrossing-free surface at energies below bifurcation thresholds.23-25 Reactive trajectories are mediated by stable and unstable manifolds (reaction pathways) attached to the NHIM. These pathways persist even in reactions whose state-to-state transitions are not dictated by purely configurational changes.26

In systems subjected to time-varying external forcing, the characterization of the NHIM as a hypersphere of constant energy breaks down. For example, field-matter interactions constitute processes in which energy is exchanged with a reacting system. These interactions lead to emergent and controllable behavior in assembly phenomena,27-30 organic synthesis,31 protein folding,32 the detection of DNA,33 and photodissociation.34 Knowledge of the mechanism by which these interactions mediate reactive flow provides a methodological tool in the design of molecular devices with unique functionality.35-37

Materials that undergo conformational changes in response to an external trigger offer examples of such emergent technology.36,38-41 Stimuli such as thermal variations, electric fields, and photoinduction have been used as triggers for the conversion of chemical energy into mechanical work.37,42-44 Assemblies that convert chemical energy into directional

0021-9606/2015/142(7)/074108/13/$30.00

142, 074108-1

©2015 AIP Publishing LLC

motion can be achieved through isomerization reactions which are induced either from light or applied electric fields.35'45'46 In these responsive materials, controlling the rate and pathway at which reactants transform to products is fundamental to harness mechanical actions for applicative purposes.

The aim of this paper is to develop a rate theory for reactions that are driven by periodic external fields in weak thermal environments. In the absence of noise, a dissipative system that is periodically driven admits a DS that is free of recrossings.47 This structure differs from the canonical view of the TS wherein the TS is a structure fixed in time at a saddle point on the potential energy surface. Here, we develop a rate theory based on reactive flux through this recrossing-free DS and the stability of the corresponding TS. In Ref. 48, we found that the stability of the moving TS directly determines the reaction rate for single-mode sinusoidal driving. In conservative systems, stability analysis is known to characterize molecular motions as well as determine the rate of configurational transitions.24,49-52 Building on our previous work, we test the viability of stability analysis to determine reaction rates in systems driven by multi-mode waveforms with no thermal driving. The extent to which the accuracy of the rate theory relying on the noise-free geometry persists in systems that are coupled to a thermal bath is also verified through inclusion and variation of the thermal driving strength.

The outline of this paper is as follows: In Sec. II, a dynamical system is introduced to model barrier crossings in chemical reactions forced by periodic external fields. In Sec. III, a dividing surface that is recrossing-free is constructed for this model in the absence of thermal driving. Section IV contains analytical theories to predict the reaction rates of driven reactions by calculation of the reactive flux through this dividing surface for both globally non-linear and locally linear dynamics. Comparison to the computational rates, computed from numerical integration of large ensembles of trajectories, is presented in Sec. V. Although not considered earlier, the effect of noise on the rate of these driven systems has also been addressed in Sec. VI. We find that the rates computed from the noise-free geometry are accurate up to relatively large values of the friction and sometimes even in the thermal regime.

II. MODEL DETAILS

The interaction of an external field with a reactant species can strongly influence the mechanism and rate of a reaction.3'53,54 As a paradigmatic example of a chemical reaction driven under kinetic control, we consider a particle of unit mass moving along a reaction coordinate x. The trajectory of the particle begins at a position x0 on the reactant side of an energy barrier that is moving in space under the influence of a time-dependent external field E(t). The chosen potential surface is the quartic form

U(x) = - 2- E(t))2 - 4 e(x - E(t))4.

With the inclusion of additional non-conservative dissipation as well as stochastic driving forces, a particle at a phase space point r = (x, v ) moving according to the potential (1) can be described by the Langevin equations of motion

x = v,

V = -yv + w¡2(x - E(t)) + e(x - E(t))3 + V2^£a(t),

where y > 0 is a dissipation parameter, wb is the barrier frequency, and e is an anharmonic coefficient. Thus, for e + 0 the coordinate of the particle is non-linearly coupled to the moving barrier. By restricting the anharmonic coefficient to values e > 0, there is a single maximum in the potential located at the BT. The random fluctuating force Ça(t) is Gaussian white noise obeying the statistical properties

<&(t)> = 0 and <&(t)&(t')> = S(t - t'),

where a denotes a specific noise sequence. The strength of the noise is varied through the parameter < > 0.

Depending on the geometry of U(x), initial conditions, as well as the specific realization of the thermal environment and the external field, a trajectory will either surmount the energy barrier and form product or remain on the reactant side. By calculation of the normalized flux of reactive trajectories through the TS, the classical reaction rate for a system evolving through (2) can be obtained.6

In this article, we consider periodic external driving of the form

E(t) = a [ sin(w t + p),

The time dependent, instantaneous position of the moving barrier top (BT) is specified by E(t).

where c R is a finite set of frequencies. The waveforms consist of a fundamental frequency Q, and convolutions of this fundamental with higher order partial frequencies. Three frequency sets are considered: the single fundamental frequency ^ = {Q}, the fundamental and the second partial frequencies = {Q, 2Q}, and the fundamental, second, and third partial frequencies = {Q, 2Q, 3Q}. The fundamental driving frequency Qf is Q for ^ and and 2Q for The products in Eq. (4) for the three sets can be written as finite sums

E1(t) = a sin(Qt + p),

E2(t) = 2 (cos(Qt) + cos(3Q t + 2p)), (5)

E3(t) = a (sin(2Qt + p) + sin(4Qt + p) + sin(6Qt + 3p) + sin(p)),

where the leading order terms exhibit the characteristic fundamental frequency. The maximum amplitude of each waveform is set to unity by adjusting the value of the parameter a accordingly. For the and sets, a = 1, a ~ 1.299,

and a ~ 1.822, respectively. The functional forms of (5) are shown in Fig. 1.

III. THE TRANSITION STATE TRAJECTORY

The construction and existence of a structure whose configuration space projection is free of recrossings is

FIG. 1. Functional forms of the periodic driving E (t) for the frequency sets (left), (middle), and (right) with fundamental driving frequency Qf = 4 are shown in the top row. The corresponding TS trajectory for each frequency set and anharmonic parameter values e e {0,4,8} are shown in the bottom row. Various extrema of E (t) are denoted by circles with arguments shown as dashed vertical lines. The corresponding extrema of x*(t) are denoted by circles and colored according to the respective e value. Parameters for all panels are r = 0, y = 1, and f = 0 in dimensionless units.

dependent on the mechanism and geometry of a given reaction. In an ion-pair dissociation, for example, Mullen et al.55,56 claimed that a recrossing-free DS does not exist, whereas Truhlar and Garrett57 had earlier found that recrossings could be eliminated through variation of the DS into an optimized orientation. Although, we have not specifically addressed that mechanism, we have previously shown that for a periodically driven system with no thermal driving, an optimal (recrossing-free) DS can be readily obtained. It is associated with an unstable PO in the region of the BT.47 Moreover, a DS that is free of recrossings is known to exist in thermally driven systems for the case of a harmonic barrier.58,59 The time evolution of the configuration space projection of this DS has been termed the transition state trajectory.47,48,58-62 It has not yet been proven that this is the only object which is free of recrossings over an arbitrary finite time interval. Nevertheless, all that is needed here is its existence, and the configuration projection of the TS trajectory defines a DS that is recrossing-free.

The TS trajectory is a specific trajectory that never descends into either the product or reactant regions, remaining bounded to BT for all time. For the system (2), it is a moving saddle point to which stable and unstable manifolds can be attached. All trajectories that exponentially approach the TS trajectory as t ^ ^ are contained on the stable manifold. These trajectories will never descend from the BT region and therefore separate reactive from nonreactive trajectories in phase space. The unstable manifold is formed from trajectories that approach the TS trajectory as t ^ The role of the unstable manifold is less important for the purposes considered here.

For an arbitrary driving E(t) of a harmonic (e = 0) potential, the equations of motion can be solved exactly and an exact form of the TS trajectory can be obtained. The eigenvalues of (2),

= -1 (y ±ylY2 + 4wb) '

correspond to the stable and unstable manifolds. The S

functionals54,59

St [ß,g; t] =

g(r) exp(ß(t - r)) dr : Re ß > 0,

+ / g(r) exp(ß(t - r)) dr : Re ß < 0,

obtained as a Green's function solution, suppress the transient exponential factor in the solution and return only the equilibrium portion. In the absence of thermal driving (<r = 0), the TS trajectory for a harmonic barrier can therefore be expressed

as47,48,61,62

x*(0 = t—V E; t] - S[Xu,E; t]),

v*(t) =

Xu - Xs

(XsS[Xs, E; t] - XuS[Xu, E; t]).

For the case of T -periodic motion of a harmonic barrier, the TS trajectory can be identified more easily by looking for a bounded solution to the equations of motion. For the single frequency case, the Ansatz

x\(t) = A sin(Qt + p) + B cos(Qi + p) (9)

yields the solution

A = Ai, B = Bi,

Ak = a

Bk = a

wb2(wb2 + (k ft)2) (y k Q)2 + (wg + (k Q)2)2,

w2hy k Q (y k Q)2 + (wb + (k Q)2)2.

In the absence of friction, y = 0, this simplifies to

A1 = a-

Bi = 0.

1 + ft2/wb'

In this case, the TS trajectory will oscillate in phase with the barrier, but with smaller amplitude A1 < a.

For the and cases, Ansätze can be constructed through Fourier series expansion of Eq. (4) yielding the solutions

x2(t) = 1 (A1 cos(Qt) - B1 sin(Qt)

- A3 cos(3Qt + 2p) + B3 sin(3Qt + 2p)) (13)

x3(t ) = 4 ( A2 sin(2Qt + p) + B2 cos(2Qt + p) + A4 sin(4Qt + p) + B4 cos(4Qt + p) - A6 sin(6Qt + 3p) - B6 cos(6Qt + 3p)

+ a sin(p)),

respectively.

For periodically driven anharmonic barriers (e + 0), the TS trajectory remains an unstable PO47 but it does not admit to an exact solution of the system of Eq. (2). Nevertheless, the phase space vector of the TS trajectory r* = (x*(t), v *(t)) is a bounded solution to the equations of motion. To find this bounded solution to arbitrary accuracy, numerical Newton-Raphson root finding methods were applied.

The dynamics of x*(t), shown in Fig. 1, illustrate the result that the instantaneous position of the TS trajectory does not correspond to the energetic maximum of the potential surface. For dissipative systems (y + 0), x*(t) will either lag behind in phase, as is the case for both the and sets, or advance in phase as is the case for the Q2 set, with respect to motion defined by E(t). Also note that x*(t) oscillates with a smaller amplitude than E(t). Thus, even for in-phase oscillations, e.g., when y = 0, it will not correspond to the location of an energetic saddle point. Figure 1 also shows the dependence of x*(t) on the anharmonic parameter e. As e is increased, the curvature of the energy barrier increases. Non-intuitively, this results in a larger amplitude of oscillation for x*(t) to remain bounded to the BT. This trend persists for all

For dynamical analysis, it is advantageous to introduce a coordinate system which has a fixed point at the origin. In relative coordinates

Ax = x - x*(t), Av = v - v*(t), the equations of motion read Ax = Av,

Av = -yAv - U'(Ax + x*(t)) + U'(x*(t)).

The relative equations of motion have a fixed point Ar* at Ax = Av = 0, i.e., on the TS trajectory, and the surrounding vector field itself will now oscillate with period T, the same period as the driving. The TS trajectory has both a stable and an unstable manifold attached. In relative coordinates, the directions of these manifolds will depend on time.

IV. REACTION RATE THEORY

In the TST formalism, the rate of a chemical reaction is

given by the time-dependence of the conversion process from reactant to product (R ^ P) where a DS in either configuration

space or phase space separates the reactive constituents. The reaction rate can be obtained from the dynamics of the normalized reactive population (PR ^ PP) either through analytical propagation of the phase space density of initial conditions or by treating large numbers of trajectories as discrete sets, and integrating the equations of motion.

Consider a set of trajectories evolving through (2) that all have initial positions x0 < x*(0) on the reactant side of the moving surface. The initial position distribution at time t = 0 is ¿(x - x0) and the initial phase space density is

po(x,v) = ¿(x - xo) q(v ),

where q(v) is a Boltzmann distribution. The initial velocity v0 of each trajectory is sampled from q(v), although at later times due to dissipation and driving this distribution will

not be conserved. A fraction of this initial density contains reactive trajectories. From the survival probability of PR, the reaction rate can be expressed as the instantaneous flux-overpopulation. The flux calculation is formally exact because the DS attached to the TS trajectory is recrossing-free.

A. Harmonic barriers

When the barrier is harmonic (e = 0), reactive trajectories will cross the moving DS at a time60

_ln I Avp - XuAxp

Xu - Xs \ Avo - ^sAxo

The crossing time is a monotonically decreasing function of the initial velocity Av0: fast trajectories cross earlier. It diverges as Av0 ^ XsAx0 approaches the stable manifold, and it tends to zero as Av0 ^ to.

At any time t > 0, the product region Ax > 0, to the right of the moving surface, will contain all those trajectories that cross the surface at a time t* < t. These are the trajectories that have an initial velocity of at least vmin = v *(0) + Avmin, where t*(Avmin) = t. From this condition, we obtain

Xue-Xut - . Avmin = -r--r— Axo.

The population of the product region at time t is therefore

/ q(v) d v, (20)

PP(t) = / q(v) d v,

vmin(t)

(16) and the flux across the moving surface is

FM(t ) = dP^

= -q(vmin(t )) = -q(vmin(t ))

dAVmin dt

= -q(Vmin(t)) Axo (Xu - Xs)

(Xu+Xs)t

(eXuf - eXsf)2

. (21)

This result is positive because Ax0 < 0.

Alternatively, the flux can be calculated directly from the flux integral

FM(t) = dAv Av pt(Ax = o, Av ),

where pt(Ax, Av) is the density of trajectories in phase space at time t. Initially, this density is

P0(Ax, Av) = ¿(Ax - Ax0) q(v*(0) + Av). (23)

At later times, it can be obtained from

pt(Ax, Av ) = eYt po(Ax(-t), Av (-t )).

Here, Ax(-t) and Av(-t) denote the phase space point reached from Ax, Av by propagating backwards to -t, i.e., it is the initial condition that has reached Ax, Av at time t. The exponential prefactor accounts for the shrinkage of phase space volume: The relative dynamics stretches distances at a rate in the u direction and by a rate < 0 in the s

X(t) OL-,{t) X(t) X{t)

FIG. 2. Phase space plots for a swarm of trajectories following the equations of motion (2) with various driving frequency sets and parameter values: with Qf = 5 and e = 1 for (a) r = 0 (noiseless) and (b) r = 0.025 (thermal) driving, with Qf = 4 and e = 3 for (c) r = 0 and (d) r = 0.025. The initial position for every trajectory, x0 = -0.1, is shown as a vertical solid white line. Reactive trajectories are colored in cyan and nonreactive trajectories are colored in orange. The TS trajectory r * is shown as a solid white curve. The critical velocity V* is indicated by a circle at the intersection of the dashed horizontal line and the line of initial conditions. The solid red curve is the critical curve Vj and the solid blue curve is the harmonic critical curve Vj |e=0. Parameters for all panels are y = 1 and i = 0 in dimensionless units.

direction. Volumes therefore are "stretched" at a constant rate Xu + Xs = -y < 0, and densities must increase accordingly.

Because the relative dynamics is linear, the equations of motion can be solved explicitly. The result is

Ax(-t) = ax Ax + av Au, Au (-t) = bx Ax + bv Au,

Xu e-Xst - Xs e-Xut

e-_e-Xst

bx = -

Xu - Xs

XuXs(e

-Xuf _,

Xu - Xs

Xu - Xs ,

Xu e-Xut - Xs e-Xst

Xu - Xs

We thus obtain the flux integral

FM(t) = eyt f dAu Au ô(av Au - x0) q(u*(0) + bv Au) Jo

= eyt H dAu Au ^(Au - xo/av) q(u*(0) + bv Au) Jo |av I

= _e_ Ol q (v* (0) + bv/av xq) , -a, a,

which can be shown to agree with Eq. (21).

In the limit t ^ the minimum velocity (19) is approximately

Aumin = Xs Axe - (Xu - Xs)Axc e"(x""xs)t + O (e"2(x""xs)t).

As expected, it tends to Xs Ax0, which is the location of the stable manifold. Therefore,

VminM = v * (0) + Xs A Xq = V *.

The critical velocity V* is determined by the point of intersection between the stable manifold and the line x = x0 of initial conditions.61,62 The identification of V* allows the separation of reactive (u0 > V*) and nonreactive trajectories (u0 < V* ) from initial conditions. The stable manifold at t = 0 can be calculated through extension of this point to all values of x0 and defines a critical curve Vc* . As illustrated in Figs. 2(a) and 2(c), V* is a time-invariant phase space object which separates the reactive and nonreactive basins.

FIG. 3. The asymptotic product population Pp(m) of the harmonic potential (e = 0) as a function of driving frequency Q for the frequency set. The curves are colored with respect to the value of the initial phase shift: i = 0 (blue), f = n/2 (cyan), and f = n (red). For each value of f, the dependency of the asymptotic population on the friction parameter y is shown by varying the linestyle: y = 0 (solid), y = 1 (dashed), and y = 2 (dotted).

The product population in the long-time limit is

Pp(to) = q(u ) d u.

J v*(0)+^s Ax0

As shown in Fig. 3, the asymptotic population of the product region depends strongly on the frequency of the barrier motion Q, the initial phase p, and the friction y. The asymptotic value is approached according to

Pp(t) = Pp(to) -

vmin(t )

q(v ) d v

= PP(«0 + qtUminM) (Xu - Xs) Axe e-(Xu-Xs)t

+ o (e-2(X»-Xs)t) .

The rate of approach, i.e., the barrier crossing rate is

Xu - Xs = ^y2 + 4^2.

It depends only on the damping and the shape of the barrier, but not on the details of the barrier motion or the distribution of initial conditions (unless q(umin(TO)) happens to vanish).

B. Anharmonic barriers

Analogous to the harmonic case, for an anharmonic barrier, we assume that there is an unstable PO with the period of the driving. This is similar to the POs used by Lehmann et al.63-65 for the case of thermal activation with additive periodic driving. Note that this TS trajectory r* is an exact solution to the equations of motion. As r* is an unstable PO, it has stable and unstable manifolds attached. The manifolds are uniquely defined and can be calculated perturbatively61,62 or with a numerical scheme. The dependence of these manifolds on e and the corresponding phase space reaction dynamics are shown in Figs. 2(a) and 2(c). With increasing anhar-monicity, V* also increases due to curvature in the stable manifold. This results in a decrease in fraction of trajectories that surmount the barrier leading to products.

The nonlinear equations of motion (16) cannot, in general, be solved exactly. Let

OCT,, to; t ) =

y x(ro, to; t) y„(ro, to; t )

represent the phase space point that is reached at time t by a trajectory that starts at r0 at time to. Because of the external driving, it depends on t and t0 separately not only on the difference t - t0. The Jacobian matrix of this trajectory with respect to the initial conditions is

J( ro, to; t) =

( dyx dyx

d xo dvo

dy„ dy„

1 d xo dvo /

All derivatives on the right hand side of (32) are to be evaluated at (ro, to; t).

Reactive trajectories are those that have an initial velocity Vi > V*, for a critical velocity V* measured at x0. Each reactive trajectory will cross the moving DS at a time tc(Avi) and with a velocity vc. If the crossing time decays monotonically from tc(AV*) = to to tc(to) = 0, the inverse function Avi(tc) or vi(tc) = v*(0) + ¿vi(tc) can be obtained. For any crossing time tc > 0, there is a unique initial velocity vi that will lead to a crossing at the given time.

The population of the product region Ax > 0 at time tf is therefore, as in Eq. (20),

Pp(tc) = / q(v) d V,

Jvi(tc)

and the flux across the moving surface is

FM(tc) = dPP dt c

= -q(v'(t c)) dTc

= -q(vi(t c))

dt c '

This result is positive because the initial velocity is a decreasing function of the crossing time.

The flux can also be evaluated directly from the flux integral (22)

FM(tc) = dAv Av Ptc(Ax = 0, Av), (35)

where pt(Ax, Av) is the density of trajectories in phase space at time t. Initially, this density is (Eq. (23))

p0(Ax, Av) = ¿(Ax - Ax0) q(v*(0) + Av). (36)

At later times, it can be obtained from Eq. (24)

Pt (Ax, Av) = eYt p0(y*(Ax, Av, t ;0),yv(Ax, Av, t;0)). (37)

Here, we have used the general notation for the flow of the equation of motion. The exponential accounts for the shrinkage of phase space volume and the corresponding increase in density. It is the same as in the harmonic case: In general, the flow of a differential equation ii = f (a) leads to a stretching of volume whose rate is the divergence of the vector field f. For Eq. (16), this rate is constant -y, so that over time t all volumes will shrink by a factor e"rt.

The flux integral formula (22) now reads

FM(tc) = tc

dAv Avd(yx(o, Av, tc;Q) - Axo)

x q(v *(o) + y„(o, Av, t c;o)).

The ô function requires that the trajectory that reaches Ax = 0, Av at time tc must have started at Ax0 at time 0. It produces a single contribution to the integral at velocity Avc(tc), so that

FM(tc) = eytc q(vi(o) + Av,(tc))

Avc(tc)

where yv(0, Av,tc;0) = Av,(tc) and the subscript tc indicates that the derivative is to be evaluated at (0,Avc(tc), tc;0). Similarly, a subscript 0 will be used to require evaluation at (Ax0, Av,(tc), 0; tc). These subscripts indicate derivatives of the flow taken along the trajectory from (Ax0, Av,(tc)) at t = 0 forward in time to (0, Avc(tc)) at t = tc (subscript 0) and along the same trajectory backward in time (subscript tc).

To verify that the flux integral (39) gives the same result as (34) that was obtained from the product population, it must be shown that

Avc(tc)

To this end, first note that Av,(tc) is defined by the condition

y>x(Ax0, Av,(tc), 0; t c ) = 0. Differentiating this condition with respect to tc gives dAv,-

dyx dAv

dyx dt

Now 5yx/5t is the velocity of the trajectory at the end point. The second term in Eq. (41) is therefore Avc(tc). With this result, the condition (40) simplifies to

dy x dAv

dy x dAv

Under the given assumptions on the geometry, the derivative on the left hand side is negative: A trajectory that arrives at the DS with larger velocity must have started further away, i.e., at smaller Ax(0).

dAvo ltc

= ey c

dAvo ltc

= ey c

The derivatives occurring in Eq. (42) are elements of the Jacobian matrices

J|fc = J(0, Avc(tc), ic;0) =

dy* dy* \

d x tc dv tc

dy„ dy„

V d x tc dv tc /

J|o = J(Axo, Av,-(t c), 0; t c) =

d x dy„

dy* dv dy„ dv

respectively. Because these matrices describe variations around the same trajectory, taken forwards and backwards in time, they are inverse to each other. Formally, this can be shown by taking derivatives of the flow property

®(® (Jo, 0; tc), tc; 0) = To for all ro,

which says that propagating an arbitrary phase space point r0 forward in time by tc and back again will return the original point.

By the well known formula for the inverse of a 2 x 2 matrix, it follows that

Jltc = (J|0)-1

det J|0

dv 0 dy„

dy* \ dv 0 dy

so that

dy * dAv

V d x 0 d x 0 /

detJ|0 dAv

- n-ytc

det J|o = e

is the factor by which phase space volumes shrink during time tc. This proves the condition (42) and therefore the equality of the two flux formulas.

25 20 15 10 5

20 15 10 5

20 15 10 5 0

o2 = {5,10}

h-—.—

i23 = {2.5,5,7.5}

FIG. 4. Time dependence of the scaled logarithm of the reactant population, -ln[PR(t)-PR(M)],for (top), n2 (middle), and n3 (bottom) with Hf = 5 for all panels. Values of the anharmonic parameter are e e {1,2,4,6,8,10}. The slope of each dashed line is the barrier crossing rate kf. The color of each line corresponds to the respective e value. In all panels, parameters are y = 1 and f = 0.

C. Dynamics near the TS

The TS trajectory is a moving saddle point and thus trajectories in the neighborhood of r* can be described by a linearization of the equations of motion. In the phase space vector relative coordinate Ar = (Ax, Au), this linearization is given by

J(t) =

Air = J(t) Ar,

+ 3e(x*(t) - E(t))2 -y

is the Jacobian of Eq. (16). The asymptotic decay rate of PR(t) is determined by the behavior of trajectories with initial conditions close to the stable manifold. For an ensemble of trajectories constituting an initial phase space density p0, trajectories that emanate close to vC (the stable manifold at t = 0) will persist in the neighborhood where (43) is valid

for long times. The decay of these trajectories determines the reaction rate.

The stretching and compression of phase space about a PO is known to dictate escape rates in conservative49-51 and dissipative systems.48 When J(t) is periodic in systems of the form of Eq. (43), the rate of deformation in the linearized phase space can be quantified through calculation of the Floquet exponents.66

To classify the stability of Ar*, we consider the dynamics of a perturbation vector A^(t). The equation of motion (43) is linear in A^(t) and thus it satisfies

Act = J(t ) Ac, Ac(0) = I,

where I is the 2 x 2 identity matrix. The principal fundamental matrix solution over one period of the driving is the monodromy matrix

= /ac(1)çt) ac(2)(T)\ m = |A<t(1)(T) act(2)(T)) •

(a) 8 7 6

(c) 8 7

(e) 8 7 6

: iil ii = 10 <

: 7=1 *"o " °Q=7o 1

---- """" fi=5 -

Sr °fi—3

•___i---

- —-- — • n=i

harmonic

: n2 : 7=1

■ °

o o — *Tf==2 "

x harmonic

(b) 8 7 6

(d) 8 7

(f) 8 7

: fix " 7=2 ^r ............C.... 0 _ °fi=7o ' O 0fl=50 j o °fi =

- —~..... o ____—•—* __•--"TT=2 .

■ m— __ m —*■ ii=1

¿r^ harmonic

. max \E(t)\ » 0.77

- o Q— o — g -! —•-~~*il=2 .

- ^^ harmonic

FIG. 5. The barrier crossing rates of systems following the equations of motion (2) as a function of the anharmonic parameter e for various frequency sets driving frequencies and values of friction y, as denoted in each panel. The circles denote the rates kf calculated from the time evolution of PR(t) through numerical simulation and correspond to the dashed lines in Fig. 4. The solid curves are the rates predicted by the difference in the characteristic Floquet exponents ¡ia -¡is of the corresponding TS trajectory.

A fundamental matrix solution Aa(t) of (43) at some later time t + kT, for k = 1,2,3, ..., can be obtained as

A a(t + kT) = Mk A a(t), (47)

through repeated operation by the monodromy matrix.

The eigenvalues ms,u of M are the Floquet multipliers. The Floquet exponents

Ms, u = T ln |ms,u I (48)

quantify the stability of Ar* and give the rate of expansion or contraction of the perturbation of per unit time.67-69 The TS trajectory has both an unstable mu > 0 and a stable ms < 0 exponent which correspond to stretching and contraction of the initial perturbation in the directions of the unstable and stable manifolds, respectively.

For an arbitrary time interval of length T, trajectories that cross the DS in this interval form a strip in the phase plane. Trajectories that cross the DS in the next following time interval T form a similar strip that is closer to the stable manifold. In the region where the linearized system is valid, the phase space density is constant. The flux of trajectories

through the DS in a given time interval is proportional to the width of the strip that contains these trajectories. During sequential periods, this width decreases by a factor e~(pu~Ms)t. From this it follows that, up to periodic modulation, the flux must decay as e~(Pu~Ms)t and the barrier crossing rate is

kf = Mu - Ms, (49)

which expresses the reaction rate in terms of the characteristic Floquet exponents of the TS trajectory. Equation (49) generalizes Eq. (30) for the case of an anharmonic barrier.

V. NUMERICAL RESULTS AND COMPARISON WITH THEORY

The reaction rate of (2) was calculated by simulating ensembles of n = 108-109 trajectories for various sets of parameters {Qy,e,<r} and following the survival probability of PR as a function of time. A Runge-Kutta-Maruyama scheme70 was implemented to perform the integration. In the absence of noise = 0), this algorithm is the well-known fourth-order Runge-Kutta method. For all numerical

Ö • i—i

s-i O CD

• 1—I

CD S-i

35 30 25 20 15 10 5 0 35 30 25 20 15 10 5 0 10

(a) - £=0 - fi = 10 - --- n=5 ; ■■■■■ 0=2 _ j 7=0.01 ~

..........................

......... ........

■ 7=0.1

................. .................

35 faß 30 •d 25 CG 03 O 20 ■(c) ; 7=1 Sir >V - » <

1 tr O 15 CD Sh 10

........1 ........1 ........

FIG. 6. The percentage of trajectories that recross the moving dividing surface attached to the DTS trajectory as a function of noise strength t with (a) y = 0.01, (b) y = 0.1, (c) y = 1, and (d) y = 3 for single-frequency (fti) driving and various values of e and Q. The black vertical line (solid) marks the noise strength when the fluctuation-dissipation theorem is obeyed.

simulations, non-dimensional parameters were used by choosing units such that the barrier frequency wb and driving amplitude are unity. Each trajectory was given an initial position x0 = -0.1 (in the reactant region) and v0 was sampled from a Boltzmann distribution with kBT = 1. The choice of initial conditions is arbitrary as the asymptotic decay rate of PR(t) is independent of the choice of initial distribution, suffice that there is enough density about the stable manifold such that a rate exists.48

The ensemble of n trajectories was evolved through the equations of motion (2). The normalized reactant population was calculated at each time step in the integration scheme. An indicator function was employed to follow the state evolution of each trajectory,

Hr[X(Î )] =

(0, x(t) > x\t), 1, x(t) < x\t),

where x*(t) is the configuration space projection of the TS trajectory. If for a specific trajectory i, x¡(t) > x*(t) that trajectory is in the product state and is not counted in the reactant population at time t. The instantaneous normalized population of the reactant region can be found by summing over all n trajectories and then normalizing by a factor 1/n,

PR(t) = - Y! är[xi(t)].

n i—i

Trajectories can only exist in one of two states, reactant or product, and so the normalized product population PP = 1 - PR.

As shown in Fig. 4, the scaled logarithm of the normalized reactant population, -ln [PR(t) - PR(ro)], is approximately linear in time after an initial transient section implying a

first-order rate process. The asymptotic reaction rate kf can thus be found as the slope of the scaled logarithmic curve in the long-time limit. Periodic modulation in the decay of PR(t) was found to become more prominent for low frequency driving (Qf < 2). In these cases, the global exponential rate was calculated as an average over these modulations.

A comparison between the rates calculated from numerical simulation and rates predicted by Eq. (49) is shown in Fig. 5. For all frequency sets and parameter values, agreement is observed. Underdamped (y < 2), overdamped (y > 2), and critically damped (y = 2) regimes of a corresponding harmonic well were considered. Agreement between the rates persists over all ranges of damping. For high frequency driving (Qf > wb), the exponential rate can be averaged over several periods of driving and modulations in the decay are minimal, as illustrated in Fig. 4. Periodic modulations in the decay of PR(t) are prominent for low driving frequencies (Qf ~ wb) and the integration of n = 108 trajectories resulted in reaching the numerical asymptote Pr(to) at times less than the period of the external driving. In those cases for which the integration time was insufficient to sample the asymptotic region, a larger number of trajectories (n = 109) were integrated. Each trajectory was integrated to a final time of at least tf = 15 and, as shown in Fig. 4, PR(^) is reached well before the end of this sampling window. Increasing the number of trajectories by an order of magnitude resulted in improved convergence of the scaled logarithmic population and marginally better agreement between the compared methodologies, as shown in Fig. 5(e) for Q = 1. The agreement between the two methods is illustrated in Fig. 5(d) for the smaller, non-unity, driving amplitude case of The decreased driving amplitude leads to a decrease in the reaction rate.

FIG. 7. Time dependence of the scaled logarithm of the reactant population for systems with single-frequency (ft1) periodic and thermal driving for y = 0.01 (top), y = 0.1 (middle), and y = 1 (bottom). The color of each line corresponds to a specific r value. The decay for systems with various anhar-monicities e e {1,3,10} is shown and denoted in each panel. The fundamental driving frequency is £f = 5 for all panels. For visual clarity, each curve is truncated at a point where the data became noisy.

VI. CHARACTERIZING NOISY REACTIONS WITH THE NOISE-FREE GEOMETRY

In systems in which the strength of an external driving force dominates over that of the thermal driving, statistical quantities can be approximated by those of a corresponding purely deterministically driven system. For thermally induced reactions, Lehmann, Reimann, and Hanggi63-65 have shown that in the overdamped (large-y) regime, when a chemical reaction is forced by a periodic field the reaction rate is determined in part by the geometry of periodic trajectories in the purely deterministic phase space. This work was later extended to cases with different scaling behaviors between the strength of thermal activation and the strength of the external

field.71-73

Our goal here is to develop a minimalist theory, applicable at the limit where the magnitude V2< of a noise sequence %a(t) is a small enough perturbation to the periodic driving E(t) that the TS trajectory of the noiseless system (the periodic

orbit) gives rise to a DS with minimal recrossings. This deterministic TS trajectory (DTS trajectory) does not solve the equations of motion (2) with a non-zero value of <. We therefore distinguish the DTS trajectory from the true TS trajectory of the noisy system (that we do not compute in this work).

A principal assumption for the use of the noise-free geometry is that the phase space density of the thermal system, and its time-dependence, is approximately that of the deterministic system, i.e., pt(Axa, Ava) « pt(Ax, Av). As shown in Fig. 2, for small values of <, the geometry of the thermal system is similar to that of its deterministic counterpart. The rate theory developed in Sec. IV C for the deterministic system can therefore be applied. This is advantageous in applications such as in comparisons with experiments in which the exact noise sequence is not known.

Thermal systems in which the fluctuation-dissipation theorem (FDR) is not obeyed due to energy dissipation constitute non-equilibrium processes. Formal treatments of fluctuation-response in periodically forced systems by Teramoto, Harada, and Sasa74,75 provide insight into the rate of energy dissipation in such systems. Green et al.76 have shown that the rate of energy dissipation is directly related to the dynamical entropy of the system. To realize non-equilibrium conditions in the present model reaction, the damping constant y is held constant and the strength of the thermal fluctuations < is increased up to the point where the FDR is satisfied. If the initial velocities are drawn from a Boltzmann ensemble with kBT = 1 (in dimensionless units), this is the case at < = y. If < < y, the thermal bath is at a lower temperature than that of the distribution of initial velocities. In this case, the temperature of the ensemble of reactants will be continuously cooled by the thermal bath. Interestingly, although the flux over population rate is well defined numerically, as exhibited in Fig. 7, the ensemble in the reactant region does not satisfy the equilibrium distribution. This requires a reconsideration of TST in the context of nonequilibrium processes as below.

The percentage of thermal trajectories that recross the DS attached to the DTS trajectory is shown in Fig. 6 for varying noise strengths < and constant dissipation rates. As shown in Figs. 6(a) and 6(b), a minimal number of recrossings occur below and up to the FDR threshold for small values of y. For the y = 1 case, shown in Fig. 6(c), trajectories persist around the BT for long times, leading to a larger number of recrossings than observed for smaller dissipation rates. For the overdamped dynamics (y = 3), shown in Fig. 6(d), the deterministic DS identifies reactive trajectories adequately only for weak thermal driving (small < ) and strong anharmonic-ity. As the harmonic limit is approached or in equilibrium systems the superimposed DS becomes very poor.

The decay of the scaled logarithm of the normalized reactant population, as calculated with the superimposed deterministic DS, is shown in Fig. 7 for various parameter values. Over all friction regimes, the population decay of the systems with additional thermal driving follows that of its deterministic counterpart if the noise strength < is sufficiently low. For y = 1, when the strength of the thermal driving approaches that of the FDR, a decrease in the reaction rate is observed. The data presented in Fig. 7 become highly

FIG. 8. The barrier crossing rates of systems following the equations of motion (2) as a function noise strength c. The rates calculated using the DS attached to the DTS trajectory for single-frequency (ftj) driving and various values of e, y, and Q are shown as circles. The horizontal lines (solid) denote the rates given by the Floquet exponents of the corresponding DTS trajectory and are colored according to a respective e value. The black vertical lines (solid) denote the noise strength where the fluctuation-dissipation theorem is obeyed.

(a) 9 8

7=0.01 0=5 ;

- e = 10 ' O ( ; e=5 1 -

. e = l

; 7=0.1 0=5

: ° e=5 0 ' - ° .=3 0 * ! '

. * e=l "

■ 7 = 1 0=5 -

' o n -

- E ° e=5 --_

- e—3 o ' ------

■ "........^........

(b) 9 8

7=0.01 0 = 10 ;

- * e=10 -

' ° Q

" * £=1

7=0.1 0=10

- e=10 " . o o < e=5 ) .

; e=3 0

; ' £=1 1 ------- "........................ ........1 1 ■ "

_......1 ........1 ........1 1 1 1 111II ' ' _ 7 = 1 ii = 10 ; ■ ° O _ ■

e=5 e=3

oscillatory at long times due to recrossings of the DS. For visual clarity, each data series has been truncated to remove this noisy tail. As observed in Figs. 4 and 7, for short times (t g 0.3), the decay of PR is non-exponential, signifying temporally global non-RRKM kinetic behavior. (Note that RRKM refers to Rice-Ramsperger-Kassel-Marcus theory.) We obtain the rate from the long-time asymptotic decay of PR, which is representative of kinetic experiments in which the concentration of a reactant species is directly measured over time.77

The thermal rates calculated using the DTS trajectory are shown in Fig. 8. As expected by the minimal number of recrossings shown in Fig. 6, stability analysis of the DTS can produce an excellent approximation to the rate in thermal environments. Through calculation of the error between the numerically calculated rate with included noise and the rate given by the Floquet exponents of the DTS trajectory, the extent of applicability of the noise-free geometry can be quantified. This error is <3% at y = 0.01 over all parameter values. It is < 1% for e e {5,10}. Increasing the dissipation by an order of magnitude (y = 0.1) results in the same general trends, with all errors generally less than 5%. The exceptions

occur at the noise strength where the FDR is obeyed (a = 0.1) at e e {1,3} and Q = 5 for which the error « 20%. For y = 1.0 and Q = 10, all calculated errors are less than or on the order of 20%, increasing monotonically as a function of a. As illustrated in Fig. 8(e), at lower-frequency driving (Q = 5) and large noise (a = 1), the error is between 30%-50%. This suggests a practical upper bound to the applicability of the noise-free geometry in estimating the reaction rates in the presence of noise. Although not shown, for overdamped dynamics, stability analysis of the DTS gives an accurate approximation to the rate only in non-equilibrium small noise regimes.

The calculated errors are on the order of the error expected from application of variational transition state theory (VTST).5 The presented methodology is advantageous over VTST as it does not require the integration of large numbers of trajectories or a flux minimization procedure. Thus, stability analysis of the DTS trajectory offers a simple rate calculation methodology that can be readily applied, in weak thermal environments, to driven chemical reactions with only prerequisite knowledge of the geometry of the energy surface and the functional shape of the driving waveform.

VII. CONCLUSIONS

We have shown that in a model chemical reaction subjected to the influence of forcing from a temporally periodic external field, a recrossing-free dividing surface can be constructed over an unstable periodic orbit in the region of a moving energetic barrier top. This forms the basis for future work on specific driven chemical reactions that can be represented by a single collective variable for the reaction coordinate under nonequilibrium conditions. Potential targets include both substitution and isomerization reactions in which the governing multi-dimensional energy surface can be parameterized by a single collective degree of freedom. Other possible targets include mechanochemical reactions and stimuli-responsive assembly mechanisms when the reaction rate is dictated predominantly by geometric properties about the moving dividing surface. Generally, force-modified and temperature-modulated energy surfaces, deforming under the influence of temporally varying forces, motivate the development of rate methodologies that go beyond the simplistic equilibrium arguments intrinsic to classical equilibrium transition state theory.

The no-recrossing surface constructed here has been shown to persist for strongly anharmonic barriers subjected to single-mode and multi-mode driving waveforms. A formally exact rate theory has been developed based on the flux of reactive trajectories through this recrossing-free surface. It rectifies the principal criterion of transition state theory for periodically driven chemical reactions.

To circumvent computationally taxing numerical calculations of the reactive flux through this surface, a rate theory has been developed based on the stability of the dividing surface. Strong agreement was observed between the rate predicted by the Floquet exponents of a trajectory defining the phase space evolution on the dividing surface, and the rate calculated from simulation of a large ensemble of trajectories. Thus, in a periodically driven chemical reaction, the asymptotic decay rate of an initial distribution of reactants can be extracted directly from the stability of the time-varying dividing surface irrespective of the dynamics of the reactive population.

Use of the noise-free geometry to approximate the corresponding structure of a driven thermal system has been shown to give an excellent approximation to the optimal dividing surface if the magnitude of the oscillating force is large compared with that from the thermal environment. For thermally activated processes, the stability exponents of the purely periodically driven system can thus be used to predict the reaction rates without an explicit treatment of the thermal dynamics. Extensions of this work to include an explicit treatment of the noise, systems with structured solvents environments,78,79 and systems displaying fluctuating rates,80 are possible next steps, and ones which we are currently pursuing.

ACKNOWLEDGMENTS

This work has been partially supported by the Air Force Office of Scientific Research through Grant No. FA9550-12-1-0483. Travel between partners was partially supported through the People Programme (Marie Curie Actions) of the

European Union's Seventh Framework Programme FP7/2007-2013/ under REA Grant Agreement No. 294974.

!K. Yamanouchi, "The next frontier," Science 295, 1659-1660 (2002).

2B. J. Sussman, D. Townsend, M. Y. Ivanov, and A. Stolow, "Dynamic Stark control of photochemical processes," Science 314, 278-281 (2006).

3S. Kawai and T. Komatsuzaki, "Quantum reaction boundary to mediate reactions in laser fields," J. Chem. Phys. 134, 024317 (2011).

4D. G. Truhlar, W. L. Hase, and J. T. Hynes, "Current status of transition-state theory," J. Phys. Chem. 87, 2664-2682 (1983).

5D. G. Truhlar and B. C. Garrett, "Variational transition state theory," Annu. Rev. Phys. Chem. 35, 159-189 (1984).

6W. H. Miller, "Beyond transition-state theory: A rigorous quantum theory of chemical reaction rates," Acc. Chem. Res. 26, 174 (1993).

7D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, "Current status of transition-state theory," J. Phys. Chem. 100, 12771-12800 (1996).

8H. Waalkens, R. Schubert, and S. Wiggins, "Wigner's dynamical transition state theory in phase space: Classical and quantum," Nonlinearity 21, R1 (2008).

9R. Hernandez, T. Bartsch, and T. Uzer, "Transition state theory in liquids beyond planar dividing surfaces," Chem. Phys. 370, 270-276 (2010).

10E. Pollak and P. Pechukas, "Transition states, trapped trajectories, and classical bound states embedded in the continuum," J. Chem. Phys. 69, 1218 (1978).

"E. Pollak and P. Pechukas, "Unified statistical model for 'complex' and 'direct' reaction mechanisms: A test on the collinear H + H2 exchange reaction," J. Chem. Phys. 70, 325-333 (1979).

12P. Pechukas and E. Pollak, "Classical transition state theory is exact if the transition state is unique," J. Chem. Phys. 71, 2062 (1979).

13E. Pollak, M. S. Child, and P. Pechukas, "Classical transition state theory: A lower bound to the reaction probability," J. Chem. Phys. 72, 1669 (1980).

14R. Hernandez and W. H. Miller, "Semiclassical transition state theory. A new perspective," Chem. Phys. Lett. 214, 129-136 (1993).

15R. Hernandez, "A combined use of perturbation theory and diagonalization: Application to bound energy levels and semiclassical rate theory," J. Chem. Phys. 101, 9534-9547 (1994).

16T. Uzer, C. Jaffe, J. Palacian, P. Yanguas, and S. Wiggins, "The geometry of reaction dynamics," Nonlinearity 15, 957 (2002).

17N. De Leon, M. A. Mehta, and R. Q. Topper, "Cylindrical manifolds in phase space as mediators of chemical reaction dynamics and kinetics. I. Theory," J. Chem. Phys. 94, 8310-8328 (1991).

18C. Li, A. Shoujiguchi, M. Toda, and T. Komatsuzaki, "Definability of no-return transition states in the high-energy regime above the reaction threshold," Phys. Rev. Lett. 97, 028302 (2006).

19H. Waalkens and S. Wiggins, "Direct construction of a dividing surface of minimal flux for multi-degree-of-freedom systems that cannot be recrossed," J. Phys. A: Math. Gen. 37, L435-L445 (2004).

20G. S. Ezra, H. Waalkens, and S. Wiggins, "Microcanonical rates, gap times, and phase space dividing surfaces," J. Chem. Phys. 130, 164118 (2009).

21G. S. Ezra and S. Wiggins, "Phase-space geometry and reaction dynamics near index 2 saddles," J. Phys. A: Math. Theor. 42, 205101 (2009).

22H. Teramoto, M. Toda, and T. Komatsuzaki, "Dynamical switching of a reaction coordinate to carry the system through to a different product state at high energies," Phys. Rev. Lett. 106, 054101 (2011).

23M. Inarrea, J. F. Palacian, A. I. Pascual, and J. P. Salas, "Bifurcations of dividing surfaces in chemical reactions," J. Chem. Phys. 135,014110 (2011).

24A. Allahem and T. Bartsch, "Chaotic dynamics in multidimensional transition states," J. Chem. Phys. 137, 214310 (2012).

25R. S. MacKay and D. C. Strub, "Bifurcations of transition states: Morse bifurcations," Nonlinearity 27, 859 (2014).

26U. giftgi and H. Waalkens, "Reaction dynamics through kinetic transition states," Phys. Rev. Lett. 110, 233201 (2013).

27N. Elsner, C. P. Royall, B. Vincent, and D. R. E. Snoswell, "Simple models for two-dimensional tunable colloidal crystals in rotating ac electric fields," J. Chem. Phys. 130, 154901 (2009).

28S. Jäger and S. H. L. Klapp, "Pattern formation of dipolar colloids in rotating fields: Layering and synchronization," Soft Matter 7, 6606 (2011).

29A. Prokop, J. Vacek, and J. Michl, "Friction in carborane-based molecular rotors driven by gas flow or electric field: Classical molecular dynamics," ACS Nano 6, 1901-1914 (2012).

30F. Ma, D. T. Wu, and N. Wu, "Formation of colloidal molecules induced by alternating-current electric fields," J. Am. Chem. Soc. 135, 7839-7842 (2013).

31P. Lidström, J. Tierney, B. Wathey, and J. Westman, "Microwave assisted organic synthesis—A review," Tetrahedron 57, 9225-9283 (2001).

32M. Platkov and M. Gruebele, "Periodic and stochastic thermal modulation of protein folding kinetics," J. Chem. Phys. 141, 035103 (2014).

33G. Loget and A. Kuhn, "Electric field-induced chemical locomotion of conducting objects," Nat. Commun. 2, 535 (2011).

34M. E. Corrales, J. González-Vázquez, G. Balerdi, I. R. Solá, R. de Nalda, and L. Bañares, "Control of ultrafast molecular photodissociation by laser-field-induced potentials," Nat. Chem. 6, 785-790 (2014).

35S. Saha and J. F. Stoddart, "Photo-driven molecular devices," Chem. Soc. Rev. 36, 77-92 (2007).

36J. Michl andE. C. H. Sykes, "Molecular rotors and motors: Recent advances and future challenges," ACS Nano 3, 1042-1048 (2009).

37C. Zazza, G. Mancini, G. Brancato, and V. Barone, "In silico study of molecular-engineered nanodevices: A lockable light-driven motor in dichloromethane solution," J. Phys. Chem. Lett. 4, 3885-3890 (2013).

38W. R. Browne and B. L. Feringa, "Making molecular machines work," Nat. Nanotechnol. 1, 25-35 (2006).

39E. R. Kay, D. A. Leigh, and F. Zerbetto, "Synthetic molecular motors and mechanical machines," Angew. Chem., Ind. Ed. 46, 72-191 (2007).

40V. Balzani, A. Credi, and M. Venturi, "Molecular devices and machines," Nano Today 2, 18-25 (2007).

41H. Meng and G. Li, "Reversible switching transitions of stimuli-responsive shape changing polymers," J. Mater. Chem. A 1, 7838-7865 (2013).

42B. L. Feringa, "Nanotechnology: In control of molecular motion," Nature 408, 151-154 (2000).

43D. A. Leigh, J. K. Wong, F. Dehez, and F. Zerbetto, "Unidirectional rotation in a mechanically interlocked molecular rotor," Nature 424,174-179 (2003).

44S. P. Fletcher, F. Dumur, M. M. Pollard, and B. L. Feringa, "A reversible, unidirectional molecular rotary motor driven by chemical energy," Science 310, 80-82 (2005).

45D. Horinek and J. Michl, "Surface-mounted altitudinal molecular rotors in alternating electric field: Single-molecule parametric oscillator molecular dynamics," Proc. Natl. Acad. Sci. U. S. A. 102, 14175-14180 (2005).

46M. Klok, N. Boyle, M. T. Pryce, A. Meetsma, W. R. Browne, and B. L. Feringa, "MHz unidirectional rotation of molecular rotary motors," J. Am. Chem. Soc. 130, 10484-10485 (2008).

47G. 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).

48G. T. Craven, T. Bartsch, and R. Hernandez, "Communication: Transition state trajectory stability determines barrier crossing rates in chemical reactions induced by time-dependent oscillating fields," J. Chem. Phys. 141, 041106(2014).

49L. P. Kadanoff and C. Tang, "Escape from strange repellers," Proc. Natl. Acad. Sci. U. S. A. 81, 1276-1279 (1984).

50R. T. Skodje and M. J. Davis, "Statistical rate theory for transient chemical species: Classical lifetimes from periodic orbits," Chem. Phys. Lett. 175, 92-100 (1990).

51P. Gaspard, Chaos, Scattering and Statistical Mechanics (Cambridge University Press, 1998), Vol. 9.

52J. R. Green, T. S. Hofer, R. S. Berry, and D. J. Wales, "Characterizing molecular motion in H2O and H3O+ with dynamical instability statistics," J. Chem. Phys. 135, 184307 (2011).

53A. E. Orel and W. H. Miller, "Collision induced absorption spectra for gas phase chemical reactions in a high power IR laser field," J. Chem. Phys. 72, 5139-5144(1980).

54S. Kawai, A. D. Bandrauk, C. Jaffé, T. Bartsch, J. Palacián, and T. Uzer, "Transition state theory for laser-driven reactions," J. Chem. Phys. 126, 164306 (2007).

55R. G. Mullen, J.-E. Shea, and B. Peters, "Transmission coefficients, committors, and solvent coordinates in ion-pair dissociation," J. Chem. Theory Comput. 10, 659-667 (2014).

56R. G. Mullen, J.-E. Shea, and B. Peters, "Communication: An existence test for dividing surfaces without recrossing," J. Chem. Phys. 140, 041104 (2014).

57D. G. Truhlar and B. C. Garrett, "Multidimensional transition state theory and the validity of Grote-Hynes theory," J. Phys. Chem. B 104, 1069-1072 (2000).

58T. Bartsch, R. Hernandez, and T. Uzer, "Transition state in a noisy environment," Phys. Rev. Lett. 95, 058301-01—058301-04 (2005).

59T. Bartsch, T. Uzer, and R. Hernandez, "Stochastic transition states: Reaction geometry amidst noise," J. Chem. Phys. 123, 204102 (2005).

60T. Bartsch, T. Uzer, J. M. Moix, and R. Hernandez, "Identifying reactive trajectories using a moving transition state," J. Chem. Phys. 124, 24431001—244310-13 (2006).

61F. Revuelta, T. Bartsch, R. M. Benito, and F. Borondo, "Communication: Transition state theory for dissipative systems without a dividing surface," J. Chem. Phys. 136, 091102 (2012).

62T. Bartsch, F. Revuelta, R. M. Benito, and F. Borondo, "Reaction rate calculation with time-dependent invariant manifolds," J. Chem. Phys. 136, 224510 (2012).

63J. Lehmann, P. Reimann, and P. Hanggi, "Surmounting oscillating barriers," Phys. Rev. Lett. 84, 1639-1642 (2000).

64J. Lehmann, P. Reimann, and P. Hanggi, "Surmounting oscillating barriers: Path-integral approach for weak noise," Phys. Rev. E 62,6282-6303 (2000).

65J. Lehmann, P. Reimann, and P. Hanggi, "Activated escape over oscillating barriers: The case of many dimensions," Phys. Status Solidi B 237, 53-71 (2003).

66C. Skokos, "On the stability of periodic orbits of high dimensional autonomous Hamiltonian systems," Physica D 159, 155-179 (2001).

67P. Cvitanovic, R. Artuso, R. Mainieri, G. Tanner, and G. Vattay, Chaos: Classical and Quantum (ChaosBook.org, Niels Bohr Institute, Copenhagen, 2012).

68R. P. Boland, T. Galla, and A. J. McKane, "Limit cycles, complex Floquet multipliers, and intrinsic noise," Phys. Rev. E 79, 051131 (2009).

69F. L. Traversa, M. Di Ventra, and F. Bonani, "Generalized Floquet theory: Application to dynamical systems with memory and Bloch's theorem for nonlocal potentials," Phys. Rev. Lett. 110, 170602 (2013).

70A. Naess and V. Moe, "Efficient path integration methods for nonlinear dynamic systems," Probab. Eng. Mech. 15, 221-231 (2000).

71R. S. Maier and D. L. Stein, "Noise-activated escape from a sloshing potential well," Phys. Rev. Lett. 86, 3942-3945 (2001).

72M. I. Dykman, B. Golding, and D. Ryvkine, "Critical exponent crossovers in escape near a bifurcation point," Phys. Rev. Lett. 92, 080602 (2004).

73M. I. Dykman and D. Ryvkine, "Activated escape of periodically modulated systems," Phys. Rev. Lett. 94, 070602 (2005).

74T. Harada and S.-I. Sasa, "Equality connecting energy dissipation with a violation of the fluctuation-response relation," Phys. Rev. Lett. 95, 130602 (2005).

75H. Teramoto and S.-I. Sasa, "Microscopic description of the equality between violation of fluctuation-dissipation relation and energy dissipation," Phys. Rev. E 72, 060102 (2005).

76J. R. Green, A. B. Costa, B. A. Grzybowski, and I. Szleifer, "Relationship between dynamical entropy and energy dissipation far from thermodynamic equilibrium," Proc. Natl. Acad. Sci. U. S. A. 110, 16339-16343 (2013).

77I. N. Levine, Physical Chemistry (McGraw-Hill, 2002).

78G. T. Craven, A. V. Popov, and R. Hernandez, "Stochastic dynamics of penetrable rods in one dimension: Occupied volume and spatial order," J. Chem. Phys. 138, 244901 (2013).

79G. T. Craven, A. V. Popov, and R. Hernandez, "Structure of a tractable stochastic mimic of soft particles," Soft Matter 10, 5350-5361 (2014).

80S. W. Flynn, H. C. Zhao, and J. R. Green, "Measuring disorder in irreversible decay processes," J. Chem. Phys. 141, 104107 (2014).