Invariant Manifolds and Rate Constants in Driven Chemical Reactions

TL;DR

This paper develops advanced methods for accurately constructing invariant manifolds in multidimensional driven chemical reactions, utilizing machine learning techniques like neural networks and Gaussian process regression to improve rate constant calculations.

Contribution

It introduces the use of Gaussian process regression alongside neural networks for constructing the NHIM in multidimensional systems, enhancing previous methodologies.

Findings

Gaussian process regression effectively constructs smooth NHIMs.

Methods accurately predict rate constants in complex multidimensional models.

Comparison shows machine learning approaches outperform traditional techniques.

Abstract

Reaction rates of chemical reactions under nonequilibrium conditions can be determined through the construction of the normally hyperbolic invariant manifold (NHIM) [and moving dividing surface (DS)] associated with the transition state trajectory. Here, we extend our recent methods by constructing points on the NHIM accurately even for multidimensional cases. We also advance the implementation of machine learning approaches to construct smooth versions of the NHIM from a known high-accuracy set of its points. That is, we expand on our earlier use of neural nets, and introduce the use of Gaussian process regression for the determination of the NHIM. Finally, we compare and contrast all of these methods for a challenging two-dimensional model barrier case so as to illustrate their accuracy and general applicability.

Figures

Tables

Table 1. Table 1: Nomenclature of the different reactive (II, IV) and non-reactive (I, III) areas in phase space separated by the two manifolds. The particle leaves the neighborhood of the saddle at the exit point x exit b , f superscript subscript 𝑥 exit b f x_{\mathrm{exit}}^{\mathrm{b,f}} when propagated in backward (b) and forward (f) time, respectively.

Area	$x_{\mathrm{exit}}^b$	$x_{\mathrm{exit}}^f$
I	$x_P$	$x_P$
II	$x_R$	$x_P$
III	$x_R$	$x_R$
IV	$x_P$	$x_R$

Table 2. Table 2: Training and function call times for the use of NN and GPR to learn and access the same NHIM for the model problem discussed in the text.

	NN	GPR
obtaining training data	$3636s$	$545s$
training	$4220s$	$9s$
optimize hyper-parameters	–	$2238s$
total training time	$7856s$	$2792s$
function call time	$14\mu s$	$253\mu s$

Equations

$$\begin{array}{rl}x& <x^{\mathrm{DS}}({\mathbf{y}}^{\mathrm{p}},{{\mathbf{v}}_{\mathbf{y}}}^{\mathrm{p}},t)\to \mathrm{reactant},\mathrm{or}\\ x& >x^{\mathrm{DS}}({\mathbf{y}}^{\mathrm{p}},{{\mathbf{v}}_{\mathbf{y}}}^{\mathrm{p}},t)\to \mathrm{product}.\end{array}$$

$$N_{\mathrm{react}}(t)-N_{\infty }\propto \exp (-kt),$$

$$V(x,t)$$

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

$$\tilde{\mathcal{L}}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)$$

$$x\left(t_{\text{b,f}}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)\right)$$

$${\mathcal{T}}_{\text{d}}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)$$

$${\mathcal{T}}_{\text{d,f}}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)$$

$${\mathcal{T}}_{\text{d,b}}({\mathbf{x}}_0,{\mathbf{v}}_0,t_0)$$

$$V(x,y,t)$$

$$c=a\left(\sum _{i=1}^N{w}_i{c}_i+b\right),$$

$$k({\mathbf{z}}_1,{\mathbf{z}}_2)=\exp \left(-\frac{{\left|{\mathbf{z}}_1-{\mathbf{z}}_2\right|}^2}{2l^2}\right),$$

$$\Sigma =\left(\begin{array}{cc}K& N^{\text{T}}\\ N& M\end{array}\right)$$

$$K_{ij}$$

$$M_{ij}$$

$$N_{ij}$$

$$\bar{\mathbf{\mu }}=N{K}^{-1}{\mathbf{f}}_{\text{train}}$$

$$\bar{\Sigma }=M-N{K}^{-1}{N}^{\text{T}}.$$

$$\Sigma =M-\tilde{N}{\tilde{N}}^{\text{T}},$$

Full text

Here, we extend our recent methods by constructing points on the NHIM accurately even for multidimensional cases. We also advance the implementation of machine learning approaches to construct smooth versions of the NHIM from a known high-accuracy set of its points. That is, we expand on our earlier use of neural nets, and introduce the use of Gaussian process regression for the determination of the NHIM. Finally, we compare and contrast all of these methods for a challenging two-dimensional model barrier case so as to illustrate their accuracy and general applicability.

1 Introduction

The framework of transition state theory (TST) provides a powerful basis for the qualitative and quantitative description of chemical reactions. When their dynamics can be described by a Born-Oppenheimer potential driven by a classical equation of motion, then TST reduces the dynamical calculation to a geometric one involving the identification of a barrier region that separates reactants and products. The position of the barrier is typically marked by a saddle point of rank 1 that has exactly one unstable direction which coincides with the reaction coordinate at that point. The remaining degrees of freedom are locally stable and are associated with the bath coordinates. In this paper, we aim to address activated processes with arbitrary dimensionality but restricted to reactions that can be characterized by a one-dimensional reaction coordinate. We do not solve this class of problems in fullest generality, but do make progress in treating systems with bath coordinates of dimensions higher than one, addressing both theoretical and numerical challenges to the computation of the dividing surface.

TST rests on the identification of a dividing surface (DS) in the barrier region which separates reactants from products.1, 2, 3, 4, 5 It is exact, if the DS is crossed once and only once by each reactive trajectory. Recrossings of the DS lead to an overestimation of the rate. Advances in the determination of a recrossing-free DS will therefore impact a broad range of problems beyond those of chemical reactions in which the overall process can be described as an activated process that is primarily characterized by a single, albeit curved, reaction coordinate in the presence of a bath. Examples abound in atomic physics 6, solid state physics 7, cluster formation 8, 9, diffusion dynamics 10, 11, cosmology 12, celestial mechanics 13, 14, and Bose-Einstein condensates,15, 16, 17, 18, 19 to name a few.

Of course, the primary aim for the current work is the construction of accurate rates in chemical reactions that take place in solution, perhaps under nonequilibrium conditions, and for which we ultimately want to achieve some level of control. The basic (or naive) theory uses a planar dividing surface to determine the flux, but it is approximate because it is invariably recrossed by trajectories that go from reactants to products.1, 2, 3, 4 It is also an upper bound, and leads to a variational transition state theory20, 21 which optimizes the dividing surface by way of minimizing recrossings. The development of semiclassical transition state theory22, 23, 24, 25, 26 taking advantage of Keck’s phase space representation of TST27 hinged on the development of good action-angle variables associated with the DS. The use of perturbation theory to construct these objects was found in parallel to work in dynamical systems theory addressing activated escape.28, 29 The theory was then generalized for chemical reactions under time-dependent conditions —arising from driving, noise, or both.30, 31 Applications to chemical reactions have included H $+$ H2,32, 33 LiCN,34 and ketene isomerization.35 While these are interesting cases, they do not represent the bulk of chemical reactions which are generally higher dimensional and which take place in more complex environments. This work is thus focused on advancing the theoretical and computational machinery that have the potential of addressing such challenging systems.

In a multidimensional autonomous Hamiltonian system, a recrossing-free DS is attached to the normally hyperbolic invariant manifold (NHIM). The latter contains all trajectories that are trapped in the saddle region both forward and backward in time. It can be constructed using a normal form expansion.36, 37, 24, 32, 38, 13, 39, 40, 41, 42 For a multidimensional time-dependent system that is driven by an external field or subject to thermal noise, the question arises whether the NHIM can be generalized in such a way that it can still be used to construct a recrossing-free DS. Recent successful constructions43, 44, 45 of such a DS suggests that it is indeed possible to do so. The purpose of this paper is to clarify the geometry that underlies this construction, and to present a readily applicable procedure to calculate rate constants of time-dependent, driven and multidimensional systems with a rank-1 saddle and an arbitrary number of bath modes. As this work is focused on advancing the methods so as to ultimately treat higher dimensional chemical reactions, numerical results are restricted to model systems —inspired by chemical systems— useful for verifying the theory.

