# Probing photoinduced proton coupled electron transfer process by means of
two-dimensional resonant electronic-vibrational spectroscopy

Jiaji Zhang<EMAIL_ADDRESS>Department of Chemistry,
Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan  Raffaele
Borrelli<EMAIL_ADDRESS>DISAFA, University of Torino, Largo Paolo
Braccini 2, I-10095 Grugliasco, Italy  Yoshitaka Tanimura
<EMAIL_ADDRESS>Department of Chemistry, Graduate School of
Science, Kyoto University, Kyoto 606-8502, Japan

###### Abstract

We develop a detailed theoretical model of photo-induced proton-coupled
electron transfer (PPCET) processes, which are at the basis of solar energy
harvesting in biological systems and photovoltaic materials. Our model enables
to analyze the dynamics and the efficiency of a PPCET reaction under the
influence of a thermal environment by disentangling the contribution of the
fundamental electron transfer (ET) and proton transfer (PT) steps. In order to
study quantum dynamics of the PPCET process under an interaction with non-
Markovian environment we employ the hierarchical equations of motion (HEOM).
We calculate transient absorption spectroscopy (TAS) and a newly defined two-
dimensional resonant electronic-vibrational spectroscopy (2DREVS) signals in
order to study the nonequilibrium reaction dynamics. Our results show that
different transition pathways can be separated by TAS and 2DREVS.

Electron coupled proton transfer, Two-dimensional electronic-vibrational
spectroscopies, Hierarchical Equations of motion

##  I Introduction

The simultaneous transfer of protons and electrons plays an important role in
many natural and artificial energy conversion processes. A typical example is
the oxygen evolving complex (OEC) of natural photosynthetic system, where the
oxygen generation consists of four stepwise proton-coupled electron transfer
(PCET) catalyzed reactions. Migliore  _et al._ (2014); Megiatto Jr  _et al._
(2014); Huynh  _et al._ (2017); Offenbacher and Barry (2020) specific pathways
taken by electrons and protons, can lead to step-wise (consecutive) or
concerted type reactions (CEPT). Unravelling the detailed mechanistic aspects
of the PCET process is fundamental for the design of artificial solar energy
utilization systems, for example, dye sensitized photo-electrochemical cell
(DS-PEC) and many other bio-mimetic systems which have been developed for
solar energy utilization and hydrogen reduction.Gagliardi  _et al._ (2010);
Becker  _et al._ (2016); Yamaguchi  _et al._ (2014); Allen (2002); Papa,
Capitanio, and Luca Martino (2006).

Various approximated quantum dynamical theories, mostly based on the
determination of reaction rate constants, have been derived for isolated
systems on the basis of the golden rule expression, linear response theory,
and on Marcus’s theory of electron-transfer (ET) processes. Cukier (1996);
Goldsmith, Soudackov, and Hammes-Schiffer (2019) Their applications have been
extended to condensed phase systems by further assuming a perturbative system-
bath interaction and a classical treatment of an environment representing, for
example, solvent.Mayer (2011); Soudackov and Hammes-Schiffer (2016, 2000,
2015); Cukier (1994, 1996); Hazra, Soudackov, and Hammes-Schiffer (2010) Rate
constants for several PCET systems in thermal equilibrium conditions have also
been computed with the aids of molecular dynamics simulations and quantum
chemistry calculations.Kretchmer and Miller (2016); Kobrak and Hammes-Schiffer
(2001); Grimminger and Schmickler (2007); Auer, Fernandez, and Hammes-Schiffer
(2011); Peluso, Di Donato, and Saracino (2000)

Yet, the sole computation of reaction rates does not provide enough
information to fully disentangle different ET and PT pathways and can hide
important information about the role of the environment. Ultrafast nonlinear
spectroscopy can be a powerful tool for unravelling the mechanistic aspects of
PPCET reactions and of photosynthesis in general. Cho (2019) For example,
infrared (IR) transient absorption spectroscopy (TAS) has been applied to
excited-state proton transfer and chemical bond cleavage, and can provide a
versatile tool to determine the relaxation mechanism after an initial
photoexcitation. Pfeiffer  _et al._ (1999); Hogle, Cunningham, and Tucker
(2018) The results of luminescence TAS have indicated that the quantum effect
of donor-acceptor (D-A) vibrations on PPCET is important for a full quantum
treatment of the total reaction system. Kuss-Petermann  _et al._ (2012); Giret
_et al._ (2020); Zheng  _et al._ (2011); Song and Shi (2017); Barbara, Walker,
and Smith (1992)

These spectroscopic techniques have also been extended to multi-dimensional
cases. Two-dimensional (2D) vibrational spectroscopy (2DVS)Tanimura and
Mukamel (1993); Mukamel (1995); Hamm and Zanni (2011) and 2D electronic
spectroscopy (2DES) have been applied to condense phase transition and
succeeded in investigating the electronic excitation dynamics and a structural
change of molecules.Gaynor  _et al._ (2016); Oliver (2018); Fox, Blair, and
Khalil (2020); Mukamel (2000); Dong  _et al._ (2015) Their combination, 2D
electronic-vibrational spectroscopy (2DEVS), has also been successfully
applied to the photo-isomerization reactions, mental-to-ligands transitions,
conical intersection wavepacket dynamics, and ultrafast excitonic
photosynthetic energy transfer reactions. Wu  _et al._ (2019a); Ikeda and
Tanimura (2018); Lewis and Fleming (2016); Courtney  _et al._ (2015); Wu  _et
al._ (2019b) By utilizing the UV-vis and IR pulses, we are now able to measure
the coupling strength and coherence between the electronic and vibrational
transitions as the off-diagonal peaks of 2D spectroscopy. Cho (2019); Mukamel
(1995) These features are useful for the investigation of PPCET reaction
dynamics.

In this paper, we present a model of a PPCET reaction and provide a detailed
analysis of its dynamics by computing the signals of TAS and a newly defined
2D resonant electronic-vibrational spectroscopy (2DREVS). The 2DREVS is an
extension of 2DEV for a strong resonant reaction system, and is useful for
investigating the dynamics of PPCET reaction, as described below. We describe
the coupled proton-electron dynamics using two-dimensional potential energy
surfaces (PESs), and complex system-bath interactions to simulate a system in
realistic conditions. We employ the numerically “exact” hierarchical equations
of motion (HEOM) approach to study the reduced system dynamics under non-
perturbative and non-Markovian system-bath interactions at finite temperature.
Tanimura and Kubo (1989); Tanimura (1990); Ishizaki and Tanimura (2005);
Tanimura (2006, 2014, 2015, 2020) The paper is organized as follows. In Sec.
II, we derive a system-bath model for a prototypical PPCET process and
introduce the HEOM approach for numerical simulation. The theory of nonlinear
response functions is also briefly sketched in this section. In Sec. III, we
present the calculated TAS and 2DREVS results and analyze their profiles.

##  II Theory

###  II.1 Model Hamiltonian

The system considered in the present work is depicted in Fig. 1. In the ground
electronic state the proton is localized at bond distance from donor $D$, and
the $D$-$H$ moiety is hydrogen bonded to the acceptor $A$.

Figure 1:  Model PPCET system with a hydrogen bridge. Here, $D(A)$ is the
donor (acceptor), $H$ is the transferring proton, $x$ is the proton coordinate
describing its distance from the center of the $D$ and $A$ units, and $Q$ is
the distance between the heavy atoms connected by the hydrogen bond.

The $x$ coordinate describes the position of the proton between $D$ and $A$,
while $Q$ is the distance between the heavy atoms which is also referred to as
the reaction promoting mode. We wish to describe the dynamics of the system
resulting from the photo-excitation of $D$, which is followed by a coupled
transfer of an electron and a proton to the $A$ moiety. As a result of the
process an hydrogen atom is transferred from $D$-$H$ to $A$, i.e. $A$ is
reduced to $A$-$H$ and $D$-$H$ is oxidized $D$.

In order to model the coupled PT and ET processes we consider an electronic
active space comprising the ground state $|\phi_{g}\rangle$ of the system, the
localized excited state $|\phi_{LE}\rangle$, in which only the moiety $D$ is
in the first excited electronic state, while $A$ is in the ground electronic
state, and the charge-transfer state $|\phi_{CT}\rangle$, in which $D$ has
transferred an electron to $A$. The diabatic representation of system is shown
in Fig. 2. The motion along the $x$ and $Q$ coordinates is described employing
realistic two-dimensional potential energy surfaces. Furthermore, we assume
that the system interacts with a condensed phase environment which can be
either a solvent or a protein scaffold. The overall Hamiltonian can therefore
be expressed as

$\displaystyle\hat{H}=\sum_{i}\hat{H}_{i}(\hat{x},\hat{Q},\\{q_{a}\\})\big{|}i\rangle\langle
i\big{|}+\sum_{i\neq j}\Delta_{ij}\big{|}i\rangle\langle
j\big{|}+\hat{H}_{B},$ (1)

where $\hat{H}_{i}(\hat{x},\hat{Q},\\{q_{a}\\})$ is the Hamiltonian for the
electronic states $i=g,CT$, and $LE$, $\Delta_{ij}$ are the electronic
couplings among different electronic states, and $\hat{H}_{B}$ is the
Hamiltonian of the thermal bath which is modeled as a collection of harmonic
oscillators

$\displaystyle\hat{H}_{B}=\sum_{a}\left(\frac{\hat{p}_{a}^{2}}{2m_{a}}+\frac{1}{2}m_{a}\omega_{a}^{2}\hat{q}_{a}^{2}\right),$
(2)

where $\hat{p}_{a}$, $\hat{q}_{a}$, $m_{a}$ and $\omega_{a}$ are the momentum,
position, mass and frequency of ${a}^{th}$ bath oscillator, respectively.

Figure 2:  Diabatic representation of the PES for the reduced system. The
black curve represent the ground state ${\left|\phi_{g}\right\rangle}$. The
red, and blue curves represent the local excited state
${\left|\phi_{LE}\right\rangle}$, and charge transfer state
${\left|\phi_{CT}\right\rangle}$, respectively.

The operators $\hat{H}_{i}$ can be explicitly written in the form

$\displaystyle\hat{H}_{i}(\hat{x},\hat{Q},\\{q_{a}\\})=\frac{\hat{p}_{x}^{2}}{2m_{x}}+\frac{\hat{p}_{Q}^{2}}{2m_{Q}}+U_{i}(\hat{x},\hat{Q},\\{q_{a}\\})+\varepsilon_{i}.$
(3)

Here, $\hat{x}$, $\hat{p}_{x}$, $m_{x}$, and $\hat{Q}$, $\hat{p}_{Q}$ $m_{Q}$
are the coordinate, momentum and mass of the proton and of the D-A vibration,
respectively, and $\varepsilon_{i}$ is the the energy of the electronic state
for $i$. The potential $U_{i}(\hat{x},\hat{Q},\\{q_{a}\\})$ describes the
variation of the electronic energy as a function of the coordinates $x,Q$ and,
furthermore, explicitly includes the interaction between these coordinates and
the bath modes $\\{q_{a}\\}$. Following previous work, Zhang, Borrelli, and
Tanimura (2020) we use an asymmetric double well Morse potential for the
proton mode, and a harmonic potential for the D-A mode

$\displaystyle U_{i}(\hat{x},\hat{Q},\\{q_{a}\\})$
$\displaystyle=D_{i}^{l}\left[1-e^{-\alpha(\hat{Q}/2+\hat{x}-{x_{e}(\\{q_{a}\\})})}\right]^{2}+D_{i}^{r}\left[1-e^{-\alpha(\hat{Q}/2-\hat{x}-x_{e}(\\{q_{a}\\}))}\right]^{2}$
$\displaystyle+\frac{1}{2}D_{k}{(\hat{Q}-{Q_{e}})}^{2},$ (4)

where $D_{i}^{l}$ and $D_{i}^{r}$ are the dissociation energy of donor (left
well) and acceptor (right well), ${x_{e}}$ and ${Q_{e}}$ are the equilibrium
distance of the proton and D-A vibrations, $\alpha$ represents the curvature
of the Morse potentials, and $D_{k}$ is the force constant of the D-A
vibration. In our model the role of bath modes $\\{q_{a}\\}$ is to dynamically
perturb the equilibrium position of the proton via a linear interaction, that
is

