Influence of external driving on decays in the geometry of the LiCN isomerization

The motion of atoms in chemical reactions can often be described by classical mechanics on a Born–Oppenheimer potential separating reactants from products through a rank-1 saddle. The corresponding barrier has one unstable direction along the reaction coordinate, and the remaining orthogonal modes along the stable direction. The framework of transition state theory (TST) provides both qualitative and quantitative description of reaction rates using the sum of the flux through a particular dividing surface (DS).^1–7 In the most naive case, the DS is a plane located at the saddle and parallel to the orthogonal modes, providing the usual qualitative, but not exact, Arrhenius-like rates.^1,2,4,8,9 More generally, the DS can be extended to a fully recrossing-free surface in phase space giving rise to exact rates.^10–21 

In the phase space associated with the barrier region of such systems, the most relevant object is the normally hyperbolic invariant manifold (NHIM),^15,22–25 which is a natural generalization of the time-independent two-dimensional periodic orbit dividing surface (PODS).^3,4 The NHIM contains all trajectories that are mathematically bound to the saddle forever and will, therefore, never leave to either the reactant or product sides. Hence, this intermediate structure lying between reactants and products is also often referred to as the transition state (TS). In multidimensional autonomous Hamiltonian systems, the NHIM can be constructed approximately up do a desired order using normal form expansions^26–35 or numerically by application up to desired accuracy of Lagrangian descriptors,^{24,25,36–38} the binary contraction method (BCM),³⁹ and machine learning algorithms.^25,40,41 Many of these methods also allow for the construction of a time-dependent NHIM when the system is driven through time-dependent potentials. Given the NHIM, a DS can then be attached^24,25,42,43 to it in the sense that it is anchored (or rooted) at the NHIM and lifted vertically in the momentum space. The resulting time-dependent DS is recrossing-free, at least, in its local neighborhood. Such recent advances are necessary for us to address the challenge of a driven chemical reaction such as done here on LiCN through external coupling of its dipole moment.

The dynamics transverse to the NHIM is unstable. Hence, any trajectory having a small deviation from it will depart to either the reactant or the product side. This departure of trajectories in a close neighborhood of the NHIM is associated with a rate that describes the decay of reactant population close to the TS.⁴⁴ One approach for obtaining these decay rates in a driven reaction was pursued using Floquet exponents and demonstrated for a one-dimensional system.⁴⁵ In this model, the NHIM contains only a single trajectory—viz., the TS trajectory^15,17,46—to which a DS can be attached. The decay rates obtained here correspond to the decay of the reactant population close to the TS in the sense of the usual flux that crosses it in the forward direction. We effectively impose absorbing boundary conditions by neglecting the long-time return of the trajectories after they have reached the product or reactant and restrict our use of the term global dynamics to refer to all motion before its crossing of this boundary to product or reactants. Thus, in the limiting case that the DS is non-recrossing in this global sense, these decay rates are the forward reaction rates, as was the case in the paradigmatic system of Ref. 45. To the extent, however, that the DS that we construct is local, then the decay rates are the direct rates associated with this particular barrier. When this barrier is rate determining, then this local decay is once again the overall rate. Even when it is not, the decay rate remains a useful quantity to describe the flow along the reaction pathway near an identified TS.

