key: cord-0604925-urptg6jx
authors: Pandey, Priyanka; Naik, Shibabrat; Keshavamurthy, Srihari
title: Classical and quantum dynamical manifestations of index-2 saddles: concerted versus sequential reaction mechanisms
date: 2020-09-16
journal: nan
DOI: nan
sha: 1e592032bc1120ddffb6104e3d0ea21b9d0f8b15
doc_id: 604925
cord_uid: urptg6jx

The presence of higher index saddles on a multidimensional potential energy surface is usually assumed to be of little significance in chemical reaction dynamics. Such a viewpoint requires careful reconsideration, thanks to elegant experiments and novel theoretical approaches that have come about in recent years. In this work, we perform a detailed classical and quantum dynamical study of a model two degree of freedom Hamiltonian, which captures the essence of the debate regarding the dominance of a concerted or a stepwise reaction mechanism. We show that the ultrafast shift of the mechanism from a concerted to a stepwise one is essentially a classical dynamical effect. In addition, due to the classical phase space being a mixture of regular and chaotic dynamics, it is possible to have a rich variety of dynamical behaviour, including a Murrell-Laidler type of mechanism, even at energies sufficiently above that of the index-2 saddle. We rationalize the dynamical results using an explicit construction of the classical invariant manifolds in the phase space.

Transition state theory (TST) is a cornerstone for reaction dynamics [1] . The simple idea of calculating the rate of a reaction as a flux through a "bottleneck" region (TS) has proved to be of immense utility in chemistry for nearly a century now. The notion of associating the TS with a saddle point on the multidimensional potential energy surface (PES) has led to the development of several powerful algorithms to determine the so called intrinsic reaction coordinate (IRC) or minimum energy path (MEP) and accurate ab initio force fields in the vicinity of the various TS. Although the IRCs and their connectivity across the high dimensional PES are useful starting points for analyzing the reactions [2] , the fact remains that MEPs and IRCs are non-dynamical. Consequently, depending on the energy and other parameters of interest, significant deviations from IRCs may be observed [3] . There is little doubt that the true TS is to be found in the full phase space of the system. This dynamical "Wignerian" view was brought out beautifully by the early work by Pollak and Pechukas [4, 5] and more recently generalized to higher degrees of freedom by Wiggins and coworkers [6, 7] .

Although, traditionally TS have been associated with index-1 saddles on the PES, in systems with a large number of degrees of freedom the true (dynamical) TS rarely coincides with the saddle points on the PES. Moreover, in such multidimensional systems there are higher index saddles on the PES that are expected to be dynamically relevant for reactions. Indeed, reactions at sufficiently high energies can occur via pathways that involve both the index-1 saddles as well as the higher index ones. Competition between the different pathways, apart from complicating the mechanistic understanding, can lead to unexpected products. For example, in a recent experimental study, Lu et al. found that a short and intense extreme ultraviolet (XUV) excited CO 2 molecule can result in the production of O 2 molecules [8] . Similarly, the ring opening reaction of cyclobutene to 1, 3butadiene can occur via a disrotatory pathway, as opposed to the conrotatory pathway predicted based on the Woodward-Hoffman rules, due to the involvement of a index-2 saddle. [9, 10] It is crucial to note here that the implication of higher index saddles to reactivity does not necessarily come under the ambit of the Murrell-Laidler theorem [11] [12] [13] or the McIver-Stanton rules [14, 15] . The reason being that these rules are based on the IRC perspective, and hence non-dynamical. Note, however, that even within the IRC perspective it is possible for a index-1 saddle to be linked to a higher index saddle which can then lead to multiple products. An example for this comes from the recent work of Harabuchi et al [16] resulting in the so called nontotally symmetric trifurcation of a reaction pathway. Similarly, modulating reactivity using external forces need not be constrained by IRC-based rules and can possibly involve pathways that utilize higher index saddles. Thus in mechanochemistry [17] [18] [19] [20] [21] the energetics of the various index saddles alone is no longer a useful criteria and it is quite feasible for the "forced" dynamics to prefer pathways that explore the various high index saddle regions.

