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

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.^1–8 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)^9–21 provides a powerful, though usually approximate, framework to calculate the rate from the reactive flux through 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 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.^22–34 By contrast, in time-dependent systems, the DS is, in general, also time-dependent and the transition state (TS) trajectory,^35–41 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,^24,36 but the approach does not have an obvious zeroth-order reference in barrier reactions⁴² 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^43,44 (LDs) provides a general and powerful construction scheme to obtain the TS trajectory in such difficult cases.^40,42,45,46 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.⁴⁷ 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 minimum 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. This time-dependent structure carries the surface that divides the reactant and product regions by identifying the nonrecrossing trajectories useful for obtaining exact rates through the reactive flux time correlation function. The latter technique has been used routinely in the chemistry community and is beautifully imposed in the recent work of Mesele and Thompson⁴⁸ 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 naïve LD-DS approach.

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

that has already been the subject of the work in Ref. 47. This potential consists of two Morse potentials $VMorse±$ providing the reactant (−) and product (+) wells and a time-dependent Gaussian potential $VGauss$ serving as the time-dependent barrier. Specifically, these dimensionless potential terms are

where $xb$ is the position of the barrier, the factor D determines the energy scale, and $±x0$ are the positions of the wells. As in Ref. 47, we use D = 1 for the energy scale, and specify the spacial scale of the system by setting x₀ = 3 as well as a = b = 1. Moreover, we apply a sinusoidal driving of the barrier according to $xb(t)=0.4sin(πt)$ so that the barrier oscillation is $0.4∕3≈13%$ of the distance to the well’s minimum. For a centered barrier top at $xb(t)=0$⁠, the barrier height is

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₀ of the well minima and the barrier height $E‡$ as a reference for the energy scale, this defines the dimensionless units of the system.

As shown in Refs. 42, 45, and 46, a nonperturbative approach for 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

where $v$ is the velocity of a certain trajectory and therefore $L$ is a measure of the trajectory’s arc length over the time interval $[t0−τ;t0+τ]$⁠.

The LD is of special importance for the reaction dynamics because the stable and unstable manifolds $Ws,u$ attached to the (time-dependent) barrier top are directly related to LD’s forward (f: $t0≤t≤t0+τ$⁠) and backward (b: $t0−τ≤t≤t0$⁠) contributions, respectively, according to^42,45,46

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 “$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 $T(t0)$ at time t₀ are located on the intersection of these manifolds,^42,45,46 $T(t0)=Ws(t0)∩Wu(t0)$⁠, so that it is directly related to the two-sided LD (4) via

As we demonstrate below, definition (4) of the LD with fixed integration time $τ$ in a closed system with finite reactant and product wells has the disadvantage that the LD, in general, exhibits not only a single minimum but also 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 the reactant and product. Thus, one needs to define the reaction region through which trajectories traverse from the reactant to product regions in order to obtain the direct rate 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 $τ$ with a variable depending on the underlying trajectory. For this purpose, we redefine the time interval over which the trajectory is integrated through the map

Here, the values $τ±[x(t)]$ are trajectory-dependent integration times which we define according to

with an appropriate “size” $x¯$ of the TS region. The redefinition of the integration time in Eq. (8) is thereby limited by both $τ$ and the time for the particle to leave the barrier region (when |x| is greater than $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 recurrences in the trajectories being globally reflected to the barrier region. By construction, our redefinition precisely retains only the trajectories that contribute to the direct rate. The value of $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 $τ±$ may also differ in the forward and backward time directions.

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 $τ=20$ and 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 start and end either 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 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^46,47 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 to avoid overloading the figure. The minima are located either on one of the manifolds $Ws,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.

Not only taking into account relation (6) but also looking at the actual LD values of the minima can help 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 $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 larger.

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 here it is 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 the 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 arc length of the trajectory and if we denote by $Lbo$ the arc length of one oscillation at the barrier as well as by $Lwo$ one oscillation in the well, the complete LD value corresponding to one of its minima is

with $m,n=0,1,2,…$⁠. The “approximate equal” sign here is 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 $Lbo$ and $Lwo$ 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 it 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 larger.

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 $τ→∞$⁠), 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 are 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 $τ$ is increased (left to right for the values $τ=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 $Ws,u$⁠. Moreover, their intersection is clearly visible and not accompanied by additional, local minima. This observation especially holds independently of the maximum integration time $τ$ so that the LD surface does not get more complicated for increasing integration time. Note that it is a consequence of the trajectory cutoff 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, 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 cutoff 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 cutoff 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 show 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≤x≤0.4$⁠) which verifies the success of the procedure.

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 $τ→∞$⁠). 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 naïvely 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 echoes those seen in atom-diatom reactions by Tiyapan and Jaffé.^49,50 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. 44, 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 the use of the LD for obtaining fixed dividing surfaces when the latter LD is beset by the existence of many local minima due to recurrences.

A.J. acknowledges the Alexander von Humboldt Foundation, Germany, for support through a Feodor Lynen Fellowship. R.H.’s contribution to this work was supported by the National Science Foundation (NSF) through Grant No. CHE-1700749. This collaboration has also benefited from the support by the people mobility programs and most recently by the European Unions Horizon 2020 research and innovation programme under Grant Agreement No. 734557. We thank the reviewer for pointing out the approach by Tiyapan and Jaffé^49,50 to resolve the scaling laws of the fractal tiling of the dividing surface as a possible future direction for resolving rates using the LD-DS method.