#  Coherent control from quantum committment probabilities

Michelle C. Anderson  Department of Chemistry, University of California,
Berkeley  Amro Dodin  Department of Chemistry, University of California,
Berkeley  Chemical Sciences Division, Lawrence Berkeley National Laboratory
Thomas P. Fay  Department of Chemistry, University of California, Berkeley
David T. Limmer<EMAIL_ADDRESS>Department of Chemistry, University of
California, Berkeley  Kavli Energy NanoSciences Institute, University of
California, Berkeley  Chemical Sciences Division, Lawrence Berkeley National
Laboratory  Materials Sciences Division, Lawrence Berkeley National Laboratory

###### Abstract

We introduce a general definition of a quantum committor in order to clarify
reaction mechanisms and facilitate control in processes where coherent effects
are important. With a quantum committor, we generalize the notion of a
transition state to quantum superpositions and quantify the effect of
interference on the progress of the reaction. The formalism is applicable to
any linear quantum master equation supporting metastability for which
absorbing boundary conditions designating the reactant and product states can
be applied. We use this formalism to determine the dependence of the quantum
transition state on coherences in a polaritonic system and optimize the
initialization state of a conical intersection model to control reactive
outcomes, achieving yields of the desired state approaching 100%. In addition
to providing a practical tool, the quantum committor provides a conceptual
framework for understanding reactions in cases when classical intuitions fail.

††preprint: 1

## Introduction

Understanding mechanisms of reactions evolving quantum mechanically is
difficult, as delocalization, tunneling and interference preclude equating a
reaction with the classical idea of motion over a barrier. The committor,
which is the probability of the system to complete a reaction,1; 2 defines an
ideal reaction coordinate, and is used to identify transition states and
mechanisms in classical reactions. Previous work has recently extended the
notion of the committor to systems in which transition states are delocalized
and the dynamics are non-adiabatic,3; 4 but has stopped short of dealing with
the dynamical consequences of coherences, due to the complications that arise
in defining trajectories in such cases.5 Here we overcome these previous
limitations by introducing a definition of a quantum committor consistent with
any quantum dynamical map in which the propagator can be formed as a linear
superoperator acting on the system density matrix with certain states
designated as absorbing boundary conditions. This allows us to explore how
mechanisms and transition states are influenced by interference, in addition
to quantum uncertainty. We have applied this framework to investigate how the
transition state in a thermal barrier crossing reaction in a polaritonic
system is affected by coherences and to coherently control the outcome of
vertical relaxation of a system with a conical intersection.

Defining a quantum committor, or splitting probability, requires the existence
of multiple metastable states. The meaning of metastability in quantum systems
remains a topic of active investigation6; 7. In cases where a system is in
contact with a decohering environment, quantum metastability can be associated
with a spectral gap in the generator of the dynamics, evident in a separation
of relaxation timescales8; 9; 10. Imposing definitions of reactants and
products on the system, quantum incarnations of transition path sampling have
been applied using unravelled quantum master equations,4 and surface hopping
models.11; 12 Ensembles of quantum trajectories have been analyzed through
extensions of transition path theory,3 however such analysis requires the
dynamics to be secular.13; 4 This assumption decouples coherences from
populations in the energy eigenbasis, mapping the system dynamics to a
classical Markov process. Such a mapping eschews any influence of quantum
coherences by assigning each energy eigenstate a committor probability.
Interference effects significantly impact reaction rates and pathways in many
systems of chemical interest, and the inability to address coherent effects
seriously limits the applicability of these quantum committor methods.14; 15;
16; 17 This is especially true for systems with conical intersections17; 18;
19; 20; 21; 22; 23 or strongly coupled to light24; 25; 26; 27; 28 for which
the secular approximation breaks down due to gaps in the system eigenspectrum
being small compared to the rate of relaxation of the bath.

To investigate the effects of coherences on reactive behavior, we have defined
a quantum mechanical committor that can be assigned to any quantum state or
superposition. This idea is developed in the context of a partially
secularized Redfield29; 30 master equation, but is general provided an
accurate Markovian description of thermalization and decoherence. With this
generalization of the committor, we have characterized coherent quantum
effects on the transition state in a polaritonic system, finding profound
alterations in strong light-matter coupling cases. We have further employed
knowledge of the committor to engineer ideal initial states in which
interference effects guarantee relaxation into a preferred product in a
conical intersection model. The ease of selecting these ideal initial
conditions opens intriguing possibilities. Given recent advances in laser
technology which allow monitoring and manipulation of systems on timescales
relevant to electronic dynamics,31; 32; 33 the ability to determine initial
states for desired reactive outcomes has potential applications in coherent
quantum control where biasing photoisomerization results or controlling
relaxation pathways is a common goal. 34; 35; 36

## Partial Secular Master Equations

To preserve the influence of coherences while also guaranteeing complete
positivity of the density matrix, we employ a partial secularization of the
Redfield master equation. The Redfield master equation29; 37 describes the
evolution of a quantum system with total Hamiltonian $H=H_{S}+H_{B}+H_{I}$,
where $H_{S}$ is the Hamiltonian of the system, $H_{B}$ is the Hamiltonian of
the thermal bath, and $H_{I}=\sum_{k}S_{k}\otimes G_{k}$ represents the weak,
bilinear coupling between the system and bath. Operator $S_{k}$ is an operator
on the system and $G_{k}$ is the $k$’th statistically independent operator on
the bath.

The Redfield master equation describing the evolution of system density
matrix, $\rho(t)$, where the bath degrees of freedom have been traced out, can
be formulated as13

$\displaystyle\dot{\rho}(t)$ $\displaystyle=$
$\displaystyle-i\left[H_{S}+{H}_{LS},{\rho}(t)\right]+\sum_{\omega,\omega^{\prime},k}\gamma_{k}(\omega^{\prime},\omega)e^{i(\omega^{\prime}-\omega)t}$
(1)
$\displaystyle\left(S_{k}(\omega){\rho}(t)S^{*}_{k}(\omega^{\prime})-\frac{1}{2}\\{S^{*}_{k}(\omega)S_{k}(\omega),\;{\rho}(t)\\}\right)$

where $\hbar=1$ and $\dot{\rho}(t)$ indicates the time derivative of
$\rho(t)$. The operator $S_{k}(\omega)$ is the portion of the operator $S_{k}$
which produces density transfer between any eigenstates $j$ and $i$ such that
$\omega=\epsilon_{i}-\epsilon_{j}$ are the Bohr frequencies of the system
Hamiltonian, with $\epsilon_{i}$ being the $i$’th energy eigenvalue with
eigenvector $|i\rangle$. Specifically,

$S_{k}(\omega)=\sum_{\epsilon_{i}-\epsilon_{j}=\omega}|i\rangle\langle
i|S_{k}|j\rangle\langle j|,$ (2)

is the operator producing density transfer between eigenstates separated by an
energy gap of $\omega$. We define $S^{*}_{k}(\omega)$ as the conjugate
transpose of $S_{k}(\omega)$. The factors
$\gamma_{k}({\omega^{\prime},\omega})$ are given by

$\gamma_{k}({\omega^{\prime},\omega})=\Gamma_{k}({\omega})+\Gamma^{*}_{k}({\omega^{\prime}}),$
(3)

and the one sided Fourier transform of the bath correlation function define

$\Gamma_{k}({\omega})=\int_{0}^{\infty}\mathrm{Tr_{B}}[\tilde{G}^{*}_{k}(s)\tilde{G}_{k}(0)\sigma_{B}]e^{i{\omega}s}\;ds,$
(4)

with $\sigma_{B}$ being the density matrix of the thermal equilibrium state of
the bath. The interaction picture bath operators, $\tilde{G}_{k}$, are
$\tilde{G}=e^{iH_{B}t}Ge^{-iH_{B}t}$. The Lamb shift Hamiltonian,