Recently there have been a number of studies that focus on the dynamical aspects of the higher index saddles. A key motivation is to generalize the dynamical perspective of TST to construct locally recrossing free dividing surfaces. Efforts along these lines have been made by Collins et al. on a model two degrees of freedom potential [22] , by Haller et al. on the double ionization of helium in an external electric field [23] , and by Nagahata et al. on the proton exchange reaction involving the H + 5 moiety [24] . The latter two examples being three degrees of freedom system. We note that typically the higher index saddles occur along with the usual index-1 saddles due to topological constraints on the potential energy surface. An early work by Mann and Hase highlighted the role of dynamics in the electrocyclic ring opening reaction of cyclopropane radical to form the allyl radical [25] . More recently, Pradhan and Lourderaj performed extensive ab initio dynamics calculations to emphasize the key role of a index-2 saddle point in the denitrogenation reaction of 1-pyrazoline [26] . Interestingly, in the study by Mann and Hase the trajectories were not propagated for a sufficiently long time to ascertain if the reaction was indeed statistical or not. On the other hand, Pradhan and Lourderaj did a much more extensive computation, and established significant deviations from the minimum energy path (MEP). They suggested the possibility of slow intramolecular vibrational energy redistribution (IVR) leading to nonstatistical dynamics. One phenomenon where the issue of competing pathways, possibly mediated by higher index saddles, continues to be of interest is that of intramolecular double proton transfer (DPT). A number of studies have addressed various aspects of the mechanism of the DPT phenomenon in a variety of molecular systems. Of particular interest in such studies is the issue of whether the mechanism of DPT is concerted or stepwise. There are indications that the mechanism is sensitive to a variety of factors including temperature and the extent of coupling between the two local hydrogen stretching modes [27] . Suggestions for the two hydrogens to be quantum mechanically entangled [28, 29] and the role of nuclear quantum effects at temperatures below 100 K [30] have been made in porphycene and other structurally related molecules. Careful and rather detailed NMR experiments on DPT in the azophenine led Rumpel and Limbach [31] to observe the breakdown of the so called "rule of geometric mean". According to this rule, the rate constants k HH , k HD , and k DD for DPT reaction with single and double isotopic substitutions are expected to obey the relation

The significant breakdown of the above relation, observed in azophenine and other systems [32] , implies that the mechanism is stepwise rather than concerted. An elegant analysis of a simple model for the isotope effect in DPT can be found in the article by Albery [33] . However, the fact that the relation Eq. (1.1) is a consequence of equilibrium and traditional TST based arguments implies that dynamical effects are expected to be important in these class of molecules.

Despite many studies, there are still doubts as to whether the mechanism of DPT can be purely concerted or stepwise. In fact there are hints that the dynamics underlying DPT reactions is fairly complicated and quite distinct from single proton transfer reactions. For instance, Accardi et al. in their quantum wavepacket dynamical studies on a model two degree of freedom system showed that mechanisms can switch on an ultrafast timescales [34] . Dynamical studies in full dimensionality, both ab initio [35] and Car-Parinello [36] , also indicate the complexity of the reaction mechanism. Indeed, in the context of Diels-Alder reactions, Houk and coworkers [37] have argued that the mechanism can be described as "dynamically concerted" since the time lag between the transfer of the two protons is less than the characteristic vibrational period. The wavepacket dynamical results, according to Accardi et al., is therefore possibly a "quantization" of the notion of classical synchronicity proposed by Houk and coworker. We also refer the reader to the work [38] of Takeuchi and Tahara for an illuminating discussion of the concerted versus stepwise debate in the context of DPT in the 7-azaindole dimer.

From the above discussion it is clear that the mechanistic implications of higher index saddles even in gas phase reactions require careful dynamical considerations. In particular, the nature of IVR and consequently the possibility of emergence of new pathways and mechanisms due to the coupling between modes deserves close attention. An important issue is whether the classically predicted dynamical pathways and preferences will be respected by the quantum dynamics. For instance, can quantum tunneling lead to the switching of a stepwise isomerization pathway to a concerted one? Moreover, to date, the signatures of the index-2 saddles on the quantum eigenstates and dynamics has not been explained clearly enough. We note that although the work of Accardi et al. [34] is suggestive of a quantum manifestation of the classical synchronicity, a confirmation by performing full dimensional quantum dynamics on the molecules studied by Houk and coworkers [37] is not practical yet. In this article, we perform extensive classical and quantum dynamical studies on the model system [22] proposed by Collins et al., which is closely related to the earlier models of DPT, to address the above issues. We demonstrate a striking classical-quantum dynamical correspondence and argue that much of the mechanism of DPT at high energies can be understood from a classical phase space perspective.

To study the dynamics in the proximity of the index-2 saddle, we consider a classical two degrees of freedom Hamiltonian of the form

where (P, Q) are momentum and position variables. The zeroth-order Hamiltonian term is given by

with the symmetric double-well potential term

The coupling term

preserves the symmetry and γ is the coupling strength.

The model Hamiltonian in Eq.(2.1) has been used to describe the double proton transfer (DPT) in a number of studies. In particular, the coordinates (Q 1 , Q 2 ) are essentially the collective coordinates introduced in a detailed study of DPT by Smedarchina et al [39, 40] . The mode-mode coupling in Eq. (2.4) may appear to be a special choice. However, as elucidated by Smedarchina et al, the coupling term γQ 2 1 Q 2 2 in the collective coordinates accounts for couplings up to fourth order in the original local proton transfer coordinates. This can be seen by transforming the potential in terms of the local proton transfer coordinates q 1,2 ≡ (Q 1 ± Q 2 )/2, resulting in the potential V (q) = V 0 (q) + V coup (q) with the zeroth-order part

