Chemical dynamics between wells across a time-dependent barrier: Self-similarity in the Lagrangian descriptor and reactive basins

TL;DR

This paper investigates the complex behavior of chemical reaction dynamics across time-dependent barriers, revealing self-similarity in Lagrangian descriptors and proposing modifications to improve the identification of transition state trajectories.

Contribution

It introduces a novel analysis of self-similar structures in Lagrangian descriptors for time-dependent potentials and proposes a modified LD approach to better identify transition states.

Findings

Self-similar structure of LD in bounded potentials

Reflections and side-minima affect LD-based TS trajectory detection

Modified LD approach improves reaction rate calculations

Abstract

In chemical or physical reaction dynamics, it is essential to distinguish precisely between reactants and products for all time. This task is especially demanding in time-dependent or driven systems because therein the dividing surface (DS) between these states often exhibits a nontrivial time-dependence. The so-called transition state (TS) trajectory has been seen to define a DS which is free of recrossings in a large number of one-dimensional reactions across time-dependent barriers, and, thus, allows one to determine exact reaction rates. A fundamental challenge to applying this method is the construction of the TS trajectory itself. The minimization of Lagrangian descriptors (LDs) provides a general and powerful scheme to obtain that trajectory even when perturbation theory fails. Both approaches encounter possible breakdowns when the overall potential is bounded, admitting the possibility of returns to the barrier long after trajectories have reached the product or reactant wells. Such global dynamics cannot be captured by perturbation theory. Meanwhile, in the LD-DS approach, it leads to the emergence of additional local minima which make it difficult to extract the optimal branch associated with the desired TS trajectory. In this work, we illustrate this behavior for a time-dependent double-well potential revealing a self-similar structure of the LD, and we demonstrate how the reflections and side-minima can be addressed by an appropriate modification of the LD associated with the direct rate across the barrier.

Figures

Equations

$$V(x,t)=V_{\text{Morse}}^{-}(x)+V_{\text{Gauss}}(x,t)+V_{\text{Morse}}^{+}(x)$$

$$V_{\text{Morse}}^{\pm }(x)$$

$$V_{\text{Gauss}}(x,t)$$

$$E_{x_{\text{b}}(t)=0}^{‡}=D(3+2{\mathrm{e}}^{-2b{x}_0}-4{\mathrm{e}}^{-b{x}_0})\approx 2.806,$$

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

$${\mathcal{W}}_{\text{s}}(t_0)$$

$${\mathcal{W}}_{\text{u}}(t_0)$$

$$\mathcal{T}(t_0)=\mathrm{argmin}\mathcal{L}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0).$$

$$\Bigl{[}t_{0}-\tau,~{}t_{0}+\tau\Bigr{]}~{}\longrightarrow~{}\Bigl{[}t_{0}-\tau^{-}[\boldsymbol{x}(t)],~{}t_{0}+\tau^{+}[\boldsymbol{x}(t)]\Bigr{]}\,.$$

$$\tau^{\pm}[\boldsymbol{x}(t)]=\min\Bigl{(}\tau,~{}t\bigl{.}\bigr{|}_{|\boldsymbol{x}(t)|>\bar{x}}\Bigr{)}$$

$$\mathcal{L}\approx m{\mathcal{L}}_{\text{bo}}+n{\mathcal{L}}_{\text{wo}}$$

Full text

I Introduction

