# Unidirectional coupling of a quantum emitter to a subwavelength grating
waveguide with engineered stationary inflection point: Supporting Information

Ilya. A. Volkov  Department of Physics and Engineering, ITMO University,
St.-Petersburg 197101, Russia  Roman S. Savelev  Department of Physics and
Engineering, ITMO University, St.-Petersburg 197101, Russia
<EMAIL_ADDRESS>
###### Contents

  1. I A. Properties of the isolated waveguides
  2. I.1 General formulation
  3. I.2 Application of the CMT to the first design
  4. I.3 Application of the CMT to the second design

##  I A. Properties of the isolated waveguides

The structures shown schematically in Figs. 3(a,d) in the main text are based
on periodic waveguides that support two quasi-degenerate modes, coupled with a
waveguide with a single mode in the spectral range of interest. The
possibility to achieve a dispersion crossing of two modes with orthogonal
dominant polarizations was demonstrated previously in the Refs. Savelev  _et
al._ (2019); Savelev and Gorlach (2020). Such modes can be considered as two
constituent elements that can be combined in an arbitrary way allowing for
controllable tailoring of local polarization of the electric field in the
plane parallel to the substrate.

Figure S1: (a) Dispersion of a periodic waveguide with two intersecting modes
with positive group velocities (top view), $\lambda$ is the wavelength, $k$ is
the Bloch wavenumber. Thickness of the waveguide membrane $h=220$ nm, period
is $a=270$ nm, dimensions of the rectangular hole are $100\times 310$ nm,
thickness of the waveguide $w_{1}=390$ nm. (c-f) Intensities of $y$ (c,e) and
$z$ (d,f) electric field in the cross section of the unit cell in the plane
$x=90$ nm for the (c,d) first modes, (e,f) second mode.

In the first design we consider a periodic waveguide with two forward quasi-
degenerate modes. Their dispersion is shown with red and blue curves in Fig.
S1. The dispersion of the modes can be approximated as following:

$\displaystyle\omega_{1}\;\text{[rad/s]}\approx 2\pi\cdot 308.5\cdot
10^{12}+5.1\cdot 10^{7}\delta k-21\delta k^{2},$ (S1)
$\displaystyle\omega_{2}\;\text{[rad/s]}\approx 2\pi\cdot 308.5\cdot
10^{12}+2.3\cdot 10^{7}\delta k-19\delta k^{2},$

where $\delta k$ is the detuning of the wavvenumber from the crossing point.

Both of these modes possess two planes of symmetry: electric (tangential
components of electric field is zero) plane of symmetry $x=0$, and electric
(magnetic) plane of symmetry $y=0$ for the first (second) mode. Consequently,
in the plane $x=0$ electric field in both modes has only $x$ component.
However, when shifted by any distance in $x$ direction, one mode acquires
(dominant) $y$ component of electric field, while another one $z$ component.
Distributions of the $y$ and $z$ components of electric field at the distance
$x=90$ nm are shown for both modes in Fig. S1.

Figure S2: Dispersion of a periodic waveguide with two intersecting modes “1”
and “2” with positive group velocities (red and blue curves) and a dispersion
of a side-coupled strip waveguide (dashed green curve) with mode “3”. Due to
the dispersion folding, dispersion of the third mode reflects from the edge of
the Brillouin zone ($ka/\pi=1$), as shown with the solid black curve, changes
the sign of the group velocity and intersects with the first and second modes.

In order to achieve a zero group velocity of the supermode of the coupled
waveguides, one has to compensate the positive group velocities of these modes
by another mode with negative group velocity. To this end we employ a
homogeneous strip waveguide with a single guided mode in the operational
frequency range. Although dispersion $\omega(k)$ of a homogeneous waveguide is
a monotonic function with a positive group velocity, introducing arbitrarily
small periodicity $a$ in the direction of the waveguide axis $z$ makes the
dispersion $\omega(k)$ a periodic function with the period equal to reciprocal
lattice vector $2\pi/a$. Consequently, due to dispersion folding, i.e.
reflection of dispersion branch from the edge of the Brillouin zone, some
parts of the dispersion curves become effectively backward, as illustrated in
Fig. S2. In other words, there is a contradirectional coupling between the
mode of the strip waveguide that propagates in the negative $z$ direction and
the modes of the periodic waveguide that propagate in the positive $z$
direction. By tuning the width of the strip waveguide, we were able to achieve
the crossing of the dispersion of a strip waveguide and the intersection point
of the modes of the periodic waveguide, as shown in Fig. S2. The interaction
between these modes is further described within the framework of the coupled
mode theory, see Sec. B.

In the second design the modes of the periodic waveguide have different signs
of the group velocities. The scheme of the waveguide, its dispersion and the
field distributions at the crossing point are shown in Fig. S3. The
dispersions of the modes can be approximated as following:

$\displaystyle\omega_{1}\;\text{[rad/s]}\approx 2\pi\cdot 317.7\cdot
10^{12}-5.51\cdot 10^{7}\delta k-2.3\delta k^{2},$ (S2)
$\displaystyle\omega_{2}\;\text{[rad/s]}\approx 2\pi\cdot 317.7\cdot
10^{12}+2.26\cdot 10^{7}\delta k-13.5\delta k^{2}.$

Again, in the plane $x=90$ nm (parallel to substrate), one of the modes has
dominant $y$ component of the electric field, while the second one – $z$
component.

Figure S3: (a) Dispersion of a periodic waveguide with two intersecting modes
with group velocities of opposite sign (top view), $\lambda$ is the
wavelength, $k$ is the Bloch wavenumber. Thickness of the waveguide membrane
$h=245$ nm, period is $a=310$ nm, radius of the holes is $90$ nm, thickness of
the waveguide $w_{1}=250$ nm. (c-f) Intensities of $y$ (c,e) and $z$ (d,f)
components of electric field in the cross section of the unit cell in the
plane $x=90$ nm for the (c,d) first mode, (e,f) second mode.

## B. Coupled mode theory

###  I.1 General formulation

To analyze the relation between the geometrical parameters of the structure
and the main characteristics of the coupled modes, such as group velocity and
local and average polarization of electric field, we have applied a coupled
mode theory (CMT) Yariv and Yeh (2003); Sukhorukov  _et al._ (2008).
Throughout the text we numerate modes in such a way that indices “1” and “2”
correspond to the modes of the periodic waveguide and index “3” – to the mode
of the homogeneous waveguide. We assume, that the field distribution of a mode
of the coupled waveguides $\mathbf{E}$ can be represented as a superposition
of the modes of the isolated waveguides $\mathbf{E}_{i}$, $i=1,2,3$ in the
following way:

$\mathbf{E}(\mathbf{r},k)=\sum_{i}a_{i}(k)\mathbf{E}_{i}(\mathbf{r},k),$ (S3)

where $\mathbf{r}$ is the radius vector, $k$ is the Bloch wavenumber, and
$a_{i}(k)$ is the complex amplitude of the $i$th mode. Then, the
eigenfrequencies $\omega$ of the coupled modes as a function of $k$ can be
found from the following set of equations:

$\mathbf{M}(k,\omega)\mathbf{A}(k)=0,$ (S4)

where $\mathbf{A}=(a_{1}(k),a_{2}(k),a_{3}(k))^{T}$, and the elements of the
coupling matrix $\mathbf{M}(k)$ are defined via the overlap integrals:

$\displaystyle
M_{ij}(k,\omega)=(\omega^{2}-\omega_{j}^{2}(k))A_{ij}(k)+\omega^{2}B_{ij}(k),$
(S5) $\displaystyle
A_{ij}(k)=\langle\varepsilon_{j}(\mathbf{r})\mathbf{E}_{i}^{*}(\mathbf{r},k)\cdot\mathbf{E}_{j}(\mathbf{r},k)\rangle,$
(S6) $\displaystyle
B_{ij}(k)=\langle(\varepsilon(\mathbf{r})-\varepsilon_{j}(\mathbf{r}))\mathbf{E}_{i}^{*}(\mathbf{r},k)\cdot\mathbf{E}_{j}(\mathbf{r},k)\rangle,$
(S7)

where $i,j=1,2,3$, $\omega_{j}(k)$ defines the dispersion of the $j$th mode,
averaging is performed over the unit cell
$\langle\ldots\rangle=\int_{0}^{a}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\ldots
dxdydz$, $a$ is the period of the structure in $z$ direction,
$\varepsilon({\bf r})$ is the permittivity of the coupled waveguides and
$\varepsilon_{j}({\bf r})$ is the permittivity of the isolated waveguide
supporting $j$th mode.

When the periodic and homogeneous waveguide are brought together, the coupling
between the modes of the periodic waveguide ($M_{12}$) occurs only via the
homogeneous waveguide perturbation. In the considered structures it is much
weaker than the coupling between the modes of the periodic waveguide and the
mode of a homogeneous waveguide ($M_{13}$ or $M_{23}$). However, as it is
shown further, to obtain the desired characteristics of the guided structure
we need to have comparable strength of coupling between all of the modes. To
achieve that, we introduce a local perturbation in the periodic waveguide
(schematics are shown in Sec. C), which almost does not affect the elements
$M_{13},M_{23},M_{31},M_{32}$, while the elements $M_{12}$ and $M_{21}$ become
of the same order of magnitude as the other non-diagonal elements. In order to
correctly find the coupling between the modes “1” and “2” caused by this
perturbation, one need to modify the integral (S7) as follows Johnson  _et
al._ (2002):

$\tilde{B}_{ij}(k)=\langle(\varepsilon(\mathbf{r})-\varepsilon_{j}(\mathbf{r}))\mathbf{E}_{i\|}^{*}(\mathbf{r},k)\cdot\mathbf{E}_{j\|}(\mathbf{r},k)-(1/\varepsilon(\mathbf{r})-1/\varepsilon_{j}(\mathbf{r}))D_{i\perp}^{*}(\mathbf{r},k)\cdot
D_{j\perp}(\mathbf{r},k)\rangle,$ (S8)

where $\mathbf{E}_{\|}$ is the vector $\mathbf{E}$ with zero component normal
to the perturbation interface, and $D_{i\perp}$ is the component of the
electric displacement field vector normal to the perturbation interface. After
estimating the order of magnitude of all elements in (S5) and neglecting all
terms larger than the first order, we rewrite the matrix equation (S4) as