and the coupling term

Thus, the potential in the collective coordinates V (Q) with a single coupling term γQ 2 1 Q 2 2 accounts for the various couplings in V (q). However, note that in the context of DPT the last term in Eq. (2.6) is not allowed since symmetry considerations demand that the couplings must be symmetric in the two coordinates and sensitive to their signs, leading to the symmetry between Q 1 and Q 2 to be broken. In the present study, the symmetry between Q 1 and Q 2 is broken due to the unequal choice of the parameters corresponding to the two double well potentials (cf. Table 1 ) in V 0 (Q). An additional difference between the models has to do with the fact that in our case all four minima are degenerate while the Smedarchina model has two pairs of isoenergetic minima. In any case, our aim here is to study the correspondence between classical and quantum dynamics of the Hamiltonian model in Eq. (2.1) without necessarily being restricted to a DPT process. Nevertheless, as seen below, our model shares many of the dynamical features observed in the earlier studies on DPT.

The numerical values of potential parameters used in this study are listed in Table 1 . The values of masses are kept equal to the mass m H of a hydrogen atom. Fig. 1 represents a contour plot of the potential energy surface (PES). The various critical points of the PES can be summarized as follows:

1. Four isoenergetic minima (wells) located at

Following Collins et al, the minima are denoted as (−−), (+−), (++), and (−+) with energy V M . In this work we consider the (−−) well to be the reactant (R) and the (++) well as the product (P). However, given that the minima are isoenergetic, one can also choose any pair of minima as R and P.

2. The minima are connected by four index-1 saddles, of which the two saddles denoted SP

have energy V † 1 lower than the ones with energy V † 2 located at

and denoted by SP

3. An index-2 saddle denoted SP 2 with energy V ‡ = 0 at (Q 1 , Q 2 ) = (0, 0).

The position of the various critical points of the PES along with their energies for the specific choice of parameters in Table 1 are given in Table 2 . For future reference and for an idea of the timescales involved, we note that an approximate harmonic approximation around the (−−) well yields the vibrational time periods of about 0.06 ps (≡ 6 × 10 −14 s ≡ 60 fs) and 0.08 ps for the Q 1 and Q 2 modes respectively. 

We begin our dynamical studies by highlighting the key features of the quantum dynamics at energies above the index-2 saddle energy. Note that this is motivated by the several quantum dynamical studies that have been done earlier on similar model Hamiltonians in the context of DPT. Our purpose here is to point out the subtleties involved, as emphasized [34] by Accardi et al and others [41] , in deciphering the reaction mechanism. 

The dynamics of the system can be investigated in terms of propagating several types of initial quantum states. However, since our interest lies in exploring the correspondence between the classical and quantum dynamics, we focus here on quantum initial states that are minimum uncertainty wavepackets. Thus, we choose normalized initial states of the form

corresponding to a wavepacket centered at (P 0 , Q 0 ) with the normalization factor N j = (2Re(β j )/π) 1/4 . For the wavepacket above the position and momentum uncertainties ∆Q j and ∆P j are equal to 1/(2 β j ) and β j respectively. Thus, ∆Q j ∆P j = /2, corresponding to a minimum uncertainty state. In what follows, we set the squeezing parameter β j = 1.0, for j = 1, 2 and vary the . The mean energy associated with the initial state is evaluated as

where C αΨ ≡ α|Ψ with |α and E α being the quantum eigenstates and eigenvalues, respectively, of the full HamiltonianĤ. The time dependent Schrödinger equation for a specific initial quantum state Ψ(Q, 0)

is solved numerically using the split operator technique [42, 43] . As this is a well known technique, we provide a brief description of the method. In Eq. (3.4),K andV represent the kinetic and potential operators, respectively. Given a specific initial quantum state Ψ(Q, 0), the quantum state Ψ(Q, t) at time t can be obtained using the unitary time evolution operatorÛ (t)

We can representÛ (t) as a product of n consecutive time evolution operators over short time intervals ∆t = t/n.Û

We note that due to the non-commutativity ofK andV , we have

The basic idea of the split operator approach is to introduce a symmetrized product ofK andV

leading to reduced error, and hence faster time evolution. By using Eq. (3.8) in Eq. (3.6), we obtain

(3.9) We note that the matrix forV is diagonal in the position representation, whereas the matrix forK is diagonal in the momentum representation. Consequently, in order to determine the time-evolved state Ψ(Q, t), we use the Fast Fourier transform and its inverse to efficiently propagate the initial state. Although it is possible to use more sophisticated splitting algorithms [44] [45] [46] , Eq. 

Following Accardi et al [34] , we define different domains D in the configuration space. As shown in Fig. 2(c) , we define five domains that include the reactant (−−), product (++), "intermediates" (+−)/(−+), and the index-2 saddle regions. As noted earlier, the precise geometric definitions of the domains is not particularly important, and the resulting domain probabilities