${H}_{LS}=\sum_{k\omega\omega^{\prime}}\Pi_{k}(\omega^{\prime},\omega)S^{*}_{k}(\omega^{\prime})S_{k}(\omega),$
(5)

includes factors,

$\Pi_{k}({\omega^{\prime}},\omega)=\frac{1}{2i}\left[\Gamma_{k}({\omega})-\Gamma^{*}_{k}({\omega^{\prime}})\right],$
(6)

which are also built from one-sided Fourier transforms. This expression for
the evolution of the reduced density matrix assumes only weak system-bath
coupling and a Markovian bath.38

Due to the cross terms generated by operators $S_{k}(\omega)$ and
$S_{k}(\omega^{\prime})$ in the summation, the terms which preserve the
influence of coherences on the evolution, this equation cannot be written in
the form of a dynamical semigroup. The commonly applied secular Redfield
master equation13; 37; 39 achieves dynamical semigroup form by arguing that,
in the weak coupling limit, the unitary oscillation of the coherent phase is
much faster than population transfer. Thus the complex exponential under the
summation in Eq. 1, will oscillate many times before significant population is
transferred, averaging to zero unless $\omega^{\prime}-\omega=0$. The
nonsecular terms, any which do not fulfill this requirement, can be neglected,
decoupling populations from coherences and yielding the fully secular quantum
master equation which is in the form of a dynamical semigroup and thus
guarantees positivity of $\rho$.13 Any coherences in the initial density
matrix at time zero decay exponentially with no influence on the populations.
When a near degeneracy occurs in the system eigenspectrum, the assumption that
we may neglect the average influence of the quickly oscillating exponential
term no longer holds and there is significant influence of quantum coherences
on the evolution of populations.40; 41

An alternative partial secularization which preserves important coherences30;
41 can be derived under the piecewise flat secular approximation. This
approximation demands that the spectral density describing the bath does not
change over energy ranges encompassing all eigenstates in which coherent
effects are relevant. If the coherence between eigenstates $i$ and $j$ is
important, then for any $\omega=\epsilon_{i}-\epsilon_{l}\;l\neq i,j$ and
$\omega^{\prime}=\epsilon_{j}-\epsilon_{l}\;l\neq\;i,j$ the Fourier transforms
of the bath correlation function are equal,
$\Gamma_{k}(\omega)=\Gamma_{k}(\omega^{\prime})$. Because coherent effects are
generally large when energy differences between states are very small relative
to system timescales, this approximation is usually not demanding. To apply
the approximation, the energy eigenspectrum of $H_{S}$ is organized into
nearly degenerate nonsecular blocks by splitting the Hamiltonian into two
parts $H_{S}=H_{S}^{(0)}+\delta H_{S}$, where $H_{S}^{(0)}$ collapses all
nearly degenerate states to the same energy $\bar{\epsilon_{i}}$ and $\delta
H_{S}$ breaks the degeneracy to recover the initial Hamiltonian. If $||\delta
H_{S}||\ll||H_{I}||$, then we may perform the secular approximation on
$H_{S}^{(0)}$ that retains coherences between degenerate states (i.e. states
in the same nonsecular block) and then add in the degeneracy breaking as a
perturbative correction. The average energy of the eigenstates in a nonsecular
block is designated $\bar{\epsilon_{i}}$ for nonsecular block $i$. The energy
difference between each block is denoted
$\bar{\omega}_{ij}=\bar{\epsilon_{i}}-\bar{\epsilon_{j}}$. A block is assigned
a first eigenstate entry and if the next higher energy eigenstate is less than
a limit, $\hbar\gamma$, higher in energy, it is added to the block. The
difference in energies between the newly added eigenstate and the next higher
eigenstate is tested and eigenstates added until one exceeds $\hbar\gamma$.
This highier eigenstate is the first entry in the next block. The factor,
$\gamma$, should formally be equal to the total rate of population exit from
the eigenstate under consideration.30

The operators, $S_{k}(\bar{\omega}),$ in the partial secular quantum master
equation are formed from the system-bath interaction Hamiltonian such that
each operator is either the dephasing operator with $\bar{\omega}=0$ or
defines a jump process between nonsecular blocks $i$ and $j$, not a single
pair of eigenstates, such that

$S_{k}(\bar{\omega})=\sum_{\bar{\epsilon}_{BL(i)}-\bar{\epsilon}_{BL(j)}=\bar{\omega}}S_{k}(\omega)$
(7)

where $\bar{\epsilon}_{BL(i)}$ refers to the average energy of the nonsecular
block to which eigenstate $i$ belongs. After defining nonsecular blocks and
relevant jump operators, the partial secular master equation is identical to
the secular equation save that all $\omega$ are replaced by $\bar{\omega}$.
This reduces Eq. 1 to

$\displaystyle\dot{\rho}(t)$ $\displaystyle=$
$\displaystyle-i\left[H_{S}+H_{LS}\right]+\sum_{\bar{\omega},k}\gamma_{k}(\bar{\omega},\bar{\omega})$
$\displaystyle\left(S_{k}(\bar{\omega})\rho(t)S_{k}^{*}(\bar{\omega})-\frac{1}{2}\\{S_{k}^{*}(\bar{\omega})S_{k}(\bar{\omega}),\rho(t)\\}\right),$

which is trace preserving and completely positive. The Lamb shift Hamiltonian
becomes block diagonal, its structure corresponding with that of the selected
nonsecular blocks. Coherent evolution occurs within the nonsecular blocks,
facilitated both by the Lambshift Hamiltonian and the action of
$S_{k}(\bar{\omega})$ operators under the summation, which produce cross terms
linking coherences and populations within the nonsecular blocks only.
Coherences between eigenstates in different nonsecular blocks still decay
exponentially without influencing population dynamics, but those within the
blocks can increase or decrease under the influence of these partially
secularized operators and may have significant impacts on population dynamics.
In the case that each block contains a single eigenstate, meaning all
coherences are decoupled from all populations, the partial secular equation
reduces to the full secular approximation.

## Coherent quantum committor

The committor is conventionally defined as the probability, given sets of
system states $A$ and $B$, that the system will visit a state in $B$ before a
state in $A$, denoted $P_{B|A}$. In the context of chemical reactions set $A$
can be considered the reactants and $B$ the products, both of which are
expected to be metastable as defined by the separation of timescales evident
in the eigenspectrum of the propagator.8; 10 In the context of the partial
secular Redfield dynamics we consider, this requires that no state of $A$ or
$B$ be in a nonsecular block with any state not also in $A$ or $B$, such that
the projective measurement implied by the conditioning in the reactive
trajectory ensemble4 does not destroy coherences relevant to the reactive
dynamics.

Under these assumptions, the committor can be evaluated by constructing a
superoperator in the energy eigenbasis such that
$\dot{\rho}(t)=\mathcal{L}\rho(t)$, modified to enforce absorbing boundary
conditions in all states within the sets of states in $A$ and $B$.13; 42 For
our partial secular equation, absorbing boundary conditions are imposed by
changing the relevant rates such that
$\gamma_{k}(\bar{\omega})=0\;\forall\;\bar{\omega}=\bar{\epsilon_{j}}-\bar{\epsilon_{i}}\;\;j\in
A\cup B$. If this approach were applied to a master equation in which
coherences linked $A$ and $B$ to the rest of the system, these would have to
be approximated as zero, which would not always be a justifiable assumption.
The modified $\mathcal{L}$ with absorbing states is denoted by
$\mathcal{L}^{\prime}$. The time integrated flux into any eigenstate $b\in B$
is then