In multidimensional systems, infinitely many trajectories are located on the NHIM of a time-dependent driven barrier. Reactive trajectories may pierce the DS close to any of them. The full dynamics of bound trajectories inside the NHIM then needs to be considered when obtaining the decay rates of a reactant population close to the TS. This problem was recently addressed in Ref. 44 through three different approaches: (i) propagation of an appropriate ensemble of reactive trajectories and observing their piercing through the DS as a function of time, (ii) analysis of the structure of the stable and the unstable manifold close to the NHIM, or (iii) a multidimensional extension of the so-called Floquet method of Ref. 45 in which the Floquet exponents are used to obtain the decay rate. These approaches enabled the analysis of a two-dimensional model reaction with numerical precision.⁴⁴ It also confirmed that in a multidimensional system, the driving potential has a large influence on the dynamics of trajectories inside the NHIM and the corresponding decay of the reactant population close to the TS. The aim of this paper is the application of these methods to the case of a model with high fidelity to a real chemical system. Specifically, we address the (periodically driven) LiCN ⇌ LiNC isomerization reaction. The heavy mass of the Li allows us to consider this reaction in the classical limit, although there has been a lot of interest in the quantum versions of TST,^28,47,48 which can be revealed in the case of the corresponding HCN reaction,⁴⁹ whose classical reaction geometry was also addressed previously in the time-independent case.⁵⁰ The important 1:2 resonance in that system⁵¹ does not show up as strongly in the present LiCN system because of the difference in masses. This difference allows us to simplify the current analysis by constraining the CN vibration to be fixed.

For the LiCN isomerization reaction, analytical expressions for both the potential energy surface⁵² and the dipole surface⁵³ are known. The LiCN ⇌ LiNC isomerization reaction has received significant attention in the literature—e. g., in molecular dynamics simulations of LiCN in an argon bath,^54–57 an analysis of the geometric structures underlying the LiCN reaction dynamics,^58,59 and as a paradigmatic example for the control of other reactions using external fields.^60–62 Here, we apply a periodically driven external field on the molecule. This introduces a complexity not captured by the PODS in two-dimensional reactions because we now require a corresponding time-dependent DS to distinguish reactant and product regions, and it allows us to reveal the influence of driving on the phase-space resolved decay rates obtained for reactant populations close to the TS. These rates describe how trajectories deviate from the time-dependent NHIM—which is the multidimensional generalization of the TS-trajectory—in analogy to the characterization of the deviation from unstable periodic trajectories through the use of Lyapunov exponents.

This paper is organized as follows: In Sec. II A, the LiCN ⇌ LiNC isomerization reaction is introduced. The static and driven equations of motion in Secs. II A 1 and II A 2, respectively, are recapitulated using the known literature forms and parameters. Details of the theory and implementation of three different approaches⁴⁴ introduced recently for the numerical evaluation of the decay rates are included in the supplementary material. Because of their complementary advantages, each is used to numerically determine the decay rates of the reactant population close to the TS. Associated structures for the periodically driven LiCN ⇌ LiNC isomerization reaction are provided in Sec. III and compared with the corresponding results for the non-driven isomerization reaction. The dynamics of trajectories on the NHIM and the associated phase-space resolved decay rates are examined.

The three-atom molecule LiCN consists of carbon (C) and nitrogen (N) atoms forming a strongly bound cyanide anion, which is weakly bound to the lithium (Li) cation regiospecifically. That is, it is an isomeric molecule with two stable conformations, LiCN and LiNC,

Its motion can be described quasiclassically in the Born–Oppenheimer approximation, and analytical approximations for both the energy surface⁵² and the dipole surface⁵³ are known from the literature. It can be represented in a body-fixed reference frame (x′, z′), illustrated in Fig. 1 using Jacobi coordinates, where the z′-axis lies on the vector R pointing from the center of mass of the cyanide toward the lithium atom. Here, the relative distance between the nitrogen and the carbon atom is labeled r_CN and the angle between R and r_CN is referred to as ϑ = ∡(R, r_CN). Consequently, the regions near ϑ = 0 and ϑ = π correspond to the LiCN and LiNC isomers, respectively. When described in a body-fixed reference frame, the angle $α=∡(e^z,R)$ yields the overall orientation of the molecule relative to a space-fixed coordinate system (x, z), as also indicated in Fig. 1. Here, $e^z$ is the unit vector in the z-direction of the space-fixed coordinate system.

In the absence of time-dependent external fields, the energy is conserved in the isomerization process and the corresponding potential energy surface of the non-driven LiCN ⇌ LiNC isomerization reaction is independent of the overall orientation α of the molecule. It can be approximated by a potential energy surface V(R, ϑ) that only depends on R = |R| and ϑ, i. e., the intrinsic degrees of freedom of the molecule. The bond distance of the cyanide anion is held fixed at |r_CN| = 2.186 a₀, with a₀ being the Bohr radius, in keeping with earlier work showing that it has little effect on the dynamics.^63–65 

