# $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.

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)$

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}$$`

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 {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_sys*bath"}. $$\begin{aligned}
U*{xp_x}^- &= \left{(x,y,x_i,y_i,p_x,p*y,p*{x*i},p*{y_i}) \; | \; y = y_w, x_i = 0, y*i = 0, p*{x*i} = 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

- 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - R. Zwanzig, “Nonlinear generalized Langevin equations,”
*Journal of Statistical Physics*, vol. 9, no. 3, pp. 215–220, Nov. 1973. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - 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. - S. Wiggins,
*Normally hyperbolic invariant manifolds in dynamical systems*, vol. 105. Springer Science & Business Media, 2013. - 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. - 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. - 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. - 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. - 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. - T. Komatsuzaki and M. Nagaoka, “Study on ‘regularity’ of barrier recrossing motion,”
*J. Chem. Phys.*, vol. 105, pp. 10838–10848, 1996. - 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. - 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. - 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. - 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. - A. Allahem and T. Bartsch, “Chaotic dynamics in multidimensional transition states,”
*J. Chem. Phys.*, vol. 137, p. 214310, 2012. - 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. - R. S. MacKay and D. C. Strub, “Bifurcations of transition states: Morse bifurcations,”
*Nonlinearity*, vol. 27, no. 5, p. 859, 2014. - 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. - 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. - 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. - 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. - I. E. Craig and D. E. Manolopoulos, “Chemical reaction rates from ring polymer molecular dynamics,”
*J. Chem. Phys.*, vol. 122, p. 084106, 2005. - 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. - 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. - 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. - 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. - 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. - 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. - 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.