${\bf M}(k){\mathbf{A}(k)}=\omega(k){\mathbf{A}(k)},$ (S9)

where the elements of the newly defined coupling matrix are simplified to:

$\displaystyle M_{ii}(k)=\omega_{i}(k),$ (S10) $\displaystyle
M_{ij}(k)=-\dfrac{\omega_{j}^{2}(k)}{\omega_{i}(k)}\dfrac{\tilde{B}_{ij}(k)}{2A_{ii}(k)},\;ij=\\{12,21\\},$
$\displaystyle
M_{ij}(k)=-\dfrac{\omega_{j}^{2}(k)}{\omega_{i}(k)}\dfrac{B_{ij}(k)}{2A_{ii}(k)},\;ij=\\{13,23,31,32\\}.$

When calculating $M_{12}$ and $M_{21}$ elements, integration (S8) is performed
only over the perturbation introduced in the periodic waveguide, while for
$M_{13},M_{23},M_{31},M_{32}$ elements this perturbation can be neglected in
the integrals (S7).

To derive the requirements to the elements of the coupling matrix that allow
to achieve the dispersion with a stationary inflection point (SIP) and
circular polarization of the local electric field in the periodic waveguide,
we impose the following conditions on the eigenmode of the Eq. (S9) for a
wavenumber $k_{0}$:

$\displaystyle\dfrac{d\omega}{dk}\Bigr{|}_{k=k_{0}}=0,\;\;\dfrac{d^{2}\omega}{dk^{2}}\Bigr{|}_{k=k_{0}}=0,$
(S11) $\displaystyle\mathbf{A}(k_{0})=(1,i\alpha,\beta)^{T},$ (S12)

where $\alpha$ and $\beta$ are arbitrary real numbers (the phase of the third
mode can be chosen at will). Note that $k_{0}$ is not a wavenumber at the
crossing point of the modes “1” and “2” but rather an arbitrary wavenumber.
The Eq. (S11) specifies the existence of the SIP,
$\omega(k)\approx\omega(k_{0})+\gamma(\delta k)^{3}$, at the point $k_{0}$,
while the eigenvector (S12) ensures the desired polarization of the field in
the periodic waveguide. The phase difference between the components is fixed
as $\pm\pi/2$, while the amplitude ratio $|\alpha|$ depends on the field
distribution of the first and second modes and the choice of the
normalization. For instance, $|\alpha|$ can be chosen in such a way that the
amplitudes of the local electric field of two constituent modes at a certain
point are equal to each other, thus realizing local circularly polarized
electric field.

The relation between the elements of the coupling matrix that satisfy the
condition (S12) can be found by straight forward substitution of the required
eigenvector at a fixed value of $k_{0}$ into Eq. (S9). The conditions (S11),
on the other hand, require the knowledge of the decomposition of the matrix
elements up to second order over the wavenumber detuning $\delta k=k-k_{0}$:

$\displaystyle\omega_{i}(k)=\omega_{i}(k_{0})+v_{gi}\delta k+D_{i}\delta
k^{2},$ (S13)
$\displaystyle|M_{ij}|^{2}=|M_{ij}|^{2}_{0}+|M_{ij}|^{2}_{1}\delta
k+|M_{ij}|^{2}_{2}\delta k^{2},\;i\neq j,$

where $v_{gi}$ and $2D_{i}$ are the group velocity and group velocity
dispersion of the $i$th mode, respectively. The large amount of the variables
emerging after such decomposition significantly complicates the analysis and
makes it practically useless. However, in realistic structure one can more or
less control only the zeroth order terms, i.e. the frequency shift of the
isolated modes and the strength of coupling between the modes, but not the
first and second order terms of the decomposition. Further numerical
estimations performed for the realistic designs have shown that the
qualitative as well as quantitative requirements to the parameters of the
system can be identified with a good level of approximation even after
neglecting some of the elements. Namely, the matrix $\mathbf{M}$ can be
simplified to:

$\displaystyle M_{ii}(\delta k)=\omega_{i}(k_{0})+v_{gi}\delta k,$ (S14)
$\displaystyle
M_{ij}=-\dfrac{\omega_{j}^{2}(k_{0})}{\omega_{i}(k_{0})}\dfrac{B_{ij}(k_{0})}{2A_{ii}(k_{0})},\;i\neq
j,$

where $v_{gi}$ is the group velocity of the $i$th mode. By introducing the
frequency detuning $\widetilde{\omega}(k)=\omega(k)-\widetilde{\omega}_{0}$,
where
$\widetilde{\omega}_{0}=\widetilde{\omega}_{1}(k_{0})=\widetilde{\omega}_{2}(k_{0})$,
and normalizing the fields in such a way that $M_{ij}=M_{ji}^{*}$, we rewrite
the Eq. (S9) in the following form:

$\begin{pmatrix}v_{g1}\delta k&M_{12}&M_{13}\\\ M_{12}^{*}&v_{g2}\delta
k&M_{23}\\\ M_{13}^{*}&M_{23}^{*}&\Delta+v_{g3}\delta
k\end{pmatrix}{\mathbf{A}(k)}=\widetilde{\omega}(k){\mathbf{A}(k)},$ (S15)

where $\Delta=\widetilde{\omega}_{3}(k_{0})-\widetilde{\omega}_{0}$. The
simplified form of the CMT (S15) is mainly used in the further analysis.

###  I.2 Application of the CMT to the first design

First, we note that the quasi-degenerate modes of the periodic waveguide have
orthogonal $y$ and $z$ dominant components of the fields, while the evanescent
field of the third mode in the region, where periodic waveguide is placed, has
both $y$ and $z$ components with comparable amplitudes and with $\pm\pi/2$
phase difference. Consequently, the elements of the coupling matrix $M_{13}$
and $M_{23}$ are expected to have a fixed phase difference $\pm\pi/2$. The
absolute values of these constants also depend on the field profiles of the
modes “1” and “2”. More specifically, on the amplitude of the term with $n=-1$
in the Fourier decomposition, according to the phase matching condition. The
ratio $M_{13}/M_{23}$ extracted from the numerical simulations was equal to
-0.85i. Taking this into account, one can notice that there is a simple
solution that satisfies the conditions (S11-S12) when all three modes
intersect at the same point $k_{0}=0$, i.e. when $\Delta=0$. Then, for a zero
wavenumber detuning $\delta k$ the system has an eigenvector:

$\mathbf{A}=(1,M_{13}/M_{23},-M_{12}/M_{23})^{T},$ (S16)

with the eigenfrequency $\omega=0$, if the coupling constants $M_{12},M_{23}$
are real and $M_{13}$ is imaginary. Furthermore, the condition

$v_{g1}|M_{23}|^{2}+v_{g2}|M_{13}|^{2}+v_{g3}|M_{12}|^{2}=0.$ (S17)

ensures that the group velocity of the coupled mode vanishes, while the group
velocity dispersion vanishes automatically.

The condition (S17) can be satisfied by choosing an appropriate ratio of the
amplitudes of the $M_{13}$, $M_{23}$ and $M_{12}$ elements. The Eq. (S16)
shows that the amplitude of the second mode $\alpha$ is equal to the ratio
$M_{13}/M_{23}$, which is close to the desired value $i$ in the realistic
systems. Since all group velocities are of the same order and $v_{g1}$ and
$v_{g2}$ have the same sign, the matrix element $M_{12}$ has to be of the same
order as $M_{13}$ or $M_{23}$. This can be achieved by introducing the
deformation in the periodic waveguide that couples its modes strongly enough.
The symmetry of the deformation that provides the real valued $M_{12}$ is
discussed in the Sec. C.

The parameters and characteristics of the considered structure (Fig. 3(a) in
the main text) that were determined according to the above-mentioned recipe,
turned out to be slightly different than those expected from the CMT. Further
analysis revealed that such discrepancy is caused by the non-negligible
coupling of the modes of the periodic waveguide with another mode of the
homogeneous waveguide (lower branch of the third mode in Fig. S2). Although
this mode is well separated in frequency, the interaction with it can be
stronger than that with the backward branch due to direct coupling (not
contradirectional coupling via reciprocal lattice vector). Although in order
to find the characteristics of the system with a very fine level of accuracy
one need to include this mode into consideration, the general concept remains
the same and with high level of accuracy even when considering the coupling
only between three main modes.

Note also, that if one takes into account all first and second order terms in
the decomposition of the matrix elements over $\delta k$ in Eqs. (S13), the
group velocity dispersion of the middle branch at $k_{0}$ is modified in the
following way:

$D=\dfrac{D_{1}+D_{2}|\alpha|^{2}+D_{3}|\beta|^{2}+v_{g1}\dfrac{|M_{23}|^{2}_{1}}{|M_{23}|^{2}_{0}}+v_{g2}\dfrac{|M_{13}|^{2}_{1}}{|M_{23}|^{2}_{0}}+v_{g3}\dfrac{|M_{12}|^{2}_{1}}{|M_{23}|^{2}_{0}}}{1+|\alpha|^{2}+|\beta|^{2}}.$
(S18)

For the wavenumber $k_{0}=0$ group velocity $v_{g}$ is still equals to zero,
but according to (S18) $D$ is non-zero, which breaks the condition of the SIP
existence. However, SIP can be restored by a very slight change of the $k_{0}$
and e.g. $M_{12}$. In realistic structures the values of $D_{i}$ and
$|M_{ij}|^{2}_{1,2}$ are so small, that the required changes in $k_{0}$ and
$M_{12}$ are much less than 1%, therefore, the effect of higher order terms in
the decomposition (S13) can be neglected.

###  I.3 Application of the CMT to the second design

In the second design we assume that the second mode of the periodic waveguide
has negative velocity, while the first mode of the periodic waveguide and the
mode of the homogeneous waveguide (third mode) have positive group velocities.
Because the coupling between the modes “1” and “3” is codirectional, and
between the modes “2” and “3” it is contradirectional, typically the coupling
constant $M_{23}$ is much less than $M_{13}$. Therefore, for the second design
for simplicity we have neglected the coupling $M_{23}$ and considered the
coupling matrix in the following form:

$\mathbf{M}=\begin{pmatrix}v_{g1}k&M_{12}&M_{13}\\\ M_{12}^{*}&v_{g2}k&0\\\
M_{13}^{*}&0&v_{g3}k+\Delta\end{pmatrix}.$ (S19)

Applying the conditions (S11-S12) provides us with the relation between the
parameters of the system:

$\displaystyle
M_{12}=-i\dfrac{k_{0}\alpha^{3}(v_{1}-v_{2})v_{2}^{2}}{(v_{1}^{2}+2\alpha^{2}v_{1}v_{2}+\alpha^{2}v_{2}^{2})},$
(S20) $\displaystyle
M_{13}=\pm\dfrac{k_{0}(v_{1}-v_{2})\sqrt{-(v_{1}+\alpha^{2}v_{2})^{3}v_{3}}}{(v_{1}^{2}+2\alpha^{2}v_{1}v_{2}+\alpha^{2}v_{2}^{2})},$
$\displaystyle\Delta=\dfrac{k_{0}v_{1}[v_{1}(v_{2}-2v_{3})+v_{2}(3\alpha^{2}(v_{2}-v_{3})+v_{3})]}{(v_{1}^{2}+2\alpha^{2}v_{1}v_{2}+\alpha^{2}v_{2}^{2})},$
$\displaystyle\omega_{0}=\dfrac{k_{0}v_{1}v_{2}(v_{1}+3\alpha^{2}v_{2})}{(v_{1}^{2}+2\alpha^{2}v_{1}v_{2}+\alpha^{2}v_{2}^{2})},$
$\displaystyle\beta=-\mathop{\mathrm{sign}}\nolimits(M_{13}){\sqrt{-\dfrac{v_{1}+\alpha^{2}v_{2}}{v_{3}}}}.$

Relations (LABEL:eq:Solution_2des) show that a SIP with a desired polarization
can be realized for an arbitrary value of $k$, except the trivial case of
$k=0$. Since $M_{12}$, $M_{13}$ are proportional to $k_{0}$, while they should
be relatively small, we also aim at small values of $k_{0}$. According to
(LABEL:eq:Solution_2des) the constant $M_{12}$ in the second design should be
imaginary. This can be achieved e.g. by shifting the holes in transverse
direction, see Sec. C for details. Note, that the shift of the holes is
accompanied by the frequency shift of the modes, i.e. by appearance of the
non-zero diagonal elements $M_{11}$ and $M_{22}$. Such modification, however,
only requires to measure the frequency and the wavenumber detuning from the
new crossing point of the shifted modes “1” and “2”.

One can also understand the dispersion, shown in Fig. 2(b) in the main text
and obtained with the use of the conditions (LABEL:eq:Solution_2des), in a
different way, by considering the coupling between the modes in two steps. At
the first step one we consider only the coupling between the modes of the
periodic waveguide. Introducing perturbation to the waveguide couples two
modes with the phase difference dictated by the coupling constant $M_{12}$.
The dispersion of the modes of the unperturbed and perturbed periodic
waveguide are shown in Fig. S4(a) with thin gray and thick blue curves,
respectively. As was mentioned earlier, due to the non-zero diagonal terms of
the coupling matrix along with the formation of the anticrossing there is a
slight shift of the dispersion. By choosing the imaginary constant $M_{12}$,
upper and lower dispersion curves acquire elliptical polarization with either
$\pi/2$ or $-\pi/2$ phase difference between the constituent modes. At the
second step, we choose the frequency shift of the third mode $\Delta\approx
M_{12}$, as shown with the green line in Fig. S4(b). Then, the third mode
couples to the upper branch of the perturbed periodic waveguide at the point
with negative group velocity. Consequently, by choosing the appropriate value
of the $M_{13}$ constant we obtain the SIP with the polarization defined
mainly by the constant $M_{12}$ at the first step, see purple curves in Fig.
S4(c).

Figure S4: (a) Dispersion of the modes of a periodic waveguide with
symmetrically placed holes (gray curves) and shifted holes (blue curves). (b)
Blue curves – the same as in panel (a), green curve shows the dispersion of
the mode of an isolated homogeneous waveguide. (c) Gray curves – the same as
blue and green curves in panel (b), purple curves show the dispersion of the
coupled periodic and homogenous waveguides. Black dot indicates a SIP.

## C. Relation between the geometrical parameters of the nanostructure and
coupling constants in the CMT

Although all variables in the CMT depend on the geometrical parameters of the
structure, the main parameters that can be controllably tuned by the
geometrical adjustment are the frequency shifts of the modes of isolated
waveguides and coupling constants at a given wavenumber.

Since the operational frequency is always close to the crossing point of the
dispersion curves of two modes of the periodic waveguide, we measure the
frequency shift $\Delta$ of the third mode from this frequency. This is done
straight forwardly by changing the width of the strip waveguide $w$. Decrease
(increase) of the waveguide thickness results in the increase (decrease) of
the frequency of the waveguide mode for a given wavenumber. For small values
of the detuning, relevant to the current designs, $\Delta$ is linearly
proportional to the $\delta w$, thus allowing for simple tuning of the
diagonal element $M_{33}$ of the coupling matrix.

The elements $M_{13}$ and $M_{23}$ describe the coupling between the modes of
the periodic waveguide and the mode of the homogeneous waveguide. The phase
difference between these elements is fixed as $\pi/2$, while their absolute
values are determined by two factors. First, the less is the separation
between the waveguides, the larger are the coupling constants. Second, since
in some cases the coupling between the modes is contradirectional, the
efficiency of coupling also depends on the $n=-1$ component in the Fourier
decomposition of the field profiles of the periodic waveguide modes.
Typically, a shift of the operational wavenumber closer to the edge of the
Brillouin zone increases the amplitude of the minus first Fourier component of
the mode and thus increases the coupling strength.

Figure S5: Schematic illustration of the two types of asymmetric perturbation
in the unit cell of a periodic waveguide. Dashed curve shows the boundary of a
hole in the unperturbed waveguide, solid curve shows the modified position of
the hole. (a) Perturbation is symmetric with respect to rotation around the
$x$ axis by 180 degrees. (b) Symmetry with respect to the $z=0$ plane is
preserved.

The element $M_{12}$ can be controlled by introducing a local perturbation in
the periodic waveguide. The value of $M_{12}$ is proportional to the overlap
integral $\int_{\text{unit
cell}}\Delta\varepsilon(\mathbf{r})C(\mathbf{r})d\mathbf{r}$, (S8), where
$\Delta\varepsilon$ describes the change of the permittivity as compared to
the unperturbed waveguide. Taking into account that $y=0$ is magnetic
(electric) symmetry plane of the first (second) mode, and that rotation around
the $x$ axis transforms the mode into its time reversal partner with complex
conjugate field, one gets the following properties of the function $C$:
$C(-y)=-C(y)$, $C(-z)=-C^{*}(z)$. Therefore, non-zero coupling is achieved by
the perturbation that breaks the symmetry with respect to the $y=0$ plane. In
general, the coupling constant $M_{12}$ is complex-valued. However, if the
symmetry of the perturbation $\Delta\varepsilon$ with respect to $z=0$ plane
is preserved, the element $M_{12}$ becomes purely imaginary. And if instead
the perturbation has a rotation symmetry around $x$ axis, the coupling is
purely real. This is illustrated in Fig. S5.

Besides the coupling between the modes, such perturbation also induces the
spectral shift of the modes, i.e. non-zero diagonal elements $M_{11}$ and
$M_{22}$. These can be partially compensated if along with the positive change
of the permittivity in some region (dielectric instead of vacuum,
$\Delta\varepsilon>0$), there is also negative change in another region
(vacuum instead of dielectric $\Delta\varepsilon<0$). For instance, in Fig.
S5(a), the domains $D_{1}$ provide negative contribution to the overlap
integral, while $D_{2}$ – positive contribution. Note, that such compensation
is not ideal, because of the slightly different field distributions in the
domains $D_{1}$ and $D_{2}$ in the unperturbed waveguide.

## References

  * Savelev  _et al._ (2019) R. S. Savelev, D. F. Kornovan, V. V. Yaroshenko, and M. I. Petrov, Journal of Applied Physics 125, 123104 (2019).
  * Savelev and Gorlach (2020) R. S. Savelev and M. A. Gorlach, Phys. Rev. B 102, 161112 (2020).
  * Yariv and Yeh (2003) A. Yariv and P. Yeh,  _Optical Waves in Crystals: Propagation and Control of Laser Radiation_ (John Wiley and Sons, New York, 2003).
  * Sukhorukov  _et al._ (2008) A. A. Sukhorukov, A. V. Lavrinenko, D. N. Chigrin, D. E. Pelinovsky, and Y. S. Kivshar, J. Opt. Soc. Am. B 25, C65 (2008).
  * Johnson  _et al._ (2002) S. G. Johnson, M. Ibanescu, M. A. Skorobogatiy, O. Weisberg, J. D. Joannopoulos, and Y. Fink, Phys. Rev. E 65, 066611 (2002).


# Nonadiabatic simulations of photoisomerization and dissociation in ethylene
using ab initio classical trajectories

K. Miyazaki  N. Ananth<EMAIL_ADDRESS>Department of Chemistry and
Chemical Biology, Cornell University, Ithaca, New York, 14853, United States

###### Abstract

We simulate the nonadiabatic dynamics of photo-induced isomerization and
dissociation in ethylene using ab initio classical trajectories in an extended
phase space of nuclear and electronic variables. This is achieved by employing
the Linearized Semiclassical Initial Value Representation (LSC-IVR) method for
nonadiabatic dynamics where discrete electronic states are mapped to
continuous classical variables using either the Meyer-Miller Stock-Thoss
representation or a more recently introduced spin mapping approach. Trajectory
initial conditions are sampled by constraining electronic state variables to a
single initial excited state, and by drawing nuclear phase space
configurations from a Wigner distribution at finite temperature. An ensemble
of classical ab initio trajectories are then generated to compute thermal
population correlation functions and to analyze the mechanisms of
isomerization and dissociation. Our results serve as a demonstration that this
parameter-free semiclassical approach is computationally efficient and
accurate, identifying mechanistic pathways in agreement with previous
theoretical studies, and also uncovering dissociation pathways observed
experimentally.

