# Lighting up the Photon Wigner Distribution via Dilepton Productions

Yu Shi<EMAIL_ADDRESS>Key Laboratory of Particle Physics and
Particle Irradiation (MOE), Institute of frontier and interdisciplinary
science, Shandong University, Qingdao, Shandong 266237, China  CPHT, CNRS,
École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
Lin Chen<EMAIL_ADDRESS>School of Science and Engineering, The Chinese
University of Hong Kong, Shenzhen, Guangdong, 518172, P.R. China  University
of Science and Technology of China, Hefei, Anhui, 230026, P.R.China  Instituto
Galego de Física de Altas Enerxías IGFAE, Universidade de Santiago de
Compostela, E-15782 Galicia-Spain  Shu-Yi Wei<EMAIL_ADDRESS>Key Laboratory
of Particle Physics and Particle Irradiation (MOE), Institute of frontier and
interdisciplinary science, Shandong University, Qingdao, Shandong 266237,
China  Bo-Wen Xiao<EMAIL_ADDRESS>School of Science and Engineering,
The Chinese University of Hong Kong, Shenzhen, Guangdong, 518172, P.R. China
Southern Center for Nuclear-Science Theory (SCNT), Institute of Modern
Physics, Chinese Academy of Sciences, Huizhou 516000, Guangdong Province,
China

###### Abstract

We present a systematic investigation of lepton pair production through
photon-photon fusion processes in heavy-ion collisions. It is demonstrated
that the dilepton production at a given impact parameter ($b_{\perp}$) with a
fixed transverse momentum imbalance ($q_{\perp}$) can be factorized into a
unified formula in terms of the Wigner photon distribution of heavy nuclei. We
show that this framework provides a comprehensive description of all the
relevant data from RHIC to the LHC, with a strong evidence that the quasi-real
photon can be radiated not only from the nucleus as a whole, standing for the
coherent contribution, but also from the sub-structures inside the nucleus,
representing the incoherent contribution. Further predictions are made for the
anisotropies in the correlations between $q_{\perp}$, $b_{\perp}$, and the
dilepton transverse momentum ($P_{\perp}$). This will help us to constrain the
photon Wigner distribution which plays a crucial role to study the gluonic
matter of nucleus at small-$x$ through the diffractive photoproduction
processes in heavy ion collision.

Introduction Dilepton production in photon-induced events measured in ultra-
peripheral collisions (UPCs) has become a rich physics program at the
Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC).
Due to a large photon flux surrounding the relativistic heavy nuclei, one can
measure the dilepton production at an abundant rate and provide detailed
information for this process at various transverse momentum ($P_{\perp}$), the
transverse momentum imbalance ($q_{\perp}$), and the impact parameter
($b_{\perp}$). In recent years, abundant high-precision experimental data on
dilepton production through photon-photon fusion process
($\gamma\gamma\rightarrow l^{+}l^{-}$) has been collected at RHIC STAR:2004bzo
; STAR:2018ldd ; STAR:2019wlg  and the LHC ATLAS:2016vdy ; ATLAS:2018pfw ;
ATLAS:2020epq ; CMS:2020skx ; ALICE:2022hvk ; ATLAS:2022yad ; ATLAS:2022srr ;
ATLAS:2022ryk ; CMS:2022arf . A large number of diverse observables from the
experimental measurement have been provided to understand and study its
production mechanism. As it provides a clean environment to study the photon
distributions from heavy nuclei, this process has attracted much attention
which leads to a series of theoretical researches Brodsky:1971ud ;
Terazawa:1973tb ; Budnev:1975poe ; Bertulani:1987tz ; Cahn:1990jk ;
Bottcher:1991nd ; Vidovic:1992ik ; Hencken:1994my ; Krauss:1997vr ;
Guclu:1999fx ; Baur:2001jj ; Baltz:2002pp ; Baur:2003ar ; Hencken:2004td ;
Bertulani:2005ru ; Baur:2007zz ; Baltz:2007kq ; Baltz:2009jk ; Klein:2018cjh ;
Klein:2018fmp ; Klusek-Gawenda:2018zfz ; Vysotskii:2018eic ; Zha:2018tlq ;
Zha:2018tlq ; Azevedo:2019fyz ; Li:2019yzy ; Li:2019sin ; Klein:2020jom ;
Klein:2020fmr ; Xiao:2020ddm ; Karadag:2019gvc ; Klusek-Gawenda:2020eja ;
Brandenburg:2020ozx ; Zha:2020cst ; Hatta:2021jcd ; Zha:2021jhf ; Wang:2021kxm
; Wang:2022ihj ; Wang:2022gkd ; Brandenburg:2022tna ; Lin:2022flv ;
Shao:2022stc ; Luo:2023syp ; Shao:2023zge ; Dai:2024imb ; Zhang:2024mql  and
Monte-Carlo studies Klein:2016yzr ; Harland-Lang:2018iur ; Burmasov:2021phy ;
Shao:2022cly ; Harland-Lang:2023ohq . Additionally, this process improves our
understanding of photo-nuclear interaction and the structure of the nucleus,
and is important to explore new physics beyond the standard model
DELPHI:2003nah ; ATLAS:2022ryk ; Knapen:2016moh ; CMS:2018erd ; Ellis:2017edi
; Xu:2022qme ; Shao:2023bga  and study QED properties under extreme conditions
Baltz:1997di ; Baur:1998ay ; Klein:2020fmr ; Steinberg:2021lfm ;
Hattori:2020htm ; Copinger:2020nyx ; Brandenburg:2021lnj .

Currently, different high-precision physical observables from RHIC and the LHC
present a challenge to theorists: to develop a unified factorized framework
incorporating both impact parameter dependence and transverse momentum
dependence, which can account for all experimental observables, such as
dilepton anisotropies and $q_{\perp}$ broadening. Based on previous
studiesLi:2019yzy ; Li:2019sin ; Klein:2020jom ; Klein:2020fmr ; Xiao:2020ddm
, we derive a generalized factorization framework on dilepton production in
terms of the photon fusion hard part and the photon Wigner distributions,
which provide five-dimensional photon distribution within the nuclei. We will
show that this framework not only reproduces the results from previous
studiesLi:2019yzy ; Li:2019sin ; Klein:2020jom ; Xiao:2020ddm , but also
reveal a range of new anisotropic correlations among the vectors
$\vec{P}_{\perp}$, $\vec{q}_{\perp}$, and $\vec{b}_{\perp}$. These
anisotropies are useful for probing the properties of linearly polarized
photons in heavy nuclei.

With a detailed analysis, we also find that the quasi-real photons not only
can be emitted coherently from the nucleus as a whole but also incoherently
from the substructures within the nucleus. Specifically, the incoherent
contribution is one of the production mechanisms responsible for $q_{\perp}$
broadening, and it predominantly contributes to events with neutron emissions
in the large-$q_{\perp}$ region. Thus, we provide a thorough description of
existing experimental measurements from RHIC and LHC within this unified
framework.

More recently, diffractive photo-production processes in heavy ion collisions,
such as heavy quarkonium production, have been shown to provide a unique probe
of gluon tomography in nuclei at small-$x$ Klein:2019qfb ; Mantysaari:2020axf
; Mantysaari:2017slo ; Morreale:2021pnn ; Mantysaari:2023prg ;
Mantysaari:2023xcu ; Kovchegov:2023bvy ; Cepila:2023dxn ; ALICE:2021gpt ;
CMS:2023snh ; LHCb:2022ahs ; ALICE:2023jgu ; STAR:2023nos . Because these
processes also depend on the incoming photon Wigner distribution, the
formalism derived in our paper will help solidify the theoretical foundation
for these studies. In other words, a comparative study between dilepton
production through photon-photon fusion and diffractive photo-production will
assist in building a rigorous framework to explore gluonic matter in cold
nuclei during heavy ion collisions.

The generalized factorization Let us describe the generalized factorization
for the dilepton production in the two-photon fusion process as follows

$\gamma(k_{1})+\gamma(k_{2})\rightarrow l^{+}(p_{1})+l^{-}(p_{2}).$ (1)

In heavy-ion collisions, two high-energy nuclei, separated by a transverse
distance $b_{\perp}$, emit a significant number of quasi-real photons with low
virtuality. Two photons from different sources collide to produce a lepton
pair with transverse momentum $p_{1\perp}$ and $p_{2\perp}$. The collision
point is distanced at $b_{1\perp}$ and $b_{2\perp}$ from the two heavy nuclei,
respectively. In both RHIC and the LHC experiments, measurements have shown
that the transverse momentum of the final lepton is often much larger than the
total momentum of the lepton pair,
$\vec{q}_{\perp}=\vec{p}_{1\perp}+\vec{p}_{2\perp}$. Consequently, we can
assume that the relative momentum of the lepton pair,
$\vec{P}_{\perp}=(\vec{p}_{1\perp}-\vec{p}_{2\perp})/2$ with $P_{\perp}\gg
q_{\perp}$, which is widely referred to as the back-to-back correlation limit.
In this situation, one often needs to use the generalized equivalent photon
approximation and the photon Wigner distribution to describe these quasi-real
photons. The photon Wigner distribution, a quasi probability distribution
function, includes the five-dimensional information of longitudinal momentum
fraction, transverse momentum, and impart parameter of the photon, which is
defined as vonWeizsacker:1934nji ; Williams:1934ad ; Ji:1996ek ;
Belitsky:2000jk ; Diehl:2001pm ; Ji:2003ak ; Belitsky:2003nz ; Hatta:2016dxp ;
Klein:2020jom

$\displaystyle\Gamma^{ij}(x,\vec{k}_{T},\vec{b}_{\perp})=\int\frac{d^{2}\Delta_{\perp}}{(2\pi)^{2}}e^{i\Delta_{\perp}\cdot
b_{\perp}}\Gamma^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp}),$ (2)

and $\Gamma^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$ represents the so-called
Generalized Transverse Momentum Dependent (GTMD) photon distributions in the
momentum space. It reads

$\displaystyle\Gamma^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=\int\frac{d\xi^{-}d^{2}r_{\perp}}{(2\pi)^{3}}e^{ixP^{+}\xi^{-}-ik_{T}\cdot
r_{\perp}}$ $\displaystyle\times\langle
A,-\frac{\Delta_{\perp}}{2}|F^{+i}(0,\frac{r_{\perp}}{2})F^{+j}(\xi^{-},-\frac{r_{\perp}}{2})|\frac{\Delta_{\perp}}{2},A\rangle,$
(3)

with the two dimensional polarization indices $i$ for the incoming photon in
the scattering amplitude and the polarization indices $j$ for the photon in
the conjugate amplitude. $k_{T}$ represents the transverse momentum of the
photon, and $\Delta_{\perp}$ corresponds to the Fourier transform of
$b_{\perp}$. Meanwhile, $x$ stands for the longitudinal momentum fractions of
the photon per nucleon. Conventionally, the GTMD photon distributions can be
parametrized into two parts: the unpolarized photon distribution
$xf_{\gamma}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$ and the linearly-polarized
photon distribution $xh_{\gamma}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$. It is
defined as

$\displaystyle\Gamma^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$
$\displaystyle=$
$\displaystyle\frac{\delta^{ij}}{2}xf_{\gamma}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$
(4) $\displaystyle+$
$\displaystyle\left(\frac{k^{i}_{+}k^{j}_{-}}{\vec{k}_{-}\cdot\vec{k}_{+}}-\frac{\delta^{ij}}{2}\right)xh_{\gamma}(x,\vec{k}_{T},\vec{\Delta}_{\perp}),$

with $k_{\pm}=k_{T}\pm\Delta_{\perp}/2$. In the correlation limit assumption,
the dilepton production from the two-photon fusion process can be separated
into two parts: the Wigner distribution and the hard scattering part. Thus,
the factorization formula reads as follows

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}q_{\perp}d^{2}P_{\perp}d^{2}b_{\perp}}=\int\Gamma_{ij}(x_{1},\vec{k}_{1T},\vec{b}_{1\perp})\otimes$
$\displaystyle\quad\quad\quad\Gamma_{kl}(x_{2},\vec{k}_{2T},\vec{b}_{2\perp})\otimes\mathcal{H}^{ijkl}(\vec{P}_{\perp}).$
(5)

The two photon Wigner distributions, $\Gamma_{ij}$ and $\Gamma_{kl}$, contain
the photon distribution information with the momentum scales significantly
less than $P_{\perp}$. The hard factor, $\mathcal{H}^{ijkl}$, represents the
hard scale part and depends on the momentum imbalance of the lepton pair
$P_{\perp}$ and its angle $\phi_{P}$. Since the lepton mass $m_{l}$ is much
smaller than $P_{\perp}$, we have neglected the lepton mass in the hard
factor. The hard factor can be expressed as follows

$\displaystyle\mathcal{H}^{ijkl}(\vec{P}_{\perp})$ $\displaystyle=$
$\displaystyle\frac{\alpha_{{\rm
em}}^{2}}{\hat{s}^{2}}\left[O^{ijkl}(\vec{P}_{\perp})-4\Omega^{ijkl}(\vec{P}_{\perp})\right],$
(6)

with the projection hard parts $O^{ijkl}$ and $\Omega^{ijkl}$ defined as

$\displaystyle O^{ijkl}(\vec{P}_{\perp})$ $\displaystyle=$ $\displaystyle
2\left(\frac{\hat{u}}{\hat{t}}+\frac{\hat{t}}{\hat{u}}\right)\left(\delta^{ij}\delta^{kl}-\delta^{ik}\delta^{jl}+\delta^{il}\delta^{jk}\right),$
(7) $\displaystyle\Omega^{ijkl}(\vec{P}_{\perp})$ $\displaystyle=$
$\displaystyle 2\Pi^{ij}(\vec{P}_{\perp})\Pi^{kl}(\vec{P}_{\perp})$ (8)
$\displaystyle-\left(\delta^{il}\delta^{jk}+\delta^{ik}\delta^{jl}-\delta^{ij}\delta^{kl}\right),$

