# Non-adiabatic Dynamics in a Continuous Circularly Polarized Laser Field with
Floquet Phase-space Surface Hopping

Zeyu Zhou  Yanze Wu  Xuezhi Bian  Joseph Eli Subotnik<EMAIL_ADDRESS>[

###### Abstract

Non-adiabatic chemical reactions involving continuous circularly polarized
light (cw CPL) have not attracted as much attention as dynamics in
unpolarized/linearly polarized light. However, including circularly (in
contrast to linearly) polarized light allows one to effectively introduce a
complex-valued time-dependent Hamiltonian, which offers a new path for control
or exploration through the introduction of Berry forces. Here, we investigate
several inexpensive semiclassical approaches for modeling such nonadiabatic
dynamics in the presence of a time-dependent complex-valued Hamiltonian,
beginning with a straightforward instantaneous adiabatic fewest-switches
surface hopping (IA-FSSH) approach (where the electronic states depend on
position and time), continuing to a standard Floquet fewest switches surface
hopping (F-FSSH) approach (where the electronic states depend on position and
frequency), and ending with an exotic Floquet phase-space surface hopping
(F-PSSH) approach (where the electronic states depend on position, frequency,
and momentum). Using a set of model systems with time-dependent complex-valued
Hamiltonians, we show that the Floquet phase-space adiabats are the optimal
choice of basis as far as accounting for Berry phase effects and delivering
accuracy. Thus, the F-PSSH algorithm sets the stage for modeling nonadiabatic
dynamics under strong externally pumped circular polarization in the future.

penn]Department of Chemistry, University of Pennsylvania, Philadelphia,
Pennsylvania 19104, United States

##  1 Introduction

Non-adiabatic transitions between electronic states typically arise in two
different contexts. First, transitions occur naturally through vibronic
interactions when molecules visit regions of configuration space where the
Born-Oppenheimer approximation is violated; second, transitions can be induced
by photo-excitations when an external incident light is coupled to the
transition dipole moment between these electronic states. Both processes are
very important in the field of photochemistry and spectroscopy,1, 2, 3, 4, 5
and both processes need not occur exclusively (i.e. both can occur at the same
time). One important difference between vibronic couplings and light-induced
couplings is that the latter is time-dependent; a typical light source
contains a central frequency $\omega$ such that the coupling contains
$\cos\omega t$. Experiments have shown that strong monochromatic continuous
wave (cw) light can change the landscape of potential energy surfaces and
reaction channels by introducing light-induced states, or Floquet states in
several different ways; e.g., there is now experimental evidence of light-
induced conical intersections.6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19

During the past thirty years, various semiclassical formalisms for studying
non-adiabatic phenomena have been demonstrated as effective and reasonably
accurate. More recently, many surface hopping formalisms 20 have been
generalized to incorporate time-dependent radiative couplings for light-
induced non-adiabatic processes.21, 22, 23, 24, 25, 26, 27, 28, 1, 29 One of
the most intuitive formalisms is to instantaneously calculate the adiabatic
potential energy surfaces, whereby one instantaneously diagonalizes the light-
matter Hamiltonian. As a result, time-derivative coupling matrix elements
$\mathbf{T}_{jk}=\left\langle\psi_{j}(\mathbf{R}(t),t)\right|\frac{d}{dt}\left|\psi_{k}(\mathbf{R}(t),t)\right\rangle$
contains contributions that arise from both nuclear motion and from the
explicit time-dependence of the external field
$\Bigl{\\{}\left|\psi_{k}(\mathbf{R}(t),t)\right\rangle\Bigr{\\}}$. In some
cases, these resulting dynamics can perform well24, 26, 1, 29, but the
algorithm faces a difficult choice when simulating energy absorption/emission
from the external light; usually, when deciding whether a hop is accepted or
frustrated, one compares the bandwidth and the energy differences during the
nonadiabatic transitions. Energy conservation as a function of photon number
is difficult to implement and so the algorithm can lose accuracy.

Apart from IA-FSSH, another possible generalization to the time-dependent
nonadiabatic problem is to apply Floquet theory30, 31, 27, 28, 29, which
transforms a time-periodic Hamiltonian into a time-independent one with larger
dimensions, leading to Floquet fewest switches surface hopping algorithm
(F-FSSH). In our own experience, we have found that, in the presence of
monochromatic light (and provided the frequency of the light is not very
small), F-FSSH performs better in most cases than IA-FSSH insofar as the
algorithm better captures energy absorption/emission, which is manifested as
transitions between Floquet states with different Fourier indices;
interestingly, when the frequency of light is small, IA-FSSH performs better,
as the Hamiltonian approaches a time-independent form. For the most part, one
can use these two algorithms to reasonably capture most standard nonadiabatic
dynamics in a linearly polarized light field.

Now, if we wish to study dynamics under a circular polarized light field (CPL)
with large frequency, and if we adopt an electric dipole Hamiltonian, it is
fairly easy to conclude that, in a Floquet representation, the light-matter
coupling terms become inherently complex-valued (and this complex-valued
nature cannot be eliminated by any simple gauge transformation). For instance,
for a two-level electronic Hamiltonian subject to an external circularly
polarized laser (e.g. $\mathbf{E}(t)=E_{x}\cos\omega t\hat{x}+E_{y}\sin\omega
t\hat{y}$), in the bare electronic basis and under electric dipole
approximation, the light-matter coupling becomes

$\left\langle\mu\right|H(t)\left|\nu\right\rangle=\mathbf{D}\cdot\mathbf{E}(t)=\mu_{x}E_{x}\cos\omega
t+\mu_{y}E_{y}\sin\omega t$ (1)

Here, $\mathbf{D}$ is the transition dipole moment between electronic states
$\left|\mu\right\rangle$ and $\left|\nu\right\rangle$. Note that this coupling
term is completely real and time-dependent. Next, if we apply a Fourier
transform, the resulting light-matter couplings in Floquet basis will be of
the form

$\left\langle
m\mu\right|H(t)\left|(m\pm\tilde{1})\nu\right\rangle=\mu_{x}E_{x}/2\pm
i\mu_{y}E_{y}/2$ (2)

Here, $m$ is the Fourier index. Clearly, these light-matter couplings between
Floquet states with $\pm\tilde{1}$ Fourier indices difference are complex-
valued. Note that the transition dipole moment will usually depend on the
nuclear configuration and hence, the phases of these coupling terms are not
identical for different nuclear configurations.

Unfortunately, the introduction of a complex-valued Hamiltonian renders most
of surface hopping schemes inapplicable. On the one hand, it is widely
acknowledged that during a hopping event between two multidimension potential
energy surfaces, to enforce energy conservation, the momentum of the
trajectory is rescaled along the direction of derivative coupling
$\mathbf{d}$.32 Thus, it is not obvious what direction to choose when the
derivative coupling $\mathbf{d}$ is complex-valued. 33, 34 On the other hand,
complex-valued derivative coupling $\mathbf{d}$ can alter the nuclear motion
by Berry phase effects. As is well-known, Berry curvature emerges when the
derivative couplings become complex-valued, resulting in an effective magnetic
field for the nuclear degrees of freedom.

Over the last few years, several attempts have been made to extend the
standard fewest switches surface hopping (FSSH) approach so as to treat a
complex-valued Hamiltonian by finding a good rescaling direction and
incorporating Berry phase effects33, 34, 35, 36, 37, 38, 39. Within our group,
we have developed two such algorithms:

  1. 1.

FSSH with ad hoc Berry forces and "$\mathbf{h}$ \+ $\mathbf{k}$" rescaling
direction in ref 38 (FSSH h+k). The basic premise of FSSH h+k is to take into
account the complex nature of the derivative couplings $\mathbf{d}$, which
leads to a local effective magnetic field known as the Berry curvature. For
this algorithm, the rescaling direction after a hop depends lies in the plane
spanned by the derivative of the norm of the diabatic Hamiltonian elements
($\mathbf{h}$) and the the derivative of the phase $\mathbf{k}$ (see section
2.2.3 for definitions).

  2. 2.

Phase-space surface hopping algorithm in ref 39 (PSSH). The basic premise is
to transform the complex-valued Hamiltonian into a real-valued one locally by
introducing a phase factor that induces a momentum shift during an electronic
transition.

To date, neither of these algorithms (or any other FSSH algorithm, as far as
we are aware) has been successfully extended so as to model the nonadiabatic
dynamics of light under a CPL laser field. The goal of this article is to
create and benchmark such extensions.

Finally, before concluding this section, we note that CPL is not the only
means by which one can introduce complex-valued Hamiltonians (and Berry
forces) into nonadiabatic light-matter systems. More generally, Berry forces
arise when there is degeneracy of the Hamiltonian, and if one models dynamics
without external light – but with spin degrees of freedom and spin-orbit
couplings – one will also find that a complex-valued Hamiltonian arises.
Moreover, recently there has been speculation that chiral induced spin
selectivity (CISS) effects might arise precisely through such coupled nuclear-
spin motion. 40, 41, 42 Thus, it is important to emphasize that all of the
theory presented below for modeling the dynamics of nuclei and electrons in a
circularly polarized light field can be equally applied to modeling the
dynamics of nuclei and electrons and spin in a strong linearly polarized light
field.

With this background in mind, an outline of this article is as follows: In
section 2, we first review in detail the IA-FSSH and the Floquet-FSSH
algorithms which were designed for a system periodically driven by a linearly
polarized light. Secondly, we discuss how existing extensions of FSSH to treat
complex-valued Hamiltonian can potentially be incorporated into F-FSSH. In
section 3, we describe the model that we will use to differentiate these
different FSSH approaches, and we will offer some visual intuition. In section
4, the results of the formalisms discussed above are compared with exact
quantum calculations. We conclude and discuss several intriguing questions
regarding F-PSSH in section 5.

##  2  Methods

###  2.1  Model Hamiltonian and exact solution

Let us consider a molecule illuminated by continuous wave circularly polarized
light (cw-CPL). As discussed in section 1, the incoming light, in the electric
dipole approximation, becomes a real-valued time-periodic coupling
$\hat{V}(t)=\hat{V}(t+T_{0})$ between electronic states with period $T_{0}$ in
the original electronic basis (see eq 1). As shown in eq 2, in the Floquet
basis, the coupling becomes time-independent and complex-valued.

As discussed above, if spin degrees of freedom are present, the vibronic
coupling can also be complex-valued – due to spin-orbit coupling or under the
influence of an external static magnetic field. Hence, for the most general
two level system, let us consider a Hamiltonian of the form:

$\displaystyle\hat{H}_{\text{tot}}(t)=\hat{\mathbb{T}}_{\mathbf{R}}+\hat{H}_{el}^{0}+\hat{V}(t)=\hat{\mathbb{T}}_{\mathbf{R}}+\hat{H}_{el}(t)$
(3)

Here, $\hat{\mathbb{T}}_{\mathbf{R}}$ is the kinetic operator and
$\hat{H}_{el}^{0}$ is the time-independent part of the complex-valued
electronic Hamiltonian $\hat{H}_{el}$ which we write as