††preprint: AIP/123-QED

##  I Introduction

Nonadiabatic effects that result from the coupling of nuclear motion to
electronic transitions are crucial to understanding the mechanism of various
interesting photophysical and photochemical phenomena including internal
conversion, intersystem crossing, and photodissociation. [1] Unfortunately,
the high dimensionality of complex molecular systems makes real-time dynamic
simulations on a pre-computed potential energy surface challenging and
motivates the development of on-the-fly dynamic methods.

The search for computationally efficient and accurate on-the-fly quantum
dynamic simulations remain an outstanding challenge, particularly in systems
where nonadiabatic effects play a central role. Although several methods have
been developed for the simulation of nonadiabatic processes, only a handful
lend themselves to on-the-fly implementations that require not just routine
force calculations but also evaluation of the nonadiabatic coupling vector at
each time step. Multiconfigurational Time-dependent Hartree (MCTDH) is,
arguably, the most exact of these techniques relying on a weighted sum of
Gaussian nuclear functions centered on a grid while solving the Schrödinger
equation exactly for the electronic degrees of freedom. [2] Variational
Gaussian methods are extensions of MCTDH where the nuclear wavefunction is
replaced by a Gaussian wavepacket. [3] In multiple spawning and related
methods, [4, 5] the number of nuclear basis functions increase through
‘spawning’ events that occur when electronic states become near-degenerate,
mimicking the quantum dynamical bifurcation of wavepackets. Formally, multiple
spawning can converge to exact quantum dynamics if a sufficiently large number
of wavepackets are included, and the transition matrix elements are evaluated
exactly at each step. [6] However, in practice, it is necessary to limit the
number of wavepackets spawned and approximate matrix elements as is done in
approaches like ab initio multiple spawning (AIMS). [3, 7, 8]

Mixed quantum-classical methods like surface hopping and Ehrenfest dynamics
further approximate nuclear motion by treating it classically. In surface
hopping [9, 10] nuclear wavepackets with finite widths in the nuclear
coordinates are replaced by independent classical trajectories. Each
trajectory is allowed to ‘hop’ between surfaces according to a stochastic
algorithm based on the amplitudes of each electronic state. Ehrenfest is a
mean-field dynamics, where the potential energy for classical nuclear
propagation is obtained by the weighted average of electronic state energies
obtained at an instantaneous nuclear geometry. [11] While the non-interacting
classical trajectories simplify the dynamics, both surface hopping and
Ehrenfest are known to suffer from strong coherence and slow decoherence,
necessitating additional ad hoc modifications. [12, 13, 14, 15, 16, 17, 18]

The methods discussed thus far focus on approximately time evolving
wavepackets that can then be used to compute real-time correlation functions
necessary to compare with experimental observables. Semiclassical (SC) and
path-integral based methods [19, 20, 21] take an alternate approach, directly
approximating the real-time quantum correlation function. As reviewed
recently, there is a hierarchy of SC methods that include quantum effects and
more specifically nonadiabatic effects to differing extents. [22] However,
only the more classical-limit SC methods are sufficiently efficient to allow
for on-the-fly dynamic simulations. In particular, the Linearized
Semiclassical Initial Value Representation (LSC-IVR) approach [23, 24] is a
classical-limit method that approximates the quantum correlation functions
between two operators, $\hat{A}$ and $\hat{B}$, as a single phase space
integral over a product of Wigner functions evaluated at time zero and
classically time-evolved phase space configurations at later time $t$
respectively. LSC-IVR is exact at time zero and for harmonic potentials, and
offers a simple implementation. Recently, the symmetrical quasiclassical (SQC)
method has been introduced [25, 26, 27] where the Wigner functions in LSC-IVR
are replaced by ‘window’ functions that have been successfully used to study a
range of model chemical systems, [26, 28, 27, 29] ab initio dynamics using the
so-called quasi-diabatic electronic states [30, 31] as well as using adiabatic
electronic states with an approximate integration scheme. [29] Unfortunately,
unlike LSC-IVR, SQC requires the choice of an ad hoc window function that can
significantly affect simulation accuracy. [32, 30, 31]

In this paper, we implement a first on-the-fly nonadiabatic dynamic simulation
with LSC-IVR using ab initio classical trajectories in an extended phase space
of nuclear and adiabatic electronic state variables. The classical Hamiltonian
employed by these semiclassical methods is obtained by mapping discrete
electronic state variables to continuous Cartesian variables. We explore two
mapping protocols: the well-established Meyer-Miller-Stock-Thoss (MMST)
approach [33, 34, 35, 36] and the more recently introduced spin mapping that
shows great promise in model system studies. [37, 38] To enable on-the-fly
simulations, we use an adiabatic electronic state representation [35, 39, 40]
in conjunction with the ‘kinematic’ momentum integration scheme. [40]
Trajectory initial conditions for the nuclei are sampled from an initial
Wigner distribution, while electronic state variables are sampled to ensure
that only a single excited electronic state is populated at time zero. We
calculate the real-time population correlation function by generating
classical trajectories from the mapping Hamiltonian, with the potential
gradients, energies, and the nonadiabatic coupling vector calculated at each
time step from a CASSCF(2o,2e) electronic structure calculation. We then
analyze the resulting ensemble of trajectories to identify the mechanisms of
photoisomerization and dissociation pathways in ethylene, and compare our
results against previous theoretical and experimental studies. We conclude
with a detailed discussion highlighting the significant advantages of the LSC-
IVR approach for on-the-fly nonadiabatic simulations while also outlining
outstanding challenges.

##  II Theory

###  II.1 LSC-IVR approximation to quantum correlation functions

In the path integral representation of quantum mechanics, the real-time
quantum correlation function,

$\displaystyle
C_{AB}(t)=\Tr\left[\hat{A}e^{i\hat{H}t/\hbar}\hat{B}e^{-i\hat{H}t/\hbar}\right],$
(1)

is expressed as a double sum over all possible forward and backward paths in
coordinate space with an overall phase corresponding to the difference in
action between the forward and backward paths. Truncating the difference in
action to first order yields the LSC-IVR approximation [23, 41] to the
correlation function in Eq. (1),

$\displaystyle C_{AB}^{LSC}(t)=\cfrac{1}{(2\pi\hbar)^{N}}\int
d\mathbf{X}_{0}\int
d\mathbf{P}_{0}\,A_{\mathcal{W}}(\mathbf{X}_{0},\mathbf{P}_{0})B_{\mathcal{W}}(\mathbf{X}_{t},\mathbf{P}_{t}).$
(2)

In Eq. (2), $O_{\mathcal{W}}(\mathbf{X},\mathbf{P})$ represents the Wigner
transform of $\hat{O}$ defined as

$\displaystyle O_{\mathcal{W}}(\mathbf{X},\mathbf{P})=\int
d\bm{\Delta}\matrixelement{\mathbf{X}+\frac{\bm{\Delta}}{2}}{\hat{O}}{\mathbf{X}-\frac{\bm{\Delta}}{2}}e^{-i\mathbf{P}\bm{\Delta}/\hbar},$
(3)

where $N$ is the number of system degrees of freedom, and
$(\mathbf{X},\mathbf{P})$ are phase space vectors. An ensemble of trajectories
are generated by sampling initial conditions from
$A_{\mathcal{W}}(\mathbf{X}_{0},\mathbf{P}_{0})$ and propagated for time $t$
according to classical equations of motion generated by the Hamiltonian,
$H(\mathbf{X},\mathbf{P})$. The function $B_{\mathcal{W}}$ is then evaluated
at the time-evolved phase space variables $(\mathbf{X}_{t},\mathbf{P}_{t})$.
Calculating the LSC-IVR correlation function in Eq. 2 is generally efficient,
incorporating quantum effects like nuclear tunneling and zero-point energy at
a computational cost similar to a classical simulation.

###  II.2 Thermal Correlation Functions and Nonadiabatic Dynamics

The LSC-IVR expression for a real-time thermal correlation function requires
the Wigner transform of the density operator,
$\hat{A}\equiv\hat{\rho}=e^{-\beta\hat{H}}$ (Eq. 3), where $\beta=1/k_{B}T$.
For a multi-state system, we assume the initial density operator is separable,
$\hat{\rho}=\hat{\rho}_{e}\hat{\rho}_{n}$, where $\hat{\rho}_{e}$ and
$\hat{\rho}_{n}$ are the electronic and nuclear density operators,
respectively, and
$\Tr_{e}\left[\hat{\rho}_{e}\right]=\Tr_{n}\left[\hat{\rho}_{n}\right]=1$.
Under this assumption, the Wigner transform of $\hat{A}$ can be written as the
product of separate transforms:
$A_{\mathcal{W}}(\mathbf{R},\mathbf{P},\mathbf{x},\mathbf{p})=\left[\hat{\rho}_{e}\right]_{\mathcal{W}}(\mathbf{x},\mathbf{p})\left[\hat{\rho}_{n}\right]_{\mathcal{W}}(\mathbf{R},\mathbf{P})$
where $\left(\mathbf{R},\,\mathbf{P}\right)$ are nuclear phase space vectors
and the $\left(\mathbf{x},\,\mathbf{p}\right)$ are the electronic state phase
space variables. Assuming a thermal distribution of harmonic normal modes, it
is possible to express the Wigner transform of nuclear density as

$\displaystyle\left[\hat{\rho}_{n}\right]_{\mathcal{W}}$
$\displaystyle=\left[e^{-\beta\hat{H}_{n}}\right]_{\mathcal{W}}$
$\displaystyle=\prod_{i=1}^{N-F}\cfrac{1}{2\pi}\,\exp\left[-\tanh\left(\cfrac{\beta\omega_{i}}{2}\right)\left(\cfrac{1}{\mu_{i}\omega_{i}}P_{i}^{2}+\mu_{i}\omega_{i}R_{i}^{2}\right)\right],$
(4)

where $\mu_{i}$ and $\omega_{i}$ are the reduced mass and frequency of the
$i$-th vibrational mode, $N$ is the total number of degrees of freedom in the
system (electronic and nuclear) and $F$ is the number of electronic states.