One of the grand challenges in the field of driven reaction dynamics is the complete characterization of the rates and pathways so as to allow for control. One possible route for such driving is the perturbation of the underlying potential energy surface by time-dependent, external fields. Yamanouchi (2002); Hershkovits and Hernandez (2005); Sussman et al. (2006); Kawai et al. (2007); Kawai and Komatsuzaki (2011); Sethi and Keshavamurthy (2009); Patra and Keshavamurthy (2015); Revuelta, Chacón, and Borondo (2015) The configurational change of the reactive system is typically mediated by an energy barrier separating reactant and product basins which must be somehow affected by the external control mechanism. In the limit of no driving, transition state theory (TST) Pitzer, Smith, and Eyring (1962); Pechukas (1981); Garrett and Truhlar (1979); Truhlar, Issacson, and Garrett (1985); Natanson et al. (1991); Truhlar, Garrett, and Klippenstein (1996); Truhlar and Garrett (2000); Komatsuzaki and Berry (2001); Waalkens, Schubert, and Wiggins (2008); Bartsch et al. (2008); Kawai and Komatsuzaki (2010); Hernandez, Bartsch, and Uzer (2010); Sharia and Henkelman (2016) provides a powerful, though usually approximate, framework to calculate the rate from the reactive flux though the dividing surface (DS) separating the reactant and product regions. Such rates are exact if the DS is crossed by each reactive trajectory exactly once. Thus a central task for applying TST is the determination of a DS with this no-recrossing property. In time-independent systems with a two-dimensional configuration space, the DS is associated with an unstable periodic orbit at the barrier top, and in higher-dimensional systems it can be constructed using a normally hyperbolic invariant manifold. Pollak and Pechukas (1978); Pechukas and Pollak (1979); Hernandez and Miller (1993); Hernandez (1994); Jaffé, Farrelly, and Uzer (2000); Koon et al. (2000); Jaffé et al. (2002); Uzer et al. (2002); Waalkens and Wiggins (2004); Jaffé et al. (2005); Li et al. (2006); Teramoto, Toda, and Komatsuzaki (2011); Çiftçi and Waalkens (2013) By contrast, in time-dependent systems, the DS is, in general, also time-dependent and the transition state (TS) trajectory, Bartsch, Uzer, and Hernandez (2005); Bartsch, Hernandez, and Uzer (2005); Bartsch et al. (2006); Craven, Bartsch, and Hernandez (2014a, b, 2015); Kawai and Komatsuzaki (2009) which is a hyperbolic trajectory close to the barrier top, has proven to give rise to an associated non-recrossing time-dependent DS.

The TS trajectory can be constructed through perturbation theory in several limiting cases,Hernandez and Miller (1993); Bartsch, Hernandez, and Uzer (2005) but the approach does not have an obvious zeroth-order reference in barrier reactionsJunginger and Hernandez (2016) and can break down when the dynamics is affected by features on the potential energy surface far from the barrier region. It has recently been shown that the minimization of Lagrangian descriptors Mendoza and Mancho (2010); Mancho et al. (2013) (LDs) provides a general and powerful construction scheme to obtain the TS trajectory in such difficult cases.Craven, Bartsch, and Hernandez (2015); Craven and Hernandez (2015); Junginger and Hernandez (2016); Craven and Hernandez (2016) In simple terms, the initial condition for the TS trajectory is the one for which the LD, integrated for some sufficiently long time, is a minimum over the domain of the underlying phase space coordinates.

The present work revisits the prototypical reaction in a time-dependent double-well potential with oscillating barrier position from our earlier work.Junginger et al. (2016) In this paper, we address a possible challenge to the minimization in the LD-DS procedure arising from the reflections of particles when the overall potential has characteristic bound reactant and product wells. Without such wells, particles leaving sufficiently far from the reaction region escape into a deep exit channel from which they never return. When the overall potential is bounded, however, the global motion gives rise to (uncorrelated) returns back to the reaction region despite the fact that the particles relaxed into a given basin in the intervening time. This effect leads, in general, to several or formally an infinite number of local LD minima, making it difficult to identify the optimal LD minima associated with the TS trajectory. We demonstrate in this paper that the LD surface in phase space exhibits a huge number of local minima which are related to trajectories leaving and reentering the barrier region. Equally important, we also show how to identify the primary minimum needed to construct the TS trajectory for calculating direct rates, and thereby resolve the use of the LD-DS method for typical (bounded) chemical reactions.

