|  Surface hopping molecular dynamics simulation of ultrafast methyl iodide
photodissociation mapped by Coulomb explosion imaging
---|---
|  Yijue Ding∗a, Loren Greenmana and Daniel Rollesa
|  We present a highly efficient method to directly simulate the
photodissociation followed by Coulomb explosion of methyl iodide. In order to
achieve statistical reliability, more than 40,000 trajectories are calculated
on accurate potential energy surfaces of both the neutral molecule and the
doubly charged cation. Non-adiabatic effects during photodissociation are
treated using a Landau-Zener surface hopping algorithm. The simulation is
performed analogous to a recent pump-probe experiment using coincident ion
momentum imaging [Ziaee et al., Phys. Chem. Chem. Phys. 2023, 25, 9999]. At
large pump-probe delays, the simulated delay-dependent kinetic energy release
signals show overall good agreement with the experiment, with two major
dissociation channels leading to $\text{I}(^{2}\text{P}_{3/2})$ and
$\text{I}^{*}(^{2}\text{P}_{1/2})$ products. At short pump-probe delays, the
simulated kinetic energy release shows a clear bifurcation near 12 fs, owing
to non-adiabatic transitions through a conical intersection. The developed
method is particularly suitable and efficient in simulating processes that
highly rely on statistics or for identifying rare reaction channels.

##

††footnotetext: a J. R. Macdonald Laboratory, Department of Physics, Kansas
State University, Manhattan, KS 66506, USA. E-mail<EMAIL_ADDRESS>
##  1 Introduction

Photochemical reactions are vital and ubiquitous in our daily lives.
Representative examples include photosynthesis, human vitamin D formation and
light-harvesting devices1. Photochemical reactions under low temperature also
play an important role in quantum control and quantum information processing2,
3, 4, 5. Following molecular structures in real time during photochemical
reactions is vital to understand and control the reactions for specific
purposes. To this end, ultrafast imaging techniques are developed to study
reaction dynamics in a time-resolved manner. These techniques include
attosecond transient absorption spectroscopy (ATAS)6 and time-resolved
photoelectron spectroscopy (TRPES)7, which provides indirect information on
nuclear geometries by measuring the electronic energy levels or the kinetic
energy of photoelectrons with extremely high temporal resolution. Other
techniques such as ultrafast electron diffraction (UED)8 and Coulomb explosion
imaging (CEI)9 provides direct information of the nuclear geometry by
interpreting the diffraction pattern or measuring the momenta of ionic
fragments.

In particular, CEI is a powerful technique that is equally sensitive to both
light and heavy atoms. CEI retrieves molecular structure by rapidly removing
multiple electrons and measuring the asymptotic momenta of ionic fragments.
Thanks to decades of development, CEI can be applied to molecules consisting
of a couple of atoms to over ten atoms in a molecule10, 11, 12, 13. In a time-
resolved manner, CEI has been used to identify isomerization processes14 or
isolate a rare roaming channel15.

Methyl iodide ($\text{CH}_{3}\text{I}$) is a prototype molecule to study
photochemical reaction dynamics due to its simplicity and rich variety of
photochemical phenomena, and thereby has attracted research interests for a
long time16, 17, 18, 19, 20, 21. CEI experiments have been conducted on this
molecule both near the ground state equilibrium22, 23 and during
photodissociation in a time-resolved manner24, 25, 26, 27. One intriguing
feature is the nondiabatic dynamics during photodissociation owing to a
conical intersection (CI) induced by the strong spin-orbit coupling effect.
Non-adiabtic transitions near such a CI result in two dissociation channels
corresponding to $\text{I}(^{2}\text{P}_{3/2})$ and
$\text{I}^{*}(^{2}\text{P}_{1/2})$ products, respectively.

In the theoretical side, several trajectory calculations have been performed
on potential energy surfaces (PES) of $\text{CH}_{3}\text{I}$ excited states
over time28, 29, 30, as well as ab initio molecular dynamics (AIMD)
calculations31, 32 and quantum mechanical calculations33, 34, 35 that have
been conducted to study the A-band photodissociation. However, accurate
simulations of both the photodissociation and the Coulomb explosion process to
quantitatively simulate the observvables from the pump-probe experiments have
not been performed so far. Particularly, no molecular dynamics simulations
have been conducted for the doubly charged ionic state, which is important for
the Coulomb explosion process.

In this work, we study the UV induced dissociation of $\text{CH}_{3}\text{I}$
probed by CEI in a time-resolved manner through trajectory calculations on ab
initio PESs of both the neutral molecule and the doubly charged cation. Since
CEI is a statistical method that relies on a large number of events,
theoretical methods like AIMD that calculate electronic structures on-the-fly
are of limited use as they can typically only predicted a few exemplary
trajectories due to high computational cost. To address this issue, we develop
an efficient program that performs fast trajectory calculations on both the
neutral and ionic surfaces with the inclusion of non-adiabatic effects.

##  2 Theory and computational methods

###  2.1 Ab initio electronic structure calculation

All electronic structures are computed at the multi-reference configuration
interaction (MRCI)36, 37 level of theory. The active space is 6 electrons in 5
orbitals for both the neutral $\text{CH}_{3}\text{I}$ molecule and the
$\text{CH}_{3}\text{I}^{2+}$ cation. The active space includes two nearly
degenerate lone pair ($n$) orbitals of the I atom, the C-I bonding ($\sigma$)
and antibonding ($\sigma^{*}$) orbitals, for the purpose of accurately
describing the A-band excited states that are dominated by
$n\rightarrow\sigma^{*}$ transitions. One additional orbital is added to the
active space to improve convergence.

The state-averaged multi-configuration self-consistent field (SA-MCSCF)38, 39,
40 calculation is first carried out to obtain the reference wavefunctions. For
the neutral $\text{CH}_{3}\text{I}$ molecule, 3 lowest singlet states and 3
lowest triplet states are computed and optimized with equal weights in the SA-
MCSCF procedure to cover all A-band excited states. For the
$\text{CH}_{3}\text{I}^{2+}$ cation, 7 singlets and 4 triplets are computed
and optimized with equal weights in the SA-MCSCF procedure to cover all
possible atomic terms of the $\text{I}^{+}$ cation. The MRCI calculation is
then performed by adding single and double excitations to external
configurations from the MCSCF reference wavefunctions.

The spin-orbit (SO) coupling effect is treated using a state-interaction (SI)
method. The SO eigen-energies are calculated by diagonalizing the SO coupled
Hamiltonian, $\hat{H}_{el}+\hat{H}_{SO}$, under the basis of MCSCF wave
functions, with the diagonal elements replaced by MRCI energies.

The Dunning type correlation-consistent basis sets41 are employed in our
electronic structure calculations. Specifically, we use cc-pVTZ basis sets for
hydrogen and carbon atoms. We also use cc-pVTZ-PP basis set for the iodine
atom, where the inner 28 core electrons are described by a relativistic
pseudo-potential42, 43.

All ab initio calculations are performed using the MOLPRO quantum chemistry
package44, 45.

###  2.2 Analytic potential energy surfaces

In order to obtain an accurate description of the electronic energies and the
corresponding gradients, we calculate over 40,000 geometries for
$\text{CH}_{3}\text{I}$ and over 20,000 geometries for
$\text{CH}_{3}\text{I}^{2+}$. The PES is constructed under
$\\{R,\alpha,\theta\\}$ degrees of freedom, where R is the C-I distance,
$\alpha$ is the angle between the C-H vector and the $C_{3}$ axis of symmetry,
and $\theta$ is the angle between the C-I vector and the $C_{3}$ axis of
symmetry. Other degrees of freedom are kept at equilibrium geometry.
Particularly, ab initio electronic structure calculations are performed on a
very dense grid ($\Delta R=0.01\text{\AA}$) near the conical intersection
between ${}^{3}\text{Q}_{0}$ and ${}^{1}\text{Q}_{1}$ states in order to
accurately characterize the nonadiabatic effects. Analytic representations of
the adiabatic PESs are generated using a multi-dimensional cubic spline
interpolation algorithm. Because the grid is sufficiently dense, such an
direct interpolation is able to characterize the conical intersection
precisely.