$V_{b}(\rho_{ij})=\left(\lim_{t\rightarrow\infty}e^{\mathcal{L^{\prime}}t}\rho_{ij}\right)_{b,b}$
(9)

where $\rho_{ij}=|i\rangle\langle j|$ designates an initial condition and the
subscript $b,b$ indicates that the $b$’th population of the density matrix has
been extracted after propagation. The time integrated flux for the entirety of
set of states in $B$ is given by the sum $V_{B|A}(\rho_{ij})=\sum_{b\in
B}V_{b}(\rho_{ij})$. Together, $\rho_{ij}$ for all $i$ and $j$ define a basis
for any initial density matrix. By calculating $V_{B|A}$ for all basis
entries, we can easily calculate $V_{B|A}$ for any density matrix
$\rho=\sum_{ij}c_{ij}\rho_{ij}$ by summing over the contributions of each
basis state. Because coherences which are not contained in a nonsecular block
cannot have any impact on population dynamics, $V_{b}$ only needs to be
calculated for the subset of the basis where it is not zero by definition. The
final summation,

$P_{B|A}(\rho)=\sum_{i,j}c_{ij}V_{B|A}(\rho_{ij}),$ (10)

is the probability that a system initialized in state $\rho$ will reach state
$B$ before any other state in $A$.

## Transition States under vibrational strong coupling

In our first application of the general quantum committor method, we
investigated a polariton model similar to that used by Lindoy and coworkers26.
Vibrational polaritons form when the vacuum photon frequency of a microcavity
is near resonance with the frequency of a vibrational mode of a molecule in
the cavity. Experimental investigations have indicated that resonance effects
in cavity confined polariton systems can influence ground state reactivity,
inverting the preferred formation of products in bond breaking reactions.43;
44; 45 Theoretical work has studied potential origins for the resonance
effects both by classical46; 47; 48; 49 and quantum mechanical26; 50; 51; 52
means, with evidence indicating that a fully quantum mechanical approach to
the dynamics is likely necessary to understand this phenomenon26; 28; 27; 53.

We have studied the quantum transition state in a polariton model as a
function of light-matter coupling using the generalized quantum committor.
This model includes a single photon and proton degree of freedom under a Shin-
Metiu54 model employing the Pauli-Fierz55 Hamiltonian which treats light and
matter coordinates on the same quantum footing. The Hamiltonian is given by

$H_{S}=\frac{P^{2}}{2M}+U(R)+\frac{p_{c}^{2}}{2}+\frac{\omega_{c}^{2}}{2}\left(q_{c}+\sqrt{\frac{2}{\omega_{c}}}\eta_{c}\mu(R)\right)^{2},$
(11)

where $P$ is the momentum operator of the proton coordinate with position $R$,
$M$ is the proton mass, $U(R)$ describes a quartic proton potential energy
surface, $q_{c}$ is the position operator for the photon, $p_{c}$ is the
momentum operator for the photon, $\omega_{c}$ is the photon frequency,
$\mu(R)$ is the proton dipole operator, and $\eta_{c}$ is the light-matter
coupling strength. The system-bath coupling is bilinear and independent for
both the proton and cavity modes

$H_{I}=R\otimes\sum_{k}c_{k,1}r_{k,1}+q_{c}\otimes\sum_{k}c_{k,2}r_{k,2},$
(12)

and the spectral density describing the coupling strengths, $c_{k,\alpha}$,
for each bath $\alpha=\\{1,2\\}$ is,

$J_{\alpha}(\omega)=\sum_{k}{c_{k,\alpha}^{2}}\delta(\omega-\omega_{k,\alpha})=\eta\omega
e^{-|\omega|/\omega_{b}},$ (13)

an Ohmic exponential form, where $\omega_{b}$ is the cutoff frequency and
$\eta$ is the coupling strength for both baths. Appendix A compiles the
parameters used for this study.

We treat the four lowest energy eigenstates only, with the first and second
eigenstates identified as $|g_{1}\rangle$ the reactant $A$ state and
$|g_{2}\rangle$ the product $B$ state. The third and fourth eigenstates,
$|e_{1}\rangle$ and $|e_{2}\rangle$, form a nonsecular block. The system
setup, including the four relevant eigenstates, is shown in Fig. 1 a). The
excited state set spanned by $|e_{1}\rangle$ and $e_{2}\rangle$ can be
represented by a Bloch sphere with the eigenstates at its poles. We define the
committor surface as the locations in the Bloch sphere where the probabilities
to arrive in state $|g_{2}\rangle$ or $|g_{1}\rangle$ are both $0.5$, a
generalized transition state for quantum systems. Due to the linear nature of
the propagator, the committor surface is guaranteed to be a plane. The
committor plane is easily found by determining $V_{g_{1}|g_{2}}(\rho_{ij})$
for a basis of the density matrix in the nonsecular block, solving for
conditions at which $P_{g_{1}|g_{2}}(\rho)=0.5$ and transferring the obtained
solution into Bloch sphere coordinates. To inspect the influence of quantum
coherent effects, the committor plane can be defined for the fully secular
master equation as well, in which case the plane must be defined by a single
$z$ coordinate.

Figure 1: a) The potential energy surface, $U(R)$ and lowest four eigenstates
plotted on the adabiatic potential energy surface beside an illustration of
the $|e_{i}\rangle$ states on the Bloch sphere. The partial secular (upper)
and fully secular (lower) committor planes through the Bloch sphere for b)
weak $\eta_{c}=0.01$, c) intermediate $\eta_{c}=0.1$, and d) strong
$\eta_{c}=0.2$ light-matter coupling. e) The population of state
$|g_{1}\rangle$ as a function of time given initialization in $\rho(0)$ equal
to each of the three marked density states in panel d). f) Population in
$|g_{1}\rangle$ given initialization in $\rho(0)$ equal to $\rho_{2}$ as
indicated in panel d) for master equations that are partial secular, secular,
and partial secular with the coherences manually set to zero after each time
step.

We solved for the committor planes for the partially and fully secular master
equation for three different values of light-matter coupling strength. In the
weak light-matter coupling case, Fig. 1 b), the partial secular plane is
nearly flat. However, it lies above the fully secular committor plane. This is
due to the piecewise flat spectral density approximation in the partial
secular master equation, which alters individual jump rates from those found
for equivalent jumps in the fully secular master equation. As the light-matter
coupling increases, the fully secular committor plane remains largely
unchanged whereas the partial secular committor plane acquires significant
tilt, indicating the increasing importance of quantum coherences. This is
evident for moderate coupling in Fig. 1 c) but especially obvious for the
strong coupling case in Fig. 1 d). The angle along the $x$ direction indicates
a significant influence exerted by the real part of the coherence between
$|e_{1}\rangle$ and $|e_{2}\rangle$ on determining the outcome of the
reaction. This implies that the transition state is represented by a
particular set of coherent states at the top of the barrier.

To demonstrate the influence of coherences on dynamics in the strongest
coupling case, three points on the surface of Bloch sphere, two on the fully
secular committor plane and one on the partial secular committor plane, were
selected as initial conditions to run partial secular dynamics, displayed in
Fig. 1 e). Initial conditions designated by $\rho_{1}$ and $\rho_{2}$ have
identical $z$ values, meaning identical populations in $|e_{1}\rangle$ and
$|e_{2}\rangle$. However, the initial state coherences are different and
initialization in $\rho_{2}$ results in a final population in $|g_{1}\rangle$
below $0.5$ whereas initialization in $\rho_{1}$ results in a final population
of $|g_{1}\rangle$ above $0.5$, demonstrating the strong influence of quantum
coherent effects on dynamics. By contrast the initial condition $\rho_{3}$
begins on the partial secular committor plane and thus reaches a final state
with population $0.5$, representing the equal likelihood of evolving to the
product state or back to the reactant state. The particular importance of the
coherence on the population dynamics is confirmed in Fig. 1 f), in which we
ran partial secular, fully secular, and zeroed coherence dynamics starting
from initial $\rho_{2}$. Zeroed coherence propagation indicates that the
partial secular master equation was employed but coherences in the system were
manually set to zero at each timestep, removing their influence. This process
artificially removes the population-coherence coupling while making the
piecewise flat spectral density approximation. Without effects from
coherences, the relaxation dynamics are nearly identical to the fully secular
relaxation dynamics.