$\displaystyle\hat{H}_{el}(t)$ $\displaystyle=\left[\begin{array}[]{c
c}H^{el}_{00}(\mathbf{R})&H^{el}_{01}(\mathbf{R},t)\\\
H^{el}_{10}(\mathbf{R},t)&H^{el}_{11}(\mathbf{R})\\\ \end{array}\right]$ (6)

Here, the couplings between the two electronic states contain two
contribution: a time-independent vibronic coupling term and a time-periodic
light-induced coupling term:

$\displaystyle H^{el}_{10}(\mathbf{R},t)$
$\displaystyle=D_{a}(\mathbf{R})\exp(i\phi_{a}(\mathbf{R}))+D_{b}(\mathbf{R})\exp(i\phi_{b}(\mathbf{R}))\cos\omega
t$ (7)

Here, $D_{a}$ and $D_{b}$ are the effective time-independent and time-
dependent coupling strengths, and $\phi_{a}(\mathbf{R})$ and
$\phi_{b}(\mathbf{R})$ are the phases of the couplings.

The exact solution of such a time-dependent nonadiabatic problem can be
obtained by propagating the Schrödinger equation on a grid using short time
steps and exponentiating $e^{-i\hat{H}_{tot}(t)dt}$ at each time step, where
$\hat{H}_{tot}(t)$ is the time-dependent total Hamiltonian. The goal of this
article is to assess cheap, semiclassical approaches to such propagation.

###  2.2 Four semiclassical methods for time-dependent coupled nuclear-
electronic dynamics

####  2.2.1 Instantaneous Adiabatic Fewest Switches Surface Hopping

Let us begin by reviewing the simplest extension to original FSSH for model
problems with time-dependent Hamiltonian. As proposed by Gonzalez and
Marquetand,24, 26, 43, 1 the basic premise is to use the instantaneous
adiabatic potential energy surfaces (that are explicitly time-dependent) to
replace the adiabatic potential energy surfaces (that are only parameterized
by nuclear configuration $\mathbf{R}$). The nuclear degrees of freedom are
evolved by Newton’s equations of motion

$\displaystyle\dot{\mathbf{R}}$ $\displaystyle=\frac{\mathbf{P}}{M}$ (8)
$\displaystyle\dot{\mathbf{P}}$
$\displaystyle=-\frac{\nabla_{\mathbf{R}}E_{\lambda}(\mathbf{R},t)}{M}.$ (9)

Here, $E_{\lambda}$ is the instantaneous, active adiabatic state of the
electronic Hamiltonian eq 6. The electronic degrees of freedom are evolved by
the electronic time-dependent Schrödinger equations.

$\displaystyle\dot{c}_{j}=-\frac{i}{\hbar}E_{j}(\mathbf{R},t)c_{j}-i\hbar\sum_{k}\mathbf{P}\cdot\mathbf{d}_{jk}c_{k}/M$
(10)

Here, $\mathbf{d}_{jk}$ is the derivative coupling matrix element between
instantaneous adiabatic electronic states $\left|\psi_{j}\right\rangle$ and
$\left|\psi_{k}\right\rangle$. In practice, we evaluate the time-derivative
matrix instead

$\mathbf{T}_{jk}=\mathbf{P}\cdot\mathbf{d}_{jk}/M=\left\langle\psi_{j}(\mathbf{R},t)\middle|\frac{d\psi_{k}(\mathbf{R},t)}{dt}\right\rangle,$
(11)

which contains contributions from both the time-dependent part
$\left\langle\psi_{j}(\mathbf{R},t)\middle|\frac{\partial\psi_{k}(\mathbf{R},t)}{\partial
t}\right\rangle$ and the nuclear motion part
$\left\langle\psi_{j}(\mathbf{R},t)\middle|\nabla_{\mathbf{R}}\psi_{k}(\mathbf{R},t)\right\rangle\cdot\frac{\mathbf{P}}{M}$.

Similar to FSSH, the hopping probability from active instantaneous adiabatic
state $\lambda$ to state $j$ is

$\displaystyle g_{\lambda
j}=\frac{-2\text{Re}(c_{\lambda}c_{j}^{*}T_{k\lambda})dt}{|c_{\lambda}|^{2}}$
(12)

When a hop from state $\lambda$ to state $j$ occurs, the momentum is only
rescaled along the direction of $\mathbf{d}_{\lambda j}$, if the energy
difference between state $\lambda$ and state $j$ is not within the bandwidth
of the time-dependent driving.

####  2.2.2 Floquet Theory and Floquet Fewest Switches Surface Hopping

In this subsection, we review Floquet theory and Floquet fewest switches
surface hopping (F-FSSH).44, 30, 25, 27, 29 For any problem with real-valued
periodic Hamiltonian $\hat{H}(t)=\hat{H}(t+T_{0})$, the time-dependent
electronic Schrödinger equation is

$\displaystyle i\hbar\frac{\partial}{\partial
t}\left|\Psi(t)\right\rangle=\hat{H}(t)\left|\Psi(t)\right\rangle$ (13)

We define the Floquet Hamiltonian as

$\displaystyle\hat{H}_{F}(t)\equiv\hat{H}(t)-i\hbar\frac{\partial}{\partial
t}$ (14)

and Floquet diabatic basis as

$\displaystyle\left|m\mu\right\rangle=\exp(im\omega t)\left|\mu\right\rangle$
(15)

Here, $\left|\mu\right\rangle=\\{\left|0\right\rangle,\left|1\right\rangle\\}$
belongs to a set of orthonormal properly chosen diabatic electronic basis. $m$
is the Fourier basis index, which represents the number of photons dressed by
the Floquet state $\left|m\mu\right\rangle$. Note that here and below, we will
use a superscript tilde to differentiate the “Floquet photon” index
($m=\tilde{0},\pm\tilde{1},...$) from the diabatic state index
($\left|\mu\right\rangle=\left|0\right\rangle,\left|1\right\rangle$).

In the basis $\left\\{\left|m\mu\right\rangle\right\\}$, the elements of the
Floquet Hamiltonian $\hat{H}_{F}(t)$ become time-independent:

$\displaystyle[\hat{\cal H}_{F}]_{(n\nu)(m\mu)}$
$\displaystyle=\frac{1}{T_{0}}\int_{0}^{T_{0}}dt\left\langle\nu\right|\hat{H}_{F}\left|m\mu\right\rangle\exp[-in\omega
t]$
$\displaystyle=\frac{1}{T_{0}}\int_{0}^{T_{0}}dt\left\langle\nu\right|\hat{H}(t)\left|\mu\right\rangle\exp[-i(n-m)\omega
t]+\delta_{\mu\nu}\delta_{mn}n\hbar\omega$ (16)

Hence, another possible way to propagate exact dynamics is to project the
initial wavefunction onto the Floquet states with $m=\tilde{0}$ and evolve the
system with the propagator for the time-independent Floquet Hamiltonian. For
the complex-valued model Hamiltonians in this paper, the explicit matrix form
is written in Appendix A.

Now, let us briefly review F-FSSH algorithm.27, 29 The nuclear degrees of
freedom $\mathbf{R},\mathbf{P}$ are evolved by Newton’s equations of motion
(see eq 8 and 9) along the active Floquet adiabatic potential energy surface
$E^{F}_{\lambda}(\mathbf{R})$ (eigenvalues of the Floquet Hamiltonian in eq
16, see also Appendix A) of the trajectory. Note that the Floquet adiabatic
potential energy surfaces do not explicitly depend on time (which is different
from IA-FSSH). Similar to standard FSSH, the electronic degrees of freedom
follow eq 10 except that (i) the propagation follows the Floquet adiabatic
quasi-energy $E_{\lambda}^{F}(\mathbf{R})$ (rather than instantaneous active
adiabatic energy $E_{\lambda}(\mathbf{R},t)$) and (ii) the relevant time-
derivative matrix element
$\mathbf{T}_{jk}^{F}=\left\langle\psi_{j}^{F}\middle|\frac{d\psi_{k}^{F}}{dt}\right\rangle$
is between Floquet adiabatic states (rather than $\mathbf{T}_{jk}$ in eq 11).

Lastly, the hopping probability from active Floquet state $\lambda$ to state
$k$ is given by the analogue of eq 12. That being said, there is the question
of how to rescale momenta after a hopping event because $\mathbf{d}_{\lambda
j}$ is complex-valued. To that end, consider a general two-level system of the
form:

$\displaystyle\hat{H}_{el}$ $\displaystyle=V\left[\begin{array}[]{c
c}-\cos\theta&e^{i\phi}\sin\theta\\\
e^{-i\phi}\sin\theta&\cos\theta\end{array}\right]$ (19)

For this Hamiltonian, the derivative couplings lie in the vector space spanned
by the two directions $\mathbf{h}$ and $\mathbf{k}$

$\displaystyle\mathbf{h}$ $\displaystyle=\nabla_{\mathbf{R}}\theta$ (20)
$\displaystyle\mathbf{k}$
$\displaystyle=\nabla_{\mathbf{R}}\phi-\frac{\nabla_{\mathbf{R}}\theta(\nabla_{\mathbf{R}}\theta\cdot\nabla_{\mathbf{R}}\phi)}{|\nabla_{\mathbf{R}}\theta|^{2}}$
(21)

We will follow the convention in ref 38, 39 and rescale momenta in the
direction $\mathbf{h}=\nabla_{R}\theta$.

Finally, at the end of the calculation, there is always the question of how to
calculate electronic observables, e.g. the population on a given electronic
state. For standard FSSH, there is no unique means of calculating such an
observable, but from both a theoretical and practical perspective, a density
matrix approach usually performs best.45 Below, we will avoid such nuances and
calculated electronic populations only in the asymptotic limit where the
diabats and adiabats are equal, and where there is no light-matter coupling.
In such a case, the final electronic population can be evaluated (to a good
approximation) by summing up the populations of all Floquet states that
correspond to a given electronic state:46

$\displaystyle\text{Prob}_{\mu}^{\text{pop}}=\frac{\sum_{m}N_{m\mu}^{\text{traj}}}{N^{\text{traj}}_{tot}}$
(22)

####  2.2.3 F-FSSH algorithm with Berry force (or “FSSH with h+k rescaling”)

The third algorithm that we have tested aims to improve the F-FSSH algorithm
by explicitly including Berry forces along dynamics on one surface 38 and
making sure that momentum rescaling yield the correct asymptotic values. To
motivate such an approach and explain how the algorithm works in practice, let
us treat these two effects separately:

Berry Forces: Recall that Berry forces are effective magnetic fields that
emerge with complex-valued Hamiltonians. The explicit form of the Berry force
is:

$\displaystyle
F_{\lambda}^{B}=\frac{2\hbar}{M}\sum_{j}\text{Im}(\mathbf{P}_{\lambda}\cdot\mathbf{d}_{j\lambda})\mathbf{d}_{\lambda
j}.$ (23)

As discussed in ref 38, one can calculate the Berry force (eq 23) and apply it
along a given FSSH trajectory.