are not significantly different for alternative definitions of D. In particular, the mechanistic information is robust to slight changes in the definitions of the domains D. In Fig. 2(b) , we show the results for the quantum domain probabilities as a function of time for a wavepacket initiated at (Q 1 , Q 2 ) with a mean energy of E wp ≈ 0.044 au. Two observations are worth noting. First, the reaction is ultrafast and the product regions are populated within a timescale of about 50 fs, which is of the order of a harmonic vibration period in the reactant well; a similar ultrashort timescale is observed for the average residence time in the intermediate and index-2 domains. Second, as noted by Accardi et al [34] , the first forward reaction is clearly concerted since the maximum probability is associated with the index-2 saddle region. However, at later times (∼ 150 fs) the forward reaction is mostly due to a sequential mechanism. In fact, even the first backward reaction from products to reactants is non-concerted. Clearly, Fig. 2(b) indicates that the role of the index-2 saddle for the overall reaction is gradually reduced with time.

Given the fact that the mean energy of initial wavepacket is much above the index-2 saddle energy, the switching of the reaction mechanism from a concerted to a sequential one is perhaps puzzling. We note here that due to the ultrashort residence times in the (+−)/(−+) and index-2 domains, it has been argued that the mechanism can be considered as effectively synchronous. Nevertheless, rationalizing the observations in Fig. 2(b) is crucial for two main reasons. Firstly, such an insight may offer interesting clues for implementing rational control strategies that exploit the presence of higher index saddles. Secondly, the dynamical switching of the mechanism might pose a challenge for TST-based rate computations in terms of fluxes across appropriately defined locally recrossing free dividing surfaces.

In order to understand the mechanism switch, in Fig. 3 we show the time evolution of the wavepacket which results in the domain probabilities shown in Fig. 2(b) . Starting from the initial time the evolution is shown in time steps of 0.012 ps until a final time of 0.168 ps. Note that this time range, as evident from Fig. 2(b) , covers the first forward and backward reaction and the second forward reaction. The various aspects of the time evolution, from the initial flux across the index-2 region to the subsequent wavepacket dispersion and relief reflections leading to the domain densities in Fig. 2(b) , are in accordance with the earlier studies. Consequently, as analyzed and argued earlier in detail [34] , the switch from a concerted to the stepwise mechanism is due to the dispersion and the subsequent relief reflections of the wavepacket from the steep repulsive wall of the PES in the product domain into the (−+)/(+−) domains.

Given the results in Fig. 2 (b) and the wavepacket dynamics shown in Fig. 3 along with the rationale provided above, it is reasonable to expect that the switching of the mechanism is of quantum origins. However, confirming such expectations requires investigating the corresponding classical dynamics. Thus, we compute the classical analog of the domain probabilities in Eq. (3.10) as

with the classical phase space density computed based on the formula [47] ρ cm (Q, P; t) = dQ dP δ[Q − Q t (Q , P )]δ[P − P t (Q , P )]ρ cm (Q , P ; 0) (3.12)

In the above formal solution of the Liouville equation, (Q t , P t ) is the classical trajectory with the initial condition (Q , P ). The classical density at t = 0 is chosen as Gaussians in the phase space with position and momentum widths consistent with the initial quantum wavepacket density |Ψ(Q, 0)| 2 . In order to compute P cm D (t) we initiate an ensemble of 20000 initial conditions sampled according to the initial density ρ cm (Q , P ; 0) and integrate their equation of motion forward in time for 300 fs (0.3 ps) using a fourth order Runge-Kutta algorithm. The resulting classical domain probabilities, shown in Fig. 2(a) , clearly indicate that the mechanism switch occurs classically as well. Moreover, apart from the expected quantitative differences, the classical and quantum results agree with regards to the timescales and the overall qualitative dynamics. This is perhaps not surprising, since the dynamics of minimum uncertainty wavepackets on short timescales are generally expected to exhibit excellent classical-quantum correspondence. Nevertheless, it behooves us to identify the origins of the mechanism switch from a classical phase space perspective.

Note that several other wavepackets, with mean energy greater than that of the index-2 saddle, exhibit similar switching dynamics. However, it is crucial to point out that there do exist initial wavepackets that undergo very different dynamics. In particular, they completely ignore the index-2 region, instead utilizing the (−+)/(+−) domains i.e., a stepwise mechanism. Such examples can be labeled as "dynamical" Murrell-Laidler (DML) cases -dynamical, since the trajectories are not on the IRC, and Murrell-Laidler since they access only the domains corresponding to the index-1 saddles. In Fig. 4 , we show one such example for DML for an initial state with a mean energyĒ ≈ 0.045 au. Again, the classical and quantum domain probabilities are in good qualitative agreement and the avoidance of the index-2 region is clear. The snapshots of the wavepacket shown in Fig. 4 indicate that the wavepacket center essentially follows the classical trajectory.