We use the potential energy surface V(R, ϑ) of Ref. 52 for our calculations. As shown in Fig. 2(a), there exist two local minima on the energy surface corresponding to the linear structures LiCN (ϑ = 0) and LiNC (ϑ = π). In between these two minima, at ϑ = 0.292π, a rank-1 saddle represents the bottleneck of the isomerization reaction, visualized by the equipotential lines in Fig. 2(a). At the energies above the barrier typical for reaction, the motion of the Li atom between the isomers appears as orbits around the cyanide.

The corresponding classical rotationless two degrees of freedom Hamiltonian, as applied to the LiCN reaction,^58,65 is

under the assumption that the canonical momentum p_α = 0 is fixed for non-rotating molecules. Using Hamilton’s formalism, the dynamics of the LiCN ⇌ LiNC isomerization reaction in the absence of external fields can be obtained for the state $γ(t)=𝜗,p𝜗,R,pRT$ by numerically integrating a set of first order differential equations for $γ̇(t)$⁠, i.e.,

In the present work, we adopt a fourth order Runge–Kutta algorithm with a fixed step size.⁶⁶ The first derivatives of the potential in Eq. (3) have been implemented analytically in C++.

The aim of this work is to reveal how external driving influences the decay rates of the reactant population close to the transition state (TS). Such external driving can be induced by a time-dependent homogeneous electric field along a fixed direction. We call this the z direction and write the driving term as

which couples to the molecule’s dipole moment μ. This makes the Hamiltonian of Eq. (3) time-dependent through the added term

to the potential energy.⁶⁷ To complete the equations, however, we must now also specify a continuous representation for the molecule’s dipole moment μ, the so-called dipole surface, accounting for the fact that the underlying electronic wave functions vary as a function of the Born–Oppenheimer coordinates.

Individual points of the dipole surface have earlier been obtained using self-consistent field (SCF) methods.⁵² Motivated by the success of Wormer and co-workers^52,64 in representing the SCF potential energy surface, Brocks et al.⁵³ constructed analytical expressions [μ^(x′)(R, ϑ), μ^(z′)(R, ϑ)] for the dipole moment of LiCN in the body-fixed reference frame (see Fig. 1). We use this dipole surface with the corrections involving sign errors, noted recently by Murgida et al.⁶² The z′-part of this dipole surface, which is at least 15 times larger than the dipole moment μ^(x′) in the x′-direction, is shown in Fig. 2(b). The body-fixed and space-fixed coordinate systems differ by a rotation with angle α, and thus, the dipole moment in the space-fixed coordinate system reads

Neglecting the small component μ^(x′)(R, ϑ) in Eq. (5), the dipole potential now reads

We then further approximate cos α ≈ 1 removing the corrections from the oscillations in α around the minimum at α = 0 and allowing us to reduce the dimensionality to only the two remaining degrees of freedom, R and ϑ. Although we have not fully explored the most general conditions for which this approximation will be valid, at the very least, they will be satisfied when the oscillations in α are faster than the other motion. In this limit, the potential on (R, ϑ) results from the effective field obtained from the average over α and reduces to a form with a renormalized prefactor, E₀, and no α dependence. We can thus focus on the reduced-dimensional Hamiltonian

where $H$ is the non-driven Hamiltonian of Eq. (2) and the time-dependent driving is included via the dipole potential

The equations of motion take the same form as for the non-driven case given in Eq. (3) with V replaced by the (time-dependent) potential

Again, solutions are found numerically using a fourth order Runge–Kutta algorithm⁶⁶ and a C++ implementation of the derivatives of V_driven with respect to R and ϑ.

