Search
N-DoF System-Bath Model 2

$n$-Degree-of-Freedom System-Bath Model

Introduction and Development of the Problem

(REMOVE LD DESCRIPTION)

Isomerization is one of the most prevalent reaction when studying atmospheric, medical, and industrial chemistry [^1], and thus have garnered interest from both theoretical and applied chemistry [1], [2], [3], [4]. From an applied perspective, the influence of various solvents on the rate constant of conformational isomerization has been pursued for specific molecules [5], [6], [7], [8], [9]. From a theoretical perspective, as is adopted in this chapter, it is instructive to formulate a systematic approach for identifying reactive trajectories in the spirit of transition path sampling [10] and reactive islands [11], [12], [13], [14].

Modeling reaction in solvents

Following the derivation of the Langevin equations [15] and its generalized system dynamics in a heat bath [16], the Kramers turnover of reaction rates (escape rates from a potential well) with increasing viscosity of the solvent (intensity of noise) has been used to obtain reaction rate expressions in gas and condensed phase reactions by Grote, Hynes, Pollak, Grabert, Hänggi [17], [18]. The generalized Langevin equation of motion for a particle trapped in a one-dimensional well with a barrier height $\mathcal{V}^{\ddagger}$ and coupled to a medium such as gas or liquid is modeled by a bath with a viscosity parameter. This framework has received much attention in the literature where the system dynamics can be obtained explicitly. The Langevin dynamics represents the motion of the particle in a bath which is modeled using a large number of harmonic oscillators. In this set-up, the bath coordinates are coupled with the system coordinates, and the Hamiltonian is simply a sum of system and bath components. However, the system dynamics of a two or more degrees of freedom coupled with bath modes has not received a global analysis from a dynamical systems perspective of reactions. In this direction, the first step would be to consider a two degrees of freedom reaction with well-understood chemical observables such as reaction coordinate, and where the harmonic bath modes can be coupled to represent the reaction in a condensed phase. In this work, we will adopt the two degrees of freedom isomerization model of De Leon and Berne who have studied the chemical reaction dynamics extensively using a dynamical systems perspective [19], [14], [20], [21], [22]. This model of isomerization describes the conformational change by the motion of an internal angle where the isomers are represented by the well in a double potential well separated by a barrier. The two degrees of freedom in the model correspond to the bond that undergoes structural change and to the bond that breaks above dissociation energy. Typically, this isomerization barrier height is lower than the dissociation energy of the molecule, and the activation energy for the isomerization is imparted by molecular collisions or photoexcitation.

Effect of solvent on reaction rates and role of phase space structures

Traditionally, the construction of a dividing surface (DS) was focused on critical points of the potential energy surface (PES), that is, in the configuration space describing the molecular system [23], [24]. Critical points on the PES do have significance in phase space; they are the equilibrium points for zero momentum. But they continue to have influence for nonzero momentum for a range of energies above the energy of the equilibrium point. The construction of a DS separating the phase space into two parts, reactants and products, has been a focus from the dynamical systems point of view in recent years [25], [26], [27]. In phase space, that is for nonzero momentum, the role of the saddle point is played by an invariant manifold of saddle stability type, the normally hyperbolic invariant manifold (NHIM) [28], [29]. In order to fully appreciate the NHIM and its role in reaction rate theory, it is useful to begin with a precursor concept the periodic orbit dividing surface or PODS. For systems with two DoF described by a natural Hamiltonian, kinetic plus potential energy, the problem of constructing the DS in phase space was solved during the 1970s by McLafferty, Pechukas and Pollak [30], [31], [32], [33]. They demonstrated that the DS at a specific energy is related to an invariant phase space structure, an unstable periodic orbit (UPO) which defines (it is the boundary of) the bottleneck in phase space through which the reaction occurs. The DS which intersects trajectories evolving from reactants to products can then be shown to have the geometry of a hemisphere in phase space whose boundary is the UPO. The same construction can be carried out for a DS intersecting trajectories crossing from products to reactants and these two hemispheres form a sphere for which the UPO is the equator. Generalisation of this construction of DS to high dimensional systems has been a central question in reaction dynamics and has only received a satisfactory answer in recent years [25], [26], [27]. The key difficulty being the high dimensional analogue of the UPO used in the two DoF system for the construction of the DS and which is resolved by considering the NHIM, which has the appropriate dimensionality for anchoring the dividing surface in high dimensional phase space [28]. Normal hyperbolicity of these invariant manifolds means that their stability, in a precise sense, is of saddle type in the transverse direction, which implies that they possess stable and unstable invariant manifolds that are impenetrable barriers and mediate reactive trajectories in phase space. These invariant manifolds of the NHIM are structurally stable, that is, stable under perturbation [29]. For two DoF systems, the NHIM is an unstable PO, and for an $N > 2$ DoF system at a fixed energy, the NHIM has the topology of a $(2N-3)$-dimensional sphere and is the equator of a $(2N-2)$-dimensional sphere which is the DS. This DS can then be used to divide the $(2N-1)$-dimensional energy surface into two parts, reactants and products [34], [35], [23], [24], [36]. An elementary description of the role of the NHIM in reaction dynamics is given in [37] along with description of their geometry using quadratic normal form Hamiltonians. Fundamental theorems assure the existence of the phase space structures  NHIM and its invariant manifolds for a range of energies above that of the saddle [29]. However, the precise extent of this range, as well as the nature and consequences of any bifurcations of the phase space structures that might occur as energy is increased, is not known and is a topic of continuing research[38], [39], [40], [41], [42], [43].

