Obtaining time-dependent multi-dimensional dividing surfaces using Lagrangian descriptors

TL;DR

This paper extends a Lagrangian descriptor-based method to construct time-dependent, multi-dimensional dividing surfaces that eliminate recrossings, enabling exact transition state theory rates in complex, time-dependent systems.

Contribution

It introduces a formalism for constructing recrossing-free, time-dependent dividing surfaces in systems with multiple degrees of freedom using Lagrangian descriptors.

Findings

Dividing surfaces are recrossing-free in numerical simulations.

The method yields exact transition state theory rates.

Extension to systems with bath degrees of freedom is demonstrated.

Abstract

Dynamics between reactants and products are often mediated by a rate-determining barrier and an associated dividing surface leading to the transition state theory rate. This framework is challenged when the barrier is time-dependent because its motion can give rise to recrossings across the fixed dividing surface. A non-recrossing time-dependent dividing surface can nevertheless be attached to the TS trajectory resulting in recrossing-free dynamics. We extend the formalism---constructed using Lagrangian Descriptors---to systems with additional bath degrees of freedom. The propagation of reactant ensembles provides a numerical demonstration that our dividing surface is recrossing-free and leads to exact TST rates.

Figures

Equations

$$\begin{array}{rl}V(x,y,t)=& E_{\mathrm{b}}\exp \left(-a{\left[x-\hat{x}\sin \left({\omega }_{x}t\right)\right]}^2\right)+\frac{{\omega }_{\mathrm{y}}^2}{2}{\left[y-\frac{2}{\pi }\arctan \left(2x\right)\right]}^2.\end{array}$$

$$\mathcal{L}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)={\int }_{t_0-\tau }^{t_0+\tau }||\mathbf{v}(t)||\mathrm{d}t.$$

$${\mathcal{W}}_{\mathrm{s}}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)$$

$${\mathcal{W}}_{\mathrm{u}}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)$$

$$\mathcal{T}(x_0,v_0,t)\equiv {\mathcal{W}}_{\mathrm{s}}(x_0,v_0,t)\cap {\mathcal{W}}_{\mathrm{u}}(x_0,v_0,t),$$

$$\rho (\mathbf{x},\mathbf{v})={\rho }_{\mathrm{therm}}\delta (x+2)\Theta (v_x),$$

$$p_{\mathrm{r}}(t)=p_{\mathrm{r},0}{\mathrm{e}}^{-kt}+c$$

$$k(k_{\mathrm{B}}T)=k_{\infty }\exp \left(-\frac{\Delta E_{\mathrm{eff}}}{k_{\mathrm{B}}T}\right),$$

Full text

1 Introduction

The accuracy in the determination of reaction rates relies on the precision with which reactants and products can be distinguished in the underlying state space. Usually, the boundary between these regions contains an energetic saddle point in phase space to which an appropriate dividing surface (DS) can be attached. Transition state theory (TST) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18] then provides a powerful basis for the qualitative and quantitative description of the reaction. The rate is obtained from the flux through the DS and it is exact if and only if the DS is free of recrossings. Advances in the determination of this fundamental quantity can impact a broad range of problems in atomic physics [19], solid state physics [20], cluster formation [21, 22], diffusion dynamics [23, 24], cosmology [25], celestial mechanics [26, 27], and Bose-Einstein condensates [28, 29, 30, 31, 32], to name a few.

In autonomous systems, the recrossing-free DS is attached to a normally hyperbolic invariant manifold that can be constructed using e. g. normal form expansions [33, 34, 35, 36, 37, 26, 38, 39, 40, 41]. The situation becomes fundamentally different if the system is time-dependent, e. g. if it is driven by an external field or subject to thermal noise. In one-dimensional time-dependent systems, a DS with the desired property is given by the transition state (TS) trajectory [42, 43, 44, 45, 46, 47, 48, 49, 50, 51] which is a unique trajectory bound to the vicinity of the saddle for all time.

In systems with dimension greater than one, the reacting particle can simply bypass the TS trajectory (point) by having a non-zero velocity perpendicular to the reaction coordinate. Thus one must attach a multi-dimensional surface to the TS trajectory that separates reactants and products. The use of perturbation theory in multi-dimensional cases provides both the TS trajectory and the associated geometry on which this dividing surface can be constructed. The challenge, addressed in this Letter, is how to obtain this multi-dimensional structure without perturbation theory. One possible approach lies in the use of the Lagrangian descriptor (LD) [52, 53] used recently by Hernandez and Craven [54, 55] to obtain the TS trajectory without resolving the DS at higher dimension. This alternate framework is necessary when there is no useful reference such as in barrierless reactions [49], and more generally to avoid the convergence issues that invariably plague a perturbation expansion far from the reference. In the case of field-induced ketene isomerization [55], the LD was computed across the entire phase space. It not only revealed the structure of the DS, but also coincided with the final state basins for each initial condition in phase space for both 1-dimensional and 2-dimensional representations. However, while the approach is formally applicable to arbitrary dimension, we have found that it is difficult to perform the minimization of the naive LD, even in dimensions as low as two.