The ab initio potential energies are calculated up to $R_{c}=5$ Å for
$\text{CH}_{3}\text{I}$ and up to $R^{ion}_{c}=10$ Å for
$\text{CH}_{3}\text{I}^{2+}$, beyond which the energy gradient of
$\text{CH}_{3}\text{I}$ is extremely
small($|dV_{\text{CH}_{3}\text{I}}/dR|<0.01$ eV/Å), and the potential of
$\text{CH}_{3}\text{I}^{2+}$ is almost purely
Coulombic($|d(V_{\text{CH}_{3}\text{I}^{2+}}-1/R)/dR|<0.01$ eV/Å). In order to
simulate the the pump-probe process at large delays (up to 500 fs) and to
obtain the asymptotic momenta of the ionic fragments after Coulomb explosion,
it is natural to extend the PES to infinite C-I distance.

The interaction between the neutral I and $\text{CH}_{3}$ fragments is
dominated by the van der Waals force at large internuclear distances, which is
in the order of $R^{-6}$. This high-order interaction is negligible, and
therefore the $\text{CH}_{3}\text{I}$ potential beyond $R_{c}$ is written as

$V_{\text{CH}_{3}\text{I}}\approx V_{CH_{3}}+V_{I}\qquad(R>R_{c}),$ (1)

where $V_{\text{CH}_{3}}$ only depends on the internal degrees of freedom of
the $\text{CH}_{3}$ fragment, and $V_{\text{I}}$ is the internal energy of the
I atom, which also defines the thresholds of the different reaction channels.

The interaction between the ionic $\text{I}^{+}$ and $\text{CH}_{3}^{+}$
fragments is dominated by the Coulomb force at large internuclear distances.
We also include a $1/R^{4}$ term to describe the charge-dipole interaction.
Thus, the $\text{CH}_{3}\text{I}^{2+}$ potential beyond $R^{ion}_{c}$ is
written as

$V_{\text{CH}_{3}\text{I}^{2+}}\approx
V_{\text{CH}_{3}^{+}}+V_{\text{I}^{+}}+(\frac{1}{R}-\frac{C_{4}}{R^{4}})\qquad(R>R^{ion}_{c}),$
(2)

where the coefficient $C_{4}$ is obtained by fitting the long range tails of
the calculated potential energies.

###  2.3 Photodissociation-Coulomb explosion trajectory calculation

In order to directly simulate the observables from time-resolved CEI
experiments, we develop a molecular dynamics program to simulate
photodissociation followed by Coulomb explosion in a single trajectory
calculation on pre-built analytic PESs of $\text{CH}_{3}\text{I}$ and
$\text{CH}_{3}\text{I}^{2+}$. At the first step, the initial conditions are
sampled from a classical Boltzmann distribution of normal mode coordinates of
the ground electronic state at a certain effective temperature. For UV photon
absorption near 266 nm, transition to the ${}^{3}\text{Q}_{0}$ state (mostly
parallel) is much stronger than to other states in the A-band20, 46. The
momenta and positions of each atom do not change during the UV absorption
(vertical transition). Thus, similar to other AIMD simulations performed on
this molecule32, 31, 26, we start all trajectories on the ${}^{3}\text{Q}_{0}$
electronic state with initial conditions sampled on the ground electronic
state. In our simulation, the equation of motion (EOM) is described by the
Euler-Lagrange equation $\frac{d}{dt}\frac{\partial
L}{\partial\dot{\textbf{q}}}-\frac{\partial L}{\partial\textbf{q}}=0$, where
$\textbf{q}=\\{\Psi,\Theta,\Phi,R,\alpha,\theta\\}$. The Euler angles
$\\{\Psi,\Theta,\Phi\\}$ are included in the EOM to describe the overall
rotation of the entire molecule. The molecule is initially prepared at its
rotational ground state with random orientation. To propagate the trajectory,
the EOM is integrated numerically using an adaptive Runge-Kutta algorithm with
a desired relative error of $10^{-7}$. Non-adiabatic effects during the
propagation are taken into account using a Landau-Zener surface hopping
algorithm, which will be discussed in the following subsection. When the
trajectory is calculated up to a pre-specified target time, the probe laser
pulse is applied to doubly ionize the molecule and induces the Coulomb
explosion. We refer to this target time as the pump-probe delay. The pulse
duration and its effects on the nuclear dynamics are not considered in our
simulation. So the molecular dynamics are always conducted on field-free PESs.
The photoionization event itself is assumed to occur instantaneously on the
timescales considered in our simulation. Therefore, the trajectory
instantaneously hops to the $\text{CH}_{3}\text{I}^{2+}$ PES and continues to
propagate until the Coulomb energy between ionic fragments is below 0.01 eV.
Surface hoppings are not considered on the $\text{CH}_{3}\text{I}^{2+}$ PES.
We call this type of trajectories photodissociation-Coulomb explosion (PD-CE)
trajectories.

Since CEI is highly dependent on statistics, and one trajectory corresponds to
one single event in the coincident ion momentum measurement, a large number of
trajectories are required to obtain a statistically reliable interpretation of
the CEI results. In our simulation, more than 40,000 PD-CE trajectories are
calculated in total, specifically, about 20,000 in the asymptotic regime and
about 20,000 for short pump-probe delays.

###  2.4 Landau-Zener surface hopping algorithm

In the present study, we employ a Landau-Zener model to treat the non-
adiabatic effects in the trajectory calculation. The Landau-Zener surface
hopping (LZSH) algorithm has been elaborated in previous studies, and was
shown to yields reliable results47, 48, 49, 50. Although Tully’s fewest-switch
surface hopping (FSSH) algorithm51 is commonly used for modeling non-adiabatic
dynamics in a classical manner, it cannot be rigorously derived from quantum
mechanics, and thus the time-evolution of the electron wave function is
inconsistent with the classical motion of nuclei. One advantage of LZSH is
that it can be readily implemented on adiabatic PESs. Moreover, frustrated
hops are less likely to occur in the Landau-Zener model.

The Landau-Zener non-adiabatic transition probability in the adiabatic picture
is written as52

$P_{i\rightarrow
j}=\text{exp}\left(-\frac{\pi}{2\hbar}\sqrt{\frac{Z_{ij}(\textbf{q}(t_{0}))^{3}}{\frac{d^{2}}{dt^{2}}Z_{ij}(\textbf{q}(t))|_{t=t_{0}}}}\right)$
(3)

where $\textbf{q}(t)$ is the nuclear coordinate vector as a function of time,
and $Z_{ij}=|V_{i}-V_{j}|$ is the absolute value of the potential energy
difference between $i$ and $j$ adiabatic states. During the trajectory
propagation, the value of $Z_{ij}$ is monitored at each time step. When a
local minimum of $Z_{ij}$ as a function of time is attained, the non-adiabatic
transition probability is calculated using Eq. (3). Then, a random number is
generated from a uniform distribution between 0 and 1, and is compared with
the probability $P_{i\rightarrow j}$ to determine whether the trajectory
should hop to a different adiabatic state. If a hop is confirmed, the velocity
of the system is rescaled by a factor $\eta$ to maintain the conservation of
energy. This requirement can be written as

$\frac{1}{2}\dot{\textbf{q}}^{T}\textbf{M}\dot{\textbf{q}}+V_{i}(\textbf{q})=\frac{1}{2}\eta\dot{\textbf{q}}^{T}\textbf{M}\eta\dot{\textbf{q}}+V_{j}(\textbf{q}),$
(4)

where $\dot{\textbf{q}}=d\textbf{q}/dt$ is the velocity vector of the nuclear
coordinates, and M is the mass matrix.

##  3 Results and discussions

###  3.1 Adiabatic potential energies of $\text{CH}_{3}\text{I}$ and
$\text{CH}_{3}\text{I}^{2+}$