The computational task is challenging because of the dimensionality of the objects: In an $n$-dimensional system with one reaction coordinate and $n-1$ bath coordinates, the phase space has dimension $2n$. The DS is a time-dependent and non-trivially curved hypersurface of dimension $2n-1$, and the NHIM is a surface of dimension $2n-2$. Previous work has shown that attaching a DS to the NHIM, whose dimension is lower, can yield a recrossing-free DS if the reaction coordinate is unbound.44, 45

An individual point on the NHIM can be computed with the algorithm discussed in Ref. 46. However, to construct a DS that covers a sufficiently large region in phase space and time, a vast number of points is required. Even with an efficient algorithm, computing the required number of points on the NHIM is a daunting task. Consequently, a continuous interpolation between the finite collection of high-accuracy —and therefore computationally expensive— points on the DS in space and time is needed to determine the reactivity of a given trajectory. We present two tools developed in the context of machine learning to approximate any function between isolated points in arbitrary high dimensions. It is notable that Pozun et al.47 earlier implemented a support-vector machine to identify a DS between reactants and products on a time-independent potential energy surface. Here, we focus on time-dependent potentials, often coupled to a bath, in which the associated NHIM is consequently time-dependent. First, a neural network (NN)48 can be trained on an arbitrary set of points located on the NHIM. Once trained, these networks provide fast access to any point needed for calculating rates in time-dependent driven systems.45 We also present an alternative machine learning algorithm based on Gaussian process regression (GPR).49, 50 We compare these two algorithms and discuss the relative advantages and disadvantages in their application to TST.

Perturbative51, 52, 53, 54, 55 and nonperturbative43, 35, 56, 57, 58 constructions of a recrossing-free DS in time-dependent systems rely on the concept of the transition state (TS) trajectory:30, 59 a unique trajectory bound to the vicinity of the saddle for all time. In systems with more than one degree of freedom, such a TS trajectory is unique only if the system is dissipative with the oscillations of the bath modes damped asymptotically. If the reactive system is Hamiltonian —that is, when in the absence of damping,— all trajectories on the NHIM are trapped near the barrier for all time and can therefore serve as TS trajectories. The consequences of the non-uniqueness of the TS trajectory were avoided in Ref. 60 by an arbitrary choice of one of the trapped trajectories as “the” TS trajectory. The DS constructed there and in Refs. 44, 45 also avoids such arbitrary choices because it considers the NHIM itself, rather than a particular trajectory on it. For multidimensional time-periodically driven Hamiltonian systems, we suggest identification of the TS trajectory to be restricted to the one trajectory, among all the bound trajectories on the NHIM, that is periodic. This reduces to the usual result in one dimension, and is no longer arbitrary in higher dimensions as it is now selected by a characteristic property of the system.

The outline of the paper is as follows: In Sec. 2, we start with a brief overview describing how rates in time-dependent systems with a rank-1 saddle are calculated using the NHIM and the associated DS. This suggests a restricted definition of the TS trajectory stated above. In Sec. 3, we present a fast and robust algorithm for calculating single points on the NHIM. Sec. 4 summarizes two interpolation algorithms based on machine learning, and provides an analysis of the respective advantages and disadvantages. In Sec. 5, we apply both methods to calculate rate constants for a two-dimensional system.

2 TST in time-dependent systems

A possible first step in obtaining the time-dependent representation of the DS is the determination of individual points on this high-dimensional object. As outlined in Ref. 44, in time-dependent systems it suffices to calculate the normally hyperbolic invariant manifold (NHIM) because the DS in time-dependent systems is attached to it in a prescribed way. In this section, we describe the former as a fundamental geometric object and demonstrate how it can be used to obtain a DS. We start by considering a one-dimensional system in which there are no bath coordinates. We then present the binary contraction method for finding the NHIM based on the geometric properties of phase space in the vicinity of a rank-1 saddle. Furthermore, we describe how the insights gained in the one-dimensional setting can be generalized to higher dimension. Finally, we will discuss how a unique TS trajectory can be defined in a multidimensional system.

2.1 The role of the NHIM in TST

In a one-dimensional system, a rank-1 saddle is simply a maximum of the potential. The relevant phase space structures in its neighborhood are shown in Fig. 1: Trajectories that approach the energy barrier with sufficiently high energy (from either side) will cross the barrier, and others will not. Its phase space can be separated into four different regions as characterized by where particles came from and where they go to. The dividing lines between these regions are the stable $\mathcal{W}_{\mathrm{s}}$ and unstable $\mathcal{W}_{\mathrm{u}}$ manifolds. They contain trajectories that are trapped near the saddle in forward ($\mathcal{W}_{\mathrm{s}}$) or backward ($\mathcal{W}_{\mathrm{u}}$) time, respectively.

A trajectory starting at a given phase space point is propagated in forward (f) and backward (b) time until it leaves the saddle region given by $x_{\mathrm{R}}\leq x\leq x_{\mathrm{P}}$ with a suitably chosen $x_{\mathrm{R}}$ at the reactant side and $x_{\mathrm{P}}$ at the product side of the reaction coordinate. The initial point is then assigned to a region, depending on whether it leaves the saddle region at $x_{\mathrm{R}}$ or $x_{\mathrm{P}}$ in forward or backward time, as explained in Table 1. Note that numerically, the trapping cannot, in general, be observed because the dynamics near the saddle is unstable. Small numerical errors cause an otherwise trapped trajectory to leave the saddle region in either time direction.

The closures of stable and unstable manifolds intersect in a point (see black cross in Fig. 1). When the potential is also time-independent, this intersection is a fixed point. In time-dependent (or driven) cases, it becomes a trajectory that is trapped in both forward and backward time. This trajectory is referred to as the TS trajectory 30, 59, 51, 52, 53, 54, 55, 56, 57, 58. In one-dimensional systems, it is unique. In multidimensional systems, the initial conditions of trajectories trapped in both forward and backward time form the time-dependent multidimensional NHIM61, 62, 63. The latter reduces to a single trajectory in one-dimensional systems, and so the unique TS trajectory in one dimension can also be called the NHIM.

For a time-independent, one-dimensional system, the NHIM is located exactly at the barrier top. In a time-dependent system, all structures shown in Fig. 1 depend on the initial time at which the trajectories are started, i.e., the stable and unstable manifolds, as well as the NHIM, are themselves time-dependent. The NHIM can be expected to reside close to the saddle, but it will in general not coincide with it. If the external driving is periodic in time, the NHIM is a periodically moving manifold with the same period as the external driving force (see blue, dashed circle in Fig. 3(d), calculated for the example of Eq. (3)).

So far, we have assumed open reactant and product basins. In the more typical case of closed basins —as in e. g. the model potential of Ref. 64, or potentials of real systems, for instance LiCN 65, 34, 66 or ketene 67, 68, 69, 35— the requirement for trajectories to be recrossing-free often needs to be restricted to mean locally recrossing-free, i.e., a particle that enters the neighborhood of the saddle will cross the DS no more than once before it leaves the saddle region. It may, however, reenter this region and recross the DS at a later time. For further information about these global recrossings and the degree to which TST can remain exact despite them, see Ref. 64.

2.2 Recrossing-free DS in time-dependent systems

For a system with $n$ degrees of freedom, the separation of the $2n$-dimensional phase space into reactant and product regions requires a surface of dimension $2n-1$. To distinguish reactant and product regions for any rank-1 saddle in an $n$-dimensional energy surface we use an appropriate coordinate system where $x$ approximates the reaction coordinate, i.e., the unstable direction of the saddle, and $\boldsymbol{y}$ the $n-1$ remaining (bath) coordinates, given by the stable directions of the saddle. The corresponding velocities are denoted by $v_{x}$ and $\boldsymbol{v_{y}}$, respectively. We will construct a dividing surface at an instantaneous reaction coordinate $x^{\mathrm{DS}}(\boldsymbol{y}^{\mathrm{p}},\boldsymbol{v_{y}}^{\mathrm{p}},t)$ that depends on the bath coordinates, their velocities, and on time. A phase space point $(x,\boldsymbol{y},v_{x},\boldsymbol{v_{y}})$ is classified as lying on the reactant or product side of the barrier at time $t$ depending on whether it fulfills

[TABLE]

While propagating an ensemble of particles, each particle can be classified as reactant or product at any time according to Eq. (1). When it crosses the moving DS, it reacts and contributes to the reaction of the time-dependent system.