## Coherent Control of a Conical Intersection

Conical intersections are commonplace in large molecules where they offer
ultrafast, nonradiative relaxation pathways,56; 57; 58 making them key in
isomerization of molecular photoswitches,59; 60 light harvesting mechanisms in
photosynthetic organisms,61; 62 and damage mitigation in DNA. 63 At a conical
intersection, the splitting between states is expected to become small,
enhancing the role of coherent and interference effects in relaxation
following photoexcitation. With the ability to probabilistically understand
reaction outcomes with the committor, we show that this knowledge can be
translated into a control strategy by preparing specific states to optimize
the relaxation outcome following photoexcitation.

The model we consider is similar to a minimal pyrazine model, with parameters
modified to produce two metastable wells in the excited state.64; 22; 17; 65;
18 The system Hamiltonian is given by,

$H_{S}=\sum_{i=1,2}|\phi_{i}\rangle
h_{i}\langle\phi_{i}|+\left(|\phi_{1}\rangle\langle\phi_{2}|+|\phi_{2}\rangle\langle\phi_{1}|\right)\lambda
Q_{c},$ (14)

with $i$ indexing the diabatic electronic states, $|\phi_{i}\rangle$, and

$h_{i}=1/2\sum_{j=c,t}\hbar\omega_{j}\left\\{P_{j}^{2}+Q_{j}^{2}\right\\}+E_{k}+\kappa_{k}Q_{t}$
(15)

where $\lambda$ is the diabatic coupling strength, $Q_{j}$ is the
dimensionless coupling ($j=c$) or tuning ($j=t$) coordinate, $P_{j}$ is the
dimensionless momentum for $Q_{j}$, $\kappa_{i}$ is the displacement along
$Q_{t}$ for each diabatic electronic state, $\omega_{j}$ is the coupling
($j=c$) or tuning frequency ($j=t$) frequency and $E_{i}$ is the diabatic
electronic state energy. The system-bath coupling Hamiltonian is bilinear in
the two modes of the conical intersection and bath modes,

$H_{I}=\left(|\phi_{1}\rangle\langle\phi_{1}|+|\phi_{2}\rangle\langle\phi_{2}|\right)\sum_{\alpha}\sum_{j=c,t}c_{\alpha,j}r_{\alpha,j}Q_{j},$
(16)

where $r_{\alpha,j}$ are dimensionless bath oscillator coordinates, and
diagonal in the diabatic electronic states. The spectral density again takes
Ohmic form,

$J_{j}(\omega)=\eta e^{-|\omega|/\omega_{b}}\;\quad j=c,t,$ (17)

where $\omega_{b}$ is the cutoff frequency and $\eta$ is the system-bath
coupling strength. The parameters of the model employed are summarized in
Appendix A.

Figure 2: a) A cut through the diabatic electronic potential energy surfaces
of the conical intersection (black) at $Q_{c}=0$ and the energies of the
fifteen eigenstates of the conical intersection, color coded by nonsecular
block, with eigenstates 13 and 14 marked with dashed and dotted lines
respectively. b) The real parts in both diabatic electronic states,
$|\phi_{1}\rangle$ and $|\phi_{2}\rangle$, of the optimal initial wavefunction
for arrival in eigenstate 1, $|\Psi_{A}\rangle$, and the optimial initial
wavefunction for arrival in eigenstate 3, $|\Psi_{B}\rangle$. c) Populations
in the lowest three eigenstates as a function of time following initialization
of the system in $\rho(0)=|\Psi_{A}\rangle\langle\Psi_{A}|$. d) Populations in
the lowest three eigenstates as a function of time following initialization of
the system such that $\rho(0)=|\Psi_{B}\rangle\langle\Psi_{B}|$.

An illustration of the structure of the system, including individual
eigenstates, non-secular blocks, and diabatic potential energy surfaces is
shown in Fig. 2 a). The lowest two eigenstates are localized in a state at
positive $Q_{t}$, and we group them into the reactant state $A$, while the
third lowest eigenstate is localized at negative $Q_{t}$ and considered the
product state, $B$. Each of these states is its own nonsecular block. There
are four other nonsecular blocks, with eigenstates 13 and 14 in block $6$ at
the top of the barrier highlighted with dashes. These states are extremely
close in energy, which suggests that quantum coherences between them are
likely to be important.

We have calculated $V_{A|B}$ for the populations and coherences of and between
eigenstates 13 and 14. We find that $V_{A|B}(\rho_{14,14})\approx 0.85$ and
that $V_{A|B}(\rho_{13,13})\approx 0.15$, with the time integrated fluxes
corresponding well with the diabatic character of these eigenstates. We also
find that $V_{A|B}(\rho_{13,14})$ is very large. To illustrate how
understanding the committment probability informs coherent quantum control, we
optimized population and coherence to direct evolution of the system into
either $A$ or $B$ by optimizing the committor value over all possible coherent
superpositions of eigenstates $13$ and $14$. When inspecting the optimized
wavefunctions displayed in Fig. 2 b), the wavefunction optimized to arrive in
$A$, $|\Psi_{A}\rangle$, with
$P_{A|B}(\rho={|\Psi_{A}\rangle\langle\Psi_{A}|})$ approaching 1.0, is nearly
entirely localized in the right well. Similarly, the wavefunction optimized to
arrive in $B$, $|\Psi_{B}\rangle$, for which the committor
$P_{B|A}(\rho={|\Psi_{B}\rangle\langle\Psi_{B}|})$ approaches 1.0, is fully
localized in the left well. In other words, this system behaves diabatically
despite being highly excited. Thus we can optimize the reactive outcomes by
using quantum interference to engineer an initial condition. This optimization
is also robust to small changes of the model as illustrated in Appendix B.

The differences in population dynamics following excitation to
$\rho(0)=|\Psi_{A}\rangle\langle\Psi_{A}|$ in Fig. 2 c) is strikingly
different from dynamics following excitation to
$\rho(0)=|\Psi_{B}\rangle\langle\Psi_{B}|$ in Fig. 2 d), which is not
unexpected given that the initial populations are very different, but reflects
the selectivity that could not be achieved by incoherent initialization into a
combination of eigenstates 13 and 14. Population in $B$, $i=3$, in Fig. 2 c)
remains at nearly zero for the entire 50000 au period whereas virtually all
population is in $B$ over that period in Fig. 2 d). Note that in Fig. 2 d)
some population has become trapped in what is effectively a dark state in
nonsecular block 4 and has not relaxed down to eigenstate 3 over the time
period displayed, but will do so in the long time limit.

