1 INTRODUCTION

Transition state theory (TST) is a cornerstone for reaction dynamics [1]. The simple idea of calculating the rate of a reaction as 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 nondynamical. 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 in the early work of 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 and the higher-index saddles. 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,3\)-butadiene can occur via a disrotatory pathway, as opposed to the conrotatory pathway predicted based on the Woodward – Hoffman rules, due to the involvement of an 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 is that these rules are based on the IRC perspective, and hence nondynamical. Note, however, that even within the IRC perspective, it is possible for an index-\(1\) saddle to be linked to a higher-index saddle which can then lead to multiple products. An example of 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 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 criterion and it is quite feasible for the “forced” dynamics to prefer pathways that explore the various high-index saddle regions.

Recently there have been several 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-degree-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 last two examples being three-degree-of-freedom systems. 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 an 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). Several 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_{\rm HH},k_{\rm HD}\), and \(k_{\rm DD}\) for DPT reaction with single and double isotopic substitutions are expected to obey the relation

$$k_{\rm HD}=\sqrt{k_{\rm HH}k_{\rm DD}}\implies\frac{k_{\rm HH}}{k_{\rm HD}}=\frac{k_{\rm HD}}{k_{\rm DD}}.$$
(1.1)

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 (1.1) is a consequence of equilibrium and traditional TST-based arguments implies that dynamical effects are expected to be important in these classes of molecules.

Despite many studies, there are still doubts as to whether the mechanism of DPT can be purely concerted or stepwise. 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 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 of Takeuchi and Tahara [38] 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 the 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 have not been explained clearly. 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 proposed by Collins et al. [22], 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.

In Section 2, we describe the model Hamiltonian for the correlated proton transfer reaction in a dimensional form. In Section 3, we study the dynamics of the quantum wavepackets and in Section 4, we study the classical dynamics in the context of the switch between the stepwise and concerted mechanisms. We also illustrate the presence of the dynamical Murrell – Laidler mechanism which ignores the concerted mechanism. In Section 5, we present our analysis of the mechanism switch using the phase space structures that underlie the quantum and classical dynamics. In Section 6, we conclude with our summary and outlook on future work.

2 MODEL HAMILTONIAN

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

$$H(\mathbf{P},\mathbf{Q})=H_{0}(\mathbf{P},\mathbf{Q})+V_{coup}(\mathbf{Q}),$$
(2.1)

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

$$H_{0}(\mathbf{P},\mathbf{Q})=\sum_{j=1,2}\left[\frac{1}{2m_{\rm H}}P_{j}^{2}+V_{j}(Q_{j})\right]$$
(2.2)

with the symmetric double-well potential term

$$V_{j}(Q_{j})=-a_{j}Q_{j}^{2}+b_{j}Q_{j}^{4}.$$
(2.3)

The coupling term

$$V_{coup}(\mathbf{Q})=\gamma Q_{1}^{2}Q_{2}^{2}$$
(2.4)

preserves the symmetry and \(\gamma\) 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 \(\gamma Q_{1}^{2}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}\equiv(Q_{1}\pm Q_{2})/2\), resulting in the potential \(V({\bf q})=V_{0}({\bf q})+V_{\rm coup}(\bf q)\) with the zeroth-order part

$$V_{0}({\bf q})=-\left(\frac{a_{1}+a_{2}}{4}\right)(q_{1}^{2}+q_{2}^{2})+\left(\frac{b_{1}+b_{2}}{16}\right)(q_{1}^{4}+q_{2}^{4})$$
(2.5)

and the coupling term

$$V_{\rm coup}({\bf q})=-\left(\frac{a_{1}+a_{2}}{2}\right)q_{1}q_{2}+\left(\frac{b_{1}-b_{2}}{4}\right)(q_{1}^{3}q_{2}+q_{1}q_{2}^{3})+\left(\frac{3(b_{1}+b_{2})}{8}\right)q_{1}^{2}q_{2}^{2}.$$
(2.6)