We now turn our attention to a detailed investigation of the manifolds in the classical phase space that are responsible for the observed reaction dynamics.

The model Hamiltonian in Eq. (2.1) has two degrees of freedom and hence one can use the Poincaré surface of section to visualize the global phase space structures. However, due to the multiwell potential, a given sectioning condition can only reveal part of the dynamical information. For example, the Q 1 = 0 section will show the trajectories that isomerize between the (−−) and (+−) wells but not the transitions between the (−−) and (−+) wells. Clearly, any sectioning condition will only reveal two among the four wells. Therefore, some care is required in interpreting the surface of sections.

Since in this work we are interested in the isomerization between the (−−) (reactant) and the (++) (product) wells, we choose the Q 1 = Q 2 sectioning condition. In Fig. 5 , we show the resulting Poincaré sections for three values of the total energy that are above that of the index-2 saddle and for coupling strengths ranging from γ = 0 to γ = 1.0 × 10 −4 . We show example configuration space representation of the different dynamics in Fig. 6 and note several observations at this stage. Firstly, for a fixed energy, increasing the coupling leads to a mixed regular-chaotic phase space. Secondly, one can observe several classes of regular trajectories, particularly at the higher energies. Clearly, the phase space is far from being completely chaotic even for the highest energy and coupling strength considered. In turn this observation, apart from illustrating the danger of naively assuming ergodicity at sufficiently high energies and couplings, implies that whether the mechanism is concerted or not depends crucially on the nature of the initially prepared state. As evident from the configuration space representation of the various trajectories in Fig. 6 , even at E = 0.05 au there is a coexistence of both concerted and non-concerted dynamics. . The outer black line in the plots is the boundary of energetically accessible energy surface, that is, the zero velocity curve or Hills region. In addition, we note that symmetric island regions with different color mean that they are classically disconnected.

Note that, amongst the regular trajectories in Fig. 5 , several nonlinear resonances can be observed. Such nonlinear resonances correspond to facile energy exchange between the two modes. A particularly interesting example is the prominent resonance shown in yellow in Fig. 5 whose configuration space representation (cf. Fig. 6(f) ) shows that they correspond to the DML type. In particular, the trajectory shown in Fig. 4 and the associated wavepacket evolution establish that purely non-concerted mechanism can manifest well above the index-2 saddle energy, depending on the initial state preparation. It is important to note that the DML mechanism occurs due to the coupling γ = 0. This fact emphasizes the crucial role of IVR in the reaction mechanism. We note that the various nonlinear resonances that may manifest in the phase space can be predicted based on the zeroth-order Hamiltonian. This is possible since the exact zeroth-order nonlinear frequencies can be computed for a quartic oscillator. However, as is well known, different sets of action-angle variables and frequencies are associated with energies below and above the index-1 saddles. Thus, as is relevant to the DML example, focusing on the case wherein the unperturbed energies for both the modes are greater than the respective index-1 saddle energies, the ratio of the exact unperturbed frequencies is given by the expression

.

(4.1)

In the above equation, E 1 and E 2 correspond to the zeroth-order energies in the two modes. The K(k) are complete elliptic integrals of the first kind with

As usual, for Ω 1 /Ω 2 = r/s with integers (r, s) one predicts the specific nonlinear resonances, visible in the surface of section as resonance islands. Note that the analytic solution q j ∝ cn(2K(k j )θ j /π, k j ), withθ j ≡ Ω j , allows one to determine the various possible nonlinear resonances. Indeed, using the standard Fourier representation of the Jacobi elliptic function it is possible to show that the primary resonances are of the form Ω 1 : Ω 2 = 2r : 2s. Consequently, the lowest order term is a 2 : 2 resonance which corresponds to the DML dynamics. Similar results can be derived for E 1 and E 2 being below the index-2 saddle energies, only now all resonances of the form Ω 1 : Ω 2 = r : s are possible.

For the model Hamiltonian of interest further analysis in terms of the zeroth-order actionangle variables is not straightforward since the zeroth-order Hamiltonian cannot be written down explicitly in terms of the action variables. However, numerically, the technique of time-frequency analysis is a powerful approach to determine the various frequency lockings. Moreover, the timefrequency analysis can be applied rather generally to the fully coupled system in order to analyze the non-integrable dynamics. Here we use the continuous wavelet transform approach [48] which has been employed earlier in several studies of IVR and reaction dynamics [49] . Using this approach, in Fig. 6 we show the frequency ratios for the various classes of dynamics that are possible in the model system. As expected, chaotic trajectories have frequencies that vary considerably with time whereas the regular trajectories have constant frequencies. Nevertheless, an important aspect to note from the time-frequency analysis is that the chaotic trajectories exhibit stickiness on fairly long timescales. Interestingly, from the lone example of chaotic case shown in Fig. 6 , the stickiness can occur near concerted as well as non-concerted phase space features.