where $\hat{s}$, $\hat{u}$ and $\hat{t}$ are the Mandelstam variables and
$\Pi^{ij}(\vec{P}_{\perp})=2\frac{P_{\perp}^{i}P_{\perp}^{j}}{P_{\perp}^{2}}-\delta^{ij}$.
The first term $O^{ijkl}$ is isotropic in terms of the angle $\phi_{P}$ as
derived in Refs. Vidovic:1992ik ; Klein:2020jom . The second term,
$\Omega^{ijkl}$, can generate azimuthal asymmetries between two linearly
polarized photons. Given an arbitrary two-dimensional unit vector $\hat{n}$,
it is interesting to note the following tensor contractions:
$\hat{n}_{i}\hat{n}_{j}\Pi^{ij}=\cos(2\Delta\phi)$ and
$\hat{n}_{i}\hat{n}_{j}\hat{n}_{k}\hat{n}_{l}\Omega^{ijkl}=\cos(4\Delta\phi)$,
where $\Delta\phi$ represents the azimuthal angle between $\vec{P}_{\perp}$
and $\hat{n}$. Integrating over the angle of impact parameter $\phi_{b}$, we
find that the $\Omega^{ijkl}$ term projects the cross-section onto the
$\cos(4\phi_{P}-4\phi_{q})$ harmonics, reproducing the result derived in Ref.
Li:2019sin . Similarly, after integrating over $q_{\perp}$, one obtains the
$\cos(4\phi_{P}-4\phi_{b})$ azimuthal asymmetry proposed in Ref. Xiao:2020ddm
from the $\Omega^{ijkl}$ tensor. In general, various azimuthal angle
asymmetries can arise from inner products of the three vectors
$\vec{P}_{\perp}$, $\vec{q}_{\perp}$, and $\vec{b}_{\perp}$.

As discussed in Ref. Klein:2020jom , it is vital to demonstrate that this
dilepton production cross-section always remains positive definite although
the photon Wigner distribution may become negative. Let us rewrite the photon
GTMD distribution as

$\Gamma^{ij}(x_{1},\vec{k}_{1T},\vec{\Delta}_{\perp})=\frac{Z^{2}\alpha_{\rm
em}}{\pi^{2}}k^{i}_{1\perp}{k^{\prime}}^{j}_{1\perp}\frac{F_{A}(k_{1})}{k_{1}^{2}}\frac{F_{A}(k_{1}^{\prime})}{{k^{\prime}}_{1}^{2}},\,$
(9)

where $F_{A}(k_{1})$ is the elastic charge form factor, $Z$ is charge number
and $m_{p}$ is the mass of proton with
$k_{1}^{2}=x_{1}^{2}m_{p}^{2}+k_{1\perp}^{2}$ and
${k^{\prime}}_{1}^{2}=x_{1}^{2}m_{p}^{2}+{k^{\prime}}_{1\perp}^{2}$. As
detailed in the supplemental material, one finds

$d\sigma\geq\left(\frac{Z^{2}}{\pi^{2}}\right)^{2}\frac{8\alpha_{{\rm
em}}^{4}}{(2\pi)^{2}\hat{s}^{2}}(G^{12}+G^{21})(G^{12*}+G^{21*}),$ (10)

where $G^{ik}\equiv\int d^{2}k_{1\perp}e^{ik_{1\perp}\cdot
b_{\perp}}F_{A}(k_{1})F_{A}(k_{2})|k_{1\perp}||k_{2\perp}|\hat{l}_{1}^{i}\hat{l}_{2}^{k}$.
The $\hat{l}_{i}$ here refers to a unit vector, pointing the direction with
the angle $\phi_{l_{i}}=\phi_{k_{i\perp}}-\phi_{P}$. Thus, it is clear that
the cross-section is ensured to be positive definite.

Phenomenology Let us employ the unified factorization formalism to calculate
various relevant observables, with a comprehensive phenomenological study
available in the supplemental material.

First, the cross-section contains abundant angular correlation between the
impact parameter and the transverse momentum of the final particles, and it
can be decomposed as follows

$\displaystyle\sigma$ $\displaystyle=$
$\displaystyle\sigma_{0}+\sigma_{2qb}\cos(2\phi_{q}-2\phi_{b})+\sigma_{4qb}\cos(4\phi_{q}-4\phi_{b})$
(11)
$\displaystyle+\sigma_{4qP}\cos(4\phi_{q}-4\phi_{P})+\sigma_{4bP}\cos(4\phi_{b}-4\phi_{P})$
$\displaystyle+\sigma_{qbP}\cos(2\phi_{q}+2\phi_{b}-4\phi_{P})+...,$

where the coefficients $\sigma_{2qb},\sigma_{4qb},...$ encode various
information of the photon flux. As shown in Refs. Li:2019yzy ; Li:2019sin ;
Xiao:2020ddm ; Klein:2020jom ; Klein:2020fmr ; Klusek-Gawenda:2020eja ;
Brandenburg:2020ozx ; Zha:2020cst ; Shao:2022stc ; Shao:2023zge , the linearly
polarized photons can generate a range of angular correlations either between
$P_{\perp}$ and $q_{\perp}$ (including $\Delta\phi$ or $\alpha$ distributions
at fixed $b_{\perp}$), or between $P_{\perp}$ and $b_{\perp}$. From the
unified framework, one sees that there are a series of non-trivial angular
correlations among $q_{\perp}$, $b_{\perp}$, and $P_{\perp}$. We find that
these anisotropies can be divided into two categories. The first category,
which does not depend on $P_{\perp}$, arises from the contraction between the
$O^{ijkl}$ term of hard factor and the photon Wigner distributions. It leads
to two sizable anisotropies, $2\langle\cos(2\phi_{b}-2\phi_{q})\rangle$ and
$2\langle\cos(4\phi_{b}-4\phi_{q})\rangle$. The other category, which involves
the $\cos(4\phi_{P})$ correlation, is due to the contraction between the
$\Omega^{ijkl}$ term and the photon Wigner distributions. We can reproduce the
two anisotropies, $2\langle\cos(4\phi_{P}-4\phi_{q})\rangle$ and
$2\langle\cos(4\phi_{P}-4\phi_{b})\rangle$, as proposed earlier in Refs.
Li:2019yzy ; Li:2019sin ; Xiao:2020ddm . Also, we find that there are
additional novel angular correlations, such as
$2\langle\cos(2\phi_{q}+2\phi_{b}-4\phi_{P})\rangle$,
$2\langle\cos(6\phi_{q}-2\phi_{b}-4\phi_{P})\rangle$ and
$2\langle\cos(6\phi_{b}-2\phi_{q}-4\phi_{P})\rangle$, etc. It is interesting
to note that the impact parameter vector $\vec{b}_{\perp}$ can be determined
through the reaction plane of the hadronic events in peripheral collisions.
Therefore, combined with the additional information provided by the momentum
anisotropy of the produced hadrons in the same event, it is possible to
measure all the aforementioned dilepton angular correlations among
$q_{\perp}$, $P_{\perp}$, and $b_{\perp}$ in peripheral collisions at RHIC and
the LHC.

To provide a complete description of the most up-to-date data, we also need to
include the incoherent contributions from the interior of the colliding
nuclei. When the transverse momentum of incoming photons is less than
$1/R_{A}$ with $R_{A}$ the nuclear radius, one can view the nucleus as a
coherent charge source of photons. While incoming photons emerge incoherently
from the interior of the heavy nuclei, they often carry much larger transverse
momentum as they resolve the nucleus’ internal sub-structure. This scenario is
often accompanied by the dissociation of these nuclei. The STAR, CMS and ATLAS
Collaborations have recently published a substantial amount of data related to
neutron emissions and nuclear break-ups in UPCs STAR:2019wlg ; ATLAS:2020epq ;
CMS:2020skx ; ATLAS:2022srr . In the following, we examine the role of
incoherent contributions in events involving neutron emissions.

In the $0n0n$ scenario in which no neutron is detected, the heavy nucleus
remains intact. The primary contribution comes from the coherent photons. As
substantiated in previous studies Klein:2018fmp ; Klein:2020jom ; Shao:2023zge
, the large transverse momentum imbalance of the lepton pair induced by the
final-state soft multi-photon radiation can be effectively addressed through
the Sudakov resummation. The detailed implementation is provided in the
supplemental material. On the other hand, for the $XnXn$ cases where neutron
emissions are observed, the heavy nucleus is considered disintegrated,
indicating the presence of not only the final-state soft multi-photon
radiation but also incoherent interactions in these events. Therefore, taking
into account photons come from different sources, we modified the GTMD photon
distribution $\Gamma^{ij}(x,k_{T},\Delta)$ as follows

$\displaystyle\Gamma^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$
$\displaystyle=$
$\displaystyle\Gamma_{A}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})+\Gamma_{p}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$
(12) $\displaystyle+\Gamma_{q}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp}).$

The first term $\Gamma_{A}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$
represents the coherent contribution from the heavy nucleus, acting as a
whole, as the photon source. The second term
$\Gamma_{p}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})$ signifies that the
proton inside the nucleus plays the role of the individual photon source,
which is given by

$\Gamma_{p}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=\frac{Z\alpha_{\rm
em}}{\pi^{2}}k^{i}_{\perp}{k^{\prime}}^{j}_{\perp}\frac{F_{p}(k)}{k^{2}}\frac{F_{p}(k^{\prime})}{{k^{\prime}}^{2}},$
(13)

where $F_{p}(k)$ stands for the proton charge form factor. The last term
corresponds to the incoherent contribution originated from quarks inside
nucleons. By treating quarks as point-like particles, then the corresponding
GTMD photon distribution
$\Gamma_{q}^{ij}(x_{,}\vec{k}_{T},\vec{\Delta}_{\perp})$ reads

$\Gamma_{q}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=\alpha_{\rm
em}\frac{\sum_{q}e_{q}^{2}}{\pi^{2}}k^{i}_{\perp}{k^{\prime}}^{j}_{\perp}\frac{1}{k^{2}}\frac{1}{{k^{\prime}}^{2}},$
(14)

where $e_{q}$ is the quark charge number. In the experimental measurements,
the maximum value of transverse momentum imbalance is approximately $1.2$ GeV,
which is already much larger than the typical confinement scale. Thus, at a
scale around $1$ GeV, dilepton production is capable of resolving quarks
within the proton, although fluctuations inside the nucleon remain relatively
insignificant. Consequently, for the sake of simplifying the calculation, we
utilize the valence quark model to outline the quark distribution within both
the proton and the neutron. We assume that these three quarks, two up quarks
and one down quark, share the same momentum within the proton. Similarly, two
down quarks and one up quark share the same momentum within the neutron. The
above approximation can be improved by following a recent study of the photon
parton distribution in the proton Manohar:2016nzj ; Manohar:2017eqh , where a
rigorous connection to the proton structure functions played an important
role. We expect a similar treatment can be made for the above incoherent
contribution to the photon Wigner distribution. We will leave this issue for a
later publication.

Figure 1: Comparison the numerical results with the ATLAS measurement of the
dimuon pair acoplanarity $\alpha$ distribution in central collisions
ATLAS:2022yad .

Figure 2: Comparisons of the numerical results with the CMS $XnXn$ data
CMS:2020skx  and the ATLAS $MnMn$ data ATLAS:2022srr .

Comparing with data In our numerical calculation, we follow the STARlight
Klein:2016yzr  and use the nucleus form factorDavies:1976zzb

$F_{A}(k)=\frac{3}{(kR_{A})^{3}}\frac{\left[\sin(kR_{A})-kR_{A}\cos(kR_{A})\right]}{1+a^{2}k^{2}},$
(15)

with $a=0.71$ fm, $R_{A}$ nucleus radius and $A$ the number of nucleons in the
nucleus. We set $R_{A}=6.6$ fm for the lead nucleus. For the proton, we employ
the dipole form factor Andivahis:1994rq ; Klein:2003vd

$F_{p}(k)=\frac{1}{[(k^{2}/0.71\text{ GeV}^{2})+1]^{2}}.$ (16)

As depicted in Fig. 1, we compare our results with the $\alpha$ distribution
data from the ATLAS Collaboration ATLAS:2022yad  in the $0-5\%$ centrality
interval. The acoplanarity $\alpha$ is defined by $\alpha=1-\Delta\phi/\pi$,
where $\Delta\phi$ denotes the angular difference between two leptons. The
magenta dashed line represents the computation derived from the first term
$O^{ijkl}$ of the hard factor, whereas the black solid line indicates the
calculation resulting from both the first term $O^{ijkl}$ and the second term
$\Omega^{ijkl}(\vec{P}_{\perp})$. In the small $\alpha$ region, there is a
notable gap between the result computed from the first term $O^{ijkl}$ of the
hard factor and the experimental data. However, after including the second
term $\Omega^{ijkl}(\vec{P}_{\perp})$, our calculation can reproduce the dip
shape observed in the experiment and perfectly aligns with the ATLAS
experimental data, thereby highlighting the importance of the
$\Omega^{ijkl}(\vec{P}_{\perp})$ term in describing the dip in the extremely
small $\alpha$ region. Similar numerical results can also be found in the
calculations in Refs. Zha:2018tlq ; Li:2019yzy ; Li:2019sin ; Wang:2021kxm ;
Wang:2022gkd ; Shao:2023zge . As $\alpha$ increases, we find that the
contribution of the $\Omega^{ijkl}(\vec{P}_{\perp})$ term tends to diminish,
and the $O^{ijkl}$ term takes the dominant contribution.

We compare our numerical calculations for the $\alpha$ distribution involving
neutron emissions with the CMS $XnXn$ and ATLAS $MnMn$ data CMS:2020skx ;
ATLAS:2022srr  in Fig. 2. In the CMS experiment, $Xn$ represents the emission
of two or more neutrons, whereas $Mn$ denotes the emission of one or more
neutrons in the ATLAS experiment. The magenta dash-dotted line stands for the
coherent calculation without the Sudakov effect, and the blue dashed line
corresponds to the coherent calculation with the Sudakov effect. In addition,
the orange dotted line represents the sum of the coherent and incoherent
contribution without the Sudakov effect. At last, the solid line includes both
incoherent contribution and the Sudakov effect.

All four calculations agree in the small-$\alpha$ region, meaning that the
coherent term originating from the nucleus acting as the photon source is the
dominant contribution, and the incoherent contribution and Sudakov effect are
negligible in this region. In the large-$\alpha$ region, both the Sudakov
effect and the incoherent contribution are important mechanisms leading to the
final $q_{\perp}$ broadening. However, neither of these mechanisms alone can
accurately describe the data with neutron emissions. When both mechanisms are
considered, we can obtain the enhanced broadening effect and find good
agreements with the CMS and ATLAS data. This implies that the final
$q_{\perp}$ broadening in the large-$\alpha$ region receive contributions from
both the Sudakov effect and the incoherent contribution. In the small-$\alpha$
region, the corresponding transverse momentum imbalance,
$q_{\perp}\sim\pi\alpha P_{\perp}$, is less than $30$ MeV. As a result,
dileptons cannot penetrate or “see” the inner structure of heavy nuclei,
allowing the heavy nucleus to be perceived as a whole entity radiating
photons. However, as $\alpha$ increases, $q_{\perp}$ can rise up to $1\sim 3$
GeV. This enables the dilepton to resolve the inner structure of heavy nuclei,
such as protons and quarks. Consequently, quarks and protons begin to act as
incoherent photon sources in the large $\alpha$ region.