Momentum Rescaling: For a complex-valued Hamiltonian, as compared to a real-
valued one, another relevant difference is that the derivative coupling matrix
$\mathbf{d}_{jk}$ becomes complex-valued. Hence, the rescaling direction is no
longer well-defined. To that end, we will follow the ansatz in ref 38. First,
even though we have many Floquet states, we calculate the rescaling direction
only during a possible hopping event between the active Floquet state and the
target Floquet state during a trajectory. Second, as discussed in ref 38, when
the trajectory has sufficient kinetic energy to hop, we add a component of the
momentum along the $\mathbf{k}$ (eq 21) direction, and then we rescale the
momentum (for energy conservation) along $\mathbf{h}$ (eq 20):

$\displaystyle\mathbf{P}_{k}=\mathbf{P}_{\lambda}+(|\left\langle\zeta_{0}\middle|k\right\rangle|^{2}-|\left\langle\zeta_{0}\middle|\lambda\right\rangle|^{2})\hbar\mathbf{k}+\alpha\mathbf{h}$
(24)

Here, $\zeta_{0}$ is the first diabatic state in eq 19, $\left|k\right\rangle$
is the adiabat that the trajectory is hopping to, $\alpha$ is the prefactor
determined by energy conservation. If there is no real solution to $\alpha$,
that is, the energy is not sufficient to supplement the momentum change, we
use a test momentum $\mathbf{P}^{\text{test}}$ to see if the kinetic energy is
consumed because of the Berry force (eq 23)

$\displaystyle\mathbf{P}^{\text{test}}=\mathbf{P}_{\lambda}-(|\left\langle\zeta_{0}\middle|\lambda\right\rangle|^{2}-|\left\langle\zeta_{0}\middle|\zeta\right\rangle|^{2})\hbar\mathbf{k}+\beta\mathbf{h}$
(25)

Here, $\zeta$ is the initial diabat. If there exist a real solution to
$\beta$, we set

$\displaystyle\mathbf{P}_{k}=\mathbf{P}_{\lambda}-\frac{\mathbf{P}\cdot\mathbf{h}}{|\mathbf{h}|^{2}}\mathbf{h}+\gamma\hbar\mathbf{k}.$
(26)

Again, $\gamma$ is determined by energy conservation.

The algorithm above is complicated and was developed empirically out of
necessity in order to match a set of data. The intuition is that one ensures
the correct asymptotic momentum along $\mathbf{k}$ by depleting the momentum
along $\mathbf{h}$. That being said, if the energy does not allow for a real-
valued $\beta$, the hop is frustrated/ (See section 3 for information about
velocity reversal.) Finally, at the end of the calculation, we evaluate the
final electronic population in the same fashion as in eq 22.

####  2.2.4 Floquet Phase-Space Surface Hopping Algorithm

In this subsection, we review our fourth candidate algorithm, a combination of
Phase-Space surface hopping (PSSH) together with a Floquet FSSH formalism. 39

Let us begin by reviewing the PSSH algorithm for two electronic states and
consider the problem in the form of eq 19. We will perform a local gauge
transformation where we introduce a phase on each electronic state (but
without rotating the states explicitly):47

$\displaystyle\left[\begin{array}[]{c}\left|\zeta_{0}^{\prime}\right\rangle\\\
\left|\zeta_{1}^{\prime}\right\rangle\end{array}\right]$
$\displaystyle=U\left[\begin{array}[]{c}\left|\zeta_{0}\right\rangle\\\
\left|\zeta_{1}\right\rangle\end{array}\right]=\left[\begin{array}[]{c
c}1&0\\\
0&\exp(-i\phi)\end{array}\right]\left[\begin{array}[]{c}\left|\zeta_{0}\right\rangle\\\
\left|\zeta_{1}\right\rangle\end{array}\right]$ (35)

The electronic Hamiltonian (eq 19) then becomes real-valued:

$\displaystyle\hat{H}_{el}$ $\displaystyle=V\left[\begin{array}[]{c
c}-\cos\theta&\sin\theta\\\ \sin\theta&\cos\theta\end{array}\right]$ (38)

and the kinetic energy part of the total Hamiltonian becomes

$\displaystyle\hat{\mathbb{T}}_{\mathbf{R}}=\frac{(\hat{\mathbf{P}}-i\hbar\hat{\mathbf{D}})^{2}}{2M}$
(39)

Here,
$\hat{\mathbf{D}}=\frac{U^{\dagger}\hat{P}U-\hat{P}}{i\hbar}=-i\nabla_{\mathbf{R}}\phi$$\left|\zeta_{1}^{\prime}\right\rangle$$\left\langle\zeta_{1}^{\prime}\right|$
where $\zeta_{1}^{\prime}$ is the basis after the gauge transform in eq 35.
The basic idea of PSSH is that one diagonalizes the Hamiltonian
$H(R,P)=\hat{H}_{el}+\mathbb{T}_{\bf R}$ that depends on both position and
momentum, and one then moves along the corresponding eigenstates.

Now, let us consider the problem with time-dependent couplings (in eq 7) and
discuss how to construct the corresponding F-PSSH equations of motion. Once
the Hamiltonian (in eq 6) is turned into a time-independent one as in eq 16,
we construct a basis set with dimensionality $N_{elec}N_{Fourier}$, where
$N_{elec}=2$ corresponds to the two electronic states indices and
$N_{Fourier}$ is the truncated number of Fourier indices at which the
calculation converges (in our calculation $N_{Fourier}=9$ as
$m=\tilde{0},\pm\tilde{1},\dots,\pm\tilde{4}$). We seek a local gauge
transformation similar to eq 35 and a new basis:

$\displaystyle\Bigl{\\{}\dots;\left|\xi^{{}^{\prime}(-\tilde{1})0}\right\rangle;\left|\xi^{{}^{\prime}(-\tilde{1})1}\right\rangle;\left|\xi^{{}^{\prime}\tilde{0}0}\right\rangle;\left|\xi^{{}^{\prime}\tilde{0}1}\right\rangle;\left|\xi^{{}^{\prime}\tilde{1}0}\right\rangle;\left|\xi^{{}^{\prime}\tilde{1}1}\right\rangle;\dots\Bigl{\\}}$
(40)

where the Floquet Hamiltonian will be strictly real-valued. Here, the ′ labels
the states after the local gauge transformation.

Unfortunately, for this model Hamiltonian (and presumably most Hamiltonians),
there is no local gauge transformation under which the Floquet Hamiltonian is
strictly real-valued. One possible approximation then is to focus on the most
important couplings (i.e. those the couplings mixing populated states that are
relevant during a surface hopping calculation) and render those couplings
real-valued (or as real-valued as possible). Framed mathematically, we seek
diagonal matrices $\hat{\cal U}_{F}$ and $\hat{\cal D}_{F}$ with the same
dimension ($N_{elec}N_{Fourier}$) as the time-independent Floquet Hamiltonian
that make the original Floquet Hamiltonian as close to real-valued as
possible.

To construct such matrices explicitly, let us return to the original Floquet
diabatic basis, where the Floquet Hamiltonian is of the form (same as eq 16):

$\displaystyle[\hat{\cal H}_{F}]_{(n\nu)(m\mu)}$
$\displaystyle=\frac{1}{T_{0}}\int_{0}^{T_{0}}dt\left\langle\nu\right|\hat{H}_{el}(t)\left|\mu\right\rangle\exp[-i(n-m)\omega
t]+\delta_{\mu\nu}\delta_{mn}n\hbar\omega$ (41)

with some complex-valued matrix elements arising potentially from both
vibronic couplings and light-matter couplings. Eq 16 is written out in matrix
form in Appendix A. Let us suppose (without loss of generality) that the
initial state corresponds to $\left|\tilde{0}0\right\rangle$. To make the
Hamiltonian as real-valued as possible, we will conjugate this Hamiltonian by
a diagonal matrix:

$\displaystyle\text{diag}(\hat{\cal U}_{F})_{\text{initial
state}=\left|\tilde{0}0\right\rangle}$
$\displaystyle=[\dots,e^{i\phi_{a}+\phi_{b}},e^{i\phi_{b}},1,e^{-i\phi_{a}},e^{-i\phi_{a}-i\phi_{b}},e^{-i\phi_{b}},\dots]$
(42) $\displaystyle\text{diag}(\hat{\cal D}_{F})_{\text{initial
state}=\left|\tilde{0}0\right\rangle}$
$\displaystyle=[\dots;-\mathbf{w}_{a}-\mathbf{w}_{b};-\mathbf{w}_{b};0;\mathbf{w}_{a};\mathbf{w}_{a}+\mathbf{w}_{b};\mathbf{w}_{b};\dots]$
(43)

Here, $\mathbf{w}_{a}=-i\nabla_{\mathbf{R}}\phi_{a}$ and
$\mathbf{w}_{b}=-i\nabla_{\mathbf{R}}\phi_{b}$ arise from the phase factors in
the vibronic coupling in eq 7. Note that, if the initial state were different
(e.g. $\left|\tilde{0}1\right\rangle$), a similar local gauge transformation
($(\hat{\cal U}_{F})_{\text{initial state}=\left|\tilde{0}1\right\rangle}$)
could also be defined. For more details about and an explicit representation
of eqs 41 \- 43, see Appendix B.

At this point, let us assume the Floquet phase-space Hamiltonian is close to
real-valued but depends on both nuclear coordinates $\mathbf{R}$ and momenta
$\mathbf{P}$,

$\displaystyle[\hat{\cal H}_{F}]^{{}^{\prime}}_{(n\nu)(m\mu)}$
$\displaystyle=\frac{1}{T_{0}}\int_{0}^{T_{0}}dt\left\langle\nu\right|\hat{H}_{el}(t)\left|\mu\right\rangle\exp[-i(n-m)\omega
t]+\delta_{\mu\nu}\delta_{mn}n\hbar\omega-\frac{\hbar^{2}\hat{\cal
D}_{F}^{2}}{2m}-\frac{i\hbar\hat{P}\cdot\hat{\cal D}_{F}}{m}$ (44)

F-PSSH then follows the same procedures as one would expect when combining
F-FSSH and PSSH. In a basis of boosted (momentum-dependent) Floquet diabatic
states (from eq 44), one diagonalizes the momentum-dependent Hamiltonian and
obtains momentum-dependent Floquet adiabatic states, or Floquet phase-space
adiabats. All subsequent dynamics move along these phase-space adiabats. In
section 3 below, we will present both the model problem and visualize the
difference between Floquet phase-space adiabats and the original Floquet
adiabats.

##  3 Simulation details

In this paper, we will focus on nonadiabtic models with two nuclear
dimensions, one light-induced avoided crossing and one vibronic coupling
induced avoided crossing. All parameters are in arbitrary units. One last
small note is now in order. In all of the four algorithms discussed above, if
a frustrated hop is encountered, the trajectory is reversed along the
rescaling direction $\mathbf{h}$ if
$(\mathbf{P}\cdot\mathbf{h})(\nabla_{\mathbf{R}}E^{F}_{j}\cdot\mathbf{h})>0$.
32, 48