The relative movement of individual atoms in the non-driven LiCN ⇌ LiNC isomerization reaction is described by trajectories obtained via the propagation of an initial state according to Eq. (3). The evolving state is identifiable as reactant $R$ if ϑ < 0.15π and product $P$ if ϑ > 0.5π when it is found in the corresponding regions highlighted in Fig. 3(a). The contours of the potential are also shown, providing a view of the underlying channel for trajectories to go between $R$ and $P$⁠.

In Fig. 3(a), representative trajectories are initialized at (R = 4.3, p_R = 0) for four different combinations of ϑ and p_ϑ as also highlighted by the corresponding marker in Fig. 3(b). Each full trajectory, displayed with a different line style, is obtained by propagating the initial point in phase space according to Eq. (3) forward and backward in time until the reactant or the product state is reached. Two of these trajectories are reactive (thick green solid with initial circle and thin white dashed one with initial square) and the other two (thin yellow solid with initial triangle and thin purple dashed one with initial diamond) are not since their energy is too low to cross the barrier.

The shaded regions of Fig. 3(b) are labeled the reactive (⁠$R→P$ and $P→R$⁠) and non-reactive (⁠$R→R$ and $P→P$⁠) initial points in the (ϑ, p_ϑ) subspace. These regions are separated by the stable and the unstable manifolds, which intersect at a particular point $(𝜗,p𝜗)NHIM$ of the NHIM for (R = 4.3, p_R = 0). The white dashed square trajectory is initialized closest to one of the manifolds and stays in the barrier region for nearly three oscillations in the stable direction of the barrier. Hence, it crosses a corresponding dividing surface (DS) separating reactants from products closest to the normally hyperbolic invariant manifold (NHIM). While the structures illustrated in Fig. 3 were previously seen in a generic model system with a driven rank-1 saddle,^24,25,39,44 the results here show that they are realized also in a specific model of a chemical reaction. In both cases, the geometry is described by the stable and unstable manifolds of the barrier, and the nature of the reactivity is determined by the associated DS attached to the NHIM.

When periodically driving the LiCN ⇌ LiNC isomerization reaction, the potential V_dip(R, ϑ, t) in Eq. (9) is explicitly time-dependent and the energy of the system is no longer conserved. Still, a two-dimensional NHIM exists in the barrier region, which contains all trajectories that never leave this region, neither forward nor backward in time. In contrast to a static system, however, the NHIM becomes time-dependent. For periodic driving, such a NHIM oscillates with the same frequency as that of the driving of the barrier because of the symmetry of the system with respect to time.

The Poincaré surface of section (PSOS) is a useful tool for resolving the dynamics of this time-dependent NHIM in a periodically driven system. Therein, the positions (R, p_R) of several trajectories on the NHIM are marked after integer multiples of the driving period for a time much longer than a single period of the external driving. When propagating trajectories for such a relatively long time, the instability of the dynamics in the NHIM can become problematic as errors produced by any numerical propagator are exponentially increasing. Consequently, trajectories initialized as precisely as possible in the NHIM nevertheless “fall down” from the barrier to either the reactant or the product side. This problem can be addressed through the use of a stabilized propagator,^41,44 which successively projects unstable trajectories back into the NHIM after an appropriately chosen time step using the binary contraction method (BCM).

In the static case according to Eq. (2), the energy is conserved, so each trajectory in Fig. 4(a) is periodic. The central fixed point corresponds to the trajectory resting at the saddle point of the barrier—viz., at R = 4.2626 and ϑ = 0.2800π. The frequency of oscillations in the stable direction of the static barrier decreases monotonously from about ω_orth = 0.043π for the trajectories close to the fixed point to about ω_orth = 0.035π for a trajectory with energy E = −0.19. The driving frequencies, chosen below, range between 2 and 4 times slower than the oscillation frequencies of the orthogonal modes and thus provide a non-trivial perturbation to the system. The PSOSs at these driving frequencies indeed show significant changes compared to the static case, as shown in Fig. 4.