To satisfy the recrossing-free requirement, this high-dimensional DS is generally non-trivially curved in a manner that depends on the details of the dynamics. The representation (1) of the DS assumes that no curvature in the $v_{x}$ direction is required. This assumption is inspired by the phase space plot in Fig. 3, where such choice of a DS appears as a vertical line going up from the NHIM. Reactive trajectories cross the DS transversely and exactly once, whereas non-reactive trajectories do not cross it. This statement is trivially true in a harmonic system. It describes a qualitative property that is robust under perturbations of either the system or the dividing surface. Many other surfaces that are sufficiently close to the vertical line will also be recrossing-free; the exact choice of DS is therefore not critical. We take advantage of this non-uniqueness by postulating that the DS has the form (1) and that its location $x^{\mathrm{DS}}=x^{\mathrm{NHIM}}$ coincides with that of the NHIM. These assumptions single out a unique DS and also provide a means to compute it because the NHIM is amenable to numerical computation, as discussed below in Sec. 3.

With this simplification, our choice of a recrossing-free DS, independent on the reactive velocity $v_{x}$, can be justified via a confirmation of the absence of unwanted recrossings in dynamical simulations, as shown in Sec. 5 for a particular two-dimensional model system. However, in more general cases when particles cross the DS far away from the NHIM, this simplification might lead to unwanted recrossings (see, e. g. Refs. 70, 71, 72, 40), as the assumption embodied in Eq. (1) has its limits. There is no general proof that a recrossing-free DS of the form (1) should exist, and exceptions are known.70 If the invariant manifolds are strongly deformed, e. g. such that they cross the chosen DS, then the fate of a trajectory starting on the DS $x=x^{\mathrm{DS}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ will depend on $v_{x}$, and the representation (1) becomes impossible. However, this is expected only for strongly non-linear systems or very high velocities, so that the representation (1) is adequate for many systems of practical relevance. For a given system, one can always numerically confirm whether a DS is truly recrossing-free. This may be done by propagating a large number of trajectories with initial conditions in the vicinity of the moving saddle, and subsequent monitoring of the crossings across the constructed DS, as has been done successfully for the model system discussed in Refs. 44, 45.

If the underlying potential has open reactant and product basins, a particle that crosses a recrossing-free DS will never return because it escapes the saddle region forever. Assuming that our choice of a DS independent on the reactive velocity $v_{x}$ is sufficiently accurate The propagation of a reacting particle can be stopped when a particle crosses a space determined by $x$ = $x^{\mathrm{DS}}$ while calculating rates in open systems. In systems with closed reactant or product basins, particles can be reflected globally at the boundaries of the basins and may therefore cross the space determined by $x=x^{\mathrm{DS}}$ infinitely many times.64, 65, 34, 66, 67, 68, 69, 35. In this case, the propagation cannot be stopped and particles can react from one side of the saddle to the other multiple times. Still, the time-dependent classification in reactant and product is always given due to the local recrossing-free property of the DS presented here (see Sec. 2.2 for further discussion on how to construct the DS and its properties).

Reacting particles can, in principle, cross the DS at some previously unknown values of the bath coordinates and velocities in between the pre-calculated points, and at a previously unknown time $t$ (schematically visualized as trajectory in Fig. 2), while the DS itself is a time-dependently moving object. To obtain a rate from the simulation of a large number of trajectories, one therefore needs to evaluate the location $x^{\mathrm{DS}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ of the DS for a wide range of the parameters $\boldsymbol{y}$, $\boldsymbol{v_{y}}$, $t$, and the parameter space will typically have high dimension. However, even the calculation of $x^{\mathrm{DS}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ at single points is challenging and numerically expensive. Indeed, it is hopeless to calculate $x^{\mathrm{DS}}(\boldsymbol{y}=\boldsymbol{y}^{p},\boldsymbol{v_{y}}=\boldsymbol{v_{y}}^{p},t)$ from scratch for every position $(\boldsymbol{y}^{p},\boldsymbol{v_{y}}^{p})$ on the trajectory of every particle $p$ at the respective time $t$, while propagating a large ensemble. Consequently, it is critical to obtain a representation of the function $x^{\mathrm{DS}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$, that describes the dependence of the DS on all its arguments over a sufficiently large range, and which is easy to evaluate.

We obtain a computationally efficient representation of the DS in two steps: First, a high-accuracy algorithm is used to obtain a small, but incomplete, number of representative points on the DS. In Sec. 3, we present an algorithm that is able to efficiently calculate these points, based on the geometric properties of phase-space in the saddle region. Second, a smooth representation of the function $x^{\mathrm{DS}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ is obtained through an efficient extrapolation of the finite points obtained in the first step. This requirement can sometimes be achieved through interpolation and fitting of expected forms of the potential. To avoid such numerics and to take advantage of the limited data set presently available, we have recently implemented machine learning approaches which should, in principle, be applicable to arbitrarily high dimension 45. The use of these methods to obtain a time-dependent description of the full DS is summarized, extended and assessed in Sec. 4.

2.3 Obtaining rate constants from reactant decay

Through TST one aims to acquire rate constants for a chemical reaction from a reactant state over an energy barrier to a product state. A typical approach for obtaining rate constants 73 is to propagate an ensemble of particles, identify only those trajectories from reactants to products that cross the DS such as that illustrated in Fig. 2, sum their flux, and divide by the population of reactant particles. In the TST approximation the reacting flux is obtained by summing the instantaneous flux across the entire DS. Rates thus obtained are exact if and only if the DS is free of recrossings.

To obtain rate constants, we monitor the number of trajectories that remain reactants over time $N_{\mathrm{react}}(t)$ before crossing the moving DS. For large times, when all trajectories have left the barrier region, it approaches the number of non-reactive particles $N_{\infty}$. After the initial transient behavior has decayed, this limit is approached exponentially

[TABLE]

where $k$ is the rate constant of this reaction.74, 3, 75 Deviations from exponential behavior also arise at long times when the number of remaining trajectories becomes too small to be statistically meaningful.

3 Obtaining single points on the NHIM

3.1 Revealing phase space structures

The phase space structures discussed in Sec. 2.2 can be uncovered through an increasing array of tools, including the ones reviewed in this section. For simplicity, we will illustrate them using a one-dimensional model system for a simple chemical reaction with a periodically moving barrier and open reactant and product channels. Specifically, we consider the motion of a particle of unit mass under the influence of the time-dependent potential

[TABLE]

Here, $E_{\mathrm{b}}$ is the height of a Gaussian barrier oscillating along the $x$ axis with frequency $\omega_{x}$ and amplitude $\hat{x}$. We use dimensionless units in which we set $E_{\mathrm{b}}=2$, $\hat{x}=0.4$, and $\omega_{x}=\pi$.

A tool to resolve the manifolds in phase space is the Lagrangian descriptor (LD).76 In the context of TST, the LD at position $\boldsymbol{x}_{0}$ in the coordinate space, velocity $\boldsymbol{v}_{0}$, and time $t_{0}$, is defined as 76, 43, 56, 58

[TABLE]

The LD measures the arc length of the trajectory $\boldsymbol{x}(t)$ that passes through $(\boldsymbol{x}_{0},\boldsymbol{v}_{0})$ at time $t=t_{0}$. The trajectory is considered forward and backward in time over the interval $[t_{0}-\tau;~{}t_{0}+\tau]$, and the parameter $\tau$ is chosen such that the relevant time scale of the system is covered. In the case of open reactant and product basins, the stable $\mathcal{W}_{\mathrm{s}}$ and unstable $\mathcal{W}_{\mathrm{u}}$ manifolds associated with the energy barrier should be revealed by local minima of the LD in phase space, but are sometimes obfuscated by ditch-like structures, as seen in Fig. 3(a). Trajectories starting near the stable and unstable manifolds generally have low LDs because they are bound to the saddle region for a finite time. One ditch —that is, a deep minimum— is found in the forward time component, indicating the location of the stable manifold $\mathcal{W}_{\mathrm{s}}$ and the other ditch is found in the backward time component, indicating the unstable manifold $\mathcal{W}_{\mathrm{u}}$.

In the case of closed reactant and product basins, the LD defined in Eq. (4) results in a fractal-like structure in phase space that also contains both manifolds, but obfuscates them even more 64. Nevertheless, $\mathcal{W}_{\mathrm{s}}$ and $\mathcal{W}_{\mathrm{u}}$ can be revealed if LDs are computed not for a fixed time $\tau$, but rather until they leave the saddle region $x_{\text{R}}\leq x\leq x_{\text{P}}$:

[TABLE]

As shown in Fig. 3(b), the manifolds are then clearly revealed by local maxima in $\tilde{\mathcal{L}}$.

The shift from minima in $\mathcal{L}$ to maxima in the modified $\tilde{\mathcal{L}}$ is attributed to the fact that trajectories close to the stable and unstable manifold are integrated over longer times, while the respective particles are moving in the saddle region. With this observation, one can find a simpler and arguably more reliable approach to reveal the manifolds. By simply tracking the time a particle spends in the saddle region, we alleviate the need to compute the arc length of the trajectory. Following the nomenclature of the LD, we refer to these times as time descriptor (TD):

[TABLE]

The closer a particle is to $\mathcal{W}_{\mathrm{s}}$ in forward time or to $\mathcal{W}_{\mathrm{u}}$ in backward time, the longer it will stay in the saddle region in the respective time component. Consequently, the stable and the unstable manifold are again revealed by the maxima of the TD, as seen in Fig. 3(c). Note that since trajectories on the stable and unstable manifolds are trapped in the barrier region for infinite time, both $\tilde{\mathcal{L}}$ and $\mathcal{T}_{\text{d}}$ may diverge, although it is unlikely when propagating the trajectories numerically.

To calculate the $v_{x}$-position of the two manifolds in phase space for a given $x$-position (reaction coordinate), one can search for the appropriate extremum of the LD, the modified LD, or the TD. Using what we refer to as the standard LD is perhaps the least robust of these methods because the LD has several local minima in addition to the primary minimum indicating the manifold, see Fig. 4(a). The modified LD and the TD do not show such substructure and have clear maxima for both stable and unstable manifolds, see Fig. 4(b) – (d). However, because of its simpler definition the TD seems to be the method of choice as long as one can clearly define a saddle region. With all of these methods, the integration in forward and backward time can also be performed separately to reveal the two manifolds individually, see Fig. 4(d). Thus, numerical ambiguities between stable and unstable manifolds near their intersection can be avoided.

3.2 1D: Finding the NHIM

With the stable and unstable manifold revealed by the LDs or TDs, it is now possible to find the NHIM. Since the NHIM is given by the intersection of the stable and unstable manifold, it can be computed by finding coordinates where their distance is zero. This can be accomplished with a root search that finds the value of the reaction coordinate $x$ for which the difference between the velocities $v_{x}$ of the stable and the unstable manifold is zero. However, this procedure is very costly as it nests an extreme value search of the LDs or TDs within a root search for the distance between these extreme values. Moreover, the computation of each of the many LD or TD values underlying this search is itself expensive because it requires the integration of a full trajectory.

An efficient alternative for obtaining the NHIM is the binary contraction method,46 which is based on the classification of trajectories in Fig. 3(d). The algorithm starts with a quadrangle whose vertices represent the four different regions, see Fig. 5. A candidate vertex is then associated with the midpoint of one of the edges of the quadrangle, and a trajectory is initiated to determine which region it belongs to. The corresponding vertex is then replaced, and the area of the quadrangle is reduced. This process is repeated for all vertices until the quadrangle converges to the NHIM with the desired accuracy —that is, sufficiently small area. Additional details of the algorithm, including an explanation on how to find good initial conditions and the treatment of exceptional cases, can be found in Ref. 46.

Thus the binary contraction method can be used to find the NHIM (or the TS trajectory) with high precision, while propagating considerably fewer trajectories than the optimization method described above. For example, in the one-dimensional model potential of Eq. (3) with periodic external driving, this leads to a periodic orbit (the TS trajectory), as demonstrated by the blue dashed circles in Fig. 3.

3.3 The NHIM in higher dimensions

In systems with $n>1$ degrees of freedom, the situation is qualitatively similar, but it is more challenging because all phase space structures have higher dimensions. As already established in Sec. 2.2, the NHIM has dimension $2n-2$ in phase space. Consequently, its associated stable and unstable manifolds both have dimension $2n-1$. The time-dependence in the potential, either because of the external driving or the noise, adds additional complexity which makes these manifolds time-dependent.

A point in phase space can be characterized by the reaction coordinates $x$ and $v_{x}$ and the bath coordinates $\boldsymbol{y}$ and $\boldsymbol{v_{y}}$. For any fixed values of the bath coordinates $\boldsymbol{y}$ and $\boldsymbol{v_{y}}$, we obtain a two-dimensional section through phase space with coordinates $x$ and $v_{x}$. We assume that in each of these sections the stable and unstable manifold intersect as shown in Fig. 1. In particular, each section contains a single point of the NHIM that lies at the intersection of the two reactive and two nonreactive regions. This condition provides a unique reaction coordinate and velocity $(x^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t),v_{x}^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t))$ on the NHIM for each bath coordinate and time, which can be found by the binary contraction algorithm. Because these coordinates depend on $2n-2$ variables $\boldsymbol{y}$ and $\boldsymbol{v_{y}}$ in phase space, the NHIM has dimension $2n-2$.

A possible counterexample arises in the simple limiting case of an uncoupled harmonic system. Therein, the point in the reactive phase space, $x^{\mathrm{NHIM}}$ and $v_{x}^{\mathrm{NHIM}}$ will be constant, independent of the bath coordinates $\boldsymbol{y}$ and $\boldsymbol{v_{y}}$. In this case, it would not be possible to represent any of the bath modes as functions of $x$ or $v_{x}$. Nevertheless, the construction above provides the correct phase space point because it does not rely on the association being injective.

A non-trivial multidimensional example is a two-dimensional model with a time-periodically moving barrier separating open reactant and product basins. It can be represented using an extension of the one-dimensional potential of Eq. (3) to

[TABLE]

This potential describes a two-dimensional potential energy landscape that has minima in the bath coordinate along the curve $y=({2}/{\pi})\arctan\left(2x\right)$. A moving Gaussian barrier with height $E_{\mathrm{b}}$ is added over this stationary potential. It is oscillating along the $x$ axis with frequency $\omega_{x}$ and amplitude $\hat{x}$. $\omega_{y}$ is the frequency of oscillations in the bath mode for a particle of unit mass. For simplicity, we again use dimensionless units in which the parameter are $E_{\mathrm{b}}=2$, $\omega_{x}=\pi$, $\omega_{\mathrm{y}}=2$, and $\hat{x}=0.4$.

We choose $x$ as the approximate reaction coordinate and $y$ as the bath coordinate. Consequently, we calculate an $(x,v_{x})$-slice for each set of the remaining coordinates $(y,v_{y})$ and for any time. Using the binary contraction method, we find the intersection ${(x,v_{x})}^{\mathrm{NHIM}}(y,v_{y},t)$. Repeating this procedure for different values $(y,v_{y})$ on an equidistant grid for any time $t$ and using an appropriate interpolation method (see below), we obtain the full NHIM of the system. It is a two-dimensional, time-dependently moving surface embedded in the four-dimensional phase space. A representation of this NHIM at $t=0.5$ is shown in Fig. 6 as a gray surface.

For certain points $(y,v_{y})$, Fig. 6 shows sections of the stable and unstable manifolds that are attached to the NHIM. To obtain them, we consider a two-dimensional section of phase space with fixed values of $y$ and $v_{x}=v_{x}^{\mathrm{NHIM}}(y,v_{y})$. The dynamical structure of this section is similar to that in Fig. 1, except that the axes are now $x$ and $v_{y}$. The manifolds are obtained by finding the value of $v_{y}$ at which, for given $x$, the TD is maximal. The stable and unstable manifolds intersect on the NHIM, as expected. Because they were obtained independently of the calculation that yields the NHIM, this observation also serves to confirm the reliability of the computation and the uniqueness of the NHIM.

3.4 TS-trajectory in higher dimensions

In a one-dimensional system, the NHIM is a single trajectory that can be identified as the TS trajectory. A higher-dimensional NHIM, by contrast, is made up of infinitely many trajectories over time, all of which are trapped in the saddle region for all times. Thus a unique identification of the TS trajectory requires additional constraints beyond the minimal requirement that it is trapped within the saddle region.

In an autonomous system, the location of the DS is often marked by the naive TS —that is the saddle at the top of the barrier.

The TS is the point on some selected reaction path that crosses the DS and it is not necessarily the same as the naive TS but is often in its proximity. We further note that there exist exceptions to cases in which the TS is near the naive TS. This includes the TS structures far from the naive TS found by Gray and Davis77 for the autonomous HeI2 system and, more recently, roaming reactions that are associated to a TS that is not necessarily close to a saddle.78, 68, 79, 80, 81 Such cases are not addressed in the present work.

The TS trajectory is thus intended to generalize the fixed TS for driven systems. In a one-dimensional, time-periodically driven system, the TS trajectory is a periodic orbit. Following this property, we suggest that for a time-periodically driven, multidimensional Hamiltonian system the period-1 orbit trapped in the barrier region, if uniquely identified, should also be called the TS trajectory. Here, period-1 qualifies the orbit as one which is periodic and for which the ratio of its period to that of the external driving is 1.

To make this definition precise, we determine the Euclidean distance $\Delta\gamma$ of two phase space coordinates $\boldsymbol{\gamma}=(\boldsymbol{x},\boldsymbol{v})^{\mathrm{T}}$ on the same trajectory, one time period $T$ of the potential apart. In practice, it can be computed by propagating initial conditions given for a certain time $t_{0}$ forward and backward in time by half a potential period $T/2$. We perform this computation in such a way to avoid errors accumulated by propagating trajectories near unstable manifolds over long times. A trajectory with $\Delta\gamma=\left|\boldsymbol{\gamma}(t_{0}+T/2)-\boldsymbol{\gamma}(t_{0}-T/2)\right|=0$ is a periodic orbit. We identify such a periodic orbit as the TS trajectory if it is the only period-1 orbit bound to the saddle region with respect to the periodicity of the external driving in a time-periodically driven system.

Again, the same assumptions as discussed in Sec. 2.2 are valid, meaning that the system should be not too non-linear, and the velocities of reacting particles not too high. Such high non-linearities lead in general to a break-down of our assumptions as the motion on the NHIM of systems with many degrees of freedom becomes highly chaotic. Examples for such chaotic regimes can be found in the high-energy limit of autonomous systems described in Refs. 70, 71. In such highly non-linear systems, the NHIM and the TS can even bifurcate, as found in Refs. 72, 40.

The distance $\Delta\gamma$ shown in Fig. 7 has been computed on the NHIM at $t_{0}=0$ for the potential (7). The color of any point on the NHIM indicates the distance $\Delta\gamma$ of the phase space vectors of two trajectories, integrated for $T/2=1$ in forward and backward time. A pronounced minimum, and thus by our proposed definition, a point on the TS trajectory, can clearly be distinguished. Using this point as the initial condition for a trajectory propagation, the periodic TS trajectory can be obtained. Since this propagation will fail after a certain time because the NHIM is unstable, one cannot compute a fully periodic trajectory numerically. Nevertheless, with sufficiently accurate initial conditions on the NHIM, one can compute a periodic trajectory up to a certain error tolerance. Alternatively, various points on the TS trajectory can be computed from scratch searching for the minimum of $\Delta\gamma$ on the NHIM calculated at different initial times $t_{0}$.

4 Machine-learning based description of the NHIM

Using the methods described in Sec. 3, a single point on the NHIM can be calculated for arbitrary values of the bath modes and time, as is done here for the illustrative case of system (7). However, for a rate calculation, the position of the DS must be compared to the instantaneous position of each propagated particle at each time step. To compute the required reaction coordinate of the DS from scratch whenever it is needed would be prohibitively expensive because even the efficient binary contraction method described above requires the propagation of dozens of trajectories. To reduce the numerical effort, it is critical to obtain a continuous representation of the DS as a function of bath modes and time that is based on the knowledge of a comparatively small number of points on the NHIM.

Machine learning methods have already been applied within the field of theoretical chemistry for the interpolation of high-dimensional potential energy surfaces.82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92 Following our recent work,45 we use machine learning techniques —and specifically through the use of NNs or GPR— to characterize an entirely different set of surfaces: the NHIM and the DS attached to it. After training on an arbitrary set of points on the NHIM, a neural network yields a continuous representation that approximates the NHIM at the given points and elsewhere as shown in Sec. 4.1. In principle, the method can be applied for arbitrarily high dimension, though the complexity of the network may increase as the number of bath coordinates grows. Alternatively, the interpolation of a given set of points on the NHIM can be performed using GPR (see, e. g. Ref. 49 and references therein) as discussed in Sec. 4.2. This method requires a priori knowledge of the length and time scales along which the position of the NHIM is expected to vary. The comparison between these two methods and a discussion of their advantages and disadvantages in a given context will follow in Sec. 4.3.

4.1 Feed-forward neural networks

Neural networks (NNs) are a powerful tool for approximating functions of arbitrary complexity in any dimension.48 Inspired by our current understanding of NNs in the brain, an artificial NN consists of a set of model neurons with complex connections between them. It transforms a given input into some output. Typically, a NN is trained in the sense that a set of free parameters is adjusted on a set of input points for which the output is known. During training, the net learns to find patterns and dependencies in the given data without requiring any user input beyond the training data. When properly trained, the network will provide a good approximation to the correct output for input similar to the points it was trained on. In our case the training data consists of a set of points on the NHIM. By training on these points, we can create two separate NNs that represent the functions $x^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ and $v_{x}^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ separately, or use a more complex network to represent the combined function $(x,v_{x})^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$. This representation will be reliable in the range of inputs that a suitably chosen set of training data covers, as long as the accuracy of the output was verified.

Here, we provide a brief overview on how we set up and use a NN to represent the DS.45 For a more general discussion on NNs, we refer the interested reader to the literature, e. g. Refs. 93, 94 and references therein. A feed-forward NN consists of several layers through which information is processed, see illustration in the schematic of Fig. 8. As we create a net to approximate ${(x,v_{x})}^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$, the input layer is given by the bath coordinates $\boldsymbol{y}$, velocities $\boldsymbol{v_{y}}$, and time $t$. Based on its input, each neuron computes an output value that is passed to the next layer of the network, from that layer to the next, and so on until the final output is obtained.

In typical neural network applications, a user must decide on the number of hidden layers, the number of neurons in each layer, and the function by which a neuron converts its input into its output. Each neuron —such as the one shown in Fig. 8(b)— of a given layer receives inputs $c_{i}$ from the $N$ neurons of the previous layer. Its output,

[TABLE]

depends on the input values and on the activation function $a(x)$ whose nonlinearity ensures that the NN is not simply a linear transformation of the input variables and hence capable of capturing the necessary complexity. In this work, we choose the activation function $a(x)=\tanh\left(x\right)$ for all neurons except those on the output layer as it is smooth and bound for $x\to\pm\infty$ while providing the requisite nonlinearity. The different weights $w_{i}$ that the neuron assigns to the output of each neuron of the previous layer and the bias $b$ are optimized individually during the training process. We note, that although a single hidden layer NN is a universal approximator, it may not be the most efficient NN to a specific problem. In the following, we tried several different NN configurations, varying the number of hidden layers and neurons. The examples presented below provide the best results among the different approaches we sampled, but may not necessarily be fully optimized.

At the onset of training, a NN is initialized with the weights and biases of the neurons chosen randomly within a reasonable range of values. After successful training, a NN should return a good approximation of the outputs over the domain of inputs. A cost function is selected to measure the quality of the approximation. In the present application, we average the mean square difference between the output of the network and the known values over all training points. This mean square cost function has units of length squared whose values are denoted with respect to the dimensionless units introduced in Sec. 3.3. To minimize the cost function, there are several optimization routines that adjust the weights and biases of the neurons iteratively. This includes the back propagation algorithm which is a gradient-based method that successively adjusts all free parameters.95, 45 Usually, a NN is exposed to the entire set of training data many times before it is optimized, and each such exposure is called an epoch.

The training of a NN is usually not performed with all the available data. Instead, the data set is split into two disjoint sets: a training set and a verification set. The NN is trained only on the training set. The cost function on the verification set is used to assess the reliability of the NN for data that it was not trained for. If the cost function on the training set is significantly smaller than the cost function on the verification set, then the NN adapted very specifically to the given training set. This kind of overfitting is, in general, an undesirable effect as it indicates that the NN is no longer learning to make more accurate predictions.

To illustrate the ease in training a NN to approximate the position $x^{\mathrm{NHIM}}$ for the NHIM in the two-dimensional model system of Eq. (7), we wrote and tested a simple C++ code without using the more generally available libraries that we use in the production runs described below. The network was constructed with two hidden layers: the first and second with $40$ and $10$ neurons, respectively. It was trained on 2000 points randomly distributed over the NHIM and obtained using the binary contraction method. The remaining data points —numbering approximately 8000— were used for verification. Fig. 9 shows the training and verification costs as a function of the epochs in the training. During training, both costs decrease monotonically at varying speeds and in concert with each other for most of the epochs. After more than about $10^{5}$ epochs, the verification cost ceases to decrease while the training cost decays further. It is at this point (when additional training ceases to improve the verification) that the training of a NN should cease. Additional training simply provides better agreement with the training data while possibly distorting the general accuracy of the NN. In this example, we see that significant accuracy of the NHIM can be obtained (down to a cost of $10^{-6}$, measured in units of dimensionless length squared as described in Sec. 3.3), with a relatively small NN.

4.2 Gaussian process regression

Gaussian process regression (GPR) offers an alternative to the NN as a supervised machine learning approach for fitting a multidimensional hypersurface such as the NHIM to isolated training points. This approach interpolates given values of a function $f(\boldsymbol{z})$ that depends on an arbitrary number of input variables $\boldsymbol{z}$. In our application, the input variables are the bath coordinates $\boldsymbol{y}$, velocities $\boldsymbol{v_{y}}$ and time $t$. The unknown functions ${x}^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ and ${v_{x}}^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ will be interpolated independently.

GPR is based on the assumption that the values of a function $f(\boldsymbol{z})$ are Gaussian random variables that are a priori unknown, apart from their statistical properties. After some values of the function have been learned, other values are inferred from the given information. For this inference to be possible, some a priori knowledge is required about the length scales along which the values of the function vary. Specifically, we assume that we know the a priori mean $\mu(\boldsymbol{z})$ of the function $f(\boldsymbol{z})$ and the covariance $k(\boldsymbol{z}_{1},\boldsymbol{z}_{2})$ between the values $f(\boldsymbol{z}_{1})$ and $f(\boldsymbol{z}_{2})$. In an unbiased a priori distribution these means would be set to zero, i.e., $\mu(\boldsymbol{z})=0$. A typical choice for $k(\boldsymbol{z}_{1},\boldsymbol{z}_{2})$ is the squared exponential kernel

[TABLE]

where the hyper-parameter $l$ gives the typical length scale along which the values of $f$ are expected to vary. It yields high covariance between closely located points $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ and a lower covariance between distant points. In an application to TST, the characteristic length $l$ could differ for time and space coordinates, or even for different bath modes.

Assume now that the values $\boldsymbol{f}_{\text{train}}={\left(f(\boldsymbol{z}^{(1)}_{\text{train}}),\dots,f(\boldsymbol{z}^{(n)}_{\text{train}})\right)}^{\text{T}}$ of the function $f$ are known in $n$ points $\boldsymbol{z}^{(i)}_{\text{train}}$ and the values $\boldsymbol{f}_{\text{test}}={\left(f(\boldsymbol{z}^{(1)}_{\text{test}}),\dots,f(\boldsymbol{z}^{(m)}_{\text{test}})\right)}^{\text{T}}$ in $m$ other points $\boldsymbol{z}^{(j)}_{\text{test}}$ are sought. Initially, the values $\boldsymbol{f}_{\text{train}}$ and $\boldsymbol{f}_{\text{test}}$ follow a multidimensional Gaussian distribution with zero mean and covariance matrix

[TABLE]

that consists of the blocks

[TABLE]

The a posteriori conditional distribution of the unknown values $\boldsymbol{f}_{\text{test}}$, given the known values $\boldsymbol{f}_{\text{train}}$, is a Gaussian function with mean 49, 50

[TABLE]

and covariance matrix

[TABLE]

The components of the conditional mean $\bar{\boldsymbol{\mu}}$ are taken to be the interpolated values of the function. If desired, the diagonal elements of the conditional covariance matrix $\bar{\Sigma}$ can be used to provide an error estimate for the interpolation. If the standard deviation $\sqrt{\bar{\Sigma}_{ii}}$ is sufficiently small, the interpolated value $f(\boldsymbol{z}^{(i)}_{\text{test}})$ is known essentially with certainty.

In practice, we compute a Cholesky decomposition $K=LL^{\text{T}}$ of the matrix $K$, where $L$ is a lower triangular matrix. It is then easy to pre-compute the vector $K^{-1}\boldsymbol{f}_{\text{train}}=L^{-\text{T}}L^{-1}\boldsymbol{f}_{\text{train}}$ that depends only on the training points, but not on the points $\boldsymbol{z}^{(i)}_{\text{test}}$. After that, the computation of one interpolated function value according to Eq. (14) requires the computation of the corresponding row of the matrix $N$ and the scalar product with the precomputed vector. The effort required is proportional to the number $n$ of training points. The conditional covariance matrix (15) can be written as

[TABLE]

where the matrix $\tilde{N}$ is the solution of the equation $N=\tilde{N}L^{\text{T}}$. Because $L$ is triangular, the computation of one row of $\tilde{N}$, and therefore of an error estimate for one data point, requires an effort of order $n^{2}$.

The training and use of a GPR machine can be illustrated relative to a sine function shown in Fig. 10 (b). We use the squared exponential kernel (9) with a characteristic length scale $l=1$. The conditional mean that serves as the interpolated function is indicated by the dashed red line, a 95 % confidence interval based on the conditional standard deviation by the gray area.

The difference between the GPR machine’s approximation and the original sine-function is small near the training points and increases elsewhere. Whereas for the prior, the mean $\mu$ and the $95\,\%$ confidence region are independent of $z$, the learning of several training points alters both. At the training points, this $95\,\%$ region is zero, as these points are assumed to be correct. The further away from training points, the more the confidence region grows. This increase happens on a length-scale compared to the characteristic length $l=1$ of the assumed squared exponential kernel (9). At points far from where no training points are given, e. g. for $z<-3$, the confidence region of the posterior realigns with the Gaussian prior. As a result, the accuracy of the approximation increases with the number of training points.

To illustrate that the posterior of the GPR machine gives a distribution over admissible functions, we plot three functions randomly sampled from the posterior distribution in Fig. 10 (b) and compare to three samples of the prior distribution in Fig. 10 (a). The samples of the posterior distribution match the training points and are mostly located within the altered $95\,\%$ confidence region around the approximation $\mu(z)$ to the original function $f(z)$. Therefore, the posterior ‘chooses’ which functions of the prior are suitable to interpolate between the given training points.

4.3 NNs vs. GPR

We now compare the advantages and disadvantages of the two machine learning methods introduced in Secs. 4.1 and 4.2 for the interpolation of the NHIM and construction of the DS in the calculation of rate constants. The aim is to provide a guideline on when to use which method.

For the training of a NN on a given set of pre-calculated points of the NHIM, little previous knowledge on the interpolated object is needed. For example, in a typical NN that approximates ${(x,v_{x})}^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$, the number of neurons in each of the input and output layers are fixed as illustrated in Fig. 8. The only information that is prescribed is the structure of the hidden layers, and it merely needs to contain as much complexity as the system that it is re-generating. Its efficacy is confirmed if the optimization of the neurons leads to weights and biases that are sufficiently small compared to the target accuracy with respect to the cost function.

When using GPR, the success and the speed of training highly depends on the chosen hyper-parameters of the kernel function. In choosing these parameters, one can include previous knowledge of the object to be interpolated, e. g. the expectation that the NHIM of a given system may be a smooth surface with small curvature. Consequently, the characteristic length in the direction of the bath coordinates can be higher than the characteristic length of the variation in time, if the external driving leads to a fast moving NHIM whose curvature does not vary heavily. If such previous knowledge is available, training can be accelerated significantly by choosing an appropriate kernel function or by a corresponding restriction of the parameter ranges in the optimization.

Another difference between both methods lies in the procedure of training itself. For NNs, training with the back propagation algorithm is an iterative process. It iterates a possibly high number of times —that is, epochs— over the given training set gradually adjusting weights and biases to achieve a certain level of convergence in the training and verification cost. This procedure typically takes a lot of computation time. Using GPR, training can be effectively broken down to basic routines in linear algebra without need for iteration. Therefore, training to a certain accuracy in the context of GPR is usually much faster compared to the use of NNs.

We compare the training of the neural network already discussed in Sec. 4.1 to the training done by GPR on the NHIM of the same two-dimensional model potential according to Eq. (7). The NN is set up with two hidden layers with $40$ and $10$ neurons as before. For GPR, we use the squared exponential kernel (Eq. (9)) with $l=2.0$ as hyper-parameter. In the left panel of Fig. 11, the drop of the cost function in relation to the needed computational time is shown for the training for each of the NN and GPR machines. Whereas training a NN converges after several thousand seconds, the GPR approach reaches the same level of accuracy within seconds.

Another difference in the implementation of the two methods lies in the relative time each machine learning approach requires to evaluate the function ${(x,v_{x})}^{\mathrm{NHIM}}(\boldsymbol{y},\boldsymbol{v_{y}},t)$ for some input of the bath coordinates. For simplicity, we report the internal run time reported by the same computer and call it the function call time. Precision for GPR is achieved by adding more training points and therefore decreasing the $95\,\%$ region of the posterior. But adding more training points leads to an increasing dimension of the matrices discussed in Sec. 4.2. Consequently, the function call time required to evaluate Eq. (14) increases linearly with the number of training points as shown by the red curve in the right panel of Fig. 11. In contrast, increasing the number of training points does not increase the evaluation time for the NN, assuming its layer geometry does not change. Thus evaluating a small NN can be significantly faster than GPR, although it might have taken a long time to train the net.

Based on these three observations, we now have guidelines for choosing between the NN and GPR machines for continuously interpolating the time-dependent NHIM or the related DS: (i) In cases when the function call time exceeds the additional training time of the neural network compared to GPR, it is computationally more efficient to use the neural network approach. A simple example would be a system with time-periodic external driving (and without friction). If the focus is on calculating rates for different initial ensembles of particles under the influence of the same external force, the NHIM and therefore the DS is always the same. Consequently, it is more efficient to use NNs so as to benefit from the short function call time. (ii) On the other hand, in cases when the same ensemble is used to calculate reaction rates for a system under influence of a different external driving or in cases with non-periodic driving, the training time for each NN can easily be higher than the overall function call time to propagate the ensemble and calculating reaction rates. It is then more efficient to use GPR. Another example is the future application of the methods presented here on thermal systems which experience nontrivial friction and fluctuations. In these systems, there will be a different NHIM (and DS) for each sequence of random numbers modeling the fluctuations in a Langevin-type approach, and each sequence will require separate training. The training time can then easily exceed the function call time. (iii) A third guideline can be formulated in terms of how many training points are needed for a given method. Using NNs, millions of training points can be fit easily, whereas with GPR, the function call time is rather high even for just thousands of training points, and its use to fit millions of training points is hopeless. This comparison can also be seen exactly in the opposite way, as prior knowledge on the object to be interpolated—here the NHIM—can be included in the GPR method by choosing an appropriate kernel function. This usually reduces the amount of training data needed significantly and provides an advantage of GPR over the NN approach. If only a few training points are available (for example if training data is hard to obtain), GPR is still capable of interpolating these few points when provided a suitable kernel function, at least one with respect to optimization of the cost function.

5 Application to a two-dimensional system

In this section, we assess the performance of the methods summarized above for the central application of this work. That is, the determination of the time-dependent NHIM and the associated rate for the driven model potential (7) for the same parameters as used in Sec. 3.3. In the previous section, we focused on the accuracy as determined by the cost function as one would typically use in training a NN or GPR machine. However, in our specific application, we also have an external criterion at our disposal that allows us to assess the accuracy of the NHIM. Namely, the DS attached to the NHIM should by construction be recrossing-free. For a numerical test of this criterion, a large ensemble of trajectories can be launched in the neighborhood of the barrier and propagated forward in time until it leaves this neighborhood. In the process, the number of crossings and recrossings of the DS can be monitored. Beyond the decrease in the cost function, the demand of the DS to be as recrossing free as possible is another—and in our case more relevant—indicator of a successful training process.

To obtain a stringent test of the DS, it is advantageous to launch trajectories close to the stable manifolds. These trajectories will remain in the neighborhood of the barrier for a long time and therefore have plenty of opportunity to recross a DS of low quality. A detailed description of such tests for the NHIM constructed by other approaches, together with numerical results, can be found in Refs. 44, 45 for various test systems. In what follows, we test the NN and GPR machines relative to this criteria.

5.1 Representation of the NHIM by NNs or GPR

For both NNs and GPR, training points are computed in the range $y\in[-4,4]$ and $v_{y}\in[-8,8]$. The NN is trained on an equidistant grid of $25$ points in $y$, $40$ points in $v_{y}$ range and $40$ points in $t$, sampled over one period, i.e., in the interval $[0,2]$. The structure of the net uses three input neurons for the variables $y$, $v_{y}$, and $t$, four hidden layers with $40$, $40$, $40$, and $10$ neurons, respectively, and one output neuron representing $x^{\mathrm{NHIM}}$. While the NN discussed in Sec. 4.3 was written from scratch in our group, here we resort to an efficient implementation using TensorFlow.96 To improve the periodicity of the resulting function $x^{\mathrm{NHIM}}$, the grid of training points is tripled to cover a full time span from $t\in[0,6]$. All in all, we use $120\,000$ training points taking advantage of the periodicity of the driven system, so just $40\,000$ points have to be obtained using the binary contraction method explained in Sec. 3.3. After training for $50\,000$ epochs, the training cost has decreased to $2.06\times 10^{-7}$ and the verification cost to $2.15\times 10^{-7}$. Although using so many training points might seem to be a bit conservative here, we did so as to ensure, that the likelihood of single, localized artifacts in the representation of the DS is reduced. Such artifacts, which are hidden by averaging over a large amount of accurately inferred data points in the loss function, will cause spurious recrossings to the trajectories. This would lead to errors in the corresponding rates which goes against our primary reason for calculating the non-recrossing DS. Thus the large number of points—though not fully optimized—was chosen in order to ensure not only the proper construction of the manifold but also to satisfy our more stringent criterion for a reduced number of recrossings.

For the GPR representation of the NHIM, we used $6000$ training points in total on equidistant grids chosen at $30$ time steps in the range $t\in[-0.5,2.4]$ and in the same $y$ and $v_{y}$ ranges as for the NN. The time range is more than the full period of $T=2$ of the time-dependent potential (7) in order to improve periodicity. Training was done using a squared exponential kernel like in Eq. (9), but with different characteristic length scales $l$ in any direction $(y,v_{y},t)$. Training was carried out for various values of the hyperparameters on a logarithmic grid, and the set that led to the smallest verification cost was retained. The final values of the hyperparameters that best reflect the length and time scales on which the NHIM moves, were $\sigma_{y}\approx 0.705$, $\sigma_{v_{y}}\approx 0.790$ and $\sigma_{t}\approx 0.505$. Note, that the equidistant grids are chosen in order to provide no initial bias, and thus focus on whether the different approaches are feasible. Choosing non-equidistant grids may be a subject for future optimization of the training process.

To assess the accuracy of the representation obtained from either the NN or the GPR, we computed the position of the NHIM at a fixed time for a grid of $50\times 50$ additional points in the ranges $y\in[-4,4]$ and $v_{y}\in[-8,8]$ with an accuracy of $10^{-12}$ using the binary contraction method. Afterwards, the values $x_{\mathrm{NN}}(y,v_{y},0)$ and $x_{\mathrm{GPR}}(y,v_{y},0)$ were obtained for both trained machines of the NHIM at all of these $50\times 50$ verification points. Fig. 12 shows the location of the NHIM at time $t=0$ (in a three-dimensional projection of the four-dimensional phase space). Very similar figures have been obtained for different times. The color of the surfaces represents the error of the machine learning representation of the NHIM compared to the numerically exact result.

The maximum error of the NN representation is $6.1\times 10^{-3}$ as seen in Fig. 12(a). It is nearly a factor of three smaller than the maximum error for GPR in Fig. 12(b). Large errors of the NN occur just in a limited region near $y=v_{y}=0$, whereas the GPR has large errors in a band-like structure near $v_{y}=0$. The maximum difference could be lowered by e. g. using more training points. But this will come with the burden of longer training times for neural networks or longer function call times for GPR. A different approach, to be pursued in the future, might be the use of a non-equidistant grid of training points. By increasing the density of points in the regions where the quality of the representation is poor and decreasing it elsewhere, it should be possible to lower the maximum error while keeping the total number of training points constant.

As noted above, the GPR interpolates the given training points exactly, whereas the NN attempts to minimize the training cost, i.e., the mean square error, averaged over all training points, without bringing it to zero. It now appears that the NN produces a smoother representation of the NHIM, with moderate errors everywhere, whereas the GPR yields very small errors near the training points at the expense of larger errors elsewhere.

Note that a NN typically has intrinsic hyper-parameters like the activation function or the learning rate. Although adjusting these parameters is not free, we neglect this adjustment time in our discussion, since in our experience, these parameters have to be adjusted only once for the general problem, but can be kept constant while just varying potential parameters as, e. g., the amplitude of the moving saddle.

Finally, we comment on the numerical effort required by both methods, as represented by the computing times listed in Table 2. These times depend critically on, e. g., the hardware used or the desired accuracy, and should be understood as a rough guidance only. Computing all $40\,000$ training points for the neural network took $3636\,\mathrm{s}$, whereas obtaining $6\,000$ points for the GPR took only $545\,\mathrm{s}$. (This corresponds to $0.091\,\mathrm{s}$ on average per training point.) The process of training the neural net for $50\,000$ epochs took $4220\,\mathrm{s}$. GPR, by contrast, could be trained in only $9\,\mathrm{s}$. However, this training had to be repeated many times in order to optimize the hyperparameters. In total, the process took $2238\,\mathrm{s}$. In summary, the neural network was was ready to use after $7856\,\mathrm{s}$ (about $131\,\mathrm{min}$) the GPR about $2792\,\mathrm{s}$ (about $47\,\mathrm{min}$). Afterwards, the evaluation of a single value of $x^{\textrm{NHIM}}$ took $14\,\upmu\mathrm{s}$ for the neural network and $253\,\upmu\mathrm{s}$ for GPR. Thus, although training in GPR is much faster, its function call time is by a factor of $18$ higher. As we concluded at the end of Sec. 4, depending on how many function calls are required after training, either method can be more computationally efficient: The NN becomes advantageous if more than about 21 million evaluations are made.

5.2 Relative accuracy of the rate constants:

NNs vs. GPR

To obtain a rate constant for a reaction over the time-dependent barrier of potential (7), we propagate an ensemble of particles initialized on the reactant side and monitor the reaction of particles to the product side. The ensemble is set up with $N_{\mathrm{react}}(t=0)=10^{6}$ particles at $x_{0}=-0.2$ with $y$ on the minimum energy path of potential (7) at time $t_{0}=0$. The velocities $(v_{x},v_{y})$ of these particles are sampled from a thermal distribution with $k_{\mathrm{B}}\,T=0.5$.

Propagating the full ensemble using a standard velocity-Verlet integrator, we monitor possible crossings of the DS and the related crossing times of each particle. Ideally, each trajectory should either descend on the reaction side of the barrier without crossing the DS or cross it at exactly once and then descend on the product side. In practice, because of inaccuracies in how the DS is represented, we expect a small number of trajectories to recross the DS.

A histogram with the number of crossings is shown in Fig. 13 (a) for the NN and in Fig. 13 (b) for the GPR approach. In both cases approx. $16\,\%$ of all particles are reactive and cross the DS once from the reactant to the product side. Nearly all other particles are non-reactive. The number of trajectories that cross the DS more than once is very small in both cases, though the GPR still produces a factor $22$ more recrossings than the NN. This difference arises because the NN represents the NHIM more accurately than GPR as shown in Fig. 12.

Rates are obtained in dimensionless units using the procedure described in Sec. 2.3. For the NN approach in Fig. 13 (c), we obtain a rate constant of $k_{\mathrm{NN}}=3.1$ compared to the rate constant of $k_{\mathrm{GPR}}=3.0$ for the GPR approach in Fig. 13 (d). As discussed in Secs. 2.3 and 3.3, these rate constants are given in units of a dimensionless inverse time. The time-scale is given by a period with $T=2$ of the oscillation with frequency $\omega_{x}=\pi$ in Eq. (7). Although these constants are very similar, $k_{\mathrm{NN}}$ is presumed to be more precise because the DS of the NN approach exhibits fewer recrossings. Note that differences between the reactant populations computed by the two methods arise once the populations have decayed to approximately the number of recrossing trajectories in GPR.

6 Conclusion and outlook

In this paper, we have presented several techniques to obtain the underlying geometric structures needed to obtain rate constants for reactions over time-dependent driven barriers in the context of TST. Specifically, we focused on the determination of the NHIM which is generally a multidimensional, time-dependent object located in phase space, to which a recrossing-free DS can be attached, and with respect to which the TST rate is obtained. As we restrict ourselves to rank-1 saddles, the procedures presented here can in principle be carried out for any multidimensional system with one unstable degree of freedom and an arbitrary number of stable bath degrees of freedom. Here, one needs to keep in mind that in higher-dimensional systems not only does the propagation of trajectories become more complex and computationally more challenging, but also the NHIM itself scales modestly with the dimensionality of the equations of motion. Depending on how one chooses to construct the DS, the number of points needed to characterize the underlying NHIM may scale exponentially with increasing dimensionality.

A central result of this work is the extension of recent approaches to address reactions of increasing dimensionality. We elaborate on the use of our recently developed binary contraction method46 for the accurate non-perturbative determination of points on the NHIM. In the case of time-periodically driven Hamiltonian systems, this allowed us to extend the one-dimensional TS trajectory to multi dimensions by defining it to be the unique periodic orbit that remains in the vicinity of the saddle.

A second central result of this work is the elaboration and extension of our recent work on NNs.45 Here we show the applicability of the NN to construct the NHIM and rates for a two-dimensional system. We also demonstrated that a different machine learning approach, the GPR, is also effective. These two approaches have competing efficiencies in training and use. Together, they offer a powerful combination for obtaining high-accuracy smooth NHIMs for barriers with rank-1 saddles in higher dimensions. Indeed, Kamath et al97 recently compared both of these methods for the first time in the context of constructing potential energy surfaces.92 Although their work has a different objective than the NHIMs obtained here, we likewise obtained similar conclusions regarding the advantages and disadvantages of the two methods with respect to the determination of the target surface. Carrington and coworkers92 employed a secondary criterion for evaluating the quality of their computed surface. Namely, the GPR outperforms NN when fitting potential energy surfaces to obtain vibrational spectra. In our case, we see in Fig. 13 that the total number of recrossings for the DS obtained by GPR is significantly higher than that found by the NN. Following this additional —and in our case more relevant—requirement to the DS that it be recrossing-free, we conclude that for our applications the NNs are more useful in the sense that they provide higher accuracy for an indirect metric not used as part of the cost function while training. We emphasize that we do not claim our setup to be fully optimized in terms of net structure or the training process, but it is the best we found using comparable approaches.

For the determination of the rate, the computation of the required reaction coordinates of the DS from scratch (using e. g. the binary contraction method) whenever it is needed would be prohibitively expensive, as the position of the DS must be compared to the instantaneous position of each propagated particle. Having a continuous representation of the DS as a function of bath modes and time reduces this numerical effort drastically. This is where the use of NN and GPR machine-learning techniques represent a key advance by way of providing a suitable smooth surfaces.

It remains to show that these approaches are effective for the determination of rate constants in systems with more complex energy landscapes or external potentials. This includes aperiodic external driving, potentials with multiple saddles, dissipative terms and thermal noise. In addition, although these methods are developed for in principle arbitrary dimensional problems, everything done so far was an application to two- and three-dimensional model systems 44, 45. So one aspect for a future work will be to test these methods on really high dimensional systems with tens of degrees of freedom to answer the question if for increasing dimensionality really an exponential increase in the obtained number of points on the NHIM is necessary or if these points can be chosen more efficiently by improving the machine learning procedures.

An application for future work will lead away from model systems towards real systems like the LiCN $\leftrightarrow$ LiNC isomerization reaction, to which an analytical approximation to the molecular potential 98 as well as to the molecular dipole moment 99 is known. Another possible application includes the ketene isomerization100, 67, 101 which has received increased attention since the initial work by some of us35 using the Lagrangian descriptor to obtain the NHIM and the associated rate constants. The determination of the full-dimensional NHIM in both of these cases remains a challenge that may be addressable using the emerging approaches presented here.

{acknowledgement}

The German portion of this collaborative work was partially supported by Deutsche Forschungsgemeinschaft (DFG) through Grant No. MA1639/14-1. The US portion was partially supported by the National Science Foundation (NSF) through Grant No. CHE 1700749. AJ acknowledges the Alexander von Humboldt Foundation, Germany, for support through a Feodor Lynen Fellowship. MF is grateful for support from the Landesgraduiertenförderung of the Land Baden-Württemberg. This collaboration has also benefited from support by the European Union’s Horizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie Grant Agreement No. 734557. Surface plots in Figs. 6 and 7 have been made with the Mayavi software package 102.