The extended LD-DS method needed to account for global recrossings is developed in Sec. II. Specifically, the LD is modified so that only a single minimum remains, and it is precisely the one corresponding to the TS trajectory. It corresponds to the imposition of reactant and product surfaces defining the entry into the corresponding regions in the reactive flux time correlation function for the rate, and which are imposed to avoid recurrences from the reverse reaction. Such an approach has been used routinely in the chemistry community, and is beautifully imposed in the recent work of Mesele and ThompsonMesele and Thompsona (2016) in obtaining expressions for the activation energy. This section also includes the structure and parameters of the model double-well used in this work to illustrate the LD-DS method. In Sec. III, we investigate the self-similar LD structure in detail with respect to the model system. The trajectories associated with the local minima are identified clearly using the modified LD, and the optimal one leads to a numerical construction of the TS trajectory. This central result is essential to the extension of the LD-DS method towards higher-dimensional, finite reactive systems. The latter necessarily includes global recrossings which would bedevil perturbation theory and the naive LD-DS approach.

II Theory

In this paper, we investigate the time-dependent double-well potential (see Fig. 1)

[TABLE]

that has already been subject of the work in Ref. Junginger et al., 2016. This potential consists of two Morse potentials $V_{\text{Morse}}^{\pm}$ providing the reactant ($-$) and product ($+$) wells and a time-dependent Gaussian potential $V_{\text{Gauss}}$ serving as the time-dependent barrier. Specifically, these dimensionless potential terms are

[TABLE]

where $x_{\text{b}}$ is the position of the barrier, the factor $D$ determines the energy scale, and $\pm x_{0}$ are the positions of the wells. As in Ref. Junginger et al., 2016, we use $D=1$ for the energy scale, and specify the spacial scale of the system by setting $x_{0}=3$ as well as $a=b=1$. Moreover, we apply a sinusoidal driving of the barrier according to $x_{\text{b}}(t)=0.4\sin(\pi t)$, so that the barrier oscillation is $0.4/3\approx 13\%$ of the distance to the well’s minimum. For a centered barrier top at $x_{\text{b}}(t)=0$, the barrier height is

[TABLE]

which varies only slightly during the oscillation of the barrier top. For simplicity, we use mass-weighted coordinates throughout this paper. Together with the length scale given by the position $x_{0}$ of the well minima and the barrier height $E^{\ddagger}$ as a reference for the energy scale, this defines the dimensionless units of the system.

As shown in Refs. Craven and Hernandez, 2015; Junginger and Hernandez, 2016; Craven and Hernandez, 2016, a nonperturbative approach to constructing a time-dependent, recrossing-free DS is given by a minimization procedure of LDs with respect to the underlying phase space coordinates. In the context of TST, the LD is defined by the integral

[TABLE]

where $\boldsymbol{v}$ is the velocity of a certain trajectory and therefore $\mathcal{L}$ is a measure of the trajectory’s arc length over the time interval $[t_{0}-\tau;t_{0}+\tau]$.

The LD is of special importance for the reaction dynamics, because the stable and unstable manifolds $\mathcal{W}_{\text{s,u}}$ attached to the (time-dependent) barrier top are directly related to the LD’s forward ($f$: $t_{0}\leq t\leq t_{0}+\tau$) and backward ($b$: $t_{0}-\tau\leq t\leq t_{0}$) contribution, respectively, according to Craven and Hernandez (2015); Junginger and Hernandez (2016); Craven and Hernandez (2016)

[TABLE]

This is the case because the trajectories on these manifolds approach the barrier top in either forward or backward time yielding extremal properties of the LD. In Eq. (5), the function ‘$\operatorname{arg\,min}$’ denotes the value of the LD argument, i. e. the respective phase space coordinates at the (local) minimum.

The coordinates of the TS trajectory $\mathcal{T}(t_{0})$ at time $t_{0}$ are located on the intersection of these manifolds, Craven and Hernandez (2015); Junginger and Hernandez (2016); Craven and Hernandez (2016) $\mathcal{T}(t_{0})=\mathcal{W}_{\text{s}}(t_{0})\cap\mathcal{W}_{\text{u}}(t_{0})$, so that it is directly related to the two-sided LD (4) via