In this section, we present the phase space structure governing the mechanism of concerted and sequential isomerization. As discussed in Sect. 3, the switch from concerted to sequential is accessible only when the energy is above the energy of the index-2 saddle. We have shown that using quantum wavepacket and classical density calculations and here we will present the underlying phase space structures for total energy, E > E S 2 , where E S 2 denotes the energy of the index-2 saddle.

As the energy is increased above the energy of the index-2 saddle, the unstable periodic orbits associated with the index-1 saddles with Q 2 = 0 coordinates (at energy −0.0225 au) coalesce to form one unstable periodic orbit, Γ 13 . Similarly, another unstable periodic orbit, Γ 24 , is obtained due to the coalescence of index-1 saddles with Q 1 = 0 (at energy −0.0162 au) coordinates. The unstable periodic orbits are obtained using a differential correction of guess initial conditions and numerical continuation to obtain the orbit at the desired energy. The unstable periodic orbits have associated invariant manifolds of geometry S × R 1 (cylindrical geometry or tubes), as shown in Fig. 7 . These are computed using globalization of the initial conditions displaced along the corresponding eigenvectors, which is used to generate the intersections with the surface of section (Eqn. 4.3) and the resulting homoclinic tangle. To justify these phase space structures as underlying the reaction mechanism switch, we present the results for E = 0.045 au and γ = 1 × 10 −5 au to compare with the wavepacket calculations. We note that our discussion holds as long as the stability type of the periodic orbits does not change with total energy and coupling strength. A detailed analysis of the changes in the stability and related manifestation of the reaction mechanism is beyond the scope of this article and will be the focus of future work. The methods used for computing the phase space structures are described in the supplemental material and we show the homoclinic tangles for sample values of the total energy and coupling strength. Let us define the surface of section:

where the ± indicates direction of crossing of the surface and in this article, we use the + direction. We note that a detailed analysis of phase space transport in this system requires combining the two diagonal sections with appropriate crossing direction. Due to the geometry of the unstable periodic orbits (UPOs) and the location of the surface of section, Γ 24 appears as a point (Q 1 = 0, P 1 = 0) and Γ 13 appears as two points (located near the energy boundary at Q 1 = 0) on the surface of section. In addition, the pair of invariant manifolds associated with Γ 13 and Γ 24 partition the energy surface into dynamically distinct regions. However, due to their co-existence at the same energy, their combined role in phase space transport of classical trajectories is to be expected [50] . An easier way to explain and visualize their interplay is via the intersection of the manifolds with the surface of section. The stable and unstable invariant manifolds of these UPOs form the homoclinic tangle as shown in Fig. 8 (c) and 9(c).

The above construction of the invariant manifolds, in light of the good classical-quantum correspondence observed for the domain probabilities in Fig. 2 , raises the intriguing possibility of modulating the extent to which the index-2 saddle influences the reaction dynamics. We note that the stable and unstable invariant manifolds of the UPOs form the homoclinic tangle, leading to large scale chaos on sufficiently long timescales. However, the ultrashort dynamics implies early time Fig. 7 . Phase space structures underlying the reaction mechanism switch at energies above the energy of the index-2 saddle. Shown here are unstable periodic orbits, Γ13 and Γ24, and its associated invariant manifolds in (a) and (b), respectively, and which are co-dimension 1 in the energy surface. The pair of cylindrical manifolds form the skeleton of the reaction mechanism switch. The manifolds have been globalized for t = 0.5 ps and at E = 0.045 au which corresponds to the mean energy of the wavepacket used in quantum dynamics calculations. structure in the homoclinic tangle and, related to the early work on reactive islands [51] , specific initial states may exhibit widely different mechanisms. We note, however, that there is a distinct possibility that the quantum wavepacket tunneling can compromise the classical "impenetrable" barriers. Apart from the domain probabilities, we perform a more detailed comparison between the classical and quantum dynamics in terms of computing the Husimi distribution [52] , which is a phase space distribution associated with the quantum state Ψ(Q, t)

where P 0 , Q 0 |φ is a minimum uncertainty state of the form in Eq.( 3.1) centered at (P 0 , Q 0 ). Here we calculate the projection of the Husimi distribution on to the Q 1 = Q 2 plane, in order to be consistent with the classical surface of section.

