It is always useful to start with the basics to make the ideas concrete. So, the following one DoF system is the simplest model of a chemical reaction over a potential energy barrier between reactants and products. This is the archetypical one DoF Hamiltonian system modelling reaction dynamics associated with a saddle point. The only configuration space coordinate, q, is the reaction coordinate. The configuration space is one dimensional and the phase space is two dimensional.
The Hamiltonian for a (linear) one DoF saddle point is given by:
H(q,p)=λ2(p2−q2)=λ2p2⏟kineticenergy−λ2q2⏟potentialenergy,λ>0.In this simple system ''reaction'' corresponds to trajectories that change sign in q, which requires H>0 (as shown in Fig. fig:1 b)) and non-reacting trajectories have H<0.
The associated Hamilton's equations of motion are:
˙q=∂H∂p=λp,˙p=−∂H∂q=λq.In Fig. fig:1 a) we show the potential energy, V(q)=−λ2q2. In Fig. fig:1 b) we show the contour lines of the Hamiltonian for some values of the total energy, typically referred to as the phase portrait corresponding to (1).
Now we discuss the phase space structures -- NHIM, its stable and unstable manifolds -- and their role in constructing the DS. All of these notions are trivial in this simple setting, but they will serve to focus the ideas when we consider more DoF.
For this case the NHIM is the saddle point at the origin (a single point is a trivial example of a manifold). It only exists on the H=0 energy surface (this is very different when we go to two, and more, DoF) and its stable and unstable manifolds are the diagonal lines (also on the H=0 energy surface--the stable and unstable manifolds of a NHIM have the same energy as the NHIM).
The non-isoenergetic DS can be taken as the line q=0. Clearly, it has the ''no-recrossing'' properties and all reacting trajectories must cross this line. The DS at a fixed (positive) energy is given by
λ2p2=H=constant,or
p=±√2λH.So for a fixed energy H > 0 the DS consists of two distinct points: p=+√2λH (the dividing surface for forward reactions) and p=−√2λH (the dividing surface for backward reactions). These points are just the intersections of the reacting trajectories with q=0.
The following two DoF system is the one DoF system (reactant) coupled to a harmonic oscillator (bath). The configuration space coordinates are q1 and q2. The configuration space is two dimensional and the phase space is four dimensional.
The system is described by a quadratic 2 DoF Hamiltonian:
H(q1,q2,p1,p2)=λ2(p21−q21)⏟H1+ω2(p22+q22)⏟H2,λ,ω>0,where H1 represents the energy in the reactive DoF and H2 corresponds to the energy associated with the bath DoF. The corresponding Hamilton's equations are given by:
˙q1=∂H∂p1=λp1,˙p1=−∂H∂q1=λq1,˙q2=∂H∂p2=ωp2,˙p2=−∂H∂q2=−ωq2,These equations have an equilibrium point of saddle-center equilibrium type (index one saddle) at the origin. In Fig. fig:2 a) we show contours of the potential energy and in Fig. fig:2 b) we show the phase portrait corresponding to (5). Since the Hamiltonians H1 and H2 are uncoupled we can sketch the phase portraits for each separately, and discuss the distribution of total energy between each DoF in a simple manner.
Notice that this system can be easily solved. Given the initial condition
x0=x(t0)=(q01,q02,p01,p02),t0=0the general solution to Eq. (6) has the form:
q1(t)=q01cosh(λt)+p01sinh(λt),p1(t)=q01λsinh(λt)+p01λcosh(λt)q2(t)=q02cos(ωt)+p02sin(ωt),p2(t)=−q02ωsin(ωt)+p02ωcos(ωt)If the energy of the bath DoF is:
H2=ω2((p02)2+(p02)2)then the trajectory in the q2-p2 plane satisfies the equation:
2H2ω=(q02)2+(p02)2giving rise to a periodic orbit (a circle) in the q2-p2 plane with radius √2H2/ω.
Note that trajectories corresponding to H1 can become unbounded and trajectories corresponding to H2 are bounded. Hence in this system reaction occurs when the q1 coordinate of a trajectory changes sign. Therefore, a ''natural'' dividing surface would be q1=0. This is a three dimensional surface in the four dimensional phase space. We want to examine it's structure more closely and, in particular, its' intersection with a fixed energy surface. We will also utilize terminology from chemistry by referring to H1 as the ''reactive mode'' and H2 is the ''bath mode''.
First, note that for reaction to occur we must have H1>0, since the q1 component of reacting trajectories changes sign. Also, it is clear from the form of H2 that H2≥0. Therefore, for reaction we must have H=H1+H2>0. The energy surface is given by:
S(H)={(q1,q2,p1,p2)|H=H1+H2=λ2(p21−q21)+ω2(p22+q22)}The intersection of q1=0 with this energy surface is given by:
D(H)=S(H)∩{q1=0}={(q1,q2,p1,p2)|H=λ2p21+ω2(p22+q22),q1=0}This is the isoenergetic DS. It has the form of a 2-sphere in the four dimensional (q1,p1,q2,p2) space. It has two ''halves'' corresponding to the forward and backward reactions, respectively:
Df(H)={(q1,q2,p1,p2)|p1=√2λ√H+λ2q21−ω2(q22+p22),q1=0}and backward reaction:
Db(H)={(q1,q2,p1,p2)|p1=−√2λ√H+λ2q21−ω2(q22+p22),q1=0}Observe that the DS has the no-recrossing property, meaning that the Hamiltonian vector field is never zero or tangent to the DS. Moreover, since Hamilton's equations impose that ˙q1=λp1, by construction any point starting on the forward DS will leave the DS in the same direction, and the same applies for the backward DS. The two hemispheres of the DS meet at the NHIM, which is the equator of the DS and is described by the set:
N(H)=S(H)∩{q1=p1=0}={(q1,q2,p1,p2)|H=ω2(p22+q22),q1=p1=0}which is a circle (S1). Furthermore, it is an invariant manifold because the phase space trajectory of an initial point on the NHIM, which has q1=p1=0, due to Hamilton's equations ˙q1=λp1 and ˙p1=λq1 that imply ˙q1=˙p1=0, remains on the NHIM forever. The NHIM is normally hyperbolic since it has saddle-type stability as the q1-p1 coordinates transverse to it have exponential growth and decay dynamics. This can be explicitly seen from Eq. (8). This behavior allows for the construction of the stable and unstable invariant manifolds of the NHIM:
WuN(H)={(q1,q2,p1,p2)|H=ω2(p22+q22),p1=q1}which are known in the dynamical systems literature as spherical cylinders because they have the topological structure R×S1, that is, the cartesian product of a line (p1=q1 or p1=−q1) and a one-dimensional sphere. It is straightforward to show that they are invariant and that any phase space point starting on the unstable (stable) manifold will approach the NHIM as t→−∞ (t→+∞). Moreover, if an initial condition is inside the spherical cylinder correspinitialonding to the stable (unstable) manifold, then it will cross the phase space bottleneck that opens in the neighborhood of the index-1 saddle at the origin in forward (backward) time. Consequently, the trajectories of these initial conditions can be classified as reactive trajectories that go from the reactants to the products region. For this reason, the stable and unstable manifolds of the NHIM are known in the Chemistry literature as reactive cylinders, giving rise to reactive islands when they intersect a Poincare section of the phase space. On the other hand, if an initial condition is outside the spherical cylinders, the trajectory will not cross from reactants to products, and that trajectory would be classified as non-reactive. For a visual illustration of all these concept see Figs. fig:3 , fig:4 and fig:5.
The following three DoF system is the one DoF system (reactant) coupled to two harmonic oscillators (bath modes). The configuration space coordinates are q1, q2 and q3. The configuration space is three dimensional and the phase space is six dimensional.
We consider a quadratic 3 DoF Hamiltonian:
H=λ2(p21−q21)⏟H1+ω22(p22−q22)⏟H2+ω32(p23−q23)⏟H3,λ,ω2,ω3>0with the corresponding Hamilton's equations given by:
˙q1=∂H∂p1=λp1,˙p1=−∂H∂q1=λq1,˙q2=∂H∂p2=ω2p2,˙p2=−∂H∂q2=−ω2q2,˙q3=∂H∂p3=ω3p3,˙p3=−∂H∂q3=−ω3q3,These equations have an equilibrium point of saddle-center-center equilibrium type (index one saddle) at the origin. Since the Hamiltonians H1, H2 and H3 are uncoupled we can analyze the phase portraits for each separately. As in the previous examples, H1 corresponds to the ''reactive mode'' (trajectories can become unbounded) and H2 and H3 are ''bath modes'' (trajectories are bounded).
In this system reaction occurs when the q1 coordinate of a trajectory changes sign. Hence, as in the 2 DoF example, a ''natural'' dividing surface would be q1=0. This is a five dimensional surface in the six dimensional phase space. We want to examine its' structure more closely and, in particular, its intersection with a fixed 5 dimensional energy surface.
First, note that for reaction to occur we must have H1>0. Also, it is clear that H2,H3≥0. Therefore, for reaction we must have H=H1+H2+H3>0. The energy surface is given by:
λ2(p21−q21)+ω22(p22+q22)+ω32(p23+q23)=H1+H2+H3=H>0,H1>0,H2,H3≥0.The intersection of q1=0 with this energy surface is given by:
λ2p21+ω22(p22+q22)+ω32(p23+q23)=H1+H2+H3=H>0,H1>0,H2,H3≥0.This is the isoenergetic DS. It has the form of a 3-sphere in the four dimensional (q1,p1,q2,p2,q3,p3) space. It has two ''halves'' corresponding to the forward and backward reactions, respectively:
p1=+√2λ√H1+H2+H3−ω2(p22+q22)−ω2(p23+q23),forward DS,The forward and backward DS ''meet'' at p1=0:
ω22(p22+q22)+ω32(p23+q23)=H2+H3≥0,NHIM,which is a normally hyperbolic invariant 3 sphere. It is {\em invariant} because on this set q1=p1=0 and, from (19), if q1=p1=0 the ˙q1=˙p1=0. Hence, q1 and p1 always remain zero, and therefore trajectories with these initial conditions always remain on (24). In other words, it is invariant. It is normally hyperbolic for the same reasons as for our 2 DoF example. The directions normal to (24), i.e. q1−p1, are linearized saddle like dynamics.
The stable and unstable manifolds of the NHIM are given by:
Wu(M(h))={(q1,p1,q2,p2,q3,p3)|q1=p1,ω22(p22+q22)+ω32(p23+q23)=H2+H3>0},Ws(M(h))={(q1,p1,q2,p2,q3,p3)|q1=−p1,ω22(p22+q22)+ω32(p23+q23)=H2+H3>0},and they are four dimensional on a fixed five dimensional energy surface. The topology of the manifolds is the product of a line (q1=p1 or q1=−p1) with a 3 sphere (ω22(p22+q22)+ω32(p23+q23)=H2+H3>0), that is R×S3 and sometimes such geometry is referred to as spherical cylinders. Since the lines q1=p1 and q1=−p1 correspond to the contour H1=0, in the six dimensional phase space, these manifolds have energy H=H1+H2+H3=0+H2+H3>0.
We consider isoenergetic two-dimensional surfaces parametrized by two coordinates and compute the Lagrangian descriptor in a square domain of size 2 units around the origin. We discretize the coordinates of the two dimensional surface and pick constant values for three of the four remaining coordinates, and use the total energy equation to solve for the sixth coordinate. Due to the form of the Hamiltonian (18), obtaining the coordinate from the constant energy condition is simply solving a quadratic equation.
- Isoenergetic two-dimensional surface parametrized by (q1,p1) On the constant energy surface, H(q1,p1,q2,p2,q3,p3)=h, we compute Lagrangian descriptor on a two-dimensional surface parametrized by (q1,p1) coordinates by defining
where
p3(q1,p1,q2=0,p2=0,q3=0;h)=√2ω3(h−λ2(p21−q21))The intersection of the two-dimensional surface U+q1p1 with the NHIM (24) becomes
M(h)∩U+q1p1={(q1,p1,q2,p2,q3,p3)|q1=0,p1=0,q2=0,p2=0,q3=0,˙q3>0:p3(q1,p1,q2,p2,q3;h)>0}.Thus, the NHIM is located at the origin (0,0) and marked by a red cross in the LD plot (Fig. fig:6).
Next, the intersection of the two-dimensional surface with the unstable (25) and stable manifolds (26) is given by
Wu(M(h))∩U+q1p1={(q1,p1,q2,p2,q3,p3)|q1=p1,q2=0,p2=0,q3=0,˙q3>0:p3(q1,p1,q2,p2,q3;h)>0},Ws(M(h))∩U+q1p1={(q1,p1,q2,p2,q3,p3)|q1=−p1,q2=0,p2=0,q3=0,˙q3>0:p3(q1,p1,q2,p2,q3;h)>0},which are one-dimensional for a fixed energy, and represent lines passing through the origin shown as dashed red (unstable) and white (stable) lines, respectively, in Fig. fig:6. The only points of local minima in the LD plot (Fig. fig:6) also lie along the lines passing through the origin and correspond to the manifolds of the NHIM.
- Isoenergetic two-dimensional surface parametrized by (q2,p2) On the constant energy surface, H(q1,p1,q2,p2,q3,p3)=h, we compute Lagrangian descriptor on a two-dimensional surface parametrized by (q2,p2) coordinates by defining
where
p1(q1=0,q2,p2,q3=0,p3=0;h)=√2λ(h−ω22(p22+q22))The intersection of the two-dimensional surface U+q2p2 with the NHIM (24) becomes
M(h)∩U+q2p2={(q1,p1,q2,p2,q3,p3)|q1=0,p1=0,p3=0,q3=0,˙q3>0:ω22(p22+q22)=h}.Thus, the NHIM is the circle of radius √2h/ω2 and marked by a dashed line in the LD plot (Fig. fig:7).
Next, the intersection of the two-dimensional surface with the unstable (25) and stable manifolds (26) is given by
Wu(M(h))∩U+q2p2={(q1,p1,q2,p2,q3,p3)|q1=p1,q1=0,q3=0,p3=0,˙q3>0:ω22(p22+q22)=h},Ws(M(h))∩U+q2p2={(q1,p1,q2,p2,q3,p3)|q1=−p1,q1=0,q3=0,p3=0,˙q3>0:ω22(p22+q22)=h},which are one-dimensional for a fixed energy, and marked by dashed red (unstable) and white (stable) lines, respectively, in Fig. fig:7. The only points of local minima in the LD plot (Fig. fig:7) also lie along the lines passing through the origin and correspond to the manifolds of the NHIM.
- Isoenergetic two-dimensional surface parametrized by (q3,p3)
On the constant energy surface, H(q1,p1,q2,p2,q3,p3)=h, we compute Lagrangian descriptor on a two-dimensional surface parametrized by (q3,p3) coordinates by defining
U+q3p3={(q1,p1,q2,p2,q3,p3)|p1=0,q2=0,p2=0,˙p1>0:q1(p1,q2,p2,q3,p3;h)>0}where
q1(p1=0,q2=0,p2=0,q3,p3;h)=√2λ(ω32(p23+q23)−h)The intersection of the two-dimensional surface U+q3p3 with the NHIM (24) becomes
M(h)∩U+q3p3={(q1,p1,q2,p2,q3,p3)|q1=0,p1=0,p2=0,q2=0,˙p1>0:ω32(p23+q23)=h}.Thus, the NHIM is the circle of radius √2h/ω3 and marked by a dashed line in the LD plot (Fig.fig:8).
Next, the intersection of the two-dimensional surface with the unstable (25) and stable manifolds (26) is given by
Wu(M(h))∩U+q3p3={(q1,p1,q2,p2,q3,p3)|q1=p1,p1=0,q2=0,p2=0,˙p1>0:ω32(p23+q23)=h},Ws(M(h))∩U+q3p3={(q1,p1,q2,p2,q3,p3)|q1=−p1,p1=0,q2=0,p2=0,˙p1>0:ω32(p23+q23)=h},which are one-dimensional for a fixed energy and represent circles of radius √2h/ω3, and marked by dashed red (unstable) and white (stable) lines, respectively, in Fig. fig:8. The only points of minima and singularity in the LD plot (Fig. fig:8) is along the circle and thus identify the manifolds of the NHIM.
The influence of index two saddle points on reaction dynamics has been studied in (missing reference). The construction of a dividing surface for general index k saddles was given in [1]. However, in this section we describe a ''global dividing surface'' that is associated with an index two saddle point and two index one saddle points. Such a structure was constructed in analyzing the isomerization dynamics of a buckled nanobeam in [2]. Intriguingly, a similar geometrical structure arises in the study of the ''roaming'' phenomenon in (missing reference). Consequently, this type of geometrical structure could be more widespread in reaction dynamics so we believe it may be useful to give an analytically tractable example of such a situation.
Unlike our linear examples above, for the system to have multiple saddle points it must be nonlinear. We will consider two identical uncoupled two well potential systems, as shown in Fig. fig:9.
The Hamiltonian for this system is given by:
H=p212−q212+q414⏟H1+p222−q222+q424⏟H2,with the associated Hamiltonian vector field:
˙q1=∂H∂p1=p1,˙p1=−∂H∂q1=q1−q31,˙q2=∂H∂p2=p2,˙p2=−∂H∂q2=q2−q32.The potential energy is given by:
V(q1,q2)=−q212+q414−q222+q424.The potential energy has nine critical points, which we list below, along with their stability type and total energy:
(0,0),index two saddle,total energy0,(0,1),(0,−1),(1,0),(−1,0),index one saddles,total energy−1/4,(1,1),(1,−1),(−1,1),(−1,−1),minima,total energy−1/2,We illustrate the critical points of (44) in Fig. fig:10.
In order to consider a surface to be a dividing surface we need to understand what the surface is dividing. In the language of chemical reactions, index one saddles give rise to dividing surfaces that divide reactants and products. The examples above were designed only to illustrate the geometry associated with the passage of trajectories through a dividing surface. In particular, in the examples ''reactants'' corresponded to a region have a particular sign of the q1 coordinate and ''products'' corresponded to a region corresponding to the opposite sign of the q1 coordinate. and it was arbitrary which region was considered reactants and which products.
References
- P. Collins, G. S. Ezra, and S. Wiggins, “Index k saddles and dividing surfaces in phase space with applications to isomerization dynamics,” The Journal of Chemical Physics, vol. 134, no. 24, p. 244105, 2011.
- P. Collins, G. S. Ezra, and S. Wiggins, “Isomerization dynamics of a buckled nanobeam,” Phys. Rev. E, vol. 86, p. 056218, 2012.