For the driven systems, the stabilized trajectories are initiated on the NHIM at t₀ = 0. They are propagated for 400 periods of the external driving with amplitude E₀ = 0.01 and frequency ω = 0.01π in atomic units. Specifically, the values are chosen to correspond to a laser field with wave length λ = 1.43 µm and amplitude E₀ = 5.14 × 10⁷ V/cm. After each period of T = 200, the instantaneous position (R, p_R) of each trajectory in the NHIM is marked. The stabilization of each trajectory into the NHIM was performed using the BCM with an error tolerance of 10⁻¹². Finally, we observe that the PSOSs for the driving cases look quite different from the static case because the driven NHIM is time-dependent.

The driven trajectories are generally no longer periodic and energy is not conserved, leading to a more complex geometric structure. Nevertheless, at first glance, the structure of the trajectories of the driven barrier in Fig. 4(b) looks very similar to those of the non-driven barrier in Fig. 4(a), having a central fixed point and tori (trajectories in the non-driven case) around it. On the right-hand side of Fig. 4(b) at about R = 4.8 and p_R = −0.2, an unstable fixed point is visible, interrupting the regular arrangement of the tori. According to the Poincaré–Birkhoff theorem,⁶⁸ only an even number of fixed points may occur in a perturbed system (if regarded not precisely at the bifurcation values), and consequently, the appearance of the unstable fixed point is accompanied by the emergence of a stable fixed point. The latter is located approximately at R = 3.9 and p_R = 0 and is not directly visible here.

To magnify the influence of the moving barrier, the PSOS in Fig. 4(c) is obtained in the same way as in Fig. 4(b) but with twice the frequency of the external driving ω = 0.02π. Physically, this corresponds to laser light with wave length λ = 716 nm. Now, the external driving has a strong impact on the dynamics of the trajectories on the NHIM. An additional unstable and an additional stable fixed point are clearly visible. Both stable fixed points belong to separate period-1 trajectories that are encircled by many quasi-periodic trajectories on various tori. These structures, not seen in the non-driven barrier and only partially seen with weak driving, arise solely due to the external driving.

A three-dimensional representation (R, ϑ, p_R) of the three periodic trajectories corresponding to the stable fixed points clearly shown in Figs. 4(b) and 4(c) is given in Fig. 5(a). Here, the inner trajectory corresponding to the visible stable fixed point in Fig. 4(b) at initial position R = 4.2579 and p_R = −0.0117 shows just very little movement in the stable direction of the barrier. However, the trajectory in between, corresponding to the stable fixed point at R = 4.2040 and p_R = 0.0313, and especially the outer trajectory, corresponding to the stable fixed point at R = 5.0096 and p_R = −0.0378 of Fig. 4(c), are subject to significant movement in the stable direction of the barrier. Note that these two trajectories show two oscillations in the direction of the orthogonal mode of the barrier although they are both period-1 trajectories with respect to the external driving.

Using the methods introduced in Ref. 44 and discussed in the supplementary material, the instantaneous reactant decay rates associated with trajectories on the driven NHIM can be obtained. First, we analyze the three period-1 trajectories displayed in Fig. 5(a). The inner trajectory—shown in blue in panel (a)—is obtained for a driving frequency of ω = 0.01π, while the two outer trajectories—shown in blue and green in panel (a)—are obtained for a system with a larger driving frequency of ω = 0.02π, as labeled.

Figure 5(b) presents the instantaneous decay rates associated with trajectories close to the inner period-1 trajectory with an initial point R = 4.2579 and p_R = −0.0117 at time t = 0 on the NHIM of a driven barrier with E₀ = 0.01 and ω = 0.01π. To obtain decay rates with the ensemble method, the trajectory is divided into 20 segments. For each segment, an ensemble of 200 reactive trajectories is initialized close to the NHIM with a distance of Δϑ = 10⁻³. Patching together these individual segments yields the thin black line in Fig. 5(b). The corresponding mean ensemble decay rate is obtained by averaging over a full period of the external driving, yielding $k¯=0.0724$ and displayed as a vertical thick solid line. This result can be verified using the Floquet method (see the supplementary material), which yields a mean Floquet rate of $k¯F=0.0724$ for this trajectory. This mean Floquet rate, shown as a vertical orange dotted line in Fig. 5(b), is in perfect agreement with the mean ensemble rate. Using the local manifold analysis (LMA) (see the supplementary material), a third verification of these decay rates can be obtained. The thin red dashed line in Fig. 5(b) marks the instantaneous manifold rate k_M(t) and the vertical thick red dashed line marks the mean manifold rate $k¯M=0.0724$⁠, averaged over a full period of the external driving. Both lines fit the results of the other two methods.