To consolidate our claim, we visualize the manifold intersections on the surface of section along with the evolution of the Husimi distribution over the same time interval for two different initial conditions. In Fig. 8 and 9 (c) we show the location of the center of the quantum wavepacket with respect to the homoclinic tangle formed by the intersection of the stable and unstable invariant manifolds associated with the UPOs. A closer look at the location of the initial wavepackets can be seen in Fig. 8 and 9 (d). It is interesting to note that despite both the initial conditions being in the homoclinic tangle and at the same energy of E = 0.045 au, the reaction mechanisms are quite different. Moreover both the initial conditions also exhibit mechanisms distinct from the example shown in Fig. 2 in terms of the extent of involvement of the index-2 saddle. Specifically, the initial wavepacket in Fig. 9 shows that the index-2 saddle starts to play a role at a much later time t ∼ 110 fs. In this case the early time mechanism is dominantly sequential or stepwise (−−) → (−+) → (++). In contrast, the dynamics in Fig. 8 clearly involves the index-2 saddle at a much earlier time. However, there is an important difference between the mechanism for the initial wavepacket in Fig. 8 (b) and the one in Fig. 2(b) . Although both do involve the index-2 saddle dominantly at early times, the former case utilizes the concerted dynamics for the (+−) → (−+) transfer and then sequentially to the product whereas in the latter case the index-2 saddle's role is predominantly to lead to the (−−) → (++) population transfer.

The snapshots of the Husimi distributions shown in the bottom panels of Fig. 8 and 9 indicates that the evolution is mainly along the "track" formed by the homoclinic tangle of the Γ 13 UPO. This indicates, apart from justifying the close agreement between the classical and quantum domain probabilities, a fairly detailed classical-quantum correspondence when it comes to the transport in the phase space. In other words, much of the observed mechanisms are justifiable from a purely classical dynamical perspective. This observation further bolsters the claim by Accardi et al [34] that the wavepacket mechanism switching is essentially a quantized form of the mechanisms observed by Black et al [37] . However, it is worth pointing out that despite this reasonably good classicalquantum correspondence, one can clearly see key differences between the classical and quantum dynamics. Firstly, as gleaned from the various domain probabilities shown, the role of the index-2 saddle is more dominant classically when compared to the quantum dynamics. Secondly, signatures of quantum tunneling can be observed both in the domain probabilities and in the time evolution of the Husimi distributions. In case of the domain probabilities the tunneling effects are, for instance, seen at the various transition times -the classical probabilities exhibit a near "step" function drop when compared to the more smoother quantum variations. In terms of the Husimi a clear indication of quantum tunneling can be seen from the snapshot in Fig. 8 at t ∼ 0.104 ps wherein lower density patches can be observed in the reactant well. Note that around this time the classical and the quantum domain probabilities are also not in a good correspondence.

Our dynamical studies on a two degrees of freedom system model for double proton transfer in molecular systems has clearly established a good classical-quantum correspondence with regards to the mechanism of the reaction. More importantly, our studies corroborate the expectation that the dynamical complexity precludes a clean separation of the mechanism into purely sequential or concerted one. We also show that the phase space of the system, despite strong couplings, continues to exhibit mixed regular-chaotic dynamics. As a consequence, the mechanism is expected to be sensitive to the specific molecule and the prepared initial state.

As we noted in the discussion in the last subsection, the detailed classical-quantum correspondence and the ultrafast nature of the dynamics may allow for modulating the influence of the index-2 saddle to the overall reaction mechanism. Although such an expectation is borne out by Figs. 8 and 9 , extracting the detailed phase space mechanisms warrant further studies. In particular, the differences in the reaction mechanism via crossing the index-2 saddle between the two initial conditions sampled from the homoclinic tangle are expected to be due to the result of the underlying structure of the lobes [53] . We suspect that this, and the reaction mechanism switch observed in Fig. 2 , can be investigated by using lobe dynamics on the section (4.3) along with the section passing through the intermediate wells. In this context a crucial question is whether experiments can prepare such exquisitely tailored wavepackets in complex systems. An equally relevant question is whether the predictions of the lobe dynamics are robust in the presence of noise and or dissipation. Further work is required to address both of the above issues.

Finally, there are indications that coupling of at least another mode to the two degree of freedom DPT model is crucial [36] . The results of the current study suggest that the classical-quantum correspondence should hold in the higher degree of freedom system as well. However, identifying the invariant manifolds and transport barriers in phase space will require techniques such as that of Lagrangian descriptors on suitably chosen sections [55] [56] [57] [58] and the recently introduced asymptotic trajectory indicator method [59] . We anticipate that a detailed study of the far richer intramolecular vibrational energy flow pathways along with the invariant manifolds that regulate the phase space transport will lead to control strategies based on exploiting the dynamics in the vicinity of the higher index saddles.

The activated complex in chemical reactions

Systematic exploration of the mechanism of chemical reactions: the global reaction route mapping (GRRM) strategy using the ADDF and AFIR methods

Visualization of dynamics effect: projection of on-the-fly trajectories to the subspace spanned by the static reaction path network

Transition states, trapped trajectories, and classical bound states embedded in the continuum

Classical transition state theory is exact if the transition state is unique

Wigner's dynamical transition state theory in phase space: classical and quantum

The geometry of reaction dynamics

Evidence for direct molecular oxygen production in CO2 photodissociation