Conclusion In this paper, we present a unified factorized framework, involving
the photon Wigner distribution, to study the lepton pair production via the
photon fusion process in heavy-ion collisions. Furthermore, we show that
various correlations exist between lepton’s transverse momentum $P_{\perp}$,
the imbalance of transverse momentum $q_{\perp}$, and the impact parameter
$b_{\perp}$. Future experimental studies can explore these anisotropies in
heavy-ion collisions at RHIC and the LHC. Additionally, our study shows that
the photon source is not solely from the heavy nucleus as a whole but also
from the protons inside the nucleus and the substructure of the nucleons (the
quarks). Considering events accompanied by neutron emissions measured at the
LHC, we show that the incoherent contributions are significant in the
large-$\alpha$ region. Our unified factorized formula has successfully
described a wide range of relevant data from the STAR, CMS, and ATLAS
Collaborations. Finally, this unified factorization formula, developed in this
study, can offer novel perspectives for exploring the factorization and
production mechanisms of other processes in UPCs.

Acknowledgments We thank Feng Yuan, Jian Zhou, Cheng Zhang, Ya-Jin Zhou and
Shi Pu for useful inputs and discussions. This work is partly supported by the
CUHK-Shenzhen under grant No. UDF01001859. Y. Shi is supported by the China
Postdoctoral Science Foundation under Grant No. 2022M720082. S.Y. Wei is
supported by the Shandong Province Natural Science Foundation under grant No.
2023HWYQ-011 and the Taishan fellowship of Shandong Province for junior
scientists.

## References

  * (1) J. Adams et al. [STAR], Phys. Rev. C 70, 031902 (2004) [arXiv:nucl-ex/0404012 [nucl-ex]].
  * (2) J. Adam et al. [STAR], Phys. Rev. Lett. 121, no.13, 132301 (2018) [arXiv:1806.02295 [hep-ex]].
  * (3) J. Adam et al. [STAR], Phys. Rev. Lett. 127, no.5, 052302 (2021) [arXiv:1910.12400 [nucl-ex]].
  * (4) [ATLAS], ATLAS-CONF-2016-025.
  * (5) M. Aaboud et al. [ATLAS], Phys. Rev. Lett. 121, no.21, 212301 (2018) [arXiv:1806.08708 [nucl-ex]].
  * (6) G. Aad et al. [ATLAS], Phys. Rev. C 104, 024906 (2021) [arXiv:2011.12211 [nucl-ex]].
  * (7) A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 127, no.12, 122001 (2021) [arXiv:2011.05239 [hep-ex]].
  * (8) S. Acharya et al. [ALICE], JHEP 06, 024 (2023) [arXiv:2204.11732 [nucl-ex]].
  * (9) G. Aad et al. [ATLAS], Phys. Rev. C 107, no.5, 054907 (2023) [arXiv:2206.12594 [nucl-ex]].
  * (10) G. Aad et al. [ATLAS], JHEP 2306, 182 (2023) [arXiv:2207.12781 [nucl-ex]].
  * (11) G. Aad et al. [ATLAS], Phys. Rev. Lett. 131, no.15, 151802 (2023) [arXiv:2204.13478 [hep-ex]].
  * (12) A. Tumasyan et al. [CMS], Phys. Rev. Lett. 131, 151803 (2023) [arXiv:2206.05192 [nucl-ex]].
  * (13) S. J. Brodsky, T. Kinoshita and H. Terazawa, Phys. Rev. D 4, 1532-1557 (1971).
  * (14) H. Terazawa, Rev. Mod. Phys. 45, 615-662 (1973).
  * (15) V. M. Budnev, I. F. Ginzburg, G. V. Meledin and V. G. Serbo, Phys. Rept. 15, 181-281 (1975).
  * (16) C. A. Bertulani and G. Baur, Phys. Rept. 163, 299 (1988).
  * (17) R. N. Cahn and J. D. Jackson, Phys. Rev. D 42, 3690-3695 (1990).
  * (18) C. Bottcher, M. R. Strayer and M. J. Rhoades- Brown, Phys. Rev. A 44, 4769-4770 (1991).
  * (19) M. Vidovic, M. Greiner, C. Best and G. Soff, Phys. Rev. C 47, 2308-2319 (1993).
  * (20) K. Hencken, D. Trautmann and G. Baur, Phys. Rev. A 51, 1874-1882 (1995) [arXiv:nucl-th/9410014 [nucl-th]].
  * (21) F. Krauss, M. Greiner and G. Soff, Prog. Part. Nucl. Phys. 39, 503-564 (1997)
  * (22) M. C. Guclu, J. Li, A. S. Umar, D. J. Ernst and M. R. Strayer, Annals Phys. 272, 7-48 (1999)
  * (23) G. Baur, K. Hencken, D. Trautmann, S. Sadovsky and Y. Kharlov, Phys. Rept. 364, 359-450 (2002) [arXiv:hep-ph/0112211 [hep-ph]].
  * (24) A. J. Baltz, S. R. Klein and J. Nystrand, Phys. Rev. Lett. 89, 012301 (2002) [arXiv:nucl-th/0205031 [nucl-th]].
  * (25) G. Baur, K. Hencken, A. Aste, D. Trautmann and S. R. Klein, Nucl. Phys. A 729, 787-808 (2003) [arXiv:nucl-th/0307031 [nucl-th]].
  * (26) K. Hencken, G. Baur and D. Trautmann, Phys. Rev. C 69, 054902 (2004) [arXiv:nucl-th/0402061 [nucl-th]].
  * (27) C. A. Bertulani, S. R. Klein and J. Nystrand, Ann. Rev. Nucl. Part. Sci. 55, 271-310 (2005) [arXiv:nucl-ex/0502005 [nucl-ex]].
  * (28) G. Baur, K. Hencken and D. Trautmann, Phys. Rept. 453, 1-27 (2007) [arXiv:0706.0654 [nucl-th]].
  * (29) A. J. Baltz, G. Baur, D. d’Enterria, L. Frankfurt, F. Gelis, V. Guzey, K. Hencken, Y. Kharlov, M. Klasen and S. R. Klein, et al. Phys. Rept. 458, 1-171 (2008) [arXiv:0706.3356 [nucl-ex]].
  * (30) A. J. Baltz, Y. Gorbunov, S. R. Klein and J. Nystrand, Phys. Rev. C 80, 044902 (2009) [arXiv:0907.1214 [nucl-ex]].
  * (31) S. R. Klein, Phys. Rev. C 97, no.5, 054903 (2018) [arXiv:1801.04320 [nucl-th]].
  * (32) S. Klein, A. H. Mueller, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 122, no.13, 132301 (2019) [arXiv:1811.05519 [hep-ph]].
  * (33) M. Kłusek-Gawenda, R. Rapp, W. Schäfer and A. Szczurek, Phys. Lett. B 790, 339-344 (2019) [arXiv:1809.07049 [nucl-th]].
  * (34) M. I. Vysotskii and E. Zhemchugov, Phys. Usp. 62, no.9, 910-919 (2019) [arXiv:1806.07238 [hep-ph]].
  * (35) W. Zha, J. D. Brandenburg, Z. Tang and Z. Xu, Phys. Lett. B 800, 135089 (2020) [arXiv:1812.02820 [nucl-th]].
  * (36) C. Azevedo, V. P. Gonçalves and B. D. Moreira, Eur. Phys. J. C 79, no.5, 432 (2019) [arXiv:1902.00268 [hep-ph]].
  * (37) C. Li, J. Zhou and Y. J. Zhou, Phys. Lett. B 795, 576-580 (2019) [arXiv:1903.10084 [hep-ph]].
  * (38) C. Li, J. Zhou and Y. J. Zhou, Phys. Rev. D 101, no.3, 034015 (2020) [arXiv:1911.00237 [hep-ph]].
  * (39) S. Klein, A. H. Mueller, B. W. Xiao and F. Yuan, Phys. Rev. D 102, no.9, 094013 (2020) [arXiv:2003.02947 [hep-ph]].
  * (40) S. Klein and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 70, 323-354 (2020) [arXiv:2005.01872 [nucl-ex]].
  * (41) B. W. Xiao, F. Yuan and J. Zhou, Phys. Rev. Lett. 125, no.23, 232301 (2020) [arXiv:2003.06352 [hep-ph]].
  * (42) S. Karadağ and M. C. Güçlü, Phys. Rev. C 102, no.1, 014904 (2020) [arXiv:1911.08507 [hep-ph]].
  * (43) M. Kłusek-Gawenda, W. Schäfer and A. Szczurek, Phys. Lett. B 814, 136114 (2021) [arXiv:2012.11973 [hep-ph]].
  * (44) J. D. Brandenburg, W. Li, L. Ruan, Z. Tang, Z. Xu, S. Yang and W. Zha, [arXiv:2006.07365 [hep-ph]].
  * (45) W. Zha, J. D. Brandenburg, L. Ruan, Z. Tang and Z. Xu, Phys. Rev. D 103, no.3, 033007 (2021) [arXiv:2006.12099 [hep-ph]].
  * (46) Y. Hatta, B. W. Xiao, F. Yuan and J. Zhou, Phys. Rev. D 104, no.5, 054037 (2021) [arXiv:2106.05307 [hep-ph]].
  * (47) W. Zha and Z. Tang, JHEP 08, 083 (2021) [arXiv:2103.04605 [hep-ph]].
  * (48) R. j. Wang, S. Pu and Q. Wang, Phys. Rev. D 104, no.5, 056011 (2021) [arXiv:2106.05462 [hep-ph]].
  * (49) X. Wang, J. D. Brandenburg, L. Ruan, F. Shao, Z. Xu, C. Yang and W. Zha, Phys. Rev. C 107, no.4, 044906 (2023) [arXiv:2207.05595 [nucl-th]].
  * (50) R. j. Wang, S. Lin, S. Pu, Y. f. Zhang and Q. Wang, Phys. Rev. D 106, no.3, 034025 (2022) [arXiv:2204.02761 [hep-ph]].
  * (51) J. D. Brandenburg, J. Seger, Z. Xu and W. Zha, Rept. Prog. Phys. 86, no.8, 083901 (2023) [arXiv:2208.14943 [hep-ph]].
  * (52) S. Lin, R. J. Wang, J. F. Wang, H. J. Xu, S. Pu and Q. Wang, Phys. Rev. D 107, no.5, 054004 (2023) [arXiv:2210.05106 [hep-ph]].
  * (53) D. Y. Shao, C. Zhang, J. Zhou and Y. J. Zhou, Phys. Rev. D 107, no.3, 036020 (2023) [arXiv:2212.05775 [hep-ph]].
  * (54) J. Luo, X. Li, Z. Tang, X. Wu and W. Zha, Phys. Rev. C 108, no.5, 054906 (2023) [arXiv:2308.03070 [hep-ph]].
  * (55) D. Y. Shao, C. Zhang, J. Zhou and Y. j. Zhou, Phys. Rev. D 108, no.11, 116015 (2023) [arXiv:2306.02337 [hep-ph]].
  * (56) J. P. Dai and S. Zhao, Phys. Rev. D 109, no.5, 054022 (2024) [arXiv:2401.04681 [hep-ph]].
  * (57) C. Zhang, L. M. Zhang and D. Y. Shao, [arXiv:2406.05618 [hep-ph]].
  * (58) S. R. Klein, J. Nystrand, J. Seger, Y. Gorbunov and J. Butterworth, Comput. Phys. Commun. 212, 258-268 (2017) [arXiv:1607.03838 [hep-ph]].
  * (59) L. A. Harland-Lang, V. A. Khoze and M. G. Ryskin, Eur. Phys. J. C 79, no.1, 39 (2019) [arXiv:1810.06567 [hep-ph]].
  * (60) N. Burmasov, E. Kryshen, P. Buehler and R. Lavicka, Comput. Phys. Commun. 277, 108388 (2022) [arXiv:2111.11383 [hep-ph]].
  * (61) H. S. Shao and D. d’Enterria, JHEP 09, 248 (2022) [arXiv:2207.03012 [hep-ph]].
  * (62) L. A. Harland-Lang, Phys. Rev. D 107, no.9, 093004 (2023) [arXiv:2303.04826 [hep-ph]].
  * (63) J. Abdallah et al. [DELPHI], Eur. Phys. J. C 35, 159-170 (2004) [arXiv:hep-ex/0406010 [hep-ex]].
  * (64) S. Knapen, T. Lin, H. K. Lou and T. Melia, Phys. Rev. Lett. 118, no.17, 171801 (2017) [arXiv:1607.06083 [hep-ph]].
  * (65) A. M. Sirunyan et al. [CMS], Phys. Lett. B 797, 134826 (2019) [arXiv:1810.04602 [hep-ex]].
  * (66) J. Ellis, N. E. Mavromatos and T. You, Phys. Rev. Lett. 118, no.26, 261802 (2017) [arXiv:1703.08450 [hep-ph]].
  * (67) I. Xu, N. Lewis, X. Wang, J. D. Brandenburg and L. Ruan, [arXiv:2211.02132 [hep-ex]].
  * (68) D. Y. Shao, B. Yan, S. R. Yuan and C. Zhang, [arXiv:2310.14153 [hep-ph]].
  * (69) G. Baur, K. Hencken and D. Trautmann, J. Phys. G 24, 1657-1692 (1998) [arXiv:hep-ph/9804348 [hep-ph]].
  * (70) A. J. Baltz and M. Strikman, Phys. Rev. D 57, 548-549 (1998) [arXiv:hep-ph/9705220 [hep-ph]].
  * (71) P. A. Steinberg [ALICE, ATLAS, CMS, LHCb and STAR], Nucl. Phys. A 1005, 122007 (2021)
  * (72) K. Hattori, H. Taya and S. Yoshida, JHEP 01, 093 (2021) [arXiv:2010.13492 [hep-ph]].
  * (73) P. Copinger and S. Pu, Int. J. Mod. Phys. A 35, no.28, 2030015 (2020) [arXiv:2008.03635 [hep-ph]].
  * (74) J. D. Brandenburg, W. Zha and Z. Xu, Eur. Phys. J. A 57, no.10, 299 (2021) [arXiv:2103.16623 [hep-ph]].
  * (75) S. R. Klein and H. Mäntysaari, Nature Rev. Phys. 1, no.11, 662-674 (2019) [arXiv:1910.10858 [hep-ex]].
  * (76) H. Mäntysaari, Rept. Prog. Phys. 83, no.8, 082201 (2020) [arXiv:2001.10705 [hep-ph]].
  * (77) H. Mäntysaari and R. Venugopalan, Phys. Lett. B 781, 664-671 (2018) [arXiv:1712.02508 [nucl-th]].
  * (78) A. Morreale and F. Salazar, Universe 7, no.8, 312 (2021) [arXiv:2108.08254 [hep-ph]].
  * (79) H. Mäntysaari, F. Salazar, B. Schenke, C. Shen and W. Zhao, Phys. Rev. C 109, no.2, 024908 (2024) [arXiv:2310.15300 [nucl-th]].
  * (80) H. Mäntysaari, F. Salazar and B. Schenke, [arXiv:2312.04194 [hep-ph]].
  * (81) Y. V. Kovchegov, H. Sun and Z. Tu, [arXiv:2311.12208 [hep-ph]].
  * (82) J. Cepila, J. G. Contreras, M. Matas and A. Ridzikova, [arXiv:2312.11320 [hep-ph]].
  * (83) S. Acharya et al. [ALICE], Eur. Phys. J. C 81, no.8, 712 (2021) [arXiv:2101.04577 [nucl-ex]].
  * (84) A. Tumasyan et al. [CMS], Phys. Rev. Lett. 131, no.26, 262301 (2023) [arXiv:2303.16984 [nucl-ex]].
  * (85) R. Aaij et al. [LHCb], JHEP 06, 146 (2023) [arXiv:2206.08221 [hep-ex]].
  * (86) S. Acharya et al. [ALICE], JHEP 10, 119 (2023) [arXiv:2305.19060 [nucl-ex]].
  * (87) [STAR], [arXiv:2311.13637 [nucl-ex]].
  * (88) C. F. von Weizsacker, Z. Phys. 88, 612-625 (1934).
  * (89) E. J. Williams, Phys. Rev. 45, 729-730 (1934).
  * (90) X. D. Ji, Phys. Rev. Lett. 78, 610-613 (1997) [arXiv:hep-ph/9603249 [hep-ph]].
  * (91) A. V. Belitsky and D. Mueller, Phys. Lett. B 486, 369-377 (2000) [arXiv:hep-ph/0005028 [hep-ph]].
  * (92) M. Diehl, Eur. Phys. J. C 19, 485-492 (2001) [arXiv:hep-ph/0101335 [hep-ph]].
  * (93) X. d. Ji, Phys. Rev. Lett. 91, 062001 (2003) [arXiv:hep-ph/0304037 [hep-ph]].
  * (94) A. V. Belitsky, X. d. Ji and F. Yuan, Phys. Rev. D 69, 074014 (2004) [arXiv:hep-ph/0307383 [hep-ph]].
  * (95) Y. Hatta, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 116, no.20, 202301 (2016) [arXiv:1601.01585 [hep-ph]].
  * (96) A. Manohar, P. Nason, G. P. Salam and G. Zanderighi, Phys. Rev. Lett. 117, no.24, 242002 (2016) [arXiv:1607.04266 [hep-ph]].
  * (97) A. V. Manohar, P. Nason, G. P. Salam and G. Zanderighi, JHEP 12, 046 (2017) [arXiv:1708.01256 [hep-ph]].
  * (98) K. T. R. Davies and J. R. Nix, Phys. Rev. C 14, 1977-1994 (1976)
  * (99) L. Andivahis, P. E. Bosted, A. Lung, L. M. Stuart, J. Alster, R. G. Arnold, C. C. Chang, F. S. Dietrich, W. R. Dodge and R. Gearhart, et al. Phys. Rev. D 50, 5491-5517 (1994)
  * (100) S. R. Klein and J. Nystrand, Phys. Rev. Lett. 92, 142003 (2004) [arXiv:hep-ph/0311164 [hep-ph]].