The same procedure used to obtain Fig. 5(b) is repeated for the two remaining trajectories of Fig. 5(a). Figures 5(c) and 5(d) present the results for the intermediate trajectory (labeled c, red) and the outer trajectory (labeled d, green) in Fig. 5(a), respectively. In both cases, the instantaneous decay rates obtained via the LMA correspond perfectly to the instantaneous decay rates of the ensemble method, and their mean rates are also in perfect agreement with the obtained Floquet rates. Hence, $k¯=k¯M=k¯F=0.0730$ for the intermediate trajectory, and $k¯=k¯M=k¯F=0.0778$ for the outer trajectory.

The decay of the reactant population close to the three different period-1 trajectories in Figs. 5(b)–5(d) is represented by very different instantaneous rates and consequently also mean rates. In addition, the amount of variation of the instantaneous rates varies strongly for each trajectory. For the inner trajectory that, according to Fig. 5(a), has the smallest movement in the orthogonal mode, the relative change in the instantaneous rate, Δk ≈ 10⁻³, is very small. Consequently, the effect of the external driving is barely noticeable. The intermediate trajectory of Fig. 5(a) has significantly more motion in the stable direction of the barrier. The relative change, Δk ≈ 10⁻², in the instantaneous rate is considerably larger, as seen in Fig. 5(c). The outer trajectory of Fig. 5(a) has by far the largest movement in the stable direction of the barrier. The relative change, Δk ≈ 0.5 × 10⁻¹, in the instantaneous rates in Fig. 5(d) is also the largest. Thus, we can conclude that the influence of the external driving is high if a trajectory has significant motion in the direction of the orthogonal modes. Further evidence for this effect is given in Sec. III C in the context of the decay rates associated with the numerous quasi-periodic trajectories on the driven NHIM.

We now obtain the phase-space resolved decay rates of the reactant population close to arbitrary trajectories on the NHIM of the periodically driven LiCN ⇌ LiNC isomerization reaction. As seen above, all three methods to calculate the decay of reactant population close to the transition state (TS) result in the same values for a specific period-1 trajectory, and hence, we choose only one of these for the present calculation. Specifically, the Floquet method is employed because it is relatively easy to implement and computationally fast to evaluate.

In all the cases of Fig. 4, the Floquet rates are overlayed on top of the PSOSs. They are obtained on equidistant grids using ∼10 000 stabilized trajectories, and they are displayed through a shaded (or colored in color) encoding. All trajectories are propagated for a total time up to t = 10 000 as needed to converge. The individual Floquet rates are interpolated using bicubic splines to smooth the discrete points. We found that the decay rates of all the trajectories located on the corresponding regular tori are indeed the same. This is expected since any quasi-periodic trajectory on such a torus, in general, covers it in full if propagated for long enough. Thus, a unique decay rate is associated with each torus. This observation is true for the static system using a simpler argument. Specifically, the Floquet rate reduces to a property of any of the periodic trajectories on the static NHIM because the tori is necessarily periodic.

When comparing the obtained mean Floquet rates of the driven system with E₀ = 0.01 and ω = 0.01π in Fig. 4(b) to the mean Floquet rates of the static system according to Fig. 4(a), the influence of the external driving is small. In Fig. 4(b), the mean Floquet rate obtained at the clearly visible elliptic fixed point at the center of the tori corresponds to the rates already obtained in Fig. 5(b) with all three methods introduced in Ref. 44 and provided as the supplementary material. In the two cases of Figs. 4(a) and 4(b), the mean Floquet rates at the central fixed point are similar—that is, they are $k¯F=0.0721$ and $k¯F=0.0724$⁠, respectively. The primary difference is the emergence of a structure that is barely visible here but more clearly visible in Fig. 4(c), as discussed below.