Thus, the potential in the collective coordinates \(V(\bf Q)\) with a single coupling term \(\gamma Q_{1}^{2}Q_{2}^{2}\) accounts for the various couplings in \(V({\bf 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}({\bf 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.

Fig. 1
figure 1

Contour plot of the model potential energy surface. Blue, black, and red points denote the positions of the index-\(1\) saddles of lower energy, index-\(1\) saddles of higher energy, and the index-\(2\) saddle, respectively. Symbolic descriptions of the minima refer to the sign of the \((Q_{1},Q_{2})\) coordinates, respectively.

Full size image

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 mH of a hydrogen atom. Figure 1 represents a contour plot of the potential energy surface (PES). The various critical points of the PES can be summarized as follows:

  1. 1)

    Four isoenergetic minima (wells) located at

    $$(Q_{1},Q_{2})=\left(\pm\sqrt{\frac{2a_{1}b_{2}-\gamma a_{2}}{4b_{1}b_{2}-\gamma^{2}}},\pm\sqrt{\frac{2a_{2}b_{1}-\gamma a_{1}}{4b_{1}b_{2}-\gamma^{2}}}\right).$$
    (2.7)

    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. 2)

    The minima are connected by four index-\(1\) saddles, of which the two saddles denoted SP\({}_{1}^{(l)}\) located at

    $$(Q_{1},Q_{2})=\left(\pm\sqrt{\frac{a_{1}}{2b_{1}}},0\right)$$
    (2.8)

    have energy \(V_{1}^{\dagger}\) lower than the ones with energy \(V_{2}^{\dagger}\) located at

    $$(Q_{1},Q_{2})=\left(0,\pm\sqrt{\frac{a_{2}}{2b_{2}}}\right)$$
    (2.9)

    and denoted by SP\({}_{1}^{(h)}\).

  3. 3)

    An index-\(2\) saddle denoted SP\({}_{2}\) with energy \(V^{\ddagger}=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 (\(\equiv 6\times 10^{-14}\) s \(\equiv 60\) fs) and \(0.08\) ps for the \(Q_{1}\) and \(Q_{2}\) modes, respectively.

Table 1 Parameters used in the potential of the system
Full size table
Table 2 Critical points and the corresponding energies (in atomic units (au)) for the model potential
Full size table

3 QUANTUM WAVEPACKET DYNAMICS

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 by Accardi et al. [34] and others [41], in deciphering the reaction mechanism.

3.1 Choosing the Initial State and Time Propagation

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

$$\Psi({\bf Q},0)=\prod_{j=1,2}N_{j}\exp{\left[-\beta_{j}(Q_{j}-Q_{j_{0}})^{2}+\frac{i}{\hbar}P_{j_{0}}(Q_{j}-Q_{j_{0}})\right]}$$
(3.1)

corresponding to a wavepacket centered at (\({\bf P}_{0},{\bf Q}_{0}\)) with the normalization factor \(N_{j}=(2{Re}(\beta_{j})/\pi)^{1/4}\). For the wavepacket above, the position and momentum uncertainties \(\Delta Q_{j}\) and \(\Delta P_{j}\) are equal to \(1/(2\sqrt{\beta_{j}})\) and \(\hbar\sqrt{\beta_{j}}\), respectively. Thus, \(\Delta Q_{j}\Delta P_{j}=\hbar/2\), corresponding to a minimum uncertainty state. In what follows, we set the squeezing parameter \(\beta_{j}=1.0\), for \(j=1,2\) and vary the \(\hbar\). The mean energy associated with the initial state is evaluated as

$$\bar{E}=\langle\Psi|\hat{H}|\Psi\rangle=\sum_{\alpha}|C_{\alpha\Psi}|^{2}E_{\alpha},$$
(3.2)

where \(C_{\alpha\Psi}\equiv\langle\alpha|\Psi\rangle\) with \(|\alpha\rangle\) and \(E_{\alpha}\) are the quantum eigenstates and eigenvalues, respectively, of the full Hamiltonian \(\hat{H}\).

The time-dependent Schrödinger equation for a specific initial quantum state \(\Psi({\bf Q},0)\)

$$\displaystyle i\hbar\frac{\partial\Psi({\bf Q},t)}{\partial t}=\left[\sum_{j=1,2}\left(-\frac{\hbar^{2}}{2m_{\rm H}}\frac{\partial^{2}}{\partial Q_{j}^{2}}+V_{j}(Q_{j})\right)+\gamma Q_{1}^{2}Q_{2}^{2}\right]\Psi({\bf Q},t)$$
(3.3)
$$\displaystyle\equiv[\hat{K}+\hat{V}]\Psi({\bf Q},t)$$
(3.4)

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), \(\hat{K}\) and \(\hat{V}\) represent the kinetic and potential operators, respectively. Given a specific initial quantum state \(\Psi({\bf Q},0)\), the quantum state \(\Psi({\bf Q},t)\) at time \(t\) can be obtained using the unitary time evolution operator \(\hat{U}(t)\)

$$\Psi({\bf Q},t)=\hat{U}(t)\Psi({\bf Q},0)\equiv\exp[-i\hat{H}t/\hbar]\Psi({\bf Q},0).$$
(3.5)

We can represent \(\hat{U}(t)\) as a product of \(n\) consecutive time evolution operators over short time intervals \(\Delta{t}=t/n\):

$$\hat{U}(t)=\underbrace{e^{-i\hat{H}\Delta{t}/\hbar}e^{-i\hat{H}\Delta{t}/\hbar}\ldots{e^{-i\hat{H}\Delta{t}/\hbar}}}_{\text{n times}}.$$
(3.6)

We note that, due to the noncommutativity of \(\hat{K}\) and \(\hat{V}\), we have

$${e^{-i[\hat{K}+\hat{V}]\Delta{t}/\hbar}}\approx{e^{-i\hat{K}\Delta{t}/\hbar}}{e^{-i\hat{V}\Delta{t}/\hbar}}+O(\Delta{t}^{2}).$$
(3.7)

The basic idea of the split operator approach is to introduce a symmetrized product of \(\hat{K}\) and \(\hat{V}\)

$${e^{-i\hat{H}\Delta{t}/\hbar}}\approx{e^{-i\hat{V}\Delta{t}/2\hbar}}{e^{-i\hat{K}\Delta{t}/\hbar}}{e^{-i\hat{V}\Delta{t}/2\hbar}}+O(\Delta{t}^{3})$$
(3.8)

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

$${e^{-i\hat{H}t/\hbar}}={e^{-i\hat{V}\Delta{t}/2\hbar}}\underbrace{{e^{-i\hat{K}\Delta{t}/\hbar}}{e^{-i\hat{V}\Delta{t}/\hbar}}{e^{-i\hat{K}\Delta{t}/\hbar}}{e^{-i\hat{V}\Delta{t}/\hbar}}\ldots{e^{-i\hat{K}\Delta{t}/\hbar}}{e^{-i\hat{V}\Delta{t}/\hbar}}}_{\text{n-1 times}}{e^{-i\hat{K}\Delta{t}/\hbar}}{e^{-i\hat{V}\Delta{t}/2\hbar}}.$$
(3.9)

We note that the matrix for \(\hat{V}\) is diagonal in the position representation, whereas the matrix for \(\hat{K}\) is diagonal in the momentum representation. Consequently, in order to determine the time-evolved state \(\Psi({\bf 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. (3.8) is sufficient for the results reported in this study due to the ultrafast nature of the reaction dynamics. For the calculations reported here, we chose a \(512\times 512\) spatial grid. The range for both sets of the spatial and momentum coordinates is \([-8.0,8.0]\) and \([-20.0,20.0]\).

3.2 Switch from Concerted to Nonconcerted Mechanism

Following Accardi et al. [34], we define different domains \(D\) in the configuration space. As shown in Fig. 2c, 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 are not particularly important, and the resulting domain probabilities

$$P_{D}(t)\equiv\int_{D}dQ_{1}dQ_{2}|\Psi({\bf Q},t)|^{2}$$
(3.10)

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\).

Fig. 2
figure 2

(b) Quantum domain probabilities (\(\hbar=1.0\)) as a function of time for a wavepacket with mean energy of \(\approx 0.044\) au. Reactant well \((--)\) (black), product well \((++)\) (red), index-\(2\) region, defined by the square region in (c) (magenta filled curve), the \((+-)\) well (green), and the \((-+)\) well (blue). Note the switch in the mechanism from a concerted to a sequential one at \(t\sim 70\) fs (\(\equiv 0.07\) ps). (a) The corresponding classical probabilities, showing a similar mechanism switch. (c) Definition of the various domains in the configuration space \((Q_{1},Q_{2})\), adapted [34] from an earlier work. The center of the wavepacket is at \((Q_{1},P_{1},Q_{2},P_{2})=(-2.50,9.57,-4.37,14.19)\).

Full size image

In Fig. 2b, 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} \approx 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 (\(\sim 150\) fs), the forward reaction is mostly due to a sequential mechanism. In fact, even the first backward reaction from products to reactants is nonconcerted. Clearly, Fig. 2b 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 the 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. 2b 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.

Fig. 3
figure 3

Configuration space snapshots of the time evolving wavepacket corresponding to Fig. 2b. The top column left panel corresponds to time \(t=0\) ps and the right column bottom panel corresponds to \(t=0.168\) ps. The rest of the panels are shown at intervals of \(0.012\) ps.

Full size image

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. 2b. 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. 2b, 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. 2b, are in agreement with the earlier studies. Consequently, as analyzed and argued earlier in detail in [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.

Fig. 4
figure 4

Initial quantum wavepacket with mean energy \(\bar{E}\approx 0.045\) au exhibiting Murrell – Laidler type dynamics. The associated domain probabilities (a) classical and (b) quantum are shown. Configuration space snapshots of the wavepacket at indicated times are also shown. For reference, the classical trajectory corresponding to the center of the wavepacket is shown as a bold black line. The center of the wavepacket is at \((Q_{1},P_{1},Q_{2},P_{2})=(-3.50,-12.07,-3.77,12.97)\).

Full size image

3.3 Is the Mechanism Switch a Quantum Effect?

Given the results in Fig. 2b 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

$$P_{D}^{cm}(t)=\int_{D}d{\bf Q}\int d{\bf P}\rho_{cm}({\bf Q},{\bf P};t)$$
(3.11)

with the classical phase space density computed based on the formula [47]

$$\rho_{cm}({\bf Q},{\bf P};t)=\int d{\bf Q}^{\prime}d{\bf P}^{\prime}\delta[{\bf Q}-{\bf Q}_{t}({\bf Q}^{\prime},{\bf P}^{\prime})]\delta[{\bf P}-{\bf P}_{t}({\bf Q}^{\prime},{\bf P}^{\prime})]\rho_{cm}({\bf Q}^{\prime},{\bf P}^{\prime};0).$$
(3.12)

In the above formal solution of the Liouville equation, \(({\bf Q}_{t},{\bf P}_{t})\) is the classical trajectory with the initial condition \(({\bf Q}^{\prime},{\bf P}^{\prime})\). 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 \(|\Psi({\bf Q},0)|^{2}\). In order to compute \(P_{D}^{cm}(t)\), we initiate an ensemble of \(20000\) initial conditions sampled according to the initial density \(\rho_{cm}({\bf Q}^{\prime},{\bf P}^{\prime};0)\) and integrate their equation of motion forward in time for \(300\) fs (\(0.3\) ps) using a fourth-order Runge – Kutta method.

The resulting classical domain probabilities, shown in Fig. 2a, 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 of 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 \(\bar{E}\approx 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.

4 CLASSICAL DYNAMICS: DYNAMICAL MURRELL – LAIDLER MECHANISMAND THE COMPETITION BETWEEN STEPWISE AND CONCERTED MECHANISMS

4.1 A Survey of the Classical Phase Space

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 multi-well 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 given by \(U_{1}^{+}\) in Eq. (4.3). 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 \(\gamma=0\) to \(\gamma=1.0\times 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 nonconcerted dynamics.

Fig. 5
figure 5

Poincaré surface of sections along \(Q_{1}=Q_{2}\) for various coupling values (indicated) at a total energy of (a) \(0.001\) au (top row), (b) \(0.01\) au (middle row), (c) \(0.05\) au (bottom row). The outer black line in the plots is the boundary of energetically accessible energy surface, that is, the zero velocity curve or Hill’s region. In addition, we note that symmetric island regions with different color mean that they are classically disconnected.

Full size image

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. 6f) shows that they correspond to the DML type. In particular, the trajectory shown in Fig. 4 and the associated wavepacket evolution establish that purely nonconcerted 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 \(\gamma\neq 0\). This fact emphasizes the crucial role of IVR in the reaction mechanism.

Fig. 6
figure 6

(a) Time-frequency analysis for the trajectories at a total energy of \(0.05\) au and a coupling strength of \(\gamma=1\times 10^{-5}\). Configuration space representation of (b)–(e) concerted and (f)–(h) nonconcerted trajectories. Note that the color code is the same as indicated in the corresponding surface of section shown in Fig. 5 and the teal colored circle represents the initial state of the trajectory.

Full size image

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 where 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

$$\frac{\Omega_{1}}{\Omega_{2}}=\left[\frac{a_{1}^{2}+4b_{1}E_{1}}{a_{2}^{2}+4b_{2}E_{2}}\right]^{1/4}\frac{\mathbb{K}(k_{1})}{\mathbb{K}(k_{2})}.$$
(4.1)

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

$$k_{j}^{2}(E_{j})=\frac{1}{2}\left[1+\frac{a_{j}}{\sqrt{a_{j}^{2}+4b_{j}E_{j}}}\right].$$
(4.2)

As usual, for \(\Omega_{1}/\Omega_{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}\propto\mathop{\rm cn}(2\mathbb{K}(k_{j})\theta_{j}/\pi,k_{j})\), with \(\dot{\theta}_{j}\equiv\Omega_{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 \(\Omega_{1}:\Omega_{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}\), which are below the index-\(2\) saddle energies, only now all resonances of the form \(\Omega_{1}:\Omega_{2}=r^{\prime}:s^{\prime}\) are possible.

For the model Hamiltonian of interest, further analysis in terms of the zeroth-order action-angle 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 time-frequency analysis can be applied rather generally to the fully coupled system in order to analyze the nonintegrable 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 the chaotic case shown in Fig. 6, the stickiness can occur near concerted as well as nonconcerted phase space features.

4.2 Unstable Periodic Orbits and Invariant Manifolds that Regulate the Reaction Dynamics

In this section, we present the phase space structure governing the mechanism of concerted and sequential isomerization. As discussed in Section 3, the switch from concerted to sequential is accessible only when the energy is above the energy of the index-2 saddle. In the previous section, we have shown the classical nature of the mechanism switch using quantum wavepacket and classical phase space density dynamics. Here we will present the underlying phase space structures that possibly regulate the mechanism switch for total energy, \(E>E_{\rm S_{2}}\), where \(E_{\rm 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, \(\Gamma_{13}\). Similarly, another unstable periodic orbit, \(\Gamma_{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 \(\mathbb{S}\times\mathbb{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 (Eq. (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 \(\gamma=1\times 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 (see [60] for more details on a similar PES). We show the homoclinic tangles for sample values of the total energy and coupling strength. Let us define the surface of section:

$$U_{1}^{\pm}=\{(Q_{1},P_{1})|Q_{2}-Q_{1}=0,-{\rm sign}(P_{2}-P_{1})=\pm 1\},$$
(4.3)

where the \(\pm\) 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.

Fig. 7
figure 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, \(\Gamma_{13}\) and \(\Gamma_{24}\), and its associated invariant manifolds in (a) and (b), respectively, and which are codimension 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.

Full size image

Due to the geometry of the unstable periodic orbits (UPOs) and the location of the surface of section, \(\Gamma_{24}\) appears as a point (\(Q_{1}=0,P_{1}=0\)) and \(\Gamma_{13}\) appears as a point (located near the energy boundary at \(Q_{1}=0\)) on the surface of section. In addition, the pair of invariant manifolds associated with \(\Gamma_{13}\) and \(\Gamma_{24}\) partition the energy surface into dynamically distinct regions. However, due to their coexistence 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. 8c and 9c.

Fig. 8
figure 8

(a) Classical and (b) quantum (\(\hbar=1.0\)) domain probabilities for an initial condition chosen on the homoclinic tangle. (c) Location of the Husimi distribution at \(t=0\) superposed on the stable (blue) and unstable (magenta) manifolds of the unstable periodic orbit \(\Gamma_{24}\). (d) A zoomed version of (c) indicating the precise location of the initial wavepacket. The bottom row shows four different time snapshots of the evolving Husimi distributions. The center of the wavepacket is at \((Q_{1},P_{1},Q_{2},P_{2})=(-5.00,-5.00,-5.00,-14.62)\).

Full size image

4.3 Modulating the Role of the Index-\(2\) Saddle

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 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 \(\Psi({\bf Q},t)\)

$$\rho_{H}({\bf P}_{0},{\bf Q}_{0};t)=\frac{1}{(2\pi\hbar)^{2}}|\langle\phi|\Psi(t)\rangle|^{2},$$
(4.4)

where \(\langle{\bf P}_{0},{\bf Q}_{0}|\phi\rangle\) is a minimum uncertainty state of the form in Eq. (3.1) centered at \(({\bf P}_{0},{\bf Q}_{0})\). Here we calculate the projection of the Husimi distribution onto the \(Q_{1}=Q_{2}\) plane in order for it to be consistent with the classical surface of section.

Fig. 9
figure 9

(a) and (b) Classical and the quantum (\(\hbar=1.0\)) domain probabilities, respectively, for another initial condition chosen on the homoclinic tangle. (c) Location of the Husimi distribution at \(t=0\) superposed on the stable (blue) and unstable (magenta) manifolds of the unstable periodic orbit \(\Gamma_{24}\). (d) A zoomed version of (c) indicating the precise location of the initial wavepacket. The bottom row shows four different time snapshots of the evolving Husimi distributions. The center of the wavepacket is at \((Q_{1},P_{1},Q_{2},P_{2})=(-3.00,7.50,-3.00,14.84)\).

Full size image

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 Figs. 8 and 9c, 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 Figs. 8 and 9d. 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 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\sim 110\) fs. In this case, the early time mechanism is dominantly sequential or stepwise \((--)\rightarrow(-+)\rightarrow(++)\). 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. 8b and the one in Fig. 2b. Although both do involve the index-\(2\) saddle dominantly at early times, the former case utilizes the concerted dynamics for the \((+-)\rightarrow(-+)\) transfer and then sequentially to the product, whereas in the latter case, the index-\(2\) saddle’s role is predominantly to lead to the \((--)\rightarrow(++)\) population transfer.

The snapshots of the Husimi distributions shown in the bottom panels of Figs. 8 and 9 indicate that the evolution is mainly along the “track” formed by the homoclinic tangle of the \(\Gamma_{24}\) 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 mechanism is 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 classical-quantum 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 the 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 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\sim 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.

5 CONCLUSIONS AND OUTLOOK

Our dynamical studies on a two-degree-of-freedom system model for double proton transfer in molecular systems have 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 a 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 on the overall reaction mechanism. Although such an expectation is borne out by Figs. 8 and 9, extracting the detailed phase space mechanisms warrants 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 studying lobe dynamics on the section passing through the intermediate wells along with the results in Section 4.3. 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 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.