Thus, calculation of reaction rate (or flux) based on the geometry of phase space structures requires identifying trajectories that start in the reactant well, cross the dividing surface constructed from the NHIM, and reach the product well. This dividing surface has been shown to be the appropriate (locally no-recrossing) surface that reactive trajectories must cross since the calculated reaction rates do not need correction due to recrossings [27]. This construction is in contrast to the "standard" transition-state-theory (TST) for constructing the dividing surface which is only exact in gas phase unimolecular reactions and when ergodicity of trajectories in the phase space holds [18]. As is now established, the no recrossing (locally) property of a dividing surface is a contribution of the phase space perspective of chemical reactions [26]. While the standard TST relies on recrossing free surface for calculating reaction flux, a dividing surface constructed in the configuration space violates this condition in the case of solvent, and the TST based reaction rate is not exact. This violation of the recrossing property when DS is constructed in the configuration space of a reaction in a high viscosity solvent also follows from the Kramers' diffusion model, Langevin equation, of chemical reactions [18]. Thus, finding the reactive trajectories, and the changes in the DS and NHIM due to a solvent presents a worthwhile step towards development of the role of phase space structures in reaction dynamics.

The geometry of unimolecular reactions dynamics has been developed using a $N$ degrees of freedom Hamiltonian where the coordinates represent intermolecular bonds. As a natural step in studying unimolecular reaction dynamics in solvents, we adopt a model where the reaction coordinates (modeled as a system Hamiltonian) are coupled with a set of harmonic bath modes (modeled as a bath Hamiltonian) [44], [45], [46], [47]. This is with the intention of parametrizing the effects of a bath (or solvent) on the system dynamics. This formulation of coupling harmonic bath modes with system dynamics also serves as a preliminary step in assessing the capabilities of a trajectory diagnostic called Lagrangian descriptors (LDs) [48] in realistic (high dimensional) chemical systems. This system-bath model will also serve as a test bed for illustrating the use of LDs in directly computing the chemical observables. But before doing that, we would like to identify the reactive trajectories and develop a systematic approach for transition path sampling (in the sense of rare event sampling) for reactions in solvent. Thus, can we identify the qualitative changes in the isomerization of a molecule in the presence of a solvent? In particular, we would like to understand the influence of the coupling strength, the solvent's viscosity, and the number of bath modes on the reactive trajectories. We will answer this using LDs which have been shown recently to be of use in detecting phase space structures that mediate reactive trajectories in dissipative, time-dependent models of chemical reactions [49], [50], and for transition path sampling in two degrees of freedom models of chemical reactions [51]. In this article, we will present an approach for sampling reactive trajectories and identifying reactive islands in high dimensional phase space of an isomerization reaction (system) in a solvent (bath coupled to each of the system configuration coordinates) which is following our work on quadratic normal form Hamiltonian systems [52].

Two DoF de Leon-Berne system coupled with bath modes

In this section we describe the system-bath model for a system with two degrees-of-freedom (DoF). While system-bath models have received a great deal of attention in the chemistry and physics community, models with systems having more than 1 DoF have received much less attention.

We will consider the 2 degree-of-freedom Hamiltonian [eqn:ham_db]{reference-type="eqref" reference="eqn:ham_db"} coupled with harmonic oscillators for a system-bath model of the form [44], [45], [46]:

\begin{align} \mathcal{H}(x,y,x_j,y_j,p_x,p_y,p_{x_j},p_{y_j}) = \underbrace{\frac{p_x^2}{2 m_s} + \frac{p_y^2}{2 m_s} + V_{DB}(x,y)}_{\text{System Hamiltonian}} + & \underbrace{\sum_{j=1}^{N_B} \frac{1}{2} \left[ \frac{p_{x_j}^2}{m_j} + \left(\omega_j x_j -\frac{c_{x,j} x}{\omega_j} \right)^2 \right]}_{\text{Coupling of the bath to $x$}} \\ & + \underbrace{\sum_{j=1}^{N_B} \frac{1}{2} \left[ \frac{p_{y_j}^2}{m_j} + \left(\omega_j y_j - \frac{c_{y,j} y}{\omega_j} \right)^2 \right]}_{\text{Coupling of the bath to $y$}} \\ %& = \mathcal{H}_{\text{system}} + \mathcal{H}_{\text{bath coupled to $x$}} + \mathcal{H}_{\text{bath coupled to $y$}} \label{eqn:sb_ham} \end{align}

where $x$ and $y$ denote the configuration space coordinates of the system, $p_x$ and $p_y$ are the associated conjugate momenta, $p_{x_j}$, $x_j$ denote the $j^{\rm th}$ bath phase space coordinates associated with the system configuration space variable $x$, and $p_{y_j}$, $y_j$ denote the $j^{\rm th}$ bath phase space coordinates associated with the system configuration space variable $y$. We assume that the frequencies, $\omega_j$, are the same for each bath, and the coupling constants for each configuration space variable to the bath are given by $c_{x,j}$ and $c_{y,j}$. We explicitly carry out the discretization that gives us the coupling constants $c_{x,j}$ and $c_{y,j}$ and the frequencies $\omega_j$.