## Supplemental material

###### Contents

  1. .1 The Generalized Factorization formalism
  2. .2 Phenomenological Studies
    1. .2.1 Anisotropies
    2. .2.2 The Sudakov Effect
    3. .2.3 Dilepton Events with Neutrons Emission
    4. .2.4 The Incoherent Contribution to the Photon Distribution
  3. .3 Comparisons with the RHIC and LHC data
    1. .3.1 Acoplanarity and Transverse Momentum Broadening
    2. .3.2 The Neutron Emission effect in the Small $\alpha$ Region
    3. .3.3 Enhancements from the Incoherent Contributions in the Large $\alpha$ Region

We provide additional details in the supplemental material. In Sec. .1, the
generalized factorization framework for the dilepton production via the
photon-photon fusion process in ultra-peripheral heavy-ion collisions is
presented. Sec. .2 discusses several topics concerning the phenomenological
studies. In Sec. .3, we present additional numerical results compared with the
experimental data.

###  .1 The Generalized Factorization formalism

Figure 3: Left: The cartoon depiction of the lepton pair production from the
photon-photon fusion process in heavy-ion collisions. Right: The Feynman
diagram illustrating the generalized TMD factorization for the
$\gamma(k_{1})+\gamma(k_{2})\rightarrow l^{+}(p_{1})+l^{-}(p_{2})$ process in
heavy-ion collisions.

The left figure in Fig. 3 shows the cartoon of the dilepton production from
electromagnetic interactions in heavy-ion collisions: two quasi-real photons,
each radiated from one of the two fast-moving nuclei at a relative distance
$b_{\perp}$, undergo fusion to generate a pair of leptons. These leptons are
almost back-to-back, each with transverse momentum $p_{1\perp}$ and
$p_{2\perp}$, and are produced at the rapidity $y_{1}$ and $y_{2}$,
respectively. The collision point, where these leptons are produced, is
situated at distances $b_{1\perp}$ and $b_{2\perp}$ from the two heavy nuclei.

In the following calculations, we assume that the relative transverse momentum
of the lepton pair $\vec{P}_{\perp}=(\vec{p}_{1\perp}-\vec{p}_{2\perp})/2\sim
p_{1\perp}\sim p_{2\perp}$ is much large than the total transverse momentum of
lepton pair $\vec{q}_{\perp}=\vec{p}_{1\perp}+\vec{p}_{2\perp}$. Thus, we can
use the relative transverse momentum of the lepton pair $P_{\perp}$ as the
largest scale to factorize the process into two different parts, the hard part
and the photon distribution part. The hard part corresponds to the hard
interaction of lepton pair production, which depends on the hard momentum
$P_{\perp}$ and its azimuthal angle $\phi_{P}$. The photon distributions
involve soft momentum scales, which are of the order of $q_{\perp}$ much less
than $P_{\perp}$.

The study of photon distributions within the nucleus has a long history.
Initially, the equivalent photon approximation (EPA), also known as the
Weizsäcker-Williams (WW) approximation method vonWeizsacker:1934nji ;
Williams:1934ad , is the most commonly used method for addressing photon
distribution. In this method, photons are treated as quasi-real particles, and
by calculating the lowest-order contributions in QED, one can obtain a photon
distribution with low virtuality and transverse momentum. However, the
original approximation does not account for both the transverse momentum and
the impact parameter of the photon. Recently, inspired by the implementation
of the Wigner distribution in QCD studies Ji:2003ak ; Belitsky:2003nz ;
Hatta:2016dxp , the generalized equivalent photon approximation including the
transverse momentum and the impact parameter dependence has been proposed to
parameterize the photon distribution from high energy nuclei, which is also
known as the photon Wigner distribution Li:2019sin ; Klein:2020jom . Also, one
can define the Generalized Transverse Momentum Dependent (GTMD) distribution
by Fourier transforming the impact parameter dependence into the momentum
space.

As illustrated in the right figure of Fig. 3, the generalized factorization
formula for the process following the GTMD in the correlation limit can be
cast into

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}$
$\displaystyle=$ $\displaystyle\frac{1}{(2\pi)^{4}}\int
d^{2}b_{1\perp}d^{2}b_{2\perp}\delta^{2}(\vec{b}_{\perp}-\vec{b}_{1\perp}+\vec{b}_{2\perp})\int
d^{2}k_{1T}d^{2}k_{2T}\int
d^{2}\Delta_{1\perp}d^{2}\Delta_{2\perp}e^{-ib_{1\perp}\cdot\Delta_{1\perp}-ib_{2\perp}\cdot\Delta_{2\perp}}$
(17)
$\displaystyle\times\Gamma_{ij}(x_{1},\vec{k}_{1T},\vec{\Delta}_{1\perp})\Gamma_{kl}(x_{2},\vec{k}_{2T},\vec{\Delta}_{2\perp})\mathcal{H}^{ijkl}(\vec{P}_{\perp})\delta^{(2)}(\vec{q}_{\perp}-\vec{k}_{1T}-\vec{k}_{2T}),$

where $\mathcal{H}^{ijkl}$ stands for the hard factor, $\Gamma_{ij}$ and
$\Gamma_{kl}$ represent the 5-dimensional GTMD distribution with the
polarization indices $i$ and $k$ of incoming two photons in the scattering
amplitude and the polarization indices $j$ and $l$ of photons in conjugate
amplitude. These two-dimensional indices run from $1$ to $2$. To reconstruct
the kinematics of the incoming photons, we write the longitudinal momentum
fraction of photon per nucleon as
$x_{1}=\text{max}[p_{1\perp},p_{2\perp}](e^{y_{1}}+e^{y_{2}})/\sqrt{s}$ and
$x_{2}=\text{max}[p_{1\perp},p_{2\perp}](e^{-y_{1}}+e^{-y_{2}})/\sqrt{s}$.
$b_{1\perp}$ and $b_{2\perp}$ stand for the distance of the photon and its
corresponding nucleus.
$\vec{k}_{1T}=(\vec{k}_{1\perp}+\vec{k}_{1\perp}^{\prime})/2$ and
$\vec{k}_{2T}=(\vec{k}_{2\perp}+\vec{k}_{2\perp}^{\prime})/2$ are transverse
momenta of two photons, where $k_{i\perp}$ and $k_{i\perp}^{\prime}$ are the
transverse momenta of incoming photons in the amplitude and conjugate
amplitude, respectively. According the momentum conservation, these momenta
satisfy
$\vec{q}_{\perp}=\vec{k}_{1\perp}+\vec{k}_{2\perp}=\vec{k}_{1\perp}^{\prime}+\vec{k}_{2\perp}^{\prime}$.
$\vec{\Delta}_{1\perp}=\vec{k}_{1\perp}^{\prime}-\vec{k}_{1\perp}$ and
$\vec{\Delta}_{2\perp}=\vec{k}^{\prime}_{2\perp}-\vec{k}_{2\perp}$ are defined
as the momentum difference for the nucleus state in the non-forward
scatterings and the momentum conservation implies
$\vec{\Delta}_{1\perp}+\vec{\Delta}_{2\perp}=0$. Thus we define
$\vec{\Delta}_{\perp}=\vec{\Delta}_{1\perp}=-\vec{\Delta}_{2\perp}$ for
convenience. It is important to note that the momenta from the scattering
amplitude and conjugate amplitude only become identical in forward scattering
with $\Delta_{1\perp}=\Delta_{2\perp}=0$. By integrating over the two delta
functions, the cross-section can be rewritten as

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}$
$\displaystyle=$ $\displaystyle\frac{1}{(2\pi)^{2}}\int
d^{2}k_{1\perp}d^{2}k_{1\perp}^{\prime}e^{ib_{\perp}\cdot\Delta_{\perp}}\Gamma_{ij}(x_{1},\vec{k}_{1T},\vec{\Delta}_{\perp})\Gamma_{kl}(x_{2},\vec{k}_{2T},-\vec{\Delta}_{\perp})\mathcal{H}^{ijkl}(\vec{P}_{\perp}),$
(18)

with
$\vec{\Delta}_{\perp}=\vec{k}_{1\perp}^{\prime}-\vec{k}_{1\perp}=\vec{k}_{2\perp}-\vec{k}_{2\perp}^{\prime}$.
By using the definition

$\Gamma_{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=\int
d^{2}b_{\perp}e^{-i\vec{b}_{\perp}\cdot\vec{\Delta}_{\perp}}\Gamma_{ij}(x,\vec{k}_{T},\vec{b}_{\perp}),$
(19)

where $b_{\perp}$ is the distance of the photon and center of the nucleus, one
can convert the GTMD photon distribution into the photon Wigner distribution
via the Fourier transform. Then, the generalized factorization formula
involving the Winger photon distribution becomes

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}p_{1\perp}d^{2}p_{2\perp}d^{2}b_{\perp}}$
$\displaystyle=$ $\displaystyle\int
d^{2}b_{1\perp}d^{2}b_{2\perp}\delta^{2}(\vec{b}_{\perp}-\vec{b}_{1\perp}+\vec{b}_{2\perp})\int
d^{2}k_{1T}d^{2}k_{2T}$ (20)
$\displaystyle\times\Gamma_{ij}(x_{1},\vec{k}_{1T},\vec{b}_{1\perp})\Gamma_{kl}(x_{2},\vec{k}_{2T},\vec{b}_{2\perp})\mathcal{H}^{ijkl}(\vec{P}_{\perp})\delta^{(2)}(\vec{p}_{1\perp}+\vec{p}_{2\perp}-\vec{k}_{1T}-\vec{k}_{2T}).$

In the meantime, with the definition

$\Gamma_{ij}(x,\vec{k}_{T},\vec{b}_{\perp})=\int
d^{2}r_{\perp}e^{-i\vec{k}_{T}\cdot\vec{r}_{\perp}}\Gamma_{ij}(x,\vec{r}_{\perp},\vec{b}_{\perp}),$
(21)

where $r_{\perp}$ represents photon position difference between the amplitude
and conjugate amplitude. The cross-section in the full coordinate space can be
written as

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}=\int
d^{2}b_{1\perp}d^{2}r_{\perp}e^{-iq_{\perp}\cdot
r_{\perp}}\Gamma_{ij}(x_{1},\vec{r}_{\perp},\vec{b}_{1\perp})\Gamma_{kl}(x_{2},\vec{r}_{\perp},-\vec{b}_{\perp}+\vec{b}_{1\perp})\mathcal{H}^{ijkl}(\vec{P}_{\perp}).$
(22)

In order to derive the hard factor $\mathcal{H}^{ijkl}$, one needs to carry
out the Feynman graph calculation for $\gamma\gamma\to l^{+}l^{-}$ with open
indices for all four photon polarizations. After some tedious tensor
manipulations, one can eventually find