###  3.1 Model problems

We choose the following diabatic states

$\displaystyle H^{el}_{00}(\mathbf{R},t)$ $\displaystyle=A\tanh(Bx)$ (45)
$\displaystyle H^{el}_{11}(\mathbf{R},t)$ $\displaystyle=-A(\tanh(Bx)+C)$ (46)

Here, $A=0.1$, $B=0.35$, $C=-0.9$. as they cross near $x=1.5$ and the energy
difference between them matches external driving $\hbar\omega=0.18$ near
$x=-1.5$. For better convergence, the time-dependent couplings (eq 7) need to
be localized near each type of crossings. Thus, we choose $D_{a}(\mathbf{R})$
and $D_{b}(\mathbf{R})$ to be gaussian functions centered around the two
avoided crossings respectively

$\displaystyle D_{a}(\mathbf{R})=E\exp(-(x-1.5)^{2}/2)$ (47) $\displaystyle
D_{b}(\mathbf{R})=E\exp(-(x+1.5)^{2}/2).$ (48)

In this paper, $E=0.02$. We pick simple phase factors as follows:

$\displaystyle\phi_{a}$ $\displaystyle=W_{a}y,\ W_{a}=6$ (49)
$\displaystyle\phi_{b}$ $\displaystyle=W_{b}y,\ W_{b}=\pm 5$ (50)

These phase factors yield two reasonably strong regions of effective magnetic
fields. Note our notation: while lowercase $\mathbf{w}_{a/b}$ (defined below
eq 43) represents the gradient of the phase $\phi_{a/b}$, here we study a
problem where both $\mathbf{w}_{a}$ and $\mathbf{w}_{b}$ depend on only one
nuclear coordinate, $y$, and so $\phi_{a/b}$ has a gradient along only one
direction $y$; above and below, we have denoted this number as $W_{a/b}$.

The relevant potential energy surfaces are presented in Figure 1. The standard
F-FSSH algorithm and F-FSSH with Berry force algorithms evolve trajectories
along the standard Floquet states. These Floquet states are parallel and
equally spaced. In comparison, for the F-PSSH algorithm, the trajectories run
along Floquet phase-space adiabatic states. These states are also parallel,
but are not equally spaced, as shown by the dashed lines.

Figure 1:  Floquet adiabatic potential quasi-energy surfaces and Floquet
phase-space adiabatic potential quasi-energy surfaces for our model. The
original F-FSSH and F-FSSH with Berry force algorithms run trajectories along
the Floquet adiabatic quasi-energy surfaces; here, the surfaces belonging to
the same electronic state are parallel and have equal quasi-energy spacing.
The Floquet phase-space adiabatic potential quasi-energy surfaces are parallel
but are not equally-spaced (because of the $\hat{\cal D}_{F}^{2}$ terms in eq
44, resulting in $W_{a}^{2}$ and $W_{b}^{2}$ terms, see also Appendix B). Note
that the states drawn on the right are only for normal incidence
($P_{y}(t=0)=0$). For oblique incidence, the resulting potential quasi-energy
surfaces will also depend on $P_{y}(t=0)$.

###  3.2 Initial Conditions

For the exact, F-FSSH, F-FSSH with Berry force and F-PSSH calculations, we
choose $x(t=0)=-6.0$ which is far enough from $x=-1.5$ such that the initial
diabats and adiabats have a one-to-one correspondence. For all surface hopping
algorithms, the initial coordinates and momenta are sampled from the Wigner
distribution of the two-dimensional Gaussian wavepacket.

For our two-dimensional model, the effective magnetic fields are parallel to
$y$ axis. Thus, we will discuss two possible incident angles, normal incidence
($y(t=0)=0$, $P_{y}(t=0)=0$) and oblique incidence at $y(t=0)=x_{0}$ and
$P_{y}(t=0)=P_{x}(t=0)=P_{init}$ with respect to the $y$ axis in Section 4.

###  3.3 Common techniques

Before presenting our results, here we will list a few useful tricks and
techniques that we used so as to simulate our calculations reliably and
efficiently.

  * •

Separation of classical and quantum time steps as in ref 49. We used a nuclear
time step of 0.5 au and an electronic time step is chosen dynamically.

  * •

Matrix logarithm to calculate the time-derivative matrix instead of derivative
couplings.50, 51, 52, 53, 54, 55, 56, 57, 49, 58, 59 This approach allows for
a larger time step without sacrificing any accuracy at all.

  * •

Parallel transport by maximal phase alignments as in ref 60. Surface hopping
within Floquet theory is especially sensitive to the phases of the adiabatic
states, and these phases must be chosen accurately.

##  4 RESULTS

In this section, we first present the transmission and reflection
probabilities obtained by exact calculations, F-FSSH, IA-FSSH and F-PSSH.
Second, we compare results between F-PSSH and F-FSSH with Berry force.

###  4.1 F-PSSH vs standard F-FSSH and IA-FSSH algorithms

####  4.1.1 Normal incidence

In Fig 2, we present the transmission and reflection probabilities on
electronic states $0$ and $1$ for different $W_{b}$ and initial electronic
states. Several observations can be made.

First, in all of the subfigures, the F-PSSH results (blue line with crosses)
agree perfectly with the exact results (black solid line). Second, standard
F-FSSH (green line, labeled as F-FSSH h rescaling) cannot capture the correct
transmission probabilities for the case $W_{b}=-5$ but does give a reasonably
good answer for the case $W_{b}=+5$. Interestingly (and incorrectly), F-FSSH
predicts almost the same results for $W_{b}=-5$ vs $W_{b}=+5$. In truth,
however, exact scattering results are different for these cases because the
presence of a crossing around $x=1.5$ with complex-valued vibronic couplings
($W_{a}=+6$) breaks any symmetry in $W$. (If we were to set $W_{a}=0$, then we
would recover the same results for $\pm W_{b}$). Third, the results obtained
by IA-FSSH (magenta line) algorithm are almost the same as standard F-FSSH. In
the low momentum regime with $W_{b}=+5$, IA-FSSH is less accurate than F-FSSH.
As with F-FSSH, IA-FSSH incorrectly predicts identical results for $W_{b}=-5$
vs $W_{b}=+5$, which suggests that the algorithm will have similar
difficulties with complex-valued Hamiltonians.

Figure 2:  Transmission and reflection probabilities at normal incidence
($P_{y}(t=0)=0$) for different choices of initial diabatic electronic states
(init $0$ or init $1$) and $W_{b}=\pm 5$. As shown in all subfigures, the
F-PSSH algorithm (blue line) yields the most accurate results as compared
against the exact results (black line). Note that for the cases with opposite
time-dependent phase factors ($W_{b}=-5$ vs $W_{b}=+5$), the standard F-FSSH
scheme (green dotted line) predicts the same results (as does IA-FSSH).
However, as shown by the exact results, these opposite $W_{b}$ cases are not
identical (because of the crossing at $x=1.5$) and admit different scattering
probabilities. Thus, this model problem highlights why one requires a
nonadiabatic algorithm that can correctly treat complex-valued Hamiltonians.

####  4.1.2 Oblique incidence

Next, in Fig 3, we present the transmission and reflection probabilities for
the case of oblique incidence ($P_{y}(t=0)=P_{x}(t=0)=P_{init}$). This regime
is a more difficult test of a surface hopping algorithm. Similar to Figure 2,
F-PSSH does accurately recover the exact results in almost all scenarios. In
contrast, standard F-FSSH and IA-FSSH cannot give satisfying results except in
the high momentum regime ($P_{init}>20\ a.u.$), where effectively $P_{y}(t)$
is not sensitive to the relatively small oscillations in phase
($W_{a/b}<<P_{y}$) that each trajectory experiences when moving along the
initial electronic state. Note that standard F-FSSH and IA-FSSH algorithms
still yield approximately the same results for opposite $W_{b}$ ($W_{b}=-5$ vs
$W_{b}=+5$), as they lack the ability to capture the effects of complex-valued
couplings. Overall, F-PSSH algorithm outperforms IA-FSSH and standard F-FSSH.

Figure 3:  Transmission and reflection probabilities for the case of oblique
incidence ($P_{y}(t=0)=P_{x}(t=0)$) for different light-induced phase factors
($W_{b}=\pm 5$) and different initial electronic states (init $0$ or init
$1$). The F-PSSH results overlap with the exact results in all figures. As in
Figure 2, standard F-FSSH (with $h$-rescaling) and IA-FSSH cannot capture the
differences between $W_{b}=-5$ vs $W_{b}=+5$, and the results are thus
inaccurate. Moreover, these two algorithms yield incorrect results in the low
momentum regime $P_{init}<20\ a.u.$. Overall, F-PSSH algorithm is clearly the
most reliable algorithm.

###  4.2 F-PSSH vs F-FSSH with Berry force

Figure 4:  Transmission and reflection probabilities for $W=-5$ with normal
($\mathbf{P}(t=0)=(P_{init},0)$) and oblique incidence
($\mathbf{P}(t=0)=(P_{init},P_{init})$). In this figure, we compare F-FSSH
with $\mathbf{h}+\mathbf{k}$ rescaling vs exact and F-PSSH results. Both of
these algorithms include Berry force in some fashion and both are expected to
perform well. In all figures here, we find that the F-FSSH with
$\mathbf{h}+\mathbf{k}$ rescaling algorithm can capture the correct
probabilities at high momentum ($P_{init}>20$) and is reasonably good at low
momentum as well; nevertheless, the algorithm is outperformed by F-PSSH in the
low momentum regime where predicting the final $\mathbf{P}$ is crucial for
capturing scattering probabilities. Nevertheless, if we were to plot F-FSSH
$\mathbf{h}+\mathbf{k}$ rescaling data within Fig 2 and Fig 3, we would find
that the algorithm predicts results just as accurately as does F-PSSH.

So far, we have demonstrated that F-PSSH results agree well with the exact
results whereas the F-FSSH and IA-FSSH algorithms do not–presumably because
the latter two algorithms do not include an Berry force effects. To better
probe the value of a phase-space surface hopping code, in this subsection, we
will present a detailed comparison of the F-PSSH algorithm with the F-FSSH
with Berry force algorithm – in other words, we compare two semiclassical
algorithms that both do account for Berry force.

In most scenarios treated within Fig 2 and Fig 3 (and in most other scenarios
we tested), F-PSSH and F-FSSH with h+k rescaling (and Berry force) yield
equally accurate results. However, for the two specific scenarios shown in Fig
4, the F-FSSH with Berry force algorithm (red dashed line, labeled as F-FSSH
h+k rescaling) is able to give accurate predictions only for high initial
momentum ($P_{init}>20\ a.u.$)38. In the low initial momentum regime, however,
the F-FSSH with Berry force algorithm allows for an erroneous $P_{y}$ during
the simulation run, and this error eventually leads to major deviation from
the exact results. While we are convinced that the F-PSSH algorithm is likely
the most accurate approach in general, F-FSSH with Berry force may be a
satisfactory ansatz in many cases.

