Hamiltonian pitchfork bifurcation in transition across index-1 saddles
HH AMILTONIAN PITCHFORK BIFURCATION IN TRANSITIONACROSS INDEX -1 SADDLES
Wenyang Lyu ∗ , Shibabrat Naik † , Stephen Wiggins ‡ School of Mathematics, University of BristolFry building, Woodland Road, Bristol BS8 1UG, UK A BSTRACT
We study the effect of changes in the parameters of a two dimensional potential energysurface on the phase space structures relevant for chemical reaction dynamics. The changesin the potential energy are representative of chemical reactions such as isomerizationbetween two structural conformations or dissociation of a molecule with an intermediate. Wepresent a two degrees of freedom quartic Hamiltonian that shows pitchfork bifurcation whenthe parameters are varied and we derive the bifurcation criteria relating the parameters.Next, we describe the phase space structures — unstable periodic orbits and its associatedinvariant manifolds, and phase space dividing surfaces — for the systems that can showtrajectories undergo reaction defined as crossing of a potential energy barrier. Finally, wequantify the reaction dynamics for these systems by obtaining the directional flux and gaptime distribution to illustrate the dependence on total energy and coupling strength betweenthe two degrees of freedom. K eywords Hamiltonian pitchfork bifurcation, Bifurcation of a potential energy surface, Gap times oftransitions, Phase space dividing surface, Normally hyperbolic invariant manifold, Global stable/unstablemanifold computation
Understanding, predicting, and controlling transition dynamics has significance in many physical, chemical,biological, and engineering systems. For example, ionization of a hydrogen atom under an electromagneticfield in atomic physics [14], transport of defects in solid state and semiconductor physics [8], chemical reactionrates and pathways in chemical physics [17, 40], dynamic buckling or collapse in structural mechanics [4, 43,28], ship capsize [31, 29, 24], escape and recapture of comets and asteroids in celestial mechanics [15, 7, 27],and escape into inflation or re-collapse to singularity in cosmology [6]. In these systems, understandingtransitions between stable configurations constitutes finding the phase space structures and analysingtheir stability. These phase space structures are equilibrium points, periodic orbits, and invariant manifolds(surfaces or curves in the phase space on which an initial condition once started stays forever) in the phasespace of the governing equations of motion. In a typical scenario of transition dynamics, these phase spacestructures form the skeleton for the transition by guiding the trajectories between stable regions in the phase ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ m a t h . D S ] F e b pace [40, 30]. This indicates that studying the changes in the stability of these phase space structures, thatis bifurcation analysis, gives an understanding of the dependence of the transitions on the parameters.In the context of reaction dynamics, transition between stable regions in the phase space represents reactionby breaking and forming of chemical bonds. The dynamical systems perspective of chemical reactionsentails describing the breaking and forming of chemical bonds using the geometry of phase space structures.The starting point for analysing the geometry and tracking any changes in the stability of these structures isthe bifurcation of equilibrium points of the vector field. In the setting of Hamiltonian models where the totalenergy (or the Hamiltonian) is a sum of kinetic plus potential energy, the equilibrium points of the vector fieldare also the critical points on the potential energy surface (PES). The potential energy in the Hamiltonian isobtained from the electronic structure calculations performed under the Born-Oppenheimer approximationwhich separates the low frequency atomic nuclei motion from the high frequency electronic motion [35]. Thetopographical features of a given potential energy function such as the mountain top, valley floor, valleyridge, and saddles influence the reaction in a dynamical sense, that is their presence is felt for non-zeromomentum (or kinetic energy) even though they are configuration space (coordinates that define the potentialenergy) features. Thus, it is natural to expect that the transition dynamics is affected by the changes in thefeatures of the PES when the parameters of the Hamiltonian model are varied; for example solvent mass insolution-phase reaction [12], energy of the excitation and total energy of the molecule [11].The issue of bifurcation in Hamiltonian models of chemical reactions has received some attention in thecontext of specific reactions such as collinear HgI → HgI + I reaction [1], collinear and non-collinear H + H exchange reaction [19, 13], bimolecular reactions [22]. In these studies, the phase space dividingsurface [41, 42, 36] undergoes bifurcation as the total energy is increased, aside from Ref. [19] who alsoconsidered the role of bath modes on the location of the phase space dividing surface. Ref. [22] studied thechanges in the topology of the transition state (which in general is called the NHIM for N-DOF systems) asthe total energy in increased. From the dynamical systems point of view and as also present in Wigner’sformulation for two DOF systems in the phase space [41, 42], the correct transition state is the unstableperiodic orbit (UPO) which depends on the total energy of the system, and thus its stability can undergobifurcation. Furthermore, as the total energy is increased the region around the saddle equilibrium pointon the energy surface can deviate from the bottleneck geometry. This can manifest in loss of the locallynon-recrossing [32] property required by the phase space dividing surface and thus leads to breakdownof the transition state theory [1, 13]. So studying the effect of total energy on reaction dynamics usingcompletely integrable (where the degrees-of-freedom are separable or uncoupled) Hamiltonian systemsexplains the changes in the geometry of the phase space structures relevant to the reactions [23]. Here,we will consider a two DOF Hamiltonian as a simple model of a reaction in a bath (liquid solvent) where thebifurcation of the unstable periodic orbit, called Hamiltonian pitchfork bifurcation, occurs due to changes inthe parameters.The description of a reacting molecular system is highly influenced by the topography of the potential energysurface. Potential wells describe “products” and saddles are the “transition points” between wells. This isa general mathematical structure of the transition dynamics, and therefore it would seem reasonable totransfer the language of “chemical transformation” to other fields with a similar potential energy surfacedescription. One such physical system and which is of engineering interest is the buckling of beams andcolumns under compression. In modeling and analysing these structural elements, transition dynamicsarises in studying the motion between stable configurations that represent mode shapes of these structuralelements [18, 2, 3, 4, 28, 43]. For example, the analysis of a dynamical model of a nanobeam undercompression in Ref. [4] considers the classical Euler-Bernoulli beam equation subject to compressive stressapplied at both ends. The authors considered a two-mode truncation of the forced beam equation. Moreprecisely, the authors considered parameter regimes where the first mode is unstable and the secondmode can be either stable or unstable, and the remaining (neglected) modes are always stable. Material arameters were used corresponding to a silicon nanobeam. The two-mode model Hamiltonian is the sumof a (diagonal) kinetic energy term and a potential energy term. The form of the potential energy functionthat the authors obtained was analogous with isomerization reactions in chemistry, where “isomerization” inthis case corresponded to a transition between two stable beam configurations. Thus, the dynamics of thebuckled beam was studied using the conceptual framework for the theory of isomerization reactions. Whenthe second mode is stable the potential energy surface has an index one saddle, and when the second modeis unstable the potential energy surface has an index two saddle and two index one saddles. Symmetry ofthe system allowed the authors to readily construct a phase space dividing surface between the two “isomers”(buckled states); and they rigorously proved that, in a specific energy range, it is a normally hyperbolicinvariant manifold [38]. The energy range is sufficiently wide that we were able to treat the effects of the indexone and index two saddles on the isomerization dynamics in a unified fashion. We computed reactive fluxes,mean gap times, and reactant phase space volumes for three stress values at several different energies. Inall cases the phase space volume swept out by isomerizing trajectories is considerably less than the reactantdensity of states, proving that the dynamics is highly non-ergodic. This was supported by the associatedgap time distributions with one or more “pulses” of trajectories. Computation of the reactive flux correlationfunction showed no sign of a plateau region; rather, the flux exhibited oscillatory decay, indicating that, for thetwo-mode model Hamiltonian in the physical regime considered, a rate constant for “isomerization” (buckling)does not exist. These insights are possible by bringing the developments in methods for studying chemicalreactions to dynamical buckling of structures. In this article, we analyse the Hamiltonian pitchfork bifurcationrelevant for transition dynamics where changes in the parameters of the potential energy shows dynamics ona closed double-well, a closed single-well, and an open single-well system.In Sect. 2, we describe the Hamiltonian model and analyse the linear stability of the equilibrium points to findthe bifurcation criteria relating the parameters. In Sect. 3, we discuss the phase space structures governingthe transition dynamics and obtain quantitative measures of the transition using the phase space structures.Finally, in Sect. 4, we present our conclusions of this work and outlook for related future work. We consider a two degrees of freedom (DOF) Hamiltonian where a quartic potential is coupled to a harmonicoscillator with a bilinear coupling. H ( x, y, p x , p y ) = T ( p x , p y ) + V ( x, y ) = 12 (cid:32) p x m s + p y m b (cid:33) − α x + β x + ω y + ε x − y ) (1)This form of a model Hamiltonian is borrowed from the solution-phase reactions literature where a finite sumof harmonic oscillators are referred to as the “bath” coupled with the quartic Hamiltonian referred to as the“system”. Furthermore, we define the forward transition (or forward reaction) as a trajectory crossing from x > to x < , and the backward transition (or backward reaction) is defined as crossing from x < to x > .The Hamilton’s equations are given by ˙ x = ∂ H ∂p x = p x m s ˙ y = ∂ H ∂p y = p y m b ˙ p x = − ∂ H ∂x = α x − β x + ε ( y − x )˙ p y = − ∂ H ∂y = − ωy + ε ( x − y ) (2) he equilibria for this system are located at x e = (0 , , , and ± x e = ± ( x e , y e , , . where x e = (cid:115) β (cid:18) α − ωεω + ε (cid:19) > , y e = x e (cid:18) εω + ε (cid:19) > . (3)where we require either β < and α < ωε/ ( ω + ε ) or β > and α > ωε/ ( ω + ε ) for the system to have threeequilibria. In the case when β < and α > ωε/ ( ω + ε ) or β > and α < ωε/ ( ω + ε ) , the system admits oneequilibrium point at x e = (0 , , , . The Jacobian of the Hamiltonian vector field is given by J ( x, y, p x , p y ) = ∂ H ∂x∂p x ∂ H ∂y∂p x ∂ H ∂p x ∂ H ∂p y ∂p x ∂ H ∂x∂p y ∂ H ∂y∂p y ∂ H ∂p x ∂p y ∂ H ∂p y − ∂ H ∂x − ∂ H ∂y∂x − ∂ H ∂p x ∂x − ∂ H ∂p y ∂x − ∂ H ∂x∂y − ∂ H ∂y − ∂ H ∂p x ∂y − ∂ H ∂p y ∂y = α − ε − βx ε ε − ω − ε (4)The eigenvalues of J ( x, y, p x , p y ) evaluated at x e are ±√ λ , ± i √− λ , where λ , = α − ω − ε ± (cid:112) ( ω + α ) + 4 ε . (5)The term inside the square root is positive which means that λ is negative for α ≤ . λ is also negative for α > as the term inside the square root is larger than or equal to α when α > . The critical value of α suchthat λ = 0 occurs when α = ωεω + ε . (6)The eigenvalues of J ( x, y, p x , p y ) evaluated at ± x e are ±√ λ , ± i √− λ , where λ , = − α − ω − ε + 3 ωε/ ( ω + ε ) ± (cid:112) [2 α − ( ω + 3 ωε/ ( ω + ε ))] + 4 ε (7) = − α − ( ω + 2 ε ) / ( ω + ε ) ± (cid:112) [2 α − ( ω + 3 ωε/ ( ω + ε ))] + 4 ε . (8)The term inside the square root is positive which means that λ is negative for α ≥ . λ is also negative for α < as the term inside the square root is larger than or equal to − α when α < . The critical value of α such that λ = 0 occurs when α = ωεω + ε . (9)We note that the existence of the critical value of α depends on the sign of α and will be discussed in thefollowing sections. We also get the expected result that the eigenvalues are independent of β since thelinearization of the vector field at the equilibria only depends on the quadratic terms in the Hamiltonian. he total energy of the equilibrium points are H ( x e ) =0 , (10) H ( ± x e ) =( x e ) (cid:18) − α β x e ) + ωε ω + ε ) (cid:19) = − β (cid:18) α − ωεω + ε (cid:19) . (11)We note that the depth of the potential well is simply the difference between the saddle equilibrium point andthe center equilibrium point. In the system considered here, the depth is the energy β (cid:16) α − ωεω + ε (cid:17) . Wealso denote the excess energy ∆ E as the energy difference between the total energy of the system, e andthe energy of the saddle equilibrium point. Now, we classify the linear stability type of the equilibria. Case α < εω/ ( ε + ω ) : β < . Three eqilibrium points, one at x e = (0 , , , . λ < , λ < which means that x e a centre xcentre equilibrium point. The other two equilibrium points are at ± x e = ± ( x e , y e , , . λ > , λ < which means that ± x e are saddle-centre equilibrium points.2. β ≥ . One equilibrium point at x e = (0 , , , . λ < , λ < which means that x e is a centre-centre equilibrium point. Case α = εω/ ( ε + ω ) : β = 0 . Line of equilibrium points at y = x ( εω + ε ) β (cid:54) = 0 . One equilibrium point at x e = (0 , , , . λ = 0 , λ < which means that x e is anonhyperbolic equilibrium point. Case α > εω/ ( ε + ω ) : β ≤ . One equilibrium point at x e = (0 , , , . λ > , λ < which means that x e is a saddle-centre equilibrium point.2. β > . Three eqilibrium points, one at x e = (0 , , , . λ > , λ < which means that x e is asaddle-centre equilibrium point. The other two equilibrium points are at ± x e = ± ( x e , y e , , , with λ < , λ < which means that ± x e are centre-centre equilibrium points.The summary of the bifurcation analysis of the equilibrium points associated to λ or λ for the coupledsystem is shown in Fig. 1. The linearised stability of the equilibrium points associated to λ or λ is notshown as it is guaranteed to have centre type stability. A preliminary step in studying the global dynamics of the double well system shown inFig. 2(b) is to compute the Poincaré section of trajectories by defining the two dimensional surface U + xp x = (cid:8) ( x, y, p x , p y ) ∈ R | y = 0 , ˙ y ( x, y, p x ; E ) > (cid:9) (12)and the resulting Poincaré surface of sections are shown in Fig. 3. We can observe a mixture of regular andchaotic dynamics for the illustrative values of parameters and for the same total energy, E = 0 . , used inthe Poincaré sections. Although, the Poincaré sections reveal the structure of the global dynamics, it can notbe applied to open potential wells since the trajectories that leave the potential well (either in short or long = ωεω + ε α β I IIIII IV (a) α x e I II (b) α x e III IV (c)
Figure 1: (a) Bifurcation diagram for the coupled system with ω = 1 , ε = 0 . . Each subfigure correspond tothe phase space geometry for 4 representative choices(I: α < εω/ ( ε + ω ) , β > , II: α > εω/ ( ε + ω ) , β > ,III: α < εω/ ( ε + ω ) , β < , IV: α > εω/ ( ε + ω ) , β < ) of α, β in the x − p x plane. Bifurcation diagram for thecoupled system with (b) β = 1 , (c) β = − . Other parameters ω = 1 , ε = 0 . . (a) (b) (c) (d) Figure 2: Changes in the topology of the potential energy surface with changes in the parameters, α, β .The green, black, and magenta curves represent the equipotential contours at total energy . , , − . ,respectively. (a-d) Typical potential energy surfaces for Cases I-IV , respectively.times) become indistinguishable due to white regions in the Poincaré surface of section. In addition, whenthe dynamics occurs on a surface of dimension greater than two, a two-dimensional Poincaré section cannot reveal the global dynamics because trajectories may not cross the surface of section. In this case, acodimension-1 Poincaré section becomes impossible to visualize completely without taking further sections orprojections. For these settings of open potential wells and high-dimensional phase space, another trajectorydiagnostic is to be relied upon and related results and discussions can be found in Refs. [16, 5, 25, 26].
Phase space dividing surface and invariant manifolds associated with the unstable periodic orbit.
For the uncoupled system, ε = 0 , the bifurcation condition becomes α = 0 . In addition, the uncoupledsystem is integrable and we can derive an explicit analytical formula for the phase space dividing surface(PSDS). We note here that the dividing surface constructed in the phase space has the locally no-recrossingproperty [32] and in general, trajectories will show global recrossings of the PODS due to the Poincarérecurrence theorem [37] when the energy surface is bounded.1. Case II and IV : When α > and β ≷ , the dividing surface anchored at the UPO (which is anormally hyperbolic invariant manifold (NHIM) [39] in a two DOF system) is defined by the condition x = 0 , p x = 0 on the energy surface. The dividing surface which is a 2-sphere ( S ) has forward a) (1 , , . , . (b) (1 , , , . (c) (1 , , , . (d) (1 , , . , . Figure 3: Poincaré surface of sections for the double-well system with different system, bath, and interactioncoefficients to illustrate the coexistence of regular and chaotic trajectories in an isomerizing system. All thesystems are at same excess energy, ∆ E = 0 . , m s = m b = 1 . with ( α, β, ω, (cid:15) ) value shown as tuples inthe caption for the surface defined by Eqn. (12). We note that even though the total energy is constant in allthe cases shown here the bottleneck width for these parameters is altered because of the changes in thepotential energy surface.( p x ( y, p y ; e ) < ) and backward ( p x ( y, p y ; e ) > ) hemispheres given byforward/backward DS: (cid:8) ( x, y, p x , p y ) ∈ R | x = 0 , p x ( y, p y ; e ) ≶ (cid:9) (13)where e is the total energy of the system and p x ( y, p y ; e ) = ± (cid:113) e − p y − ωy is calculated from thefixed energy constraint. The forward reaction occurs when the trajectory crosses the forward DSand the backward reaction occurs when the trajectory crosses the backward DS. In case IV, which y p y (a) x y p y (b) Figure 4: DS (blue, red disks), NHIM(black) on the energy surface for the coupled system with (a) α = β =1 , ε = 0 . (b) negative α = β = − , ε = 0 . . Other parameters are chosen as ∆ E = 0 . , ω = 1 .is simply a two DOF saddle, the energy surface is unbounded (open potential well), trajectoriesgo-off to infinity and do not return to recross the PODS. However, trajectories will recross the PODSin case II due to Poincaré recurrence theorem and we will characterize this using a quantitativemeasure of this recrossing called gap times [10]. We note here that the global recrossing property isindependent of the existence of chaos in a dynamical system. The uncoupled system is an exampleof that since it is integrable, and hence can not exhibit chaotic dynamics.2. Case III : When α < and β < , there are two index-1 saddles, and hence two corresponding UPOs( Γ , ) and two PSDSs. The UPOs are defined by the condition x = ± (cid:112) α/β, p x = 0 on the energysurface and the PSDSs are given by:forward/backward DS − : (cid:110) ( x, y, p x , p y ) ∈ R | x = (cid:112) α/β, p x ( y, p y ; e ) ≶ (cid:111) (14)forward/backward DS + : (cid:110) ( x, y, p x , p y ) ∈ R | x = − (cid:112) α/β, p x ( y, p y ; e ) ≶ (cid:111) (15)where e is the total energy of the system and p x ( y, p y ; e ) = ± (cid:113) e − p y − ωy is calculated from thefixed energy constraint.We note that the case I α < , β > only has one equilibrium point of centre-centre stability, so the systemdoes not show transition (or reaction) dynamics as defined in for this model Hamiltonian.For the coupled system, the phase space dividing surface is sampled using the algorithm in Ref. [9] afterobtaining the unstable periodic orbit for a given energy, Γ( E ) using the open-source python package,uposham [20]. Let us denote the unstable periodic orbit by Γ( E ) := ( x, y, p x , p y ) ∈ R or by Γ − , + when thereare two unstable periodic orbits at x < and x > .1. Case II and IV : When α > ωε/ ( ω + ε ) , the PODS is given byforward/backward DS: (cid:8) ( x, y, p x , p y ) ∈ R | p x ( x, y, p y ; e ) ≶ (cid:9) (16)where p x = ± (cid:113) e − p y − V ( x, y )) is the fixed energy constraint and x, y ∈ UPO . . Case III : When α < ωε/ ( ω + ε ) , β < , the PSDSs are given byforward/backward DS − : (cid:8) ( x, y, p x , p y ) ∈ R | p x ( x − , y − , p y ; e ) ≶ (cid:9) (17)forward/backward DS + : (cid:8) ( x, y, p x , p y ) ∈ R | p x ( x + , y + , p y ; e ) ≶ (cid:9) (18)where p x = ± (cid:113) e − p y − V ( x − , + , y − , + )) is the fixed energy constraint and x − , + , y − , + ∈ Γ − , + .The phase space dividing surfaces on the energy surface for the coupled systems are shown in Fig. 4 . Wenote that the PODS have a geometry of discs in the two configuration space coordinates with one momentumcoordinate.From a phase space perspective of transition dynamics across index-1 saddles, the phase space dividingsurface can be viewed as a local barrier between reactive and non-reactive trajectories. However, the globalbarriers between the two type of trajectories are the invariant (stable and unstable) manifolds associatedwith the UPO. These stable and unstable invariant manifolds have a cylindrical geometry, that is R × S , andform impenetrable barriers between the reactive and non-reactive trajectories and shown in Fig. 5. We use the definition of the gap time which is defined as the time interval between two successive crossingsof the phase space DS [10]. For case II , that is when β > , α > ( ωε ) / ( ω + ε ) , the gap time distribution iscalculated by launching trajectories from the DS with p x > such that the trajectory enters the reactantsregion. Then, the time interval until the next crossing of the DS represents the gap time or the first passagetime. We note that local recrossings are impossible since the initial conditions are on a phase space dividingsurface. For this case with double-well and a bounded energy surface, when a trajectory recrosses the DSand enters the products region, the x − momentum changes sign and x < . This happens when the trajectorycrosses the forward DS given by p x < in the Eqn. (16). We perform this calculation for an ensemble ofinitial conditions at a fixed energy and we record the gap time for all the trajectories starting on the PSDSwith positive p x within 100 time units (same for all the case discussed here), and call this the microcanonicalgap time distribution [10].For case III , that is when β < , α < ( ωε ) / ( ω + ε ) , there are two phase space dividing surfaces associatedwith the two index-1 saddles located at x < and x > . Thus, for each PSDS there are two gap timedistributions: one for crossing the PSDS where the ensemble gets initialized and one for crossing the otherPSDS. A comparison of these two gap times is beyond the scope of this article but interested readers can seean examples of a similar analysis for a three degrees of freedom model Hamiltonian of HCN isomerization inRef. [33]. We note here that out of the two phase space dividing surfaces, we choose DS − to initialize anensemble at a fixed energy with p x > and integrate until they either cross DS − or DS + . We record only thetime to cross the PSDS DS + within 100 time units.We use the algorithm in Ref. [9] to sample points on the phase space DS (which is a hemisphere of the2-sphere, S ) to calculate the gap time distribution. Firstly , we obtain the UPO associated with the index-1saddle (or saddles in case III) equilibrium point for a given energy e . Second , we use the configuration spaceprojection of the UPO { ( x i , y i ) } . For each point on the configuration space { ( x i , y i ) } , the ( p x , p y ) satisfiesthe circle given by
12 ( p x + p y ) = e − V ( x i , y i ) . (19)Therefore, the maximum value of p x is p maxx = (cid:112) e − V ( x i , y i )) (20) a)(b) (c) Figure 5:
Cylindrical invariant manifolds and energy surface for the three reactive systems due to thepitchfork bifurcation. Stable (blue) and unstable (red) manifolds of the unstable periodic orbit associated withthe index-1 saddle (NHIM for two DOF systems) in the bottleneck while the energy surface is shown as thegrey surface. (a)
Case II: α = 1 . , β = 1 . (b) Case III: α = 0 . , β = − . (c) Case IV: α = 1 . , β = − . .Other parameters, ω = ε = 1 . , and excess energy, ∆ E = 0 . , are same for all cases.and we select p x,i uniformly from the interval [ − p maxx , p maxx ] . Using the definition of the Hamiltonian, wecalculate the value of p y,i , which is either positive or negative. We note that p maxy = p maxx , p y ∈ (cid:2) − p maxy , p maxy (cid:3) . Thirdly , we generate a set of microcanonical ensemble { ( x i , y i , p x,i , p y,i ) } e on the phase space DS at atotal energy e . We select the initial conditions with positive p x and integrate so that the trajectories enterthe reactants ( x > in case II or between the index-1 saddles in case III) region, and leave the reactantsregion by either crossing the same phase space DS (as in case II) or the other phase space DS (as in caseIII). Finally, the time interval between entering and exiting the reactants region is the gap time of the initialcondition. For the coupled systems, the UPO is computed using an open-source python library for computingunstable periodic orbits in two DOF Hamiltonian systems [20]. Thus, for case II, to simplify this detection ofcrossing the forward DS, we use an alternate boundary at x = − line which ensures that the trajectory hasleft the reactants region and is in the products region. The gap time distributions for ∆ E = 0 . , . , . areshown in Fig. 6(a) and 6(b) for the uncoupled and coupled systems, respectively. For the case III, we use analternate boundary at x = 1 . which ensures that the trajectory has left the reactants region and is in the roducts region. The gap time distributions for ∆ E = 0 . , . , . are shown in Fig. 7(a) and 7(b) for theuncoupled and coupled systems, respectively. Gap time F r a c t i on o f t r a j e c t o r i e s ∆ E=0.01 ∆ E=0.05 ∆ E=0.1 (a)
Gap time F r a c t i on o f t r a j e c t o r i e s ∆ E=0.01 ∆ E=0.05 ∆ E=0.1 (b)
Figure 6:
Gap time distributions for case II for the (a) uncoupled (cid:15) = 0 and the (b) coupled (cid:15) = 0 . system.Parameters in the potential energy are chosen as α = 1 , β = 1 , ω = 1 . Gap time F r a c t i on o f t r a j e c t o r i e s ∆ E=0.01 ∆ E=0.05 ∆ E=0.1 (a)
Gap time F r a c t i on o f t r a j e c t o r i e s ∆ E=0.01 ∆ E=0.05 ∆ E=0.1 (b)
Figure 7:
Gap time distributions for case III for the (a) uncoupled (cid:15) = 0 (b) and the coupled (cid:15) = 0 . system.Parameters in the potential energy are chosen as α = − , β = − , ω = 1 .We observe some common features in the gap time distributions for both the cases. We see from the figuresthat the gap time distributions have a unimodal shape for the uncoupled system except for a time shift of thepeak of the distribution. For the coupled systems, the gap time distributions have multiple pulses and thetime of the first pulse appears at time instant that decreases with increasing excess energy. Thus, as ∆ E increases, most trajectories cross the phase space DS to the products region at short timescales comparedto the timescale in which all the trajectories exit the reactant region. This is observed across systems withsingle and multiple phase space dividing surfaces [10, 33]. .3 Directional flux of the transition The directional flux of transition trajectories through the phase space dividing surface is given by theaction of the normally hyperbolic invariant manifold (NHIM) which is of dimension S N − in a N degrees offreedom [32, 34]. In the two DOF system, the NHIM is an unstable periodic orbit (UPO) and the directionalflux given by the action simplifies to the line integral: Q = (cid:90) UPO p · d q . (21)For the uncoupled systems, the UPO can be expressed implicitly by substituting the condition x = 0 , p x = 0 in the Hamiltonian.1. Case II and IV:
UPO: (cid:26) ( x, y, p x , p y ) ∈ R | x = 0 , p x = 0 , p y + ω y = e (cid:27) (22)Thus, the directional flux, Q , is given by Q = (cid:90) UPO p · d q = (cid:90) T p · d q dt dt = (cid:90) T ( p x ˙ x + p y ˙ y ) dt = (cid:90) T p y ˙ y dt = (cid:90) T p y dt = T e = 2 πe √ ω (23)where T is the time period of the UPO.2. Case III:
When α < and β < , the directional flux through either DS, is given by Q = (cid:90) UPO p · d q = (cid:90) T p · d q dt dt = (cid:90) T ( p x ˙ x + p y ˙ y ) dt = (cid:90) T p y ˙ y dt = (cid:90) T p y dt = T ∆ E = 2 π ∆ E √ ω (24)where T is the time period of the UPO.In Eqns. (23), (24), p x ˙ x vanishes because the vector field is tangential on the UPO. We calculate the integralby expressing p y in terms of t , which can be solved explicitly for the uncoupled Hamilton’s equations ofmotion (2). We note that this expression is the N = 2 case for the flux formula in Ref. [32].We note the directional flux is π ∆ E √ ω all three uncoupled systems. For Case II and IV, ∆ E = e − H ( x e ) = e where e is the total energy of the system and H ( x e ) = 0 is the energy of the saddle equilibrium point.For Case III, ∆ E = e − H ( ± x e ) = e + 14 β (cid:16) α − ωεω + ε (cid:17) where e is the total energy of the system and H ( ± x e ) = − β (cid:16) α − ωεω + ε (cid:17) is the energy of the saddle equilibrium points. Case II-IV:
For the coupled systems, the directional flux through the phase space DS is computed byevaluating the integral Q = (cid:90) UPO p · d q = (cid:90) T p d q dt dt (25) Q = (cid:90) T ( p x ˙ x + p y ˙ y ) dt = (cid:90) T ( p x ˙ x + p y ˙ y ) dt (26) Q = (cid:90) T p x + p y dt (27)In Fig. 8(a) we present the directional flux, Q through the DS for α = 1 , β = 1 and ω = 1 with variation of ∆ E . We can see that as ∆ E increases, Q increases linearly. At a fixed value of ∆ E , Q decreases as weincrease the strength of the coupling ε . In Fig. 8(b) we present the directional flux, Q through the DS for .00 0.02 0.04 0.06 0.08 0.10 ∆ E D i r e c t i ona l f l u x , Q † = 0 † = 0 . † = 1 (a) ∆ E D i r e c t i ona l f l u x , Q † = 0 † = 0 . † = 0 . (b) Figure 8: Directional flux for the transition through the phase space dividing surface for (a) case II: α =1 , β = 1 , ω = 1 , (b) case III and DS − : α = − , β = − , ω = 1 . α = − , β = − and ω = 1 with variation of ∆ E . We can also see that as ∆ E increases, Q increases linearly.At a fixed value of ∆ E , Q decreases as we increase the strength of the coupling ε . In this article, we analyzed the Hamiltonian pitchfork bifurcation due to changes in parameters of the potentialenergy. We related the qualitative changes in the topology of the potential energy surface with the changes inthe linear stability of the equilibrium points. We obtained the phase space structures relevant for the transitionacross an index-1 saddle by computing the invariant manifolds. These are essentially the unstable periodicorbit (a normally hyperbolic invariant manifold in two DOF systems) associated with the index-1 saddle thatanchors both the phase space dividing surface and the phase space conduits of transition across the index-1saddle/s. Finally, we quantified the transition dynamics for the sub-systems (typical PES for the three casesII, III, IV) that arise due to the bifurcation by obtaining the gap time distributions and directional flux.In future, related work will focus on correlating the statistical properties of the transition with geometry ofthe phase space structures. One such question involves characterizing the pulses in gap time distributionswith the homoclinic and heteroclinic tangle due to the stable and unstable manifolds associated with theunstable periodic orbit. Furthermore, the quartic potential coupled with a harmonic oscillator admits a simpleparametrization of the depth of the potential well which can be of use for a qualitative and quantitativeanalysis of the influence of depth of the PES on the reaction dynamics quantities; similar to Ref. [21].
ACKNOWLEDGMENTS
We acknowledge the support of EPSRC Grant No. EP/P021123/1 and Office of Naval Research (Grant No.N00014-01-1-0769). eferences [1] I. Burghardt and P. Gaspard. The Molecular Transition State: From Regular to Chaotic Dynamics. The Journal of Physical Chemistry , 99(9):2732–2752, Mar. 1995. ISSN 0022-3654, 1541-5740. doi: .[2] S. M. Carr, W. E. Lawrence, and M. N. Wybourne. Buckling cascade of free-standing mesoscopicbeams.
Europhysics Letters (EPL) , 69(6):952–958, Mar. 2005. ISSN 0295-5075, 1286-4854. doi: .[3] A. Chakraborty. Buckled nano rod – a two state system and quantum effects on its dynamics usingsystem plus reservoir model.
Molecular Physics , 109(4):517–526, Feb. 2011. doi: .[4] P. Collins, G. S. Ezra, and S. Wiggins. Isomerization dynamics of a buckled nanobeam.
Phys. Rev. E ,86:056218, Nov 2012.[5] G. T. Craven, A. Junginger, and R. Hernandez. Lagrangian descriptors of driven chemical reactionmanifolds.
Physical Review E , 96(2):022222, 2017.[6] H. P. de Oliveira, A. M. Ozorio de Almeida, I. Damião Soares, and E. V. Tonini. Homoclinic chaos in thedynamics of a general Bianchi type-IX model.
Phys. Rev. D , 65(8):9, 2002.[7] M. Dellnitz, O. Junge, M. W. Lo, J. E. Marsden, K. Padberg, R. Preis, S. D. Ross, and B. Thiere.Transport of Mars-crossing asteroids from the quasi-Hilda region.
Physical Review Letters , 94:231102,2005.[8] B. Eckhardt. Transition state theory for ballistic electrons.
J. Phys. A-Math. Gen. , 28(12):3469, 1995.[9] G. S. Ezra and S. Wiggins. Sampling phase space dividing surfaces constructed from normallyhyperbolic invariant manifolds (nhims).
The Journal of Physical Chemistry A , 122(42):8354–8362, 2018.doi: .[10] G. S. Ezra, H. Waalkens, and S. Wiggins. Microcanonical rates, gap times, and phase space dividingsurfaces.
The Journal of Chemical Physics , 130(16):164118, 2009. doi: .[11] S. C. Farantos, R. Schinke, H. Guo, and M. Joyeux. Energy Localization in Molecules, BifurcationPhenomena, and Their Spectroscopic Signatures: The Global View.
Chemical Reviews , 109(9):4248–4271, Sept. 2009. ISSN 0009-2665, 1520-6890. doi: . URL https://pubs.acs.org/doi/10.1021/cr900069m .[12] R. Garcia-Meseguer, B. Carpenter, and S. Wiggins. The influence of the solvent’s mass on the locationof the dividing surface for a model Hamiltonian.
Chemical Physics Letters: X , 3:100030, July 2019.ISSN 25901419. doi: .[13] M. Inarrea, J. F. Palacian, A. I. Pascual, and J. P. Salas. Bifurcations of dividing surfaces in chemicalreactions.
J. Chem. Phys. , 135:014110, 2011.[14] C. Jaffé, D. Farrelly, and T. Uzer. Transition state theory without time-reversal symmetry: chaoticionization of the hydrogen atom.
Phys. Rev. Lett. , 84:610–613, 2000.[15] C. Jaffé, S. D. Ross, M. W. Lo, J. E. Marsden, D. Farrelly, and T. Uzer. Statistical theory of asteroidescape rates.
Physical Review Letters , 89(1):011101, 2002.[16] A. Junginger, L. Duvenbeck, M. Feldmaier, J. Main, G. Wunner, and R. Hernandez. Chemical dynamicsbetween wells across a time-dependent barrier: Self-similarity in the lagrangian descriptor and reactivebasins.
The Journal of chemical physics , 147(6):064101, 2017.[17] T. Komatsuzaki and R. S. Berry. Regularity in chaotic reaction paths. I. Ar6.
The Journal of ChemicalPhysics , 110(18):9160–9173, 1999.
18] W. Lawrence. Phonon description and the Euler buckling instability of a mesoscopic bar at fixedstrain.
Physica B: Condensed Matter , 316-317:448–451, May 2002. ISSN 09214526. doi: .[19] C.-B. Li, M. Toda, and T. Komatsuzaki. Bifurcation of no-return transition states in many-body chemicalreactions.
The Journal of Chemical Physics , 130(12):124116, Mar. 2009. doi: .[20] W. Lyu, S. Naik, and S. Wiggins. UPOsHam: A Python package for computing unstable periodic orbitsin two-degree-of-freedom Hamiltonian systems.
Journal of Open Source Software , 5(45):1684, 2020.doi: .[21] W. Lyu, S. Naik, and S. Wiggins. The Role of Depth and Flatness of a Potential Energy Surface inChemical Reaction Dynamics.
Regular and Chaotic Dynamics , 25(5):453–475, Sept. 2020. ISSN1560-3547, 1468-4845. doi: .[22] R. S. MacKay and D. C. Strub. Bifurcations of transition states: Morse bifurcations.
Nonlinearity , 27(5):859, 2014.[23] F. A. L. Mauguière, P. Collins, G. S. Ezra, and S. Wiggins. Bifurcations of normally hyperbolic invariantmanifolds in analytically tractable models and consequences for reaction dynamics.
Int. J. Bifurcat.Chaos , 23(12), 2013.[24] S. Naik and S. D. Ross. Geometry of escaping dynamics in nonlinear ship motion.
Communications inNonlinear Science and Numerical Simulation , 47:48–70, June 2017. doi: .[25] S. Naik and S. Wiggins. Finding normally hyperbolic invariant manifolds in two and three degrees offreedom with Hénon-Heiles-type potential.
Phys. Rev. E , 100:022204, 2019.[26] S. Naik and S. Wiggins. Detecting reactive islands in a system-bath model of isomerization.
PhysicalChemistry Chemical Physics , 22(32):17890–17912, 2020. ISSN 1463-9076, 1463-9084. doi: .[27] S. D. Ross. Statistical theory of interior-exterior transition and collision probabilities for minor bodiesin the solar system. In G. Gómez, M. W. Lo, and J. J. Masdemont, editors,
Libration Point Orbits andApplications , pages 637–652. World Scientific, 2003.[28] J. Sieber, J. W. Hutchinson, and J. M. T. Thompson. Nonlinear dynamics of spherical shells bucklingunder step pressure.
Proceedings of the Royal Society A: Mathematical, Physical and EngineeringSciences , 475(2223):20180884, Mar. 2019. ISSN 1364-5021, 1471-2946. doi: .[29] J. M. T. Thompson and J. R. de Souza. Suppression of escape by resonant modal interactions: in shellvibration and heave-roll capsize.
Proc. R. Soc. Lond. A , 452:2527–2550, 1996.[30] T. Uzer, C. Jaffé, J. Palacián, P. Yanguas, and S. Wiggins. The geometry of reaction dynamics.
Nonlinearity , 15:957–992, 2002.[31] L. N. Virgin. Approximate criterion for capsize based on deterministic dynamics.
Dynamics and Stabilityof Systems , 4(1):56–70, 1989.[32] H. Waalkens and S. Wiggins. Direct construction of a dividing surface of minimal flux for multi-degree-of-freedom systems that cannot be recrossed.
Journal of Physics A: Mathematical and General , 37(35):L435–L445, aug 2004. doi: .[33] H. Waalkens, A. Burbanks, and S. Wiggins. Phase space conduits for reaction in multidimensionalsystems: HCN isomerization in three dimensions.
The Journal of Chemical Physics , 121(13):6207–6225,Oct. 2004. ISSN 0021-9606, 1089-7690. doi: .[34] H. Waalkens, A. Burbanks, and S. Wiggins. A formula to compute the microcanonical volume of reactiveinitial conditions in transition state theory.
Journal of Physics A: Mathematical and Theoretical , 38(45):l759, 2005. doi: .
35] D. Wales.
Energy Landscapes: Applications to Clusters, Biomolecules and Glasses . CambridgeMolecular Science. Cambridge University Press, 2004. doi: .[36] S. Wiggins. On the geometry of transport in phase space I. Transport in k degree-of-freedom Hamiltoniansystems, ≤ k < ∞ . Physica D , 44:471–501, 1990.[37] S. Wiggins.
Introduction to Applied Nonlinear Dynamical Systems and Chaos , volume 2. SpringerScience & Business Media, 2003.[38] S. Wiggins.
Normally hyperbolic invariant manifolds in dynamical systems.
Springer-Verlag, New York,2014. ISBN 978-1-4612-8734-6.[39] S. Wiggins. The role of normally hyperbolic invariant manifolds (nhims) in the context of the phasespace setting for chemical reaction dynamics.
Regular and Chaotic Dynamics , 21(6):621–638, 2016.doi: .[40] S. Wiggins, L. Wiesenfeld, C. Jaffé, and T. Uzer. Impenetrable barriers in phase space.
Phys. Rev. Lett. ,86:5478–5481, 2001.[41] E. P. Wigner. The transition state method.
Trans. Faraday Soc. , 34:29–40, 1938.[42] E. P. Wigner. Some remarks on the theory of reaction rates.
J. Chem. Phys. , 7:646–650, 1939.[43] J. Zhong, L. N. Virgin, and S. D. Ross. A tube dynamics perspective governing stability transitions: Anexample based on snap-through buckling.
Int. J. Mech. Sci. , 000(May):1–16, 2018., 000(May):1–16, 2018.