[TABLE]

As we demonstrate below, the definition (4) of the LD with fixed integration time $\tau$ in a closed system with finite reactant and product wells has the disadvantage that the LD, in general, does not only exhibit a single minimum, but a huge number of minima. This makes the identification and construction of the TS trajectory in Eq. (6) difficult or even ambiguous. The reason lies in the fact that when a particle has an energy high enough to cross the barrier at least once, it will continuously be reflected to the barrier and undergo repeated barrier crossings (that is, global recrossings or recurrences) which are not associated with the direct rate across the barrier between reactant and product. Thus one needs to define the reaction region thorugh which trajectories traverse from the reactant to product regions in order to obtain the direct rae across the barrier. A unique distinction between reactants and products can be made locally at the saddle (in terms of local crossings) by means of the TS trajectory. To overcome the issue of global recurrences in the construction of the TS trajectory we modify the definition of the LD (4) by simply replacing the fixed integration time $\tau$ with a variable one depending on the underlying trajectory. For this purpose, we redefine the time interval over which the trajectory is integrated through the map:

[TABLE]

Here, the values $\tau^{\pm}[\boldsymbol{x}(t)]$ are trajectory-dependent integration times which we define according to

[TABLE]

with an appropriate “size” $\bar{x}$ of the TS region. The redefinition of the integration time in Eq. (8) is thereby limited by both $\tau$ and the time for the particle to leave the barrier region (when $|x|$ is greater than $\bar{x}$.) The redefinition in Eq. (7) may appear ad hoc, but is motivated by the fact that all the undesired LD properties clearly result from recurrencs in the trajectories being globally reflected to the barrier region. By construction, our redefinition precisely retains only the trajectories the contribute to the direct rate. The value of $\bar{x}$ in Eq. (8) needs to be appropriately defined: On the one hand, it should not affect trajectories which do not leave the barrier; this suggests that it be a minimal value on the order of the barrier oscillation amplitude. On the other hand, it should be large enough to remove the effect of global recrossings of particles returning from the reactant or product wells; this suggests that its maximum value should be on the order of the distance between the barrier top and the well’s minimum. Finally, we note that the integration times $\tau^{\pm}$ may also differ in the forward and backward time direction.

III Results

In this section, we apply the LD formalism, Eq. (6), to the double-well potential in Eq. (1). In Fig. 2, we present the LD (top row) as well as the reactive basin portraits (bottom row) in phase space for a fixed integration time $\tau=20$ and the panels (a)–(e) are different zoom levels close to the TS. The reactive basin plots are coded through colors denoting whether or not a particle is reactive: black and white regions indicate nonreactive particles which either start and end in the reactant (black) or product (white) basin, and the red regions show reactive particles which undergo forward (dark red) or backward (light red) reactions. The rich structure of the LD landscape, in the top left panel of Fig. 2 for example, indicates the presence of chaotic and regular regions in the dynamics.

The LD portraits show recurring structures at the different zoom levels (a)–(e) in Fig. 2, revealing a self-similar structure of the LD in a regular region of the phase space (which we have further verified to be present down to the limit of numerical accuracy; not shown). Equally important, the reactive basin plots show the same recurring structure and the borders between the different reactive regions coincide with those observed in the LD plots. This verifies previous results Craven and Hernandez (2016); Junginger et al. (2016) in which such a connection has been observed and extends them to the present case including reflections of the particles at the potential walls.

In Fig. 3(a), we illustrate the LD portrait from Fig. 2(a) together with a selection of its (local) minima: Each of the white dots represents a local minimum of the LD as the result of a numerical search. We note that the minima highlighted by the white dots are only a selection of points. As can be seen in the LD portrait, there are more intersections of the manifolds related to more local minima, but not all of them are highlighted in order to not overload the figure. The minima are located either on one of the manifolds $\mathcal{W}_{\text{s,u}}$ which are equivalent to LD minimum valleys according to Eq. (5) or on their intersections, respectively. However, even for this amount of local minima, it becomes obvious that Eq. (6) is difficult to apply in order to locate the single minimum related to the TS trajectory.