The adiabatic PESs of $\text{CH}_{3}\text{I}$ and $\text{CH}_{3}\text{I}^{2+}$
are constructed based on high-level ab initio electronic structure
calculations and multi-dimensional cubic spline interpolations. Figure 1 shows
the potential energies of $\text{CH}_{3}\text{I}$ and
$\text{CH}_{3}\text{I}^{2+}$ by stretching the C-I bond length $R$. Fifteen
ionic potential curves with thresholds corresponding to different internal
states of $\text{I}^{+}$ are shown in Fig. 1(a), and are also compared with
the pure Coulomb interactions between $\text{CH}_{3}^{+}$ and $\text{I}^{+}$.
The energies of the ionic states at the Franck-Condon point correspond to the
double ionization potentials of $\text{CH}_{3}\text{I}$. These energies range
from about 26 eV to 31 eV in our calculation, in good agreement with previous
theoretical and experimental studies53, 54. The actual energies of all states
are smaller than the pure Coulomb energy, owing to the attractive covalent
interaction between $\text{CH}_{3}^{+}$ and $\text{I}^{+}$. This attractive
interaction is particularly strong for small $R$. As the C-I distance
increases, the potential becomes strongly repulsive and can be well
approximated by pure Coulomb interactions for $R>4$ Å.

Twelve neutral potential curves corresponding to the ground state and A-band
excited states are shown in Fig. 1(b). At large C-I distance, these states
dissociate to two atomic states of iodine, $\text{I}(^{2}\text{P}_{3/2})$ and
$\text{I}^{*}(^{2}\text{P}_{1/2})$, with an energy splitting of 0.9 eV in our
calculation. Using the Mulliken notation16, the ${}^{3}Q_{0}$, ${}^{1}Q_{1}$
and ${}^{3}Q_{1}$ states are dipole allowed transitions from the ground state.
An UV pulse near 260 nm mostly excites the molecule to the ${}^{3}Q_{0}$
state20 through parallel transition, with a photoabsorption cross section much
larger than the other states. It is noted that a non-adiabatic crossing
between ${}^{3}Q_{0}$ and ${}^{1}Q_{1}$ states is present at $R=2.38$ Å, which
will induce population transfer to the ${}^{1}Q_{1}$ state during
photodissociation.

Fig. 1: Ab initio potential energy curves for the $\text{CH}_{3}\text{I}^{2+}$
cation (a) and the neutral $\text{CH}_{3}\text{I}$ molecule (b) by stretching
the C-I bond length $R$. All other degrees of freedom remain the same as in
the equilibrium geometry. The energies are relative to the ground state energy
of $\text{CH}_{3}\text{I}$ at equilibrium. In panel (a), solid curves are
adiabatic potentials with different colors corresponding to different
dissociation thresholds of the $\text{I}^{+}$ ions. Dashed curves are pure
Coulomb potentials shifted to the corresponding $\text{I}^{+}$ thresholds for
comparison. The potential curves in panel (b) correspond to the ground state
and A-band excited states of the neutral molecule. The curves shown in thick
black, blue, and red are the relevant states during UV induced
photodissociation. Fig. 2: Potential energy surfaces of
$\text{CH}_{3}\text{I}^{2+}$ and $\text{CH}_{3}\text{I}$ projected onto the
$\alpha$ (angle between C-H vector and $C_{3}$ axis) and $R$ (C-I distance)
plane. All other degrees of freedom are kept the same as in the equilibrium
geometry such that the molecule(cation) is at $C_{3v}$ symmetry. The black dot
denotes the equilibrium geometry. (a) The lowest adiabatic electronic state of
$\text{CH}_{3}\text{I}^{2+}$ corresponding to the
$\text{I}^{+}(^{3}\text{P}_{2})$ dissociation threshold. (b) The 7th adiabatic
electronic state of $\text{CH}_{3}\text{I}$ where all trajectory propagations
start. The black dashed line represents the seam location where the energies
of ${}^{3}\text{Q}_{0}$ and ${}^{1}\text{Q}_{1}$ states are degenerate, i.e.
the conical intersections. The blue (hop) and red (no hop) circles indicate
the geometries at which a surface hopping decision is made for 100
trajectories. The energy difference between adjacent contours is 0.1 eV in
Panel (a) and 0.14 eV in Panel (b).

In previous CEI experiments, the measured delay-independent kinetic energy
release (KER) show two major peaks centered around 4.4 eV and 5.2 eV for the
doubly charged ion22, 27. Our simulation, as well as previous works22, 27,
reveals that the KER peak near the 4.4 eV can be attributed to Coulomb
explosions on the lowest ionic PES.444The KER curves in Ref.22, 27 are
mislabeled. The labels associated to the $\text{I}^{+}(^{3}\text{P}_{2})$
state and the $\text{I}^{+}(^{1}\text{D}_{2})$ state should be swapped.
Therefore, the molecule is highly likely ionized to the lowest
$\text{CH}_{3}\text{I}^{2+}$ ionic state corresponding to
$\text{I}^{+}(^{3}\text{P}_{2})$ threshold after IR multi-photon absorption,
and we simulate the Coulomb explosion part on this PES for most of our PD-CE
trajectories. Figure 2(a) shows the lowest adiabatic PES of
$\text{CH}_{3}\text{I}^{2+}$ as a function of the C-I distance $R$ and the
angle $\alpha$ when the cation is in $C_{3v}$ symmetry. Contrary to the 1D
potential curve shown in Fig.1(a) and in other literature22, 26, there is no
local minimum on this PES. Consequently, our simulation finds that no
metastable $\text{CH}_{3}\text{I}^{2+}$ cations can be formed for this ionic
state, which is in contrast to the results using 1D wave packet propagation
reported in Ref. 22. By further examination of all of the fifteen states in
Fig. 1(a), only the state corresponding to the
$\text{I}^{+}(^{1}\text{S}_{0})$ threshold has a shallow well near $R=2.25$ Å,
which may support metastable $\text{CH}_{3}\text{I}^{2+}$ cations as will be
discussed with more details in Section 3.4.

Since all trajectories propagate on adiabatic PESs, we also show the PES of
the 7th adiabatic state of $\text{CH}_{3}\text{I}$ in Fig. 2(b). This
adiabatic state is separated into ${}^{3}Q_{0}$ and ${}^{1}Q_{1}$ states by a
"seam" line, which indicates the geometries of the conical intersections. All
PD-CE trajectories are initially prepared on this adiabatic PES near the
Franck-Condon (FC) point, located at $R_{0}=2.16$ Å and $\alpha_{0}=107.6$°.
As shown in Fig. 2(b), when the trajectories propagate close to the seam line,
non-adiabatic transitions are likely to take place according to our Landau-
Zener model and the trajectories may hop to the next adiabatic state.

###  3.2 Simulating the kinetic energy release

In a CEI experiment, the asymptotic momentum of each ionic fragment is
measured via either velocity map imaging (VMI)14, 22, 23, 24, 25, 26 or COLd
Target Recoil-Ion Momentum Spectrometer (COLTRIMS)9, 11, 12, 13, 15, 27. In
the latter case, the measured kinetic energy release (KER) is the sum of the
kinetic energies of $\text{I}^{+}$ and $\text{CH}_{3}^{+}$ fragments detected
in coincidence. In our simulation, this corresponds to the translational
kinetic energy of the relative motion between I($\text{I}^{+}$) and
$\text{CH}_{3}$($\text{CH}_{3}^{+}$) fragments, which can be written as

$KER=T_{trans}=\frac{1}{2}\mu\dot{\textbf{R}}_{I-CH_{3}}^{2},$ (5)

where $\mu$ is the reduced mass between I($\text{I}^{+}$) and
$\text{CH}_{3}$($\text{CH}_{3}^{+}$) fragments, $\textbf{R}_{I-CH_{3}}$ is the
displacement vector between I($\text{I}^{+}$) and the center of mass of
$\text{CH}_{3}$($\text{CH}_{3}^{+}$), and
$\dot{\textbf{R}}_{I-CH_{3}}=d\textbf{R}_{I-CH_{3}}/dt$ is the velocity of the
relative motion.