For photo-initiated processes, the electronic density is defined as the
projection operator onto a single initially occupied $i^{\text{th}}$
electronic state,

$\displaystyle\left[\hat{\rho}_{e}\right]_{W}=\left[\ket{i}\bra{i}\right]_{W}.$
(5)

In the SC framework, the discrete state-space electronic density matrix and
the corresponding multi-state system Hamiltonian are treated by mapping them
to a continuous Cartesian variable representation, and we describe two such
schemes below.

####  II.2.1 MMST mapping:

The MMST approach effectively maps the occupation of an electronic state to a
single excitation quantum of a corresponding harmonic oscillator, [33]

$\displaystyle\ket{n}\bra{m}\rightarrow\hat{a}^{\dagger}_{n}\hat{a}_{m}$ (6)

where $F$ is the total number of electronic states and the creation and
annihilation operators of $k$th harmonic oscillator are
$\hat{a}^{\dagger}_{k}=(\hat{x}_{k}-i\hat{p}_{k})/\sqrt{2}$ and
$\hat{a}_{k}=(\hat{x}_{k}+i\hat{p}_{k})/\sqrt{2}$ in terms of electronic phase
space variables. Using this mapping, the multi-state Hamiltonian expressed in
the adiabatic electronic state representation can be written as, [42, 43, 40]

$\displaystyle H$
$\displaystyle=\frac{(\mathbf{P}+\Delta\mathbf{P})^{2}}{2\mu}+\sum_{n}^{F}\frac{1}{2}\left(p_{n}^{2}+x_{n}^{2}-\gamma\right)E_{n}(\mathbf{R})$
$\displaystyle=\frac{\mathbf{P}^{2}_{\text{kin}}}{2\mu}+\sum_{n}^{F}\frac{1}{2}\left(p_{n}^{2}+x_{n}^{2}-\gamma\right)E_{n}(\mathbf{R}).$
(7)

We note that the MMST mapping is exact when $\gamma=1$, and this is the value
we use in this manuscript although some semiclassical simulations treat
$\gamma$ as a zero-point energy (ZPE) parameter that can be modified to
increase numerical stability. In Eq. II.2.1, $E_{n}(\mathbf{R})$ is the
adiabatic energy of the $n$-th electronic state, the kinematic momentum
$\mathbf{P}_{\text{kin}}=\mathbf{P}+\Delta\mathbf{P}$, with

$\displaystyle\Delta\mathbf{P}=\sum_{n\neq
m}^{F}x_{n}p_{m}\mathbf{d}_{nm}(\mathbf{R}),$ (8)

where
$\mathbf{d}_{nm}(\mathbf{R})=\matrixelement{\phi_{n}}{\frac{\partial}{\partial\mathbf{R}}}{\phi_{m}}$
is the nonadiabatic coupling vector between the electronic states
$\ket{\phi_{n}}$ and $\ket{\phi_{m}}$. [40] It is further possible to
construct an equivalent, exact, symmetrized form of the MMST Hamiltonian,

$\displaystyle H$
$\displaystyle=\cfrac{\mathbf{P}_{\text{kin}}^{2}}{2\mu}+\cfrac{1}{F}\sum_{n}^{F}E_{n}(\mathbf{R})$
$\displaystyle\,\,\,\,\,\,\,\,+\cfrac{1}{F}\sum_{n,m}^{F}\cfrac{1}{4}\,\left(p_{n}^{2}+x_{n}^{2}-p_{m}^{2}-x_{m}^{2}\right)\left(E_{n}(\mathbf{R})-E_{m}(\mathbf{R})\right)$
$\displaystyle=\cfrac{\mathbf{P}_{\text{kin}}^{2}}{2\mu}+V_{\text{eff}},$ (9)

such that the equations of motion are independent of the ZPE parameter. [26,
40]

In order to evaluate the electronic density matrix as defined in Eq. 5, we
consider the Wigner transform of the projection operator. In the MMST
framework, the phase space expression obtained through Wigner transformation
differs depending on the quantum mechanical definition, [41, 44, 45] and we
consider both forms in this work. Defining the projection operator in the
singly excited oscillator (SEO) basis yields the so-called Wigner population
estimator, [46]

$\displaystyle P^{i}_{\mathcal{W}}$
$\displaystyle=\left[\ket{i}\bra{i}\right]_{\mathcal{W}}^{\text{SEO}}$
$\displaystyle=2^{F+1}\left(x_{i}^{2}+p_{i}^{2}-\cfrac{1}{2}\right)\exp\left[-\sum_{j}^{F}\left(x_{j}^{2}+p_{j}^{2}\right)\right].$
(10)

Expressing the projection operators using the mapping variable operators
yields the SC population estimator, [44]

$\displaystyle
P^{i}_{\text{SC}}=\left[\cfrac{1}{2}\left(\hat{x}_{i}^{2}+\hat{p}_{i}^{2}-1\right)\right]_{\mathcal{W}}=\cfrac{1}{2}\left(x_{i}^{2}+p_{i}^{2}-1\right).$
(11)

####  II.2.2 Spin mapping

As an alternative to MMST mapping, various schemes that map electronic states
to spin variables have been previously proposed; [47, 34, 48] here we employ a
recently introduced scheme that appears to out-perform MMST mapping in model
system studies. [37, 38, 49]

The spin mapping (SM) approach maps an $F$-level electronic system to an
$(F-1)/2$ spin system, [37] and the corresponding Hamiltonian is obtained by
recognizing that the $F$-level system Hamiltonian can be exactly expressed as
a linear combination of the $F^{2}-1$ spin angular momentum matrices that
define a spin $(F-1)/2$ system and the identity operator,

$\displaystyle\hat{H}=H_{0}\hat{\mathds{1}}+\sum_{i=1}^{F^{2}-1}H_{i}\hat{S}_{i}$
(12)

where $H_{0}$ and the $\\{H_{i}\\}$ are the expansion coefficients, and the
spin angular momentum matrices, $\\{\hat{S}_{i}\\}_{i=1}^{F^{2}-1}$, are
traceless and orthogonal. The Hamiltonian in Eq. (12) can be written in terms
of continuous phase-space Cartesian operators using the Stratonovich-Weyl
transforms. [50, 51, 52] Interestingly, after symmetrization, in the adiabatic
electronic state representation the resulting SM Hamiltonian is identical to
the symmetrized MMST Hamiltonian in Eq. (II.2.1). In this manuscript, we focus
on the $W$-representation of the Stratonovich-Weyl transform, which results in
a correlation function analogous to that of LSC-IVR (Eq. 2), and that recent
work suggests is the most successful choice of representation. [37, 38]

As before, we define the initial electronic density matrix as a projection
onto a single discrete electronic states. Since the specific form of the
projection operator depends on the number of electronic states, we consider
the case where $F=3$ in line with our study of the photoisomerization of
ethylene. Expressing the projection operator in terms of spin matrices,

$\displaystyle\ket{j}\bra{j}=C_{0}\hat{\mathds{1}}+\sum_{i}^{F^{2}-1}C_{i}\hat{S}_{i},$
(13)

and using the $W$-representation leads to expressions for the population
estimators, [38]

$\displaystyle P_{\text{SM}}^{1}=\left[\ket{1}\bra{1}\right]_{\text{SM}}$
$\displaystyle=\cfrac{1}{3}+\cfrac{1}{6}\left(2r_{1}^{2}-r_{2}^{2}-r_{3}^{2}\right)$
(14a) $\displaystyle
P_{\text{SM}}^{2}=\left[\ket{2}\bra{2}\right]_{\text{SM}}$
$\displaystyle=\cfrac{1}{3}-\cfrac{1}{6}\left(r_{1}^{2}-2r_{2}^{2}+r_{3}^{2}\right)$
(14b) $\displaystyle
P_{\text{SM}}^{3}=\left[\ket{3}\bra{3}\right]_{\text{SM}}$
$\displaystyle=\cfrac{1}{3}-\cfrac{1}{6}\left(r_{1}^{2}+r_{2}^{2}-2r_{3}^{2}\right)$
(14c)

where $r_{i}^{2}=x_{i}^{2}+p_{i}^{2}$ and
$\sum_{j}^{F}\left[\ket{j}\bra{j}\right]_{\text{SM}}=1$.

Table 1:  The radii of the sampling functions for the three electronic population estimators ($\mathcal{W}$, SC, and SM) for a system with three coupled electronic states. $M$ (Eq. 16) |  Sampling radius, $r_{k}$
---|---
|  Occupied |  Unoccupied
$\mathcal{W}$ (Eq. II.2.1) |  1.55892111The larger of the two roots is chosen following Ref. 46. |  1/2
SC (Eq. 11) |  3 |  1
SM (Eq. 14) |  8/3 |  2/3

###  II.3 Simulation details

Excited state dynamics in the LSC-IVR framework are obtained from the
population correlation function,

$\displaystyle C_{P_{j}}(t)$ $\displaystyle=\frac{1}{(2\pi\hbar)^{N}}\int
d{\bf R}_{0}\int d{\bf P}_{0}\int d{\bf x}_{0}\int d{\bf p}_{0}$
$\displaystyle\times\left[\hat{\rho}_{n}\right]_{\mathcal{W}}({\bf R}_{0},{\bf
P}_{0})\left[\hat{\rho}_{e}\right]_{M}({\bf x}_{0},{\bf
p}_{0})P^{j}_{M}\left({\bf x}_{t},{\bf p}_{t}\right),$ (15)

where the Wigner transform of the nuclear density matrix,
$\left[\hat{\rho}_{n}\right]_{\mathcal{W}}$, is defined in Eq. (II.2) and used
to sample the initial nuclear phase space variables at temperature $T=300$K,
with the frequencies and the reduced masses calculated at the equilibrium
geometry of the $S_{0}$ state. For a system initially in a single excited
state, the electronic density matrix in Eq. (15) can be expressed in terms of
electronic population estimators,

$\displaystyle\left[\hat{\rho}_{e}\right]_{M}=\delta\left(P_{M}^{1}-1\right)\prod_{j\neq
1}^{2}\delta\left(P_{M}^{j}\right),$ (16)