The location of the density trapped in a dark state can be determined, as can
the overall relaxation pathways, by constructing a hidden Markov state model
between the nonsecular blocks, with the precise quantum state in the block
being the hidden variable which changes over time due to coherent evolution.
This hidden state of the block determines when and how density will depart
from the block. This removes information about the specific history of
individual pathways that would be evident in a quantum jump approach,3 but is
comparatively very simple to implement and provides adequate information for
our needs. In Fig. 3 a) the relaxation pathway following excitation into
$\rho(0)=|\Psi_{A}\rangle\langle\Psi_{A}|$ moves between nonsecular blocks
with 90% of density arriving in eigenstate 1 over the simulation period and 8%
remaining in eigenstate 2 with the remaining density largely found in
nonsecular block 4. This is a very different pathway than that following
excitation to $\rho(0)=|\Psi_{B}\rangle\langle\Psi_{B}|$ in Fig. 3 c) where
two alternate relaxation pathways exist, one following from 6 to 5 to 3 and
the other from 6 to 4 to 3, with almost 14% of density remaining in nonsecular
block 4 after 50000 au. Despite the difference in relaxation pathways,
inspecting probability density distributions during early times of relaxation
in 3 b) and d) shows that the system remains localized in both cases, relaxing
quickly towards the bottom of the destination well. This method of flux
tracing illuminates the general pathways of relaxation, showing that a
bifurcation is involved in optimal relaxation to $B$.

Figure 3: a) The net reactive flux traveling between nonsecular blocks for the
system initialized in $\rho(0)=|\Psi_{A}\rangle\langle\Psi_{A}|$ over a
duration of 50000 au. b) Population density in electronic state 2 for the
system relaxing from $\rho(0)=|\Psi_{A}\rangle\langle\Psi_{A}|$ at
$2000\;\mathrm{au}$ (top) and $5000\;\mathrm{au}$ (bottom). c) The net
reactive flux traveling between nonsecular blocks for the system initialized
in $\rho(0)=|\Psi_{B}\rangle\langle\Psi_{B}|$ over a duration of 50000 au.
Note that not all density arrives in block 3 over this period of time, with a
significant amount remaining in block 4. d) Population density in electronic
state 1 for the system relaxing from
$\rho(0)=|\Psi_{B}\rangle\langle\Psi_{B}|$ at $2000\;\mathrm{au}$ (top) and
$5000\;\mathrm{au}$ (bottom).

## Conclusions

By generalizing the definition of the committor to quantum coherent systems,
we have quantified the impact of coherences on dynamics in polaritonic systems
and conical intersections. Following initial calculations of committors for a
basis of the density matrix, the committor value of any quantum state can be
trivially obtained. This method can be used to define committors for any
linear master equation in which absorbing boundary conditions can be defined,
allowing fast determination of initial excitations providing the desired
outcome in coherent quantum control problems and access to the committor
surface, the quantum generalization of a transition state, which is valuable
in the study of quantum reaction mechanisms.

## Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences, CPIMS Program Early Career Research Program under Award
DEFOA0002019. D.T.L was supported by the Alfred P. Sloan Foundation.

### Appendix A: Simulation Parameters

Simulations for the polariton were carried out according to parameters in
Table 1 with a DVR66 basis for $R$ and a harmonic oscillator basis for $q_{c}$
with respective dimensions of 81 and 60 with $\Delta R$ for the DVR basis of
$0.03$ au, although a significantly smaller basis was likely sufficient. The
system was then truncated to include the lowest four energy eigenstates.

The polariton model employs a quartic potential energy surface,26

$E(R)=\frac{c_{ob}^{4}}{16c_{eb}}R^{4}-\frac{c_{ob}^{2}R^{2}}{2}-c_{cu}R^{3},$
(18)

with the dipole operator approximated by

$\mu(R)=v\mathrm{tanh}(yR)+zR,$ (19)

although its exact form does not appear to change the observed trends. The
coordinates on the surface of the Bloch sphere used to define density matrix
$\rho_{1}$ in Fig. 1 d are $(x,y,z)=(-0.87,-0.17,-0.31)$. Those for $\rho_{2}$
are $(0.87,0.383,-0.31)$ and those for $\rho_{3}$ are $(0.87,0.17,0.46)$.

Simulations for the conical intersection were carried out according to
parameters in Table 2 with a harmonic oscillator basis for $Q_{c}$ and $Q_{t}$
with dimensions 40 and 90 respectively, although a significantly smaller basis
was likely sufficient. The system was then truncated to include fifteen
eigenstates.

Parameter |  atomic units unless specified
---|---
$\beta$ |  1052.584413
$\omega_{b}$ |  0.05
$\omega_{c}$ |  0.025344
$\eta$ |  0.2 (unitless)
$c_{ob}$ |  0.8
$c_{eb}$ |  0.05
$M$ |  1836
$c_{cu}$ |  0.004
$v$ |  -1.7
$y$ |  3.0 (unitless)
$z$ |  0.6

Table 1: Parameters employed during simulation of the polariton model

.

Parameter |  atomic units unless specified
---|---
$\beta$ |  1052.584413
$\omega_{b}$ |  0.01
$\eta$ |  0.1 (unitless)
$E_{1}$ |  -0.00139
$E_{2}$ |  0.00139
$\omega_{c}$ |  0.004116
$\omega_{t}$ |  0.002279
$\kappa_{1}$ |  -0.006836
$\kappa_{2}$ |  0.006836
$\lambda$ |  0.00091153
Table 2: Parameters employed during simulation of the conical intersection
model. Note that this simulation is carried out with $Q_{c,t}$ and $P_{c,t}$
in dimensionless units.

### Appendix B: Coherent Control Stability Analysis

Figure 4: The real part of the optimal wavefunctions, $\Psi_{m}$ (multiplied
by a factor for easy viewing) for arrival in eigenstates $1$ and $2$ in the
unperturbed and two perturbed systems, where the perturbation adjusts
$\omega_{t}$ relative to the default value of $0.002279\;\mathrm{au}$. The
change occurs smoothly. The optimal wavefunction is largely localized in the
lower energy well in all cases.

In order to assess the sensitivity of the optimized wavefunctions and
optimization efficiency to small changes in the Hamiltonian parameters, we
determined the optimal initial wavefunction for relaxation into $A$ obtainable
from a coherent superposition of eigenstates 13 and 14 as we modified
$\omega_{t}$. The original $\omega_{t}$ is designated $\omega_{0}$.
Modifications were small to avoid changing the energy ordering of the
eigenstates. Optimization of relaxation outcomes was effective for
significantly perturbed systems. In Fig. 4, the optimal initial wavefunctions
for the case of slightly smaller $\omega_{t}$, original $\omega_{t}$ and
slightly larger $\omega_{t}$ are seen to be very similar, although as
$\omega_{t}$ increases the amount of wavefunction density located in the
metastable well, diabatic electronic state 1, increases, and, correspondingly,
efficiency decreases somewhat.

Comparing the overlap between the optimal wavefunction for the unperturbed
system with $\omega_{t}=\omega_{0}$ and the optimal wavefunction for systems
with perturbed $\omega_{t}$ in Fig. 5, the overlap decreases in tandem with
the fall in the committor value obtained from employing the optimal
wavefunction for $\omega_{t}=\omega_{0}$ as an initial condition in the
perturbed system. However, the decrease in overlap and committor value is not
drastic indicating that small variations in the Hamiltonian will not disrupt
the process of wavefunction optimization. Similarly, in Fig. 5 an optimal
wavefunction with yield of 98% or higher could be identified for all perturbed
systems, indicating that the chosen parameters in this work are not unique in
providing opportunities for optimization.

Perturbations of other parameters produced similar results or, in the case of
$\omega_{c}$ and $E_{2}$, showed much greater stability of the optimized
wavefunction to parameter perturbation. A sufficiently large perturbation of
any parameter will, however, completely change the eigenstate ordering and,
though an optimal coherent initialization may still exist, it will no longer
involve a superposition of eigenstates 13 and 14 and direct comparisons cannot
be drawn.

Figure 5: Asterisk symbols indicating the overlap between the optimal
wavefunction in the unperturbed case and the optimal wavefunction in the
perturbed case, $|\Psi_{0}^{*}\Psi_{m}|$, together with committor to $A$,
designated $P_{A|B}(\rho)$, as $\omega_{t}$ is modified for both the optimal
initialization density matrix for the unperturbed system, $\rho_{0}$, and the
perturbed optimal initialization density matrix, $\rho_{m}$.