The translational energy has contributions from two origins: the neutral
excited state potential energy during photodissociation and the ionic state
potential energy during Coulomb explosion. Thus, we further decompose the KER
into two parts,

$KER=KER_{PD}+KER_{CE},$ (6)

where $KER_{PD}$ and $KER_{CE}$ are the kinetic energies gained from
photodissociation and Coulomb explosion, respectively.

Fig. 3: Translational kinetic energy of the relative motion between I and
$\text{CH}_{3}$ fragments as a function of propagation time for 20
trajectories during photodissociation. The blue and red curves correspond to
$\text{I}(^{2}\text{P}_{3/2})$ and $\text{I}^{*}(^{2}\text{P}_{1/2})$
dissociation channels respectively.

Figure 3 shows the translational energy as a function of propagation time
during photodissociation. During the initial 30 fs, the translational energy
increases rapidly due to the strong repulsion of the excited state potential.
Some trajectories (marked by red curves) show small jumps of the translational
energy near 12 fs, which indicates that a surface hopping takes place and the
kinetic energy is adjusted according to Eq. (4) to maintain the conservation
of total energy. The translational energy saturates after about 50 fs for both
$I$ and $I^{*}$ reaction channels, indicating the photodissociation is
complete, and any further change in KER is attributed to Coulomb explosion
contribution. This photodissociation time scale is consistent with the
calculation and measurement using XUV absorption spectroscopy 32, 55, in which
one absorption peak disappears at about 50 fs because of the dipole forbidden
I(${}^{2}\text{P}_{1/2}$) to I(${}^{2}\text{D}_{5/2}$) transition.

When a one-dimensional (1D) model is applied to the PD-CE process, that is,
considering the C-I distance $R$ as the only reaction coordinate and freezing
all other degrees of freedom at the equilibrium geometry, $KER_{PD}$ and
$KER_{CE}$ are approximated to be

$\begin{split}KER_{PD}(t)\approx
V_{\text{CH}_{3}\text{I}}(R_{0})-V_{\text{CH}_{3}\text{I}}(R(t)),\\\
KER_{CE}(t)\approx
V_{\text{CH}_{3}\text{I}^{2+}}(R(t))-V_{\text{CH}_{3}\text{I}^{2+}}(\infty),\\\
t=\int_{R_{0}}^{R}\frac{d\tilde{R}}{\sqrt{2(V_{\text{CH}_{3}\text{I}}(R_{0})-V_{\text{CH}_{3}\text{I}}(\tilde{R}))/\mu}},\end{split}$
(7)

where $t$ is the pump-probe delay. This 1D model is used for interpreting the
experimental data in Ref.27. At large pump-probe delays, the C-I distance $R$
is large and the photodissociation is complete. In that case, $KER_{PD}\approx
V_{\text{CH}_{3}\text{I}}(R_{0})-V_{\text{CH}_{3}\text{I}}^{th}$ is
approximately a constant value and $KER_{CE}\approx 1/R(t)$ is approximately
the pure Coulomb energy.

###  3.3 Coulomb explosion in the asymptotic region

Fig. 4: Kinetic energy release (KER) at different pump-probe delays ranging
from 60 fs to 500 fs. (a) Signals from our simulation of about 20,000 PD-CE
trajectories at T=200 K. The KER equals to the final translational kinetic
energy of the relative motion between $\text{CH}_{3}^{+}$ and $\text{I}^{+}$
fragments. (b) Signals adapted from the experiment in Ref.27. The KER is the
sum of the kinetic energies of $\text{CH}_{3}^{+}$ and $\text{I}^{+}$
fragments measured in a coincident manner for a UV pump intensity of $1\times
10^{13}\text{W}/\text{cm}^{2}$. The yellow and black curves are predictions
using the 1D model in Eq.(7). Signals in both (a) and (b) are counted using a
bin size $h_{t}=20$ fs and $h_{E}=0.01$ eV.

About 20,000 PD-CE trajectories are calculated at T=200 K with the pump-probe
delay ranging from 50 fs to 500 fs and with a step of $\Delta t=1$ fs. Since
previous studies have concluded $\text{CH}_{3}\text{I}$ will be highly likely
ionized to the lowest ionic states after IR multi-photon absorption22, 23, the
Coulomb explosion part is simulated on the lowest PES corresponding to the
$\text{I}^{+}(^{3}\text{P}_{2})$ threshold. The simulated delay-dependent KER
signals are shown in Fig. 4(a), along with the relation between KER and the
pump-probe delay predicted by the 1D model in Eq.(7). The signals form two
major bands, corresponding to the I and $\text{I}^{*}$ reaction channels. The
lower band ($\text{I}^{*}$ channel) is significantly stronger than the upper
band (I channel), indicating that a large fraction of trajectories undergo
non-adiabatic transitions. The branching ratio
$[\text{I}^{*}]/([\text{I}^{*}]+[\text{I}])$ is 0.82 in our simulation. This
value is in overall agreement with simulations using the FSSH algorithm32, 31,
wave packet evolutions35 and experiments56, 57, 58, 59, which range from 0.65
to 0.85.

The two bands show clear correlation between the KER and the pump-probe delay,
which agree with the prediction of the 1D model. The energy of the upper band
is slightly lower than that of the 1D model, because more potential energy is
converted to the vibrational and rotational motion of $\text{CH}_{3}$ fragment
in $\text{I}^{*}$ reaction channel. In fact, over ninety percent of the
potential energy is redistributed to the translational energy of
$\text{CH}_{3}$ and I fragments during photodissociation, while a small
fraction is redistributed to the vibrational and rotational motion of
$\text{CH}_{3}$. Previous studies reveal that higher vibrational and
rotational quantum states of $\text{CH}_{3}$ are excited in $\text{I}^{*}$
dissociation channel31, 33.

The simulated delay-dependent KER signals are also in good agreement with the
experiment, as shown in Fig.4(b), although the experiment shows slightly lower
KER than the 1D model. We only compare our simulation with the experiment for
KER below 4 eV because higher KER region in the experiment is dominated by
delay-independent signals due to multi-photon ionization of ground-state
$\text{CH}_{3}\text{I}$ by the probe pulse, which is beyond the focus of this
work. A weak contribution of lower KER in the experiment is attributed to
multi-photon Rydberg state excitation by the pump pulse, which is not
accounted for in our simulation.

Fig. 5:  Experimental kinetic energy release (KER) signals (dashed black,
adapted from Ref.27) compared with simulated $KER_{PD}$ signals (solid black).
The experimental signals are integrated over the pump-probe delay from 2.5 ps
to 3ps. The intensity of the UV pump is $1\times
10^{13}\text{W}/\text{cm}^{2}$. The simulated $KER_{PD}$ include only the
kinetic energy release during photodissociation, and are integrated over the
pump-probe delay from 50 fs to 100 fs. The blue and red curves with shadows
are separate signals corresponding to $\text{I}^{*}$ and I dissociation
channels respectively in our simulation, scaled by a factor of 1/2. The
vertical lines indicate the potential energy differences between
${}^{3}\text{Q}_{0}$ state at equilibrium geometry and the thresholds of
$\text{I}^{*}$(red) and I(blue) channels, with the relative heights matching
the $[\text{I}^{*}]/([\text{I}^{*}]+[\text{I}])$ branching ratio in our
simulation. Signals are normalized by the total number of events between 1 eV
and 3 eV.

When the pump-probe delay is sufficiently large such that the Coulomb energy
is negligible ($KER_{CE}$<0.05 eV when t>2.5 ps), the measured KER is entirely
contributed from $KER_{PD}$ during photodissociation, which has already
saturated at t=50 fs. Therefore, it is convenient to compare the simulated
$KER_{PD}$ at t>50 fs with the measured KER at t>2.5 ps, as shown in Fig. 5.
The overall shape and peak positions of the simulated $\text{KER}_{PD}$
signals shows good agreement with the experimentally measured KER. The
relative height between the two peaks in our simulation is slightly different
from the experiment, which reflects the branching ratio in our simulation is
higher than the experiment. This can be attributed to the limitation of the
present surface hopping model. Another reason is that we neglect the small
fraction of dipole transitions to the ${}^{1}\text{Q}_{1}$ state that mostly
leads to dissociations in I channel.