Discretization of the Spectral Density: Derivation of the Parameters of the Bath {#app:discrete}

The coupling of the bath of harmonic oscillators to the configuration space coordinates is described by a spectral density:

$$\begin{aligned} J_x (\omega) = \frac{\pi}{2} \sum_{i=1}^{N_B} \frac{c_{x,i}^2}{\omega_i} \delta (\omega - \omega_i), \label{x-disc_spec} \\ J_y (\omega) = \frac{\pi}{2} \sum_{i=1}^{N_B} \frac{c_{y,i}^2}{\omega_i} \delta (\omega - \omega_i), \label{y-disc_spect}\end{aligned}$$

and these result from the discretization of a continuous Ohmic (linear) form with an exponential cutoff:

$$\begin{aligned} \bar{J}_x (\omega) = \eta_x \omega e^{-\frac{\omega}{\omega_{c}}} , \label{z_1-cont_spec} \\ \bar{J}_y (\omega) = \eta_y \omega e^{-\frac{\omega}{\omega_{c}}} , \label{z_2-cont_spect}\end{aligned}$$

the discretization that gives us the coupling constants $c_{x,j}$, $c_{y,j}$, and the frequencies $\omega_j$, and are given by $$\omega_j = -\omega_c \log \left(\frac{j-\frac{1}{2}}{N_B} \right), \quad j = 1, \ldots, N_B. \label{disc_freq_1}$$ and $$\begin{aligned} c_{x,j} = \sqrt{\frac{2 \eta_{x} \omega_{c}}{\pi N_B}} \, \omega_j, \label{coupling_x} \\ c_{y,j} = \sqrt{\frac{2 \eta_{y} \omega_{c}}{\pi N_B}} \, \omega_j, \quad j = 1, \ldots, N_B. \label{coupling_y}\end{aligned}$$

In this appendix we describe a scheme for discretizing the continuous spectral density which was given in [47]. In the following we will drop the subscripts $z_1$ and $z_2$ on the various quantities for the sake of a simpler notation since we will follow the same discretization procedure for each spectral density. The subscripts can then be added back afterwards.

Re-establishing the notation, the discrete spectral density function is given by:

$$J(\omega) = \frac{\pi}{2} \sum_{i=1}^{n_b} \frac{c_i^2}{\omega_i} \delta (\omega - \omega_i), \label{sf_disc}$$

and the continuous spectral density is given by:

$$\bar{J} (\omega) = \eta \omega e^{-\frac{\omega}{\omega_c}}. \label{sf_cont}$$

Discretization is obtained by carrying out the following steps.

  1. Require

    $$\int_{0}^{\infty} J(\omega) F(\omega) d \omega \approx

    \int_{0}^{\infty} \bar{J} (\omega) F(\omega) d \omega,
\label{quad_1}$$

    for any integrable function $F(\omega)$

  2. Substitute ([sf_disc]{reference-type="ref" reference="sf_disc"}) into the left-hand side of ([quad_1]{reference-type="ref" reference="quad_1"}) to obtain:

    $$\frac{\pi}{2} \sum_{i=1}^{n_b} \frac{c_i^2}{\omega_i} F(\omega_i).

    \label{quad_2}$$

  3. Approximate the right-hand side of ([quad_1]{reference-type="ref" reference="quad_1"}) by an appropriate quadrature. This is the step that we will now carry out in detail.

The quadrature recommended in [47] is the midpoint rule, after changing variables in the integral to $x = e^{-\frac{\omega}{\omega_c}}$. The reason that they give for choosing this quadrature is that it gives a uniform distribution of grid points in the unit interval $0 < x_i < 1$ and therefore a logarithmic distribution of bath frequencies $\omega_i = -\omega_c \log x_i$. They claim that such a distribution of bath frequencies is appropriate on physical grounds for an exponentially decaying density of bath states.

The change of variables gives:

$$\begin{aligned} x = e^{-\frac{\omega}{\omega_c}} \Rightarrow \log x = -\frac{\omega}{\omega_c} \Rightarrow \omega =-\omega_c \log x, \label{cv_1} \\ dx = -\frac{1}{\omega_c} e^{-\frac{\omega}{\omega_c}} d \omega = -\frac{1}{\omega_c} x d \omega.\label{cv_2}\end{aligned}$$

Then ([sf_cont]{reference-type="ref" reference="sf_cont"}) becomes

$$\bar{J} (\omega) = -\eta x \omega_c \log x, \label{sf_cont_cv}$$

and using this expression, and the change of variables given in ([cv_1]{reference-type="ref" reference="cv_1"}) and ([cv_2]{reference-type="ref" reference="cv_2"}), the right-hand side of ([quad_1]{reference-type="ref" reference="quad_1"}) becomes:

$$-\int_{0}^{1} \eta \omega_c^2 ( \log x) F(\omega (x)) dx. \label{quad_3}$$

We discretize this integral using the midpoint rule. We partition the unit interval into $n_b$ intervals of length $\frac{1}{n_b}$ and evaluate the integrand at the midpoint, $x_i$, of each sub-interval:

$$x_i = \frac{i-\frac{1}{2}}{n_b}, \, i=1, \ldots , n_b, \label{midpoint}$$

and obtain:

$$-\int_{0}^{1} \eta \omega_c^2 ( \log x) F(\omega (x)) dx \approx -\sum_{i=1}^{n_b} \eta \omega_c^2 \log x_i F(\omega(x_i)) \frac{i}{n_b}. \label{disc_quad_1}$$

Equating each term in the sum ([quad_2]{reference-type="ref" reference="quad_2"}) to each term in the sum ([disc_quad_1]{reference-type="ref" reference="disc_quad_1"}) gives:

$$\frac{\pi}{2} \frac{c_i^2}{\omega_i} = -\frac{1}{n_b} \eta \omega_c^2 \log x_i = \frac{1}{n_b} \eta \omega_c \omega_i,$$

which gives:

$$c_i =\sqrt{\frac{2 \eta \omega_c}{\pi n_b}} \, \omega_i, \quad i=1, \ldots, n_b. \label{coupling}$$

and from ([cv_1]{reference-type="ref" reference="cv_1"}) and ([midpoint]{reference-type="ref" reference="midpoint"}) we see that:

$$\omega_i = -\omega_c \log \left(\frac{i-\frac{1}{2}}{n_b} \right), \quad i=1, \ldots, n_b. \label{disc_freq_2}$$

[47] choose $\omega_c = 500 {\rm cm}^{-1}$.

Thus, these quantities are given by

\begin{equation} \omega_j = -\omega_c \log \left(\frac{j-\frac{1}{2}}{N_B} \right), \quad j=1, \ldots, N_B \label{eqn:disc_freq_1} \end{equation}

and

\begin{align} c_{x,j} = \sqrt{\frac{2 \eta_{x} \omega_{c}}{\pi N_B}} \, \omega_j, \; c_{y,j} = \sqrt{\frac{2 \eta_{y} \omega_{c}}{\pi N_B}} \,\omega_j, \quad j=1, \ldots, N_B. \label{eqn:coupling_coeff} \end{align}

We note here that the baths are coupled to each other through the system dynamics. Hamilton's equations for the system-bath dynamics, using Eqn. [eqn:sys_pot]{reference-type="eqref" reference="eqn:sys_pot"} and Eqn. [eqn:sb_ham]{reference-type="eqref" reference="eqn:sb_ham"}, are given by

\begin{align} \dot{x} = & \frac{\partial H}{\partial p_x} = \frac{p_x}{m_s}, \\ \dot{y} = & \frac{\partial H}{\partial p_y} = \frac{p_y}{m_s}, \\ \dot{x}_j = & \frac{\partial H}{\partial p_{x_j}} = \frac{p_{x_j}}{m_j}, \qquad j=1, \ldots, N_B \\ \dot{y}_j = & \frac{\partial H}{\partial p_{y_j}} = \frac{p_{y_j}}{m_j}, \qquad j=1, \ldots, N_B, \\ \dot{p}_x = &-\frac{\partial H}{\partial x} = 2 D_x \lambda \exp(-\lambda x) (\exp(-\lambda x) - 1) + \\ & \qquad \qquad \dfrac{\mathcal{V}^{\ddagger}}{y_w^4}\zeta \lambda y^2(y^2 - 2y_w^2)\exp(-\zeta \lambda x) + \sum_{j=1}^{N_B} \frac{c_{x,j}}{\omega_j} \left(\omega_j x_j - \frac{c_{x,j} x}{\omega_j} \right), \\ \dot{p}_y = &- \frac{\partial H}{\partial y} = -4 \dfrac{\mathcal{V}^{\ddagger}}{y_w^4}y(y^2 - y_w^2)\exp(- \zeta \lambda x) + \sum_{j=1}^{N_B} \frac{c_{y,j}}{\omega_j} \left(\omega_j y_j - \frac{c_{y,j} y}{\omega_j} \right), \\ \dot{p}_{x_j} = &- \frac{\partial H}{\partial {x_j}} = - \omega_j \left(\omega_j x_j - \frac{c_{x, j} x}{\omega_j} \right) , \qquad j=1, \ldots, N_B \\ \dot{p}_{y_j} = & -\frac{\partial H}{\partial {y_j}} = - \omega_j \left(\omega_j y_j - \frac{c_{y, j} y}{\omega_j} \right) , \qquad j=1, \ldots, N_B, \\ \label{eqn:Hameq_sys_bath} \end{align}

where the bath frequencies $\omega_j$ and coupling coefficients $c_{x,j},c_{y,j}$ are given by Eqn. [eqn:disc_freq_1]{reference-type="eqref" reference="eqn:disc_freq_1"} and Eqn. [eqn:coupling_coeff]{reference-type="eqref" reference="eqn:coupling_coeff"}, respectively. The total energy will be denoted by $\mathcal{H}(x,y,x_j,y_j,p_x,p_y,p_{x_j},p_{y_j}) = E = E_{\rm saddle} + \Delta E$ where $j=1, 2, \ldots, N_B$ and $\Delta E$ is the excess energy with respect to the isomerization barrier energy. In this article, we are adopting the contraction $x_j$ to denote $x_1, x_2, \ldots, x_{N_B}$; $p_{x_j}$ to denote $p_{x_1}, p_{x_2}, \ldots, p_{x_{N_B}}$, and so forth.

In this form of coupling each bath mode is coupled with all the degrees of freedom of the system, the total number of degrees of freedom are $N_S N_B + N_S$ and as a result the dimension of phase space is $2N_S(N_B + 1)$. Thus, it gives rise to a high dimensional Hamiltonian model that allows us to incorporate a high number of degrees of freedom by increasing the number of bath modes. In this study, we will first focus on varying the number of bath modes, in particular when $N_B = 1, 2, 4, 8, 16$ which corresponds to $8, 12, 20, 36, 68$ dimensional phase space. Even though, the number of bath modes are not large enough to capture the effect of a solvent dynamics in the Langevin sense, but nonetheless this formulation can be extended as long as enough computational time is reserved to obtain trajectories and the associated diagnostic method. However, the high dimensionality of the system-bath model will be a setting to assess how the method of Lagrangian descriptors performs in more realistic chemical systems. This also gives a natural way to model the unimolecular conformational isomerization in a solvent and allows us to vary the friction of the environment from low-density gases to high-density liquids [53].

Model verification

\centering \subfigure[]{\includegraphics[width=0.24\textwidth]{./figures/fig3B2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E1-500_x-px_tau100_NB1.png}} \subfigure[]{\includegraphics[width=0.24\textwidth]{./figures/fig3B2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB1.png}} \centering \subfigure[]{\includegraphics[width=0.24\textwidth]{./figures/fig3C2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E1-500_x-px_tau100_NB1.png}} \subfigure[]{\includegraphics[width=0.24\textwidth]{./figures/fig3C2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB1.png}}

Revealing Phase Space Structures

The equilibrium points of the system-bath model [eqn:Hameq_sys_bath]{reference-type="eqref" reference="eqn:Hameq_sys_bath"} are located at $\overline{q}_s = (0,0,0_j,0_j,0,0,0_j,0_j)$ and at

$$\overline{q}_c = (x_{eq},\pm y_w,(c_{x,j}/\omega_j^2)x_{eq},(c_{y,j}/\omega_j^2)y_w,0,0,0_j,0_j)$$

where $x_{eq}$ is the $x$-coordinate of the equilibrium point in the system model and needs to be obtained using numerical nonlinear root solver. This shows the system coordinates $(x,y)$ of the equilibrium points do not depend on the bath parameters and thus the isomer wells are still at the same location. This lets us use the same surface of section to compute Lagrangian descriptor even when increasing the bath modes.

As noted in the article, the location of the equilibrium points in the system-bath model are dependent on the frequencies and coupling strength of the bath modes. Hence, we would like to analyze the stability of these equilibrium points. The Jacobian of the system-bath Hamiltonian vector field [eqn:Hameq_sys_bath]{reference-type="eqref" reference="eqn:Hameq_sys_bath"} is

$$\begin{aligned} \mathbb{J} = \begin{pmatrix} \text{\normalfont\Large\bfseries 0}& \hspace*{-\arraycolsep}\vline\hspace*{-\arraycolsep}& \mathbf{M} \\ \hline \frac{\partial^2 V_{\rm SB}}{\partial x_i x_j} & \hspace*{-\arraycolsep}\vline\hspace*{-\arraycolsep}& \text{\normalfont\Large\bfseries 0} \end{pmatrix}, \quad \text{where} \; i, j = 1, 2, \ldots , N_S(N_B + 1) \label{eqn:jacobian_sbham_vecfield}\end{aligned}$$

where each block matrix is of size $N_S(N_B + 1) \times N_S(N_B + 1)$. The matrix $\mathbf{M}$ is given by:

$$\begin{aligned} \mathbf{M} = \begin{pmatrix} \frac{1}{m_s} & 0 & 0 & 0 & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & \frac{1}{m_s} & 0 & 0 & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 0 & \frac{1}{m_1} & 0 & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 0 & 0 & \frac{1}{m_2} & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 0 & 0 & 0 & \ddots & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 0 & 0 & 0 & 0 & \cdots & \frac{1}{m_1} & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & 0 & \cdots & 0 & \frac{1}{m_2} & \cdots & 0 \\ 0 & 0 & 0 & 0 & 0 & \cdots & \cdots & \cdots & \ddots & 0 \end{pmatrix}\end{aligned}$$

which is a diagonal matrix with mass of system and bath degrees of freedom along the main diagonal.

The Hessian of the potential energy function for the system-bath model is given by:

$$\begin{aligned} \frac{\partial^2 V_{\rm SB}}{\partial x_i x_j} = \begin{pmatrix} - \frac{\partial^2 V_{\rm DB}}{\partial^2 x^2} - \sum\limits_{j} (\frac{c_{x,j}}{\omega_j})^2 & -\frac{\partial^2 V_{\rm DB}}{\partial y \partial x} & c_{x,1} & c_{x,2} & \cdots & c_{x,N_B} & 0 & 0 & \cdots & 0 \\ -\frac{\partial^2 V_{\rm DB}}{\partial y \partial x} & - \frac{\partial^2 V_{\rm DB}}{\partial^2 y^2} - \sum\limits_{j} (\frac{c_{y,j}}{\omega_j})^2 & 0 & 0 & \cdots & 0 & c_{y,1} & c_{y,2} & \cdots & c_{y,N_B} \\ c_{x,1} & 0 & -\omega_1^2 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ c_{x,2} & 0 & 0 & -\omega_2^2 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ \vdots & 0 & 0 & 0 & \ddots & 0 & 0 & 0 & \cdots & 0 \\ c_{x,N_B} & 0 & 0 & 0 & \cdots & -\omega_{N_B}^2 & 0 & 0 & \cdots & 0 \\ 0 & c_{y,1} & 0 & 0 & \cdots & 0 & -\omega_1^2 & 0 & \cdots & 0 \\ 0 & c_{y,2} & 0 & 0 & \cdots & 0 & 0 & -\omega_2^2 & \cdots & 0 \\ 0 & \vdots & 0 & 0 & \cdots & 0 & 0 & 0 & \ddots & 0 \\ 0 & c_{y,N_B} & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & -\omega_{N_B}^2 \\ \end{pmatrix}\end{aligned}$$

The Jacobian evaluated at the equilibrium points $\overline{x}$ is given by

$$\begin{aligned} \mathbb{J}(\overline{q}) = \begin{pmatrix} \text{\normalfont\Large\bfseries 0}& \hspace*{-\arraycolsep}\vline\hspace*{-\arraycolsep}& \mathbf{M} \\ \hline \frac{\partial^2 V_{\rm SB}}{\partial x_i x_j} (\overline{q}) & \hspace*{-\arraycolsep}\vline\hspace*{-\arraycolsep}& \text{\normalfont\Large\bfseries 0} \end{pmatrix}, \quad \text{where} \; i, j = 1, 2, \ldots , N_S(N_B + 1) \label{eqn:jacobian_sbham_eqpt}\end{aligned}$$

The eigenvalues are given by the solutions of the characteristic equation:

$$\rm{det} \left( \mathbb{J}(\overline{q}) - \lambda \mathbb{I} \right) = 0$$

where $\mathbb{J}(\overline{q})$ and $\mathbb{I}$ are $2N_S(N_B + 1) \times 2N_S(N_B + 1)$ matrices.

\centering Eigenvalues of the equilibrium point $\overline{q}_s$ of the system-bath model [\[eqn:Hameq\_sys\_bath\]](#eqn:Hameq_sys_bath){reference-type="eqref" reference="eqn:Hameq_sys_bath"} where the real part is shown as red dots and imaginary part is shown as blue cross. The pair of purely positive and negative eigenvalues verify the equilibrium point is index-1.{width="0.95\linewidth"}

The linear stability analysis of these equilibrium points gives their stability does not change from the system model by adding the bath modes in the form of Eqn. [eqn:sb_ham]{reference-type="eqref" reference="eqn:sb_ham"}. Next, the total energy of the index-1 equilibrium point at $\overline{q}_s$ is $\mathcal{H}(0,0,0_j,0_j,0,0,0_j,0_j) = \epsilon_s$ which is the same as in the system Hamiltonian [eqn:ham_db]{reference-type="eqref" reference="eqn:ham_db"}.

Implications for Reaction Dynamics

Detecting changes in reactive islands using Lagrangian descriptors

It has been demonstrated that the reactive islands and their hierarchical structure can be used to search for new reactive trajectories using the so-called shooting method [51]. In high dimensional molecular phase space, this can be a non-trivial task due to the available phase space volume and so has been phrased as "harvesting rare trajectories" and computed using transition path sampling methods [54]. In this section, we would like to investigate if the LD contour maps can provide a systematic approach to finding reactive trajectories and hence building the high dimensional analogs of reactive islands in a system-bath model.

$N_S(N_B + 1)$-DoF: $N_B$ bath modes coupled to each of the $N_S$ system DoF

To sample the reactive trajectories, we introduce the following isoenergetic two-dimensional slice in the $2N_S(N_B + 1)$-dimensional phase space of Eqn. [eqn:Hameq_sys_bath]{reference-type="eqref" reference="eqn:Hameq_sysbath"}. $$\begin{aligned} U{xp_x}^- &= \left{(x,y,x_i,y_i,p_x,py,p{xi},p{y_i}) \; | \; y = y_w, x_i = 0, yi = 0, p{xi} = 0, p{y_i} = 0, p_y(\cdot;E) < 0 \right} \label{eqn:sos_Uxpx_systembath} \end{aligned}$$

Influence of the number of bath modes

Let us consider the bath parameters $\eta_{x} = \eta_{y} = 0.01$ and $\omega_{c} = 1.0$ and study the effect of increasing the number of bath modes on the reactive islands of the 2 DoF isomerization.

Changes in reactive islands with increasing number of bath modes for $\tau = 100$ which we considered to be enough integration time for the isomerization given by the system [eqn:vec_field_db]{reference-type="eqref" reference="eqn:vec_field_db"}.

\centering \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3B2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E1-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3B2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3B2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E7-500_x-px_tau100_NB64.png}} \centering \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3C2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E1-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3C2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3C2_eta1e-1_omega2236e-3_tau100/slice1d_backward_lag_desc_systembath_E7-500_x-px_tau100_NB64.png}}