##  5 Conclusions and Discussions

In summary, we have considered four different surface hopping algorithms:

  * •

instantanously adiabatic (IA)-FSSH

  * •

Floquet-FSSH

  * •

Floquet-FSSH with Ad Hoc Berry Forces

  * •

Floquet-PSSH

We have benchmarked these algorithm as far as treating a time-dependent and
complex-valued Hamiltonian, as relevant to coupled nuclear-electronic problems
in a circularly polarized light field).

For our model problems, both IA-FSSH and the standard F-FSSH algorithms
completely fail to capture the resulting Berry phase effects acting on the
nuclear motion and transmission/reflection probabilities which arise from the
complex-valued nature of the Hamiltonian. Between the remaining two algorithms
which do include Berry force (at least to some extent), the algorithm with an
ad hoc Berry force and rescaling direction ($\mathbf{h}+\mathbf{k}$) can
predict accurate results in the high momentum regime; the Floquet phase-space
surface hopping algorithm can recover the exact results almost always both
qualitatively and quantitatively.

Interestingly, as a side note, the Hamiltonian in eqs 6, 7 and 50 is spatially
periodic in $Y$ (with period $2\pi$, or angular frequency $\Omega_{Y}=1$).
Thus, just as in Bloch theory, it would appear natural to construct an
electronic basis (with a phase parameterized by the nuclear position in the
$Y$-direction) of the form

$\displaystyle\left|\Xi^{pm\mu}\right\rangle=\exp(ip\Omega_{Y}Y)\exp(im\omega_{T}t)\left|\mu\right\rangle$
(51)

where $\omega_{T}$ is the temporal frequency of external periodic driving and
$p$ is a Fourier index. Exploring such a basis (and relating this basis back
to the basis in eqs 15 and 40) might be very revealing as far as simulating
systems with both time and spatial symmetry.61

Looking forward, we can find several paths for future exploration with many
open questions. A first question is, as discussed in section 2, how to
properly calculate final electronic populations. Formally, whewn working
within a Floquet picture, one should account for non-zero interference terms
of the form: 44, 29

$\text{Prob}^{\text{interference}}_{\nu}=\sum_{n\neq
m}\frac{\sum_{r=1}^{N_{n\nu}^{\text{traj}}}\sum_{s=1}^{N_{m\nu}^{\text{traj}}}\tilde{c}^{r}_{n\nu}(\tilde{c}^{s}_{m\nu})^{*}\exp(i(n-m)\omega
t)}{{N^{\text{traj}}_{n\nu}\times N^{\text{traj}}_{m\nu}}}.$ (52)

Thus, in the context of F-PSSH, one may wonder if/how we might be able to
properly calculate the interference terms between wavepackets with different
momenta. This question should be addressed in the future.

A second and equally important question relates to the ubiquitous extra
decoherence that must be included within all surface hopping algorithms, i.e.
the need to account for wave packet separation. Do existing decoherence
ansatzes for FSSH needs to be modified when applied to PSSH approaches, and in
particular to F-PSSH?62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
76, 77, 78, 79 Our intuitions is that, as in standard FSSH, for F-PSSH the
coherence of a trajectory moving along one Floquet phase-space adiabat
trajectory should also be damped as wavepackets on different Floquet phase-
space adiabat separate after passing through the crossing region. But for any
Floquet based scheme, we have an interesting scenario whereby we have two
groups of parallel surfaces, with forces $F_{0}$ and $F_{1}$, and we know that
decoherence should be proportional to $F_{0}-F_{1}$. Thus, one might wonder:
should there be a unique decoherence rate for every pair of states,
$\left|\xi^{{}^{\prime}\tilde{a}0}\right\rangle$ and
$\left|\xi^{{}^{\prime}\tilde{b}1}\right\rangle$ for all possible $a,b$ photon
indices? Or should there be just one effective decoherence rate for the
system, capturing an effective rate of wavepacket separation between all those
wavepackets on Floquet surfaces $a0$ and all those wavepackets on Floquet
surfaces $b1$?

Detailed benchmarking needs to be done before drawing any definitive
conclusions.

Third and finally, the exact functional form of phase factors $\phi$ can be
obtained only by performing ab initio calculations. As far as the implications
of a circularly polarized light field and the corresponding Berry forces are
concerned, the most essential question is: how strong are the phase
oscillation ($\nabla_{\mathbf{R}}\phi$) for realistic systems and are these
oscillations localized or delocalized across configuration space. These
questions require further investigation as well.

In the end, although there are many questions left unanswered, the possibility
of using circularly polarized light to push forward coupled nuclear-electronic
dynamics is a tantalizing prospect with a largely underexplored knob to
control the dynamics (the ellipcity of the light). We believe the Floquet PSSH
approach proposed here should offer an essential tool going forward to begin
such exploration.

This work was supported by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences, under Award No. DE-SC0019397

##  Appendix A Diabatic Floquet Hamiltonian

In this appendix, for visual ease, we write out explicitly the Floquet
Hamiltonian in matrix form. Our basis $\\{\left|m\mu\right\rangle\\}$ is:

$\displaystyle\Bigl{\\{}\dots;\left|(-\tilde{1})0\right\rangle;\left|(-\tilde{1})1\right\rangle;\left|\tilde{0}0\right\rangle;\left|\tilde{0}1\right\rangle;\left|\tilde{1}0\right\rangle;\left|\tilde{1}1\right\rangle;\dots\Bigl{\\}}$
(53)

Here again, $m=\tilde{0},\pm\tilde{1},...$ represents the Floquet photon
indices and $\mu=0,1$ represents the electronic state indices. The explict
matrix form of eq 16 is (see also eq 6 and 7):

$\displaystyle\hat{\cal H}_{F}=$ (54)
$\displaystyle\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&H^{el}_{00}-\hbar\omega&D_{a}e^{i\phi_{a}}&0&D_{b}e^{i\phi_{b}}/2&0&0&\dots\\\
\dots&D_{a}e^{-i\phi_{a}}&H^{el}_{11}-\hbar\omega&D_{b}e^{i\phi_{b}}/2&0&0&0&\dots\\\
\dots&0&D_{b}e^{-i\phi_{b}}/2&H^{el}_{00}&D_{a}e^{i\phi_{a}}&0&D_{b}e^{i\phi_{b}}/2&\dots\\\
\dots&D_{b}e^{-i\phi_{b}}/2&0&D_{a}e^{-i\phi_{a}}&H^{el}_{11}&D_{b}e^{i\phi_{b}}/2&0&\dots\\\
\dots&0&0&0&D_{b}e^{-i\phi_{b}}/2&H^{el}_{00}+\hbar\omega&D_{a}e^{i\phi_{a}}&\dots\\\
\dots&0&0&D_{b}e^{-i\phi_{b}}/2&0&D_{a}e^{-i\phi_{a}}&H^{el}_{11}+\hbar\omega&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(63)

Note that we omit here all dependence on nuclear configuration $\mathbf{R}$.

##  Appendix B Constructing A Boosted Electronic Hamiltonian that is as Real
as Possible for a Phase-space Floquet Hamiltonian Approach

In section 2.2.4, we argued that the appropriate gauge transformation
$\hat{\cal U}_{F}$ must depend on the initial electronic state of a given
problem. This scenario differs from original PSSH algorithm in Ref. 39, where
the initial state was irrelevant. The difference arises now from the fact that
for the Floquet Hamiltonian as explicitly shown in eq 54, it is impossible to
define a gauge transformations such that the entire Hamiltonian becomes real-
valued. To see this, let us attempt to make the $6\times 6$ block in eq 54
real.

Let us denote the target $6$ basis functions as

$\displaystyle\left[\begin{array}[]{c}\left|\xi^{{}^{\prime}(-\tilde{1})0}\right\rangle\\\
\left|\xi^{{}^{\prime}(-\tilde{1})1}\right\rangle\\\
\left|\xi^{{}^{\prime}\tilde{0}0}\right\rangle\\\
\left|\xi^{{}^{\prime}\tilde{0}1}\right\rangle\\\
\left|\xi^{{}^{\prime}\tilde{1}0}\right\rangle\\\
\left|\xi^{{}^{\prime}\tilde{1}1}\right\rangle\\\ \end{array}\right]=\hat{\cal
U}_{F}\left[\begin{array}[]{c}\left|(-\tilde{1})0\right\rangle\\\
\left|(-\tilde{1})1\right\rangle\\\ \left|\tilde{0}0\right\rangle\\\
\left|\tilde{0}1\right\rangle\\\ \left|\tilde{1}0\right\rangle\\\
\left|\tilde{1}1\right\rangle\end{array}\right]=\left[\begin{array}[]{c c c c
c c}e^{i\eta_{1}}&&&&&\\\ &e^{i\eta_{2}}&&&&\\\ &&e^{i\eta_{3}}&&&\\\
&&&e^{i\eta_{4}}&&\\\ &&&&e^{i\eta_{5}}&\\\
&&&&&e^{i\eta_{6}}\end{array}\right]\left[\begin{array}[]{c}\left|(-\tilde{1})0\right\rangle\\\
\left|(-\tilde{1})1\right\rangle\\\ \left|\tilde{0}0\right\rangle\\\
\left|\tilde{0}1\right\rangle\\\ \left|\tilde{1}0\right\rangle\\\
\left|\tilde{1}1\right\rangle\end{array}\right]$ (88)

In order to reach a real-valued Hamiltonian, the relative phase differences
between the basis functions need to satisfy the following overdetermined set
of algebraic equations:

$\displaystyle\left[\begin{array}[]{c c c c c c}1&-1&&&&\\\ &&1&-1&&\\\
&&&&1&-1\\\ &1&-1&&&\\\ 1&&&-1&&\\\ &&&1&-1&\\\
&&1&&&-1\end{array}\right]\left[\begin{array}[]{c}\eta_{1}\\\ \eta_{2}\\\
\eta_{3}\\\ \eta_{4}\\\ \eta_{5}\\\
\eta_{6}\end{array}\right]=\left[\begin{array}[]{c}\phi_{a}\\\ \phi_{a}\\\
\phi_{a}\\\ \phi_{b}\\\ \phi_{b}\\\ \phi_{b}\\\ \phi_{b}\end{array}\right]$
(109)

The first three rows arise from forcing the complex-valued vibronic couplings
to be real-valued, and the last four rows arise from forcing the complex-
valued light-matter couplings to be real-valued. However, there is no solution
to these algebraic equations. For instance, if we sum over rows $1,2,4$ and
compare the sum with row $5$, we find:

$\displaystyle\eta_{1}-\eta_{2}+\eta_{3}-\eta_{4}+\eta_{2}-\eta_{3}$
$\displaystyle=\eta_{1}-\eta_{4}=2\phi_{a}+\phi_{b}$ (110)
$\displaystyle\eta_{1}-\eta_{4}$ $\displaystyle=\phi_{b}$ (111)