###  3.4 Coulomb explosion at short pump-probe delays

Fig. 6: Simulated kinetic energy release at different pump-probe delays
ranging from 0 fs to 50 fs using about 20,000 PD-CE trajectories at T=100 K.
The Coulomb explosion is simulated on the lowest ionic potential energy
surface of $\text{CH}_{3}\text{I}^{2+}$ that corresponds to
$\text{I}^{+}(^{3}\text{P}_{2})$ final dissociation limit. (a) Signals
including all active states during photodissociation, which will lead to both
I and $\text{I}^{*}$ channels. (b) Signals only from 7th adiabatic state
during photodissocaition, which correspond to I channel. (c) Signals only from
9th adiabatic state during photodissociation, which corresponds to
$\text{I}^{*}$ channel at long distance. All signals are counted using a bin
size $h_{t}=0.5$ fs and $h_{E}=0.01$ eV.

At short pump-probe delays (t<50 fs), the $\text{CH}_{3}\text{I}$ molecule is
in the intermediate states of photodissociation. Non-adiabatic transitions
also take place at this time scale. To investigate the non-adiabatic dynamics
in this regime, additional 20,000 PD-CE trajectories are calculated with
initial conditions sampled at T=100 K and with a pump-probe delay step $\Delta
t=0.1$ fs up to 50 fs. Here we choose a lower effective temperature to make
the initial conditions a little more condensed in the phase space such that
the non-adiabatic features are better resolved in KER signals. The Coulomb
explosion part is still simulated on the lowest PES corresponding to
$\text{I}^{+}(^{3}\text{P}_{2})$ threshold. The resulting delay-dependent KER
signals are shown in Fig. 6. A clear bifurcation of signals is present between
10 fs and 20 fs in Fig. 6(a), which indicates the occurrence of non-adiabatic
transitions. The photodissociation separates to two channels leading to I and
I* products. Fig. 6(b) shows that only a small fraction of trajectories (about
20 percent) undergoes adiabatic passage and dissociate to I channel. Most
trajectories (about 80 percent) start to hop to the 9th adiabatic state near
12 fs and dissociate to $\text{I}^{*}$ channel, as shown in Fig. 6(c). It
should be noted that the 1D model is a poor approximation at short pump-probe
delays. The C-I stretching is highly coupled with
$\text{CH}_{3}$($\text{CH}_{3}^{+}$) umbrella motion and the C-I bending at
small C-I distances, which is not characterized by the 1D model. Not only does
the 1D model predict a fake potential well on $\text{I}^{+}(^{3}\text{P}_{2})$
ionic PES, but also it significantly underestimates the KER. Therefore, we do
not compare our simulation with the 1D model at short pump-probe delays.

The bifurcation of KER signals has not been observed in CEI experiments so
far. One reason is that the conical intersection is very close to the FC point
and $\text{CH}_{3}$ is a very light fragment. Thus, non-adiabatic transitions
occur within 20 fs after UV excitation, which requires pump and probe pulses
with very high temporal resolutions to resolve. One option to alleviate this
difficulty might be choosing a molecule with heavier fragments such as
$\text{CF}_{3}\text{I}$. An alternative method to map the non-adiabatic
transitions is using ATAS, which has been studied both theoretically and
experimentally55, 32.

Fig. 7: Simulated $\text{CH}_{3}\text{I}^{2+}$ signals at very short pump-
probe delays. (a) Potential energy surface of $\text{CH}_{3}\text{I}^{2+}$
corresponding to $\text{I}^{+}(^{1}\text{S}_{0})$ threshold, projected to
$\alpha$ and $R$ plane. The energy difference between adjacent contours is 0.1
eV. The Coulomb explosion is simulated on this surface. 4,000 PD-CE
trajectories are calculated at $T=100$ K with pump-probe delays from 0 fs to
10 fs. The black triangles, which are counted as $\text{CH}_{3}\text{I}^{2+}$
signals, indicate the final geometries of 956 PD-CE trajectories after
propagating for over 2 ps on the ionic potential energy surface. (b)
$\text{CH}_{3}\text{I}^{2+}$ signals as a function of pump-probe delay.

As mentioned before, only the ionic state associated with
$\text{I}^{+}(^{1}\text{S}_{0})$ threshold has a shallow well on the PES and
can support a metastable $\text{CH}_{3}\text{I}^{2+}$ state, as shown in Fig.
7(a). But it should be noted that this channel has a rather small transition
probability during the IR multi-photon ionization. In order to investigate the
$\text{CH}_{3}\text{I}^{2+}$ signals, we calculate 4,000 PD-CE trajectories at
very short pump-probe delays (t<10 fs), for which the Coulomb explosion is
simulated on the $\text{I}^{+}(^{1}\text{S}_{0})$ ionic PES. 956 trajectories
form $\text{CH}_{3}\text{I}^{2+}$ states and never break up into ionic
fragments, as shown on top of the PES in Fig. 7(a). The
$\text{CH}_{3}\text{I}^{2+}$ signal has a maximum yield near 1 fs, and damps
rapidly with the pump-probe delay, as shown in Fig. 7(b). Due to this shallow
well on the ionic PES, the $\text{CH}_{3}\text{I}^{2+}$ signals are detected
in the CEI experiment27, and an enhancement of $\text{CH}_{3}\text{I}^{2+}$
yield is observed near zero pump-probe delay, which is attributed to the cross
correlation between the pump and the probe pulses.

##  4 Conclusion and outlook

In this work, we have developed an effective method to directly simulate the
photodissociation followed by Coulomb explosion of methyl iodide with the
inclusion of non-adiabatic effects. Potential energy surfaces of both the
neutral molecule and the doubly charged cation are built upon high-level ab
initio electronic structure calculations. More than 40,000 trajectories are
calculated on both the neutral and ionic PESs to obtain a statistically
reliable interpretation of the delay-dependent KER signals. The non-adiabatic
effects are treated using a Landau-Zener surface hopping algorithm, which
yields a branching ratio in overall agreement with previous theories and
experiments. Our simulation of the delay-dependent KER at large pump-probe
delays shows good agreement with a recent experiment using coincident ion
momentum imaging techniques27, with the KER signals reflecting two reaction
channels. At short pump-probe delays, the simulation shows a clear bifurcation
of the KER signals near 12 fs, indicating the occurrence of non-adiabatic
transitions through a conical intersection between the ${}^{3}\text{Q}_{0}$
and ${}^{1}\text{Q}_{1}$ states. The bifurcation of the KER signals in our
simulation is statistically reliable and our simulation directly predicts
experimental KER images, which is beyond estimations using averaged AIMD
trajectories26. Moreover, metastable $\text{CH}_{3}\text{I}^{2+}$ ionic states
can be formed with small probability on one of the ionic surfaces near zero
pump-probe delay.

This work reveals that dynamics on the ionic PES plays an important role on
the delay-dependent KER signals, especially in the regime of short pump-probe
delays which cannot be described by a simple 1D model. Our work also confirms
that Coulomb explosion imaging technique can effectively be used for direct
visualization of conical intersections in coordinate space, provided
sufficient temporal resolution.

The method developed in this work is highly efficient in calculating
trajectories with non-adiabatic effects, which can be applied to molecular
dynamics simulations that highly rely on statistics. It can also be employed
to identify reaction channels with small probabilities such as roaming. In the
future, we plan to employ this method to investigate photo-isomerization
processes such as bond rearrangement in diiodomethane and ring opening in
thiophenone.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