Embedding of the Saddle Point of Index two on the PES of the Ring Opening of Cyclobutene

Conrotatory and disrotatory stationary points for the electrocyclic isomerization of cyclobutene to cis-butadiene

Symmetries of activated complexes

Limitations of the Murrell-Laidler theorem

Saddle points of index 2 on potential energy surfaces and their role in theoretical reactivity investigations

Group theoretical selection rules for the transition states of chemical reactions

Group theory and reaction mechanisms: An extension of the Mciver-Stanton rules

Nontotally symmetric trifurcation of an S N 2 reaction pathway

Toward a theory of mechanochemistry: Simple models from the very beginnings

Force-Induced Dissolution of Imaginary Mode in Mechanochemical Reaction: Dibenzophenazine Synthesis

Should the Woodward-Hoffmann Rules be Applied to Mechanochemical Reactions?

Covalent mechanochemistry: theoretical concepts and computational tools with applications to molecular nanomechanics

Bond breaking in a Morse chain under tension: Fragmentation patterns, higher index saddles, and bond healing

Index k saddles and dividing surfaces in phase space with applications to isomerization dynamics

Transition states near rank-two saddles: Correlated electron dynamics of helium

Reactivity boundaries for chemical reactions associated with higher-index and multiple saddles

Ab initio direct dynamics study of cyclopropyl radical ring-opening

Can reactions follow non-traditional second-order saddle pathways avoiding transition states?

Quantum tautomerization in porphycene and its isotopomers: Path-integral molecular dynamics simulations

Entanglement and co-tunneling of two equivalent protons in hydrogen bond pairs

Quantum entanglement and nonlocal proton transfer dynamics in dimers of formic acid and analogues

Elucidating the nuclear quantum dynamics of intramolecular double hydrogen transfer in porphycene

NMR study of kinetic HH/HD/DD isotope, solvent and solid-state effects on the double proton transfer in azophenine

Dynamic NMR study of the kinetic HH/HD/DD isotope effects on the double proton transfer in cyclic bis (p-fluorophenyl) formamidine dimers

Isotope effects in double-proton-transfer reactions

From synchronous to sequential double proton transfer: Quantum dynamics simulations for the model Porphine

Successive mechanism of double-proton transfer in formic acid dimer: A classical study

Car-Parrinello molecular dynamics study of the intramolecular vibrational mode-sensitive double proton-transfer mechanisms in porphycene

Dynamics, transition states, and timing of bond formation in Diels-Alder reactions

The Answer to Concerted Versus Step-wise controversy for the Double Proton Transfer Mechanism of 7-Azaindole Dimer in Solutiuon

Correlated double-proton transfer. I. Theory

Mechanisms of double proton transfer. Theory and applications

Conditional Born-Oppenheimer Dynamics: Quantum Dynamics Simulations for the Model Porphine

Introduction to quantum mechanics: a time-dependent perspective

Program for quantum wave-packet dynamics with timedependent potentials

Symplectic splitting operator methods for the time-dependent Schrödinger equation

Improved exponential split operator method for solving the timedependent Schrödinger equation

Efficient fourth-order split operator for solving the triatomic reactive schrodinger equation in the time-dependent wavepacket approach

Structures in classical phase space and quantum chaotic dynamics

Time-frequency analysis of classical trajectories of polyatomic molecules

Scaling Perspective on Intramolecular Vibrational Energy Flow: Analogies, Insights, and Challenges

The role of normally hyperbolic invariant manifolds (NHIMS) in the context of the phase space setting for chemical reaction dynamics

Reactive islands as essential mediators of unimolecular conformational isomerization: A dynamical study of 3-phospholene

Some formal properties of the density matrix

On the geometry of transport in phase space I. Transport in k-degree-of-freedom Hamiltonian systems, 2 k < ∞

Computational method for phase space transport with applications to lobe dynamics and rate of escape. Regul. Chaotic Dyn

Hidden geometry of ocean flows

Detecting reactive islands using Lagrangian descriptors and the relevance to transition path sampling

Lagrangian descriptors: A method for revealing phase space structures of general time dependent dynamical systems

Exploring Isomerization Dynamics on a Potential Energy Surface with an Index-2 Saddle using Lagrangian Descriptors

Identifying Reaction Pathways in Phase Space via Asymptotic Trajectories

The authors would like to acknowledge the fruitful discussions with Stephen Wiggins, the tireless support of Eleanor Machin during the COVID-19 lockdown in Bristol and IIT Kanpur High Performance Computing facility for computing resources. 

The authors declare that they have no conflicts of interest.

Priyanka Pandey performed the classical and quantum dynamical computations. Shibabrat Naik contributed in terms of the computation of the classical invariant manifolds and their intersections. All authors participated in discussing the results and were involved in writing the text of the paper.

Supplementary material with details on computation of the invariant manifolds is available for this article at and are accessible for authorized users.