For $\phi_{a}\neq k\pi$, $k\in\mathbb{Z}$, there is no solution to this set of
algebraic equations. Similarly, when we sum over rows $2,3,6$ and compare the
sum with row $7$, we obtain the same type of contradiction. In the end, there
simply is no gauge transformation ${\cal U}_{F}$ under which the (infinite-
dimensional) Floquet Hamiltonian becomes real-valued.

With this constraint in mind, a practical approach forward is to recognize
that only a few Floquet states are actually populated during a typical surface
hopping calculation. After all, for this same reason, the formally infinite
dimensional Floquet Hamiltonian can be safely truncated; for our calculations
in Fig. 1 \- 4 , the corresponding Floquet Hamiltonians are $18\times 18$
matrices. Thus, it makes sense for us to concern ourselves and make real-
valued only those states that directly couple to the initial Floquet state.

For the scenario that the initial state corresponds to
$\left|\tilde{0}0\right\rangle$, we choose to make the couplings to
$\left|\tilde{0}1\right\rangle$ and $\left|\tilde{-1}1\right\rangle$ real. If
we fix $\eta_{3}=0$, the corresponding matrix of phases is:

$\displaystyle\text{diag}(\hat{\cal U}_{F})_{\text{initial
state}=\left|\tilde{0}0\right\rangle}$
$\displaystyle=[\dots,e^{i\eta_{1}},e^{i\phi_{b}},1,e^{-i\phi_{a}},e^{i\eta_{5}},e^{i\eta_{6}},\dots]$
(112)

The choice of $\eta_{1}$, $\eta_{5}$, $\eta_{6}$ is irrelevant to the results.
Simply for the sake of concreteness, we will choose the phase $\eta_{1}$ for
states $\left|(-\tilde{1})0\right\rangle$ so as to make the coupling
$\left\langle(-\tilde{1})0\right|\hat{\cal
H}_{F}\left|(-\tilde{1})1\right\rangle$ real-valued, we will choose the phase
$\eta_{5}$ for state $\left|\tilde{1}0\right\rangle$ so as to make the
coupling $\left\langle\tilde{0}1\right|\hat{\cal
H}_{F}\left|\tilde{1}0\right\rangle$ real-valued and we choose the phase
$\eta_{6}$ for state $\left|\tilde{1}1\right\rangle$ so as to make the
coupling $\left\langle\tilde{0}0\right|\hat{\cal
H}_{F}\left|\tilde{1}1\right\rangle$ to be real-valued. The final result is:

$\displaystyle\text{diag}(\hat{\cal U}_{F})_{\text{initial
state}=\left|\tilde{0}0\right\rangle}=[\dots,e^{i\phi_{b}+i\phi_{a}},e^{i\phi_{b}},1,e^{-i\phi_{a}},e^{-i\phi_{a}-i\phi_{b}},e^{-i\phi_{b}},\dots]$
(113)

For the couplings between the rest of the states with larger Floquet photon
indices, we follow the same fashion such that the closest complex-valued
couplings to the initial state is transformed to be real-valued.

Naturally, there is a different transformation if are to simulate dynamics
with $\left|\tilde{0}1\right\rangle$ as the initial state. Now, if we fix
$\eta_{4}=0$, we obtain

$\displaystyle\text{diag}(\hat{\cal U}_{F})_{\text{initial
state}=\left|\tilde{0}1\right\rangle}$
$\displaystyle=[\dots,e^{i\eta_{1}},e^{i\eta_{2}},e^{i\phi_{a}},1,e^{-i\phi_{b}},e^{i\eta_{6}},\dots]$
$\displaystyle=[\dots,e^{i\phi_{b}},e^{i\phi_{a}+i\phi_{b}},e^{i\phi_{a}},1,e^{-i\phi_{b}},e^{-i\phi_{b}-i\phi_{a}},\dots]$
(114)

As one would expect, the final results do not depend on $\eta_{1}$,
$\eta_{2}$, $\eta_{6}$.

The matrices $\hat{\cal D}_{F}$ for these two scenarios are obviously:

$\displaystyle\text{diag}(\hat{\cal D}_{F})_{\text{initial
state}=\left|\tilde{0}0\right\rangle}$
$\displaystyle=[\dots;-\mathbf{w}_{a}-\mathbf{w}_{b};-\mathbf{w}_{b};0;\mathbf{w}_{a};\mathbf{w}_{a}+\mathbf{w}_{b};\mathbf{w}_{b};\dots]$
(115) $\displaystyle\text{diag}(\hat{\cal D}_{F})_{\text{initial
state}=\left|\tilde{0}1\right\rangle}$
$\displaystyle=[\dots-\mathbf{w}_{b};-\mathbf{w}_{a}-\mathbf{w}_{b};-\mathbf{w}_{a};0;\mathbf{w}_{b};\mathbf{w}_{a}+\mathbf{w}_{b};\dots]$
(116)

Under these gauge transformations, we obtain two diabatic Floquet Hamiltonians
respectively:

$\displaystyle(\hat{\cal H}_{F})_{\text{initial
state}=\left|\tilde{0}0\right\rangle}=$ (117)
$\displaystyle\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&H^{el}_{00}-\hbar\omega&D_{a}&0&D_{b}e^{-2i\phi_{a}}/2&0&0&\dots\\\
\dots&D_{a}&H^{el}_{11}-\hbar\omega&D_{b}/2&0&0&0&\dots\\\
\dots&0&D_{b}/2&H^{el}_{00}&D_{a}&0&D_{b}/2&\dots\\\
\dots&D_{b}e^{2i\phi_{a}}/2&0&D_{a}&H^{el}_{11}&D_{b}/2&0&\dots\\\
\dots&0&0&0&D_{b}/2&H^{el}_{00}+\hbar\omega&D_{a}e^{2i\phi_{a}}&\dots\\\
\dots&0&0&D_{b}/2&0&D_{a}e^{-2i\phi_{a}}&H^{el}_{11}+\hbar\omega&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(126)

$\displaystyle(\hat{\cal H}_{F})_{\text{initial
state}=\left|\tilde{0}1\right\rangle}=$ (127)
$\displaystyle\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&H^{el}_{00}-\hbar\omega&D_{a}e^{2i\phi_{a}}&0&D_{b}/2&0&0&\dots\\\
\dots&D_{a}e^{-2i\phi_{a}}&H^{el}_{11}-\hbar\omega&D_{b}/2&0&0&0&\dots\\\
\dots&0&D_{b}/2&H^{el}_{00}&D_{a}&0&D_{b}e^{-2i\phi_{a}}/2&\dots\\\
\dots&D_{b}/2&0&D_{a}&H^{el}_{11}&D_{b}/2&0&\dots\\\
\dots&0&0&0&D_{b}/2&H^{el}_{00}+\hbar\omega&D_{a}&\dots\\\
\dots&0&0&D_{b}e^{2i\phi_{a}}/2&0&D_{a}&H^{el}_{11}+\hbar\omega&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(136)

Note that there is a real-valued 3x3 matrix block inside of each Floquet
Hamiltonian.

Lastly, by following eq 44, if the initial state is
$\left|\tilde{0}1\right\rangle$, the final diabatic phase-space Floquet
Hamiltonian is:

$\displaystyle[\hat{\cal H}_{F}]^{{}^{\prime}}(\mathbf{R},\mathbf{P})=$ (137)
$\displaystyle\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&H^{el}_{00}-\hbar\omega&D_{a}&0&D_{b}e^{-2i\phi_{a}}/2&0&0&\dots\\\
\dots&D_{a}&H^{el}_{11}-\hbar\omega&D_{b}/2&0&0&0&\dots\\\
\dots&0&D_{b}/2&H^{el}_{00}&D_{a}&0&D_{b}/2&\dots\\\
\dots&D_{b}e^{2i\phi_{a}}/2&0&D_{a}&H^{el}_{11}&D_{b}/2&0&\dots\\\
\dots&0&0&0&D_{b}/2&H^{el}_{00}+\hbar\omega&D_{a}e^{2i\phi_{a}}&\dots\\\
\dots&0&0&D_{b}/2&0&D_{a}e^{-2i\phi_{a}}&H^{el}_{11}+\hbar\omega&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(146)
$\displaystyle-\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&\frac{\hbar^{2}(\mathbf{w}_{a}+\mathbf{w}_{b})^{2}}{2m}&0&0&0&0&0&\dots\\\
\dots&0&\frac{\hbar^{2}\mathbf{w}_{b}^{2}}{2m}&0&0&0&0&\dots\\\
\dots&0&0&0&0&0&0&\dots\\\
\dots&0&0&0&\frac{\hbar^{2}\mathbf{w}_{a}^{2}}{2m}&0&0&\dots\\\
\dots&0&0&0&0&\frac{\hbar^{2}(\mathbf{w}_{a}+\mathbf{w}_{b})^{2}}{2m}&0&\dots\\\
\dots&0&0&0&0&0&\frac{\hbar^{2}\mathbf{w}_{b}^{2}}{2m}&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(155)
$\displaystyle+\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&\frac{i\hbar\mathbf{P}\cdot(\mathbf{w}_{a}+\mathbf{w}_{b})}{m}&0&0&0&0&0&\dots\\\
\dots&0&\frac{i\hbar\mathbf{P}\cdot\mathbf{w}_{b}}{m}&0&0&0&0&\dots\\\
\dots&0&0&0&0&0&0&\dots\\\
\dots&0&0&0&-\frac{i\hbar\mathbf{P}\cdot\mathbf{w}_{a}}{m}&0&0&\dots\\\
\dots&0&0&0&0&-\frac{i\hbar\mathbf{P}\cdot(\mathbf{w}_{a}+\mathbf{w}_{b})}{m}&0&\dots\\\
\dots&0&0&0&0&0&-\frac{i\hbar\mathbf{P}\cdot\mathbf{w}_{b}}{m}&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(164)

If the initial state is $\left|\tilde{0}1\right\rangle$, then the final
diabatic phase-space Floquet Hamiltonian is:

$\displaystyle[\hat{\cal H}_{F}]^{{}^{\prime}}(\mathbf{R},\mathbf{P})=$ (165)
$\displaystyle\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&H^{el}_{00}-\hbar\omega&D_{a}e^{2i\phi_{a}}&0&D_{b}/2&0&0&\dots\\\
\dots&D_{a}e^{-2i\phi_{a}}&H^{el}_{11}-\hbar\omega&D_{b}/2&0&0&0&\dots\\\
\dots&0&D_{b}/2&H^{el}_{00}&D_{a}&0&D_{b}e^{-2i\phi_{a}}/2&\dots\\\
\dots&D_{b}/2&0&D_{a}&H^{el}_{11}&D_{b}/2&0&\dots\\\
\dots&0&0&0&D_{b}/2&H^{el}_{00}+\hbar\omega&D_{a}&\dots\\\
\dots&0&0&D_{b}e^{2i\phi_{a}}/2&0&D_{a}&H^{el}_{11}+\hbar\omega&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(174)
$\displaystyle-\left(\begin{array}[]{cccccccc}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&\frac{\hbar^{2}\mathbf{w}_{b}^{2}}{2m}&0&0&0&0&0&\dots\\\
\dots&0&\frac{\hbar^{2}(\mathbf{w}_{a}+\mathbf{w}_{b})^{2}}{2m}&0&0&0&0&\dots\\\
\dots&0&0&\frac{\hbar^{2}\mathbf{w}_{a}^{2}}{2m}&0&0&0&\dots\\\
\dots&0&0&0&0&0&0&\dots\\\
\dots&0&0&0&0&\frac{\hbar^{2}\mathbf{w}_{b}^{2}}{2m}&0&\dots\\\
\dots&0&0&0&0&0&\frac{\hbar^{2}(\mathbf{w}_{a}+\mathbf{w}_{b})^{2}}{2m}&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(183) $\displaystyle+\left(\begin{array}[]{c c c c c c c
c}\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\reflectbox{$\ddots$}\\\
\dots&\frac{i\hbar\mathbf{P}\cdot\mathbf{w}_{b}}{m}&0&0&0&0&0&\dots\\\
\dots&0&\frac{i\hbar\mathbf{P}\cdot(\mathbf{w}_{a}+\mathbf{w}_{b})}{m}&0&0&0&0&\dots\\\
\dots&0&0&\frac{i\hbar\mathbf{P}\cdot\mathbf{w}_{a}}{m}&0&0&0&\dots\\\
\dots&0&0&0&0&0&0&\dots\\\
\dots&0&0&0&0&-\frac{i\hbar\mathbf{P}\cdot\mathbf{w}_{b}}{m}&0&\dots\\\
\dots&0&0&0&0&0&-\frac{i\hbar\mathbf{P}\cdot(\mathbf{w}_{a}+\mathbf{w}_{b})}{m}&\dots\\\
\reflectbox{$\ddots$}&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$
(192)

## References

  * Mai and González 2020 Mai, S.; González, L. Molecular photochemistry: Recent developments in theory. _Angewandte Chemie International Edition_ 2020, _59_ , 16832–16846
  * Nelson et al. 2020 Nelson, T. R.; White, A. J.; Bjorgaard, J. A.; Sifain, A. E.; Zhang, Y.; Nebgen, B.; Fernandez-Alberti, S.; Mozyrsky, D.; Roitberg, A. E.; Tretiak, S. Non-adiabatic excited-state molecular dynamics: Theory and applications for modeling photophysics in extended molecular materials. _Chemical Reviews_ 2020, _120_ , 2215–2287
  * Stock and Thoss 2005 Stock, G.; Thoss, M. Classical description of nonadiabatic quantum dynamics. _Advances in chemical physics_ 2005, _131_ , 243–376
  * Levine and Martínez 2007 Levine, B. G.; Martínez, T. J. Isomerization through conical intersections. _Annu. Rev. Phys. Chem._ 2007, _58_ , 613–634
  * Penfold et al. 2018 Penfold, T. J.; Gindensperger, E.; Daniel, C.; Marian, C. M. Spin-vibronic mechanism for intersystem crossing. _Chemical reviews_ 2018, _118_ , 6975–7025
  * Moiseyev et al. 2008 Moiseyev, N.; Šindelka, M.; Cederbaum, L. S. Laser-induced conical intersections in molecular optical lattices. _Journal of Physics B: Atomic, Molecular and Optical Physics_ 2008, _41_ , 221001
  * Halász et al. 2011 Halász, G. J.; Vibók, Á.; Šindelka, M.; Moiseyev, N.; Cederbaum, L. S. Conical intersections induced by light: Berry phase and wavepacket dynamics. _Journal of Physics B: Atomic, Molecular and Optical Physics_ 2011, _44_ , 175102
  * Halász et al. 2012 Halász, G. J.; Šindelka, M.; Moiseyev, N.; Cederbaum, L. S.; Vibók, Á. Light-induced conical intersections: Topological phase, wave packet dynamics, and molecular alignment. _The Journal of Physical Chemistry A_ 2012, _116_ , 2636–2643
  * Corrales et al. 2014 Corrales, M.; González-Vázquez, J.; Balerdi, G.; Solá, I.; De Nalda, R.; Bañares, L. Control of ultrafast molecular photodissociation by laser-field-induced potentials. _Nature chemistry_ 2014, _6_ , 785–790
  * Natan et al. 2016 Natan, A.; Ware, M. R.; Prabhudesai, V. S.; Lev, U.; Bruner, B. D.; Heber, O.; Bucksbaum, P. H. Observation of quantum interferences via light-induced conical intersections in diatomic molecules. _Physical Review Letters_ 2016, _116_ , 143004
  * Csehi et al. 2017 Csehi, A.; Halász, G. J.; Cederbaum, L. S.; Vibók, Á. Intrinsic and light-induced nonadiabatic phenomena in the NaI molecule. _Physical Chemistry Chemical Physics_ 2017, _19_ , 19656–19664
  * Leclerc et al. 2017 Leclerc, A.; Viennot, D.; Jolicard, G.; Lefebvre, R.; Atabek, O. Exotic states in the strong-field control of dissociation dynamics: from exceptional points to zero-width resonances. _Journal of Physics B: Atomic, Molecular and Optical Physics_ 2017, _50_ , 234002
  * Csehi et al. 2017 Csehi, A.; Halász, G. J.; Cederbaum, L. S.; Vibók, A. Competition between light-induced and intrinsic nonadiabatic phenomena in diatomics. _The Journal of Physical Chemistry Letters_ 2017, _8_ , 1624–1630
  * Szidarovszky et al. 2018 Szidarovszky, T.; Halász, G. J.; Császár, A. G.; Cederbaum, L. S.; Vibók, A. Direct signatures of light-induced conical intersections on the field-dressed spectrum of Na2. _The Journal of Physical Chemistry Letters_ 2018, _9_ , 2739–2745
  * Triana and Sanz-Vicario 2019 Triana, J. F.; Sanz-Vicario, J. L. Revealing the presence of potential crossings in diatomics induced by quantum cavity radiation. _Physical review letters_ 2019, _122_ , 063603
  * Kübel et al. 2020 Kübel, M.; Spanner, M.; Dube, Z.; Naumov, A. Y.; Chelkowski, S.; Bandrauk, A. D.; Vrakking, M. J.; Corkum, P. B.; Villeneuve, D.; Staudte, A. Probing multiphoton light-induced molecular potentials. _Nature communications_ 2020, _11_ , 1–8
  * Zanchet et al. 2020 Zanchet, A.; García, G. A.; Nahon, L.; Bañares, L.; Poullain, S. M. Signature of a conical intersection in the dissociative photoionization of formaldehyde. _Physical Chemistry Chemical Physics_ 2020, _22_ , 12886–12893
  * Farag et al. 2021 Farag, M. H.; Mandal, A.; Huo, P. Polariton induced conical intersection and berry phase. _Physical Chemistry Chemical Physics_ 2021, _23_ , 16868–16879
  * Fábri et al. 2021 Fábri, C.; Halász, G. J.; Cederbaum, L. S.; Vibók, Á. Signatures of light-induced nonadiabaticity in the field-dressed vibronic spectrum of formaldehyde. _The Journal of Chemical Physics_ 2021, _154_ , 124308
  * Tully 1990 Tully, J. C. Molecular dynamics with electronic transitions. _The Journal of Chemical Physics_ 1990, _93_ , 1061–1071
  * Mitrić et al. 2009 Mitrić, R.; Petersen, J.; Bonačić-Kouteckỳ, V. Laser-field-induced surface-hopping method for the simulation and control of ultrafast photodynamics. _Physical Review A_ 2009, _79_ , 053416
  * Mitrić et al. 2011 Mitrić, R.; Petersen, J.; Wohlgemuth, M.; Werner, U.; Bonačić-Kouteckỳ, V. Field-induced surface hopping method for probing transition state nonadiabatic dynamics of Ag 3. _Physical Chemistry Chemical Physics_ 2011, _13_ , 8690–8696
  * Lisinetskaya and Mitrić 2011 Lisinetskaya, P. G.; Mitrić, R. Simulation of laser-induced coupled electron-nuclear dynamics and time-resolved harmonic spectra in complex systems. _Physical Review A_ 2011, _83_ , 033408
  * Richter et al. 2011 Richter, M.; Marquetand, P.; González-Vázquez, J.; Sola, I.; González, L. SHARC: ab initio molecular dynamics with surface hopping in the adiabatic representation including arbitrary couplings. _Journal of chemical theory and computation_ 2011, _7_ , 1253–1258
  * Bajo et al. 2012 Bajo, J. J.; González-Vázquez, J.; Sola, I. R.; Santamaria, J.; Richter, M.; Marquetand, P.; González, L. Mixed quantum-classical dynamics in the adiabatic representation to simulate molecules driven by strong laser pulses. _J. Phys. Chem. A_ 2012, _116_ , 2800–2807
  * Mai et al. 2015 Mai, S.; Marquetand, P.; González, L. A general method to describe intersystem crossing dynamics in trajectory surface hopping. _International Journal of Quantum Chemistry_ 2015, _115_ , 1215–1231
  * Fiedlschuster et al. 2016 Fiedlschuster, T.; Handt, J.; Schmidt, R. Floquet surface hopping: Laser-driven dissociation and ionization dynamics of H 2+. _Phys. Rev. A_ 2016, _93_ , 053409
  * Fiedlschuster et al. 2017 Fiedlschuster, T.; Handt, J.; Gross, E.; Schmidt, R. Surface hopping in laser-driven molecular dynamics. _Physical Review A_ 2017, _95_ , 063424
  * Zhou et al. 2020 Zhou, Z.; Chen, H.-T.; Nitzan, A.; Subotnik, J. E. Nonadiabatic dynamics in a laser field: Using Floquet fewest switches surface hopping to calculate electronic populations for slow nuclear velocities. _Journal of Chemical Theory and Computation_ 2020, _16_ , 821–834
  * Hone et al. 1997 Hone, D. W.; Ketzmerick, R.; Kohn, W. Time-dependent Floquet theory and absence of an adiabatic limit. _Phys. Rev. A_ 1997, _56_ , 4045–4054
  * Grifoni and Hänggi 1998 Grifoni, M.; Hänggi, P. Driven quantum tunneling. _Physics Reports_ 1998, _304_ , 229–354
  * Jasper and Truhlar 2003 Jasper, A. W.; Truhlar, D. G. Improved treatment of momentum at classically forbidden electronic transitions in trajectory surface hopping calculations. _Chem. Phys. Lett._ 2003, _369_ , 60–67
  * Miao et al. 2019 Miao, G.; Bellonzi, N.; Subotnik, J. An extension of the fewest switches surface hopping algorithm to complex Hamiltonians and photophysics in magnetic fields: Berry curvature and “magnetic” forces. _The Journal of chemical physics_ 2019, _150_ , 124101
  * Miao et al. 2020 Miao, G.; Bian, X.; Zhou, Z.; Subotnik, J. A “backtracking” correction for the fewest switches surface hopping algorithm. _The Journal of Chemical Physics_ 2020, _153_ , 111101
  * Bian et al. 2021 Bian, X.; Wu, Y.; Teh, H.-H.; Zhou, Z.; Chen, H.-T.; Subotnik, J. E. Modeling nonadiabatic dynamics with degenerate electronic states, intersystem crossing, and spin separation: A key goal for chemical physics. _The Journal of Chemical Physics_ 2021, _154_ , 110901
  * Bian et al. 2022 Bian, X.; Qiu, T.; Chen, J.; Subotnik, J. E. On the meaning of Berry force for unrestricted systems treated with mean-field electronic structure. _The Journal of Chemical Physics_ 2022, _156_ , 234107
  * Bian et al. 2022 Bian, X.; Wu, Y.; Teh, H.-H.; Subotnik, J. E. Incorporating Berry Force Effects into the Fewest Switches Surface-Hopping Algorithm: Intersystem Crossing and the Case of Electronic Degeneracy. _Journal of Chemical Theory and Computation_ 2022, _18_ , 2075–2090
  * Wu and Subotnik 2021 Wu, Y.; Subotnik, J. E. Semiclassical description of nuclear dynamics moving through complex-valued single avoided crossings of two electronic states. _The Journal of Chemical Physics_ 2021, _154_ , 234101
  * Wu et al. 2022 Wu, Y.; Bian, X.; Rawlinson, J. I.; Littlejohn, R. G.; Subotnik, J. E. A phase-space semiclassical approach for modeling nonadiabatic nuclear dynamics with electronic spin. _The Journal of chemical physics_ 2022, _157_ , 011101
  * Naaman and Waldeck 2012 Naaman, R.; Waldeck, D. H. Chiral-induced spin selectivity effect. _The journal of physical chemistry letters_ 2012, _3_ , 2178–2187
  * Fransson 2020 Fransson, J. Vibrational origin of exchange splitting and” chiral-induced spin selectivity. _Physical Review B_ 2020, _102_ , 235416
  * Hore et al. 2020 Hore, P.; Ivanov, K. L.; Wasielewski, M. R. Spin chemistry. 2020
  * Mai et al. 2018 Mai, S.; Marquetand, P.; González, L. Nonadiabatic dynamics: The SHARC approach. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2018, _8_ , e1370
  * Ho et al. 1983 Ho, T.-S.; Chu, S.-I.; Tietz, J. V. Semiclassical many-mode Floquet theory. _Chem. Phys. Lett._ 1983, _96_ , 464–471
  * Subotnik et al. 2016 Subotnik, J. E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W.; Bellonzi, N. Understanding the surface hopping view of electronic transitions and decoherence. _Annu. Rev. Phys. Chem._ 2016, _67_ , 387–417
  * 46 This formula ignores the contribution from interferences between different Floquet states that correspond to the same electronic state.29 However, this approximate formula is fairly valid for our model problem because the phases of wavepackets along different Floquet states are significantly different (due to the complex-valued Hamiltonians that lead to different asymptotic momentum). Hence, in practice, we find the contribution from interference terms to the final transmission probabilities to be very small.
  * 47 Note that the $angle\phi$ depends on the nuclear configuration $\mathbf{R}$ and there is no diabatic basis in which this complex-valued Hamiltonian becomes real-valued Hamiltonian.
  * Jain and Subotnik 2015 Jain, A.; Subotnik, J. E. Surface hopping, transition state theory, and decoherence. II. Thermal rate constants and detailed balance. _J. Chem. Phys._ 2015, _143_ , 134107
  * Jain et al. 2016 Jain, A.; Alguire, E.; Subotnik, J. E. An efficient, augmented surface hopping algorithm that includes decoherence for use in large-scale simulations. _J. Chem. Theory Comput._ 2016, _12_ , 5256–5268
  * Hammes-Schiffer and Tully 1994 Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. _J. Chem. Phys._ 1994, _101_ , 4657–4667
  * Fabiano et al. 2008 Fabiano, E.; Keal, T.; Thiel, W. Implementation of surface hopping molecular dynamics using semiempirical methods. _Chem. Phys._ 2008, _349_ , 334–347
  * Barbatti et al. 2010 Barbatti, M.; Pittner, J.; Pederzoli, M.; Werner, U.; Mitrić, R.; Bonačić-Kouteckỳ, V.; Lischka, H. Non-adiabatic dynamics of pyrrole: Dependence of deactivation mechanisms on the excitation energy. _Chem. Phys._ 2010, _375_ , 26–34
  * Plasser et al. 2012 Plasser, F.; Granucci, G.; Pittner, J.; Barbatti, M.; Persico, M.; Lischka, H. Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer. _J. Chem. Phys._ 2012, _137_ , 22A514
  * Fernandez-Alberti et al. 2012 Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules. _J. Chem. Phys._ 2012, _137_ , 014512
  * Nelson et al. 2013 Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Artifacts due to trivial unavoided crossings in the modeling of photoinduced energy transfer dynamics in extended conjugated molecules. _Chem. Phys. Lett._ 2013, _590_ , 208–213
  * Wang and Prezhdo 2014 Wang, L.; Prezhdo, O. V. A simple solution to the trivial crossing problem in surface hopping. _J. Phys. Chem. Lett._ 2014, _5_ , 713–719
  * Meek and Levine 2014 Meek, G. A.; Levine, B. G. Evaluation of the time-derivative coupling for accurate electronic state transition probabilities from numerical simulations. _J. Phys. Chem. Lett._ 2014, _5_ , 2351–2356
  * Dell’Angelo and Hanna 2017 Dell’Angelo, D.; Hanna, G. Importance of eigenvector sign consistency in computations of expectation values via mixed quantum-classical surface-hopping dynamics. _Theor. Chem. Acc._ 2017, _136_ , 75
  * Lee and Willard 2019 Lee, E. M.; Willard, A. P. Solving the trivial crossing problem while preserving the nodal symmetry of the wavefunction. _J. Chem. Theory Comput._ 2019,
  * Zhou et al. 2019 Zhou, Z.; Jin, Z.; Qiu, T.; Rappe, A. M.; Subotnik, J. E. A Robust and Unified Solution for Choosing the Phases of Adiabatic States as a Function of Geometry: Extending Parallel Transport Concepts to the Cases of Trivial and Near-Trivial Crossings. _Journal of Chemical Theory and Computation_ 2019, _16_ , 835–846
  * Fläschner et al. 2016 Fläschner, N.; Rem, B.; Tarnowski, M.; Vogel, D.; Lühmann, D.-S.; Sengstock, K.; Weitenberg, C. Experimental reconstruction of the Berry curvature in a Floquet Bloch band. _Science_ 2016, _352_ , 1091–1094
  * Schwartz and Rossky 1994 Schwartz, B. J.; Rossky, P. J. Aqueous solvation dynamics with a quantum mechanical solute: computer simulation studies of the photoexcited hydrated electron. _J. Chem. Phys._ 1994, _101_ , 6902–6916
  * Bittner and Rossky 1995 Bittner, E. R.; Rossky, P. J. Quantum decoherence in mixed quantum-classical systems: Nonadiabatic processes. _J. Chem. Phys._ 1995, _103_ , 8130–8143
  * Schwartz et al. 1996 Schwartz, B. J.; Bittner, E. R.; Prezhdo, O. V.; Rossky, P. J. Quantum decoherence and the isotope effect in condensed phase nonadiabatic molecular dynamics simulations. _J. Chem. Phys._ 1996, _104_ , 5942–5955
  * Prezhdo and Rossky 1997 Prezhdo, O. V.; Rossky, P. J. Evaluation of quantum transition rates from quantum-classical molecular dynamics simulations. _J. Chem. Phys._ 1997, _107_ , 5863–5878
  * Prezhdo and Rossky 1997 Prezhdo, O. V.; Rossky, P. J. Mean-field molecular dynamics with surface hopping. _J. Chem. Phys._ 1997, _107_ , 825–834
  * Fang and Hammes-Schiffer 1999 Fang, J.-Y.; Hammes-Schiffer, S. Improvement of the internal consistency in trajectory surface hopping. _J. Phys. Chem. A_ 1999, _103_ , 9399–9407
  * Fang and Hammes-Schiffer 1999 Fang, J.-Y.; Hammes-Schiffer, S. Comparison of surface hopping and mean field approaches for model proton transfer reactions. _J. Chem. Phys._ 1999, _110_ , 11166–11175
  * Volobuev et al. 2000 Volobuev, Y. L.; Hack, M. D.; Topaler, M. S.; Truhlar, D. G. Continuous surface switching: An improved time-dependent self-consistent-field method for nonadiabatic dynamics. _J. Chem. Phys._ 2000, _112_ , 9716–9726
  * Hack and Truhlar 2001 Hack, M. D.; Truhlar, D. G. Electronically nonadiabatic trajectories: Continuous surface switching II. _J. Chem. Phys._ 2001, _114_ , 2894–2902
  * Horenko et al. 2001 Horenko, I.; Schmidt, B.; Schütte, C. A theoretical model for molecules interacting with intense laser pulses: The Floquet-based quantum-classical Liouville equation. _The Journal of Chemical Physics_ 2001, _115_ , 5733–5743
  * Wong and Rossky 2002 Wong, K. F.; Rossky, P. J. Dissipative mixed quantum-classical simulation of the aqueous solvated electron system. _J. Chem. Phys._ 2002, _116_ , 8418–8428
  * Wong and Rossky 2002 Wong, K. F.; Rossky, P. J. Solvent-induced electronic decoherence: Configuration dependent dissipative evolution for solvated electron systems. _J. Chem. Phys._ 2002, _116_ , 8429–8438
  * Horenko et al. 2002 Horenko, I.; Salzmann, C.; Schmidt, B.; Schütte, C. Quantum-classical Liouville approach to molecular dynamics: Surface hopping Gaussian phase-space packets. _J. Chem. Phys._ 2002, _117_ , 11075–11088
  * Jasper and Truhlar 2005 Jasper, A. W.; Truhlar, D. G. Electronic decoherence time for non-Born-Oppenheimer trajectories. _J. Chem. Phys._ 2005, _123_ , 064103
  * Bedard-Hearn et al. 2005 Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. Mean-field dynamics with stochastic decoherence (MF-SD): A new algorithm for nonadiabatic mixed quantum/classical molecular-dynamics simulations with nuclear-induced decoherence. _J. Chem. Phys._ 2005, _123_ , 234106
  * Subotnik and Shenvi 2011 Subotnik, J. E.; Shenvi, N. Decoherence and surface hopping: When can averaging over initial conditions help capture the effects of wave packet separation? _J. Chem. Phys._ 2011, _134_ , 244114
  * Subotnik and Shenvi 2011 Subotnik, J. E.; Shenvi, N. A new approach to decoherence and momentum rescaling in the surface hopping algorithm. _J. Chem. Phys._ 2011, _134_ , 024105
  * Landry and Subotnik 2012 Landry, B. R.; Subotnik, J. E. How to recover Marcus theory with fewest switches surface hopping: Add just a touch of decoherence. _J. Chem. Phys._ 2012, _137_ , 22A513