YD thanks Sabre Kais for discussions at the early state of this study. YD
thanks Panwang Zhou and Runze Liu for helpful advice on the ab initio
calculations and non-adiabatic dynamics. The authors thank Farzaneh Ziaee for
providing the experimental data shown in Figures 4 and 5. The authors also
thank B. D. Esry for suggestions on improving the manuscript. This work is
supported by the Chemical Sciences, Geosciences, and Biosciences Division,
Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy,
Grant no. DE-FG02-86ER13491.

## References

  * Alharbi and Kais 2015 F. H. Alharbi and S. Kais, _Renewable and Sustainable Energy Reviews_ , 2015, 43, 1073–1089.
  * Bohn _et al._ 2017 J. L. Bohn, A. M. Rey and J. Ye, _Science_ , 2017, 357, 1002–1010.
  * Zhao and Pan 2022 B. Zhao and J.-W. Pan, _Chem. Soc. Rev._ , 2022, 51, 1685–1701.
  * Ding _et al._ 2016 Y. Ding, J. Pérez-Ríos and C. H. Greene, _ChemPhysChem_ , 2016, 17, 3756–3763.
  * Ding _et al._ 2017 Y. Ding, J. P. D’Incao and C. H. Greene, _Phys. Rev. A_ , 2017, 95, 022709.
  * Wu _et al._ 2016 M. Wu, S. Chen, S. Camp, K. J. Schafer and M. B. Gaarde, _J. Phys. B: At., Mol. Phys._ , 2016, 49, 062003.
  * Stolow _et al._ 2004 A. Stolow, A. E. Bragg and D. M. Neumark, _Chem. Rev._ , 2004, 104, 1719–1758.
  * Filippetto _et al._ 2022 D. Filippetto, P. Musumeci, R. K. Li, B. J. Siwick, M. R. Otto, M. Centurion and J. P. F. Nunes, _Rev. Mod. Phys._ , 2022, 94, 045004.
  * Vager _et al._ 1989 Z. Vager, R. Naaman and E. P. Kanter, _Science_ , 1989, 244, 426–431.
  * Schouder _et al._ 2022 C. A. Schouder, A. S. Chatterley, J. D. Pickering and H. Stapelfeldt, _Annual Review of Physical Chemistry_ , 2022, 73, 323–347.
  * Pitzer _et al._ 2013 M. Pitzer, M. Kunitski, A. S. Johnson, T. Jahnke, H. Sann, F. Sturm, L. P. H. Schmidt, H. Schmidt-Böcking, R. Dörner, J. Stohner, J. Kiedrowski, M. Reggelin, S. Marquardt, A. Schießer, R. Berger and M. S. Schöffler, _Science_ , 2013, 341, 1096–1100.
  * Boll _et al._ 2022 R. Boll, J. M. Schäfer, B. Richard, K. Fehre, G. Kastirke, Z. Jurek, M. S. Schöffler, M. M. Abdullah, N. Anders, T. M. Baumann, S. Eckart, B. Erk, A. De Fanis, R. Dörner, S. Grundmann, P. Grychtol, A. Hartung, M. Hofmann, M. Ilchen, L. Inhester, C. Janke, R. Jin, M. Kircher, K. Kubicek, M. Kunitski, X. Li, T. Mazza, S. Meister, N. Melzer, J. Montano, V. Music, G. Nalin, Y. Ovcharenko, C. Passow, A. Pier, N. Rennhack, J. Rist, D. E. Rivas, D. Rolles, I. Schlichting, L. P. H. Schmidt, P. Schmidt, J. Siebert, N. Strenger, D. Trabert, F. Trinter, I. Vela-Perez, R. Wagner, P. Walter, M. Weller, P. Ziolkowski, S.-K. Son, A. Rudenko, M. Meyer, R. Santra and T. Jahnke, _Nat. Phys._ , 2022, 18, 423–428.
  * Li _et al._ 2022 X. Li, A. Rudenko, M. S. Schöffler, N. Anders, T. M. Baumann, S. Eckart, B. Erk, A. De Fanis, K. Fehre, R. Dörner, L. Foucar, S. Grundmann, P. Grychtol, A. Hartung, M. Hofmann, M. Ilchen, C. Janke, G. Kastirke, M. Kircher, K. Kubicek, M. Kunitski, T. Mazza, S. Meister, N. Melzer, J. Montano, V. Music, G. Nalin, Y. Ovcharenko, C. Passow, A. Pier, N. Rennhack, J. Rist, D. E. Rivas, I. Schlichting, L. P. H. Schmidt, P. Schmidt, J. Siebert, N. Strenger, D. Trabert, F. Trinter, I. Vela-Perez, R. Wagner, P. Walter, M. Weller, P. Ziolkowski, A. Czasch, D. Rolles, M. Meyer, T. Jahnke and R. Boll, _Phys. Rev. Res._ , 2022, 4, 013029\.
  * Burt _et al._ 2018 M. Burt, K. Amini, J. W. L. Lee, L. Christiansen, R. R. Johansen, Y. Kobayashi, J. D. Pickering, C. Vallance, M. Brouard and H. Stapelfeldt, _J. Chem. Phys._ , 2018, 148, 091102.
  * Endo _et al._ 2020 T. Endo, S. P. Neville, V. Wanie, S. Beaulieu, C. Qu, J. Deschamps, P. Lassonde, B. E. Schmidt, H. Fujise, M. Fushitani, A. Hishikawa, P. L. Houston, J. M. Bowman, M. S. Schuurman, F. Légaré and H. Ibrahim, _Science_ , 2020, 370, 1072–1077.
  * Mulliken 2004 R. S. Mulliken, _J. Chem. Phys._ , 2004, 8, 382–395.
  * Riley and Wilson 1972 S. J. Riley and K. R. Wilson, _Faraday Discuss. Chem. Soc._ , 1972, 53, 132–146.
  * Shapiro and Bersohn 1980 M. Shapiro and R. Bersohn, _J. Chem. Phys._ , 1980, 73, 3810–3817.
  * Godwin _et al._ 1987 F. Godwin, P. Gorry, P. Hughes, D. Raybone, T. Watkinson and J. Whitehead, _Chem. Phys. Lett._ , 1987, 135, 163–169.
  * Eppink and Parker 1998 A. T. J. B. Eppink and D. H. Parker, _J. Chem. Phys._ , 1998, 109, 4758–4767.
  * Corrales _et al._ 2014 M. E. Corrales, V. Loriot, G. Balerdi, J. González-Vázquez, R. de Nalda, L. Bañares and A. H. Zewail, _Phys. Chem. Chem. Phys._ , 2014, 16, 8812–8818.
  * Corrales _et al._ 2012 M. E. Corrales, G. Gitzinger, J. González-Vázquez, V. Loriot, R. de Nalda and L. Bañares, _J. Phys. Chem. A_ , 2012, 116, 2669–2677.
  * Zhang _et al._ 2017 D. Zhang, S. Luo, H. Xu, M. Jin, F. Liu, B. Yan, Z. Wang, H. Liu, D. Jiang, A. Eppink, W. Roeterdink, S. Stolte and D. Ding, _Euro. Phys. J. D_ , 2017, 71, 148.
  * Amini _et al._ 2018 K. Amini, E. Savelyev, F. Brauße, N. Berrah, C. Bomme, M. Brouard, M. Burt, L. Christensen, S. Düsterer, B. Erk, H. Höppner, T. Kierspel, F. Krecinic, A. Lauer, J. W. L. Lee, M. Müller, E. Müller, T. Mullins, H. Redlin, N. Schirmel, J. Thøgersen, S. Techert, S. Toleikis, R. Treusch, S. Trippel, A. Ulmer, C. Vallance, J. Wiese, P. Johnsson, J. Küpper, A. Rudenko, A. Rouzée, H. Stapelfeldt, D. Rolles and R. Boll, _Structural Dynamics_ , 2018, 5, 014301.
  * Allum _et al._ 2018 F. Allum, M. Burt, K. Amini, R. Boll, H. Köckert, P. K. Olshin, S. Bari, C. Bomme, F. Brauße, B. Cunha de Miranda, S. Düsterer, B. Erk, M. Géléoc, R. Geneaux, A. S. Gentleman, G. Goldsztejn, R. Guillemin, D. M. P. Holland, I. Ismail, P. Johnsson, L. Journel, J. Küpper, J. Lahl, J. W. L. Lee, S. Maclot, S. R. Mackenzie, B. Manschwetus, A. S. Mereshchenko, R. Mason, J. Palaudoux, M. N. Piancastelli, F. Penent, D. Rompotis, A. Rouzée, T. Ruchon, A. Rudenko, E. Savelyev, M. Simon, N. Schirmel, H. Stapelfeldt, S. Techert, O. Travnikova, S. Trippel, J. G. Underwood, C. Vallance, J. Wiese, F. Ziaee, M. Brouard, T. Marchenko and D. Rolles, _J. Chem. Phys._ , 2018, 149, 204313.
  * Corrales _et al._ 2019 M. E. Corrales, J. González-Vázquez, R. de Nalda and L. Bañares, _J. Phys. Chem. Lett._ , 2019, 10, 138–143.
  * Ziaee _et al._ 2023 F. Ziaee, K. Borne, R. Forbes, K. R. P., Y. Malakar, B. Kaderiya, T. Severt, I. Ben-Itzhak, A. Rudenko and D. Rolles, _Phys. Chem. Chem. Phys._ , 2023, 25, 9999–10010.
  * Amatatsu _et al._ 1991 Y. Amatatsu, K. Morokuma and S. Yabushita, _J. Chem. Phys._ , 1991, 94, 4858–4876.
  * Amatatsu _et al._ 1996 Y. Amatatsu, S. Yabushita and K. Morokuma, _J. Chem. Phys._ , 1996, 104, 9783–9794.
  * Alekseyev _et al._ 2007 A. B. Alekseyev, H.-P. Liebermann, R. J. Buenker and S. N. Yurchenko, _J. Chem. Phys._ , 2007, 126, 234102.
  * Kamiya and Taketsugu 2019 M. Kamiya and T. Taketsugu, _J. Comput. Chem._ , 2019, 40, 456–463.
  * Wang _et al._ 2019 H. Wang, M. Odelius and D. Prendergast, _J. Chem. Phys._ , 2019, 151, 124106.
  * Guo 1992 H. Guo, _J. Chem. Phys._ , 1992, 96, 6629–6642.
  * Xie _et al._ 2000 D. Xie, H. Guo, Y. Amatatsu and R. Kosloff, _J. Phys. Chem. A_ , 2000, 104, 1009–1019.
  * Evenhuis and Manthe 2011 C. R. Evenhuis and U. Manthe, _J. Phys. Chem. A_ , 2011, 115, 5992–6001.
  * Werner and Knowles 1988 H. Werner and P. J. Knowles, _J. Chem. Phys._ , 1988, 89, 5803–5814.
  * Knowles and Werner 1988 P. J. Knowles and H.-J. Werner, _Chem. Phys. Lett._ , 1988, 145, 514–522.
  * Werner and Knowles 1985 H. Werner and P. J. Knowles, _J. Chem. Phys._ , 1985, 82, 5053–5063.
  * Knowles and Werner 1985 P. J. Knowles and H.-J. Werner, _Chem. Phys. Lett._ , 1985, 115, 259–267.
  * Kreplin _et al._ 2019 D. A. Kreplin, P. J. Knowles and H.-J. Werner, _J. Chem. Phys._ , 2019, 150, 194106.
  * Dunning 1989 J. Dunning, Thom H., _J. Chem. Phys._ , 1989, 90, 1007–1023.
  * Peterson _et al._ 2003 K. A. Peterson, D. Figgen, E. Goll, H. Stoll and M. Dolg, _J. Chem. Phys._ , 2003, 119, 11113–11123.
  * Peterson _et al._ 2006 K. A. Peterson, B. C. Shepler, D. Figgen and H. Stoll, _J. Phys. Chem. A_ , 2006, 110, 13877–13883.
  * Werner _et al._ 2020 H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, I. Miller, Thomas F., A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut and M. Sibaev, _J. Chem. Phys._ , 2020, 152, 144107.
  * Werner _et al._ 2012 H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, _WIREs Computational Molecular Science_ , 2012, 2, 242–253.
  * Alekseyev _et al._ 2007 A. B. Alekseyev, H.-P. Liebermann and R. J. Buenker, _J. Chem. Phys._ , 2007, 126, 234103.
  * Xie and Domcke 2017 W. Xie and W. Domcke, _J. Chem. Phys._ , 2017, 147, 184114.
  * Xie _et al._ 2019 W. Xie, M. Sapunar, N. Došlić, M. Sala and W. Domcke, _J. Chem. Phys._ , 2019, 150, 154119.
  * Belyaev _et al._ 2014 A. K. Belyaev, C. Lasser and G. Trigila, _J. Chem. Phys._ , 2014, 140, 224108.
  * Suchan _et al._ 2020 J. Suchan, J. Janoš and P. Slavíček, _J. Chem. Theo. Comput._ , 2020, 16, 5809–5820.
  * Tully 1990 J. C. Tully, _J. Chem. Phys._ , 1990, 93, 1061–1071.
  * Belyaev and Lebedev 2011 A. K. Belyaev and O. V. Lebedev, _Phys. Rev. A_ , 2011, 84, 014701\.
  * Pernpointner _et al._ 2012 M. Pernpointner, J. P. Zobel, E. Fasshauer and A. N. Sil, _Chem. Phys._ , 2012, 407, 39–45.
  * Hult Roos _et al._ 2017 A. Hult Roos, J. Eland, D. Koulentianos, R. Squibb, L. Karlsson and R. Feifel, _Chem. Phys._ , 2017, 491, 42–47.
  * Chang _et al._ 2021 K. F. Chang, H. Wang, S. M. Poullain, D. Prendergast, D. M. Neumark and S. R. Leone, _J. Chem. Phys._ , 2021, 154, 234301.
  * Attar _et al._ 2015 A. R. Attar, A. Bhattacherjee and S. R. Leone, _J. Phys. Chem. Lett._ , 2015, 6, 5072–5077.
  * Riley and Wilson 1972 S. J. Riley and K. R. Wilson, _Faraday Discuss. Chem. Soc._ , 1972, 53, 132–146.
  * Hess _et al._ 1986 W. P. Hess, S. J. Kohler, H. K. Haugen and S. R. Leone, _J. Chem. Phys._ , 1986, 84, 2143–2149.
  * Baklanov _et al._ 2000 A. V. Baklanov, M. Aldener, B. Lindgren and U. Sassenberg, _Chem. Phys. Lett._ , 2000, 325, 399–404.