A further increase in the frequency of the external driving to ω = 0.02π leads to a drastic change in the structure, as shown in Fig. 4(c). Two clearly separated elliptical fixed points are now visible in the PSOS with each encircled by an individual set of tori dividing the NHIM into two regions of very different Floquet rates. This means that two regions of very different stability emerge on the periodically driven NHIM. In the surrounding of the fixed point at R = 4.2040 and p_R = 0.0313, the decay rates of the reactant population are rather small and approximately correspond to the mean Floquet rate of $k¯F=0.0730$ obtained for the central period-1 trajectory according to Fig. 5(c). On the other hand, the mean Floquet rates in a region near the fixed point at R = 5.0096 and p_R = −0.0378 are rather large and approximately correspond to the mean Floquet rate $k¯F=0.0778$ of the associated central period-1 trajectory according to Fig. 5(d).

For a periodically driven system, the instantaneous energy E(t) of trajectories in the NHIM is not conserved. However, by averaging the instantaneous and non-conserved energy of a given trajectory over many oscillations of the periodic driving, the mean energy $E¯$ of a specific torus is well defined and characteristic of the trajectory and its associated initial condition. The quantitative correspondence between the mean energy of a trajectory and the Floquet rate is revealed in Fig. 6 by plotting the Floquet rates obtained in Fig. 4 with respect to the corresponding (mean) energies of the trajectories.

As a control for the driven cases, Fig. 6(a) shows the relation between the Floquet rate and the mean energy for the non-driven system of Fig. 4(a). In this limit, energy is conserved. Thus, the mean energy $E¯$ corresponds to the instantaneous energy E(t) of the respective trajectories. According to Fig. 6(a), the relation between the energy and the Floquet rate in the static case is not monotonic and has a maximum at approximately E = −0.18. A further increase in the energy of trajectories on the NHIM goes together with a decrease in the corresponding decay rate of the reactant population close to the TS. However, a maximum in these decay rates is associated with a minimum in the stability of the activated complex close to the TS. According to Fig. 6(a), there seems to be a lower bound for the stability of this activated complex in the non-driven LiCN isomerization reaction.

The results for the driven system in Fig. 6(b) for ω = 0.01π and in Fig. 6(c) for ω = 0.02π correspond to cases in Figs. 4(b) and 4(c), respectively. Despite the driving, they retain some similarity with the static case of Fig. 6(a). Again, the relation between the Floquet rate of a trajectory on a given torus and the associated mean energy is non-monotonic, and there seems to exist a maximum decay rate for the reactant population associated with trajectories on the NHIM.

In contrast to the static case, however, for the driven barrier case with ω = 0.02π, a significant gap arises in the curve of Fig. 6(c) at about $E¯=−0.22$⁠. This gap corresponds exactly to the emergence of tori around distinct fixed points seen in Fig. 4(c). At the beginning of the upper edge of the gap in Fig. 6(c), the curve initially fluctuates before becoming smooth again with increasing mean energy. The origin of this behavior is numerical because the Floquet method for non-periodic trajectories is defined only in the limit t → ∞ (see the supplementary material). However, trajectories in Fig. 4(c) are necessarily propagated for a finite time, and this leads to errors in the calculated Floquet rates, which are small and appear as fluctuations. In other words, close to the boundary between the two regions of the different decay rates on the driven NHIM, the circulation times of quasi-periodic trajectories on the regular tori increase and longer propagation times are necessary to obtain a desired accuracy using the Floquet method. This behavior can also be seen in Fig. 6(b) for the case with a driving frequency of ω = 0.01π. A second stable fixed point is barely visible in Fig. 4(b). Consequently, nearly no gap appears in Fig. 6(b). Nevertheless, fluctuations are visible at the energy of the indicated fixed point on the left-hand side of Fig. 4(b), and hence, a gap is just emerging there.