$\displaystyle\mathcal{H}^{ijkl}(\vec{P}_{\perp})$ $\displaystyle=$
$\displaystyle\frac{\alpha_{\rm
em}^{2}}{\hat{s}^{2}}\Bigg{\\{}2\frac{(m_{l}^{4}-\hat{t}\hat{u})(6m_{l}^{4}-4m_{l}^{2}\hat{t}-4m_{l}^{2}\hat{u}+\hat{t}^{2}+\hat{u}^{2})}{(m_{l}^{2}-\hat{t})^{2}(m_{l}^{2}-\hat{u})^{2}}\left(\delta^{ij}\delta^{kl}-\delta^{ik}\delta^{jl}+\delta^{il}\delta^{jk}\right)$
(23)
$\displaystyle\qquad\quad+\frac{4m_{l}^{4}(2m_{l}^{2}-\hat{t}-\hat{u})^{2}}{(m_{l}^{2}-\hat{t})^{2}(m_{l}^{2}-\hat{u})^{2}}\left(-\delta^{ij}\delta^{kl}+\delta^{ik}\delta^{jl}+\delta^{il}\delta^{jk}\right)$
$\displaystyle\qquad\quad-\frac{2m_{l}^{2}(\hat{t}+\hat{u})(2m_{l}^{2}-\hat{t}-\hat{u})^{2}}{(m_{l}^{2}-\hat{t})^{2}(m_{l}^{2}-\hat{u})^{2}}\left(\delta^{ij}\delta^{kl}+\delta^{ik}\delta^{jl}-\delta^{il}\delta^{jk}\right)$
$\displaystyle\qquad\quad-8m_{l}^{2}\frac{(2m_{l}^{2}-\hat{t}-\hat{u})(m_{l}^{4}-\hat{t}\hat{u})}{(m_{l}^{2}-\hat{t})^{2}(m_{l}^{2}-\hat{u})^{2}}\left[\delta^{ij}\Pi^{kl}(\vec{P}_{\perp})+\Pi^{ij}(\vec{P}_{\perp})\delta^{kl}\right]$
$\displaystyle\qquad\quad-4\frac{(m_{l}^{4}-\hat{t}\hat{u})^{2}}{(m_{l}^{2}-\hat{t})^{2}(m_{l}^{2}-\hat{u})^{2}}\Omega^{ijkl}(\vec{P}_{\perp})\Bigg{\\}},$

where $\hat{s}$, $\hat{u}$ and $\hat{t}$ are the Mandelstam variables. Inside
the curly brackets, the first three terms are the isotropic contributions
which are independent of the orientation of $\vec{P}_{\perp}$, the fourth term
gives rise to the $\cos(2\phi_{P})$ angular correlation, and the last term can
generate the $\cos(4\phi_{P})$ angular correlation. The tensors
$\Pi^{ij}(\vec{P}_{\perp})$ and $\Omega^{ijkl}(\vec{P}_{\perp})$ are defined
as
$\Pi^{ij}(\vec{P}_{\perp})=2\frac{P_{\perp}^{i}P_{\perp}^{j}}{P_{\perp}^{2}}-\delta^{ij}$
and
$\Omega^{ijkl}(\vec{P}_{\perp})=2\Pi^{ij}(\vec{P}_{\perp})\Pi^{kl}(\vec{P}_{\perp})-\left(\delta^{il}\delta^{jk}+\delta^{ik}\delta^{jl}-\delta^{ij}\delta^{kl}\right)$.
They can project the cross-sections onto $\cos(2\phi_{P})$ and
$\cos(4\phi_{P})$ azimuthal angular correlations, respectively. For example,
with four arbitrary two-dimensional unit vector $\hat{n}_{1,2,3,4}$, one can
obtain

$\displaystyle\Pi^{ij}\hat{n}_{1i}\hat{n}_{2j}=\cos(\phi_{1}+\phi_{2}-2\phi_{P})\,,$
(24)
$\displaystyle\Omega^{ijkl}\hat{n}_{1i}\hat{n}_{2j}\hat{n}_{3k}\hat{n}_{4l}=\cos(\phi_{1}+\phi_{2}+\phi_{3}+\phi_{4}-4\phi_{P}),$
(25)

where $\phi_{1}$ represents the azimuthal angle orientation of $\hat{n}_{1}$,
and similarly for other azimuthal angles.

If we average over the angle of impact parameter $b_{\perp}$ and project
$\cos(2\phi_{P})$ terms into $\cos(2\phi_{P}-2\phi_{q})$ and $\cos(4\phi_{P})$
term into $\cos(4\phi_{P}-4\phi_{q})$, we can reproduce the result derived in
Refs. Shao:2022stc . Since the contribution of lepton mass is power suppressed
by factors of $m_{l}^{2}/P^{2}_{\perp}$, we neglect the mass of lepton effect
in the hard factor. Then, the hard factor can be significantly simplified and
it reads

$\displaystyle\mathcal{H}^{ijkl}(\vec{P}_{\perp})$ $\displaystyle=$
$\displaystyle\frac{\alpha_{{\rm
em}}^{2}}{\hat{s}^{2}}\left[O^{ijkl}-4\Omega^{ijkl}(\vec{P}_{\perp})\right],$
(26)

with

$\displaystyle O^{ijkl}$ $\displaystyle=$ $\displaystyle
2\left(\frac{\hat{u}}{\hat{t}}+\frac{\hat{t}}{\hat{u}}\right)\left(\delta^{ij}\delta^{kl}-\delta^{ik}\delta^{jl}+\delta^{il}\delta^{jk}\right),$
(27) $\displaystyle\Omega^{ijkl}(\vec{P}_{\perp})$ $\displaystyle=$
$\displaystyle
2\Pi^{ij}(\vec{P}_{\perp})\Pi^{kl}(\vec{P}_{\perp})-\left(\delta^{il}\delta^{jk}+\delta^{ik}\delta^{jl}-\delta^{ij}\delta^{kl}\right).$
(28)

The terms related to $\cos(2\phi_{P})$ are proportional to the lepton’s mass,
thus they vanish in the massless limit. In the remaining terms, the $O^{ijkl}$
term represents the isotropic component with respect to $P_{\perp}$ and
coincides with the angle-averaged $\phi_{P}$ term as discussed in Ref.
Klein:2020jom . This term emerges from the contributions of both the
unpolarized and the linearly-polarized photon distributions. The
$\Omega^{ijkl}(\vec{P}_{\perp})$ term, which originates from the interaction
between two linearly-polarized photons, will generate the $\cos(4\phi_{P})$
angular correlations in the following calculations.

In general, the GTMD photon distribution can be parametrized as follows

$\Gamma^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=\frac{\delta^{ij}}{2}xf_{\gamma}(x,\vec{k}_{T},\vec{\Delta}_{\perp})+\left(\frac{k^{i}_{+}k^{j}_{-}}{\vec{k}_{-}\cdot\vec{k}_{+}}-\frac{\delta^{ij}}{2}\right)xh_{\gamma}(x,\vec{k}_{T},\vec{\Delta}_{\perp}),$
(29)

with $\vec{k}_{\pm}=\vec{k}_{\perp}\pm\vec{\Delta}_{\perp}/2$. The first term
represents the unpolarized photon distribution and the second term stands for
linearly-polarized photon distribution. We can assume that the unpolarized
photon distribution is the same as the linearly-polarized photon distribution,
and the GTMD photon distribution can be parametrized as

$\Gamma^{ij}(x_{1},\vec{k}_{1T},\vec{\Delta}_{\perp})=\frac{Z^{2}\alpha_{\rm
em}}{\pi^{2}}k^{i}_{1\perp}{k^{\prime}}^{j}_{1\perp}\frac{F_{A}(k_{1})}{k_{1}^{2}}\frac{F_{A}(k_{1}^{\prime})}{{k^{\prime}}_{1}^{2}},$
(30)

where $F_{A}(k_{1})$ is the normalized elastic charge form factor for the
nucleus, $Z$ is charge number of nucleus and $m_{p}$ is the mass of proton
with $k_{1}^{2}=x_{1}^{2}m_{p}^{2}+k_{1\perp}^{2}$ and
${k^{\prime}}_{1}^{2}=x_{1}^{2}m_{p}^{2}+{k^{\prime}}_{1\perp}^{2}$. Now, the
cross-section also can be written as

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}$
$\displaystyle=$ $\displaystyle\frac{1}{(2\pi)^{2}}\frac{2\alpha_{{\rm
e}}^{4}}{\hat{s}^{2}}\left(\frac{Z^{2}}{\pi^{2}}\right)^{2}\int
d^{2}k_{1\perp}d^{2}k_{1\perp}^{\prime}e^{ib_{\perp}\cdot\Delta_{\perp}}\frac{F_{A}(k_{1})}{k_{1}^{2}}\frac{F(k_{2})}{k_{2}^{2}}\frac{F_{A}(k_{1}^{\prime})}{{k^{\prime}}_{1}^{2}}\frac{F_{A}(k_{2}^{\prime})}{{k^{\prime}}_{2}^{2}}k_{1\perp}{k^{\prime}}_{1\perp}k_{2\perp}{k^{\prime}}_{2\perp}$
$\displaystyle\times\left[\left(\frac{\hat{u}}{\hat{t}}+\frac{\hat{t}}{\hat{u}}\right)\cos(\phi_{k_{1\perp}}-\phi_{{k^{\prime}}_{1\perp}}+\phi_{k_{2\perp}}-\phi_{{k^{\prime}}_{2\perp}})-2\cos(\phi_{k_{1\perp}}+\phi_{{k^{\prime}}_{1\perp}}+\phi_{k_{2\perp}}+\phi_{{k^{\prime}}_{2\perp}}-4\phi_{P})\right].$

From the above equation, by averaging over the impact parameter angle
$b_{\perp}$ and projecting the $\cos(4\phi_{P})$ term into
$\cos(4\phi_{P}-4\phi_{q})$ direction, we can reproduce the results derived in
Refs. Li:2019sin . In the coordinate space, if we integrate over $q_{\perp}$
and $\phi_{q}$ and implement the delta function track in the anisotropy term,
it is straightforward to find that the above equation is equal to the results
given by Eq. (5) in Ref. Xiao:2020ddm . One can find that all anisotropic
terms with $\phi_{q}$ disappear, and only $\cos(4\phi_{b}-4\phi_{P})$ remains.

The cross-section, indicative of the interaction probability between
particles, must be positively defined. As mentioned in Ref. Klein:2020jom ,
the photon GTMD is not always positive definite. Consequently, it is essential
to demonstrate that dilepton production via photon-photon fusion is invariably
guaranteed to yield a positive definite cross-section. Based on Eq.
(LABEL:eq:cos), it is interesting to show that

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}$
$\displaystyle\geq$
$\displaystyle\frac{1}{(2\pi)^{2}}\left(\frac{Z^{2}}{\pi^{2}}\right)^{2}\int
d^{2}k_{1\perp}d^{2}k_{1\perp}^{\prime}e^{ib_{\perp}\cdot\Delta_{\perp}}\frac{F_{A}(k_{1})}{k_{1}^{2}}\frac{F(k_{2})}{k_{2}^{2}}\frac{F_{A}(k_{1}^{\prime})}{{k^{\prime}}_{1}^{2}}\frac{F_{A}(k_{2}^{\prime})}{{k^{\prime}}_{2}^{2}}k_{1\perp}{k^{\prime}}_{1\perp}k_{2\perp}{k^{\prime}}_{2\perp}$
(32) $\displaystyle\times\frac{4\alpha_{{\rm
e}}^{4}}{\hat{s}^{2}}\left[\cos(\phi_{k_{1\perp}}-\phi_{{k^{\prime}}_{1\perp}}+\phi_{k_{2\perp}}-\phi_{{k^{\prime}}_{2\perp}})-\cos(\phi_{k_{1\perp}}+\phi_{{k^{\prime}}_{1\perp}}+\phi_{k_{2\perp}}+\phi_{{k^{\prime}}_{2\perp}}-4\phi_{P})\right].$

By using the following identities

$\displaystyle
2\sin(\phi_{k_{1\perp}}+\phi_{k_{2\perp}}-2\phi_{P})\sin(\phi_{{k^{\prime}}_{1\perp}}+\phi_{{k^{\prime}}_{2\perp}}-2\phi_{P})$
(33) $\displaystyle=$
$\displaystyle\cos(\phi_{k_{1\perp}}-\phi_{{k^{\prime}}_{1\perp}}+\phi_{k_{2\perp}}-\phi_{{k^{\prime}}_{2\perp}})-\cos(\phi_{k_{1\perp}}+\phi_{{k^{\prime}}_{1\perp}}+\phi_{k_{2\perp}}+\phi_{{k^{\prime}}_{2\perp}}-4\phi_{P}),$

and

$\sin(\phi_{k_{1\perp}}+\phi_{k_{2\perp}}-2\phi_{P})=\sin(\phi_{k_{1\perp}}-\phi_{P})\cos(\phi_{k_{2\perp}}-\phi_{P})+\cos(\phi_{k_{1\perp}}-\phi_{P})\sin(\phi_{k_{2\perp}}-\phi_{P})\,,$
(34)

one can rewrite the cross section as

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}\geq\left(\frac{Z^{2}}{\pi^{2}}\right)^{2}\frac{8\alpha_{{\rm
e}}^{4}}{(2\pi)^{2}\hat{s}^{2}}\left(G^{12}+G^{21}\right)\left(G^{12*}+G^{21*}\right),$
(35)

where we define unit vector $\hat{l}_{i}$ which has the angle
$\phi_{l_{i}}=\phi_{k_{i\perp}}-\phi_{P}$ and $G^{ik}$ is defined as
$G^{ik}=\int d^{2}k_{1\perp}e^{ik_{1\perp}\cdot
b_{\perp}}F_{A}(k_{1})F_{A}(k_{2})|k_{1\perp}||k_{2\perp}|\hat{l_{1}}^{i}\hat{l_{2}}^{k}$.
It is then obvious to see that the cross-section is positive definite, since
the right hand side is essentially related to the square of an amplitude.

###  .2 Phenomenological Studies

Based on the above unified GTMD factorized framework, we discuss
phenomenological studies of dilepton production in this subsection in detail.
First, we elaborate on the various dilepton anisotropies and show the
numerical predictions. Then, we include the Sudakov effect, which resums the
soft photon radiation emitted by the lepton pairs. In addition, we take the
neutron emission effect in ultra-peripheral collisions (UPCs) into account.
Finally, we discuss the incoherent contribution of the photon Wigner
distribution.

####  .2.1 Anisotropies

Figure 4:  The prediction of different dilepton anisotropies as a function of
$b_{\perp}$ with pair mass $M=8-45$ GeV, fixed $q_{\perp}=0.05$ GeV and
$P_{\perp}=4-20$ GeV at the LHC kinematics. The rapidities of the dilepton are
integrated over the regions $[-1,1]$.  Figure 5: Predictions of different
dilepton anisotropies are presented as a function of $q_{\perp}$ for a pair
mass $M=[8,45]$ GeV, with $b_{\perp}$ fixed at $10$ fm and $P_{\perp}$ ranging
from $4$ to $20$ GeV, within the LHC kinematics. The rapidities of the
dileptons are integrated over the interval $[1,1]$.  Figure 6: Under RHIC
kinematics, the anisotropies of dileptons are predicted as functions of
$b_{\perp}$, with the pair mass $M$ ranging from $0.4$ to $0.76$ GeV,
$q_{\perp}$ fixed at $0.05$ GeV, and $P_{\perp}$ between $0.2$ and $2$ GeV.
The rapidities of the dileptons are integrated over the interval $[-1,1]$.
Figure 7:  For RHIC kinematics, the anisotropies of dileptons are predicted as
functions of $q_{\perp}$, with the pair mass $M$ ranging from $0.4$ to $0.76$
GeV, $b_{\perp}$ fixed at $10$ fm, and $P_{\perp}$ between $0.2$ and $2$ GeV.
The rapidities of the dileptons are integrated over the interval $[-1,1]$.