$\displaystyle{x_{e}}(\\{q_{a}\\})={x_{0}}-\sum_{a}g_{a}\hat{q}_{a},$ (5)

where ${x_{0}}$ is the equilibrium distance without the heat bath, and
$\\{g_{a}\\}$ are coupling strength parameters. Finally, we simplify this
potential by expanding Eq. (4) in terms of the collective coordinate
$\hat{X}=\sum_{a}g_{a}q_{a}$ up to the first-order, which is similar to the
reaction surface approach. Tew, Handy, and Carter (2006) The system PES and
the resulting exponential-linear (EL) system-bath interaction are then
expressed as

$\displaystyle{U_{i}}(\hat{x},\hat{Q},X)=U_{i}^{0}(\hat{x},\hat{Q})+\hat{V}_{i}(\hat{x},\hat{Q})\hat{X},$
(6)

where

$\displaystyle U_{i}^{0}(\hat{x},\hat{Q})$
$\displaystyle=D_{i}^{l}\left[1-e^{-\alpha(\hat{Q}/2+\hat{x}-x_{0})}\right]^{2}+D_{i}^{r}\left[1-e^{-\alpha(\hat{Q}/2-\hat{x}-x_{0})}\right]^{2}$
$\displaystyle+\frac{1}{2}D_{K}{(\hat{Q}-Q_{0})}^{2}.$ (7)

The operator $V_{i}$ depends solely on system variables and can be explicitly
written asZhang, Borrelli, and Tanimura (2020)

$\displaystyle\hat{V}_{i}(\hat{x},\hat{Q})$ $\displaystyle=2\alpha
D_{i}^{l}\left[1-e^{-\alpha(\hat{Q}/2+\hat{x}-x_{0})}\right]e^{-\alpha(\hat{Q}/2+\hat{x}-x_{0})}$
(8) $\displaystyle+2\alpha
D_{i}^{r}\left[1-e^{-\alpha(\hat{Q}/2-\hat{x}-x_{0})}\right]e^{-\alpha(\hat{Q}/2-\hat{x}-x_{0})}.$
(9)

The structure of this rather complex form of system operator can be easily
understood once we expand it in terms of $\hat{x}$ and $\hat{Q}$ as

$\hat{V}_{i}(\hat{x},\hat{Q})=V_{i}^{(0)}+V_{i,x}^{(1)}\hat{x}+V_{i,Q}^{(1)}\hat{Q}++V_{i,x}^{(2)}\,\hat{x}^{2}+V_{i,xQ}^{(2)}\,\hat{x}\hat{Q}+...,$
(10)

where the $V_{i}^{(0)}$, $V_{i,x}^{(1)}$, etc. are constants whose analytical
expressions are given in Appendix A. Hence it is clear that the coupling of
Eq. (5) introduces linear interactions with the electronic subsystem via the
constant $V_{i}^{(0)}$, and with the nuclear coordinates $\hat{x}$ and
$\hat{Q}$ via $V_{i,x}^{(1)}$ and $V_{i,Q}^{(1)}$. As discussed in Ref.
Tanimura, 2006, the linear-linear (LL) interaction, such as
$V_{i,x}^{(1)}\hat{x}q_{a}$ contributes mainly to energy relaxation, while the
square-linear (SL) system-bath interaction, such as
$V_{i,x}^{(2)}\,\hat{x}^{2}q_{a}$ leads to vibrational dephasing in the slow
modulation case, due to the frequency fluctuation of the system vibrations.
Finally we note that the LL contribution in the proton mode vanishes for
symmetric double well potential, i.e. $D_{i}^{l}=D_{i}^{r}$. Zhang, Borrelli,
and Tanimura (2020)

For simplicity, we further assume that all of the electronic states are
coupled to the same heat bath. Then, Eq. (1) can be rewritten as

$\displaystyle\hat{H}_{\rm tot}$
$\displaystyle=\sum_{i}\hat{H}_{i}^{0}{\big{|}i\big{\rangle}\big{\langle}i\big{|}}+\sum_{i\neq
j}\Delta_{ij}{\big{|}i\big{\rangle}\big{\langle}j\big{|}}+\hat{H}_{B+I},$ (11)

where

$\displaystyle\hat{H}_{i}^{0}=\frac{\hat{p}_{x}^{2}}{2m_{x}}+\frac{\hat{p}_{Q}^{2}}{2m_{Q}}+U_{i}^{0}(\hat{x},\hat{Q})$
(12)

and

$\displaystyle\hat{H}_{B+I}$
$\displaystyle=\sum_{a}\left\\{\frac{\hat{p}_{a}}{2m_{a}}+\frac{m_{a}\omega_{a}^{2}}{2}\left[\hat{q}_{a}-\frac{g_{a}\left(\sum_{i}{\big{|}i\big{\rangle}\big{\langle}i\big{|}}\hat{V}_{i}\right)}{m_{a}}\omega_{a}^{2}\right]^{2}\right\\}.$
(13)

We also include the counter-term in the definition of $\hat{{H}}_{B+I}$ in
order to maintain the translational symmetry of the system. Tanimura (2006,
2014)

###  II.2 Hierarchical Equations of Motion Approach

Next, we briefly introduce the hierarchical equations of motion (HEOM)
approach, which is employed to investigate quantum dynamics of the PCET system
in a numerically rigorous way.Tanimura (2006, 2020) We can also employ the
multistate quantum hierarchical Fokker-Planck equations (MQHFPE), which has
been applied to both optical and nonadiabatic transition problems described by
complex PESs. Tanimura and Maruyama (1997); Ikeda and Tanimura (2017); Ikeda,
Dijkstra, and Tanimura (2019) However, here we choose the regular HEOM in the
energy eigenstate representation for both electronic and vibrational modes.
This is because the proton motion is well confined in the PESs, and the
computational cost for using MQHFPE is much higher that regular HEOM.

The heat bath is described by the spectral distribution function (SDF)

$\displaystyle
J(\omega)=\pi\sum_{a}\frac{g_{a}^{2}}{2m_{{a}}\omega_{a}}\delta(\omega-\omega_{a})$
(14)

and the inverse temperature, $\beta\hbar=1/k_{B}T$, where $k_{B}$ is the
Boltzmann constant. The overall noise effect on the system is characterized by
the correlation function

$\displaystyle
C(t)=\hbar\int_{0}^{\infty}{\rm{d}}{\omega}J(\omega)\left[{\coth}\left(\frac{\beta\hbar\omega}{2}\right){\cos}(\omega
t)-i{\sin}(\omega t)\right],$ (15)

where the notation $\langle\ldots\rangle_{B}$ represents the thermal average
taken with the canonical distribution of the bath. In this paper, we use a
Drude formed SDF,

$\displaystyle
J(\omega)=\frac{\zeta}{2\pi}\frac{\omega\gamma^{2}}{\gamma^{2}+\omega^{2}},$
(16)

where $\zeta$ represents the coupling strength, and $\gamma$ is the reciprocal
of the noise correlation time, representing the width of the spectral
distribution. Then, Eq. (15) can be expressed in terms of a combination of
linear exponential functions and of the $\delta(t)$ function, as

$\displaystyle
C(t)=\sum_{k=0}^{K}(c_{k}^{\prime}+ic_{k}^{\prime\prime})\gamma_{k}e^{-\gamma_{k}t}+2c_{\delta}\cdot\delta(t),$
(17)

where $c_{k}^{\prime}$, $c_{k}^{\prime\prime}$, $\gamma_{k}$ and $c_{\delta}$
are constants determined by the chosen decomposition method. Here we employ
the Padé decomposition method Hu  _et al._ (2011); Ding, Xu, and Yan (2012)
which is known to enhance the efficiency of numerical calculations. By
introducing the auxiliary density operators (ADO) $\hat{\rho}_{\vec{n}}$, the
HEOM can be derived as Tanimura and Kubo (1989); Tanimura (1990); Ishizaki and
Tanimura (2005); Tanimura (2006, 2014, 2015, 2020)

$\displaystyle\frac{\partial}{\partial t}\hat{\rho}_{\vec{n}}(t)=$
$\displaystyle\,-\left[\frac{i}{\hbar}\mathcal{\hat{L}}_{S}+\sum_{k}n_{k}\gamma_{k}+c_{\delta}\hat{\Phi}^{2}\right]\hat{\rho}_{\vec{n}}(t)$
$\displaystyle\,-\sum_{k}\hat{\Phi}\hat{\rho}_{\vec{n}+\vec{e}_{k}}(t)-\sum_{k}n_{k}\hat{\Theta}_{k}\hat{\rho}_{\vec{n}-\vec{e}_{k}}(t),$
(18)

where the superoperators are defined as
${\mathcal{\hat{L}}}_{S}{\hat{A}}\equiv[{{\hat{H}}_{S}},{{\hat{A}}}]$ and
$\hat{\Theta}_{k}\equiv
c_{k}^{\prime}\hat{\Phi}-c_{k}^{\prime\prime}\hat{\Psi}$ with

$\displaystyle{\hat{\Phi}}{\hat{A}}\equiv\frac{i}{\hbar}\left[\sum_{i}{\big{|}i\big{\rangle}\big{\langle}i\big{|}}\hat{V}_{i},{{\hat{A}}}\right],~{}~{}~{}{\hat{\Psi}}{\hat{A}}\equiv\frac{1}{\hbar}\left\\{\sum_{i}{\big{|}i\big{\rangle}\big{\langle}i\big{|}}\hat{V}_{i},{{\hat{A}}}\right\\},$
(19)

for any physical operator $\hat{A}$. The components of multi-index vector
$\vec{n}=(...,n_{k},...)$ are all non-negative integers, and $\vec{e}_{k}$ is
the $k^{th}$ unit vector. In HEOM formalism, only the first element,
$\vec{n}=(0,...,0)$, has a physical meaning, corresponding to the reduced
density operator of system. The others are served as the treatment of non-
perturbative and non-Markovian heat bath effect. Tanimura (2006, 2020)
Although Eq. (18) consists of infinite equations, we can truncate it at a
properly chosen large $N$ value, for $N=\sum_{k}n_{k}$. Ishizaki and Tanimura
(2005) In order to reduce the computational cost for the time integration, we
rescale the ADOs as
$\hat{\rho}_{\vec{n}}=\hat{\rho}_{\vec{n}}/{\prod_{u,k}\sqrt{n_{u,k}!}}$.
Then, Eqs. (18) are rewritten as Shi  _et al._ (2009); Ikeda and Tanimura
(2019)

$\displaystyle\frac{\partial}{\partial t}\hat{\rho}_{\vec{n}}(t)=$
$\displaystyle\,-\left[\frac{i}{\hbar}\mathcal{\hat{L}}_{S}+\sum_{k}n_{k}\gamma_{k}+c_{\delta}\hat{\Phi}^{2}\right]\hat{\rho}_{\vec{n}}(t)$
$\displaystyle\,-\sum_{k}\sqrt{n_{k}+1}\hat{\Phi}\hat{\rho}_{\vec{n}+\vec{e}_{k}}(t)-\sum_{k}\sqrt{n_{k}}\hat{\Theta}_{k}\hat{\rho}_{\vec{n}-\vec{e}_{k}}(t).$
(20)

###  II.3 Projection operators for PT and ET states

In order to analyze the PCET process, next we introduce a set of projection
operators defined as

$\displaystyle\hat{\mathcal{\theta}}_{i}^{l}=\big{|}i\rangle\langle
i\big{|}\hat{h}(-x),~{}~{}~{}\hat{\mathcal{\theta}}_{i}^{r}=\big{|}i\rangle\langle
i\big{|}\hat{h}(x),$ (21)

where $i=CT$ and $LE$, $\hat{h}(x)$ is the Heaviside step function for the
proton coordinate, and the symbols $l$ and $r$ represent the proton localized
in the left (donor) and right (acceptor) well, respectively. The corresponding
population of the superposition is
$P_{i}^{\alpha}(t)=\rm{Tr}\left\\{\hat{\mathcal{\theta}}_{i}^{\alpha}\hat{\rho}(t)\right\\}$
for $\alpha=l$ or $r$. The populations of ${\left|\phi_{LE}\right\rangle}$ and
${\left|\phi_{CT}\right\rangle}$ is then separated as
$P_{i}(t)=P_{i}^{l}(t)+P_{i}^{r}(t)$, whereas that in the left and right well
is expressed as $P^{\alpha}(t)=P_{LE}^{\alpha}(t)+P_{CT}^{\alpha}(t)$.