In this paper, we study the influence of a periodically oscillating external field on the LiCN ⇌ LiNC isomerization reaction. In the static case, the energy is conserved and all trajectories on the associated normally hyperbolic invariant manifold (NHIM) are periodic. For a periodically driven system, nearly all trajectories contained in the NHIM are quasi-periodic and located on stable tori. The dynamics on the driven NHIM has been revealed by means of a Poincaré surface of section (PSOS). Depending on the frequency of the external field, a new set of tori around a second stable fixed point emerges on the NHIM.

The associated mean decay rates of reactant population in a close neighborhood of the periodic or quasi-periodic trajectories on the NHIM are obtained by various methods.⁴⁴ Both in the static and driven cases, the decay rates differ significantly across these neighborhoods and are approximately prescribed by the mean decay rates associated with the corresponding central periodic trajectory. The instantaneous decay rates for these central trajectories are also sensitive to the driving. Specifically, the relation between the mean decay rates associated with a specific torus and the corresponding mean energy of trajectories on this torus is non-trivial. There are two clearly distinguished regions associated with different mean decay rates of the reactant population close to the transition state (TS). They are approximately prescribed by period-1 trajectories identified at the center of the corresponding tori of the PSOS. They are consequently not continuously connected and a significant gap between the corresponding decay rate emerges. The regularity of the structure of the PSOS for the trajectories near the NHIM is not chaotic as one might have expected, given the recently observed chaotic sea in the global phase space.⁶⁹ This points to the simplification arising from considering decay rates in the local neighborhood of the TS. Future work to extend these decay rates to the global rates thus needs to consider not just the possibility of additional barrier regions but also the complex dynamics that can arise from chaotic regions.

Thus, the influence of the periodic external driving on the LiCN ⇌ LiNC isomerization reaction is large when compared to the static problem. Indeed, the emergence of different regions of the reactant population decay for the driven LiCN ⇌ LiNC isomerization reaction and the existence of gaps in both energy and decay rate between these separate regions are effects solely invoked by the external driving.

The possible effects from our approximate simplification of the LiCN ⇌ LiNC isomerization reaction to a body-fixed axes may need to be addressed. This would allow the molecules to rotate with respect to the external field and increase the dimensionality of the coordinate space of the problem beyond the two-dimensional case considered here. However, there is precedent from the earlier work of Ref. 62 that the reaction is much faster than the rotation, and hence, such effects are small. One could also include a Langevin-type description of the dynamics—according to Ref. 57—to address the effects on the decay rates from a solvent as represented by the inclusion of both noise and friction. Alternatively, the effects on the decay rates from an all-atom solvent can be uncovered by inferring the stability of the NHIM from molecular dynamics trajectories. For example, the position of reactive trajectories across the dividing surface (DS) can be tracked by simulations of LiCN in an argon bath such as those of Refs. 54–57. Applying an appropriate external driving could shift these reaction positions to different regions of the NHIM and increases the reaction rate. From the results of the current work, we would expect that there would once again emerge different regions in the PSOS and the corresponding reactant decay rates, which should be interpretable as reactive channels. When applying an appropriate external field, a second and faster reactive channel with higher decay rates opens up and might serve to an increase in the reaction rate.

See the supplementary material for a detailed description of the methods used in the primary text, summary of the details for the implementation of three approaches in determining the decay rates—Floquet analysis, ensemble method, and the local manifold analysis, and a technical derivation of the last of these methods.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The German portion of this collaborative work was supported by the Deutsche Forschungsgemeinschaft (DFG) through Grant No. MA1639/14-1. R.H.’s contribution to this work was supported by the National Science Foundation (NSF) through Grant No. CHE-1700749. M.F. is grateful for support from the Landesgraduiertenförderung of the Land Baden–Württemberg. This collaboration has also benefited from the support by the European Union’s Horizon 2020 Research and Innovation Program under the Marie Sklodowska–Curie Grant Agreement No. 734557.