The time-dependent Lagrangian descriptor dividing surface (LDDS), introduced in this Letter is the natural extension to $n$ dimensions for $n>1$. We freely choose $2n-2$ phase-space coordinates for which we fix the initial condidions, and use the LD approach to identify a corresponding trajectory, which we call an anchor trajectory. It is defined by the intersection of the stable and unstable manifolds of the time-dependent Hamiltonian [56, 57]. The TS trajectory is the anchor trajectory—which necessarily remains in the vicinity of the TS region for all past and future time—with the least vibrational motion orthogonal to the reactive degree of freedom. The LDDS is attached to the family of anchor trajectories and is necessarly $(2n-1)$-dimensional. In the special case of a one-dimensional system ($n=1$), the LDDS coincides with the moving DS on the TS trajectory [54].

2 Theory and Methods

2.1 Two-Dimensional Model System

We illustrate the construction of the LDDS by modeling the dynamics of a two-dimensional chemical reaction with stationary open reactant and product basins. Hamilton’s equation of motion propagates the particle according to a non-autonomous Hamiltonian in mass-weighted coordinates with potential

[TABLE]

Here, $E_{\mathrm{b}}$ is the height of a Gaussian barrier with width $a$ oscillating along the $x$ axis with frequency $\omega_{x}$ and amplitude $\hat{x}$, $\omega_{y}$ is the frequency of the harmonic potential in the $y$ direction, and the term $({2}/{\pi})\arctan\left(2x\right)$ is the minimum energy path whose form induces a nonlinear coupling between the two degrees of freedom. For simplicity, all variables are presented in dimensionless units, where the scales in energy (and $k_{\mathrm{B}}T$), length, and time are set according to half the maximum barrier height of the potential, twice the variance of the Gaussian distribution, and the inverse of the periodic frequency, respectively. In these units, the dimensionless parameters in Eq. (1) are set to $E_{\mathrm{b}}=2$, $a=1$, $\omega_{x}=\pi$, $\omega_{\mathrm{y}}=2$, and $\hat{x}=0.4$.

2.2 Using Lagrangian Descriptors to Obtain Dividing Surfaces

As the dividing surface between reactant and product basins is in general a high-dimensional hypersurface, the stable and unstable manifold itself become high-dimensional objects. In the context of TST, the Lagrangian descriptor (LD) at position $\boldsymbol{x}_{0}$, velocity $\boldsymbol{v}_{0}$, and time $t_{0}$, is defined as the integral [54, 49, 51],

[TABLE]

It is a measure of the arc length of the unique trajectory $\boldsymbol{x}(t)$ in forward and backward time over the time interval $[t_{0}-\tau;~{}t_{0}+\tau]$, and the parameter $\tau$ is chosen such that it covers the relevant time scale of the system [in this letter, we use $\tau=10$ corresponding to five periods of the oscillating barrier in Eq. (1)]. The importance of the LD (2) naturally results from the fact that the stable and unstable manifolds $\mathcal{W}_{\mathrm{s,u}}$ which are attached to the barrier top in phase space, correspond to the minimum of the forward (f: $t_{0}\leq t\leq t_{0}+\tau$) and backward (b: $t_{0}-\tau\leq t\leq t_{0}$) contributions to the LD,

[TABLE]

Here, the function $\operatorname*{arg\,min}$ denotes the argument of the local minimum of the LD hypersurface close to the barrier top. In $n$ dimensions, we fix $2n-2$ variables freely (which would be least associated with the reactive degree of freedom) and perform the minimization in Eqs. (3). The intersection of these two manifolds

[TABLE]

is the $t=0$ value of the anchor trajectory to which a moving DS can be attached. The central result of this Letter is that the the family of these anchor trajectories $\mathcal{T}(t)$ carries the associated family of moving dividing surfaces that we call the Lagrangian descriptor dividing surface (LDDS), and that we show below to be a recrossing-free DS. The anchor surface $\mathcal{T}(t)$ is a $(2n-2)$-dimensional object embedded in the $2n$-dimensional phase space meaning in the special case of a one-dimensional system, this intersection is a single point, namely the position of the TS trajectory at given time [54].