As shown in Fig. 2, the superposition ${\left|\phi_{LE}^{l}\right\rangle}$
represents the configuration $D^{*}-H\cdots A$, and
${\left|\phi_{LE}^{r}\right\rangle}$ represents ${\rm{D^{*}}}^{-}\cdots
H-{\rm{A}}^{+}$. Similarly, ${\left|\phi_{CT}^{l}\right\rangle}$ represents
$D^{+}-H\cdots A^{-}$, and ${\left|\phi_{CT}^{r}\right\rangle}$ represents
${\rm{D}}\cdots H-{\rm{A}}$. Thus, the pure PT process corresponds to the
transitions
${\left|\phi_{LE}^{l}\right\rangle}\leftrightarrow{\left|\phi_{LE}^{r}\right\rangle}$
and
${\left|\phi_{CT}^{l}\right\rangle}\leftrightarrow{\left|\phi_{CT}^{r}\right\rangle}$.
The pure ET process corresponds to the transitions
${\left|\phi_{LE}^{l}\right\rangle}\leftrightarrow{\left|\phi_{CT}^{l}\right\rangle}$
and
${\left|\phi_{LE}^{r}\right\rangle}\leftrightarrow{\left|\phi_{CT}^{r}\right\rangle}$.
The CEPT process corresponds to the transitions
${\left|\phi_{LE}^{l}\right\rangle}\leftrightarrow{\left|\phi_{CT}^{r}\right\rangle}$
and
${\left|\phi_{LE}^{r}\right\rangle}\leftrightarrow{\left|\phi_{CT}^{l}\right\rangle}$.

###  II.4 Nonlinear Response Function

The nonlinear response functions can be calculated within the framework of the
HEOM formalism. Tanimura (2006, 2020) The third-order optical response
function can be expressed as

$\displaystyle
R^{(3)}(t_{3},t_{2},t_{1})={\left(\frac{i}{\hbar}\right)}^{3}{\rm{Tr}}\left\\{\hat{\mu}_{4}\mathcal{G}(t_{3})\hat{\mu}_{3}^{\times}\mathcal{G}(t_{2})\hat{\mu}_{2}^{\times}\mathcal{G}(t_{1})\hat{\mu}_{1}^{\times}{\hat{\rho}^{eq}}\right\\},$
(22)

where $\hat{\mu}_{k}$ is the dipole operator of the $k$-th laser interaction,
${\mathcal{G}}(t)$ is the Green’s function of the total Hamiltonian without
laser interactions, and ${\hat{\rho}^{eq}}$ is the initial state density
operator. In the HEOM approach, the density matrix is replaced by a reduced
one, and ${\mathcal{G}}(t)$ is evaluated from Eq. (18) (or Eq.(20)). Tanimura
(2006) The operator $\hat{\mu}_{k}^{\times}$ is the commutator of the dipole
operator $\hat{\mu}_{i}$. The right-hand side of Eq. (22) can be evaluated as
follows: The system is first in the initial equilibrium state
${\hat{\rho}^{eq}}$, and is excited by the first interaction
$\hat{\mu}_{1}^{\times}$ at $t=0$. The time evolution is computed by
numerically integrating Eq. (18) up to a chosen time $t_{1}$. Then, the system
is excited by the second and third interactions $\hat{\mu}_{2}^{\times}$ and
$\hat{\mu}_{3}^{\times}$ in a similar way. The final signal is computed by the
expectation value of $\hat{\mu}_{4}$. We compute $R^{(3)}(t_{3},t_{2},t_{1})$
for a set of values of $t_{1}$, $t_{2}$, and $t_{3}$.

Here we assume that the PES of $|\phi_{g}\rangle$ and $|\phi_{LE}\rangle$ have
the same equilibrium positions, their energy difference is large and the
population relaxation in the excited states is small. The direct excitation
from ${\left|\phi_{g}\right\rangle}$ to ${\left|\phi_{CT}\right\rangle}$ is
also prohibited. Thus, the initial state is described by the thermal
equilibrium distribution of the $|\phi_{g}\rangle$ as
${\hat{\rho}^{eq}}={\hat{\rho}_{g}^{eq}}$. Assuming the PPCET reaction is
initialized by a pair of impulsive pump pulses that excite the system from
${\left|\phi_{g}\right\rangle}$ to ${\left|\phi_{LE}\right\rangle}$, we set
the initial conditions as
${\hat{\rho}}^{(2)}(0)=-\hat{\mu}_{1}^{\times}\hat{\mu}_{2}^{\times}{\hat{\rho}_{g}^{eq}}/\hbar^{2}$
for further response function analysis. With the previous assumption, we can
further set ${\hat{\rho}}^{(2)}(0)={\hat{\rho}_{LE}^{eq}}$, where
${\hat{\rho}_{LE}^{eq}}$ is evaluated as the steady state solution of the HEOM
for the $|\phi_{LE}\rangle$ state without non-adiabatic coupling with the
$|\phi_{CT}\rangle$. Thus, our discussion in the following only considers the
dynamics between ${\left|\phi_{LE}\right\rangle}$ and
${\left|\phi_{CT}\right\rangle}$.

The transient absorption response function can be evaluated from Eq. (22) by
keeping $t_{1}=0$ as

$\displaystyle
R^{\mathrm{TA}}(t,t^{\prime})=\frac{i}{\hbar}{\rm{Tr}}\left\\{\hat{\mu}_{4}\mathcal{G}(t)\hat{\mu}_{3}^{\times}\mathcal{G}(t^{\prime}){\hat{\rho}}^{(2)}(0)\right\\}.$
(23)

Transient absorption spectrum (TAS) at different $t^{\prime}$ is evaluated as

$\displaystyle I^{\mathrm{TA}}(\omega,t^{\prime})$
$\displaystyle\equiv\omega\mathrm{Im}\int_{0}^{\infty}{\rm{d}}te^{i\omega
t}R^{\mathrm{TA}}(t,t^{\prime}),$ (24)

which also corresponds to linear absorption spectrum for non-equilibrium
initial conditions. For the calculations of TAS, the dipole operators
$\hat{\mu}_{k}$ for $k\geq 3$ are assumed to be either electron part
$\hat{\mu}_{e}$, or proton part $\hat{\mu}_{p}$, defined as

$\displaystyle\hat{\mu}_{e}$
$\displaystyle={\left|\phi_{LE}\right\rangle}{\left\langle\phi_{CT}\right|}+{\left|\phi_{CT}\right\rangle}{\left\langle\phi_{LE}\right|},$
$\displaystyle\hat{\mu}_{p}$
$\displaystyle=\hat{x}\cdot\left({\left|\phi_{LE}\right\rangle}{\left\langle\phi_{LE}\right|}+{\left|\phi_{CT}\right\rangle}{\left\langle\phi_{CT}\right|}\right).$
(25)

Here, $\hat{\mu}_{e}$ is for spectroscopy of the electronic subsystem, and
$\hat{\mu}_{p}$ is for spectroscopy of vibrational degrees of freedom,
respectively.

The fifth-order transient 2D spectroscopy is defined in a similar way as

$\displaystyle
R^{\mathrm{(5)}}(t_{3},t_{2},t_{1})={\left(\frac{i}{\hbar}\right)}^{3}{\rm{Tr}}\left\\{\hat{\mu}_{6}\mathcal{G}(t_{3})\hat{\mu}_{5}^{\times}{\mathcal{G}}(t_{2})\hat{\mu}_{4}^{\times}\mathcal{G}(t_{1})\hat{\mu}_{3}^{\times}{\hat{\rho}}^{(2)}(0)\right\\},$
(26)

where ${\hat{\rho}}^{(2)}(0)$ is the same as TAS. The transient 2D correlation
spectroscopy are then evaluated as

$\displaystyle I^{\mathrm{Corr}}(\omega_{3},t_{2},\omega_{1})$
$\displaystyle=I^{\mathrm{(NR)}}(\omega_{3},t_{2}^{\prime},\omega_{1})+I^{\mathrm{(R)}}(\omega_{3},t_{2},\omega_{1}),$
(27)

where the non-rephasing and rephrasing parts of the signal are expressed as

$\displaystyle I^{\mathrm{NR}}(\omega_{3},t_{2},\omega_{1})$
$\displaystyle=\mathrm{Im}\int_{0}^{\infty}{\rm{d}}t_{3}\int_{0}^{\infty}dt_{1}e^{i\omega_{3}t_{3}}e^{i\omega_{1}t_{1}}R^{\mathrm{(5)}}(t_{3},t_{2},t_{1}),$
(28) $\displaystyle I^{\mathrm{R}}(\omega_{3},t_{2},\omega_{1})$
$\displaystyle=\mathrm{Im}\int_{0}^{\infty}{\rm{d}}t_{3}\int_{0}^{\infty}dt_{1}e^{i\omega_{3}t_{3}}e^{-i\omega_{1}t_{1}}R^{\mathrm{(5)}}(t_{3},t_{2},t_{1}).$
(29)

In a typical system measured by 2DEVS, the frequency of electronic excitation
is much higher than vibrational modes and the signals only have off-diagonal
components. However, in our model, the energy levels of electron and proton
are similar, as shown in Figs. 2 and 3. Thus, it may not easy to excite either
electron or proton modes separately. Thus, here we assume the dipole operators
to be the summation of both electron and proton,
$\hat{\mu}_{k}=\hat{\mu}_{e}+\hat{\mu}_{p}$ for $k\geq 3$. The signals in this
measurement then be a summation of 2DES, 2DEV, and 2DVS and is refereed to as
2D resonant electronic vibrational spectroscopy (2DREVS).

Table 1: System parameters $\alpha$ | 2.0 $\rm{\AA}^{-1}$
---|---
$x_{0}$ | 1.0 $\rm{\AA}$
$Q_{0}$ | 3.0 $\rm{\AA}$
$D_{k}$ | 303435 $\rm{cm}^{-1}\rm{\AA}^{-2}$
$\Delta$ | 50 $\rm{cm}^{-1}$
$D_{LE}^{l}$ |  33715 $\rm{cm}^{-1}$
$D_{LE}^{r}$ |  31715 $\rm{cm}^{-1}$
$D_{CT}^{l}$ |  31715 $\rm{cm}^{-1}$
$D_{CT}^{r}$ |  33715 $\rm{cm}^{-1}$
Table 2: The lowest 10 energy eigenvalues of each electronic state as a unit of $\omega_{0}$.  Eigen numbers $(m,n)$ | $\left|\phi_{LE}^{(m,n)}\right\rangle$ | $\left|\phi_{CT}^{(m,n)}\right\rangle$
---|---|---
(0, 0) | 0.00 | -0.02
(0, 1) | 0.82 | 0.81
(0, 2) | 1.63 | 1.64
(1, 0) | 1.94 | 1.95
(0, 3) | 2.44 | 2.47
(1, 1) | 2.78 | 2.74
(2, 0) | 3.06 | 3.06
(0, 4) | 3.26 | 3.31
(1, 2) | 3.64 | 3.55
(2, 1) | 3.94 | 3.97

##  III Numerical Results