Taking into account not only the relation (6), but also looking at the actual LD values of the minima can help to single out the desired minimum: This can be understood from Fig. 3(b) where the LD value of the local minima is visualized for the different trajectories. (Note that, on the horizontal axis, the different trajectories are sorted ascending according to their LD value, so that this axis has no physical meaning). The small blue dots show the LD value of all local minima and they exhibit a step structure with some steps indicating a big increase of the LD value. From all the minima, there is one outstanding value (highlighted as a red dot with an LD value of $\mathcal{L}=2.4717286$). This value stands out because it is, by far, the smallest one with the second-smallest minimum exhibiting an LD value that is already three times as large.

To obtain a precise understanding of the background and the occurrence of the observed LD structure as well as their values at the minima, we present in Fig. 4 a selection of ten typical trajectories which are obtained from LD minima. (Note that it is here not important which of the trajectories belongs to which minimum, because the general behavior of the trajectories explains the LD structure and that is the same in all the cases.) We see in Fig. 4 that the trajectories corresponding to minima of the LD are characterized by a combination of small oscillations at the barrier (oscillations with amplitude $x<1$) in addition to large oscillations in the wells (oscillations with amplitude $x>2$). The periodic TS trajectory is also shown as a black line for comparison. Each of these oscillations corresponds to a certain arg length of the trajectory and if we denote by $\mathcal{L}_{\text{bo}}$ the arc length of one oscillation at the barrier as well as by $\mathcal{L}_{\text{wo}}$ one oscillation in the well, the complete LD value corresponding to one of its minima is

[TABLE]

with $m,n=0,1,2,\ldots$ The ‘approximate equal’ sign is here intended to take into account that small deviations from the exact combinations occur in practice, because the single oscillations need to be connected smoothly according to the underlying dynamical equations. Equation (9) directly explains the self-similar structure of the LD because any linear combination of the LD contributions $\mathcal{L}_{\text{bo}}$ and $\mathcal{L}_{\text{wo}}$ leads to a (local) minimum if those values themselves correspond to local minima. In addition, this relation also explains the step structure presented in Fig. 3(b) as a result of different integers $m,n$ while the small slope seen within some of the steps comes from fulfilling the dynamical boundary conditions between barrier and well oscillations. In addition, Fig. 4 also makes clear that we have two important time scales in our system: one of them corresponds to the period of the barrier top (fast, small-amplitude oscillations) and the other one is the recurrence time of a trajectory (slow, big-amplitude oscillations) which is roughly four times as large.

We have, so far, investigated and explained the self-similar LD structure close to the barrier top in the system with finite reactant and product wells, and we have seen that it is the properties of the underlying trajectories which result in these observations. As we have already mentioned in the introduction, it is a major purpose of the LD method to provide a construction scheme for finding the TS trajectory. With the results previously presented in this work, i. e. the occurrence of a huge number of local LD minima close to the barrier region (or formally an infinite number in the limit $\tau\to\infty$), one can easily imagine that this can result in significant problems in the application of the method. Also, the result from Fig. 3(b) that the desired local LD minimum is that with the smallest LD value, i. e. the global minimum in the barrier region, is only of minor help, because of the following reasons: First, its systematic search then requires global optimization procedures, which are not easy to apply, especially in context with the observed self-similar minimum structure. Second, only the comparison of the LD values from a selection of numerically obtained points is not sufficient, because there is no guarantee that the desired minima is part of the selection so that there cannot be another minimum with an even smaller LD value.