The algorithm used to obtain $\mathcal{T}(t)$ can be explained by means of one of the insets in Fig. 1. These insets show the LD of an $x$-$v_{x}$-section for a certain time $t$ and fixed $y$ and $v_{y}$. The LD is calculated according to Eq. (2) by integrating trajectories with the respective initial conditions ($x,y,v_{x},v_{y},t$). They are obtained using a standard (symplectic) Velocity-Verlet integrator with a sufficiently small time-step to capture the time-dependence in the potential, Eq. (1), and to ensure convergence in the final positions and velocities. The structure of the stable and unstable manifold is identified through the local minima in the LD’s $x$-$v_{x}$-section. Their intersection yields the phase space coordinates, $x$ and $v_{x}$, of the point $\mathcal{T}(y,v_{y},t)$ to which the LDDS is attached. Repeating this procedure for an equidistant grid in the $y$-$v_{y}$ space (for a fixed time $t$) results in a mesh of points of the LDDS $\mathcal{T}(t)$. The smooth surfaces shown here are constructed through spline interpolation of this mesh.

3 Results

3.1 Trajectory Analysis

In Fig. 1 (center), we present a typical reactive trajectory (red solid line) undergoing a transition from the reactants ($x\to\infty$) to products ($x\to-\infty$). Because of the oscillating barrier in the two degree of freedom system (1), the trajectory shows several loops close to the barrier top. Its dominant motion is perpendicular to the reaction coordinate, but the trajectory also shows oscillations along the latter. Such nontrivial oscillations are a general feature of particles with an energy slightly above the barrier top. As a consequence, it is not generally possible to define a recrossing-free DS in the configuration space alone.

Although the particle’s dynamics is rather complicated near the barrier top, the reaction dynamics becomes clearer by focusing on the relative motion of the particle with respect to the time-dependent manifolds. In Fig. 1, phase space portraits of the LD are displayed for eight illustrative points along the selected trajectory. The stable (unstable) manifold corresponding to the minimum valleys of the LD according to Eq. (3) is shown as a black (yellow) dashed line. The time-dependent position $x^{\times}(y,v_{y})$ where they intersect is highlighted by a vertical, black dotted line. In the first three points, the particle is on the RHS of $x^{\times}(y,v_{y})$, crosses it at point 4, and then remains on the LHS of $x^{\times}(y,v_{y})$ for the last 4 point, as noted with the corresponding symbol defined in the caption. Each of the LD plots in the insets—labeled according to the corresponding point $1,\ldots,8$—shows an $x$-$v_{x}$-cut through phase space for the instantaneous values $y,v_{y}$ at the respective times $t$. In this and every other trajectory we have sampled, the particle crosses the corresponding $x^{\times}(y,v_{y})$ not more than once satisfying the recrossing-free criteria. For a single trajectory (that fixes $y$ and $v_{y}$ as the two remaining degrees of freedom in phase space and therefore leads to an effective one-dimensional system), the intersection of the manifolds (4) thus defines a recrossing-free DS that coincides with the TS trajectory of the effective one-dimensional system.

In the full phase space description of the two-dimensional system defined in Eq. (1), we can define a family of intersections $x^{\times}(y,v_{y})$ whose values and time-dependence vary depending on the two remaining bath coordinates (here $y$, $v_{y}$) according to Eq. (4). Indeed the union of these intersections is the time-dependent, two-dimensional anchor surface to which we attach the LDDS. The $\mathcal{T}(t)$ of Eq. (4) in the $x$, $y$, $v_{y}$ subspace is displayed in Fig. 2 (panel 7) at the time of the corresponding point in the trajectory. It is located near the saddle of the potential (1), shown as a contour-surface below, and exhibits a nontrivial curvature along all the axes. Note that the calculation of the associated time-dependent LDDS via the intersections of the manifolds in the $(x,v_{x})$ frame for a given $y$ and $v_{y}$ is exemplary, and the analogous approach using a different frame, such as $(y,v_{y})$, leads to the same result.

A higher-dimensional representation of snapshots of the anchor surface for the same trajectory of Fig. 1 is shown in Fig. 2. The nontrivial motion of the curved anchor hypersurface over time emerges as one follows it in relation to the fixed black dashed line. The oscillation of the anchor surface is in phase with that of the barrier top, but its amplitude is about one order of magnitude smaller in configuration space. The recrossings of the LDDS attached to the anchor surface due to the loops in the particle’s trajectory of Fig. 1 are avoided by the motion and appropriate bending of the surface if the particle re-approaches the barrier region.

3.2 Ensemble Analysis