### Appendix C: Derivation of Quantum Committor

In this section, we will derive Eq. (10) for the committor from the main text.
This derivation will proceed by first defining the quantum jump unraveling
corresponding to the partial secular master equation, Eq. (Partial Secular
Master Equations), and using this stochastic Poisson process to define a
commitor equivalently to standard classical Poisson processes. Then, we define
an auxiliary stochastic process that implements absorbing boundary conditions
in the reactant and product non-secular blocks while leaving dynamics
unchanged outside of this manifold, which conveniently expresses the commitor.
Finally, we average over the realizations of the jump process to obtain a new
quantum master equation with absorbing boundary conditions yielding Eq. (9).
Crucially, we see that this committor, and the resulting interpretation, are
independent of the unraveling used in the derivation.

Equation (Partial Secular Master Equations) is a quantum master equation that
we can equivalently solve as the average over a stochastic Poisson process for
a density matrix, $\hat{\rho}$, called the quantum jump unraveling,

$\displaystyle d\hat{\rho}(t)$ $\displaystyle=$
$\displaystyle\left[\mathcal{L}_{0}+\sum_{\Omega}\gamma_{\Omega}\left(\mathcal{A}_{\Omega}+\langle
n_{\Omega}(\hat{\rho})\rangle\right)\right]\hat{\rho}dt$ (20a)
$\displaystyle+\sum_{\Omega}\left(\frac{\mathcal{J}_{\Omega}}{\langle
n_{\Omega}(\hat{\rho})\rangle}-\mathds{1}\right)\hat{\rho}dN_{\Omega}(\hat{\rho})$
$\mathbb{E}\left[dN_{\Omega}(\hat{\rho})\right]=\gamma_{\Omega}\langle
n_{\Omega}(\hat{\rho})\rangle dt$ (20b) $\langle
n_{\Omega}(\hat{\rho})\rangle\equiv\operatorname{Tr}\\{\mathcal{J}_{\Omega}\hat{\rho}\\}=-\operatorname{Tr}\\{\mathcal{A}_{\Omega}\hat{\rho}\\},$
(20c)

where the double index $\Omega=(k,\overline{\omega})$ has been introduced for
brevity, $\gamma_{\Omega}\equiv\gamma_{k}(\bar{\omega},\bar{\omega})$,
$\mathds{1}$ is the identity superoperator,
$\mathcal{L}_{0}=-i[H_{S}+H_{LS},\cdot]$ is the unitary system propagator,
$\mathcal{J}_{\Omega}=S_{k}(\overline{\omega})\cdot
S_{k}^{*}(\overline{\omega})$ is the jump operator describing the change in
system state due to observing a change in bath eigenstate, and
$\mathcal{A}_{\Omega}=\\{{S_{k}^{*}(\overline{\omega})S_{k}(\overline{\omega}),\cdot\\}}$
is the anti-unitary (sometimes called Zeno) drift due to not observing a
change in bath eigenstate. The differentials $dN_{\Omega}$ are the increments
of independent Poisson processes with state dependent rates given by Eq.
(20b). Each realization of the stochastic process $\hat{\rho}$ corresponds to
a different history of bath measurement outcomes. In the partial secular limit
we know that the jumps couple superpositions in one secular manifold to
superpositions in another, decoupled from intra-manifold coherences. We will
use the decorator $\hat{\cdot}$ throughout to indicate random variables
defined for single realizations of the unraveling. This stochastic
differential equation is connected to Eq. (Partial Secular Master Equations)
by $\mathbb{E}[\hat{\rho}]=\rho$, where $\mathbb{E}$ denotes an average over
realizations of the stochastic unraveling. Conditional density matrices remain
normalized but the evolution is nonlinear in $\sigma$ since $\langle
n_{\Omega}(\hat{\rho})\rangle$ is itself a linear function of $\hat{\rho}$.

We now proceed to define the committor equivalently to the case of a classical
jump process. We are able to use this approach since the quantum jump
unraveling is simply a classical jump process over an unusual state-space of
density matrices. First, we select a reactant non-secular block $A$ and a
product block $B$. The committor can be constructed from the trajectory
observable $\hat{P}_{A|B}[\hat{\rho}(t)]$ which is $1$ if the trajectory
visits $B$ before $A$ and $0$ otherwise. If we average over all trajectories
with initial state $\bm{\rho}$, we obtain the committor
$P_{A|B}(\bm{\rho})=\mathbb{E}\left[\hat{P}_{A|B}[\hat{\rho}(t)]|\hat{\rho}(0)=\bm{\rho}\right]$.
The commutator can then be computed by sampling realizations of this
stochastic process starting in state $\bm{\rho}$ and computing the average
over the trajectory ensemble. However, since the stochastic equation is
nonlinear in $\hat{\rho}$ this must be recomputed for every initial state
$\bm{\rho}$ repeatedly which can rapidly become costly since the size of the
state space grows rapidly with increasing system dimensions.

It is therefore desirable to obtain an expression for the committor using a
quantum master equation which is linear in $\bm{\rho}$. We could then compute
the contribution to the committor from each density matrix element
$V_{B|A}(\rho_{ij})$ separately. The linearity of the master equation would
then allow us to compute the commutator for any
$\bm{\rho}=\sum_{ij}c_{ij}\rho_{ij}$, a superposition of the density matrix
elements $\rho_{ij}$, by evaluating the same superposition of
$V_{B|A}(\rho_{ij})$. To accomplish this, we define an auxiliary unraveling
$\hat{\tilde{\rho}}$, identical to Eq. (20), but removing all jumps out of
blocks $A$ and $B$. This stochastic process has the same dynamics as the
original unraveling outside of $A\cup B$ but does not allow population to
leave $A$ or $B$ once it enters either block. If we ignore all dynamics within
$A$ and $B$ and consider only the total population of the blocks, this
procedure implements two absorbing boundary conditions at the reactant and
product states.

Using the absorbing auxiliary unraveling, $\hat{\tilde{\rho}}$, we can easily
rewrite the trajectory observable,

$\hat{P}_{A|B}[\hat{\rho}(t)]=\hat{P}_{A|B}[\hat{\tilde{\rho}}(t)]=\lim_{t\to\infty}\operatorname{Tr}\\{\mathcal{P}_{B}\hat{\tilde{\rho}}(t)\\}$
(21)

as the steady-state population of non-secular product block $B$, written in
terms of the $B$ projection operator $\mathcal{P}_{B}=\sum_{b\in
B}\ket{b}\bra{b}$. For each trajectory, this quantity is always exactly $1$ or
$0$ since the probability of simultaneous jumps vanishes, and in the non-
secular master equation each jump can only connect one pair of blocks.
Therefore, it is impossible for the unraveling to simultaneously enter blocks
$A$ and $B$.

The introduction of the absorbing process has now reduced the calculation of
the committor to the expectation of a steady-state population. This can be
equivalently done after averaging over the stochastic processes to obtain an
auxiliary quantum master equation with generator $\mathcal{L}^{\prime}$,
yielding Eq. (10) in the main text. Notably, this definition of the committor
is independent of the choice of unraveling used in the derivation, indicating
that it captures the properties of the underlying dynamical process
independent of any assumed measurement protocol on the bath.

## References