# Weak Memory Demands Model-based Compiler Testing

Luke Geeson<EMAIL_ADDRESS>University College London, UK

###### Abstract.

A compiler bug arises if the behaviour of a compiled concurrent program, as
allowed by its architecture memory model, is not a behaviour permitted by the
source program under its source model. One might reasonably think that most
compiler bugs have been found in the decade since the introduction of the
C/C++ memory model. We observe that processor implementations are increasingly
exploiting the behaviour of relaxed architecture models. As such, compiled
programs may exhibit bugs not seen on older hardware. To account for this we
require model-based compiler testing.

While this observation is not surprising, its implications are broad.
Compilers and their testing tools will need to be updated to follow hardware
relaxations, concurrent test generators will need to be improved, and
assumptions of prior work will need revisiting. We explore these ideas using a
compiler toolchain bug we reported in LLVM.

††booktitle: The Future of Weak Memory, January 15, 2024

##  1\. Introduction

Consider the Message Passing litmus test in Fig. 1. The outcome {P1:r0=0; y=2}
is forbidden by the C/C++ memory model (C11, ) (Fig. 2). If Fig. 1 is compiled
using LLVM (Fig. 3), the outcome {P1:r0=0; y=2} is allowed (Fig. 4) under the
Arm AArch64 model (Seal:2000:AAR:517257, ). We found (swp, ) this compiler bug
after updates (rlxswp, ) to the herd memory model tool.