The system parameters chosen to simulate our PPCET model are listed in Table.
1, based on a typical PT system. Sato and Iwata (1988) The determination of
the electronic couplings $\Delta$ is a critical point of any PPCET reaction,
in that it provides the major contribution to the discrimination between
adiabatic and non-adiabatic mechanisms. Here, we choose to study the system
under moderate non-adiabatic conditions, and set
$\Delta=50~{}\mathrm{cm^{-1}}$, which is close to previously reported studies.
Hammes-Schiffer and Stuchebrukhov (2010) The energy eigenstates of the system
${\big{|}\phi_{i}^{(m,n)}\big{\rangle}}$ are obtained by diagonalizing the
matrix representation of the system Hamiltonian. The energy eigenvalues of the
lowest several states are presented in Table. 2, and a schematic view of is
given in Fig. 3. Here, the $m$ and $n$ represent the quantum numbers of the
proton and D-A modes that are determined from the number of nodes along the
$x$ and $Q$ directions. While, ${\big{|}\phi_{i}^{(m,n)}\big{\rangle}}$ with
$m=0$ mainly correspond to the charge localized states
${\big{|}\phi_{LE}^{l}\big{\rangle}}$ and
${\big{|}\phi_{CT}^{r}\big{\rangle}}$, those with $m=1$ mainly corresponded to
the intermediate transition states ${\left|\phi_{LE}^{r}\right\rangle}$ and
${\left|\phi_{CT}^{l}\right\rangle}$, respectively. The states for $m\geqq 2$
are strongly delocalized along the proton coordinate, and provide almost no
contribution to the pure PT processes. The numerical simulations of the HEOM
were conducted using the energy eigenstates representation, and we employed
the lowest 20-40 eigenstates for each electronic state based on the value of
system-bath coupling strength $\zeta$. The time integrals were carried out
using the low-storage fourth-order Runge-Kutta (LSRK4) method. The time step
was chosen as $\delta t=0.01\omega_{0}^{-1}$, where $\omega_{0}$ is a
characteristic frequency taken as the unit for all the other physical
variables. Here, we choose $\omega_{0}=500{\rm{cm^{-1}}}$. We also fixed the
inverse correlation time as $\gamma=0.5\omega_{0}$ and the bath temperature as
$\beta\hbar\omega_{0}=2.4$ (300K). The HEOM parameters required for a
converged calculation were chosen as $N=10$ and $K=5$. In the following, we
investigate the effects of the environment on the PCET mechanism as a function
of $\zeta$ by studying both the population dynamics and the TAS and 2DREVS
signals.

Figure 3: A schematic view for ${\left|\phi_{LE}\right\rangle}$ (red curve)
and ${\left|\phi_{CT}\right\rangle}$ (blue curve) in the diabatic
representation along $\hat{x}$ at the minimum of $\hat{Q}$. The lowest several
eigenstates for each PES are also plotted. The labeled orange and green arrows
represent the corresponding proton and electron transitions appearing in
nonlinear spectroscopy. See main text for the meaning of the labels.

###  III.1 Population dynamics

Figure 4: The population dynamics that represents proton and electron
localization. The red, yellow, green and blue curves represent
$P_{LE}^{l}(t)$, $P_{LE}^{r}(t)$, $P_{CT}^{l}(t)$, and $P_{CT}^{r}(t)$,
respectively, and the corresponding configurations of model system are
illustrated in Fig. 2. The population $P^{l}(t)$ and $P_{CT}(t)$ are also
presented as the dashed purple and black curves.

First, we illustrate the time evolution of the population states for various
values of $\zeta$. The electron and proton transfer rates can be estimated
from $P_{CT}(t)$ and $P^{l}(t)$. The calculated results are depicted in Fig.
4, for (a) a weak $(\zeta=0.0005\omega_{0})$, (b) moderate
$(\zeta=0.005\omega_{0})$, and (c) strong $(\zeta=0.03\omega_{0})$ coupling
cases. Note that, as illustrated in our PT investigation, Zhang, Borrelli, and
Tanimura (2020) the effective coupling strength on the present exponential-
linear system-bath coupling model is different from the conventional linear-
linear coupling models. The strength of the coupling parameter is determined
on the basis of the relaxation dynamics of the populations and spectral line
shape of TAS, as we will show below.

In the weak coupling case, Fig. 4(a), coherent recursive oscillations of state
populations are observed. These oscillations do not affect the equilibrium
distribution and do not contribute to the population transfer rates. Although
the contribution is minor, the population exchange between $P_{CT}(t)$ and
$P^{l}(t)$ suggests the presence of a charge transfer process. For the
moderate and strong coupling cases in Fig. 4(b) and (c), the linear term of
$\hat{V}_{i}$ causes the population relaxation suppressing the coherent
oscillations. In the $\hat{x}$ direction (proton mode), the nonlinear terms of
$\hat{V}_{i}$ also lead to a decrease of the energy barrier so that proton
transfer is promoted. In the $\hat{Q}$ direction, the linear term of
$\hat{V}_{i}$ leads to a decrease of proton distance from the heavy atoms and
increase the PT efficiency. A constant term $\hat{V}_{i}(0,0)$ (see Eq. (31))
is also present, corresponding to the interaction between the electronic
states and the heat bath. As a result, for larger $\zeta$, both the electron
and the proton are equally distributed in the two wells because of the
symmetric PES. Note that we cannot disentangle the contribution from the CEPT,
ET, and PT processes only from the analysis of population dynamics because
they are mixed in the population states.

###  III.2 Transient absorption spectroscopy (TAS)

Figure 5: The contour map of TAS calculated for the vibrational excitation (
$\hat{\mu}=\hat{\mu}_{p}$) in the cases of the (a) weak, (b) moderate, and (c)
strong system-bath interactions. All the peaks are normalized with respect to
maximum intensity of $t^{\prime}=10.0~{}\omega_{0}^{-1}$. The contours are
drawn from -0.5 to 0.5. The red and blue areas represent the positive
absorption and negative emission, respectively.

Next we present the results of TAS analysis. Although TAS has the capability
to analyze the populations in the ET and PT states separately following the
position of the absorption peaks, this is not easy in the present case,
because the excitation energies of the ET and PT processes are similar, and
the abortion peaks are often overlapped. Hence, here we calculated TAS for the
electronic and vibrational modes separately to help the analysis of 2DEVS. In
TAS, the charge transition rate can be evaluated from the intensity of
corresponding transition peaks, while coherent oscillation appears as a
$\delta$-function like peak. The characteristic time scale of various
transitions can also be evaluated as a function of $t^{\prime}$.

In Fig. 5 we present TAS for the vibrational excitation of the proton mode
obtained for a waiting time up to $t^{\prime}=10.0\omega_{0}^{-1}$ and by
setting $\hat{\mu}_{3}=\hat{\mu}_{4}=\hat{\mu}_{p}$. In each figure, the
negative and positive peaks represent the emission and absorption,
respectively. Note that although the energy eigenvalues of the
${\big{|}\phi_{LE}^{(m,n)}\big{\rangle}}$ and
${\big{|}\phi_{CT}^{(m,n)}\big{\rangle}}$ in the diabatic representation are
degenerate, those in the adiabatic representation are separated by the
frequency $\Delta$ because of the diabatic coupling.

