1]Department of Chemistry, Stanford University, Stanford, CA, USA

2]Stanford PULSE Institute, SLAC National Accelerator Laboratory, Menlo Park,
CA, USA

3]Department of Chemistry, Norwegian University of Science and Technology,
Trondheim, 7491, Norway

4]Deutsches Elektronen-Synchrotron DESY, Hamburg, Germany

5]Institute of Physical Chemistry, University of Hamburg, Hamburg, Germany

# Unexpected hydrogen dissociation in thymine: predictions from a novel
coupled cluster theory

Eirik F. Kjønstad<EMAIL_ADDRESS>O. Jonathan Fajen  Alexander C. Paul
Sara Angelico  Dennis Mayer  Markus Gühr  Thomas J. A. Wolf  Todd J. Martínez
<EMAIL_ADDRESS>Henrik Koch<EMAIL_ADDRESS>[  [  [  [  [

###### Abstract

The fate of thymine upon excitation by ultraviolet radiation has been the
subject of intense debate over the past three decades. Today, it is widely
believed that its ultrafast excited state decay stems from a radiationless
transition from the bright $\pi\pi^{\ast}$ state to a dark $n\pi^{\ast}$
state. However, conflicting theoretical predictions have made the experimental
data difficult to interpret. Here we simulate the ultrafast dynamics in
thymine at the highest level of theory to date, performing wavepacket dynamics
with a new coupled cluster method. Our simulation confirms an ultrafast
$\pi\pi^{\ast}$ to $n\pi^{\ast}$ transition ($\tau=$41\pm
14\text{\,}\mathrm{fs}$$). Furthermore, the predicted oxygen-edge X-ray
absorption spectra agree quantitatively with the experimental results. Our
simulation also predicts an as-yet uncharacterized photochemical pathway: a
$\pi\sigma^{\ast}$ channel that leads to hydrogen dissociation at one of the
two N-H bonds in thymine. Similar behavior has been identified in other
heteroaromatic compounds, including adenine, and several authors have
speculated that a similar pathway may exist in thymine. However, this was
never confirmed theoretically or experimentally. This prediction calls for
renewed efforts to experimentally identify or exclude the presence of this
channel.

## Introduction

Thymine, like other nucleobases, undergoes ultrafast radiationless relaxation
back to the ground state after being excited by ultraviolet radiation. This
property has been tied to the resilience of genetic material against light-
induced damage [1]. However, the exact mechanism of this decay is not fully
understood and has been a subject of debate for several decades. Experiments
have identified at least two excited state decay channels, one with a lifetime
on the order of $\lesssim$$100\text{\,}\mathrm{fs}$, and one considerably
longer, on the order of several ps [2, 3, 4, 5]. Yet, the underlying
mechanisms have been challenging to discern, with proposed explanations
necessarily relying on simulations of the molecular dynamics. These
simulations, in turn, introduce approximations with errors that are difficult
to control. Different theoretical methods have therefore produced different
explanations, and a consensus has yet to emerge.

Most reported simulations implicate two low-lying excited states in the
relaxation: a dark $n\pi^{\ast}$ state (S1) and the bright $\pi\pi^{\ast}$
state (S2) into which the system is initially excited. Several simulations
predict a $\pi\pi^{\ast}$ trapping channel in which the initial excitation to
the $\pi\pi^{\ast}$ state is rapidly followed by relaxation into a minimum on
the $\pi\pi^{\ast}$ surface, where the wavepacket is trapped for tens or
hundreds of fs [6, 7, 8, 9, 10, 11, 12]. These simulations disagree, however,
on the amount of $\pi\pi^{\ast}$ trapping, as well as the timescale and nature
of the subsequent processes, with proposed mechanisms including
$\pi\pi^{\ast}$ to $n\pi^{\ast}$ relaxation [11] and direct $\pi\pi^{\ast}$
relaxation to the ground state [8]. Some of these studies indicate that the
amount of $\pi\pi^{\ast}$ trapping is reduced by improving the description of
dynamical correlation (instantaneous electron-electron interactions) [11, 8],
a pattern that has also been found in the closely related nucleobase uracil
[13]. In line with this, an early density functional theory (DFT) study found
significant $n\pi^{\ast}$ population within the first $50\text{\,}\mathrm{fs}$
[14] and a more recent mixed-reference spin-flip DFT study also found rapid
barrier-less $\pi\pi^{\ast}/n\pi^{\ast}$ transfer
($\tau=$30\text{\,}\mathrm{fs}$$), including subsequent $n\pi^{\ast}$ trapping
[15].

Experimental evidence has implicated the $n\pi^{\ast}$ state in the early
dynamics. Indeed, by determining the gas phase oxygen-edge time-resolved X-ray
absorption spectrum, Wolf et al. [16] found a fast component ($\tau=$60\pm
30\text{\,}\mathrm{fs}$$) which was attributed to population of the
$n\pi^{\ast}$ state. Thus, the wavepacket appears to already transfer some of
its population to the $n\pi^{\ast}$ state within the first
$100\text{\,}\mathrm{fs}$. Furthermore, Wolf et al. found that the
$n\pi^{\ast}$ signature lasts for several ps, revealing a second relaxation
mechanism. After passing through the $\pi\pi^{\ast}/n\pi^{\ast}$ conical
intersection, parts of the wavepacket appears to get trapped in a minimum on
the $n\pi^{\ast}$ surface. While experiments have shown that the $n\pi^{\ast}$
state is involved in the early sub-100 fs dynamics, it remains an open
question whether or not there is some trapping in the $\pi\pi^{\ast}$ state
[17, 18]. More accurate simulations of the dynamics are therefore essential to
unravel the precise relaxation mechanism in thymine.

Here we present the highest-level wavepacket simulation on thymine to date. To
the best of our knowledge, this is also the simulation with the highest level
of electronic structure theory performed on a molecular system of this size.
Thanks to recent developments [19, 20, 21, 22, 23, 24], we were able to
describe the electronic structure with the highly accurate coupled cluster
singles and doubles (CCSD) [25] method. This is the first time that this
method is applied in nonadiabatic dynamics. Coupled cluster (CC) theory is
well-known for effectively capturing dynamical correlation, but it has been
widely regarded as unsuited for excited state dynamics due to the presence of
numerical artifacts at conical intersections [26, 27, 19, 28]. Recent work has
shown that these problems can be removed with similarity constrained CC (SCC)
theory [20, 21, 24]. Here, we apply the SCC with singles and doubles (SCCSD)
method to perform molecular dynamics simulations. To our knowledge, the
present study is the first nonadiabatic dynamics simulation with a coupled
cluster theory that describes conical intersections correctly. This
demonstrates that the method is a viable electronic structure method for
simulating photochemical phenomena.

Figure 1: Photochemical pathways in the thymine simulation. Following
photoexcitation to the bright $S_{2}$($\pi\pi^{\ast}$) state, the simulation
predicts two channels. The main channel (87%) is the $n\pi^{\ast}$ trapping
channel. Here, the wavepacket passes through the $S_{1}$/$S_{2}$ intersection
and heads toward a minimum on the $S_{1}(n\pi^{\ast})$ surface. This minimum
is reached in two ways, either by heading to the minimum directly (solid line)
or by reaching it indirectly through a $\pi\pi^{\ast}$ region on $S_{1}$
(dashed line). The second channel (13%) is an N-H dissociation channel. Here,
the wavepacket acquires $\pi\sigma^{\ast}$ character on $S_{2}$ (with the
anti-bonding $\sigma^{\ast}$ orbital located at an N-H bond) before it passes
through the $S_{1}$/$S_{2}$ intersection while retaining its
$\pi\sigma^{\ast}$ character when it transfers to the $S_{1}$ surface. Once
the wavepacket is on $S_{1}$, N-H dissociation quickly follows as the
wavepacket heads towards an intersection with the electronic ground state.

The main photochemical pathways in our simulation are illustrated in Figure 1.
Our simulation confirms an ultrafast $\pi\pi^{\ast}$ to $n\pi^{\ast}$
conversion. In particular, the simulation predicts a time-resolved X-ray
absorption spectrum (see Figure 2) that captures the ultrafast component of
the time-resolved oxygen-edge spectrum observed experimentally by Wolf et al.
[16]. We show that this component is unequivocally associated with the
$n\pi^{\ast}$ state, demonstrating that the ultrafast component of the decay
is due to internal conversion from the $\pi\pi^{\ast}$ state. We further find
no significant $\pi\pi^{\ast}$ trapping.

Figure 2: Thymine oxygen-edge X-ray absorption spectrum. Panel A illustrates
the pump-probe scheme for thymine, including the numbering used for the
important atoms in the dynamics. Panels B, C, and D show oxygen-edge X-ray
absorption spectra for thymine, simulated (B and D) and experimental (C) [16].
Simulated spectra were computed using CC3/cc-pVDZ transition strengths and
AIMS dynamics with CCSD/cc-pVDZ and SCCSD/cc-pVDZ, together with Gaussian
broadening, $0.3\text{\,}\mathrm{eV}$/$10\text{\,}\mathrm{fs}$ FWHM (panel B)
and $0.3\text{\,}\mathrm{eV}$/$100\text{\,}\mathrm{fs}$ FWHM (panel D). For
more details, see Supporting Information S2. The most prominent feature in the
spectrum is a bright signal at $526\text{\,}\mathrm{eV}$, which is due to the
$n\pi^{\ast}$ state and reflects the rapid $\pi\pi^{\ast}$/$n\pi^{\ast}$
internal conversion after photoexcitation to the $\pi\pi^{\ast}$ state (at
time $t=0$). The oscillations in the $n\pi^{\ast}$ feature (see panel B) are
washed out by a larger time broadening (see panel D), reflecting the time
resolution of the experiment (see panel C).

The interplay of a bright $\pi\pi^{*}$ and a dark $n\pi^{*}$ state is also
central to the excited-state relaxation of a wide variety of other
chromophores, including nucleobases such as uracil and adenine, although
lifetimes and branching ratios differ. In adenine, another decay mechanism has
been proposed, one in which a dark $\pi\sigma^{*}$ state is accessed through
interconversion from the $\pi\pi^{*}$ state [29, 30]. More generally, the
importance of these low-lying $\pi\sigma^{*}$ states are well-recognized in a
wide variety of heteroaromatic systems [31, 32, 33, 34]. In adenine, the
$\pi\sigma^{*}$ state is dissociative along the N-H bond coordinate and
intersects with the ground state ($S_{0}$) at larger N-H distances. Following
excitation to the $\pi\pi^{*}$ state, the wavepacket may reach this
$\pi\pi^{*}/\pi\sigma^{*}$ conical intersection within
$20\text{\,}\mathrm{fs}$. Once on the repulsive $\pi\sigma^{*}$ state, N-H
dissociation proceeds rapidly, resulting in a characteristically anisotropic
total kinetic energy release (TKER) spectrum of H atom fragments [35].

Our high-level dynamics simulation predicts that thymine also has a similar
$\pi\sigma^{*}$-mediated N-H dissociation channel (see Figure 1). This finding
comes as a surprise because no such dissociative channel has been previously
reported. In fact, TKER spectra across a broad range of wavelengths, 270–230
nm, showed no signature of ultrafast N-H dissociation [36], and recent
experimental and theoretical investigations have not invoked this pathway to
explain the molecular dynamics. Nonetheless, in analogy to the photochemistry
of adenine and other heteroaromatics, it may be that certain regions of the
$\pi\pi^{*}$ surface of thymine—in particular, regions accessed by a highly
excited N-H stretching mode—lead to rapid relaxation through an intersection
with a $\pi\sigma^{*}$ state, followed by hydrogen dissociation at the N-H
bond.

## Results and discussion

Figure 1 illustrates the two main photochemical pathways observed in the
simulation. The first pathway is the $n\pi^{\ast}$ trapping channel (87%).
Here, the wavepacket travels directly from the Franck-Condon region to the
$S_{2}$($\pi\pi^{\ast}$)/$S_{1}$($n\pi^{\ast}$) conical intersection, where it
transfers population to the $S_{1}$($n\pi^{\ast}$) state. The transferred
population then gets trapped in a minimum on the $S_{1}$($n\pi^{\ast}$)
surface. The second pathway is a dissociative $\pi\sigma^{\ast}$ channel
(13%). Here, the wavepacket acquires $\pi\sigma^{\ast}$ character on $S_{2}$,
where the anti-bonding $\sigma^{\ast}$ orbital is localized at an N–H bond,
leading to rapid hydrogen dissociation.

Figure 3: Dynamics with CCSD (panels A and B) and SCCSD (panels C and D).
Panels A and C show a branching or $gh$ plane for an initial condition that
encounters numerical artifacts in the simulations using CCSD, where $\Delta E$
is the energy difference between $S_{1}$ and $S_{2}$. For CCSD (A), there is
an ellipse of degeneracy in the $gh$ plane with unphysical complex-valued
energies in the interior of the elliptical boundary; for SCCSD (C), we instead
see a single point of degeneracy and no unphysical energies. In CCSD
simulations, the wavepacket may approach the intersection too closely and end
up in the region with complex-valued energies (the interior of the ellipse in
the $gh$ plane). Whenever this happens, we re-run the simulation with SCCSD.
Panels B and D show the corresponding potential energy curves for the
simulation with CCSD (B) and SCCSD (D). At $71\text{\,}\mathrm{fs}$, the CCSD
simulation enters the complex-valued region and is terminated (B). The SCCSD
simulation, on the other hand, does not encounter any problems (D). The
conical intersection point used in the branching plane calculation is the
geometry with the smallest $\Delta E$ (as given by SCCSD) in the nonadiabatic
event at $71\text{\,}\mathrm{fs}$.

The $n\pi^{\ast}$ trapping channel can be experimentally identified from the
gas-phase time-resolved oxygen-edge X-ray absorption spectrum. In Figure 2, we
compare the simulated and experimental [16] X-ray absorption spectra. The most
striking feature in these spectra is a bright signal, at around
$526\text{\,}\mathrm{eV}$, that grows in intensity in the first
50–$70\text{\,}\mathrm{fs}$. The theoretical and experimental spectra are in
agreement for this feature (see panels C and D), which we find to be due to
population transfer from $S_{2}$($\pi\pi^{\ast}$) to $S_{1}$($n\pi^{\ast}$),
in agreement with the mechanism proposed by Wolf et al. [16]. By analyzing the
simulation data, we find a characteristic time of $\tau=$41\pm
14\text{\,}\mathrm{fs}$$ for this conversion (see Supporting Information S1),
which fits well with experimental estimates ($60\pm 30\text{\,}\mathrm{fs}$
[16] and $39\pm 1\text{\,}\mathrm{fs}$ [37]). In the simulated spectrum, we
also find weak but visible features that are associated with the initial
$\pi\pi^{\ast}$ dynamics: in the first $20\text{\,}\mathrm{fs}$, there are
broad and diffuse features at around 526–528 eV and at
$534\text{\,}\mathrm{eV}$, where the $534\text{\,}\mathrm{eV}$ feature is
partially hidden by the ground state bleach (see Figure 2B). These features
reflect the rapid movement of the wavepacket away from the Franck-Condon
region and towards the $S_{2}$($\pi\pi^{\ast}$)/$S_{1}$($n\pi^{\ast}$) conical
intersection. Moreover, we find frequency oscillations associated with
dynamics on the $n\pi^{\ast}$ state. These oscillations (at
$526\text{\,}\mathrm{eV}$) have a period of about $20\text{\,}\mathrm{fs}$ and
an amplitude of about $1\text{\,}\mathrm{eV}$ (see panel B). The limited
($100\text{\,}\mathrm{fs}$ FWHM time resolution of the experiment [16] washes
these oscillations out (see Figure 2C), but we predict that they would be
observed with improved time resolution. The oscillations are similarly washed
out when we apply the same time broadening to the simulated spectrum (see
Figure 2D).

An analysis of the simulation data, and in particular the growth of the
$n\pi^{\ast}$ signal in the simulated spectrum, yields a
$\pi\pi^{\ast}$/$n\pi^{\ast}$ conversion time of $\tau=$41\pm
14\text{\,}\mathrm{fs}$$. This is consistent with the rate of
$\pi\pi^{\ast}$/$n\pi^{\ast}$ conversion in the simulated dynamics, that is,
from the observed change in character of the wavepacket from $\pi\pi^{\ast}$
to $n\pi^{\ast}$. Analysis of the electronic character of the wavepacket
reveals that its character does not always correspond to what one would expect
from the population on $S_{1}$ and $S_{2}$, assuming that these correspond to
$n\pi^{\ast}$ and $\pi\pi^{\ast}$, respectively. The adiabats mostly retain
the orbital character in the Franck-Condon region, with $S_{1}$ mostly of
$n\pi^{\ast}$ character and $S_{2}$ mostly of $\pi\pi^{\ast}$ character.
However, the simulated spectrum shows that the adiabatic populations
overestimate the true rate of internal conversion through the
$\pi\pi^{\ast}$/$n\pi^{\ast}$ intersection; in particular, the adiabatic
population transfer to $S_{1}$ is faster $(\tau=$17\pm
1\text{\,}\mathrm{fs}$)$ than the growth of the $n\pi^{\ast}$ signal in the
spectrum $(\tau=$41\pm 14\text{\,}\mathrm{fs}$)$; see Supporting Information
S1. This overestimation has two distinct causes. The first is that a part of
the wavepacket that is transferred to $S_{1}$ retains its $\pi\pi^{\ast}$
character for some period of time, before eventually acquiring $n\pi^{\ast}$
character and making its way towards the $n\pi^{\ast}$ minimum. This slows
down the appearance of the $n\pi^{\ast}$ signal compared to that expected from
the adiabatic populations. The second is the $\pi\sigma^{\ast}$ character
associated with the N1-H dissociation channel. Here, the wavepacket that
initially has $\pi\pi^{\ast}$ character acquires $\pi\sigma^{\ast}$ character
on $S_{2}$ and retains it after conversion to $S_{1}$ (that is, $n\pi^{\ast}$
and $\pi\sigma^{\ast}$ flip energetic ordering). This conversion
($\pi\pi^{\ast}\rightarrow\pi\sigma^{\ast}$) contributes to the adiabatic
$S_{2}\rightarrow S_{1}$ population transfer but not to the
$\pi\pi^{\ast}\rightarrow n\pi^{\ast}$ population transfer. Indeed, when we
reassign the populations according to the electronic character of the
wavepacket, we find that the simulated spectrum is in close agreement with the
diabatic populations $(\tau=37\pm 9$), thus showing that the
$526\text{\,}\mathrm{eV}$ signal is due to the $n\pi^{\ast}$ character of the
wavepacket.

To better understand the dynamics, and in particular the $n\pi^{\ast}$
trapping channel, we identify stationary points on the $S_{1}$($n\pi^{\ast}$)
and $S_{2}$($\pi\pi^{\ast}$) surfaces in the vicinity of the Franck-Condon
region. In agreement with previous calculations, we find a minimum on the
$S_{1}$($n\pi^{\ast}$) surface at an extended C4–O8 bond length [11, 16]. This
extension is due to the anti-bonding character of the $\pi^{\ast}$ orbital
along the bond. We also find an
$S_{1}$($n\pi^{\ast}$)/$S_{2}$($\pi\pi^{\ast}$) minimum energy conical
intersection (MECI) that can be reached from the Franck-Condon region through
C5–C6 elongation. Close to the intersection, CCSD exhibits numerical artifacts
which can be removed with SCCSD. Figure 3 shows one of these artifacts as well
as how it is corrected.

By analyzing stationary points on $S_{1}$($n\pi^{\ast}$) and
$S_{2}$($\pi\pi^{\ast}$), Wolf et al. [16] suggested that the excited state
decay of thymine follows a two-step process in the C5-C6 and C4-O8 coordinates
(see Figure 2A for atom labelling): following photoexcitation, the C5–C6 bond
is first elongated, and along this stretching coordinate, the
$S_{2}$($\pi\pi^{\ast}$)/$S_{1}$($n\pi^{\ast}$) intersection seam is
accessible; then, after interconversion to the $S_{1}$($n\pi^{\ast}$) state,
the C4–O8 bond elongates as the wavepacket heads towards the
$S_{1}$($n\pi^{\ast}$) minimum. This picture is borne out by our dynamics
simulation. In Figure 4, we show the time evolution of the nuclear density in
the C5–C6 and C4–O8 bond coordinates, and we indeed see this two-step process
unfolding in real time.

The most intriguing result in the dynamics simulation is the presence of a
hydrogen dissociation channel (13%). Figure 5 characterizes this channel in
terms of the normal mode $Q$ associated with the N1-H stretch (see panel A).
When thymine is displaced along $Q$ away from the equilibrium geometry
($Q=0$), a dissociative $\pi\sigma^{\ast}$ state comes down and intersects
with $n\pi^{\ast}$ and $\pi\pi^{\ast}$ states at around $Q=0.4$ (see panel B),
which corresponds to an N1-H bond length of around $1.3$ Å. At longer N1-H
bond lengths, the character of $S_{1}$ is $\pi\sigma^{\ast}$ (see panel C).
The two dissociating initial conditions both have a short initial N1-H bond
length, and they both transfer population to $S_{1}$ at long N1-H bond lengths
(see panel D). The simplest explanation for the dissociative channel,
therefore, is that initial conditions with short N1-H bonds experience a rapid
extension of the bond in the initial phase of the nuclear dynamics, allowing
them to access the $\pi\sigma^{\ast}$ state. This suggests that initial
conditions with shorter N1-H bonds at $t=0$ would be more likely to
dissociate. This indeed appears to be the case: out of 17 additional initial
conditions specifically selected to have N1-H bond lengths shorter than $0.9$
Å, we find that 7 of the conditions (41%) have access to parts of the
potential energy surfaces with $\pi\sigma^{\ast}$ character (see Supporting
Information S6).

Figure 4: Nuclear density snapshots in the C5-C6 and C4-O8 coordinates. Three
important stationary points are shown: the ground state minimum
$S_{0}^{\mathrm{min}}$, the minimum energy conical intersection $S_{1}/S_{2}$
MECI, and the $S_{1}$ minimum $n\pi^{\ast}_{\mathrm{min}}$. The wavepacket
quickly moves away from the Franck-Condon region ($S_{0}^{\mathrm{min}}$) and
towards the minimum-energy conical intersection, where it starts transferring
population to the $n\pi^{\ast}$ surface ($19\text{\,}\mathrm{fs}$), eventually
causing the wavepacket to split ($38\text{\,}\mathrm{fs}$). At longer times,
after almost all of its population has transferred to the lower surface
($58\text{\,}\mathrm{fs}$), the wavepacket settles in the vicinity of a
minimum on the $n\pi^{\ast}$ surface (77 and $96\text{\,}\mathrm{fs}$).

In the dynamics simulation, $2$ out of $16$ initial conditions (13%) lead to
hydrogen dissociation at the N1-H bond. Inspection of the natural transition
orbitals (NTOs) of one of the conditions shows how the dissociation happens
(see Supporting Information S6). The other initial condition has similar
behavior. Already at $2.5\text{\,}\mathrm{fs}$, $S_{3}$ is mainly of
$\pi\sigma^{\ast}$ character, while $S_{2}$ is still of $\pi\pi^{\ast}$
character and $S_{1}$ of $n\pi^{\ast}$ character. As the wavepacket moves
along the N1-H stretching coordinate, the $\pi\sigma^{\ast}$ state is
stabilized, ultimately falling below the $\pi\pi^{\ast}$ and $n\pi^{\ast}$
states. The $S_{3}$/$S_{2}$ gap is about $0.5\text{\,}\mathrm{eV}$ at
$2.4\text{\,}\mathrm{fs}$, but by $3.4\text{\,}\mathrm{fs}$, this gap has
decreased to $0.1\text{\,}\mathrm{eV}$. At the same time, $S_{3}$ is of mainly
$\pi\sigma^{\ast}$ character at $2.4\text{\,}\mathrm{fs}$, but by
$4.0\text{\,}\mathrm{fs}$, the character of $S_{2}$ and $S_{3}$ has flipped.
Eventually, at around $4.1\text{\,}\mathrm{fs}$, there is an intersection
between $S_{1}$ and $S_{2}$, where most of the population is transferred to
$S_{1}$. Moving away from the intersection, $S_{1}$ is dominated by
$\pi\sigma^{\ast}$ character. N1-H dissociation occurs rapidly over the next
$15\text{\,}\mathrm{fs}$, eventually reaching a ground state intersection at
an extended N1-H bond length (greater than 2.0 Å).

Current experimental data does not allow us to verify or exclude the existence
of the dissociative channel. We do not find any clear signature of the
$\pi\sigma^{\ast}$ channel in our simulated spectrum (see Figure 2). However,
the channel appears to be small ($\sim$13%), and an oxygen-edge spectrum is
not expected to be highly sensitive to changes at the N1-H bond. A nitrogen-
edge spectrum, on the other hand, should be sensitive to these changes, as
indicated by calculations at selected geometries (see Supporting Information
S7). However, no experimental nitrogen-edge spectrum has been reported. Other
experiments have appeared to disconfirm the existence of a dissociative
$\pi\sigma^{\ast}$ channel. Previous TKER spectra showed only smooth,
isotropic H-atom kinetic energy, indicating no involvement of ultrafast,
$\pi\sigma^{\ast}$-mediated N-H dissociation [36]. However, these spectra only
scanned over excitation wavelengths from 270–230 nm, and dissociation might
become more common at shorter wavelengths. For adenine, the TKER spectra for
wavelengths from 280–234 nm also show smooth, isotropic H-atom kinetic
energies. Only for wavelengths 233 nm and shorter do anisotropic, fast H-atom
peaks consistent with ultrafast $\pi\sigma^{\ast}$-mediated N-H dissociation
appear, and this signature becomes more intense with increasing excitation
energy [35]. This behavior for adenine, coupled with the presence of the
$\pi\sigma^{\ast}$ N-H dissociation in our dynamics, leads us to suggest that
if the TKER experiments for thymine described by Schneider and coworkers [36]
are carried out using excitation wavelengths in the range 230–200 nm, then one
might observe anisotropic, fast H-atom peaks consistent with this dissociative
channel. However, wavelengths in the range 230–200 nm may not be directly
comparable to our simulations as it would also excite the system to states
above the $S_{2}(\pi\pi^{\ast})$ state, in particular, to the $\pi\pi^{\ast}$
band that lies about 1.0 eV above $S_{2}(\pi\pi^{\ast})$ (see Figure 5E and
Supporting Information S9). This is similar to adenine, where several states
may contribute to the $\pi\sigma^{\ast}$ dissociation channel [35]. The
utility of TKER experiments for confirming/denying the presence of the
$\pi\sigma^{\ast}$ N-H dissociation channel in thymine has previously been
suggested by others [38, 34].

Figure 5: The $\pi\sigma^{\ast}$-mediated N-H dissociation pathway. Panel A
illustrates the normal mode associated with the N-H stretch, where $Q$ denotes
the displacement along this mode relative to the ground state equilibrium
geometry ($Q=0$ corresponds to the $S_{0}$ minimum geometry). Panel B shows
potential energy curves along this normal mode, with the associated N-H bond
length $R($N1-H) given on the secondary axis. The shaded area denotes the
values of $Q$ for which we see spawning events in the dynamics. Panel C shows
the natural transition orbitals of $S_{1}$ for $Q=1.0$, where the state is of
$\pi\sigma^{\ast}$ character. Panel D shows the initial geometries (at time
$t=0$) and the spawning geometries (where population is transferred between
$S_{2}$ and $S_{1}$), distinguishing between the dissociative geometries and
the non-dissociative geometries. For each geometry, we give the N-H bond
length and the projected value of the mode displacement $Q$. Dissociating
trajectories have negative initial values of $Q$ (short N-H distances) and
spawn at high positive values of $Q$ (long N-H distances), corresponding to
regions of the potential energy surfaces where the $\pi\sigma^{\ast}$ state
becomes accessible. Panel E shows the recorded ultraviolet gas-phase spectrum
for thymine with a comparison to a simulated CC3 spectrum (see Supporting
Information S9). The shaded area corresponds to wavelengths where TKER spectra
have not been recorded.

The ultimate fate of the dissociative $\pi\sigma^{\ast}$ channel cannot be
resolved in this work, as CC theory is known to break down when the
$\pi\sigma^{\ast}$ state becomes degenerate with the ground state (the
dissociation limit). Because of this, there may also be a pathway back to the
ground state through an intersection with the ground state, which might
compete with N-H bond fission [34]. However, the presence and accessibility of
the $\pi\sigma^{\ast}$ channel is well-described in our simulations. Both the
early-time involvement of the $\pi\sigma^{\ast}$ state, and the initial rapid
N-H stretch, are treated correctly with CC theory, and this aspect of the
simulated dynamics paints a picture in line with both theoretical and
experimental studies of $\pi\sigma^{\ast}$ states in other heteroaromatics
[33, 34].

## Conclusion

The photorelaxation pathways of thymine are still under debate and a clear
consensus has yet to emerge, despite numerous theoretical and experimental
investigations. Here, we have simulated the ultrafast dynamics of thymine
using high-level nonadiabatic dynamics simulations (_ab initio_ multiple
spawning) combined with coupled cluster theory (with single and double
excitations) for the electronic structure. This represents, to the best of our
knowledge, the highest level simulation performed on thymine to date, and the
simulation with the highest level of electronic structure on a molecular
system of this size. From the data, we simulated the oxygen-edge X-ray
absorption spectrum of Wolf et al. [16], obtaining excellent agreement on the
appearance and position of a signal unequivocally associated with the
$n\pi^{\ast}$ state, as well as for the associated $\pi\pi^{\ast}/n\pi^{\ast}$
conversion time: our theoretical estimate of $41\pm 14\text{\,}\mathrm{fs}$
agrees quantitatively with experimental estimates ($60\pm
30\text{\,}\mathrm{fs}$ and $39\pm 1\text{\,}\mathrm{fs}$ [16, 37]).
Furthermore, we find no significant $\pi\pi^{\ast}$ trapping.

Interestingly, our simulation predicts an additional, minor channel in which
population is rapidly transferred to a $\pi\sigma^{\ast}$ state, leading to
N-H dissociation. This type of channel has been implicated in the ultrafast
excited state dynamics of other heteroaromatic systems, including adenine [33,
34]. In view of this surprising finding, we believe further experiments—for
example, measuring the nitrogen-edge X-ray absorption spectrum—are warranted
to confirm or disprove the existence of this dissociative channel as one of
the relaxation pathways in thymine.

This work was made possible by developments in coupled cluster theory. It is
the first nonadiabatic dynamics simulation based on this electronic structure
method (CCSD), and, furthermore, the first simulation based a coupled cluster
method that correctly describes conical intersections. Earlier work by some of
the authors [19, 20, 21] had already indicated that the method could be
modified to correctly describe conical intersections and therefore also
nonadiabatic dynamics. Here, we have shown that the method introduced in these
papers, the similarity constrained coupled cluster method, can in fact be
applied in nonadiabatic dynamics simulations, opening up a range of
applications that may now be studied with the hierarchy of coupled cluster
methods. Given their ability to capture dynamical correlation, we expect that
these methods will shed new light on the photochemistry of a variety of
systems.

## Methods

### Theoretical

The excited state dynamics simulation was performed with the _ab initio_
multiple spawning (AIMS) method [39, 40, 41]. We prepared 16 initial
conditions (ICs) and simulated the dynamics for a total of
$4000\text{\,}\mathrm{au}$ ($\sim 100\text{\,}\mathrm{fs}$). The initial
trajectory basis functions (TBFs) were sampled from a $0\text{\,}\mathrm{K}$
harmonic Wigner distribution obtained from the ground state equilibrium
geometry at the CCSD/aug-cc-pVDZ level. Both the geometry and the frequencies
were obtained with CCSD/aug-cc-pVDZ. The dynamics simulation was performed at
the CCSD/cc-pVDZ level. For the ICs, we adopt the independent first generation
approximation, that is, we average results over 16 independent AIMS
simulations. Of these, 12 run normally for the whole simulation. Two ICs reach
intersections between $S_{0}$ and $S_{1}$ around $20\text{\,}\mathrm{fs}$ and
are terminated at this point because CCSD cannot describe $S_{0}$/$S_{1}$
intersections (termination is caused by convergence problems). A further two
ICs enter a defective region surrounding the $S_{1}$/$S_{2}$ intersection, at
which point the energies become complex-valued and the simulations are
terminated. These two ICs were rerun with similarity constrained CC theory
(SCCSD) using the $\mathscr{E}$ with $T=0$ projection [24]. For both CCSD and
SCCSD, the coupling elements were evaluated with the nuclear derivative acting
on the right vector without normalization; for more details, see Refs. 23, 24.
In the time-integration, we have used a default timestep of
$20\text{\,}\mathrm{au}$ ($0.5\text{\,}\mathrm{fs}$) and a smaller timestep of
$5\text{\,}\mathrm{au}$ ($0.1\text{\,}\mathrm{fs}$) in regions of high
coupling. The spawning threshold was set to 0.05, where the spawning criterion
is given as the norm of the coupling times the velocity. This produced a total
of 67 TBFs over the course of the simulation.

Using the data from the dynamics simulation, we simulated the oxygen-edge
X-ray absorption spectrum using the CC with approximate triples (CC3) [42]
method with the cc-pVDZ basis. Core excited states were obtained with the
core-valence separation approximation [43]. To simulate the spectrum, we have
applied the incoherent approximation [44], where the spectrum is calculated as
an average of the spectra computed at the centers of the TBFs, weighted with
the corresponding TBF amplitude. To avoid state assignment ambiguities, we
have used CCSD for the valence excited states and CC3 for the core excited
states, where the strengths are evaluated with CCSD (excluding the approximate
triples in the CC3 states). Similarly, absorption energies are evaluated from
CCSD-to-CC3 energy differences. All X-ray absorption spectra were shifted by
$-0.5\text{\,}\mathrm{eV}$, corresponding to the required shift needed to
align the first ground state peak with the experimental value at
$531.4\text{\,}\mathrm{eV}$. The UV-vis absorption spectrum was calculated at
the CC3/aug-cc-pVDZ+KBJ(3-4) level for the ground state equilibrium geometry,
determined with CCSD/aug-cc-pVDZ, see Supporting Information S9.

All electronic structure calculations were performed using development
versions of the eT program [45], version 1.8 for the dynamics and version 1.7
for the spectra. The AIMS dynamics was run with the FMS program and an
interface to the eT program. See Supporting Information S8 for a description
of the interface.

### Experimental

The experimental UV spectrum of thymine was taken with a Cary 5E UV-Vis-NIR
spectrometer using an in-house developed gas cell described in Ref. [46]. The
sample was purchased from Sigma-Aldrich with ¿99% purity and used without
further refinement. The cell was heated up to
$150\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to obtain sufficient absorption.
Spectra were recorded over a range of 200–$400\text{\,}\mathrm{nm}$ with a
step size of $0.5\text{\,}\mathrm{nm}$ and an integration time of
$0.5\text{\,}\mathrm{s}$ per data point. A background spectrum of the empty
cell has been recorded at the same temperature and settings and subtracted
from the spectra including the sample.

Acknowledgments This work was supported by the Norwegian Research Council
through FRINATEK project 275506, the European Research Council (ERC) under the
European Union’s Horizon 2020 Research and Innovation Program (Grant No.
101020016), and the AMOS program within the U.S. Department of Energy (DOE),
Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and
Biosciences Division. OJF is a U.S. Department of Energy Computational Science
Graduate Fellow (Grant No. DE-SC0023112). We acknowledge computing resources
through UNINETT Sigma2–the National Infrastructure for High Performance
Computing and Data Storage in Norway, project NN2962k.

## References

  *   * [1] Crespo-Hernandez, C. E., Cohen, B., Hare, P. M. & Kohler, B.  Ultrafast excited-state dynamics in nucleic acid.  _Chem. Rev._ 104, 1977–2020 (2004).
  * [2] Canuel, C. _et al._ Excited states dynamics of DNA and RNA bases: Characterization of a stepwise deactivation pathway in the gas phase.  _J. Chem. Phys._ 122, 074316 (2005).
  * [3] Ullrich, S., Schultz, T., Zgierski, M. Z. & Stolow, A.  Electronic relaxation dynamics in DNA and RNA bases studied by time-resolved photoelectron spectroscopy.  _Phys. Chem. Chem. Phys._ 6, 2796–2801 (2004).
  * [4] McFarland, B. K. _et al._ Ultrafast x-ray auger probing of photoexcited molecular dynamics.  _Nat. Commun._ 5, 4235 (2014).
  * [5] Kang, H., Lee, K. T., Jung, B., Ko, Y. J. & Kim, S. K.  Intrinsic lifetimes of the excited state of DNA and RNA bases.  _J. Am. Chem. Soc._ 124, 12958–12959 (2002).
  * [6] Hudock, H. R. _et al._ Ab initio molecular dynamics and time-resolved photoelectron spectroscopy of electronically excited uracil and thymine.  _J. Phys. Chem. A_ 111, 8500–8508 (2007).
  * [7] Szymczak, J. J. _et al._ Photodynamics simulations of thymine: Relaxation into the first excited singlet state.  _J. Phys. Chem. A_ 113, 12686–12693 (2009).
  * [8] Asturiol, D., Lasorne, B., Robb, M. A. & Blancafort, L.  Photophysics of the $\pi,\pi^{*}$ and $n,\pi^{*}$ states of thymine: MS-CASPT2 minimum-energy paths and CASSCF on-the-fly dynamics.  _J. Phys. Chem. A_ 113, 10211–10218 (2009).
  * [9] Asturiol, D., Lasorne, B., Worth, G. A., Robb, M. A. & Blancafort, L.  Exploring the sloped-to-peaked $S_{2}/S_{1}$ seam of intersection of thymine with electronic structure and direct quantum dynamics calculations.  _Phys. Chem. Chem. Phys._ 12, 4949–4958 (2010).
  * [10] Barbatti, M. _et al._ Relaxation mechanisms of UV-photoexcited DNA and RNA nucleobases.  _Proc. Natl. Acad. Sci. U.S.A._ 107, 21453–21458 (2010).
  * [11] Stojanović, L. _et al._ New insights into the state trapping of UV-excited thymine.  _Molecules_ 21, 1603 (2016).
  * [12] Mai, S., Richter, M., Marquetand, P. & González, L.  The DNA nucleobase thymine in motion – intersystem crossing simulated with surface hopping.  _Chem. Phys._ 482, 9–15 (2017).
  * [13] Chakraborty, P., Liu, Y., Weinacht, T. & Matsika, S.  Effect of dynamic correlation on the ultrafast relaxation of uracil in the gas phase.  _Farad. Disc._ 228, 266–285 (2021).
  * [14] Picconi, D., Barone, V., Lami, A., Santoro, F. & Improta, R.  The interplay between $\pi\pi^{*}/n\pi^{*}$ excited states in gas-phase thymine: A quantum dynamical study.  _ChemPhysChem_ 12, 1957–1968 (2011).
  * [15] Park, W., Lee, S., Huix-Rotllant, M., Filatov, M. & Choi, C. H.  Impact of the dynamic electron correlation on the unusually long excited-state lifetime of thymine.  _J. Phys. Chem. Lett._ 12, 4339–4346 (2021).
  * [16] Wolf, T. J. A. _et al._ Probing ultrafast $\pi\pi^{*}/n\pi^{*}$ internal conversion in organic chromophores via K-edge resonant absorption.  _Nat. Commun._ 8 (2017).
  * [17] Wolf, T. J. A. & Gühr, M.  Photochemical pathways in nucleobases measured with an X-ray FEL.  _Philos. Trans. R. Soc. A_ 377, 20170473 (2019).
  * [18] Mayer, D., Lever, F. & Gühr, M.  Time-resolved x-ray spectroscopy of nucleobases and their thionated analogs.  _J. Photochem. Photobiol._ (2024).
  * [19] Kjønstad, E. F., Myhre, R. H., Martínez, T. J. & Koch, H.  Crossing conditions in coupled cluster theory.  _J. Chem. Phys._ 147 (2017).
  * [20] Kjønstad, E. F. & Koch, H.  Resolving the notorious case of conical intersections for coupled cluster dynamics.  _J. Phys. Chem. Lett._ 8, 4801–4807 (2017).
  * [21] Kjønstad, E. F. & Koch, H.  An orbital invariant similarity constrained coupled cluster model.  _J. Chem. Theory Comput._ 15, 5386–5397 (2019).
  * [22] Schnack-Petersen, A. K., Koch, H., Coriani, S. & Kjønstad, E. F.  Efficient implementation of molecular CCSD gradients with Cholesky-decomposed electron repulsion integrals.  _J. Chem. Phys._ 156, 244111 (2022).
  * [23] Kjønstad, E. F. & Koch, H.  Communication: Non-adiabatic derivative coupling elements for the coupled cluster singles and doubles model.  _J. Chem. Phys._ 158, 161106 (2023).
  * [24] Kjønstad, E. F., Angelico, S. & Koch, H.  Coupled cluster theory for nonadiabatic dynamics: nuclear gradients and nonadiabatic couplings in similarity constrained coupled cluster theory.  _arXiv:2403.01007 [physics.chem-ph]_ (2024).
  * [25] Purvis III, G. D. & Bartlett, R. J.  A full coupled-cluster singles and doubles model: The inclusion of disconnected triples.  _J. Chem. Phys._ 76, 1910–1918 (1982).
  * [26] Hättig, C.  Structure optimizations for excited states with correlated second-order methods: CC2 and ADC(2).  _Adv. Quantum Chem._ 50, 37–60 (2005).
  * [27] Köhn, A. & Tajti, A.  Can coupled-cluster theory treat conical intersections?  _J. Chem. Phys._ 127 (2007).
  * [28] Williams, D. M. G., Kjønstad, E. F. & Martínez, T. J.  Geometric phase in coupled cluster theory.  _J. Chem. Phys._ 158, 214122 (2023).
  * [29] Sobolewski, A. L. & Domcke, W.  On the mechanism of nonradiative deca of DNA bases: ab initio and TDDFT results for the excited states of 9H-adenine.  _Eur. Phys. J. D_ 20, 369–374 (2002).
  * [30] Perun, S., Sobolewski, A. L. & Domcke, W.  Ab initio studies on the radiationless decay mechanisms of the lowest excited singlet states of 9H-adenine.  _J. Am. Chem. Soc._ 127, 6257–6265 (2005).
  * [31] Kang, H. _et al._ Photo-induced dissociation of protonated tryptophan trph+: A direct dissociation channel in the excited states controls the hydrogen atom loss.  _Phys. Chem. Chem. Phys._ 6, 2628–2632 (2004).
  * [32] Gregoire, G., Jouvet, C., Dedonder, C. & Sobolewski, A.  On the role of dissociative $\pi$$\sigma$* states in the photochemistry of protonated tryptamine and tryptophan: An ab initio study.  _Chem. Phys._ 324, 398–404 (2006).
  * [33] Ashfold, M., Cronin, B., Devine, A., Dixon, R. & Nix, M.  The role of $\pi$$\sigma$* excited states in the photodissociation of heteroaromatic molecules.  _Science_ 312, 1637–1640 (2006).
  * [34] Roberts, G. M. & Stavros, V. G.  The role of $\pi$$\sigma$* states in the photochemistry of heteroaromatic biomolecules and their subunits: insights from gas-phase femtosecond spectroscopy.  _Chem. Sci._ 5, 1698–1722 (2014).
  * [35] Nix, M. G. D., Devine, A. L., Cronin, B. & Ashfold, M. N. R.  Ultraviolet photolysis of adenine: Dissociation via the ${}^{1}\pi\sigma^{*}$ state.  _J. Chem. Phys._ 126, 124312 (2007).
  * [36] Schneider, M. _et al._ Photodissociation of thymine.  _Phys. Chem. Chem. Phys._ 8, 3017–3021 (2002).
  * [37] Miura, Y. _et al._ Formation of long-lived dark states during electronic relaxation of pyrimidine nucleobases studied using extreme ultraviolet time-resolved photoelectron spectroscopy.  _J. Am. Chem. Soc._ 145, 3369–3381 (2023).
  * [38] Ashfold, M. N. R. _et al._ $\pi\sigma^{\ast}$ excited states in molecular photochemistry.  _Phys. Chem. Chem. Phys._ 12, 1218–1238 (2010).
  * [39] Ben-Nun, M., Quenneville, J. & Martínez, T. J.  Ab initio multiple spawning: Photochemistry from first principles quantum molecular dynamics.  _J. Phys. Chem. A_ 104, 5161–5175 (2000).
  * [40] Ben-Nun, M. & Martínez, T. J.  Ab initio quantum molecular dynamics.  _Adv. Chem. Phys._ 121, 439–512 (2002).
  * [41] Curchod, B. F. E. & Martínez, T. J.  Ab initio nonadiabatic quantum molecular dynamics.  _Chem. Rev._ 118, 3305–3336 (2018).
  * [42] Koch, H., Christiansen, O., Jørgensen, P., Sanchez de Merás, A. M. & Helgaker, T.  The CC3 model: An iterative coupled cluster approach including connected triples.  _J. Chem. Phys._ 106, 1808–1818 (1997).
  * [43] Coriani, S. & Koch, H.  Communication: X-ray absorption spectra and core-ionization potentials within a core-valence separated coupled cluster framework.  _J. Chem. Phys._ 143 (2015).
  * [44] List, N. H., Dempwolff, A. L., Dreuw, A., Norman, P. & Martínez, T. J.  Probing competing relaxation pathways in malonaldehyde with transient x-ray absorption spectroscopy.  _Chem. Sci._ 11, 4180–4193 (2020).
  * [45] Folkestad, S. D. _et al._ eT 1.0: An open source electronic structure program with emphasis on coupled cluster and multilevel methods.  _J. Chem. Phys._ 152 (2020).
  * [46] Mayer, D., Picconi, D., Robinson, M. S. & Gühr, M.  Experimental and theoretical gas-phase absorption spectra of thionated uracils.  _Chem. Phys._ 558, 111500 (2022).

See pages - of SI.pdf


# Drell Yan production in the parton branching method and Powheg matching

A. Bermudez Martinez  L.I. Estevez Banos  Deutsches Elektronen-Synchrotron
DESY, Germany  M. Fernández Moreira  Higher Institute of Technology and
Applied Sciences, University of Havana, InSTEC  J. Girones Dominguez  Higher
Institute of Technology and Applied Sciences, University of Havana, InSTEC  H.
Jung  Deutsches Elektronen-Synchrotron DESY, Germany  R. Ramos Blazquez
Higher Institute of Technology and Applied Sciences, University of Havana,
InSTEC

###### Abstract

As part of the DESY summer student program 2021, transverse momentum dependent
(TMD) parton distributions obtained from the Parton Branching (PB) method were
combined with next-to-leading-order (NLO) using the Powheg method.
Computations of the resulting Drell-Yan (DY) transverse momentum spectrum were
performed. A good agreement of the theoretical predictions with the
measurement performed by the CMS experiment at the center-of-mass energy of 13
TeV is found, at low and intermediate DY $p_{\rm T}$. The new scale choice
option for the matching has been included in the Cascade3 event generator.

DESY-21-221

##  1 Introduction

The precise description of the Drell-Yan [1] (DY) $p_{\rm T}$ spectrum in
$\mathup{{{p}}}\mathup{{{p}}}$ collisions requires higher-order perturbative
calculations at large transverse momentum $p_{\rm T}$and the resummation of
soft gluons at all orders at small $p_{\rm T}$ [2, 3, 4, 5, 6] as it is done
in analytical resummation [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
and collinear parton shower approaches [20, 21, 22, 23] matched to higher-
order calculations [24, 25, 26, 27, 28, 29]. In [30, 31] the proposed
combination of the Parton Branching formulation [32, 33] of TMD evolution with
next-to-leading order (NLO) calculations using the mc@nlo method provided a
good description of the DY $p_{\rm T}$ spectrum as measured by the CMS [34],
ATLAS [35], NuSea [36, 37], R209 [38], and PHENIX [39] experiments,
corresponding to a wide range of DY masses and center-of-mass energies. The
sensitivity of the predictions to non-perturbative contribution from the TMD
was found small. In addition, in [40] a TMD multi-jet merging approach at
leading-order was developed which provides a very good description over the
whole DY $p_{\rm T}$ spectrum.

An alternative NLO matching method is Powheg [26]. The Powheg formulation
produces positive weighted events, and as opposed to mc@nlo, the parton level
events it provides do not depend on the specific parton shower used to
simulate the parton radiation. In this paper we apply the PB approach to
calculate the DY $p_{\rm T}$ spectrum at NLO using the Powheg method. We use
TMDs obtained from a fit to inclusive deep-inelastic scattering (DIS)
precision data [41] from HERA. This work was done as part of the DESY summer
student program 2021 were young scientists connected regularly from remote,
working in a virtual environment to achieve the results hereby presented.

###  1.1 Matching PB-TMDs with Powheg NLO calculations

In the following we will describe how to use the PB-TMDs together with higher
order perturbative calculations using the Powheg method. We make use of the
Powheg-BOX [26, 42, 43] (version 2). Powheg generates one emission similar to
the first step of a parton shower, with the difference that the evolution
kernel is given by the ratio of the real emission and Born matrix element,
rather than the splitting functions. Events generated by Powheg are weighted
with a Sudakov factor determined by the exponentiation of the Powheg evolution
kernel.

The NLO generator Powheg-BOX produces events in the LHE format [44] which can
then be read by parton shower event generators. We use the events output by
Powheg-BOX and apply the PB-TMDs to modify the kinematics of the initial-state
partons (and correspondingly the final state) according to the transverse
momentum distributions given by the PB-TMDs. The addition of the transverse
momentum to the system leads to shifts in the longitudinal momentum of the
incoming partons while the invariant mass and rapidity of the partonic system
are required to remain unchanged.

The DY transverse momentum spectrum calculated at pure NLO within the Powheg-
BOX is shown in Fig. 1 (left). In order to obtained the pure NLO prediction we
suppressed the action of the Powheg Sudakov by chosing a very small value for
the hdamp parameter (${\texttt{hdamp}}=0.001$). In Fig. 1 (right) the pure
Powheg calculation (LHE level) without any suppresion of the Powheg Sudakov is
shown. The predictions are overlaid with the DY $p_{\rm T}$ spectrum
measurement performed by the CMS collaboration at 13 TeV of center-of-mass
energy [34].

Figure 1: (Left) DY $p_{\rm T}$ spectrum as obtained from a pure NLO
calculation. (Right) DY $p_{\rm T}$ spectrum as obtained from a pure Powheg
(LHE level). The calculations are shown together with the measurement of the
DY $p_{\rm T}$ spectrum by the CMS experiment at the center-of-mass energy of
13 TeV [34].

As observed in Fig. 1, the Powheg method provides a positive spectrum at the
LHE level compared to pure NLO, as well as compared to mc@nlo [30]. The
potential of a combined PB-TMD + Powheg calculation can be seen in Fig. 1
(right), where a poor descrition by pure Powheg of low and intermediate DY
$p_{\rm T}$ is observed.

In order to match the PB-TMDs to the Powheg calculation we set the
$p_{\perp,\text{min}}$ parameter in Powheg by changing the input variable
ptsqmin. Below $p_{\perp,\text{min}}$ Powheg produces n–parton kinematics
while above $p_{\perp,\text{min}}$ it produces (n+1)–parton events. In Fig. 2
we show the result at the LHE level of the Powheg calculation using
$\texttt{ptsqmin}=\text{361 GeV}^{2}$.

Figure 2: DY $p_{\rm T}$ spectrum as obtained from a pure Powheg (LHE level)
setting $\texttt{ptsqmin}=\text{361 GeV}^{2}$. The calculation is shown
together with the measurement of the DY $p_{\rm T}$ spectrum by the CMS
experiment at the center-of-mass energy of 13 TeV [34].

In order to avoid double counting between the real emission produced by Powheg
and emissions from the PB-TMD parton shower the matching scale is set equal to
the parameter SCALUP (included in the LHE Powheg output). The Powheg-BOX
generator by default attributes a value
$\texttt{SCALUP}=\sqrt{\texttt{ptsqmin}}$ which is then used as the upper
scale for the transverse momentum generated by the TMD.

The inclusion of transverse momenta from the PB-TMD and the parton shower is
performed using the Cascade3 event generator [45]. The factorization scale
$\mu$ choice was chosen according to the performance of the matched PB-TMD +
Powheg calculation on the differential jet rate $d_{01}$, which represents the
square of the energy scale at which a 0–jet event is resolved as an 1–jet
event. A smoother $d_{01}$ implies a better matching between the PB-TMD shower
contribution to the 0–parton calculation and the real emission contribution.
The choice $\mu=\sum_{i}\sqrt{m_{i}^{2}+p_{t,i}^{2}}$ for 1–parton events and
$\mu=\sqrt{\hat{s}}$ for 0–parton events (where $\hat{s}$ is the partonic
center-of-mass energy) showed the smoothest $d_{01}$ distribution. This option
has been added to the Cascade3 event generator.

In Fig. 3 the described matching procedure is applied to DY production. In
Fig. 3 (left) the $d_{01}$ distribution is shown together with the
contributions from the 0– and 1–jet calculations matched to the PB-TMD parton
shower. In Fig. 3 (right) the computation of the DY $p_{\rm T}$ spectrum is
depicted.

Figure 3: (Left) The $d_{01}$ spectrum representing the energy-square scale at
which an 0–jet event is resolved as an 1–jet event. The dotted curves
represent the contributions of the single-multiplicity samples. (Right) DY
$p_{\rm T}$ spectrum as obtained from a matched PB-TMD + Powheg computation
compared to the measurement of the DY $p_{\rm T}$ spectrum by the CMS
experiment at the center-of-mass energy of 13 TeV [34].

As observed in Fig. 3, a good matching of the PB-TMD with Powheg NLO
calculation of DY production is achieved, resulting in a good description of
the DY $p_{\rm T}$ spectrum at low and intermediate $p_{\rm T}$ . The results
are in agreement with the previous results obtained using mc@nlo matching [30,
31]. The difference observed between the data and the calculation at high
$p_{\rm T}$ is due to missing higher orders.

We should point out that further refinements of the method are under
investigation. Among these is to avoid the posibility of having parton shower
emissions that are are harder than the Powheg real emission even though the
resulting transverse momentum from the TMD is not. This effect is due to the
angular ordering of the partonic emissions in the PB-TMD parton shower.

##  2 Conclusion

The DESY summer student program ran exclusively online in its 2021 edition. In
a dedicated effort, the students worked on the matching of PB-TMDs determined
from NLO fits to DIS data, and the resulting PB-TMD parton shower, to NLO
calculation using the Powheg method. The computation has been compared to the
DY $p_{\rm T}$ spectrum measured by the CMS collaboration at the center-of-
mass energy of 13 TeV. A very good agreement at small and intermediate $p_{\rm
T}$ was observed, which is also compatible with results obtained using the
mc@nlo method [30, 31]. We should poin out that this performance is achived
without further adjustment of parton shower parameters as opposed to collinear
parton shower approaches. The new scale choice option for the matching has
been included in the Cascade3 event generator.

Although a smooth $d_{01}$ distribution distribution has been observed,
further refinement of the calculation might be needed in order to avoid
remaining double counting of emissions, which has its origin in the angular
ordering of the PB-TMD parton emissions.

## References

  * [1] S. Drell and T.-M. Yan, Phys.Rev.Lett. 25, 316 (1970)
  * [2] Y. L. Dokshitzer, D. Diakonov and S. I. Troian, Phys. Lett. B 79 (1978), 269-272
  * [3] G. Parisi and R. Petronzio, Nucl. Phys. B 154 (1979), 427-440
  * [4] G. Curci and M. Greco, Phys. Lett. B 92 (1980), 175-178
  * [5] G. Altarelli, R. K. Ellis and G. Martinelli, Z. Phys. C 27 (1985), 617
  * [6] J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys. B 250 (1985), 199-224
  * [7] W. Bizon, X. Chen, A. Gehrmann-De Ridder, T. Gehrmann, N. Glover, A. Huss, P. F. Monni, E. Re, L. Rottoli and P. Torrielli, JHEP 12 (2018), 132 [arXiv:1805.05916]
  * [8] W. Bizon, A. Gehrmann-De Ridder, T. Gehrmann, N. Glover, A. Huss, P. F. Monni, E. Re, L. Rottoli and D. M. Walker, Eur. Phys. J. C 79 (2019) no.10, 868 [arXiv:1905.05171]
  * [9] S. Catani, D. de Florian, G. Ferrera and M. Grazzini, JHEP 12 (2015), 047 [arXiv:1507.06937]
  * [10] I. Scimemi and A. Vladimirov, Eur. Phys. J. C 78 (2018) no.2, 89 [arXiv:1706.01473]
  * [11] A. Bacchetta, G. Bozzi, M. Lambertsen, F. Piacenza, J. Steiglechner and W. Vogelsang, Phys. Rev. D 100 (2019) no.1, 014018 [arXiv:1901.06916]
  * [12] A. Bacchetta, G. Bozzi, M. Radici, M. Ritzmann and A. Signori, Phys. Lett. B 788 (2019), 542-545 [arXiv:1807.02101]
  * [13] G. A. Ladinsky and C. P. Yuan, Phys. Rev. D 50 (1994), R4239 [arXiv:hep-ph/9311341]
  * [14] C. Balazs and C. P. Yuan, Phys. Rev. D 56 (1997), 5558-5583 [arXiv:hep-ph/9704258]
  * [15] F. Landry, R. Brock, P. M. Nadolsky and C. P. Yuan, Phys. Rev. D 67 (2003), 073016 [arXiv:hep-ph/0212159]
  * [16] P. Nadolsky et al., http://hep.pa.msu.edu/resum/
  * [17] S. Alioli, C. W. Bauer, C. Berggren, F. J. Tackmann and J. R. Walsh, Phys. Rev. D 92 (2015) no.9, 094020 [arXiv:1508.01475]
  * [18] G. Bozzi and A. Signori, Adv. High Energy Phys. 2019 (2019), 2526897 [arXiv:1901.01162]
  * [19] S. P. Baranov, A. V. Lipatov and N. P. Zotov, Phys. Rev. D 89 (2014) no.9, 094025 [arXiv:1402.5496]
  * [20] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands, Comput. Phys. Commun. 191 (2015), 159-177 [arXiv:1410.3012]
  * [21] J. Bellm, S. Gieseke, D. Grellscheid, S. Plätzer, M. Rauch, C. Reuschle, P. Richardson, P. Schichtel, M. H. Seymour and A. Siódmok, et al. Eur. Phys. J. C 76 (2016) no.4, 196 [arXiv:1512.01178]
  * [22] M. Bahr, S. Gieseke, M. A. Gigg, D. Grellscheid, K. Hamilton, O. Latunde-Dada, S. Platzer, P. Richardson, M. H. Seymour and A. Sherstnev, et al. Eur. Phys. J. C 58 (2008), 639-707 [arXiv:0803.0883]
  * [23] T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert and J. Winter, JHEP 02 (2009), 007 [arXiv:0811.4622]
  * [24] S. Frixione, P. Nason and B. R. Webber, JHEP 08 (2003), 007 [arXiv:hep-ph/0305252]
  * [25] S. Frixione and B. R. Webber, JHEP 06 (2002), 029 [arXiv:hep-ph/0204244]
  * [26] S. Frixione, P. Nason and C. Oleari, JHEP 11 (2007), 070 [arXiv:0709.2092]
  * [27] P. Nason and B. Webber, Ann. Rev. Nucl. Part. Sci. 62 (2012), 187-213 [arXiv:1202.1251]
  * [28] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07 (2014), 079 [arXiv:1405.0301]
  * [29] R. Frederix, S. Frixione, A. Papaefstathiou, S. Prestel and P. Torrielli, JHEP 02 (2016), 131 [arXiv:1511.00847]
  * [30] A. Bermudez Martinez, P. Connor, D. Dominguez Damiani, L. I. Estevez Banos, F. Hautmann, H. Jung, J. Lidrych, M. Schmitz, S. Taheri Monfared and Q. Wang, et al. Phys. Rev. D 100 (2019) no.7, 074027 [arXiv:1906.00919]
  * [31] A. Bermudez Martinez, P. L. S. Connor, D. Dominguez Damiani, L. I. Estevez Banos, F. Hautmann, H. Jung, J. Lidrych, A. Lelek, M. Mendizabal and M. Schmitz, et al. Eur. Phys. J. C 80 (2020) no.7, 598 [arXiv:2001.06488]
  * [32] F. Hautmann, H. Jung, A. Lelek, V. Radescu and R. Zlebcik, Phys. Lett. B 772 (2017), 446-451 [arXiv:1704.01757]
  * [33] F. Hautmann, H. Jung, A. Lelek, V. Radescu and R. Zlebcik, JHEP 01 (2018), 070 [arXiv:1708.03279]
  * [34] A. M. Sirunyan et al. [CMS], JHEP 12 (2019), 061 [arXiv:1909.04133]
  * [35] G. Aad et al. [ATLAS], Eur. Phys. J. C 76 (2016) no.5, 291 [arXiv:1512.02192]
  * [36] J. C. Webb et al. [NuSea], [arXiv:hep-ex/0302019]
  * [37] J. C. Webb, [arXiv:hep-ex/0301031]
  * [38] G. Moreno, C. N. Brown, W. E. Cooper, D. Finley, Y. B. Hsiung, A. M. Jonckheere, H. Jostlein, D. M. Kaplan, L. M. Lederman and Y. Hemmi, et al. Phys. Rev. D 43 (1991), 2815-2836
  * [39] C. Aidala et al. [PHENIX], Phys. Rev. D 99 (2019) no.7, 072003 [arXiv:1805.02448]
  * [40] A. B. Martinez, F. Hautmann and M. L. Mangano, Phys. Lett. B 822 (2021), 136700 [arXiv:2107.01224]
  * [41] A. Bermudez Martinez, P. Connor, H. Jung, A. Lelek, R. Žlebčík, F. Hautmann and V. Radescu, Phys. Rev. D 99 (2019) no.7, 074008 [arXiv:1804.11152]
  * [42] S. Alioli, P. Nason, C. Oleari and E. Re, JHEP 06 (2010), 043 [arXiv:1002.2581]
  * [43] P. Nason, JHEP 11 (2004), 040 [arXiv:hep-ph/0409146]
  * [44] J. Alwall, A. Ballestrero, P. Bartalini, S. Belov, E. Boos, A. Buckley, J. M. Butterworth, L. Dudko, S. Frixione and L. Garren, et al. Comput. Phys. Commun. 176 (2007), 300-304 [arXiv:hep-ph/0609017]
  * [45] S. Baranov, A. Bermudez Martinez, L. I. Estevez Banos, F. Guzman, F. Hautmann, H. Jung, A. Lelek, J. Lidrych, A. Lipatov and M. Malyshev, et al. Eur. Phys. J. C 81 (2021) no.5, 425 [arXiv:2101.10221]