⬇ { *x = 0; *y = 0; } // fixed initial state #define relaxed
memory_order_relaxed #define release memory_order_release #define acquire
memory_order_acquire P0 (atomic_int* y,atomic_int* x) {
atomic_store_explicit(x,1,relaxed); atomic_thread_fence(release);
atomic_store_explicit(y,1,relaxed); } P1 (atomic_int* y,atomic_int* x) {
atomic_exchange_explicit(y,2,release); atomic_thread_fence(acquire); int r0 =
atomic_load_explicit (x,relaxed); } exists (P1:r0=0 /\ y=2)
---
Figure 1. Message Passing Litmus Test. We reported (swp, ) this bug in LLVM’s
implementation of atomic_exchange. We found this using the Téléchat compiler
testing tool.

The compiled counterpart of Fig. 1 did not exhibit this outcome under the Arm
AArch64 (aarch64, ) model until the herd tool was updated (rlxswp, ) to
reflect the intent of the architecture specification (Seal:2000:AAR:517257, ).
Likewise, older hardware that does not implement the relaxations allowed by
the architecture may not exhibit the behaviour. Consequently compiler testing
tools that depend on old hardware miss this bug. This motivates a need for
compiler testing techniques that are parameterised over source and
architecture models. In CGO 2024 (geeson2023compiler, ), we will present the
Téléchat tool that does exactly this. We have found several concurrency bugs
using Téléchat, numerous tooling bugs, and as far as we know it the first tool
of its kind to be deployed in industry automated testing of a production
compiler (It is the first to automatically test Arm Compiler). Compilers and
their testing tools need to account for future relaxations, as they occur.

Fig. 1 is not generated by today’s litmus test generators. Fig. 1 is a special
kind of message passing test that discards the read of the atomic_exchange and
checks y instead. Historically message passing checks whether the read of the
atomic_exchange can be reordered past the read into P1:r0. Test generators,
such as diy (Alglave:2012:FWM:2205506.2205523, ), currently implement the
historic case. Bugs similar to Fig. 1 have thus only been found by manual
inspection, until now. We conclude that improvements to test generators are
needed.

Interestingly, this bug disappears under direct observation. Using thread-
local reads int r1 = atomic_exchange... (P1 is an Observer \- see §B.2
(Seal:2000:AAR:517257, )) makes the bug disappear. The exchange is compiled to
an Arm SWP instruction and the LLVM dead register optimisation (LLVMdrd, )
rewrites the SWP destination register to the zero register (WZR), since the
read of the exchange is unused (Fig. 3). Unfortunately instructions, whose
destination register is a zero register (WZR/XZR), are not regarded as doing a
read for the purpose of a acquire barrier (§C3.2.12 of (Seal:2000:AAR:517257,
)). LLVM has introduced a bug (Fig. 4). This optimisation was allowed,
following an assumption (cmmtest, ) that “optimisations affecting only the
thread-local state cannot induce concurrency compiler bugs”. Fig. 1 challenges
this assumption with a new kind of Heisenbug. We conclude that assumptions
surrounding compiler testing with models may need to be revisited in light of
hardware relaxations.

⬇ { P1:r0=0; [y]=1; } { P1:r0=1; [y]=1; } { P1:r0=1; [y]=2; }
---
Figure 2. Outcomes of running Fig. 1 under the C/C++ memory model (C11, ). Note the absense of {P1:r0=0; y=2}. ⬇ AArch64 test { [P1_r0]=0;[x]=0;[y]=0; uint64_t %P0_x=x;uint64_t %P0_y=y; uint64_t %P1_P1_r0=P1_r0; uint64_t %P1_x=x;uint64_t %P1_y=y } P0 | P1 ; MOV W9,#1 | MOV W9,#2 ; STR W9,[X%P0_x] | SWPL W9, WZR,[X%P1_y]; DMB ISH | DMB ISHLD ; STR W9,[X%P0_y] | LDR W8,[X%P1_x] ; | STR W8,[X%P1_P1_r0] ; | ; exists (P1_r0=0 /\ [y]=2)
---
Figure 3. Compiled version (swp, ) of Fig. 1. Note the destination register of
the SWPL is the zero register (WZR).  ⬇ { P1:r0=0; [y]=2; } bug! { P1:r0=0;
[y]=1; } { P1:r0=1; [y]=1; } { P1:r0=1; [y]=2; }
---
Figure 4. Outcomes of running Fig. 3 under the Arm AArch64 model (aarch64, ). Note the presence of {P1:r0=0; y=2}. ⬇ AArch64 test { [P1_r0]=0;[x]=0;[y]=0; uint64_t %P0_x=x;uint64_t %P0_y=y; uint64_t %P1_P1_r0=P1_r0; uint64_t %P1_x=x;uint64_t %P1_y=y } P0 | P1 ; MOV W9,#1 | MOV W9,#2 ; STR W9,[X%P0_x] | SWPL W9, W15,[X%P1_y]; DMB ISH | DMB ISHLD ; STR W9,[X%P0_y] | LDR W8,[X%P1_x] ; RET | STR W8,[X%P1_P1_r0] ; | RET ; exists (P1_r0=0 /\ [y]=2)
---
Figure 5. Compiled version of Fig. 1 (fixed). Note the destination register of
the SWPL is not the zero register WZR, but rather register W15. We recommended
this fix in our bug report (swp, ). ⬇ { P1:r0=0; [y]=1; } { P1:r0=1; [y]=1; }
{ P1:r0=1; [y]=2; }
---
Figure 6. Outcomes of running Fig. 5 under the Arm AArch64 model (aarch64, ).
Note the absence of {P1:r0=0; y=2}.

## Acknowledgment

We thank supervisors James Brotherston, Earl Barr, and Lee Smith. Ana Farinha,
Richard Grisenthwaite, Alastair Donaldson, John Wickerson, Tyler Sorensen,
Wilco Dijkstra, Tomas Matheson, Arm’s Compiler Teams and Arm Architecture &
Technology Group for their feedback and assistance. This work was supported by
the Engineering and Physical Sciences Research Council [grant number
EP/V519625/1]. The views of the authors expressed in this paper are not
endorsed by Arm or any other company mentioned.

## References

  * (1) Alglave, J., Maranget, L., Sarkar, S., and Sewell, P. Fences in Weak Memory Models (Extended Version).  Form. Methods Syst. Des. 40, 2 (Apr. 2012), 170–205.
  * (2) Arm-Limited.  Arm Architecture Reference Manual.  Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 2023.
  * (3) Deacon, W., and Alglave, J. Armv8 AArch64 Memory Model.  https://github.com/herd/herdtools7/blob/master/herd/libdir/aarch64.cat, 2021\.
  * (4) Geeson, L. [AArch64]: Atomic Exchange Allows Reordering past Acquire Fence .  https://github.com/llvm/llvm-project/issues/68428, 2023.
  * (5) Geeson, L., and Smith, L. Compiler testing with relaxed memory models.  https://arxiv.org/abs/2310.12337, 2023.  To appear In Proceedings of the 2024 International Symposium on Code Generation and Optimization.
  * (6) ISO-C-Std, O. ISO/IEC 9899:201x.  https://www.open-std.org/jtc1/sc22/wg14/www/docs/n2912.pdf, 2022\.  Accessed: 2023-11-10.
  * (7) LLVM.  AArch64 Dead register definitions.  https://llvm.org/docs/doxygen/AArch64DeadRegisterDefinitionsPass_8cpp_source.html.
  * (8) Morisset, R., Pawan, P., and Zappa Nardelli, F. Compiler Testing via a Theory of Sound Optimisations in the C11/C++11 Memory Model.  In Proceedings of the 34th ACM SIGPLAN Conference on Programming Language Design and Implementation (New York, NY, USA, 2013), PLDI ’13, ACM, pp. 187–196.
  * (9) Nikoleris, N. Restore architectural intent for aarch64 swp/cas that return to zr.  https://github.com/herd/herdtools7/pull/683.