In the following, we explore the anisotropy of the dilepton production process
$\gamma+\gamma\rightarrow l^{+}+l^{-}$. The cross-section contains various
anisotropies. For example, there are a lot of studies about the angular
correlation between $P_{\perp}$ and $b_{\perp}$, and the angular correlation
between $P_{\perp}$ and $q_{\perp}$, which results in the commonly measured
$\cos(4\phi_{P}-4\phi_{q})$ and $\alpha$ distribution, as shown in Refs.
Li:2019yzy ; Li:2019sin ; Xiao:2020ddm ; Klein:2020jom ; Klein:2020fmr ;
Klusek-Gawenda:2020eja ; Brandenburg:2020ozx ; Zha:2020cst ; Shao:2022stc ;
Shao:2023zge .

Based on the expression of the hard factor
$\mathcal{H}^{ijkl}(\vec{P}_{\perp})$ in Eq. (26) and its projection
properties when contracting with the impact parameter ($\vec{b}_{\perp}$) and
the transverse momentum of the final lepton pair ($\vec{q}_{\perp}$), we can
project the cross-section in the following general Fourier harmonic series

$\displaystyle\sigma$ $\displaystyle=$
$\displaystyle\sigma_{0}+\sigma_{2qb}\cos(2\phi_{q}-2\phi_{b})+\sigma_{4qb}\cos(4\phi_{q}-4\phi_{b})+\sigma_{4qP}\cos(4\phi_{q}-4\phi_{P})$
(36)
$\displaystyle+\sigma_{4bP}\cos(4\phi_{b}-4\phi_{P})+\sigma_{qbP}\cos(2\phi_{q}+2\phi_{b}-4\phi_{P})+...,$

where the harmonic coefficients $\sigma_{2qb},\sigma_{4qb},...$ contains the
detailed quantum information of the photon flux. From the GTMD factorization
framework, we find that linearly polarized photons can generate various
anisotropies with considerable magnitudes from the correlation among azimuthal
angles of $q_{\perp}$, $P_{\perp}$, and $b_{\perp}$. Not only can we reproduce
anisotropies studied in the previous works, but we have also found several
other interesting correlations arising from the inner products of $P_{\perp}$,
$q_{\perp}$ and $b_{\perp}$. These anisotropies can be divided into two
categories.

The first category of correlation is independent of $P_{\perp}$. When the
$O^{ijkl}$ term from the hard factor $\mathcal{H}^{ijkl}(\vec{P}_{\perp})$, an
isotropic term associated with $P_{\perp}$, couples with two incoming Wigner
photon distributions, it results in two significant angular correlations:
$\cos(2\phi_{q}-2\phi_{b})$ and $\cos(4\phi_{q}-4\phi_{b})$. More precisely,
when we analyze the cross-section in momentum space, these correlations
originate from the contraction between the hard factor and the linearly
polarized photon distribution via the phase factor
$e^{i\vec{b}_{\perp}\cdot\vec{\Delta}_{\perp}}$. This second category involves
the $\cos(4\phi_{P})$ harmonics, which comes from the tensor contraction
between the $\Omega^{ijkl}(\vec{P}_{\perp})$ term in the hard factor and the
incoming linearly polarized photon distributions. Two anisotropies proposed in
Refs Li:2019yzy ; Xiao:2020ddm , namely $\cos(4\phi_{P}-4\phi_{q})$ and
$\cos(4\phi_{P}-4\phi_{b})$, originate from this mechanism and belong to this
category. In addition, we also find that the same mechanism can give rise to
several novel angular correlation of three final observables $P_{\perp}$,
$q_{\perp}$, and $b_{\perp}$, such as $\cos(4\phi_{P}-2\phi_{q}-2\phi_{b})$,
$\cos(4\phi_{P}-6\phi_{q}+2\phi_{b})$, and
$\cos(4\phi_{P}-6\phi_{b}+2\phi_{q})$, and etc.

To extract the harmonic coefficients, we can use the following definition and
compute the average value of the harmonics $\cos(...)$

$2\langle\cos(...)\rangle=\frac{2\int
d\mathcal{P.S.}\cos(...)\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}}{\int
d\mathcal{P.S.}\frac{d\sigma}{dy_{1}dy_{2}d^{2}P_{\perp}d^{2}q_{\perp}d^{2}b_{\perp}}},$
(37)

where $\int d\mathcal{P.S.}$ represents the phase space integrals. In
computing the angular correlations with given $q_{\perp}$ and $b_{\perp}$, we
define it as $\int dy_{1}dy_{2}d^{2}P_{\perp}d\phi_{q}d\phi_{b}$. We show
predictions for the anisotropies as a function of $b_{\perp}$ with fixed
$q_{\perp}=50$ MeV at the LHC in the Fig. 4. As shown in these plots, we find
that there are nine azimuthal angular correlations, which are numerically
significant. As discussed above, the first category includes two angular
correlation terms $2\langle\cos(2\phi_{q}-2\phi_{b})\rangle$ and
$2\langle\cos(4\phi_{q}-4\phi_{b})\rangle$.
$2\langle\cos(2\phi_{q}-2\phi_{b})\rangle$ has a sizeable positive
contribution. As the impact parameter $b_{\perp}$ increases, it increases in
the region $b_{\perp}<2R_{A}$, and peaks at $2R_{A}$, then decreases in the
region $b>2R_{A}$. The $2\langle\cos(4\phi_{q}-4\phi_{b})\rangle$ term is
usually negative, and its magnitude increases as the impact parameter
$b_{\perp}$ increases.

The second category of the correlations are related to the $\cos(4\phi_{P})$
harmonics, and it can be further subdivided into three different types. The
first type is the angular correlations between $P_{\perp}$ and $q_{\perp}$. It
includes one angular correlation term, i.e.,
$2\langle\cos(4\phi_{P}-4\phi_{q})\rangle$. The second type is the angular
correlation between $P_{\perp}$ and $b_{\perp}$, namely,
$2\langle\cos(4\phi_{P}-4\phi_{b})\rangle$. There are lots of related
discussions found in Refs. Li:2019yzy ; Li:2019sin ; Xiao:2020ddm ;
Shao:2022stc ; Shao:2023zge  about these two types of correlations. The last
type is the angular correlation that mixes all three vectors, $q_{\perp}$,
$b_{\perp}$, and $P_{\perp}$. Here, we show five angular correlations as
examples. Among them, the $2\langle\cos(4\phi_{P}-2\phi_{q}-2\phi_{b})\rangle$
term has the largest magnitude. This correlation is negative and its magnitude
increases as the impact parameter $b_{\perp}$ increases, and the value tends
to be constant in the region $b>2R_{A}$.
$2\langle\cos(4\phi_{P}+2\phi_{b}-6\phi_{q})\rangle$ and
$2\langle\cos(4\phi_{P}+2\phi_{q}-6\phi_{b})\rangle$ are positive, while
$2\langle\cos(4\phi_{P}+4\phi_{b}-8\phi_{q})\rangle$ and
$2\langle\cos(4\phi_{P}+4\phi_{q}-8\phi_{b})\rangle$ are negative. The
amplitudes of these four terms all show similar behaviour, and they increases
as the impact parameter $b_{\perp}$ increases. Similarly, as shown in the Fig.
5, we can find that these anisotropies show similar magnitude as the function
of $q_{\perp}$ at fixed $b=10$ fm in the LHC.

At a distance of $b=10$ fm, the colliding nuclei undergo peripheral
collisions, which also generates a large number of hadrons. Consequently, the
reaction plane of the hadronic events can be used to determine the impact
parameter vector $\vec{b}_{\perp}$ in peripheral collisions. By combining the
hadronic reaction plane information with the dilepton measurements, it becomes
possible to measure various dilepton angular correlations between $b_{\perp}$
and $q_{\perp}$, $P_{\perp}$ in peripheral collisions.

In addition, we show predictions about anisotropies as a function of
$b_{\perp}$ and $q_{\perp}$ at RHIC energy in Fig. 6 and Fig. 7. In these
predictions, we determine the values of the pair mass $M$, the relative
transverse momentum of the lepton pair $P_{\perp}$, and the rapidity region
based on the relevant kinematic region on RHIC and LHC. Although the
kinematics are different at the RHIC and LHC, the longitudinal momentum
fractions at different energies are almost identical, and the photon Wigner
distributions at different energies are comparable. Therefore, we can see that
the predictions show very similar behaviors at the RHIC and LHC. In principle,
there are infinite number of correlations in this series. We only plot the
nine examples of the angular correlations in the above plots. In the future,
it would be interesting to compare our results with RHIC and the LHC
measurements and it will help us to understand the property of the photon
Wigner distribution. We note that the mass effect has been neglected. If one
consider the mass effect, more correlation can be observed. We leave it for
future study.

####  .2.2 The Sudakov Effect

Figure 8: Comparisons our results with the ATLAS experimental data
ATLAS:2016vdy  in the dimuon pair acoplanarity $\alpha$ distribution in UPCs.

This subsection will discuss the effect of multiple soft-photon radiation,
which is also known as the Sudakov soft radiation effect. These soft
radiations give rise to large logarithms that must be resummed. The
resummation improved cross-section is given as Klein:2018fmp ; Klein:2020jom ;
Hatta:2021jcd ; Shao:2022stc ; Shao:2023zge

$\frac{d\sigma}{\pi db_{\perp}^{2}d^{2}q_{\perp}d^{2}P_{\perp}}=\int
d^{2}l_{\perp}\mathcal{S}(l_{\perp},Q,m_{l})\frac{d\sigma}{\pi
db_{\perp}^{2}d^{2}q^{\prime}_{\perp}d^{2}P_{\perp}}(q^{\prime}_{\perp}=q_{\perp}-l_{\perp}).$
(38)

As shown in Refs. Klein:2020jom ; Shao:2023zge , the Sudakov factor at the
leading double logarithm order in the center-of-mass frame of the final lepton
pair reads

$\mathcal{S}(r_{\perp},Q,m_{l})=\frac{\alpha_{\rm
em}}{2\pi}\left[\ln^{2}\frac{Q^{2}r^{2}}{c_{0}^{2}}-\ln^{2}\frac{m_{l}^{2}r^{2}}{c_{0}^{2}}\theta\left(m^{2}_{l}-r^{2}/c_{0}^{2}\right)\right].$
(39)

with the lepton mass $m_{l}$ and $\alpha_{\rm em}=1/137$. As shown in Ref.
Klein:2020jom , the resummation has an analytic form in the momentum space,
which can be expressed as follows

$\displaystyle\mathcal{S}(l_{\perp},Q,m_{l})$ $\displaystyle=$
$\displaystyle\int\frac{d^{2}r_{\perp}}{(2\pi)^{2}}e^{il_{\perp}\cdot
r_{\perp}}e^{-\frac{Q_{0}^{2}r_{\perp}^{2}}{4}}e^{-S_{\alpha}(r_{\perp},Q,m_{l})}$
(40) $\displaystyle=$ $\displaystyle\frac{e^{\gamma_{0}}}{\pi
Q_{0}^{2}}\left(\frac{c_{0}^{2}Q_{0}^{2}}{4Q^{2}}\right)^{\beta_{1}}\Gamma(1-\beta_{1})~{}_{1}F_{1}(1-\beta_{1},1,-\frac{l_{\perp}^{2}}{Q_{0}^{2}}),$
(41)

with $\beta_{1}=\frac{\alpha_{\rm
em}}{\pi}\ln\frac{Q^{2}}{l_{\perp}^{2}+m_{l}^{2}}$,
$\gamma_{0}=\frac{\alpha_{\rm em}}{2\pi}\ln^{2}\frac{Q^{2}}{m_{l}^{2}}$,
$c_{0}=2e^{-\gamma_{E}}$, the Euler constant $\gamma_{E}$, and Hypergeometric
function $~{}_{1}F_{1}$. To suppress the large $r_{\perp}$ contribution and
expedite the numerical calculations, we use a Gaussian factor with a Gaussian
width $Q_{0}=10$ MeV. As a consistent check, we find that our numerical
results is insensitive to the choice of the width as long as $Q_{0}$ is kept
small.

In Fig. 8, we present our calculations with and without the Sudakov effect
compared to the ATLAS data ATLAS:2016vdy . The Leading Order (LO) calculation
could provide a good description of the lepton pair acoplanarity in the small
$\alpha$ region and small-$q_{\perp}$ broadening. However, the cross-sections
without the Sudakov factor fall too rapidly at large acoplanarity regions.
When considering leptons with larger acoplanarity or larger $q_{\perp}$, one
must include the multiple emissions of soft photons. The calculations that
include the Sudakov effect better fit the experimental data at large
acoplanarity regions. It demonstrates that Sudakov resummation is critical for
accurately computing large acoplanarity regions in dilepton production.

####  .2.3 Dilepton Events with Neutrons Emission

In this section, we focus on the measurements of dilepton production with
neutron emissions and explore the resulting bias in the impact parameter and
its corresponding connection with the incoherent photon contribution.
Recently, there have been many measurements on dilepton production with
neutron tagging from the STAR, CMS, and ATLAS Collaborations in UPCs
STAR:2019wlg ; ATLAS:2020epq ; CMS:2020skx ; ATLAS:2022srr . When an event
involves neutron emission, these neutrons, which are almost collinear with the
beam direction, can be detected by the Zero Degree Calorimeter (ZDC).

Generally, the number of emitted neutrons depends on the impact parameter
$b_{\perp}$ between two heavy ions. As assumed in Ref. Baltz:2002pp , the
probability of having $N$ neutron emissions follows the Poisson distribution

$P(N,b_{\perp})=\frac{1}{N!}[P_{s}(b_{\perp})]^{N}\exp[-P_{s}(b_{\perp})],$
(42)

where $P_{s}(b_{\perp})=S/b_{\perp}^{2}$ is the lowest-order probability for
an excited nucleus to emit one neutron, and $P(N,b_{\perp})$ represents the
probability of emitting $N$ neutrons. Here, one uses the parameter $S$ to
indicate the strength of neutron emissions. As shown in Ref. Baur:1998ay , the
Giant-Dipole Resonance model gives an expression for $S$ as follows