Influence of friction in the bath modes

\centering \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3B2_eta1e-2_omega1e0/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3B2_eta1e-1_omega1e0/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3B2_eta1e0_omega1e0/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}} \centering \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3C2_eta1e-2_omega1e0_tau100/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3C2_eta1e-1_omega1e0_tau100/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}} \subfigure[]{\includegraphics[width=0.32\textwidth]{./figures/fig3C2_eta1e0_omega1e0_tau100/slice1d_backward_lag_desc_systembath_E4-500_x-px_tau100_NB64.png}}

References

  1. G. M. Wieder and R. A. Marcus, “Dissociation and Isomerization of Vibrationally Excited Species. II. Unimolecular Reaction Rate Theory and Its Application,” The Journal of Chemical Physics, vol. 37, no. 8, pp. 1835–1852, Oct. 1962.
  2. J. W. McIver and A. Komornicki, “Structure of transition states in organic reactions. General theory and an application to the cyclobutene-butadiene isomerization using a semiempirical molecular orbital method,” Journal of the American Chemical Society, vol. 94, no. 8, pp. 2625–2633, Apr. 1972.
  3. C. Dugave and L. Demange, “Cis−Trans Isomerization of Organic Molecules and Biomolecules: Implications and Applications ^\textrm†,” Chemical Reviews, vol. 103, no. 7, pp. 2475–2532, Jul. 2003.
  4. T. J. Donohoe, T. J. C. O’Riordan, and C. P. Rosa, “Ruthenium-Catalyzed Isomerization of Terminal Olefins: Applications to Synthesis,” Angewandte Chemie International Edition, vol. 48, no. 6, pp. 1014–1017, Jan. 2009.
  5. C. C. Price and W. H. Snyder, “SOLVENT EFFECTS IN THE BASE-CATALYZED ISOMERIZATION OF ALLYL TO PROPENYL ETHERS,” Journal of the American Chemical Society, vol. 83, no. 7, pp. 1773–1773, Apr. 1961.
  6. T. Halicioǧlu and O. Sinanoǧlu, “Solvent Effects on Cis-Trans Azobenzene Isomerization: A Detailed Application of a Theory of Solvent Effects on Molecular Association,” Annals of the New York Academy of Sciences, vol. 158, no. 1, pp. 308–317, 1969.
  7. S. R. Flom, V. Nagarajan, and P. F. Barbara, “Dynamic solvent effects on large-amplitude isomerization rates. 1. 2-Vinylanthracene,” The Journal of Physical Chemistry, vol. 90, no. 10, pp. 2085–2092, 1986.
  8. E. S. Eberhardt, S. N. Loh, A. P. Hinck, and R. T. Raines, “Solvent effects on the energetics of prolyl peptide bond isomerization,” Journal of the American Chemical Society, vol. 114, no. 13, pp. 5437–5439, 1992.
  9. E. M. Duffy, D. L. Severance, and W. L. Jorgensen, “Solvent effects on the barrier to isomerization for a tertiary amide from ab initio and Monte Carlo calculations,” Journal of the American Chemical Society, vol. 114, no. 19, pp. 7535–7542, 1992.
  10. P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, “Transition Path Sampling: Throwing Ropes Over Rough Mountain Passes, in the Dark,” Annual Review of Physical Chemistry, vol. 53, no. 1, pp. 291–318, 2002.
  11. M. J. Davis, “Bottlenecks to intramolecular energy transfer and the calculation of relaxation rates,” The Journal of Chemical Physics, vol. 83, no. 3, pp. 1016–1031, Aug. 1985.
  12. M. J. Davis, “Phase space dynamics of bimolecular reactions and the breakdown of transition state theory,” The Journal of Chemical Physics, vol. 86, no. 7, pp. 3978–4003, Apr. 1987.
  13. C. C. Marston and N. De Leon, “Reactive islands as essential mediators of unimolecular conformational isomerization: A dynamical study of 3‐phospholene,” The Journal of Chemical Physics, vol. 91, no. 6, pp. 3392–3404, Sep. 1989.
  14. N. De Leon and C. C. Marston, “Order in chaos and the dynamics and kinetics of unimolecular conformational isomerization,” The Journal of Chemical Physics, vol. 91, no. 6, pp. 3405–3425, Sep. 1989.
  15. H. A. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica, vol. 7, no. 4, pp. 284–304, Apr. 1940.
  16. R. Zwanzig, “Nonlinear generalized Langevin equations,” Journal of Statistical Physics, vol. 9, no. 3, pp. 215–220, Nov. 1973.
  17. E. Pollak, H. Grabert, and P. Hänggi, “Theory of activated rate processes for arbitrary frequency dependent friction: Solution of the turnover problem,” The Journal of chemical physics, vol. 91, no. 7, pp. 4073–4087, 1989.
  18. E. Pollak and P. Talkner, “Transition-state recrossing dynamics in activated rate processes,” Physical Review E, vol. 51, no. 3, pp. 1868–1878, Mar. 1995.
  19. N. De Leon and B. J. Berne, “Intramolecular rate process: Isomerization dynamics and the transition to chaos,” The Journal of Chemical Physics, vol. 75, no. 7, pp. 3495–3510, Oct. 1981.
  20. N. De Leon, M. A. Mehta, and R. Q. Topper, “Cylindrical manifolds in phase space as mediators of chemical reaction dynamics and kinetics. I. Theory,” The Journal of Chemical Physics, vol. 94, no. 12, pp. 8310–8328, Jun. 1991.
  21. N. De Leon, M. A. Mehta, and R. Q. Topper, “Cylindrical manifolds in phase space as mediators of chemical reaction dynamics and kinetics. II. Numerical considerations and applications to models with two degrees of freedom,” The Journal of Chemical Physics, vol. 94, no. 12, pp. 8329–8341, Jun. 1991.
  22. N. De Leon, “Cylindrical manifolds and reactive island kinetic theory in the time domain,” The Journal of Chemical Physics, vol. 96, no. 1, pp. 285–297, Jan. 1992.
  23. T. Komatsuzaki and M. Nagoaka, “A dividing surface free from a barrier recrossing motion in many-body systems,” Chem. Phys. Lett., vol. 265, pp. 91–98, 1997.
  24. T. Komatsuzaki and R. S. Berry, “Local regularity and non-recrossing path in transition state: a new strategy in chemical reaction theories,” J. Mol. Struct. THEOCHEM, vol. 506, pp. 55–70, 2000.
  25. S. Wiggins, L. Wiesenfeld, C. Jaffé, and T. Uzer, “Impenetrable Barriers in Phase-Space,” Physical Review Letters, vol. 86, no. 24, pp. 5478–5481, Jun. 2001.
  26. T. Uzer, C. Jaffé, J. Palacián, P. Yanguas, and S. Wiggins, “The geometry of reaction dynamics,” Nonlinearity, vol. 15, no. 4, pp. 957–992, Jul. 2002.
  27. H. Waalkens and S. Wiggins, “Direct construction of a dividing surface of minimal flux for multi-degree-of-freedom systems that cannot be recrossed,” J. Phys. A: Math. Gen., vol. 37, no. 35, p. L435, 2004.
  28. S. Wiggins, “On the geometry of transport in phase space I. Transport in k degree-of-freedom Hamiltonian systems, 2 \le k < ∞,” Physica D, vol. 44, pp. 471–501, 1990.
  29. S. Wiggins, Normally hyperbolic invariant manifolds in dynamical systems, vol. 105. Springer Science & Business Media, 2013.
  30. P. Pechukas and F. J. McLafferty, “On transition state theory and the classical mechanics of collinear collisions,” J. Chem. Phys., vol. 58, no. 4, pp. 1622–1625, 1973.
  31. P. Pechukas and E. Pollak, “Trapped trajectories at the boundary of reactivity bands in molecular collisions,” J. Chem. Phys., vol. 67, no. 12, pp. 5976–5977, 1977.
  32. E. Pollak and P. Pechukas, “Transition states, trapped trajectories, and bound states embedded in the continuum,” J. Chem. Phys., vol. 69, pp. 1218–1226, 1978.
  33. P. Pechukas and E. Pollak, “Classical transition state theory is exact if the transition state is unique,” J. Chem. Phys., vol. 71, no. 5, pp. 2062–2068, 1979.
  34. R. E. Gillilan and G. S. Ezra, “Transport and turnstiles in multidimensional Hamiltonian mappings for unimolecular fragmentation: Application to van der Waals predissociation,” J. Chem. Phys., vol. 94, pp. 2648–2668, 1991.
  35. T. Komatsuzaki and M. Nagaoka, “Study on ‘regularity’ of barrier recrossing motion,” J. Chem. Phys., vol. 105, pp. 10838–10848, 1996.
  36. T. Komatsuzaki and R. S. Berry, “A dynamical propensity rule for transitions in chemical reactions,” J. Phys. Chem. A, vol. 106, pp. 10945–10950, 2002.
  37. S. Wiggins, “The role of normally hyperbolic invariant manifolds (NHIMS) in the context of the phase space setting for chemical reaction dynamics,” Regular and Chaotic Dynamics, vol. 21, no. 6, pp. 621–638, 2016.
  38. C. B. Li, M. Toda, and T. Komatsuzaki, “Bifurcation of no-return transition states in many-body chemical reactions,” J. Chem. Phys., vol. 130, p. 124116, 2009.
  39. M. Inarrea, J. F. Palacian, A. I. Pascual, and J. P. Salas, “Bifurcations of dividing surfaces in chemical reactions,” J. Chem. Phys., vol. 135, p. 014110, 2011.
  40. A. Allahem and T. Bartsch, “Chaotic dynamics in multidimensional transition states,” J. Chem. Phys., vol. 137, p. 214310, 2012.
  41. F. A. L. Mauguière, P. Collins, G. S. Ezra, and S. Wiggins, “Bifurcations of normally hyperbolic invariant manifolds in analytically tractable models and consequences for reaction dynamics,” Int. J. Bifurcat. Chaos, vol. 23, no. 12, 2013.
  42. R. S. MacKay and D. C. Strub, “Bifurcations of transition states: Morse bifurcations,” Nonlinearity, vol. 27, no. 5, p. 859, 2014.
  43. R. S. MacKay and D. C. Strub, “Morse bifurcations of transition states in bimolecular reactions,” Nonlinearity, vol. 28, no. 12, pp. 4303–4329, Jan. 2015.
  44. A. M. Berezhkovskii, E. Pollak, and V. Y. Zitserman, “Activated rate processes: Generalization of the Kramers-Grote-Hynes and Langer theories,” J. Chem. Phys., vol. 97, no. 4, pp. 2422–2437, 1992.
  45. E. Hershkovitz and E. Pollak, “Multidimensional generalization of the Pollak-Grabert-Hanggi turnover theory for activated rate processes,” J. Chem. Phys., vol. 106, no. 18, pp. 7678–7699, 1997.
  46. S. K. Reese and S. C. Tucker, “Curvilinear-path based theory of the energy transfer limited rate of a two-dimensional solute in a dissipative bath,” Chemical Physics, vol. 235, pp. 171–187, 1998.
  47. I. E. Craig and D. E. Manolopoulos, “Chemical reaction rates from ring polymer molecular dynamics,” J. Chem. Phys., vol. 122, p. 084106, 2005.
  48. A. M. Mancho, S. Wiggins, J. Curbelo, and C. Mendoza, “Lagrangian Descriptors: A Method for Revealing Phase Space Structures of General Time Dependent Dynamical Systems,” Communications in Nonlinear Science and Numerical Simulation, vol. 18, no. 12, pp. 3530–3557, Dec. 2013.
  49. A. Junginger and R. Hernandez, “Uncovering the Geometry of Barrierless Reactions Using Lagrangian Descriptors,” The Journal of Physical Chemistry B, vol. 120, no. 8, pp. 1720–1725, Mar. 2016.
  50. A. Junginger et al., “Transition state geometry of driven chemical reactions on time-dependent double-well potentials,” Physical Chemistry Chemical Physics, vol. 18, no. 44, pp. 30270–30281, 2016.
  51. S. Patra and S. Keshavamurthy, “Detecting reactive islands using Lagrangian descriptors and the relevance to transition path sampling,” Physical Chemistry Chemical Physics, vol. 20, no. 7, pp. 4970–4981, 2018.
  52. S. Naik, V. J. García-Garrido, and S. Wiggins, “Finding NHIM: Identifying high dimensional phase space structures in reaction dynamics using Lagrangian descriptors,” Communications in Nonlinear Science and Numerical Simulation, vol. 79, p. 104907, 2019.
  53. S. C. Tucker, “Reaction rates in condensed phases. Perspective on ‘Brownian motion in a field of force and the diffusion model of chemical reactions,’” Theoretical Chemistry Accounts: Theory, Computation, and Modeling (Theoretica Chimica Acta), vol. 103, no. 3-4, pp. 209–211, Feb. 2000.
  54. S. Patra and S. Keshavamurthy, “Classical-quantum correspondence in a model for conformational dynamics: Connecting phase space reactive islands with rare events sampling,” Chemical Physics Letters, vol. 634, pp. 1–10, Aug. 2015.