## References

  * Onsager (1938) L. Onsager, “Initial recombination of ions,” Phys. Rev. 54, 554 (1938).
  * Hummer (2004) G. Hummer, “From transition paths to transition states and rate coefficients,” J. Chem. Phys. 120, 516–523 (2004).
  * Anderson, Schile, and Limmer (2022) M. C. Anderson, A. J. Schile, and D. T. Limmer, “Nonadiabatic transition paths from quantum jump trajectories,” J. Chem. Phys. 157, 164105 (2022).
  * Schile and Limmer (2018) A. J. Schile and D. T. Limmer, “Studying rare nonadiabatic dynamics with transition path sampling quantum jump trajectories,” J. Chem. Phys. 149, 214109 (2018).
  * Manzano, Horowitz, and Parrondo (2015) G. Manzano, J. M. Horowitz, and J. M. Parrondo, “Nonequilibrium potential and fluctuation theorems for quantum maps,” Phys. Rev. E 92, 032129 (2015).
  * Patra and Keshavamurthy (2015) S. Patra and S. Keshavamurthy, “Classical-quantum correspondence in a model for conformational dynamics: Connecting phase space reactive islands with rare events sampling,” Chem. Phys. Lett. 634, 1–10 (2015).
  * De Leon and Marston (1989) N. De Leon and C. C. Marston, “Order in chaos and the dynamics and kinetics of unimolecular conformational isomerization,” J. Chem. Phys. 91, 3405–3425 (1989).
  * Macieszczak  _et al._ (2016) K. Macieszczak, M. Guţă, I. Lesanovsky, and J. P. Garrahan, “Towards a theory of metastability in open quantum dynamics,” Phys. Rev. Lett. 116, 240404 (2016).
  * Macieszczak  _et al._ (2021) K. Macieszczak, D. C. Rose, I. Lesanovsky, and J. P. Garrahan, “Theory of classical metastability in open quantum systems,” Phys. Rev. Res. 3, 033047 (2021).
  * Brown, Macieszczak, and Jack (2024) C. A. Brown, K. Macieszczak, and R. L. Jack, “Unraveling metastable Markovian open quantum systems,” Phys. Rev. A 109, 022244 (2024).
  * Reiner  _et al._ (2023) M. M. Reiner, B. Bachmair, M. X. Tiefenbacher, S. Mai, L. González, P. Marquetand, and C. Dellago, “Nonadiabatic forward flux sampling for excited-state rare events,” J. Chem. Theory Comput. 19, 1657–1671 (2023).
  * Coffman  _et al._ (2023) A. J. Coffman, Z. Jin, J. Chen, J. E. Subotnik, and D. V. Cofer-Shabica, “Use of QM/MM surface hopping simulations to understand thermally activated rare-event nonadiabatic transitions in the condensed phase,” J. of Chem. Theory Comput. 19, 7136–7150 (2023).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press on Demand, Oxford, 2002).
  * Keefer  _et al._ (2020) D. Keefer, T. Schnappinger, R. de Vivie-Riedle, and S. Mukamel, “Visualizing conical intersection passages via vibronic coherence maps generated by stimulated ultrafast X-ray Raman signals,” Proc. Natl. Acad. Sci. U.S.A. 117, 24069–24075 (2020).
  * Schlosshauer (2019) M. Schlosshauer, “Quantum decoherence,” Phys. Rep. 831, 1–57 (2019).
  * Yang, Xiang, and Xiong (2020) Z. Yang, B. Xiang, and W. Xiong, “Controlling quantum pathways in molecular vibrational polaritons,” ACS Photonics 7, 919–924 (2020).
  * Duan, Miller, and Thorwart (2016) H.-G. Duan, R. J. D. Miller, and M. Thorwart, “Impact of vibrational coherence on the quantum yield at a conical intersection,” J. Phys. Chem. Lett. 7, 3491–3496 (2016).
  * Qi  _et al._ (2017) D.-L. Qi, H.-G. Duan, Z.-R. Sun, R. J. D. Miller, and M. Thorwart, “Tracking an electronic wave packet in the vicinity of a conical intersection,” J. Chem. Phys. 147 (2017).
  * Wang, Allodi, and Engel (2019) L. Wang, M. A. Allodi, and G. S. Engel, “Quantum coherences reveal excited-state dynamics in biophysical systems,” Nat. Rev. Chem. 3, 477–490 (2019).
  * Zhu and Yarkony (2016) X. Zhu and D. R. Yarkony, “Non-adiabaticity: the importance of conical intersections,” Mol. Phys. 114, 1983–2013 (2016).
  * Worth and Cederbaum (2004) G. A. Worth and L. S. Cederbaum, “Beyond Born-Oppenheimer: Molecular dynamics through a conical intersection,” Annu. Rev. Phys. Chem. 55, 127–158 (2004).
  * Schile and Limmer (2019) A. J. Schile and D. T. Limmer, “Simulating conical intersection dynamics in the condensed phase with hybrid quantum master equations,” J. Chem. Phys. 151, 014106 (2019).
  * Halpern and Limmer (2020) N. Y. Halpern and D. T. Limmer, “Fundamental limitations on photoisomerization from thermodynamic resource theories,” Phys. Rev. A 101, 042116 (2020).
  * Mandal and Huo (2019) A. Mandal and P. Huo, “Investigating new reactivities enabled by polariton photochemistry,” J. Phys. Chem. Lett. 10, 5519–5529 (2019).
  * Bhuyan  _et al._ (2023) R. Bhuyan, J. Mony, O. Kotov, G. W. Castellanos, J. Gómez Rivas, T. O. Shegai, and K. Börjesson, “The rise and current status of polaritonic photochemistry and photophysics,” Chem. Rev. 123, 10877–10919 (2023).
  * Lindoy, Mandal, and Reichman (2023) L. P. Lindoy, A. Mandal, and D. R. Reichman, “Quantum dynamical effects of vibrational strong coupling in chemical reactivity,” Nat. Commun. 14, 2733 (2023).
  * Fiechter  _et al._ (2023) M. R. Fiechter, J. E. Runeson, J. E. Lawrence, and J. O. Richardson, “How quantum is the resonance behavior in vibrational polariton chemistry?” J. of Phys. Chem. Lett. 14, 8261–8267 (2023).
  * Li, Nitzan, and Subotnik (2020) T. E. Li, A. Nitzan, and J. E. Subotnik, “On the origin of ground-state vacuum-field catalysis: Equilibrium consideration,” J. Chem. Phys. 152, 234107 (2020).
  * Redfield (1965) A. G. Redfield, “The theory of relaxation processes,” in  _Adv. Magn. Reson._, Advances in Magnetic and Optical Resonance, edited by J. S. Waugh (Academic Press, 1965) pp. 1–32.
  * Trushechkin (2021) A. Trushechkin, “Unified Gorini-Kossakowski-Lindblad-Sudarshan quantum master equation beyond the secular approximation,” Phys. Rev. A 103, 062226 (2021).
  * Zhu  _et al._ (1997) L. Zhu, K. Suto, J. A. Fiss, R. Wada, T. Seideman, and R. J. Gordon, “Effect of resonances on the coherent control of the photoionization and photodissociation of HI and DI,” Phys. Rev. Lett. 79, 4108–4111 (1997).
  * Krausz and Ivanov (2009) F. Krausz and M. Ivanov, “Attosecond physics,” Rev. Mod. Phys. 81, 163–234 (2009).
  * Arasaki  _et al._ (2011) Y. Arasaki, K. Wang, V. McKoy, and K. Takatsuka, “Monitoring the effect of a control pulse on a conical intersection by time-resolved photoelectron spectroscopy,” Phys. Chem. Chem. Phys. 13, 8681–8689 (2011).
  * Lucas and Hornberger (2014) F. Lucas and K. Hornberger, “Incoherent control of the retinal isomerization in rhodopsin,” Phys. Rev. Lett. 113, 058301 (2014).
  * Vogt  _et al._ (2005) G. Vogt, G. Krampert, P. Niklaus, P. Nuernberger, and G. Gerber, “Optimal control of photoisomerization,” Phys. Rev. Lett. 94, 068305 (2005).
  * Herek  _et al._ (2002) J. L. Herek, W. Wohlleben, R. J. Cogdell, D. Zeidler, and M. Motzkus, “Quantum control of energy flow in light harvesting,” Nature 417, 533–535 (2002).
  * Nitzan (2006) A. Nitzan,  _Chemical dynamics in condensed phases: Relaxation, transfer and reactions in condensed molecular systems_ (Oxford University Press, Oxford, 2006).
  * Levy, Rabani, and Limmer (2021) A. Levy, E. Rabani, and D. T. Limmer, “Response theory for nonequilibrium steady states of open quantum systems,” Phys. Rev. Res. 3, 023252 (2021).
  * Balzer and Stock (2005) B. Balzer and G. Stock, “Modeling of decoherence and dissipation in nonadiabatic photoreactions by an effective-scaling nonsecular Redfield algorithm,” Chem. Phys. 310, 33–41 (2005).
  * Egorova, Kühl, and Domcke (2001) D. Egorova, A. Kühl, and W. Domcke, “Modeling of ultrafast electron-transfer dynamics: multi-level Redfield theory and validity of approximations,” Chem. Phys. 268, 105–120 (2001).
  * Cattaneo  _et al._ (2019) M. Cattaneo, G. L. Giorgi, S. Maniscalco, and R. Zambrini, “Local versus global master equation with common and separate baths: superiority of the global approach in partial secular approximation,” New J. Phys. 21, 113045 (2019).
  * Sharpe and Wales (2021) D. J. Sharpe and D. J. Wales, “Nearly reducible finite markov chains: Theory and algorithms,” J. Chem. Phys. 155 (2021).
  * Thomas  _et al._ (2019) A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen, “Tilting a ground-state reactivity landscape by vibrational strong coupling,” Science 363, 615–619 (2019).
  * Thomas  _et al._ (2016) A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux, C. Genet, J. A. Hutchison, and T. W. Ebbesen, “Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field,” Angew. Chem. 128, 11634–11638 (2016).
  * Lather  _et al._ (2019) J. Lather, P. Bhatt, A. Thomas, T. W. Ebbesen, and J. George, “Cavity catalysis by cooperative vibrational strong coupling of reactant and solvent molecules,” Angew. Chem. Int. Edit. 58, 10635–10638 (2019).
  * Li, Mandal, and Huo (2021a) X. Li, A. Mandal, and P. Huo, “Theory of mode-selective chemistry through polaritonic vibrational strong coupling,” J. Phys. Chem. Lett. 12, 6974–6982 (2021a).
  * Philbin  _et al._ (2022) J. P. Philbin, Y. Wang, P. Narang, and W. Dou, “Chemical reactions in imperfect cavities: Enhancement, suppression, and resonance,” J. Phys. Chem. C 126, 14908–14913 (2022).
  * Galego  _et al._ (2019) J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, “Cavity Casimir-Polder forces and their effects in ground-state chemical reactivity,” Phys. Rev. X 9, 021057 (2019).
  * Li, Mandal, and Huo (2021b) X. Li, A. Mandal, and P. Huo, “Cavity frequency-dependent theory for vibrational polariton chemistry,” Nat. Comm. 12, 1315–1324 (2021b).
  * Schäfer  _et al._ (2022) C. Schäfer, J. Flick, E. Ronca, P. Narang, and A. Rubio, “Shining light on the microscopic resonant mechanism responsible for cavity-mediated chemical reactivity,” Nat. Comm. 13, 7817 (2022).
  * del Pino, Feist, and Garcia-Vidal (2015) J. del Pino, J. Feist, and F. J. Garcia-Vidal, “Quantum theory of collective strong coupling of molecular vibrations with a microcavity mode,” New J. Phys. 17, 053040–53050 (2015).
  * Yang and Cao (2021) P.-Y. Yang and J. Cao, “Quantum effects in chemical reactions under polaritonic vibrational strong coupling,” J. Phys. Chem. Lett. 12, 9531–9538 (2021).
  * Anderson  _et al._ (2023) M. C. Anderson, E. J. Woods, T. P. Fay, D. J. Wales, and D. T. Limmer, “On the mechanism of polaritonic rate suppression from quantum transition paths,” J. Phys. Chem. Lett. 14, 6888–6894 (2023).
  * Shin and Metiu (1995) S. Shin and H. Metiu, “Nonadiabatic effects on the charge transfer rate constant: A numerical study of a simple model system,” J. Chem. Phys. 102, 9285–9295 (1995).
  * Pauli and Fierz (1938) W. Pauli and M. Fierz, “Zur theorie der emission langwelliger lichtquanten,” Il Nuovo Cimento (1924-1942) 15, 167–188 (1938).
  * Neville  _et al._ (2016) S. P. Neville, Y. Wang, A. E. Boguslavskiy, A. Stolow, and M. S. Schuurman, “Substituent effects on dynamics at conical intersections: Allene and methyl allenes,” J. Chem. Phys. 144, 014305 (2016).
  * Matsika (2005) S. Matsika, “Three-state conical intersections in nucleic acid bases,” J. Phys. Chem. A 109, 7538–7545 (2005).
  * Curchod and Martínez (2018) B. F. E. Curchod and T. J. Martínez, “Ab initio nonadiabatic quantum molecular dynamics,” Chem. Rev. 118, 3305–3336 (2018).
  * Muždalo  _et al._ (2018) A. Muždalo, P. Saalfrank, J. Vreede, and M. Santer, “Cis-to-Trans isomerization of azobenzene derivatives studied with transition path sampling and quantum mechanical/molecular mechanical molecular dynamics,” J. Chem. Theory Comput. 14, 2042–2051 (2018).
  * Bull  _et al._ (2018) J. N. Bull, M. S. Scholz, E. Carrascosa, G. da Silva, and E. J. Bieske, “Double molecular photoswitch driven by light and collisions,” Phys. Rev. Lett. 120, 223002 (2018).
  * Song  _et al._ (2021) H. Song, Y. Nam, D. Keefer, M. Garavelli, S. Mukamel, and S. Tretiak, “Nonadiabatic molecular dynamics study of the relaxation pathways of photoexcited cyclooctatetraene,” J. Phys. Chem. Lett. 12, 5716–5722 (2021).
  * Savolainen  _et al._ (2008) J. Savolainen, R. Fanciulli, N. Dijkhuizen, A. L. Moore, J. Hauer, T. Buckup, M. Motzkus, and J. L. Herek, “Controlling the efficiency of an artificial light-harvesting complex,” Proc. Natl. Acad. Sci. U.S.A. 105, 7641–7646 (2008).
  * Barbatti  _et al._ (2010) M. Barbatti, A. J. A. Aquino, J. J. Szymczak, D. Nachtigallová, P. Hobza, and H. Lischka, “Relaxation mechanisms of UV-photoexcited DNA and RNA nucleobases,” Proc. Natl. Acad. Sci. U.S.A. 107, 21453–21458 (2010).
  * Chen  _et al._ (2016) L. Chen, M. F. Gelin, V. Y. Chernyak, W. Domcke, and Y. Zhao, “Dissipative dynamics at conical intersections: simulations with the hierarchy equations of motion method,” Faraday Discus. 194, 61–80 (2016).
  * Duan  _et al._ (2018) H.-G. Duan, D.-L. Qi, Z.-R. Sun, R. D. Miller, and M. Thorwart, “Signature of the geometric phase in the wave packet dynamics on hypersurfaces,” Chem. Phys. 515, 21–27 (2018).
  * Colbert and Miller (1992) D. T. Colbert and W. H. Miller, “A novel discrete variable representation for quantum mechanical reactive scattering via the s‐matrix kohn method,” J. Chem. Phys. 96, 1982–1991 (1992).