where the subscript $M=\mathcal{W},\,\text{SC},\,\text{SM}$ indicates the
specific choice of mapping framework/estimator, and $P_{M}^{i}$ is the
population estimator for the $i$-th electronic state (with $i=0,1,2$ for the
$S_{0}$, $S_{1}$, and $S_{2}$ electronic states of ethylene, respectively).
Initial values for the electronic phase space variables are sampled using the
focusing approximation [53, 46] such that the initial electronic state
population is exactly 1 for the occupied $S_{1}$ state and 0 for the
unoccupied states ($S_{0}$ and $S_{2}$). This is achieved by sampling initial
electronic mapping variables for the $i$-th electronic state from a circle
with radius $x_{i}^{2}+p_{i}^{2}=r_{i}^{2}$; The radii for occupied and
unoccupied states in each implementation are specified in Table 1, and are
obtained by solving the corresponding equations for the population estimators
provided in Eq. II.2.1 for the Wigner estimator in the MMST framework, in Eq.
11 for the SC estimators in the MMST framework, and in Eq. 14 for
$W$-representation in the SM framework. Finally, in Eq. (15), the electronic
population estimator at time $t$, $P^{i}_{M}({\bf x}_{t},{\bf p}_{t})$, is
evaluated using the time-evolved electronic positions and momenta obtained by
propagating trajectory initial conditions under the classical analog
symmetrized mapping Hamiltonian defined in Eq. (II.2.1) for all three
implementations.

Classical equations of motion are integrated using the $4^{\text{th}}$ order
Adams-Bashforth-Moulton predictor-corrector integrator with a time step of 1
a.u. time ($\sim 1/40$ fs) for a total time of $400$ fs. The necessary
classical forces at each time step are obtained from on-the-fly calculations
of the adiabatic energies of the three electronic states, their gradients, and
the nonadiabatic coupling vectors using Pople’s $6$-$31$+G∗ basis set and
CASSCF(2o2e) with three states included in state averaging in the electronic
structure package GAMESS.[54]

##  III Result and discussion

###  III.1 Photoisomerization

We discuss the results of our on-the-fly ab initio LSC study of ethylene
photoisomerization and dissociation in context with previous theoretical
efforts using AIMS, [55, 56, 57, 58] surface hopping, [59] and SQC [30, 31] as
well as experimental studies. [60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71]
As with detailed earlier theoretical simulations of this system, [55] we
include the ground state ($S_{0}$) and two excited states
($\pi\rightarrow\pi^{*}=S_{1}$ and $\pi^{*2}=S_{2}$), and exclude the low-
lying Rydberg states that are expected to be present but appear to play no
significant role in the quenching process. [56, 57, 72] In FIG. 1, we plot the
population correlation function for the three electronic states obtained using
the Wigner population estimator in the MMST mapping framework. We find that
the initial photo-excited $S_{1}$ state decays to about 50% of its original
population on a 60-80 fs timescale, in good agreement with previous theory
[59, 30, 31] and experiment. [70] This is accompanied by a rapid but small
population transfer to the higher lying $S_{2}$ state that, in turn, decays
back to the ground state in about $100$ fs.

Analysis of the trajectory ensemble employed in our LSC simulation yields
insights into the coupled nuclear motions that drive electronic state
transitions. In FIG. 2, we show the timeline on which key ethylene molecular
structures appear as trajectories enter regions with a near-degeneracy between
two electronic states, defined here as $E_{S1}-E_{S0}<0.2$ eV and
$E_{S2}-E_{S1}<0.2$ eV. We find that there are four primary structures, as
suggested by previous work. [55, 73] These structures are shown in FIG 3:
twisted, pyramidalized (pyramidalization of one or both $\text{CH}_{2}$
groups), H-migration, and ethylidene. The twisted geometry dominates at early
times ($\sim 10$ fs) weakening the $\pi$ electronic structure and enabling
population transfer from the $S_{1}$ to the doubly excited $S_{2}$ state as
seen in FIG. 1. The weakening of the $\pi$ bond also facilitates the formation
of the remaining three structures that occur at near-degeneracies between the
$S_{1}$ and $S_{0}$ states. Pyramidalization and H-migration appear at
relatively early times ($<80$ fs) while the formation of ethylidene occurs at
later times and persists, in keeping with experiments that see signatures of
this structure up to $600$ fs. [58]

Figure 1: The population of the three electronic states obtained from LSC-IVR
simulation that employ the MMST mapping with Wigner estimator are shown. The
initial photo-excited $S_{1}$ state population is shown in red, the doubly
excited $S_{2}$ state population is shown in blue, and the ground $S_{0}$
state population is shown in green. These results were obtained by averaging
the number of the isomerization trajectories shown in Table 2. The error bars
in the plot are obtained from the averages of $3$ subsets of trajectories,
each containing $48$ trajectories. Figure 2: Characterization of molecular
structures observed in $300$ trajectories from the conventional LSC-IVR
simulation when the $S_{2}$-$S_{1}$ and $S_{1}$-$S_{0}$ energy gaps are small
($<0.2$eV). The observed structures can be classified into the four structures
associated with the ethylene conical intersections (FIG. 3). The
$S_{2}$-$S_{1}$ small energy gaps are characterized mostly by twisted
structure and the $S_{1}$-$S_{0}$ gap by H-migration, pyramidalized or
ethylidene structure.

(a) (b)

(c) (d)

Figure 3: The representative geometries of conical intersections encountered
in the quenching of ethylene excited state through isomerization: (a) twisted,
(b) pyramidalized, (c) H-migration, and (d) ethylidene.

###  III.2 Photodissociation pathways

In addition to isomerization, we observe a significant number of trajectories
that describe photodissociation of ethylene upon excitation to $S_{1}$ as
documented in Table 2. This is in keeping with both previous experimental work
[60, 61, 62, 63, 67, 68, 70] and some theoretical simulations. [74, 75, 58]

After photoisomerization, we see the elimination of molecular H2 as the most
significant channel for photodissociation. In analyzing this subset of
trajectories, we can identify the structures that are the major precursors of
$\text{H}_{2}$ molecule production: about 39% of the trajectories produce H2
and acetylene (HCCH) from the ethylidene structure shown in FIG. 4a, while the
majority (52%) produce H2 and vinylidene ($\text{H}_{2}\text{CC:}$) from the
pyramidalized geometry shown in FIG. 4b. These findings reproduce the two
channels for H2 molecule elimination identified by experiments and the
intermediate structures we report have been previously characterized by
transition state calculations using electronic structure. [67, 68]

(a) (b)

Figure 4: The structures at transition states leading to $\text{H}_{2}$
elimination: (a) ethylidene resulting in a $\text{H}_{2}$ molecule and
acetylene (HCCH) and (b) pyramidalized structure resulting in a $\text{H}_{2}$
molecule and vinylidene ($\text{H}_{2}\text{CC:}$).

A smaller but statistically significant number of our trajectories as reported
in Table 2 also describe the elimination of 2 H atoms, the first leading to
formation of the vinyl radical ($\text{H}_{2}\text{CCH}$), and the subsequent
H atom elimination yielding acetylene. We note that while this mechanism is in
keeping with experimental studies, [67, 68] a very low but persistent yield
(2%) of the vinyl radical has been observed suggesting that it is possible to
stop with the elimination of a single H-atom. We find $\sim 5\%$ of
trajectories in our simulation correspond to the single H atom dissociation
channel, which is in a reasonable agreement with this experimental
observation. Although experimental measurements of the H2:H branching ratio at
$157$ nm excitation differ slightly, $0.44:0.56$ in Ref. 63 and $0.3:0.7$ in
Ref. 67, both consistently find higher yields of H atoms than $\text{H}_{2}$
molecule. The LSC-IVR simulations do not reproduce this ratio, yielding H2:H
$=0.66:0.34$, an error we attribute the level of electronic structure theory
employed rather than the dynamics as discussed extensively in the context of
C-C bond cleavage.

Finally, we note that in Table 2, we show that 20% of our trajectories result
in unphysical C-C bond cleavage for which, to our knowledge, there is no
experimental evidence other than in pump-probe or high energy ionization
experiments, [58, 71] where the cleavage is explicitly targeted. It is
important to identify the source of this error as due to either the nature of
the dynamics employed or the underlying electronic structure. To unravel this,
we start by noting that the initial excitation energy in our simulation
obtained with SA2-CASSCF(2o,2e) with $6$-$31+$G∗ basis, $\langle\Delta
E_{S_{1}S_{0}}(t=0)\rangle$, is $9.77$ eV as reported in Table. 2, a number
significantly above the C-C bond dissociation energy of $7.7$eV at this level
of theory. Notably, our $\langle\Delta E_{S_{1}S_{0}}(t=0)\rangle$ is also
significantly above the experimentally quantified Franck-Condon excitation
energy of $7.6$ eV.[76] We find that using a more extensive basis (aug-cc-
pVDZ) reduces the calculated excitation energy to $8.86$eV, and further using
second-order perturbative energy correction (XMCQDP) yields $7.63$eV, a value
in close agreement with experiment. Unfortunately, at present, we cannot
calculate the nonadiabatic coupling matrix element (NACME) in GAMESS with the
perturbation correction, so in order to test the dynamics we run an additional
50 trajectories using the aug-cc-pVDZ basis set. We find that the number of
trajectories that exhibit C-C bond cleavage drops significantly from 20% to
only 10%. This suggests that the underlying source of error leading to
unphysical trajectories in this case can indeed be attributed to the level of
electronic structure theory. The dependence of the nature of trajectories on
the initial energy can also be seen by simply removing high energy
trajectories: when we do not include trajectories with $\Delta
E_{S_{1}S_{0}}(t=0)>10$ eV, the fraction of trajectories exhibiting C-C
cleavage drops by 10%, while that of isomerization increases by 10% and other
categories remain almost unchanged.

###  III.3 Comparing the three different variants of LSC-IVR

In discussing the results of our ab initio study of photoisomerization and
dissociation in ethylene, we confined ourselves to interpreting the results
from the LSC simulation with MMST mapping and Wigner estimators. We now
motivate this choice and further provide a detailed discussion of the three
variants of LSC explored here, paying particular attention to energy
conservation, observed reaction channels, and electronic population dynamics.