We have, so far, demonstrated the recrossing-free nature of the LDDS attached to an anchor surface for a single trajectory. In the following, we extend this verification to an ensemble of trajectories. Specifically, 160 000 particles are initialized on an equidistant grid along the $x$-, $y$-, $v_{x}$-, and $v_{y}$-axes with 20 points along each axis. The grid is located close to the barrier top at $t=0$ and the grid size is chosen such that significant numbers of, both, reacting as well as nonreacting particles are observed during the time-evolution. For this ensemble, we compare in Fig. 3 different DSs with respect to their number of (re-)crossings. As can be seen in the first panel, the time-dependent LDDS attached to the anchor surface $\mathcal{T}(t)$ provides a recrossing-free DS for the whole ensemble of particles. For comparison, we introduce additional DSs for which we again calculate the number of recrossings yielding the other panels in Fig. 3. A naive way to construct a DS would be a planar surface perpendicular to the minimum energy path and attached to the time-dependently moving barrier top $x^{\ddagger}(t)$ which we refer to as $\mathcal{P}(t)$. In a second case, called $\mathcal{P}(0)$, we keep this planar surface fixed at the saddle’s initial position $x^{\ddagger}(0)$. The last two histograms are calculated using the LDDS for the fixed, time-independent potential $V(x,~{}y,~{}0)$ of Eq. (1). The resulting time-independent DS is either fixed at the saddle’s initial position [$\mathcal{T}_{\mathrm{f}}~{}@~{}x^{\ddagger}(0)$] or periodically moving with the top of the barrier [$\mathcal{T}_{\mathrm{f}}~{}@~{}x^{\ddagger}(t)$]. As Fig. 3 shows, the LDDS is the only choice for which either no or only single crossings occur while two and more crossings do not occur. By contrast, all the other choices of DSs exhibit several recrossings and will therefore lead to overestimates in the corresponding rates.

Finally, we regard a thermal ensemble of $10^{6}$ particles in the reactant well with the density distribution

[TABLE]

where $\rho_{\mathrm{therm}}$ is a Boltzmann-distribution, $\delta$ is the Dirac-delta function and $\Theta$ the Heaviside step function. The corresponding time-evolution of the reactant population $p_{\mathrm{r}}(t)$ is shown in Fig. 4. The decay corresponding to the LDDS is the only strictly monotonic one. Both of the time-independent DSs exhibit a modulation of the reactant population with a very small but nonvanishing amplitude. When fixed DSs are attached to the barrier top, $x^{\ddagger}(t)$, the oscillations are huge. Each increase in the number of reactants is due to recrossings through the respective surface. Note that for $t\rightarrow\infty$, the reactant population for all DSs is the same, as the particles have fallen down from the barrier either in the reactant or in the product basin and their classification is independent of the choice of the dividing surface, as long as it is located sufficiently close to the saddle.

The red dashed line in Fig. 4 shows an exponential fit

[TABLE]

to the long-time decay of the reactant population with fit coefficients $p_{\mathrm{r,0}},k,c$ from which we extract the reaction rate $k$. The respective rates obtained for the different DSs at various temperatures $k_{\mathrm{B}}T$ are shown in the inset of Fig. 4. Rates obtained from the LDDS are the smallest throughout as should be expected from a recrossing-free DS. The time-independent $\mathcal{P}(0)$ and $\mathcal{T}_{\mathrm{f}}(0)$ yield slightly higher rates while the rates obtained from the DSs attached to the barrier top [$\mathcal{P}(t)$ and $\mathcal{T}_{\mathrm{f}}(t)$] overestimate the rates by a factor of 2 to 5 due to recrossings. In addition, the LDDS rates exhibit a temperature dependence according to Arrhenius’ rate equation (see red fit curve in the inset of Fig. 4),

[TABLE]

where $\Delta E_{\mathrm{eff}}=0.135$ is the effective height of the potential and $k_{\infty}=3.14$ is the high-temperature limit of the rate. The effective barrier height is significantly lower than the spread of barrier heights between 1.70 to 2.00 at the naive saddle point. The fact that $\Delta E_{\mathrm{eff}}$ is not higher than these barrier heights is a good consistency check. It is lower in energy because the driving of the system maintains the effective reactant population in an activated state energetically higher than the naive reactant population near the minimum of the potential.

4 Conclusion

In this Letter, we have developed and verified the explicit construction of the time-dependent LDDS for multidimensional systems as a recrossing-free dividing surface. It reduces to the well-known TS trajectory [42] formalism in the one-dimensional limit.

The central result of this Letter is thus the justification and validation of a generalization of the LDDS method whose underlying minimization converges even in higher dimensions. Specifically, the construction is realizable because the LD (2) remains a scalar quantity regardless of the phase space dimension. The construction also results in a global DS, i. e. the recrossing-free property does not only hold close to the barrier top but for the complete hypersurface in full phase space. The method thus positions us to resolve the dividing surface (and associated reaction rates) in time-dependent molecular reactions [58, 59, 60, 61, 62, 63], and perhaps also spintronic devices [64, 65] driven by tailored external fields or thermal noise.

Acknowledgments

AJ acknowledges the Alexander von Humboldt Foundation, Germany, for support through a Feodor Lynen Fellowship. RH’s contribution to this work was supported by the National Science Foundation (NSF) through Grant No. CHE-1700749. This collaboration has also benefited from support by the people mobility programs, and most recently by the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No. 734557. Surface plots have been made with the Mayavi software package [66].