$S=5.45\times 10^{-5}\frac{Z^{3}(A-Z)}{A^{2/3}}\text{ fm}^{2},$ (43)

with $A=$ atomic mass number and $Z$ being the proton number. For Pb-Pb
collisions, the above formula gives $S=108\text{ fm}^{2}$. In addition, one
may also choose another parametrization $S=303\text{ fm}^{2}$ which includes
all additional radiation mechanism for Pb-Pb collisions Baltz:1997di . These
two different prescriptions show different dependency relations between the
probability of neutron emission and the impact parameter. They will be
referred to as PAR I and PAR II in the subsequent numerical studies,
respectively. We find that some observables are sensitive to the integration
range of the impact parameter, and we will discuss in detail the effects of
these two prescriptions in the next section. Based on the Poisson
distribution, we can write the probability of different numbers of neutron
emission as follows

$\displaystyle P_{0}(b)$ $\displaystyle=$ $\displaystyle\exp[-P_{s}(b)],$ (44)
$\displaystyle P_{1}(b)$ $\displaystyle=$ $\displaystyle
P_{s}(b)\exp[-P_{s}(b)],$ (45) $\displaystyle P_{M}(b)$ $\displaystyle=$
$\displaystyle 1-\exp[-P_{s}(b)],$ (46) $\displaystyle P_{X}(b)$
$\displaystyle=$ $\displaystyle 1-\exp[-P_{s}(b)]-P_{s}(b)\exp[-P_{s}(b)],$
(47)

where $P_{0}(b)$ is the probability of no neutron emission, $P_{1}(b)$ is the
probability of one neutron emission, $P_{M}(b)$ is the probability of at least
one neutron emission, and $P_{X}(b)$ is the probability of at least two
neutron emissions.

Furthermore, we can categorize the UPC event into different scenarios
according to the number of collinearly emitted neutrons detected in the ZDC on
both ends. For example, the scenario $1nXn$ means that exactly one neutron is
detected on one ZDC, and at least two neutrons are detected in the other ZDC.
$P_{nn}(b)$ is then defined as the combined probability of various type of
neutron emissions from both nuclei as a function of impact parameter. We
summarize various scenarios in the following table and plot the probability
distributions in the right part of Fig. 9. As shown in the plot of the
combined probabilities, the neutron emission multiplicity biases toward a
small value of the impact parameter. This means that the distance between two
colliding nuclei is, on average, smaller in the events with higher neutron
multiplicity. Therefore, neutron emissions can be used as a tool to select
events with different impact parameters.

| Neutrons emission | $P_{nn}(b)$
---|---|---
1 | $0n0n$ | $P_{00}(b)=[P_{0}(b)]^{2}$
2 | $0n1n$ | $P_{01}(b)=2P_{0}(b)P_{1}(b)$
3 | $0nMn$ | $P_{0M}(b)=2P_{0}(b)P_{M}(b)$
4 | $0nXn$ | $P_{0X}(b)=2P_{0}(b)P_{X}(b)$
5 | $1n1n$ | $P_{11}(b)=[P_{1}(b)]^{2}$
6 | $1nXn$ | $P_{1X}(b)=2P_{1}(b)P_{X}(b)$
7 | $MnMn$ | $P_{MM}(b)=[P_{M}(b)]^{2}$
8 | $XnXn$ | $P_{XX}(b)=[P_{X}(b)]^{2}$

Figure 9: Left: Table of the combined probability for various neutron
emissions. Right: The plot of the combined probabilities as functions of
$b_{\perp}$.

In addition, to include the effect of neutron emission in our framework, we
cast the cross-section into

$\displaystyle\frac{d\sigma_{nn}}{d^{2}P_{\perp}d^{2}q_{\perp}}$
$\displaystyle=$ $\displaystyle\int^{\infty}_{b_{\min}}db_{\perp}b_{\perp}\int
d\phi_{b}\frac{d\sigma}{d^{2}b_{\perp}d^{2}P_{\perp}d^{2}q_{\perp}}P_{nn}(b_{\perp}),$
(48)

where $b_{\rm min}$ is defined as the minimum impact parameter value to
distinguish between peripheral collisions and the UPCs. For the case of
$0n0n$, we note that the neutron emission probability $P_{00}$ is non-
vanishing at large $b_{\perp}$, resulting in slow numerical convergence when
combined with the phase factor. Since the probability of $P_{00}(b_{\perp})$
approaches one when the impact parameter is much larger than $R_{A}$, we can
use the following technique to calculate the $0n0n$ case numerically. We first
separate the cross-section into two terms as follows

$\displaystyle\frac{d\sigma_{00}}{d^{2}P_{\perp}d^{2}q_{\perp}}$
$\displaystyle=$ $\displaystyle\int^{\infty}_{b_{\rm
min}}db_{\perp}b_{\perp}\int
d\phi_{b}\frac{d\sigma}{d^{2}b_{\perp}d^{2}P_{\perp}d^{2}q_{\perp}}P_{00}(b_{\perp})$
(49) $\displaystyle=$ $\displaystyle\int^{\infty}_{b_{\rm
min}}db_{\perp}b_{\perp}\int
d\phi_{b}\frac{d\sigma}{d^{2}b_{\perp}d^{2}P_{\perp}d^{2}q_{\perp}}-\int^{\infty}_{b_{\rm
min}}db_{\perp}b_{\perp}\int
d\phi_{b}\frac{d\sigma}{d^{2}b_{\perp}d^{2}P_{\perp}d^{2}q_{\perp}}\left[1-P_{00}(b_{\perp})\right],$

where the second term is convergent rapidly in numerical computation, and the
first term is the full UPC contribution

$\frac{d\sigma^{\rm UPC}}{d^{2}P_{\perp}d^{2}q_{\perp}}=\int^{\infty}_{b_{\rm
min}}db_{\perp}b_{\perp}\int
d\phi_{b}\frac{d\sigma}{d^{2}b_{\perp}d^{2}P_{\perp}d^{2}q_{\perp}}.$ (50)

We can analytically perform the integration over impact parameter from
$b_{\rm\min}$ to infinity for $\sigma^{\rm UPC}$ by following the same
technique as suggested in Ref. Klein:2020jom .

####  .2.4 The Incoherent Contribution to the Photon Distribution

Figure 10: The illustration of the coherent and incoherent contributions to
the GTMD photon distribution, with the coherent nucleus represented by dashed
lines, the incoherent proton contribution by dotted lines, and the incoherent
quark contribution by dash-dotted lines. The solid lines depict the combined
distribution of coherent and incoherent contributions. We set the momentum
fraction at $x=0.0016$ and $\Delta_{\perp}$ to $0$ GeV.

Let us discuss the incoherent contributions to the photon distribution and its
connection to the neutron emission events. For dilepton production from the
photon-photon fusion process, when the momentum imbalance $q_{\perp}$ of the
lepton pair is small compared to $1/R_{A}$ (where $R_{A}$ is the radius of the
nucleus), the nucleus, as a whole, can be treated as the source of photons.
This contribution is commonly referred to as the coherent contribution.
However, as the momentum imbalance of the lepton pair increases, the
transverse momenta carried by the incoming photons should also increase
accordingly. Therefore, according to the uncertainty principle, the photon
with large transverse momentum no longer originates solely from the coherent
nuclear radiation process but can be emitted from the inner structure of the
nucleus, such as the protons or the quarks inside the protons and neutrons.
These contributions are known as the incoherent contributions.

If quasi-photons with relatively large transverse momenta are radiated from
the interior of heavy nuclei, these emissions can likely lead to the
disintegration of these nuclei. Consequently, the structural integrity of the
heavy nuclei can serve as a key indicator in determining whether this
incoherent interaction occurs or not. As discussed in the above section, lots
of experimental data with the neutron emission effect have been observed by
STAR, CMS, and ATLAS Collaborations STAR:2019wlg ; ATLAS:2020epq ; CMS:2020skx
; ATLAS:2022srr . For the $0n0n$ event, the outgoing nuclei are assumed to
still remain intact, thus the dominant contribution is from the coherent
photon emissions. In the $MnMn$ and $0nMn$ events in which ZDC can detect a
number of neutrons, the outgoing nuclei are disintegrated, and incoherent
photon interactions should take place. Therefore, we should include the
incoherent contribution in our framework to provide a complete description of
the experimental data with neutron emissions.

To account for the incoherent contributions, it is still reasonable to treat
the incoming photon as a quasi-real particle and employ the generalized
equivalent photon approximation for the calculation. The only difference is
that we are treating the nucleus as an incoherent collection of electric
charges. Consequently, considering the coherent and incoherent contributions
together, one can modify the corresponding
$\Gamma^{ij}(x,k_{T},\Delta_{\perp})$, and write it as

$\Gamma^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=\Gamma_{A}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})+\Gamma_{p}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})+\Gamma_{q}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp}),$
(51)

where the first term is the coherent contribution from the nucleus. The second
and last terms are the incoherent contributions from the protons and quarks
inside nucleons, respectively. When the protons inside the heavy nucleus are
viewed as the photon sources, the corresponding total GTMD of the nucleus is
given as the incoherent sum of the photon distributions from individual
protons

$\Gamma_{p}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=Z\frac{\alpha_{\rm
em}}{\pi^{2}}k^{i}_{\perp}{k^{\prime}}^{j}_{\perp}\frac{F_{p}(k)}{k^{2}}\frac{F_{p}(k^{\prime})}{{k^{\prime}}^{2}},$
(52)

where $F_{p}(k)$ is the proton form factor. When the photon is radiated
incoherently by the quark inside protons or neutrons, we can write the photon
GTMD by treating the quarks as point particles, which is presented as follows

$\Gamma_{q}^{ij}(x,\vec{k}_{T},\vec{\Delta}_{\perp})=\left[Z\sum
e_{q/p}^{2}+(A-Z)\sum e_{q/n}^{2}\right]\frac{\alpha_{\rm
em}}{\pi^{2}}k^{i}_{\perp}{k^{\prime}}^{j}_{\perp}\frac{1}{k^{2}}\frac{1}{{k^{\prime}}^{2}},$
(53)

where $e_{q}$ is the charge number of quark from protons and neutrons. In the
above equation, the first and second terms correspond to the contribution of
the quark inside protons and neutrons, respectively. According to the
experimental cuts used in the measurement, the transverse momentum imbalance
of the final lepton pairs $q_{\perp}$ is expected to be less than $1.2$ GeV,
which means the quantum fluctuations within nucleons are comparatively small
and the QCD evolution effect may be neglected. Therefore, we can use the
valence quark model to simplify the calculation. In this model, a proton
comprises two up quarks and one down quark, while a neutron consists of one up
quark and two down quarks, with these quarks inside nucleons taking the same
momenta.

In our numerical calculations, we employ the following parametrizations of the
nuclear form factor Davies:1976zzb
$F_{A}(k)=\frac{3}{(kR_{A})^{3}}\frac{1}{1+a^{2}k^{2}}\left[\sin(kR_{A})-kR_{A}\cos(kR_{A})\right]$,
where $a=0.71$ fm, and $R=6.6$ fm for $Pb$ and $R=6.4$ fm for $Au$ and the
proton form factor Andivahis:1994rq ; Klein:2003vd
$F_{p}(k)=\frac{1}{[(k^{2}/0.71\text{ GeV}^{2})+1]^{2}}$.

For comparison, we plot three different contributions (coherent, incoherent
proton part, and incoherent quark part) to the GTMD photon distribution in
Fig. 10. The momentum fraction $x$ in this plot is estimated with
$p_{1\perp}=p_{2\perp}=4$ GeV and $\Delta_{\perp}=0\text{ GeV}$ at midrapidity
$y_{1}=y_{2}=0$ at LHC energy. We observe that the photon distribution
originating from the coherent nucleus $\Gamma_{A}$ dominates at low $k_{T}$
values, while the incoherent contributions are a couple of magnitudes smaller
than the coherent nucleus contribution. In the large $k_{T}$ regime, the
coherent contribution $\Gamma_{A}$ falls off rather rapidly as $k_{T}$
increases, and the incoherent distribution becomes much larger than the
coherent contribution when $k_{T}>0.15$ GeV. Furthermore, the quark
contribution to the photon distribution $\Gamma_{q}$ is larger than that from
protons $\Gamma_{p}$. For the sum of coherent and incoherent contributions, in
the small $k_{\perp}$ regions, the nucleus acts as the dominant photon source,
which means the modified GTMD is consistent with the traditional GTMD in the
small $k_{\perp}$ region. In the large $k_{\perp}$ region, the incoherent
distribution from quarks takes over. It is interesting to note that the
combined photon GTMD distribution is proportional to $Z^{2}$ in the small
transverse momentum region and proportional to $(A-Z)\sum e_{q/n}^{2}+Z\sum
e_{q/p}^{2}=\frac{2}{3}(A-Z)+Z$ in the asymptotically large transverse
momentum region. The incoherent emission of photons is rare compared to the
coherent emission, yet it is more efficient in generating high $k_{T}$ photons
than the coherent emission.

For dilepton production, if the transverse momentum imbalance of the final
lepton pairs $q_{\perp}$ is small, the two colliding initial photons typically
possess low transverse momenta, which in turn suggests that the photons
originate predominantly from the coherent nucleus. Under such conditions, the
photons “see” the nucleus as a whole. Conversely, when $q_{\perp}$ is large,
either or both of the initial photons are likely to have high transverse
momenta and are typically emitted by protons or quarks. This means that these
high-momentum photons can allow us to perceive the internal structures within
the nucleus. Therefore, this incoherent process is likely to excite the
nucleus, causing it to emit one or more neutrons in the beam direction, which
can subsequently be recorded by the ZDC detector. The incoherent contribution
is expected to be the dominant production mechanism for $q_{\perp}$ broadening
in dilepton production for events with one or more neutron emissions, such as
$0nMn$ and $MnMn$. In the next section, we will compare our numerical
calculations with experimental data involving neutron tagging and show that
the inclusion of the incoherent contribution is required if one wishes to
describe both the data with and without neutron emissions.

###  .3 Comparisons with the RHIC and LHC data

In this section, we compare our numerical results with the dilepton
experimental data collected in a wide variety of measurements. Firstly, we
provide a numerical comparison with the acoplanarity and $q_{\perp}$
broadening results. Furthermore, we compare our results with the neutron-
tagging data from STAR, CMS, and ATLAS Collaborations. Also, we turn on the
incoherent contributions to account for enhancements in the large $\alpha$
region when neutron emissions are detected. To minimize various uncertainty,
we choose to normalize both the computed cross-section and the experimental
data to unity. In addition, as the measurements conducted thus far have
focused solely on the angular correlation at a fixed impact parameter
$b_{\perp}$ without specifying its orientation, we integrate over the
azimuthal angle of $b_{\perp}$ in the subsequent numerical calculations when
comparing with the data.