Figure 5: The absolute value of jumps in the total energy between two
consecutive time steps in 300 simulated trajectories are compiled in a
histogram for the three LSC variants: MMST mapping with Wigner estimator
(blue), SC estimator (red), and spin mapping (green).

Table 2 summarizes the breakdown of the trajectory types obtained in the three
different LSC implementations. We note that the energy conservation along a
trajectory is generally considered poor in ab initio implementations due to
the self-consistent field calculations at every time step 31. For this reason,
the effect of the initial conditions on the energy conservation of resulting
trajectories can be masked by that of electronic structure calculations. In
FIG. 5, we histogram the energy jump between consecutive time steps for all
300 trajectories generated in each LSC implementation. Recall that trajectory
initial conditions for the nuclei are identical in all three implementations
as is the Hamiltonian for dynamics; they only differ in the way the electronic
mapping variables are sampled, and in the form of the population estimator
employed at time $t$. Interestingly, FIG. 5 does capture the dependence on
initial conditions — we find that 92% of the trajectories in MMST mapping with
Wigner estimator exhibit energy jumps of less than 1 meV, whereas this numbers
drops to 82% for the MMST with the Semiclassical estimator and to 88% for spin
mapping. All three simulations also yield fragmention trajectories that
correspond to H2 molecule elimination and H-atom dissociation, two channels
that have been observed experimentally. All three also overestimate the
likelihood of C-C bond dissociation but as discussed we believe this error can
be attributed to electronic structure rather than the dynamics themselves.

Table 2:  Results from the 3 variants of LSC-IVR are shown. Succesful trajectories are identified as those that undergo experimentally observed pathways - photoisomerization, $\text{H}_{2}$ elimination, and H atom elimination. C-C bond cleavage is only seen at extremely high excitation energies. The 1 "other" trajectory observed in the SC + MMST mapping and spin mapping LSC implementations corresponds to an outlier trajectory that results in the molecule separating into individual atoms. We also report the average $S_{1}$ to $S_{0}$ excitation energy at $t=0$ for each subset of trajectories in eV as $\langle\Delta E_{S_{1}S_{0}}\rangle$. LSC-IVR variants | Wigner + MMST mapping | SC + MMST mapping | Spin mapping
---|---|---|---
| Count | $\langle\Delta E_{S_{1}S_{0}}\rangle$ | Count | $\langle\Delta E_{S_{1}S_{0}}\rangle$ | Count | $\langle\Delta E_{S_{1}S_{0}}\rangle$
Isomerization | $144$ | $9.65\pm 0.44$ | $47$ | $9.62\pm 0.48$ | $73$ | $9.61\pm 0.46$
$\text{H}_{2}$ elimination | $78$ | $9.78\pm 0.40$ | $128$ | $9.75\pm 0.42$ | $119$ | $9.75\pm 0.51$
H atom elimination | $22$ | $9.76\pm 0.46$ | $74$ | $9.82\pm 0.47$ | $45$ |  $9.71\pm 0.47$
C-C cleavage | $56$ | $10.06\pm 0.41$ | $50$ | $9.93\pm 0.54$ | $62$ |  $9.97\pm 0.38$
Other | $0$ | – | $1$ | – | $1$ | –
Total | $300$ | $9.77\pm 0.45$ | $300$ | $9.78\pm 0.47$ | $300$ | $9.75\pm 0.48$
Figure 6: We compare the thermal population correlation functions obtained
from LSC simulations that employ the MMST mapping with Wigner estimator
("$-$"), SC estimator ("$\blacktriangledown$"), and the spin mapping
("$\times$"). The initial photo-excited $S_{1}$ state population is shown in
red, the doubly excited $S_{2}$ state population is shown in blue, and the
ground $S_{0}$ state population is shown in green. These results were obtained
by averaging the number of the isomerization trajectories shown in Table 2.

All three LSC implementations exhibit qualitatively similar population
dynamics and rates of quenching: the $S_{2}$ state is quickly populated $\sim
10$ fs after photoexcitation, and it takes $60$ to $80$ fs for $50$% of
$S_{1}$ population to decay to the ground state, as shown in FIG. 6. In
comparing the three approaches, we pay particular attention to the issue of
negative electronic state populations: it is well known that the classical
dynamics employed here to time evolve the electronic mapping variables
preserves the sum of the individual state populations, but does not constrain
individual state populations to take on values between 0 and 1. Effectively,
the classical dynamics employed here fails to properly constrain the mapping
variables to the quantum mechanically allowed phase space. [77] Although
individual trajectories might explore unphysical values of state population in
all three LSC implementations, we use the ensemble average populations shown
in FIG. 6 to identify the ‘best’ choice; it is clear that LSC with MMST
mapping and Wigner populations as well as spin mapping LSC yield ensemble
average state populations that are between 0 and 1, whereas the SC estimator
in the MMST mapping framework yields a significant negative value for the
$S_{2}$ state population.

Based on previous studies, we expected spin mapping LSC to significantly
outperform the MMST mapping approaches, but the trajectory breakdown in Table
2 and energy conservation shown in FIG. 5 make it clear that is not
necessarily the case. We find that both in terms of energy conservation and
the overall number of C-C bond cleavage trajectories, MMST with Wigner
estimator emerges as the better implementation. We further note that careful
analysis of energy conservation in all three implementations yields no
correlation between the frequency with which individual trajectories exhibit
negative populations and trajectories that fail to conserve energy. While this
is likely something that should be analyzed on a case-by-case basis, this
finding does provide some reassurance that an individual trajectory exhibiting
negative electronic state population at a given time does not always lead to
unphysical behavior in terms of the overall system dynamics.

##  IV Conclusion

We make the case that employing ab initio classical trajectories within the
semiclassical LSC mapping framework shows great promise as an efficient and
accurate on-the-fly simulation technique for the study of nonadiabatic
processes. In the specific case of ethylene, we show that the use of an
ensemble of trajectories in LSC allows us to directly calculate correlation
functions, and to identify different, statistically significant reaction
pathways with support from previous experimental observations and theoretical
simulations.

We also explore three different implementations of LSC that differ in the
initial conditions used for the electronic mapping variables and the
electronic state population estimators. We establish that although all three
yield relatively similar population dynamics, MMST mapping with Wigner
estimators emerges as the best choice in this study on the basis of energy
conservation, positive ensemble average electronic populations, and the
relatively small number of C-C bond cleavage trajectories. Given previous
studies highlighting the favorable properties of spin mapping, however,
further case studies must be made before this observation can be generalized.

Finally, we discuss the limitations of this approach at present. Like all ab
initio dynamics, it is clear that the level of electronic structure plays a
significant role in the overall accuracy of our findings. Setting that aside,
we note that while the LSC implementation is streamlined and involves no free
parameters such as hopping probabilities, decoherence corrections, or spawning
thresholds, there are a few advances that would allow us to move towards even
more efficient implementation. Most notably, the development of improved
integrators that will allow us to implement dynamics under the symmetrized
Hamiltonian with a larger time step, and a more rigorous study of the
dependence of mapping variable dynamics on the initial conditions to minimize
the number of trajectories that yield negative electronic state populations.

###### Acknowledgements.

The authors acknowledge funding through NSF CAREER Award No. CHE-1555205. The
authors thank Prof. Ben Levine for helpful discussions about the theory and
implementation of AIMS.