In the following, we show how the problem of additional (local) minima is solved using the modified LD definition according to Eqs. (7) and (8). A comparison between the corresponding LD phase space portraits is presented in Fig. 5 using the standard LD definition in the top row and the modified definition with variable integration time in the bottom row. As can be seen in the top row and as we have discussed above, a complicated structure of the stable and unstable manifolds or the LD valleys is present with the standard definition of the LD, Eq. (4). This structure becomes more and more complicated with an increasing number of local LD minima if the integration time $\tau$ is increased (left to right for the values $\tau=10,20,30,40$). The reason for this is that longer integration times allow for more reflections (and therefore global recrossings) which, then, induce the details in the substructure.

The occurrence of the complicated LD structure is immediately suppressed if the variable integration times (8) are used as shown in the bottom row. Here, only single, minorly curved lines can be observed for the manifolds $\mathcal{W}_{\text{s,u}}$. Moreover, their intersection is clearly visible and not accompanied by additional, local minima. This observation especially holds independently of the maximum integration time $\tau$, so that the LD surface does not get more complicated for increasing integration time. Note that it is a consequence of the trajectory cut-off that the LD surface apart from the intersection of the manifolds is now characterized by very small LD values while they naturally get larger if one approaches the intersection. It therefore appears from the color map of the figure that the desired point has turned into a maximum of the LD. However, the relation (6) is still true because there is a very sharp minimum at the intersection of the manifolds (which has an extension smaller than the resolution of the figure). We emphasize that the cut-off of the integration time as defined in Eq. (8) does not have any effect on the position of the global LD minimum, i. e. the initial conditions of the TS trajectory, because this trajectory does not leave the barrier region so that the cut-off is never applied to the TS trajectory.

The TS trajectory which corresponds to the single minimum shown in the bottom row of Fig. 5 is finally presented as the solid blue line in Fig. 6 (the dashed lines shows the time evolution of the barrier top). As expected, the trajectory does not leave the barrier region but remains within the range of the barrier amplitude ($-0.4\leq x\leq 0.4$) which verifies the success of the procedure.

IV Conclusion and outlook

In this paper, we have investigated the phase space structure of a time-dependent double-well potential with special regard to the local minima of the LD and their connection to the TS trajectory. We have seen that reflections of the particle at the potential walls lead to global recrossings and therefore to a self-similar LD structure with a huge number of local minima (or formally an infinite number in the limit $\tau\to\infty$). We have demonstrated that this structure is directly related to the trajectories with each local minimum corresponding to a linear combination of barrier and well oscillations. This complicated structure naively appears to be a significant obstacle in the application of the LD formalism to the construction of TS trajectories. It can be overcome, however, by a simple modification of the LD definition in a way that the underlying trajectories are cut off as soon as they leave the barrier region. As a consequence, only local crossings of the barrier are taken into account. This improved formalism is especially important to the application of the LD-DS method to systems in high-dimensional phase space because the number of global recrossings—viz. oscillations—increases with dimensionality.

The chaotic and self-similar behavior of this driven system echo those seen in atom-diatom reactions by Tyapan and Jaffe.Tiyapan and Jaffé (1994, 1995) As in their work, the DS on the TS trajectory is a homoclinic tangle that can be characterized as a fractal tiling. Future work could and should obtain the scaling law for the tiling and use it to obtain renormalized rates from the sum of fluxes on each tile. Such an approach holds promise in obtaining the rate using the LD-DS scheme locally while addressing increasingly higher-dimensional and more complex barriers.

The LD-DS method presented in this paper can also be used to revisit the earlier uses of the LD. For example, in Ref. Mancho et al. (2013) the fractal-like structures we found in our work is already slightly visible as one cannot identify a single intersection of the manifolds near the saddle. The modification of the LD needed here to obtain the TS trajectory and clearly reveal the fractal structure should be applicable to their use of the LD for obtaining fixed dividing surfaces when the latter LD is beset by the existencs of many local minima due to recurrences.

Acknowledgements.

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 programm under Grant Agreement No. 734557. We thank reviewer 2 for pointing out the approach by Tiyapan and JafféTiyapan and Jaffé (1994, 1995) to resolve the scaling laws of the fractal tiling of the dividing surface as a possible future dirction for resolving rates using the LD-DS method.