In the weak coupling case, Fig. 5(a), the peak “A” $(0.05\omega_{0})$
predominantly arises from the CPET,
${\big{|}\phi_{i}^{(m,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{j}^{(m,n)}\big{\rangle}}$.
This transition always occurs due to the large overlap between two electronic
potential surfaces. The peak “B” $(0.5\omega_{0})$ and “E” $(2.7\omega_{0})$
arise from the pure PT with and without the participation of the D-A mode,
where “E” represents
${\big{|}\phi_{i}^{(0,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(1,n)}\big{\rangle}}$,
and “B” represents
${\big{|}\phi_{i}^{(0,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(1,n-2)}\big{\rangle}}$.
The peak “C” $(0.8\omega_{0})$ represents the excitation of the D-A mode,
${\big{|}\phi_{i}^{(m,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(m,n+1)}\big{\rangle}}$,
which arises because the proton and the D-A mode are strongly coupled. The
proton distribution varies as a function of the quantum number $n$ in the D-A
mode, even when the quantum number of the proton mode $m$ is unchanged. The
other three peaks represent the delocalization of the proton in the higher
energy states ($m\leqq 2$), which do not contribute to either PT or CEPT. The
peaks “D” $(1.1\omega_{0})$, “F” $(1.9\omega_{0})$, and “G” $(3.0\omega_{0})$
represent
${\big{|}\phi_{i}^{(1,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(2,n)}\big{\rangle}}$,
${\big{|}\phi_{i}^{(2,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(3,n)}\big{\rangle}}$,
and
${\big{|}\phi_{i}^{(0,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(2,n)}\big{\rangle}}$,
respectively. Most of these peaks consist of several small peaks because of
the participation of the D-A mode excited states $(n>0)$. A schematic view of
all the transitions is illustrated in Fig. 3.

We then analyze the effects of the system-bath coupling strength, $\zeta$,
through the peak intensities as a function of $t^{\prime}$. In the weak
coupling case, Fig. 5(a), most of the peaks are unchanged regardless of
$t^{\prime}$ except for the peak “A”, whose intensity changes sign near
$t^{\prime}=1.0$. In the moderate and strong coupling cases, Fig. 5(b) and
5(c), the peak intensity of “A” changes from almost 0 to a positive value in
the initial time period. This indicates that the CEPT process is promoted by
the system-bath interaction, and occurs in a relatively short time period. The
promotion effect can be explained by the linear term $V_{i,Q}^{(1)}$ in the
$\hat{Q}$ direction, which reduces the transfer distance and enhances the
vibronic coupling. The intensity of the peak “C” changes from positive to
negative values near $t^{\prime}=4.0$ in the case 5(b), and $t^{\prime}=2.0$
in the case 5(c). This implies that the characteristic time scale of the D-A
excitation is larger than the CEPT. The intensities of the PT peaks “B” and
“E” are almost unchanged, which indicates that the pure PT plays a minor role
in the present case.

Figure 6: The contour maps of TAS for the electronic excitation
($\hat{\mu}=\hat{\mu}_{e}$) in the cases of the (a) weak, (b) moderate, and
(c) strong system-bath interactions. The figure is depicted in the same way
with Fig. 5.

Finally we present TAS for electronic excitation case in Fig. 6, which was
computed by setting $\hat{\mu}_{3}=\hat{\mu}_{4}=\hat{\mu}_{e}$. In the weak
coupling case in Fig. 6(a), the peak labeled by “H” corresponds to the
transition
${\big{|}\phi_{i}^{(0,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{j}^{(1,n)}\big{\rangle}}$,
and “I” corresponds to
${\big{|}\phi_{i}^{(0,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{j}^{(1,n+1)}\big{\rangle}}$.
These two peaks represent the pure ET; a possible transition for different
$n^{\prime}$ increases for larger $\zeta_{x}$. In the moderate and strong
coupling cases, Figs. 6(b) and 6(c), the peaks “H” and “I” are significantly
broadened and enhanced because of the constant term $V_{i}^{(0)}$, which
introduces linear interactions between the electron subsystem and the heat
bath. Furthermore, several additional peaks appear in the range of
$0.0\leq\omega_{1}\leq 1.5\omega_{0}$. These peaks arise from the electronic
transitions, but their peak locations are the same in the vibrational
excitation case depicted in Fig. 5. This result can be ascribed to the strong
correlation between the electronic subsystem and the vibrational coordinates.
The increase of $\zeta$ has a promotion effect on the proton transfer, which
in turn opens additional transition pathways of the electron transfer.

We also find that most of the peaks are unchanged regardless of $t^{\prime}$
even in strong coupling case. Thus, pure ET process is not favored in all
cases for different $\zeta$ because of the pretty small electronic coupling
strength $\Delta$. With regard to the CEPT peak near $\omega=0.05\omega_{0}$,
the peak intensity changes sign in both weak and strong coupling cases, as can
be clearly seen from Figs. 6(a) and 6(c). Such variation becomes more
prominent in moderate coupling case around $t^{\prime}=2.0\omega_{0}^{-1}$, if
Fig. 6(b), which indicates that a turn-over feature under a strong enough
interaction occurs. According to the above results, we find that the CEPT is
the predominant process mostly because of the exact resonance conditions
between initial and final states.

Figure 7: The contour maps of 2DREVS for a weak coupling case, which
correspond to the case in Figs. 5(a) and 6(a). The intensities are normalized
with respect to the maximum value of each case in order to see the peak
profile, and the contour lines are drawn from -0.3 to 0.3. The red and blue
curves represent the absorption and emission, respectively. We also plot the
peaks along the diagonal line in the outside above, i.e.
$I(\omega_{1},t_{2},\omega_{1})$. The peaks in black and purple boxes
represent “ET-PT” and “CEPT” peaks, while the other outside peaks are “PT”
peaks.

###  III.3 Two-dimensional resonant electronic-vibrational spectroscopy
(2DREVS)

Next, we describe the 2DREVS signals as computed from Eq. (26). The contour
maps of the 2D correlation spectroscopy in the weak, moderate, and strong
coupling cases are illustrated in Figs. 7, 8, and 9, respectively, in which we
keep the system parameters the same used to obtain the TAS signals. Note that
most of the peaks along the diagonal line are relatively weak and not clearly
visible in contour maps. Therefore, we plot these peaks outside above as
$I^{Diag}(\omega_{1},t_{2})=I^{Corr}(\omega_{1},t_{2},\omega_{1})$.

The 2D correlation spectroscopy peak profiles in the weak coupling case is
presented in Fig. 7. For each peak, the positive intensity arises from the
stimulated emission (SE) or ground state bleaching (GSB), and the negative
intensity arises from the excited state absorption (ESA) for $n>0$. Using the
information obtained from TAS, we classify all the observed peaks into three
parts: 1\. “ET-PT” peaks represent the cross peaks in black boxes, which arise
from the coherent ET-PT processes. 2\. “CEPT” peaks represent the peaks in the
purple box, which arise from the CEPT process. 3\. “PT” peaks represent the
other peaks outside boxes, which are the vibrational cross peaks and represent
the coherence between the proton and D-A modes. The cross peaks for the ET
process are not visible because of the small electronic coupling $\Delta$.

We first discuss the “ET-PT” peaks. Most of them appear in the same
$\omega_{1}$ position as Fig. 6, and at the same $\omega_{3}$ position as Fig.
5, representing the corresponding ET-PT transitions. Here, we only concentrate
on “a1”, “a2” and “a3” that do not appear in TAS. In the $\omega_{3}$
direction, the peak “a1” represents the transition
${\big{|}\phi_{i}^{(0,n)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(0,n+2)}\big{\rangle}}$
that arises from to the SL kind of interaction in $\hat{V}_{i}$ in the
$\hat{Q}$ direction. The peaks “a2” and “a3” represent
${\big{|}\phi_{i}^{(1,1)}\big{\rangle}}\rightarrow{\big{|}\phi_{i}^{(0,5)}\big{\rangle}}$
that arise from the back PT process with a participation of the D-A mode.
These results indicate that we can analyze the combination of the ET and PT
transition from the cross peaks in 2DREVS, while these contributions are mixed
and appear as a single peak in TAS.

We now concentrate on the “CEPT” peaks. The diagonal peak “b1” arises from the
CEPT transition denoted as “A” in Fig. 3, and the other cross peaks represent
the combination of CEPT-PT and CEPT-ET, where “b2” and “b5” correspond to
“A”-“D” and “A”-“G”, and “b3” and “b4” correspond to “A”-“H” and “A”-“I”,
respectively. In addition, the peaks associated with the CEPT processes appear
at symmetric positions with diagonal line. Finally, we focus on the peak
profiles at different $t_{2}$. Both “ET-PT” and “CEPT” peaks stay unchanged
with $t_{2}$ because of the weak heat-bath effect. By contrast, the
intensities of “PT” peaks decrease when the excited proton reaches the
equilibrium distribution due to the linear interaction $V_{i,x}^{(1)}$.

Figure 8: The contour maps of 2DREVS for a moderate coupling case, which
correspond to the case in Figs. 5(b) and 6(b). The contour lines are drawn
from -0.5 to 0.5, while the other parameters are the same with Fig. 7.

In the moderate and strong coupling cases presented in Figs. 8 and 9, most of
the peaks that are related to the proton and D-A transitions are broadened
either in the $\omega_{1}$ or $\omega_{3}$ direction. We first discuss the
“ET-PT” peaks. The peak positions in the $\omega_{1}$ direction are almost
unchanged, which indicates that the system-bath interaction has a minor effect
on the pure ET, as observed in TAS. In the moderate coupling case, the peak
intensities increase in Figs. 8(ii) and 8(iii) in comparison with Fig. 8(i),
which indicates that the time scale of the corresponding transition is
relatively short. Most of the peaks almost disappear in the strong coupling
case. Such turn-over feature with $\zeta$ is also observed in TAS for
electronic excitation, but is more clear in the 2DREVS. This is because of the
suppression effect on ET-PT coherence and mainly comes from the square linear
interaction $V_{i,x}^{(2)}$.

For the “CEPT” peaks, the intensities increase with $t_{2}$, as evident from
Fig. 8, and become more apparent in Fig. 9 for a larger $\zeta$. This result
indicates the existence of a bath induced CEPT process, which is also observed
in TAS. For the “PT” cross peaks, most of them do not change till $t_{2}=1.0$,
as illustrated in Figs. 8(ii) and 9(ii), and almost vanish after a long
$t_{2}$ time as illustrated in Figs. 8(iii) and 9(iii). Also, the intensity of
the twisted positive and negative cross peak around
$(\omega_{1},\omega_{3})=(0.8\omega_{0},0.4\omega_{0})$ is reversed at
$t_{2}=10.0$. This peak mainly arises from the combination of the “B” and “C”
transitions (see Fig. 3), and the reverse indicates the relaxation of excited
proton. Thus, both proton and D-A motion have relatively longer time scales
compared to CPET and ET-PT, and they are always mixed.

Figure 9: The contour map of 2DREVS for a strong coupling case, which
correspond to the case in Figs. 5(c) and 6(c). The contour lines are drawn
from -0.5 to 0.5, while the other parameters are the same with Fig. 7.

At the end of this section, we discuss the peaks along the diagonal line,
$I^{Diag}(\omega_{1},t_{2})$, which represent the adiabatic transitions. For
all the coupling cases, these peaks occur in the same position found in the
TAS signals, representing the corresponding transitions. Among different
$t_{2}$ cases, the peaks representing the CEPT transition “A” play a major
role. The proton and D-A mode vibrations are only visible after $t_{2}=1.0$ in
Figs. 7(ii) and 7(iii). These features also corroborate the previous results
that the characteristic time scale of CEPT is shorter than proton and D-A mode
vibrations. The turn-over feature is also visible because the vibration peaks
become more apparent in Fig. 8(ii) than in Fig. 9(ii). The heat-bath plays a
minor role in ET processes so that the corresponding peaks are not visible in
all these cases. Although most of the results mentioned above are also
observed in TAS, 2D spectra allows a better understanding of each single
contribution.

##  IV Conclusion

In this paper we introduce a system-bath model in a multi-state two-
dimensional configuration space to describe the dynamics of PPCET process.
Using the HEOM in the eigenstate representation of the system, it is possible
to investigate the environment effects under a realistic system-bath
interaction that causes not only fluctuation and relaxation, but also
vibrational dephasing. Our results of population dynamics and TAS indicate
that CEPT is the predominant process and has a shorter time scale when
resonance conditions between initial and final states occur. Pure ET and PT
processes also take place at much longer time. The overall reaction would be a
summation of both concerted and sequential reaction mechanism. It is shown
that 2DREVS provides a wealth of information due to the coherence among the
excitation and detection periods. With the aids of the off-diagonal peaks, we
could detect the pathway of sequential ET-PT and PT-ET transition, and
concerted CEPT transition separately, whereas the diagonal peaks could
reproduce the results of TAS.

Although calculating nonlinear spectra is numerically intensive, 2DREVS with
TAS provides a valuable framework for studying PPCET processes. Since we use
the eigenstate representation of the system, it is also possible to improve
the description of the reacting system by increasing the dimension of its
configuration space, and by introducing a more complex and structured system-
bath interaction, for example, with the help of machine learning
approaches.Ueno and Tanimura (2020, 2021) This provides a powerful tool to
analyze the non-equilibrium reaction dynamics for rather complex PPCET
reactions.

###### Acknowledgements.

The financial support from The Kyoto University Foundation is acknowledged. RB
acknowledges the support of the University of Torino for the local research
funding Grant No. BORR-RILO-19-01.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Expansion of $\hat{V}_{i}$

In this Appendix, we expand the interaction function,
$\hat{V}_{i}(\hat{x},\hat{Q})$, with respect to $\hat{x}$ and $\hat{Q}$ up to
second order as,

$\displaystyle\hat{V}_{i}(\hat{x},\hat{Q})=V_{i}^{(0)}+V_{i,x}^{(1)}\,\hat{x}+V_{i,Q}^{(1)}\,\hat{Q}+V_{i,x}^{(2)}\,\hat{x}^{2}+V_{i,Q}^{(2)}\,\hat{Q}^{2}+V_{i,xQ}^{(2)}\,\hat{x}\hat{Q}+...$
(30)

where

$\displaystyle V_{i}^{(0)}$
$\displaystyle=\hat{V}_{i}(0,0)=2\alpha(D_{i}^{l}+D_{i}^{r})\left(\alpha
e^{\alpha x_{0}}-\alpha e^{2\alpha x_{0}}\right),$ (31) $\displaystyle
V_{i,x}^{(1)}$
$\displaystyle=\frac{\partial\hat{V}_{i}(\hat{x},\hat{Q})}{\partial\hat{x}}\bigg{|}_{(0,0)}=2\alpha(D_{i}^{l}-D_{i}^{r})\left(\alpha
e^{\alpha x_{0}}-2\alpha e^{2\alpha x_{0}}\right),$ (32) $\displaystyle
V_{i,Q}^{(1)}$
$\displaystyle=\frac{\partial\hat{V}_{i}(\hat{x},\hat{Q})}{\partial
Q}\bigg{|}_{(0,0)}=2\alpha(D_{i}^{l}+D_{i}^{r})\left(\alpha e^{\alpha
x_{0}}-2\alpha e^{2\alpha x_{0}}\right),$ (33) $\displaystyle V_{i,x}^{(2)}$
$\displaystyle=\frac{\partial^{2}\hat{V}_{i}(x,Q)}{\partial\hat{x}^{2}}\bigg{|}_{(0,0)}=2\alpha(D_{i}^{l}+D_{i}^{r})\left(\alpha^{2}e^{\alpha
x_{0}}-4\alpha^{2}e^{2\alpha x_{0}}\right),$ (34) $\displaystyle
V_{i,Q}^{(2)}$
$\displaystyle=\frac{\partial^{2}\hat{V}_{i}(\hat{x},\hat{Q})}{\partial
Q^{2}}\bigg{|}_{(0,0)}=\alpha(D_{i}^{l}+D_{i}^{r})\left(\alpha^{2}e^{\alpha
x_{0}}-4\alpha^{2}e^{2\alpha x_{0}}\right),$ (35) $\displaystyle
V_{i,xQ}^{(2)}$
$\displaystyle=\frac{\partial^{2}\hat{V}_{i}(\hat{x},\hat{Q})}{\partial\hat{x}\hat{Q}}\bigg{|}_{(0,0)}=2\alpha(D_{i}^{l}-D_{i}^{r})\left(-\alpha^{2}e^{\alpha
x_{0}}+4\alpha^{2}e^{2\alpha x_{0}}\right).$ (36)

## References

  * Migliore  _et al._ (2014) A. Migliore, N. F. Polizzi, M. J. Therien, and D. N. Beratan, “Biochemistry and theory of proton-coupled electron transfer,” Chemical Reviews 114, 3381–3465 (2014), pMID: 24684625, https://doi.org/10.1021/cr4006654 .
  * Megiatto Jr  _et al._ (2014) J. D. Megiatto Jr, D. D. Méndez-Hernández, M. E. Tejeda-Ferrari, A.-L. Teillout, M. J. Llansola-Portolés, G. Kodis, O. G. Poluektov, T. Rajh, V. Mujica, T. L. Groy, D. Gust, T. A. Moore, and A. L. Moore, “A bioinspired redox relay that mimics radical interactions of the tyr–his pairs of photosystem ii,” Nature Chemistry 6, 423–428 (2014).
  * Huynh  _et al._ (2017) M. T. Huynh, S. J. Mora, M. Villalba, M. E. Tejeda-Ferrari, P. A. Liddell, B. R. Cherry, A.-L. Teillout, C. W. Machan, C. P. Kubiak, D. Gust, T. A. Moore, S. Hammes-Schiffer, and A. L. Moore, “Concerted one-electron two-proton transfer processes in models inspired by the tyr-his couple of photosystem ii,” ACS Central Science 3, 372–380 (2017), https://doi.org/10.1021/acscentsci.7b00125 .
  * Offenbacher and Barry (2020) A. R. Offenbacher and B. A. Barry, “A proton wire mediates proton coupled electron transfer from hydroxyurea and other hydroxamic acids to tyrosyl radical in class ia ribonucleotide reductase,” The Journal of Physical Chemistry B 124, 345–354 (2020), pMID: 31904962, https://doi.org/10.1021/acs.jpcb.9b08587 .
  * Gagliardi  _et al._ (2010) C. J. Gagliardi, B. C. Westlake, C. A. Kent, J. J. Paul, J. M. Papanikolas, and T. J. Meyer, “Integrating proton coupled electron transfer (pcet) and excited states,” Coordination Chemistry Reviews 254, 2459 – 2471 (2010), 18th International Symposium on the Photochemistry and Photophysics of Coordination Compounds Sapporo, 2009.
  * Becker  _et al._ (2016) R. Becker, S. Amirjalayer, P. Li, S. Woutersen, and J. N. H. Reek, “An iron-iron hydrogenase mimic with appended electron reservoir for efficient proton reduction in aqueous media,” Science Advances 2 (2016), 10.1126/sciadv.1501014, https://advances.sciencemag.org/content/2/1/e1501014.full.pdf .
  * Yamaguchi  _et al._ (2014) A. Yamaguchi, R. Inuzuka, T. Takashima, T. Hayashi, K. Hashimoto, and R. Nakamura, “Regulating proton-coupled electron transfer for efficient water splitting by manganese oxides at neutral ph,” Nature Communications 5, 4256 (2014).
  * Allen (2002) J. F. Allen, “Photosynthesis of atp—electrons, proton pumps, rotors, and poise,” Cell 110, 273–276 (2002).
  * Papa, Capitanio, and Luca Martino (2006) S. Papa, G. Capitanio, and P. Luca Martino, “Concerted involvement of cooperative proton–electron linkage and water production in the proton pump of cytochrome c oxidase,” Biochimica et Biophysica Acta (BBA) - Bioenergetics 1757, 1133 – 1143 (2006), mitochondria: from Molecular Insight to Physiology and Pathology.
  * Cukier (1996) R. I. Cukier, “Proton-coupled electron transfer reactions: evaluation of rate constants,” The Journal of Physical Chemistry 100, 15428–15443 (1996), https://doi.org/10.1021/jp961025g .
  * Goldsmith, Soudackov, and Hammes-Schiffer (2019) Z. K. Goldsmith, A. V. Soudackov, and S. Hammes-Schiffer, “Theoretical analysis of the inverted region in photoinduced proton-coupled electron transfer,” Faraday Discuss. 216, 363–378 (2019).
  * Mayer (2011) J. M. Mayer, “Simple marcus-theory-type model for hydrogen-atom transfer/proton-coupled electron transfer,” The Journal of Physical Chemistry Letters 2, 1481–1489 (2011), pMID: 21686056, https://doi.org/10.1021/jz200021y .
  * Soudackov and Hammes-Schiffer (2016) A. V. Soudackov and S. Hammes-Schiffer, “Proton-coupled electron transfer reactions: analytical rate constants and case study of kinetic isotope effects in lipoxygenase,” Faraday Discuss. 195, 171–189 (2016).
  * Soudackov and Hammes-Schiffer (2000) A. Soudackov and S. Hammes-Schiffer, “Derivation of rate expressions for nonadiabatic proton-coupled electron transfer reactions in solution,” The Journal of Chemical Physics 113, 2385–2396 (2000), https://doi.org/10.1063/1.482053 .
  * Soudackov and Hammes-Schiffer (2015) A. V. Soudackov and S. Hammes-Schiffer, “Nonadiabatic rate constants for proton transfer and proton-coupled electron transfer reactions in solution: Effects of quadratic term in the vibronic coupling expansion,” The Journal of Chemical Physics 143, 194101 (2015), https://doi.org/10.1063/1.4935045 .
  * Cukier (1994) R. I. Cukier, “Mechanism for proton-coupled electron-transfer reactions,” The Journal of Physical Chemistry 98, 2377–2381 (1994), https://doi.org/10.1021/j100060a027 .
  * Hazra, Soudackov, and Hammes-Schiffer (2010) A. Hazra, A. V. Soudackov, and S. Hammes-Schiffer, “Role of solvent dynamics in ultrafast photoinduced proton-coupled electron transfer reactions in solution,” The Journal of Physical Chemistry B 114, 12319–12332 (2010), pMID: 20809583, https://doi.org/10.1021/jp1051547 .
  * Kretchmer and Miller (2016) J. S. Kretchmer and T. F. Miller, “Tipping the balance between concerted versus sequential proton-coupled electron transfer,” Inorganic Chemistry 55, 1022–1031 (2016), pMID: 26440812, https://doi.org/10.1021/acs.inorgchem.5b01821 .
  * Kobrak and Hammes-Schiffer (2001) M. N. Kobrak and S. Hammes-Schiffer, “Molecular dynamics simulation of proton-coupled electron transfer in solution,” The Journal of Physical Chemistry B 105, 10435–10445 (2001), https://doi.org/10.1021/jp012102x .
  * Grimminger and Schmickler (2007) J. Grimminger and W. Schmickler, “Dynamics of combined electron- and proton transfer at metal electrodes,” Chemical Physics 334, 8 – 17 (2007).
  * Auer, Fernandez, and Hammes-Schiffer (2011) B. Auer, L. E. Fernandez, and S. Hammes-Schiffer, “Theoretical analysis of proton relays in electrochemical proton-coupled electron transfer,” Journal of the American Chemical Society 133, 8282–8292 (2011), pMID: 21524104, https://doi.org/10.1021/ja201560v .
  * Peluso, Di Donato, and Saracino (2000) A. Peluso, M. Di Donato, and G. A. A. Saracino, “An alternative way of thinking about electron transfer in proteins: Proton assisted electron transfer between the primary and the secondary quinones in photosynthetic reaction centers,” The Journal of Chemical Physics 113, 3212–3218 (2000), https://doi.org/10.1063/1.1286918 .
  * Cho (2019) M. Cho,  _Coherent Multidimensional Spectroscopy_ (Springer, Singapore, 2019).
  * Pfeiffer  _et al._ (1999) M. Pfeiffer, C. Chudoba, A. Lau, K. Lenz, and T. Elsaesser, “Coherent motion of low frequency vibrations in ultrafast excited state proton transfer,” Laser Chemistry 19, 101–103 (1999).
  * Hogle, Cunningham, and Tucker (2018) D. G. Hogle, A. R. Cunningham, and M. J. Tucker, “Equilibrium versus nonequilibrium peptide dynamics: Insights into transient 2d ir spectroscopy,” The Journal of Physical Chemistry B 122, 8783–8795 (2018), pMID: 30040900, https://doi.org/10.1021/acs.jpcb.8b05063 .
  * Kuss-Petermann  _et al._ (2012) M. Kuss-Petermann, H. Wolf, D. Stalke, and O. S. Wenger, “Influence of donor–acceptor distance variation on photoinduced electron and proton transfer in rhenium(i)–phenol dyads,” Journal of the American Chemical Society 134, 12844–12854 (2012), pMID: 22809316, https://doi.org/10.1021/ja3053046 .
  * Giret  _et al._ (2020) Y. Giret, P. Guo, L.-F. Wang, and J. Cheng, “Theoretical study of kinetics of proton coupled electron transfer in photocatalysis,” The Journal of Chemical Physics 152, 124705 (2020), https://doi.org/10.1063/5.0001825 .
  * Zheng  _et al._ (2011) R. Zheng, Y. Jing, L. Chen, and Q. Shi, “Theory of proton coupled electron transfer reactions: Assessing the born–oppenheimer approximation for the proton motion using an analytically solvable model,” Chemical Physics 379, 39 – 45 (2011).
  * Song and Shi (2017) K. Song and Q. Shi, “Theoretical study of photoinduced proton coupled electron transfer reaction using the non-perturbative hierarchical equations of motion method,” The Journal of Chemical Physics 146, 184108 (2017), https://doi.org/10.1063/1.4982928 .
  * Barbara, Walker, and Smith (1992) P. F. Barbara, G. C. Walker, and T. P. Smith, “Vibrational modes and the dynamic solvent effect in electron and proton transfer,” Science 256, 975–981 (1992), https://science.sciencemag.org/content/256/5059/975.full.pdf .
  * Tanimura and Mukamel (1993) Y. Tanimura and S. Mukamel, “Two‐dimensional femtosecond vibrational spectroscopy of liquids,” The Journal of Chemical Physics 99, 9496–9511 (1993), https://doi.org/10.1063/1.465484 .
  * Mukamel (1995) S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_ (Oxford University Press, New York, 1995).
  * Hamm and Zanni (2011) P. Hamm and M. Zanni,  _Concepts and Methods of 2D Infrared Spectroscopy_ (Cambridge University Press, 2011).
  * Gaynor  _et al._ (2016) J. D. Gaynor, T. L. Courtney, M. Balasubramanian, and M. Khalil, “Fourier transform two-dimensional electronic-vibrational spectroscopy using an octave-spanning mid-ir probe,” Opt. Lett. 41, 2895–2898 (2016).
  * Oliver (2018) T. A. A. Oliver, “Recent advances in multidimensional ultrafast spectroscopy,” Royal Society Open Science 5, 171425 (2018), https://royalsocietypublishing.org/doi/pdf/10.1098/rsos.171425 .
  * Fox, Blair, and Khalil (2020) Z. W. Fox, T. J. Blair, and M. Khalil, “Determining the orientation and vibronic couplings between electronic and vibrational coordinates with polarization-selective two-dimensional vibrational-electronic spectroscopy,” The Journal of Physical Chemistry Letters 11, 1558–1563 (2020), pMID: 32004009, https://doi.org/10.1021/acs.jpclett.9b03752 .
  * Mukamel (2000) S. Mukamel, “Multidimensional femtosecond correlation spectroscopies of electronic and vibrational excitations,” Annual Review of Physical Chemistry 51, 691–729 (2000), pMID: 11031297, https://doi.org/10.1146/annurev.physchem.51.1.691 .
  * Dong  _et al._ (2015) H. Dong, N. H. C. Lewis, T. A. A. Oliver, and G. R. Fleming, “Determining the static electronic and vibrational energy correlations via two-dimensional electronic-vibrational spectroscopy,” The Journal of Chemical Physics 142, 174201 (2015), https://doi.org/10.1063/1.4919684 .
  * Wu  _et al._ (2019a) E. C. Wu, Q. Ge, E. A. Arsenault, N. H. C. Lewis, N. L. Gruenke, M. J. Head-Gordon, and G. R. Fleming, “Two-dimensional electronic-vibrational spectroscopic study of conical intersection dynamics: an experimental and electronic structure study,” Phys. Chem. Chem. Phys. 21, 14153–14163 (2019a).
  * Ikeda and Tanimura (2018) T. Ikeda and Y. Tanimura, “Phase-space wavepacket dynamics of internal conversion via conical intersection: Multi-state quantum fokker-planck equation approach,” Chemical Physics 515, 203 – 213 (2018), ultrafast Photoinduced Processes in Polyatomic Molecules:Electronic Structure, Dynamics and Spectroscopy (Dedicated to Wolfgang Domcke on the occasion of his 70th birthday).
  * Lewis and Fleming (2016) N. H. C. Lewis and G. R. Fleming, “Two-dimensional electronic-vibrational spectroscopy of chlorophyll a and b,” The Journal of Physical Chemistry Letters 7, 831–837 (2016), pMID: 26894783, https://doi.org/10.1021/acs.jpclett.6b00037 .
  * Courtney  _et al._ (2015) T. L. Courtney, Z. W. Fox, K. M. Slenkamp, and M. Khalil, “Two-dimensional vibrational-electronic spectroscopy,” The Journal of Chemical Physics 143, 154201 (2015), https://doi.org/10.1063/1.4932983 .
  * Wu  _et al._ (2019b) E. C. Wu, E. A. Arsenault, P. Bhattacharyya, N. H. C. Lewis, and G. R. Fleming, “Two-dimensional electronic vibrational spectroscopy and ultrafast excitonic and vibronic photosynthetic energy transfer,” Faraday Discuss. 216, 116–132 (2019b).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath,” Journal of the Physical Society of Japan 58, 101–114 (1989), https://doi.org/10.1143/JPSJ.58.101 .
  * Tanimura (1990) Y. Tanimura, “Nonperturbative expansion method for a quantum system coupled to a harmonic-oscillator bath,” Phys. Rev. A 41, 6676–6687 (1990).
  * Ishizaki and Tanimura (2005) A. Ishizaki and Y. Tanimura, “Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach,” Journal of the Physical Society of Japan 74, 3131–3134 (2005), https://doi.org/10.1143/JPSJ.74.3131 .
  * Tanimura (2006) Y. Tanimura, “Stochastic liouville, langevin, fokker–planck, and master equation approaches to quantum dissipative systems,” Journal of the Physical Society of Japan 75, 082001 (2006), https://doi.org/10.1143/JPSJ.75.082001 .
  * Tanimura (2014) Y. Tanimura, “Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities,” The Journal of Chemical Physics 141, 044114 (2014), https://doi.org/10.1063/1.4890441 .
  * Tanimura (2015) Y. Tanimura, “Real-time and imaginary-time quantum hierarchal fokker-planck equations,” The Journal of Chemical Physics 142, 144110 (2015), https://doi.org/10.1063/1.4916647 .
  * Tanimura (2020) Y. Tanimura, “Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (heom),” The Journal of Chemical Physics 153, 020901 (2020), https://doi.org/10.1063/5.0011599 .
  * Zhang, Borrelli, and Tanimura (2020) J. Zhang, R. Borrelli, and Y. Tanimura, “Proton tunneling in a two-dimensional potential energy surface with a non-linear system–bath interaction: Thermal suppression of reaction rate,” The Journal of Chemical Physics 152, 214114 (2020), https://doi.org/10.1063/5.0010580 .
  * Tew, Handy, and Carter (2006) D. P. Tew, N. C. Handy, and S. Carter, “A reaction surface hamiltonian study of malonaldehyde,” The Journal of Chemical Physics 125, 084313 (2006), https://doi.org/10.1063/1.2338891 .
  * Tanimura and Maruyama (1997) Y. Tanimura and Y. Maruyama, “Gaussian–markovian quantum fokker–planck approach to nonlinear spectroscopy of a displaced morse potentials system: Dissociation, predissociation, and optical stark effects,” The Journal of Chemical Physics 107, 1779–1793 (1997), https://doi.org/10.1063/1.474531 .
  * Ikeda and Tanimura (2017) T. Ikeda and Y. Tanimura, “Probing photoisomerization processes by means of multi-dimensional electronic spectroscopy: The multi-state quantum hierarchical fokker-planck equation approach,” The Journal of Chemical Physics 147, 014102 (2017), https://doi.org/10.1063/1.4989537 .
  * Ikeda, Dijkstra, and Tanimura (2019) T. Ikeda, A. G. Dijkstra, and Y. Tanimura, “Modeling and analyzing a photo-driven molecular motor system: Ratchet dynamics and non-linear optical spectra,” The Journal of Chemical Physics 150, 114103 (2019), https://doi.org/10.1063/1.5086948 .
  * Hu  _et al._ (2011) J. Hu, M. Luo, F. Jiang, R.-X. Xu, and Y. Yan, “Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems,” The Journal of Chemical Physics 134, 244106 (2011), https://doi.org/10.1063/1.3602466 .
  * Ding, Xu, and Yan (2012) J.-J. Ding, R.-X. Xu, and Y. Yan, “Optimizing hierarchical equations of motion for quantum dissipation and quantifying quantum bath effects on quantum transfer mechanisms,” The Journal of Chemical Physics 136, 224103 (2012), https://doi.org/10.1063/1.4724193 .
  * Shi  _et al._ (2009) Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan, “Efficient hierarchical liouville space propagator to quantum dissipative dynamics,” The Journal of Chemical Physics 130, 084105 (2009), https://doi.org/10.1063/1.3077918 .
  * Ikeda and Tanimura (2019) T. Ikeda and Y. Tanimura, “Low-temperature quantum fokker–planck and smoluchowski equations and their extension to multistate systems,” Journal of Chemical Theory and Computation 15, 2517–2534 (2019), pMID: 30776312, https://doi.org/10.1021/acs.jctc.8b01195 .
  * Sato and Iwata (1988) N. Sato and S. Iwata, “Promotion of the proton transfer reaction by the intermolecular stretching mode: Application of the two‐dimensional finite element method to the nuclear schrödinger equation,” The Journal of Chemical Physics 89, 2932–2937 (1988), https://doi.org/10.1063/1.454998 .
  * Hammes-Schiffer and Stuchebrukhov (2010) S. Hammes-Schiffer and A. A. Stuchebrukhov, “Theory of coupled electron and proton transfer reactions,” Chemical Reviews 110, 6939–6960 (2010), pMID: 21049940, https://doi.org/10.1021/cr1001436 .
  * Ueno and Tanimura (2020) S. Ueno and Y. Tanimura, “Modeling intermolecular and intramolecular modes of liquid water using multiple heat baths: Machine learning approach,” J. Chem. Theory Comput. 16, 2099–2108 (2020).
  * Ueno and Tanimura (2021) S. Ueno and Y. Tanimura, “Modeling and simulating excited-state dynamics of a system in condensed phases: Machine learning approach,” arXiv https://arxiv.org/abs/2102.02427 xx, xxxx–xxxx (2021).


# Decisions & Disruptions 2: Decide Harder
A custom cyber security incident response exercise

Benjamin Shreeve
University of Bristol  Joseph Gardiner
University of Bristol  Joseph Hallett
University of Bristol  David Humphries
City of London Police  Awais Rashid
University of Bristol

###### Abstract

Cyber incident response is critical to business continuity—we describe a new
exercise that challenges professionals to play the role of Chief Information
Security Officer (CISO) for a major financial organisation. Teams must decide
how organisational team and budget resources should be deployed across
Enterprise Architecture (EA) upgrades and cyber incidents. Every choice made
has an impact—some prevent whilst others may trigger new or continue current
attacks. We explain how the underlying platform supports these interactions
through a reactionary event mechanism that introduces events based on the
current attack surface of the organisation. We explore how our platform
manages to introduce randomness on top of triggered events to ensure that the
exercise is not deterministic and better matches incidents in the _real
world_. We conclude by describing next steps for the exercise and how we plan
to use it in the future to better understand risk decision making.

Copyright is held by the author/owner. Permission to make digital or hard
copies of all or part of this work for personal or classroom use is granted
without fee.

_Workshop on Deconstructing Gamified Approaches to Security and Privacy,
located at USENIX Symposium on Usable Privacy and Security (SOUPS)_ 2023\.
August 6–8, 2023, Anaheim, CA, United States.

##  1 Introduction

Major cyber security incidents regularly disrupt businesses, and in extreme
circumstance have even bankrupted them. We have created a major new incident
response exercise to help businesses. We have worked with specialists from law
enforcement and major financial organisations to create an exercise that
challenges teams to handle a major, responsive cyber security incident. The
aim of the exercise is to expose senior managers and incident response teams
to the time, resource and political pressures they will encounter whilst
handling a major crisis whilst at the same time gathering granular decision-
making data to inform cyber response handling research.

This work builds upon a freely available previous game released under a CC-BY-
NC license: Decisions & Disruptions (D-D) 111decisions-disruptions.org. D-D is
a highly successful tabletop exercise utilised by police forces across the UK
and businesses across the world. It was designed to explore how people make
risk decisions around cyber-physical infrastructures. Whilst D-D has provided
valuable practical and research insights [5, 13, 14, 12], we argue that it is
limited by a deterministic mechanism, fixed and tied to one sector. We,
therefore, propose a new game: Decisions & Disruptions 2: Decide Harder (D-D
2) which provides an extensible engine for risk decision making exercises that
incorporates randomness and more complicated threat relationships that can be
targeted towards any industry, rather than just for critical national
infrastructure. This paper introduces our proposed new game D-D 2 and the new
features it incorporates.

##  2 Related Work

Cyber security exercises are a popular way of raising awareness of the subject
matter. A number of exercises have been developed specifically as part of
University courses (e.g, [3, 10]). However, such exercises tend to have a
relatively narrow scope related to the content of specific University courses
and they are rarely validated with or used outside of academia.

Other existing exercises have been created to help raise awareness of cyber
security issues in industry (e.g., [4, 6, 11, 5, 2, 8]). All of these
exercises are tabletop exercises and with the exception of Frey et al. [5] all
use card-based mechanisms. Whilst the card-based mechanisms are valuable for
raising awareness they often limit how well exercises can reflect real-world
scenarios. For example, such mechanisms rely on a mixed deck of cards (or
several decks) which are then drawn from at random to provide game events. As
such it is hard to such mechanisms to capture the way that events in the real-
world may be related with one event causing another to occur. This is not to
say they are not of value, but they often emphasise learning of specific
aspects rather than emulation of scenarios. For example, Hart et al. [8]
introduce the Riskio serious game—their exercise is aimed at non-technical
participants and challenges them to consider potential threat vectors and then
identify possible countermeasures. This provides a no doubt valuable learning
experience, but is not the same as exposing participants to an emulation of
decision-making under pressure which we is the aim of our exercise. Hart et
al. [7] have taken lessons learned from the development of Riskio [8] and used
it, along with a careful analysis of a wide range of cyber security games
(including D-D), to create the MOTENS design model for serious games. This
model suggests that the most effective cyber security games include: Multiple
modes of learning—exposing players to a wide range of cyber security aspects;
Ownership Self-Learning—Providing a range of options to help meeting learning
objectives;Theory—that supports the design;Environment—creating an appropriate
environment where people feel they can learn;Negotiation—moving toward a
coaching and problem-based learning style;Self learning—enabling participants
to build upon their base knowledge.

Frey et al. [5] provide a different approach—their exercise _Decisions &
Disruptions_ provides teams with a Lego representation of a hydro-electric
company with an plant (operations) site and separate office site. Teams play
through 4 rounds, investing a budget each round to implement security controls
on these sites. At the end of each round they discover what cyber attacks have
befallen the business as a result of their choices. The game mechanism makes
it clear to players that there is a direct correlation between the investment
choices they have made and the events they have suffered. This is a
significant improvement in terms of realism, enabling more complex attack
scenarios to be developed. However, there are limitations to the mechanism
created by Frey et al. [5], the game mechanism has hard-coded paths through it
(e.g., event A will occur if security control B has not been purchased by
round 3). This means that teams cannot replay the exercise as they will always
experience the same set consequences for their actions.

We seek to build on the success of D-D by developing an exercise where the
players affect not only the landscape in which they play but also the
consequences of their choices. An exercise where events that have occurred and
choices that have been made increase the likelihood of (or even directly
cause) specific events to occur, and where events in turn can change the
landscape. We work closely with CISOs and Cyber Security Incident Response
Teams to develop a unique, replayable exercise, inspired by D-D that
challenges teams to make cyber security decisions under the time, resource and
political pressures of a series of unfolding cyber incidents.

##  3 The Exercise

The exercise challenges teams to help the CISO of a fictitious financial
organisation handle a _very bad_ cyber security week—each day of the week is
represented by a 20 minute time-limited round. Each day, a series of cyber-
related events will occur, which can be handled in various ways. Players have
to decide which tasks the CISO’s team will tackle each day to stay on top of
what’s happening to their company.

Teams have the ability to affect the overall attack-surface of the business
and therein the type and effectiveness of possible attack vectors used against
the business by updating assets through actions such as patching and staff
training. How attacks (or symptoms) are handled affect whether attack(s) are
prevented, continue or whether new related attacks are triggered. The exercise
is designed to be played by CSIRT, CISOs and senior executives.

###  3.1 A Reactionary Event Mechanism

In order to be able to represent the complexities of real-world decision-
making a complex and reactionary event mechanism was developed (see figure 1).
Teams are tasked with protecting the EA of the business (see figure 2). They
are able to affect the company’s attack surface by investing staff hours each
day into a wide range of possible upgrades (and in some cases by sacrificing
profit). The attack surface is evaluated at the start of each round and used
to identify which events can occur (out of a library of 120 possible events).
Of those that can occur, 5 are randomly selected and added to the event list
for the next round. Events can be specified to only occur given a particular
set of criteria—such as a particular EA upgrade having been purchased or a
previous event having occurred. Each choice that a team makes has an
associated cost in terms of staff hours and company profitability. Events have
been carefully designed in tandem with cyber security law enforcement officers
to provide a realistic representation of the threats facing financial
organisations. The range of events reflect the major MITRE ATT&CK
222https://attack.mitre.org/matrices/enterprise/ areas highlighted through
interviews with Chief Information Security Officers (CISOs) in major financial
organisations as problematic.

For many choices there are also consequences. These are either explicit
feedback to teams as to the impact of their choices—including additional
impact to hours/profitability/shareprice—or an in-game consequence which can
include triggering other events in the future. These performance indicators
were flagged up through interviews with CISOs and board members of major
financial organisations as vital indications of performance during a cyber
incident [15, 1]. These triggered events are then added to the event list for
the next round (possibly bypassing any criteria evaluation). Events can also
trigger other events to occur if they are ignored; for example if a team
decides an event is not important in one round it can still have an impact on
their next round by queuing up an event to penalise their neglect. Some events
are designated _‘on-draw’_ events—they affect the round as soon as they are
drawn, deliberately introducing variance to the game. For example, an event
may tell a team that a member of the CSIRT is ill that day; they therefore
will start the round 8 hours down, with only 72 hours available to them to
utilise. These _‘on-draw’_ events can force hours, profitability, shareprice
and even make EA upgrade choices for teams before a round starts.

Figure 1: Summary of the different mechanisms in the game and the outcomes
that can happen.

For example, as part of the game an event may occur where a staff member
reports finding a USB stick on the ground. The CISO can decide whether to
ignore it, to forensically analyse it (at a cost of time) or to destroy the
device. Each decision will have impacts (and could in turn lead to more events
occurring. To visualise these kinds of decisions we created an automated tool
to create attack trees from the engine’s database of events and explore what
happens. This gives a quick pictorial guide to the consequences of any
decision in the game and the threat landscape of any particular configuration.

Figure 2: Financial Organisation enterprise architecture

###  3.2 Resources

Teams have to negotiate various resource constraints. Firstly, they have to
contend with the challenge that a day has a finite length—in this case a
configurable 20 minute limit. If they fail to utilise their resources
effectively in that time then the game automatically moves forward to the next
day. The primary unit of daily resources are the number of _‘hours’_ available
to the team each day—that is, how many of staff in the CSIRT team that they
are managing are available that day (the exercise assumes that there are 10
staff members, each with a maximum 8 hours per day, 80 hours in total per
day). These hours are then allocated to event/EA actions.

The businesses performance is represented by two broad financial metrics that
are both affected by choices made and consequences. Firstly, the team have to
manage the overall long-term _‘projected profit’_ for the business—this starts
at £500,000. However, it can be reduced through loss of business or fines that
occur as a result of choices made. Teams can also choose to reinvest some of
this projected profit into the business through certain EA upgrades or event
choices which have associated financial costs. Secondly, teams have to
consider the short-term performance of the business in the form of its _‘share
price’_ —which can be affected both positively and negatively as a result of
choices made from an initial value of £100. If the players actions result in
the _‘projected profit’_ or _‘share price’_ reaching zero then the game ends.

###  3.3 Exercise Management

The exercise comprises two parts: a game master tool that tracks and evaluates
the state of play and an associated physical card deck (see figure 3). The
game master utilises a digital tool, written in Java, that reads a database
containing details on the possible events (and their relationship to one
another) as well as the assets and upgrades that possible for the EA. Each
event/choices and asset/upgrade are also printed to a deck of cards providing
session flexibility.

###  3.4 Typical Play Through

At the start of each _day_ (round) the exercise interface updates to reflect
the status of the business (see example in appendix B).

This right hand panel consists of a list of possible actions that can be
undertaken to help improve the cyber security of their assets. The left hand
panel shows a list of the most pertinent cyber events that have hit the
business in the last 24 hours and that need resolving (a minimum of 5 occur
each day) plus any that have yet to be actioned in some way that are remaining
from previous rounds.

The EA upgrades (see example in appendix B) can be purchased at any point
during a round and affect the attack surface of the game in the next round.

Each event presents presents teams with several choices about how they might
choose to resolve each of these events (see example in appendix B). All events
by default include ‘ignore this round’ which defers making a choice for that
event to the next round. All choices are final—once a choice has been made
(other than ignore), then it cannot be undone.

Teams work together to identify which events and EA upgrades they should
dedicate resources towards in the round. As they make each choice the game
provides them with instant feedback which may result in their hours,
profitability or share price being immediately hit, which may in turn then
limit the choices they had planned. When the round concludes the game utilises
the reactionary event mechanism logic (see section 3.1 and figure 1) to
evaluate which events will occur in the next round (as a result of changes the
team have made to the businesses attack surface), or that the team have
directly triggered (because of choices made in during the round).

The exercise continues this way through multiple rounds with teams starting
each round with at least 5 events to address and a starting number of hours
for the day affected by ‘on-draw’ events. The exercise concludes either when
the round limit is reached (7 rounds) or when one of the failure criteria are
triggered. The game has 3 failure criteria: the profitability of the business
can reach zero, the share price can reach zero or the team can have 10 or more
concurrent events open. We consider a team that has failed to resolve (or has
triggered) 10 (or more) events to have reached a point of event saturation
that could never be resolved in the real world. The extent to which
consequences of event choices affect the profitability and shareprice of the
organisation vary based on the magnitude of the event. Some consequences may
have a positive impact whilst others may be severe, or even sufficient to
bankrupt the organisation.

##  4 Design Choices

We have created a mechanism that can identify if an attack is viable based
upon not only the current state of the organisations attack surface, but also
on whether specific previous events have occurred and particular choices made.
This means that events that are presented are far more representative of real-
world scenarios.

We have taken this further by treating certain defensive approaches as
_‘resistances’_. Defensive mechanisms, like _firewalls_ and _antivirus_ , can
never be 100% effective. Instead, our game mechanism captures the current
success rate of these mechanisms (for example the Firewall is only 60%
effective) and then uses these values to establish if an attack that could be
prevented by a firewall will occur by incorporating chance. We also provide
teams with the ability to improve these resistances by investing in specific
EA upgrades—and punish _negligence_ for not doing so in a timely manner. In
doing so we create an exercise where there is no longer a 1:1 relationship
between events and defences and instead there are combinations of defences
that can limit the likelihood of specific attack vectors being exploited—just
like in the real world. The tool is itself also extensible. The events for a
given game are taken from a database and can be rewritten for different
sectors or markets. Penalties for different events and how events interlink is
configurable allowing for a wide range of variations of the game to be created
with relative ease.

Our design choices fulfil Hart et al.’s [7] MOTENS pedagogical design
framework for serious cyber games: Multiple Modes of Learning: The mechanism
enables participants to experience not only a range of events, but also to
explore the relationship between EA and system security and attacks. Ownership
and Self-Learning: The underlying platform is flexible, allowing sessions to
be tailored to specific learning objectives. Theory: The overarching game
principle is informed by experiential learning theory [9]. Environment: The
exercise environment is designed so teams have more forgiving opening rounds
to help them familiarise themselves with the exercise mechanics and
expectations. Negotiation: The promotion of shared decision-making amongst
participants and the immediacy of feedback through the reactionary event
mechanism move learning from static presentation to an immersive exploratory
experience. Finally, Self-Learning: the facilitation of these exercises
enables participants to ask for clarifications when needed, enabling each
participant to start from their own knowledge baseline.

##  5 Next steps and Future Work

The exercise is ready for testing with real-world CSIRTs and executives. This
stage will be used to identify UI bugs, and refine the content of the event
and asset database. Once this is complete the exercise will be made freely
available under a Creative Commons License for organisations to use. We will
continue to work with our industry partners using the exercise to gather data
around how organisations prioritise and respond to different incidents and
triggers.

Future work will focus on the creation of new asset and event databases for
different sectors, making it possible to explore how different sectors and
different technology stacks handle incident scenarios. Decisions and
Disruptions has always helped organisations better understand how they make
risk decisions: but now they have to _decide harder_.

## Acknowledgments

We would like to thank Cyber Griffin, the City of London Police and the City
of London Corporation for their funding and support.

## References

  * [1] Ioannis Agrafiotis, Jason R C Nurse, Michael Goldsmith, Sadie Creese, and David Upton.  A taxonomy of cyber-harms: Defining the impacts of cyber-attacks and understanding how they propagate.  Journal of Cybersecurity, 4(1), 10 2018.  tyy006.
  * [2] Kristian Beckers and Sebastian Pape.  A serious game for eliciting social engineering security requirements.  In 2016 IEEE 24th International Requirements Engineering Conference (RE), pages 16–25. IEEE, 2016.
  * [3] Kevin Bock, George Hughey, and Dave Levin.  King of the hill: A novel cybersecurity competition for teaching penetration testing.  In 2018 USENIX Workshop on Advances in Security Education (ASE 18). USENIX Association, 2018.
  * [4] Tamara Denning, Adam Lerner, Adam Shostack, and Tadayoshi Kohno.  Control-alt-hack: The design and evaluation of a card game for computer security awareness and education.  In CCS, pages 915–928. ACM, 2013.
  * [5] Sylvain Frey, Awais Rashid, Pauline Anthonysamy, Maria Pinto-Albuquerque, and Syed Asad Naqvi.  The good, the bad and the ugly: A study of security decisions in a cyber-physical systems game.  IEEE Transactions on Software Engineering, 45(5):521–536, 2019\.
  * [6] Mark Gondree and Zachary NJ Peterson.  Valuing security by getting [d0x3d!]: Experiences with a network security board game.  In Presented as part of the 6th Workshop on Cyber Security Experimentation and Test. USENIX, 2013.
  * [7] Stephen Hart, Basel Halak, and Vladimiro Sassone.  Motens: A pedagogical design model for serious cyber games, 2021.
  * [8] Stephen Hart, Andrea Margheri, Federica Paci, and Vladimiro Sassone.  Riskio: A serious game for cyber security awareness and education.  Computers & Security, 95:101827, 2020.
  * [9] David A Kolb.  Experiential learning: Experience as the source of learning and development.  FT press, 2014.
  * [10] John R. Morelock and Zachary Peterson.  Authenticity, ethicality, and motivation: A formal evaluation of a 10-week computer security alternate reality game for CS undergraduates.  In Wu-chang Feng and Ashley L. Podhradsky, editors, 2018 USENIX Workshop on Advances in Security Education, ASE 2018, Baltimore, MD, USA, August 13, 2018. USENIX Association, 2018.
  * [11] Adam Shostack.  Elevation of privilege: Drawing developers into threat modeling.  In 2014 USENIX Summit on Gaming, Games, and Gamification in Security Education (3GSE 14), 2014.
  * [12] Benjamin Shreeve, Catarina Gralha, Awais Rashid, João Araujo, and Miguel Goulão.  Making sense of the unknown: How managers make cyber security decisions.  ACM Trans. Softw. Eng. Methodol., aug 2022.  Just Accepted.
  * [13] Benjamin Shreeve, Joseph Hallett, Matthew Edwards, Pauline Anthonysamy, Sylvain Frey, and Awais Rashid.  “So if Mr Blue Head here clicks the link…” Risk thinking in cyber security decision making.  Transactions on Privacy and Security, 24(1):29, November 2020.
  * [14] Benjamin Shreeve, Joseph Hallett, Matthew Edwards, Kopo Marvin Ramokapane, Richard Atkins, and Awais Rashid.  The best laid plans or lack thereof: Security decision-making of different stakeholder groups.  IEEE Transactions on Software Engineering, 48(5):1515 – 1528, 2020\.
  * [15] Onur Kemal Tosun.  Cyber-attacks and stock market activity.  International Review of Financial Analysis, 76:101795, 2021.

##  Appendix A Example Physical Cards

Figure 3: Examples of physical cards used to play the game when the UI is
unavailable.

##  Appendix B Exercise UI

Figure 4: Main exercise interface.

Figure 5: Examples of an event and an EA.

  *[CISO]: Chief Information Security Officer
  *[EA]: Enterprise Architecture
  *[D-D]: Decisions & Disruptions
  *[D-D 2]: Decisions & Disruptions 2: Decide Harder
  *[CISOs]: Chief Information Security Officer
  *[CSIRT]: Cyber Security Incident Response Team
  *[CSIRTs]: Cyber Security Incident Response Team