####  .3.1 Acoplanarity and Transverse Momentum Broadening

Figure 11: Comparisons of our results with the ATLAS experimental data
ATLAS:2022yad  in the dimuon pair acoplanarity $\alpha$ distribution from
central collisions to UPC. Figure 12: Comparisons of our results with the
ATLAS experimental data ATLAS:2022yad  in the dimuon pair $k_{\perp}$
distribution from central collisions to UPC. Figure 13: Comparisons of our
results with the STAR experimental data STAR:2019wlg  in the dilepton pair
transverse momentum imbalance $q_{\perp}$ distribution at $60-80\%$ central
collisions.

In this subsection, we provide quantitative analyses on the $\alpha$
distribution and its related observable $k_{\perp}$ in comparison with the
ATLAS experimental data. We then further calculate the $q_{\perp}$
distribution and compare it with the STAR measurements for studying the
dilepton broadening effects.

The acoplanarity $\alpha$ is the physical observable describing the azimuthal
angular correlation between the lepton pairs. It is defined as

$\alpha=1-|\Delta\phi|/\pi,$ (54)

where $\Delta\phi=\phi_{1}-\phi_{2}$ is the angle difference between the
azimuthal angle of lepton and anti-lepton $\phi_{1}$ and $\phi_{2}$. In Fig.
11, we compare our results with the ATLAS experimental data ATLAS:2022yad  in
the dimuon pair acoplanarity $\alpha$ from central collisions to UPC. The
values of the impact parameter corresponding to the different centralities are
extracted using the Glauber model. Since the normalized central collision
results are almost identical for the impact parameters $b_{\perp}=0$, $1.75$,
$3.5$ fm, thus we simply use the average $b_{\perp}=1.75$ fm for the $0-5\%$
centrality bin. One observes that our results can perfectly describe the
experimental data. Although the dominant contribution is attributed to the
angle average of $\phi_{P}$ term $O^{ijkl}$, it does not capture the dip
observed in very small $\alpha$ regions at central collisions. By including
the angular correlation $\cos 4\phi_{P}$ term $\Omega^{ijkl}$, we can describe
the dip across various centrality ranges. Hence, the role of the
$\Omega^{ijkl}$ term is crucial in accurately describing the dip in central
collisions. As the impact parameter $b_{\perp}$ increases, the significance of
the $\Omega^{ijkl}$ term diminishes. Notably, in the case of UPC with large
impact parameters, the dominant contribution stems from the $O^{ijkl}$ term,
rendering the $\Omega^{ijkl}$ term negligible. In addition, we can find that
the peak of $\alpha$ distribution is shifted to smaller values of $\alpha$
from central collisions to UPC. We note that these numerical findings can also
be derived from the computations in Refs. Zha:2018tlq ; Li:2019yzy ;
Li:2019sin ; Wang:2021kxm ; Wang:2022gkd ; Shao:2023zge , with this dip
arising from the presence of the $\cos 4\Delta\phi$ term.

Besides the $\alpha$ distribution, ATLAS experiments have also used another
physical observable $k_{\perp}$ which is closely related to $\alpha$. This
observable can provide us a good estimate of the momentum imbalance and it is
defined as

$k_{\perp}=\pi\alpha P_{\perp},$ (55)

where $P_{\perp}$ is the relative transverse momentum of the lepton pair. As
shown in Fig. 12, we compare our results with the ATLAS experimental data
ATLAS:2022yad  in the dimuon pair acoplanarity $k_{\perp}$, ranging from
central collisions to UPC. Our results offer a good description of the data
across a wide range. We observe that the $k_{\perp}$ distribution shows
similar behaviours with the $\alpha$ distribution, where the dip at small
$k_{\perp}$ for small $b_{\perp}$ can be explained by including the
contribution of the anisotropic $\Omega^{ijkl}$ term.

Additionally, we present our calculation alongside the STAR experimental data
STAR:2019wlg , in the transverse momentum imbalance $q_{\perp}$ distribution
of dilepton pairs at $60-80\%$ central collisions. Our results provide a
satisfactory description of the data. Given that the momentum imbalance
$q_{\perp}$ is an angle-independent observation, only the average term
$O^{ijkl}$ of the hard factor contributes, and the $\Omega^{ijkl}$ term of the
hard factor vanishes.

Therefore, the above figures show that the dilepton acoplanarity and the
broadening of the momentum imbalance $q_{\perp}$ encode the initial photons’
transverse momentum information and they allow us to further explore the
$5$-dimensional photon Wigner distribution function and GTMD.

####  .3.2 The Neutron Emission effect in the Small $\alpha$ Region

Figure 14: Left: Comparison of our numerical results for the transverse
momentum imbalance $q_{\perp}$ distribution of dielectron pairs, with and
without triggering the neutron emission, with the STAR experimental data
STAR:2019wlg  measured in UPCs with neutron emissions; Right: Comparisons of
our results for the dimuon pair $\langle\alpha\rangle$ with the CMS
experimental data CMS:2020skx . Figure 15: Comparisons of our results with the
CMS experimental data CMS:2020skx  for the dimuon pair acoplanarity $\alpha$
distribution with the neutron emission in UPCs.

In the small $\alpha$ region or low transverse momentum imbalance regime, the
neutron multiplicity in UPCs can serve as an indicator of the impact
parameter, as more neutron emission biases toward smaller impact parameter.
Firstly, we show the impact of neutron emission on the $q_{\perp}$
distribution at RHIC energy. Secondly, we study the average acoplanarity
$\langle\alpha\rangle$ distribution measured at CMS Collaboration. The
comparisons of our numerical calculations with the acoplanarity $\alpha$
distribution from CMS Collaboration are also presented for different
configurations of neutron emission tags. Lastly, we compare our numerical
results with the $\alpha$ distribution measured at ATLAS Collaboration in the
large $\alpha$ region.

In the left plot of Fig. 14, we compare our numerical results with the
neutron-tagged UPC data from the STAR Collaboration STAR:2019wlg . The STAR
experiment uses $Mn$ to denote the emission of one or more neutrons. Here, we
present two numerical calculations, one including the corresponding
probability to account for the neutron emission effect, and the other without.
For this calculation, we adopt the prescription given in Eq. (43) (PAR I) and
choose $b_{\rm min}=2R_{A}$. The other parametrization gives similar results.
Since the neutron emission probability function for $MnMn$ is biased towards
smaller impact parameter $b_{\perp}$ values in UPC, we can see from the result
in the previous subsection that $q_{\perp}$ broadening effect is stronger in
smaller values of $b_{\perp}$. Therefore, we observe that including the
neutron emission probability shifts the peak of the $q_{\perp}$ distribution
to a higher $q_{\perp}$ value. We can find that the calculation combined with
the neutron emission probability offers an accurate description of the
neutron-tagged UPC experimental data from the STAR Collaboration. This
calculation shows that the neutron emission probability is important for
describing the experimental data with neutron tagging in the low $q_{\perp}$
regime.

In the right plot of Fig. 14, we present the comparison of the average
normalized acoplanarity, $\langle\alpha\rangle$, with the CMS data reported in
six neutron multiplicity classes. For this calculation, we utilize both
parametrizations for neutron emission probabilities: PAR I with $S=108\text{
fm}^{2}$ in Eq. (43) for lead ions, and PAR II with $S=303\text{ fm}^{2}$.
These two prescriptions yield different probabilities of neutron emission
scenarios as a function of the impact parameter $b_{\perp}$. As shown in the
figure, these two prescriptions show the right trend and magnitudes comparing
to the CMS data. Furthermore, the average acoplanarity value increases as
neutron scenarios transition from $0n0n$ to $XnXn$. We observe that the
$\langle\alpha\rangle$ depicts a strong dependence on the parameter $S$ in
different neutron scenarios, and the average value $\langle\alpha\rangle$
decreases as the parameter $S$ increases. For example, for the $XnXn$, the
peak of the neutron emission probability $P_{nn}(b_{\perp})$ shifts towards a
smaller impact parameter $b_{\perp}$ as the parameter $S$ decreases.
Furthermore, the edges of the two bands are calculated by varying $b_{\min}$
from $2R_{A}$ to $14$ fm in Eq. (48), and we notice that the average value
$\langle\alpha\rangle$ decreases as $b_{\min}$ increases. The average values
$\langle\alpha\rangle$ are sensitive to the lower limit of integration
$b_{\min}$, especially for high neutron multiplicity classes.

In Fig. 15, we present the comparison of the normalized $\alpha$ distribution
with the CMS data for six neutron multiplicity classes. In this calculation,
we employ two different parametrizations for the neutron emission probability:
PAR I and PAR II, as previously defined. Both prescriptions can provide a good
description of the CMS data. Interestingly, the peak value of the $\alpha$
distribution shifts to larger value from the $0n0n$ class to the $XnXn$ class.
In the large $\alpha$ region, we find that these distributions are insensitive
to the choice of the parameter $S$. In summary, our analysis reveals that the
neutron emission probability is essential in describing the measurement with
the neutron tagging.

####  .3.3 Enhancements from the Incoherent Contributions in the Large
$\alpha$ Region

Figure 16: Comparisons of our result with the ATLAS experimental data
ATLAS:2020epq  in the dimuon pair acoplanarity $\alpha$ distribution for
$0n0n$ with the rapidities regions of the dimuon $|y^{\mu\mu}|<2.4$. Figure
17: Comparison of our results with the CMS experimental data CMS:2020skx  in
the dimuon pair acoplanarity $\alpha$ distribution with neutron emissions in
UPCs. Figure 18: Comparison of our results with the ATLAS experimental data
ATLAS:2020epq  in the dimuon pair acoplanarity $\alpha$ distribution, taking
into account neutron emission in UPCs. Figure 19: Comparing of our results
with the ATLAS experimental measurements ATLAS:2022srr  of the $e^{+}e^{-}$
pair acoplanarity $\alpha$ distribution with neutron emissions in UPCs.

When high transverse momentum photons (i.e., incoherent photons) are emitted
from the substructure of a nucleon (proton, neutron), it often leads to the
breakup of the emitting nucleus. Using the study on events with neutron
emissions in the small $\alpha$ region discussed in the last subsection as
calibration, we further incorporate the incoherent contributions into our
numerical study and explore the production mechanism of the $q_{\perp}$
broadening effect on the data involving the neutron emission effect in the
large $\alpha$ region. As a result, we can demonstrate that the incoherent
contributions are indispensable for fully describing the $\alpha$ distribution
measurements with neutron emissions.

In the following calculations, since the rest of the observables are not very
sensitive to the choice of neutron emission prescription, we will only show
numerical results based on Eq. (43), PAR I, as the prescription to account for
the neutron emission effect.

Firstly, for the events with one or no neutron emission, the outgoing nucleus
is almost intact, and the radiated photons are mostly from the nucleus as a
whole. Therefore, the radiation of coherent photons should be the dominant
contribution in these events. As shown in the above section, the Sudakov
resummation, which resums the final-state soft-photon radiations, is the
dominant mechanism for the $q_{\perp}$ broadening in the large $\alpha$
region. In Fig. 16, we show the comparison of our calculations involving the
Sudakov effect with the ATLAS $0n0n$ data measured in the rapidities regions
$|y^{\mu\mu}|<2.4$ ATLAS:2020epq . We find that the Sudakov calculation can
describe the $0n0n$ data from the ATLAS Collaboration in the large $\alpha$
region. In upper panel of Fig. 17, we compare our results incorporating the
Sudakov effect with the CMS experimental data CMS:2020skx . For the three
scenarios ($0n0n$, $0n1n$, and $1n1n$), we observe that coherent results
incorporating the Sudakov effect accurately describe the CMS data in the large
$\alpha$ region. Similarly, as indicated in top panel of Fig. 18 and Fig. 19,
the coherent contribution with the Sudakov effect accurately describes both
the $0n0n$ and UPC data from ATLAS in dimuon ATLAS:2020epq  and dielectron
ATLAS:2022srr  pair production process. These results serve as the benchmark
for our calculation before including the incoherent contributions.

In addition, the neutron emissions indicate that the disassociation of the
nucleus when two colliding nuclei get too close to each other in UPCs. In this
scenario, the photon might not only be radiated from the nucleus as a whole
but also from the protons and quarks inside nucleons. Therefore, the
incoherent contribution naturally becomes important for events with more than
one neutron emission. Thus, we turn on the incoherent contribution to the
photon Wigner distribution. In the lower panel of Fig. 17, we present a
comparison between our results with and without the incoherent effect and the
CMS data across three scenarios ($0nXn$, $1nXn$ and $XnXn$) in dimuon pair
production. Similarly, the low panel of Figs. 18 and 19 illustrate comparisons
of our results with ATLAS data in two scenarios ($0nMn$ and $MnMn$) for
dielectron and dimuon pair production. The magenta dashed line represents only
the coherent contribution with the Sudakov effect, and the black solid line
represents the sum of coherent and incoherent contributions with the Sudakov
effect. In the small-$\alpha$ region, both curves agree with each other. It
shows that, for the event with neutron emissions, the coherent contribution,
which stands for the radiated photon from the nucleus as a whole, is still the
dominant production mechanism of the $q_{\perp}$ broadening effect in the
small $\alpha$ region. In the large $q_{\perp}$ region, these coherent
calculations with the Sudakov effect underestimate the experimental data for
tags with neutron emissions, but the sum of the coherent and incoherent
contribution with the Sudakov effect provides good descriptions of the data
from CMS and ATLAS Collaborations. It means that besides the Sudakov effect,
the incoherent contribution, which represents the radiated photon from protons
and quarks, become the dominant contribution to the $q_{\perp}$ broadening
effect in the large $\alpha$ regions for events with multiple neutron
emissions.

Let us interpret the numerical results shown in this part. It is evident that
the transverse momentum imbalance $q_{\perp}$ of the lepton pair is small in
the small-$\alpha$ region. Thus, the initial heavy nucleus as a unified entity
radiates photons with low momenta, and final dileptons do not know the inner
structure inside heavy nuclei. In this region, the dominant contribution in
the small-$\alpha$ region is from coherent photons. On the contrary, in the
large $\alpha$ region where the momentum imbalance of the dilepton pair
$q_{\perp}$ is much larger than $1/R_{A}$, the coherent contribution to the
broadening drops off very rapidly as $q_{\perp}$ increases. In this case,
incoming photons with high transverse momentum are mostly emitted from the
sub-structure in heavy nuclei. Consequently, incoherent photons originating
from quarks and protons start to play a vital role.

In the end, through the comprehensive and quantitative comparison with the
experimental data above, we demonstrate that the incoherent photon
contribution is indispensable for describing the data both with and without
neutrons emissions simultaneously.