## References

  * Bircher  _et al._ [2017] M. P. Bircher, E. Liberatore, N. J. Browning, S. Brickel, C. Hofmann, A. Patoz, O. T. Unke, T. Zimmermann, M. Chergui, P. Hamm, U. Keller, M. Meuwly, H.-J. Woerner, J. Vaníček, and U. Rothlisberger, Struct. Dyn. 4, 061510 (2017).
  * Richings and Habershon [2018] G. W. Richings and S. Habershon, J. Chem. Phys. 148, 134116 (2018).
  * Worth, Robb, and Lasorne [2008] G. A. Worth, M. A. Robb, and B. Lasorne, Molecular Physics 106, 2077–2091 (2008).
  * Martínez, Ben-Nun, and Levine [1996] T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. 100, 7884–7895 (1996).
  * Makhov  _et al._ [2014] D. V. Makhov, W. J. Glover, T. J. Martínez, and D. V. Shalashilin, J. Chem. Phys. 141, 054110 (2014).
  * Ben-Nun and Martínez [2002] M. Ben-Nun and T. J. Martínez, “Ab initio quantum molecular dynamics,” in  _Advances in Chemical Physics_ (John Wiley & Sons, Ltd, 2002) pp. 439–512, https://onlinelibrary.wiley.com/doi/pdf/10.1002/0471264318.ch7 .
  * Fedorov  _et al._ [2017] D. A. Fedorov, S. Seritan, S. B. Fales, T. J. Martínez, and B. G. Levine, J. Chem. Theory Comput. 16, 5485–5498 (2017).
  * Curchod and Martínez [2018] B. F. E. Curchod and T. J. Martínez, Chem. Rev. 118, 3305–3336 (2018).
  * Tully [1990] J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
  * Tully [1998] J. C. Tully, Faraday. Discuss. 110, 407–419 (1998).
  * Doltsinis [2002] N. L. Doltsinis, NIC Series 10, 377–397 (2002).
  * Bittner and Rossky [1995] E. R. Bittner and P. J. Rossky, J. Chem. Phys. 103, 8130 (1995).
  * Schwartz  _et al._ [1996] B. J. Schwartz, E. R. Bittner, O. V. Prezhdo, and P. J. Rossky, J. Chem. Phys. 104, 5942 (1996).
  * Prezhdo and Rossky [1997] O. V. Prezhdo and P. J. Rossky, J. Chem. Phys. 107, 5863 (1997).
  * Volobuev [2000] Y. L. Volobuev, J. Chem. Phys. 112, 9716 (2000).
  * Hack and Truhlar [2001] M. D. Hack and D. G. Truhlar, J. Chem. Phys. 114, 9305 (2001).
  * Bedard-Hearn, Larsen, and Schwartz [2005] M. J. Bedard-Hearn, R. E. Larsen, and B. J. Schwartz, J. Chem. Phys. 123, 234106 (2005).
  * Jain, Alguire, and Subotnik [2016] A. Jain, E. Alguire, and J. E. Subotnik, J. Chem. Theory Comput. 12, 5256–5268 (2016).
  * Habershon  _et al._ [2013] S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. M. III, Annu. Rev. Phys. Chem. 64, 387–413 (2013).
  * Kundu and Makri [2022] S. Kundu and N. Makri, Annu. Rev. Phys. Chem. 73, 349–75 (2022).
  * Ananth [2022] N. Ananth, Annu. Rev. Phys. Chem. 73, 299–322 (2022).
  * Malpathak, Church, and Ananth [2022] S. Malpathak, M. S. Church, and N. Ananth, J. Phys. Chem. A 126, 6359–6375 (2022).
  * Shi and Geva [2003] Q. Shi and E. Geva, J. Chem. Phys. 118, 8173 (2003).
  * Miller [2001] W. H. Miller, J. Phys. Chem. A 105, 2942–2955 (2001).
  * Cotton and Miller [2013a] S. J. Cotton and W. H. Miller, J. Phys. Chem. A 117, 7190–7194 (2013a).
  * Cotton and Miller [2013b] S. J. Cotton and W. H. Miller, J. Chem. Phys. 139, 234112 (2013b).
  * Cotton and Miller [2019] S. J. Cotton and W. H. Miller, J. Chem. Phys. 150, 104101 (2019).
  * Cotton, Igumenshchev, and Miller [2014] S. J. Cotton, K. Igumenshchev, and W. H. Miller, J. Chem. Phys. 141, 084104 (2014).
  * Talbot, Head-Gordon, and Cotton [2022] J. J. Talbot, M. Head-Gordon, and S. J. Cotton, Molecular Physics (2022), 10.1080/00268976.2022.2153761.
  * Zhou, Mandal, and Huo [2019] W. Zhou, A. Mandal, and P. Huo, J. Phys. Chem. Lett. 10, 7062–7070 (2019).
  * Weight, Mandal, and Huo [2021] B. M. Weight, A. Mandal, and P. Huo, J. Chem. Phys. 155, 084106 (2021).
  * Bellonzi, Jain, and Subotnik [2016] N. Bellonzi, A. Jain, and J. E. Subotnik, J. Chem. Phys. 144, 154110 (2016).
  * Stock and Thoss [1997] G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997).
  * Thoss and Stock [1999] M. Thoss and G. Stock, Phys. Rev. A 59, 64 (1999).
  * Meyer, Manthe, and Cederbaum [1990] H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73–78 (1990).
  * Miller and Cotton [2016] W. H. Miller and S. J. Cotton, Faraday Discuss. 195, 9 (2016).
  * Runeson and Richardson [2019] J. E. Runeson and J. O. Richardson, J. Chem. Phys. 151, 044119 (2019).
  * Runeson and Richardson [2020] J. E. Runeson and J. O. Richardson, J. Chem. Phys. 152, 084110 (2020).
  * Ananth, Venkataraman, and Miller [2007] N. Ananth, C. Venkataraman, and W. H. Miller, J. Chem. Phys. 127, 084114 (2007).
  * Cotton, Liang, and Miller [2017] S. J. Cotton, R. Liang, and W. H. Miller, J. Chem. Phys. 147, 064112 (2017).
  * Sun, Wang, and Miller [1998] X. Sun, H. Wang, and W. H. Miller, J. Chem. Phys. 109, 7064 (1998).
  * Meyer and Miller [1979a] H. D. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979a).
  * Meyer and Miller [1980] H. D. Meyer and W. H. Miller, J. Chem. Phys. 72, 2272 (1980).
  * Pierre  _et al._ [2017] S. Pierre, J. R. Duke, T. J. H. Hele, and N. Ananth, J. Chem. Phys. 147, 234103 (2017).
  * Saller, Kelly, and Richardson [2019] M. A. C. Saller, A. Kelly, and J. O. Richardson, J. Chem. Phys. 150, 071101 (2019).
  * Duke and Ananth [2015] J. R. Duke and N. Ananth, J. Phys. Chem. Lett. 6, 4219–4223 (2015).
  * Meyer and Miller [1979b] H. D. Meyer and W. H. Miller, J. Chem. Phys. 71, 2156 (1979b).
  * Cotton and Miller [2015] S. J. Cotton and W. H. Miller, J. Phys. Chem. A 119, 12138–12145 (2015).
  * Bossion  _et al._ [2022] D. Bossion, W. Ying, S. N. Chowdhury, and P. Huo, J. Chem. Phys. 157, 084105 (2022).
  * Stratonovich [1957] R. L. Stratonovich, Soviet Physics JETP 4, 891 (1957).
  * Brif and Mann [1998] C. Brif and A. Mann, J. Phys. A: Math. Gen. 31, L9–L17 (1998).
  * Brif and Mann [1999] C. Brif and A. Mann, Phys. Rev. A 59, 971 (1999).
  * Bonella and Coker [2003] S. Bonella and D. F. Coker, J. Chem. Phys. 118, 4370 (2003).
  * Barca  _et al._ [2020] G. M. J. Barca, C. Bertoni, L. Carrington, D. Datta, N. De Silva, J. E. Deustua, D. G. Fedorov, J. R. Gour, A. O. Gunina, E. Guidez, T. Harville, S. Irle, J. Ivanic, K. Kowalski, S. S. Leang, H. Li, W. Li, J. J. Lutz, I. Magoulas, J. Mato, V. Mironov, H. Nakata, B. Q. Pham, P. Piecuch, D. Poole, S. R. Pruitt, A. P. Rendell, L. B. Roskop, K. Ruedenberg, T. Sattasathuchana, M. W. Schmidt, J. Shen, L. Slipchenko, M. Sosonkina, V. Sundriyal, A. Tiwari, J. L. Galvez Vallejo, B. Westheimer, M. Wloch, P. Xu, F. Zahariev, and M. S. Gordon, “Recent developments in the general atomic and molecular electronic structure system,” The Journal of Chemical Physics 152, 154102 (2020).
  * Ben-Nun and Martínez [1998] M. Ben-Nun and T. J. Martínez, Chem. Phys. Lett. 298, 57–65 (1998).
  * Ben-Nun and Martínez [2000] M. Ben-Nun and T. J. Martínez, Chem. Phys. 259, 237–248 (2000).
  * Mori  _et al._ [2012] T. Mori, W. J. Glover, M. S. Schuurman, and T. J. Martínez, J. Phys. Chem. A 116, 2808–2818 (2012).
  * Allison  _et al._ [2012] T. K. Allison, H. Tao, W. J. Glover, T. W. Wright, A. M. Stooke, C. Khurmi, J. van Tilborg, Y. Liu, R. W. Falcone, T. J. Martínez, and A. Belkacem, J. Chem. Phys. 136, 124317 (2012).
  * Barbatti, Ruckenbauer, and Lischka [2005] M. Barbatti, M. Ruckenbauer, and H. Lischka, J. Chem. Phys. 122, 174307 (2005).
  * Satyapal  _et al._ [1990] S. Satyapal, G. W. Johnston, R. Bersohn, and I. Oref, J. Chem. Phys. 93, 6398 (1990).
  * Balko, Zhang, and Lee [1992] B. A. Balko, J. Zhang, and Y. T. Lee, J. Chem. Phys. 97, 935 (1992).
  * Lin  _et al._ [1998] J. J. Lin, D. W. Hwang, Y. T. Lee, and X. Yang, J. Chem. Phys. 109, 2979 (1998).
  * Lin  _et al._ [2000] J. J. Lin, C. C. Wang, Y. T. Lee, and X. Yang, J. Chem. Phys. 113, 9668 (2000).
  * Farmanara, Stert, and Radloff [1998] P. Farmanara, V. Stert, and W. Radloff, Chem. Phys. Lett. 288, 518–522 (1998).
  * Farmanara  _et al._ [1999] P. Farmanara, O. Steinkellner, M. Wick, M. Wittmann, G. Korn, V. Stert, and W. Radloff, J. Chem. Phys. 111, 6264 (1999).
  * Stert  _et al._ [2004] V. Stert, H. Lippert, H.-H. Ritze, and W. Radloff, Chem. Phys. Lett. 388, 144–149 (2004).
  * Lee, Lee, and Yang [2004] S.-H. Lee, Y. T. Lee, and X. Yang, J. Chem. Phys. 120, 10983 (2004).
  * Lee, Lee, and Lee [2006] S.-H. Lee, Y.-C. Lee, and Y. T. Lee, J. Phys. Chem. A 110, 2337–2344 (2006).
  * Mestdagh  _et al._ [2000] J. M. Mestdagh, J. P. Vesticot, M. Elhanine, and B. Soep, J. Phys. Chem. A 113, 237 (2000).
  * Kosma  _et al._ [2008] K. Kosma, S. A. Trushin, W. Fuss, and W. E. Schmid, J. Phys. Chem. A 112, 7514–7529 (2008).
  * Karashima  _et al._ [2022] S. Karashima, A. Humeniuk, W. J. Glover, and T. Suzuki, J. Phys. Chem. A 126, 3873–3879 (2022).
  * Champenois  _et al._ [2016] E. G. Champenois, N. H. Shivaram, T. W. Wright, C.-S. Yang, A. Belkacem, and J. P. Cryan, J. Chem. Phys. 144, 014303 (2016).
  * Levine and Martínez [2007] B. G. Levine and T. J. Martínez, Annu. Rev. Phys. Chem. 58, 613–634 (2007).
  * Chang  _et al._ [1998] A. H. H. Chang, A. M. Mebel, X.-M. Yang, S. H. Lin, and Y. T. Lee, J. Chem. Phys. 109, 2748 (1998).
  * Pe$\tilde{\text{n}}$a-Gallego, Martínez-Nú$\tilde{\text{n}}$ez, and Vázquez [2002] A. Pe$\tilde{\text{n}}$a-Gallego, E. Martínez-Nú$\tilde{\text{n}}$ez, and S. A. Vázquez, Chem. Phys. Lett. 353, 418–425 (2002).
  * Robin [1985] M. B. Robin,  _Higher Excited States of Polyatomic Molecules_ , Vol. 3 (Elsevier Inc., 1985).
  * Ananth [2013] N. Ananth, J. Chem. Phys. 139, 124102 (2013).