# Dimuon and ditau production in photon-photon collisions at next-to-leading
order in QED

Hua-Sheng Shao<EMAIL_ADDRESS>Laboratoire de Physique
Théorique et Hautes Energies (LPTHE), UMR 7589,
Sorbonne Université et CNRS, 4 place Jussieu, 75252 Paris Cedex 05, France
David d’Enterria<EMAIL_ADDRESS>CERN, EP Department, CH-1211
Geneva, Switzerland

###### Abstract

Next-to-leading-order (NLO) quantum electrodynamics (QED) corrections to the
production of muon and tau pairs in photon-photon collisions,
$\gamma\gamma\to\mu^{+}\mu^{-},\tau^{+}\tau^{-}$, are calculated in the
equivalent photon approximation. We mostly consider $\gamma\gamma$ processes
in ultraperipheral collisions of hadrons at the LHC, but the
$\gamma\gamma\to\tau^{+}\tau^{-}$ process in $\mathrm{e}^{+}\mathrm{e}^{-}$
collisions at LEP is also discussed. The NLO terms are found to modify the
total cross sections by up to 5%, increasing the tails of the dilepton
acoplanarity and transverse momentum distributions, and depleting by up to 15%
the yields at high masses, with respect to the leading-order predictions
including the very small virtuality of the colliding photons. At the LHC, the
calculations obtained with the charge form factor for protons and lead ions
including the NLO QED corrections improve the data–theory agreement for all
measured differential distributions, and prove an indispensable ingredient for
the extraction of precision quantities in photon-photon processes, such as the
anomalous magnetic moment of the tau lepton.

##  I Introduction

The study of photon-photon ($\gamma\gamma$) interactions at multi-GeV energies
was first realized in the laboratory in $\mathrm{e}^{+}\mathrm{e}^{-}$
collisions at DESY PETRA in the 1980s Morgan:1994ip  and at CERN LEP in the
1990s Whalley:2001mk ; Przybycien:2008zz , and has received a significant
experimental and theoretical boost in collisions with hadron beams in the last
twenty years thanks to the large $\gamma$ energies and luminosities accessible
at the BNL Relativistic Heavy-Ion Collider (RHIC) Bertulani:2005ru  and CERN
Large Hadron Collider (LHC) Baltz:2007kq ; deFavereaudeJeneret:2009db . In the
equivalent photon approximation (EPA) vonWeizsacker:1934nji ; Williams:1934ad
, the electric field created by a charge accelerated at high energies can be
interpreted as a flux of quasireal photons whose energies $E_{\gamma}$ and
number densities $N_{\gamma}$ grow proportionally to the Lorentz relativistic
factor ($E_{\gamma}\propto\gamma_{L}$) and squared charge ($N_{\gamma}\propto
Z^{2}$) of the beam particles, respectively Brodsky:1971ud ; Budnev:1975poe .
At the LHC, one can exploit the availability of multi-TeV beams (leading to
boosts of $\gamma_{L}=E_{\text{beam}}/m_{\mathrm{p,N}}\approx 7500,\,3000$ in
collisions with protons and ions111For proton and nucleon masses:
$m_{\mathrm{p,N}}=0.9315,0.9383$ GeV, respectively. Natural units,
$\hbar=c=1$, are used throughout the paper. at nucleon-nucleon center-of-mass
energies (c.m.) of $\sqrt{s_{{}_{\text{NN}}}}=14,\,5.5$ TeV, respectively) and
large electric charges (up to $Z=82$ for Pb ions) to produce a very large
number of different processes via photon-photon fusion Shao:2022cly .

Most of the studies of $\gamma\gamma$ processes have been carried out in
proton-proton (p-p), proton-nucleus (p-A), and nucleus-nucleus (A-A)
“ultraperipheral” collisions (UPCs) at impact parameters larger than twice
their transverse radii, i.e., without overlap of their transverse matter
profiles, and with the protons and/or ions coherently emitting quasireal
photons. Since the hadrons remain intact after their electromagnetic
interaction, the UPCs feature extremely clean topologies characterized by just
an exclusive final state produced in an otherwise empty detector. The
requirement of coherent emission from the hadron charge distribution imposes
very low maximum virtualities of the emitted photon:
$Q_{\text{max}}^{2}\approx 1/R^{2}$, where $R$ is the charge radius, i.e.,
$Q^{2}\lesssim~{}0.08$ GeV2 for protons ($R_{\mathrm{p}}\approx 0.7$ fm), and
$Q^{2}\lesssim 8\cdot 10^{-4}$ GeV2 for Pb nuclei ($R_{\mathrm{Pb}}\approx 7$
fm). More recently, $\gamma\gamma$ collisions have also been observed in
nonexclusive final states where the hadrons overlap and produce many
additional particles in parton-parton interactions concomitant with the two-
photon interaction STAR:2018ldd ; ATLAS:2018pfw ; ALICE:2022hvk , as well as
in high-luminosity p-p runs with large number of overlapping collisions
(pileup) in the same bunch crossing CMS:2024skm . To identify
$\gamma\gamma\to\mathrm{X}$ collisions in both exclusive and nonexclusive
events, one exploits the fact that the colliding photons have very small
virtualities and, therefore, that the final state is basically produced at
rest in the transverse plane, i.e., $p_{\mathrm{T}}^{X}\approx 0$, and thus
its decay products are back-to-back in azimuthal angle $\phi$, i.e., have a
very small acoplanarity $A_{\phi}=1-|\Delta\phi|/\pi\approx 0$. In the case of
p-p collisions, one can in addition detect one or both surviving proton(s) in
dedicated very forward detectors inside the LHC beamline FP420RD:2008jqg ;
CMS:2014sdw ; Tasevsky:2015xya ; CMS:2022hly , provided that the
$\gamma\gamma$ collision produces a system with large enough masses
$m_{\gamma\gamma}\gtrsim 100$ GeV Bruce:2018yzs ; CMS:2018uvs ; ATLAS:2020mve
.

Figure 1: Selection of diagrams for the photon-fusion production of a pair of
charged leptons, $\gamma\gamma\to\ell^{+}\ell^{-}$ with
$\ell^{\pm}=e^{\pm},\mu^{\pm},\tau^{\pm}$, in collisions of electrons,
protons, and/or ions: Born (leftmost), and real (second) and virtual (rest of
the diagrams) NLO corrections.

Arguably, the simplest photon-photon collision process in the standard model
(SM) is the $t$-channel production of a pair of charged leptons,
$\gamma\gamma\to\ell^{+}\ell^{-}$ with
$\ell^{\pm}=e^{\pm},\mu^{\pm},\tau^{\pm}$, whose elementary cross section is
proportional to the square of the electromagnetic coupling, $\alpha=1/137.036$
(Fig. 1, left). Among those, the Breit–Wheeler (B–W) process
$\gamma\gamma\to\mathrm{e}^{+}\mathrm{e}^{-}$ by which pure light is
transformed into matter, was the first one studied in quantum electrodynamics
(QED) Breit:1934zz . The simplicity and relatively large cross section of the
B–W process have facilitated its measurement in $\mathrm{e}^{+}\mathrm{e}^{-}$
collisions JADE:1986axr ; L3:1997acq  as well as in UPCs at fixed-target and
collider energies (by the WA93 Vane:1992ms , CERES/NA45 CERESNA45:1994cpb ,
STAR STAR:2004bzo ; STAR:2019wlg , PHENIX PHENIX:2009xtn , CDF CDF:2006apx ;
CDF:2009xwv , ALICE ALICE:2013wjo , CMS CMS:2012cve ; CMS:2018erd ;
CMS:2018uvs ; CMS:2024tfd , and ATLAS ATLAS:2015wnx ; ATLAS:2017fur ;
ATLAS:2020mve ; ATLAS:2022srr  experiments), and found in overall agreement
with the theoretical EPA predictions at leading-order (LO) accuracy in QED.
However, incorporation of theoretical contributions from final-state $\gamma$
radiation (FSR) (second diagram of Fig. 1) Klein:2020jom ; CMS:2024tfd  —as
well as properly accounting for the polarization of the incoming photons in
order to reproduce particular dielectron azimuthal distributions
Brandenburg:2022tna — are needed in order to reach accord with the latest more
precise data. A good control of the B–W process where one or both $\rm
e^{\pm}$ radiate a hard bremsstrahlung photon, is of relevance also as a
background in studies of light-by-light (LbL) scattering
$\gamma\gamma\to\gamma\gamma$ dEnterria:2013zqi , as confirmed by the
experimental measurements ATLAS:2017fur ; CMS:2018erd ; ATLAS:2019azn ;
CMS:2024tfd .

Exclusive dimuon production, $\gamma\gamma\to\mu^{+}\mu^{-}$, is even better
controlled than the B–W process for experimental and theoretical reasons.
Experimentally, high-energy muons can be often identified and reconstructed
more accurately than electrons. Theoretically, muons are 200 times more
massive than electrons, and corrections due to collinear and/or soft photon
emissions, including lepton “dressing” to account for photon FSR, are better
under control. For these two reasons, the $\gamma\gamma\to\mu^{+}\mu^{-}$
process has been often considered as a useful “standard candle” to validate
the experimental reconstruction, including its use as a “luminometer”
Khoze:2000db ; Shamov:2002yi , and the theoretical ingredients ($\gamma$
fluxes, hadron survival probabilities, radiative corrections, etc.)
Vermaseren:1982cz ; Berends:1986if  in two-photon physics. A large number of
experimental measurements of $\gamma\gamma\to\mu^{+}\mu^{-}$ exist in
$\mathrm{e}^{+}\mathrm{e}^{-}$ PLUTO:1984gmq ; TOPAZ:1992bfv ; OPAL:1993law ;
L3:1997acq ; L3:2004ful , e-p H1:2003nlr , p-$\overline{\mathrm{p}}$
CDF:2009xey ; CDF:2009xwv , p-p CMS:2011vma ; ATLAS:2015wnx ; ATLAS:2017sfe ;
CMS:2018uvs ; ATLAS:2020mve , p-Pb ALICE:2023mfc  and Pb-Pb CMS:2020skx ;
ATLAS:2020epq  collisions, with the latest, more precise, Pb-Pb UPC results
showing increasing deviations in some differential distributions with respect
to the predictions of LO QED Monte Carlo (MC) event generators Klein:2016yzr ;
Harland-Lang:2020veo .

The third dilepton process of interest is $\gamma\gamma\to\tau^{+}\tau^{-}$,
which directly probes the photon-tau vertex and is thus sensitive to the
anomalous electromagnetic couplings of the $\tau$ lepton. Given the short
$\tau$ lifetime, which precludes measuring its properties at storage rings,
the photon-fusion channel has been long proposed to probe the tau lepton
anomalous magnetic, $a_{\tau}\equiv(g-2)_{\tau}/2$, and electric dipole,
$d_{\tau}$, moments at $\mathrm{e}^{+}\mathrm{e}^{-}$ Cornet:1995pw  and
hadron delAguila:1991rm ; Atag:2010ja ; Beresford:2019gww ; Dyndal:2020yen
colliders. Two-photon production of tau pairs has been measured in
$\mathrm{e}^{+}\mathrm{e}^{-}$ OPAL:1993law ; L3:1997acq ; DELPHI:2003nah  and
UPCs at the LHC CMS:2022arf ; ATLAS:2022ryk ; CMS:2024skm . The currently best
experimental measurement of $a_{\tau}=0.0009_{-0.0031}^{+0.0032}$ has been
achieved via $\gamma\gamma\to\tau^{+}\tau^{-}$ in p-p UPCs at 13 TeV
CMS:2024skm , superseding the previous $-0.052<a_{\tau}<0.013$ limit from LEP
DELPHI:2003nah . This latest LHC measurement has a precision and accuracy
approaching those needed to probe the first “Schwinger term” quantum
correction to $(g-2)_{\tau}$, which is common to all leptons and amounts to
$a_{\tau}^{\text{NLO}}=\alpha/(2\pi)\approx 0.00116$ Schwinger:1948iu . Future
measurements will be sensitive to higher-order SM terms, which result in
$a_{\tau}=0.00117721(5)$ Eidelman:2007sb , and to any potential beyond SM
(BSM) modification. An improved measurement of $a_{\tau}$ has become an
experimental priority given the long-standing data-theory discrepancies
observed for its muon counterpart $a_{\mu}$ Aoyama:2020ynm , and the fact
that, by standard lepton mass-scaling arguments, the $\tau$ particle is more
sensitive than the muon to potential BSM physics by a factor of
$m_{\tau}^{2}/m_{\mu}^{2}\approx 280$.

The availability of higher $\gamma\gamma\to\ell^{+}\ell^{-}$ data samples
measured at the LHC, with experimental cross section uncertainties reduced to
the few percent level, calls for an improved theoretical description to match
the experimental accuracy and precision. In this context, full next-to-
leading-order (NLO) corrections accounting for virtual and real QED emissions,
such as those shown in four rightmost diagrams of Fig. 1, are now a
requirement to extract precise physics results from photon-photon collisions.
This is particularly true when aiming at comparing data with theory to extract
precision SM parameters, such as $a_{\tau}$, and/or to search for deviations
in the absolute or differential cross section from the SM predictions due to
elusive new physics contributions. To our knowledge, the first theoretical
calculations that included NLO (QCD and/or QED) corrections for photon-photon
processes in UPCs were obtained for the $\gamma\gamma\to t\overline{t}$ and
$\gamma\gamma\to\gamma\gamma$ processes employing the gamma-UPC code in Refs.
Shao:2022cly  and AH:2023kor ; AH:2023ewe , respectively. Both studies showed
that the NLO corrections augment the predicted cross sections (by about 50%
for $t\overline{t}$ and 5% for LbL scattering) and reduce their theoretical
QCD scale uncertainties, emphasizing the need to include them for the accurate
and precise calculation of $\gamma\gamma$ processes. The purpose of this work
is to extend the number of physics processes available with improved
theoretical accuracy and precision, by calculating the NLO QED corrections for
the photon-fusion production of pairs of muon and tau leptons.

The paper is organized as follows. Section II describes the theoretical setup
used to compute the LO and NLO predictions for two-photon production of dimuon
and ditau final states. Results for exclusive
$\gamma\gamma\to\mu^{+}\mu^{-},\tau^{+}\tau^{-}$ processes in UPCs at hadron
colliders —for the ditau case, a short discussion of the photon-fusion
production in $\mathrm{e}^{+}\mathrm{e}^{-}$ collisions at LEP is also
provided— are presented in Sections III and IV, respectively, including total
cross sections as well as differential distributions. The paper is closed with
a summary in Section V. Appendix A provides technical details on the
implementation of $k_{\perp}$ smearing needed to account for the small
initial-state photon virtualities, and its approximate matching to NLO
distributions for lepton pairs.

##  II Theoretical setup

The matrix elements calculations reported in this work have been performed
within the MadGraph5_aMC@NLO framework (MG5_aMC for short) Alwall:2014hca
extended with the automated computations of NLO electroweak corrections
discussed in Ref. Frederix:2018nkq . For the UPC results with protons and
ions, the initial photon fluxes as well as the hadron survival probabilities
for exclusive final states are obtained with the gamma-UPC MC code
Shao:2022cly . The gamma-UPC program provides two types of elastic photon
fluxes as a function of impact parameter $b$, $N_{\gamma/Z}(E_{\gamma},b)$,
based on the electric-dipole (EDFF) and charge (ChFF) form factors for proton
and nuclei, respectively, as well as associated survival probabilities of the
photon-emitting hadrons. At the LHC, the $\gamma\gamma$ cross sections at a
given $\gamma\gamma$ c.m. energy computed with the ChFF flux are larger than
those obtained with the EDFF flux by about 10–20% (for p-p and p-Pb UPCs) and
20–40% (for Pb-Pb UPCs), with the difference rising for increasingly heavier
final states Shao:2022cly . In general, the gamma-UPC cross sections computed
with the ChFF flux match those derived with the Superchic MC Harland-
Lang:2020veo , whereas those computed using the EDFF photon fluxes agree with
the ones obtained with the Starlight code Klein:2016yzr . The EDFF photon
number density is divergent for small impact parameters, $b\to 0$, and an
arbitrary infrared (IR) cutoff needs to be imposed in the calculations.
Usually, the cutoff is chosen as $b>R_{\mathrm{A}}$, where $R_{\mathrm{A}}$ is
the radius parameter of the hadron transverse profile. Such a limitation
prevents from using the EDFF fluxes in nonexclusive $\gamma\gamma$ collisions
with hadron overlap. Since the ChFF photon flux is well-behaved for all
hadron-hadron impact parameters, the results obtained with it are considered
more realistic, are found to agree better with the most precise B–W results to
date CMS:2024tfd , and will be our default choice hereafter. Results obtained
with the EDFF flux will be shown below in some cases for comparison purposes
given that such a flux has been often used in the past.

The ChFF flux contains an explicit dependence on the photon transverse
momentum $k_{\perp}$ through its conjugate dependence on the impact parameter
$b$, which is related to the photon virtuality via
$Q^{2}=k_{\perp}^{2}+E_{\gamma}^{2}/\gamma_{\mathrm{L}}^{2}$. The unintegrated
photon number density at a given photon energy $E_{\gamma}$ and transverse
momentum $k_{\perp}$ (see its $b$-dependent counterpart in Eq. (13) of Ref.
Shao:2022cly ) reads

$\displaystyle N_{\gamma/Z}^{\mathrm{\tiny ChFF}}(E_{\gamma},k_{\perp})$
$\displaystyle=$ $\displaystyle 2\pi
k_{\perp}\frac{Z^{2}\alpha}{\pi^{2}}\left|\int_{0}^{\infty}\frac{\mathrm{d}bk_{\perp}^{2}}{k_{\perp}^{2}+E_{\gamma}^{2}/\gamma_{\mathrm{L}}^{2}}F_{\mathrm{ch},A}\left(\sqrt{\smash[b]{k_{\perp}}^{2}+E_{\gamma}^{2}/\gamma_{\mathrm{L}}^{2}}\right)J_{1}(bk_{\perp})\;\right|^{2}$
(1) $\displaystyle=$
$\displaystyle\frac{2Z^{2}\alpha}{\pi}\frac{k_{\perp}^{3}}{\left(k_{\perp}^{2}+E_{\gamma}^{2}/\gamma_{\mathrm{L}}^{2}\right)^{2}}\left[F_{\mathrm{ch},A}\left(\sqrt{k_{\perp}^{2}+E_{\gamma}^{2}/\gamma_{\mathrm{L}}^{2}}\right)\right]^{2}.$

where $F_{\mathrm{ch},A}(Q^{2})$ is the form factor of the ion A emitting the
photon (the full analytical expressions for protons and ions can be found in
Ref. Shao:2022cly ), and $J_{1}$ is the Bessel function of the first kind.
However, the matrix element calculations in MG5_aMC expect collinear (i.e.,
$k_{\perp}$-independent) input photon distributions, and therefore the photon
fluxes need to be fully integrated over $Q$ and $k_{\perp}$ (cf. Eq. (13) in
Ref. Shao:2022cly ). Thus, when computing a given photon-photon process at LO
accuracy with the gamma-UPC plus MG5_aMC (or HELAC-Onia222Quarkonium final
states in $\gamma\gamma$ fusion can also be simulated at LO accuracy with
gamma-UPC and the HELAC-Onia event generator Shao:2012iz ; Shao:2015vga .)
codes out of the box, the central system will be produced exactly at rest in
the transverse plane. For the $\gamma\gamma\to\ell^{+}\ell^{-}$ case, this
means $p_{\text{T}}^{\ell\ell}=0$ with both leptons emitted back-to-back in
azimuth, i.e., $A_{\phi}=0$. Since the photon density follows a
$1/k_{\perp}^{2}$ dependence and its $k_{\perp}$ values are many orders-of-
magnitude smaller than their longitudinal energy $E_{\gamma}$, the EPA
assumption that both photons are real ($Q\approx 0$) is well fulfilled, and
has no actual numerical impact on the computed fully integrated cross
sections. However, in reality the coherently emitted photons can have very
small but nonzero virtualities up to about $Q_{\mathrm{max}}^{2}\approx 0.08$
GeV2 for protons and $Q^{2}_{\text{max}}\approx 8\cdot 10^{-4}$ GeV2 for Pb
nuclei, and as a result the lepton pairs feature a steeply falling spectrum
reaching up to $p_{\text{T}}^{\ell\ell}\approx 0.4$ and 0.2 GeV for p-p, and
Pb-Pb UPCs, respectively, with very small (but nonzero) acoplanarity tails. In
general, the experimental selection criteria applied to select (semi)exclusive
photon-photon events require $p_{\text{T}}^{\ell\ell}$ and $A_{\phi}$ very
close to zero, though not exactly null because of detector smearing effects
and, thus, it is more realistic to generate MC events that also account for
the small virtuality-induced variations of the $p_{\text{T}}^{\ell\ell}$ and
$A_{\phi}$ distributions. To do so, we explicitly include in our gamma-UPC
setup the initial $\gamma$ virtualities through a small extra $k_{\perp}$
implemented by directly modifying the 4-momenta of the incoming photons and
outgoing particles in the Les Houches (lhe) file output of the generated MC
events. Such a “$k_{\perp}$ smearing” of the initial and final states is
technically implemented as explained in detail in App. A. For each independent
lhe event, we sample the virtualities of the transverse momenta of the two
initial photons according to Eq. (1), and the final four-momenta of the
produced particles are reshuffled to guarantee the onshell and momentum
conservation conditions event-by-event, as done similarly in Ref.
Frixione:2019fxg . We note that, in general, such a $k_{\perp}$-smearing
procedure can only be applied directly onto the lhe dilepton events output
when computed at LO accuracy. This is so because at NLO, the different Born,
virtual, real, and infrared-divergence-subtraction counterterms are strongly
correlated, and one cannot just smear the initial photon $k_{\perp}$ of each
contribution independently. However, when a particular photon-photon final
state does not have IR-divergence cancellations between the NLO real and
virtual contributions, such as in the LbL $\gamma\gamma\to\gamma\gamma$ case
AH:2023kor ; AH:2023ewe , the NLO lhe files can indeed be smeared out
following the method outlined above. A procedure to combine, in an effective
manner, the impact of the initial photon $k_{\perp}$ distribution and the NLO
corrections in the generated $p_{\text{T}}^{\ell\ell}$ and $A_{\phi}$
distributions is also presented in App. A.

For the determination of NLO QED corrections, we proceed as follows. Since the
initial photons are almost onshell, we work in the so-called $\alpha(0)$
renormalization scheme, found suitable for final states with similar
topologies Pagani:2021iwa . In order to avoid the complications of “dressed”
massless leptons (affected by the collinear $\gamma$ emission)
Frederix:2016ost , the masses of the charged leptons are kept nonzero (since
the electron mass is very small, this prevents us from considering the NLO
corrections to the B–W process, $\gamma\gamma\to\mathrm{e}^{+}\mathrm{e}^{-}$,
in this work). Because of the large hierarchies between the hard scale given
by the invariant mass of the lepton pairs and the low scale of the lepton
masses, the code has been optimized to avoid numerical instabilities. The
relevant ultraviolet and rational $R_{2}$ counterterms of the NLO QED
amplitudes calculation have been coded in the Universal Feynman Output format
Degrande:2011ua ; Darme:2023jdn  generated with an ad-hoc Mathematica program.
In the $\alpha(0)$ scheme, thanks to the Ward identity, only external lepton
masses are relevant, while all other charged fermion masses that enter into
the renormalization constant of $\alpha$ cancel out exactly with the
corresponding terms in the wavefunction renormalization constants of the
external photons. The numerical results have been obtained with
$\alpha=1/137.036$, and $m_{\mu}=0.10566$ GeV and $m_{\tau}=1.77686$ GeV
masses ParticleDataGroup:2022pth .

The calculations of the $\gamma\gamma\to\ell^{+}\ell^{-}$ cross sections
presented here have been obtained with the latest version (v1.6) of the gamma-
UPC code gammaUPCv2 , which includes several improvements compared to the v1.0
version Shao:2022cly , such as updated photon flux of the proton to account
for corrections beyond the dipole elastic form factor, and recalculated ion
survival probabilities based on a parametrization of the nuclear overlap
functions dEnterria:2020dwq  computed with a MC Glauber code Loizides:2017ack
rather than based on analytical “optical Glauber” expressions. The associated
parametric uncertainties of the predictions (mostly due to the modeling of the
proton and ion survival probabilities) have been also estimated gammaUPCv2
and amount to 1–3% for the range of dilepton masses measured in UPCs at the
LHC, $m_{\ell\ell}\approx 5$–500 GeV, and are smaller than the LO versus NLO
differences discussed here.

##  III Two-photon dimuon production at NLO in QED

This section studies the case of exclusive dimuon production,
$\gamma\gamma\to\mu^{+}\mu^{-}$, in UPCs with protons and Pb ion beams. Such a
final state is considered as a clean standard-candle process that can be used
to calibrate our theoretical understanding of the photon fluxes and hadron
survival probabilities. At the LHC, the process has been precisely measured in
p-p CMS:2011vma ; ATLAS:2015wnx ; ATLAS:2017sfe ; CMS:2018uvs ; ATLAS:2020mve
and Pb-Pb ATLAS:2020epq ; CMS:2020skx  collisions. In Ref. Shao:2022cly , we
compared our gamma-UPC LO results with the ATLAS Pb-Pb UPC measurement at
$\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV ATLAS:2020epq , finding that the ChFF
predictions normalized to the total inclusive cross section reproduced better
the shape of the kinematic distributions compared to the EDFF case. There
exist also detailed comparisons of the LO calculations from the Superchic
Harland-Lang:2020veo  and Starlight Klein:2016yzr  programs with these ATLAS
data in Ref. Harland-Lang:2021ysd . Here, we confront our LO and NLO QED
predictions with the four existing measurements of
$\gamma\gamma\to\mu^{+}\mu^{-}$ cross sections in Pb-Pb and p-p UPCs at the
LHC performed within the fiducial phase spaces listed in Table 1.

Table 1: Exclusive dimuon measurements at the LHC (compared with the LO and NLO predictions computed in this work) with their kinematic phase space indicated. System, experiment |  $\gamma\gamma\to\mu^{+}\mu^{-}$ fiducial phase space
---|---
p-p at 7 TeV, CMS CMS:2011vma  |  $p_{\text{T}}^{\,\mu}>4$ GeV, $|\eta^{\,\mu}|<2.1$, $m_{\mu\mu}>11.5$ GeV
p-p at 7 TeV, ATLAS ATLAS:2015wnx  |  $p_{\text{T}}^{\,\mu}>10$ GeV, $|\eta^{\,\mu}|<2.4$, $m_{\mu\mu}>20$ GeV
p-p at 13 TeV, ATLAS ATLAS:2017sfe  |  $p_{\text{T}}^{\,\mu}>6$ GeV, $|\eta^{\,\mu}|<2.4$, $m_{\mu\mu}=12\mbox{--}30$ GeV
$p_{\text{T}}^{\,\mu}>10$ GeV, $|\eta^{\,\mu}|<2.4$,
$m_{\mu\mu}=30\mbox{--}70$ GeV
Pb-Pb at 5.02 TeV, ATLAS ATLAS:2020epq  |  $p_{\text{T}}^{\,\mu}>4$ GeV, $|\eta^{\,\mu}|<2.4$, $m_{\mu\mu}>10$ GeV

### Exclusive two-photon dimuon production in p-p collisions at the LHC

For the p-p UPC case, measurements have been carried out so far by CMS at
$\sqrt{s_{{}_{\text{NN}}}}=7$ TeV CMS:2011vma  and by ATLAS at
$\sqrt{s_{{}_{\text{NN}}}}=7,\,13$ TeV ATLAS:2015wnx ; ATLAS:2017sfe . The
muon pseudorapidity $\eta^{\,\mu}$ selections are relatively similar in both
experiments, but the CMS measurement is carried out at somehow lower
$p_{\text{T}}^{\mu}$ and dimuon masses. Table 2 presents our LO and NLO
fiducial cross sections ($\sigma^{\mathrm{LO}}$ and $\sigma^{\mathrm{NLO}}$,
respectively) compared with the four experimental results. Our default results
are those obtained with ChFF fluxes, but we provide for comparison in
parenthesis the EDFF values, which are about 12–20% smaller than those derived
with the former. The first observation is that the NLO corrections decrease
the theoretical cross section by 3.5–5.0% with respect to the LO results
(independently of the photon flux used). The impact of the NLO corrections
increases with $m_{\mu\mu}$, due to the mass-dependence of the quasi-collinear
$\alpha\ln{m_{\mu}^{2}/m_{\mu\mu}^{2}}$ logarithm that affects the FSR
contribution. The experimental uncertainties on the cross sections are
relatively small, of the order of 3–15% depending on the bin, so the precision
of the measurement is commensurate with the size of the NLO corrections. In
general, the NLO ChFF predictions agree with the data (as indicated by the
$\sigma^{\mathrm{data}}/\sigma^{\mathrm{NLO}}\approx 1$ ratio in the last
column of Table 2), whereas the LO ChFF results tend to overestimate the data
by 4–10%, and the NLO EDFF calculations to undershoot the measurements by
about 10%. These results further confirm that the ChFF photon flux is more
realistic than the EDFF one, which is already theoretically justified since
the latter has an arbitrary infrared cutoff on the UPC impact parameter.

Table 2: Fiducial exclusive dimuon cross sections measured in UPCs at the LHC (within the phase space defined in Table 1), compared with the theoretical LO and NLO QED results obtained with gamma-UPC using the ChFF and EDFF $\gamma$ fluxes. The last column lists the corresponding $\sigma^{\mathrm{data}}/\sigma^{\mathrm{NLO}}$ ratios.  $\gamma\gamma\to\mu^{+}\mu^{-}$ | measured $\sigma^{\mathrm{data}}$ |  gamma-UPC $\sigma^{\mathrm{LO}}$ |  gamma-UPC $\sigma^{\mathrm{NLO}}$ | ratio $\sigma^{\mathrm{data}}/\sigma^{\mathrm{NLO}}$
---|---|---|---|---
System, experiment |  | ChFF (EDFF) | ChFF (EDFF) | ChFF (EDFF)
p-p at 7 TeV, CMS CMS:2011vma  |  $3.38^{+0.62}_{-0.59}$ pb |  $3.62~{}(3.20)$ pb |  $3.50~{}(3.10)$ pb |  $0.97^{+0.18}_{-0.17}$ ($1.09^{+0.20}_{-0.19}$)
p-p at 7 TeV, ATLAS ATLAS:2015wnx  |  $0.628\pm 0.038$ pb |  $0.687~{}(0.59)$ pb |  $0.653~{}(0.56)$ pb |  $0.96\pm 0.06$ ($1.12\pm 0.07$)
p-p at 13 TeV, ATLAS ATLAS:2017sfe  |  $3.12\pm 0.16$ pb |  $3.23~{}(2.88)$ pb |  $3.09~{}(2.76)$ pb |  $1.00\pm 0.05$ ($1.13\pm 0.06$)
Pb-Pb at 5.02 TeV, ATLAS ATLAS:2020epq  |  $34.1\pm 0.8$ $\mu$b |  $39.4~{}(31.5)$ $\mu$b |  $37.5~{}(30.0)$ $\mu$b |  $0.91\pm 0.02$ ($1.14\pm 0.03$)

Figure 2: Distributions of dimuon acoplanarity $A_{\phi}^{\mu\mu}$ (left) and
$p_{\text{T}}^{\mu\mu}$ (right) in p-p UPCs at $\sqrt{s}=13$ TeV (within the
phase space defined in Table 1, third row) computed with gamma-UPC (ChFF
photon flux) at LO$\,+\,k_{\perp}$-smearing (blue dashed histograms) and NLO
QED (red histograms) accuracy, and their corresponding combined distributions
(filled grey histograms).

It is important to notice that the exclusive $\gamma\gamma\to\mu^{+}\mu^{-}$
cross sections in p-p UPCs are usually obtained from a fit of the acoplanarity
distribution measured in data over the $A_{\phi}\approx 0$–0.1 range in order
to remove the contributions from semiexclusive (single and double proton
dissociation) processes constrained by templates from MC simulations. It is,
thus, important to have a good theoretical control of the spectral shape of
the signal acoplanarity. Although the large majority of the signal is clearly
located at very small acoplanarities ($A_{\phi}\lesssim 0.02$), Fig. 2 (left)
shows that the NLO QED corrections are needed to properly describe the dimuon
acoplanarity tail beyond $A_{\phi}\approx 0.04$. Similarly, the NLO
contributions enhance also the tail of the dimuon transverse momentum, beyond
$p_{\text{T}}^{\mu\mu}\gtrsim 1$ GeV (Fig. 2, right). Any experimental
criteria applied on the dimuon $p_{\text{T}}^{\mu\mu}$ distribution has an
impact on the determination of (a fraction of) the full exclusive cross
section and, thus, requires the use of NLO predictions.

Figure 3: Differential cross section for exclusive dimuon production as a
function of invariant mass, $\mathrm{d}\sigma/\mathrm{d}m_{\mu^{+}\mu^{-}}$,
in p-p UPCs at $\sqrt{s}=13$ TeV within the fiducial cuts of Table 1 (third
row). The data (dots with orange boxes) ATLAS:2017sfe  are compared with LO
(blue dashed histogram) and NLO (red histogram) QED predictions with EDFF
(left) and ChFF (right) photon fluxes. The lower panels of each plot show the
corresponding data/NLO ratios (dots with orange uncertainties boxes) and
NLO/LO $K$ factors (violet line).

Besides the fiducial integrated cross section, the ATLAS analysis
ATLAS:2017sfe  provided measurements with small uncertainties of the
differential dimuon cross section as a function of pair invariant mass over
$m_{\mu\mu}=12$–70 GeV. Figure 3 shows the data (black dots with orange
uncertainties boxes) compared to our LO (blue dashed histogram) and NLO (red
histogram) predictions with the EDFF (left) and ChFF (right) photon fluxes.
The top panels show the data–theory comparison of the distributions, and the
bottom panels present the data/NLO ratios (dots with orange uncertainties
boxes) and the $K=$ NLO/LO factors (violet line). The $K\approx 0.945$–0.965
factors are relatively constant as a function of $m_{\mu\mu}$ and independent
of the photon flux type. Both, the NLO corrections and ChFF fluxes improve the
theoretical reproduction of the data. A data–theory goodness-of-fit yields
$\chi^{2}/N_{\mathrm{dof}}\approx 1.2,\,1.8$ (for $N_{\mathrm{dof}}=4$ degrees
of freedom) for the NLO predictions obtained with the ChFF and EDFF fluxes,
respectively. This result further confirms the optimized choice of ChFF photon
fluxes plus NLO corrections for accurate calculations of dilepton production
in $\gamma\gamma$ processes.

### Exclusive two-photon dimuon production in Pb-Pb collisions at the LHC

Exclusive dimuon production has been also measured in Pb-Pb UPCs, whose cross
section are enhanced by a $Z^{4}\approx 4.5\cdot 10^{7}$ factor with respect
to their p-p counterpart. The ATLAS experiment has measured the fiducial cross
section of $\sigma(\gamma\gamma\to\mu^{+}\mu^{-})=34.1\pm
0.3\,(\mathrm{stat}.)\pm 0.7\,(\mathrm{syst}.)~{}\mu$b in Pb-Pb at 5.02 TeV,
reaching a 2.3% precision ATLAS:2020epq . As shown in the last row of Table 2,
the LO ChFF cross section overshoots the data by 13%, whereas the LO EDFF
cross section undershoots them by 8%. The NLO QED corrections reduce the LO
cross sections by about 5% and, as a conclusion, the NLO ChFF (EDFF) deviates
from the ATLAS data by $+9\%$ ($-13\%)$. All in all, the NLO$\,+\,$ChFF result
features the closest agreement with the experimental fiducial cross section.
Figure 4 shows the impact of the NLO corrections on the dimuon acoplanarity
$A_{\phi}^{\mu\mu}$ (left plot) and pair $p_{\text{T}}^{\mu\mu}$ (right plot).
As found in the p-p case, the higher-order corrections extend the tails of
both distributions to values much above those obtained with
LO$\,+\,k_{\perp}$-smearing calculations. In the left plot of Fig. 4, the
predicted acoplanarity distributions are compared to the
$\mathrm{d}\sigma/\mathrm{d}A_{\phi}^{\mu\mu}$ measurement of Ref.
ATLAS:2020epq  (black dots with orange uncertainties), and the ratio of data
over merged NLO$\,+\,k_{\perp}$-smearing prediction is shown on the bottom
panel. As expected, the LO$\,+\,k_{\perp}$-smearing result (blue dashed
histogram) describes well the back-to-back region ($A_{\phi}^{\mu\mu}\lesssim
5\cdot 10^{-3}$), but clearly misses the experimental data in the tail above
$A_{\phi}^{\mu\mu}\approx 10^{-2}$. On the other hand, the NLO-alone
calculation (red histogram) exhibits the opposite behaviour. A proper
combination of both regimes (as discussed in App. A) can reproduce the data
within $\pm 20\%$ over the full acoplanarity range, as displayed in the bottom
left ratio. Those results emphasize once more the need to properly account for
NLO corrections in the extraction of the exclusive dimuon signal from the
experimental Pb-Pb UPC data.

Figure 4: Distributions for dimuon acoplanarity $A_{\phi}^{\mu\mu}$ (left) and
$p_{\text{T}}^{\mu\mu}$ (right) in Pb-Pb UPCs at
$\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV computed with gamma-UPC (ChFF photon
flux) at LO$\,+\,k_{\perp}$-smearing (blue dashed histograms) and NLO (red
histogram) QED accuracy, and their corresponding combined distributions
(filled grey histograms). The acoplanarity distributions are compared to the
$\mathrm{d}\sigma/\mathrm{d}A_{\phi}$ measurement of Ref. ATLAS:2020epq
(black dots with orange uncertainties), and the ratio of data over merged
NLO$\,+\,k_{\perp}$-smearing prediction is shown in the bottom left panel.

In the data-theory comparisons of differential distributions in Pb-Pb UPCs
carried out in Ref. Shao:2022cly , we concluded that the normalized shapes of
the measured spectra favoured the ChFF photon flux over the EDFF one, albeit
the overall normalization of the LO ChFF prediction appeared somewhat worse
than the LO EDFF one. We examine now how such a conclusion is altered by
including NLO QED corrections. The first differential cross section considered
is that as a function of the invariant mass of the muon pair shown in Fig. 5.
From left to right, the three panels present the data-theory comparisons for
three different pair rapidity intervals ($|y_{\mu\mu}|<0.8$,
$0.8<|y_{\mu\mu}|<1.6$, and $1.6<|y_{\mu\mu}|<2.4$, respectively), with the
upper (lower) plots showing the ChFF (EDFF) predictions. The top panels of
each plot show the differential cross sections in data and at LO and NLO
accuracies, and the bottom panels display the data/NLO ratios (black dots with
orange uncertainties) and NLO/LO $K$ factors (violet histogram). The first
observation is that the NLO ChFF results are in better agreement with data
than the EDFF ones for all kinematic ranges: the data/NLO(ChFF) ratios are
close to unity for the three measured distributions. The second observation is
that the NLO QED corrections are small in the low $m_{\mu\mu}$ range, but
become increasingly relevant and reduce the differential cross section by up
to 15% at the highest dimuon invariant mass. This can be understood from the
aforementioned $m_{\mu\mu}$ dependence of the quasi-collinear logarithms
arising in the final state radiation. This result shows that it would be
incorrect to apply a global constant $K$ factor to account for NLO corrections
in the exclusive dimuon mass distributions.

Figure 5: Differential cross section as a function of invariant mass,
$\mathrm{d}\sigma/\mathrm{d}m_{\mu^{+}\mu^{-}}$, for exclusive dimuon
production in Pb-Pb UPCs at $\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV in different
ranges of dimuon rapidity. The data (black dots with orange uncertainties)
ATLAS:2020epq  are compared with LO (blue dashed histogram) and NLO (red
histogram) QED gamma-UPC predictions with ChFF (top) and EDFF (bottom)
$\gamma$ fluxes. The lower panels of each plot show the corresponding data/NLO
ratios (black dots with orange uncertainties) and NLO/LO $K$ factors (violet
histogram).

Figure 6: Differential cross section as a function of dimuon absolute
rapidity, $\mathrm{d}\sigma/\mathrm{d}|y_{\mu^{+}\mu^{-}}|$, for exclusive
dimuon production in Pb-Pb UPCs at $\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV in
different ranges of invariant mass. The data (black dots with orange
uncertainties) ATLAS:2020epq  are compared with LO (blue dashed histogram) and
NLO (red histogram) QED gamma-UPC predictions with ChFF (top) and EDFF
(bottom) $\gamma$ fluxes. The lower panels of each plot show the corresponding
data/NLO ratios (black dots with orange uncertainties) and NLO/LO $K$ factors
(violet histogram).

The second differential cross section of interest is that as a function of
absolute pair rapidity, $\mathrm{d}\sigma/\mathrm{d}|y_{\mu\mu}|$. The ATLAS
study ATLAS:2020epq  revealed systematic deviations between the measured
distribution, in particular in the tail, and the LO Starlight predictions. Our
comparison of the data to the LO and NLO distributions with ChFF and EDFF
fluxes is shown in Fig. 6 with the same layout as in Fig. 5. The three panels
correspond to three invariant mass intervals:
$m_{\mu\mu}\in[10,20],[20,40],[40,80]$ GeV, from left to right. Again, the NLO
ChFF predictions agree nicely with the data for all kinematic ranges within
the experimental uncertainties, whereas the NLO EDFF calculations
underestimate them, especially for increasing $|y_{\mu\mu}|$. However, at
variance with the behaviour seen in the invariant mass distribution, the $K$
factors are quite flat in $|y_{\mu\mu}|$.

The last differential cross section studied is as a function of the scattering
angle $\theta^{\star}_{\mu\mu}$ of a muon in the rest frame of the dimuon
system. The absolute value of the cosine of $\theta^{\star}_{\mu\mu}$ can be
obtained from the rapidity separation of the two muons via the relationship
$|\cos{\theta^{\star}_{\mu\mu}}|=\tanh{\left(\Delta y/2\right)}$, where
$\Delta y=y^{\,\mu^{+}}-y^{\,\mu^{-}}$. Figure 7 shows this angular
distribution over the full pair rapidity coverage ($|y_{\mu\mu}|<2.4$) in
three dimuon invariant mass windows: $10<m_{\mu\mu}<20$ GeV,
$20<m_{\mu\mu}<40$ GeV, and $40<m_{\mu\mu}<80$ GeV, from left to right. The
rapid depletion of the distribution approaching
$|\cos{\theta^{\star}_{\mu\mu}}|=1$ is due to the experimental fiducial cuts,
in particular on $\eta^{\,\mu}$, that remove events with the muons flying
along the beam direction. As found in previous differential cross sections,
the NLO ChFF predictions feature a much better agreement with the data than
the NLO EDFF ones, as indicated by the data/NLO ratios for the former centered
at unity in all bottom panels. The $K$ factor has nontrivial dependencies on
both $m_{\mu\mu}$ and $|\cos{\theta^{\star}_{\mu\mu}}|$. In the lowest mass
window, the $K$ factor decreases from 0.96 at
$|\cos{\theta^{\star}_{\mu\mu}}|=0$ to 0.88 at
$|\cos{\theta^{\star}_{\mu\mu}}|\approx 0.92$. The experimental fiducial
criteria remove the population of events in the region of
$|\cos{\theta^{\star}_{\mu\mu}}|\gtrsim 0.92$. In the $20<m_{\mu\mu}<40$ GeV
range, the $K$ factor is quite flat as a function of
$|\cos{\theta^{\star}_{\mu\mu}}|$ up to
$|\cos{\theta^{\star}_{\mu\mu}}|\gtrsim 0.8$. Finally, in the highest mass
regime, the $K$ factor increases from 0.9 at
$|\cos{\theta^{\star}_{\mu\mu}}|=0$ to 0.93 for
$|\cos{\theta^{\star}_{\mu\mu}}|>0.95$. Such a nontrivial kinematic dependence
calls for a rigorous NLO QED computation as presented in this paper for any
more precise future studies of the $\theta^{\star}_{\mu\mu}$ scattering angle
in data.

Figure 7: Differential cross section as a function of the cosine of the
scattering angle in the rest frame of the dimuon system,
$\mathrm{d}\sigma/\mathrm{d}|\cos{\theta^{\star}_{\mu\mu}}|$, for exclusive
dimuon production in Pb-Pb UPCs at $\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV (over
$|y_{\mu\mu}|<2.4$). The data (black dots with orange uncertainties)
ATLAS:2020epq  are compared with LO (blue dashed histogram) and NLO (red
histogram) QED gamma-UPC predictions with ChFF (top) and EDFF (bottom)
$\gamma$ fluxes. The lower panels of each plot show the corresponding data/NLO
ratios (black dots with orange uncertainties) and NLO/LO $K$ factors (violet
histogram).

The energies $E_{\gamma,1}$ and $E_{\gamma,2}$ of the two initial photons in
the process $\gamma\gamma\to\mu^{+}\mu^{-}$ can be inferred from the four-
momenta of the final-state muons by energy conservation. At LO accuracy, they
can be obtained from the dimuon mass and rapidity as
$E_{\gamma,1}=(m_{\mu\mu}/2)\exp{(+y_{\mu\mu})}$ and
$E_{\gamma,2}=(m_{\mu\mu}/2)\exp{(-y_{\mu\mu})}$, respectively. At NLO
accuracy, both expressions are not fully exact but remain useful for most of
the cases. Following the ATLAS analysis, we define the maximum and minimum
photon energies in an event as
$E_{\gamma,\mathrm{max}}=\mathrm{max}(E_{\gamma,1},E_{\gamma,2})$ and
$E_{\gamma,\mathrm{min}}=\mathrm{min}(E_{\gamma,1},E_{\gamma,2})$,
respectively. The distributions of $E_{\gamma,\mathrm{min}}$ and
$E_{\gamma,\mathrm{max}}$ measured in data are shown in Fig. 8 left and right,
respectively, compared to our theoretical LO and NLO predictions obtained with
the ChFF (top) and EDFF (bottom) photon fluxes. The measured
$E_{\gamma,\mathrm{min,max}}$ distribution are well reproduced by the NLO ChFF
calculations except in the last bin333Note that the similar data-theory
comparison shown at LO accuracy in our previous work Shao:2022cly  (Fig. 10
left of the paper) is incorrect in the first $E_{\gamma,\mathrm{min}}$ bin, as
we had wrongly normalized the distribution by a 2-GeV (instead of 1.7 GeV) bin
size. where the theory underestimates the data although the latter suffers
from large statistical uncertainties. On the other hand, the calculations with
the EDFF $\gamma$ flux undershoot both initial $E_{\gamma,\mathrm{min}}$ and
$E_{\gamma,\mathrm{max}}$ photon energy distributions over all their measured
ranges, as observed already at LO accuracy in the comparison of the ATLAS data
to the Starlight MC predictions ATLAS:2020epq , which are based on this photon
flux. The $K$ factors for both $E_{\gamma,\mathrm{min}}$ and
$E_{\gamma,\mathrm{max}}$ spectra, shown in the lower panels of both plots,
feature a monotonic decrease with energy, from $K\approx 0.95$ down to 0.85
and 0.91 respectively, due again to the effects of the quasi-collinear FSR
logarithms.

Figure 8: Differential cross section as a function of the the minimum
($\mathrm{d}\sigma/dE_{\gamma,\mathrm{min}}$, left) and maximum
($\mathrm{d}\sigma/dE_{\gamma,\mathrm{max}}$, right) energy of the initial
photon in exclusive dimuon production in Pb-Pb UPCs at
$\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV. The data (black dots with orange
uncertainties) ATLAS:2020epq  are compared with LO (blue dashed histogram) and
NLO (red histogram) QED gamma-UPC predictions with ChFF (top) and EDFF
(bottom) $\gamma$ fluxes. The lower panels of each plot show the corresponding
data/NLO ratios (black dots with orange uncertainties) and NLO/LO $K$ factors
(violet histogram).

In summary, the detailed data-theory comparisons presented here for four
differential observables measured in $\gamma\gamma\to\mu^{+}\mu^{-}$ processes
in p-p and Pb-Pb UPCs at the LHC, confirm that the NLO QED calculations with
the ChFF photon flux can accurately describe them in all kinematic regions of
phase space probed. The measured dimuon masses (Figs. 3 and 5), absolute
rapidity (Fig. 6), $|\cos{\theta^{\star}_{\mu\mu}}|$ (Fig. 7) and the
$E_{\gamma,\mathrm{min,max}}$ (Fig. 8) distributions agree nicely with the NLO
prediction obtained with the ChFF $\gamma$ flux, whereas the calculations with
the EDFF flux systematically undershoot the experimental results in all
regions of phase space.

##  IV Two-photon ditau production at NLO in QED

Let us start by considering $\gamma\gamma\to\tau^{+}\tau^{-}$ collisions at
$\mathrm{e}^{+}\mathrm{e}^{-}$ colliders where the process was first studied.
More concretely, we consider the measurement of exclusive ditau in photon-
photon collisions performed by the DELPHI experiment from which the previous
best limits on $(g-2)_{\tau}$ had been derived DELPHI:2003nah . In this study,
the extraction of $(g-2)_{\tau}$ from the measured
$\gamma\gamma\to\tau^{+}\tau^{-}$ cross section was performed with a LO EPA
calculation, with freely floating $(g-2)_{\tau}$ values fitted to the data,
where the upper limit of the integration over 4-momenta of the emitted photons
($Q^{2}_{\mathrm{max}}$) was chosen such that the fiducial cross section
predicted by the EPA (for the SM value of the anomalous electromagnetic
moment) agreed with the more theoretically accurate radcor prediction
Berends:1986if . This latter code calculated the process at LO QED accuracy
including FSR corrections on the electron and positron lines (second diagram
in Fig. 1). We compute here the cross sections for
$\mathrm{e}^{+}\mathrm{e}^{-}\xrightarrow[]{\gamma\gamma}\mathrm{e}^{+}\mathrm{e}^{-}\tau^{+}\tau^{-}$
at LO and NLO accuracy within our setup using the default Weizsäcker-Williams
EPA photon flux of MG5_aMC. Our cross section calculations implement the
original fiducial selection criteria for the signal, mostly applied to remove
highly virtual $\gamma\gamma$ processes: $m_{\tau\tau}<40$ GeV, $\theta^{\rm
e+}<10^{\circ}$, $\theta^{\rm e-}<2^{\circ}$ (where $\theta$ is the angle
between incoming and outgoing $\mathrm{e}^{+}\mathrm{e}^{-}$, ranging from 0
to 180∘) for the various c.m. energies considered in the DELPHI measurement:
$\sqrt{s}=182.7,\,188.7,\,197.6,\,$ and 206.3 GeV. Within this phase space,
and limiting the maximum virtuality of the colliding photons to values
$Q^{2}_{\mathrm{max}}\approx(10$–15 GeV)2 as in the original EPA approach
employed by the DELPHI analysis, we obtain cross sections consistent with the
experimental collision-energy/luminosity-weighted value of
$\sigma(\gamma\gamma\to\tau^{+}\tau^{-})=429\pm 17$ pb. Interestingly, we find
that the NLO corrections are positive for this process and augment the LO
cross sections by a factor of about 0.9%. Although the DELPHI experimental
uncertainty of 4% is larger than the size of the NLO corrections, and though
probably a leading fraction of the latter were already accounted for in the
original LO$\,+\,$FSR radcor calculations, it would have been more
theoretically accurate to go beyond the real radiative corrections and take
also into account any potential cancellations between NLO real and virtual
terms as considered here. This result underscores the importance of accounting
for higher-order QED terms for future precise extractions of $(g-2)_{\tau}$
via $\gamma\gamma\to\tau^{+}\tau^{-}$ measurements at
$\mathrm{e}^{+}\mathrm{e}^{-}$ colliders such as at Belle-II Belle-II:2018jsg
and, in the longer term, FCC-ee FCC:2018evy ; Dam:2018rfz , with expected
subpercent or even subpermil experimental uncertainties.

To date, there exist two measurements of $\gamma\gamma\to\tau^{+}\tau^{-}$
production cross sections in UPCs at the LHC. First, at low ditau invariant
mass ($m_{\tau\tau}\gtrsim 5$ GeV) using PbPb data at 5.02 TeV CMS:2022arf
and, secondly, at higher masses ($m_{\tau\tau}>50$ GeV) in p-p collisions at
13 TeV CMS:2024skm . The measured fiducial cross sections have, however, still
large uncertainties in the 16%–30% range. We have transformed the experimental
fiducial cross sections (including convoluted $\eta,p_{\text{T}}$ acceptance
criteria applied at the detector-level to the different leptonic and hadronic
decay products of both tau leptons) into more inclusive cross sections that
are simpler to compare to the corresponding gamma-UPC predictions at LO and
NLO QED accuracy with ChFF and EDFF photon fluxes. The experimental and
theoretical results are listed in Table 3. Based on the data–theory
comparisons of Ref. CMS:2022arf , we find that the extrapolation of the
measured Pb-Pb fiducial result,
$\sigma^{\mathrm{CMS}}_{\mathrm{fid}}(\gamma\gamma\to\tau^{+}\tau^{-})=4.8\pm
0.8$ fb, to a fully inclusive cross section depends on whether one uses the
calculations of Ref. Beresford:2019gww  or Dyndal:2020yen . In the former
case, one obtains $570\pm 100$ $\mu$b, whereas in the latter case the result
is $850\pm 150$ $\mu$b. Given the relatively poor experimental precision and
the additional apparent extrapolation uncertainty, we quote the range span by
both values. In addition, we provide predictions for
$\gamma\gamma\to\tau^{+}\tau^{-}$ in p-p collisions at $\sqrt{s}=13.6$ TeV in
the phase space region of very high masses ($m_{\tau\tau}>300$ GeV) where both
protons can be tagged by the forward near-beam spectrometers CMS:2018uvs .

Table 3: Exclusive ditau cross sections in Pb-Pb UPCs at $\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV and in p-p collisions at $\sqrt{s}=13$ TeV ($m_{\tau\tau}>50$ GeV) and 13.6 TeV ($m_{\tau\tau}>500$ GeV) experimentally measured and obtained with gamma-UPC at LO and NLO QED accuracy using the EDFF and ChFF $\gamma$ fluxes. The first column lists the experimental measurements of Refs. CMS:2022hly ; CMS:2024skm  extrapolated to the broader phase space considered in the calculations. The last column lists the corresponding NLO-over-LO $K$ factors.  $\gamma\gamma\to\tau^{+}\tau^{-}$ | Measured $\sigma^{\mathrm{data}}$ |  gamma-UPC $\sigma^{\mathrm{LO}}$ |  gamma-UPC $\sigma^{\mathrm{NLO}}$ |  $K$ factor
---|---|---|---|---
Colliding system (kinematic cuts) | (extrapolated) | ChFF (EDFF) | ChFF (EDFF) | $\sigma^{\mathrm{NLO}}/\sigma^{\mathrm{LO}}$
Pb-Pb at 5.02 TeV (inclusive) | 580–850 $\mu$b |  $1060$ $(860)$ $\mu$b |  $1070$ $(870)$ $\mu$b | $1.01$
p-p at 13 TeV ($m_{\tau\tau}>50$ GeV) |  $855_{-215}^{+260}$ fb |  $900$ $(730)$ fb |  $895$ $(725)$ fb | $0.995$
p-p at 13.6 TeV ($m_{\tau\tau}>300$ GeV, $|y^{\tau}|<2.5$) | – |  $1.35$ $(1.02)$ fb |  $1.31$ $(0.98)$ fb | $0.965$

On the one hand, as similarly found for the $\gamma\gamma\to\mu^{+}\mu^{-}$
case (Table 2), the theoretical cross sections derived with the ChFF fluxes
are about 20% (30%) larger than those obtained with the EDFF ones in Pb-Pb
(p-p) UPCs. The difference between the cross sections obtained with both
photon fluxes rises with ditau mass because the ChFF spectrum becomes harder
than the EDFF one for increasing $E_{\gamma}$ energies Shao:2022cly . On the
other hand, the impact of the NLO corrections on the theoretical cross
sections is only of about $+1\%$ at the low $m_{\tau\tau}$ masses probed in
PbPb UPCs, about $-0.5\%$ at the intermediate $m_{\tau\tau}$ probed in p-p
UPCs, whereas at high masses it is more significant and leads to a reduction
by about $-3.5\%$ of the LO $\sigma(\gamma\gamma\to\tau^{+}\tau^{-})$
predictions. We note that for very large ditau masses, weak boson corrections
Godunov:2023myj  should be of the same size (i.e., percent level) as the QED
corrections and, therefore, would need to be also taken into account, but we
leave their calculation for future work. Unfortunately, the current accuracy
and precision of the experimental fiducial cross sections precludes any
detailed conclusion on the level of data–theory agreement, but one should
expect that the outcome of the more precise total and differential
$\gamma\gamma\to\mu^{+}\mu^{-}$ cross section studies discussed in the
previous section, i.e., that ChFF-plus-NLO constitutes the state-of-the-art
prediction, should hold also for the $\gamma\gamma\to\tau^{+}\tau^{-}$
process.

The impact of NLO corrections is more visible in the
$\gamma\gamma\to\tau^{+}\tau^{-}$ differential cross sections as a function of
pair invariant mass, $m_{\tau\tau}$, and tau $p_{\text{T}}$ as shown for Pb-Pb
UPCs at $\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV in Fig. 9, and for high-mass p-p
UPCs at $\sqrt{s_{{}_{\text{NN}}}}=13.6$ TeV in Fig. 10. We see that although
the integrated cross sections change by a few percent going from LO to NLO,
the differential distributions show in general small enhancements in the $K$
factors at low ditau masses and $p_{\text{T}}^{\tau}$ values that are then
compensated by larger reductions of the $K$-factor for increasing
$m_{\tau\tau}$ and tau $p_{\text{T}}$, by up to a factor of $-6\%$ in the tail
of the distributions. Since both kinematic distributions are the most
sensitive ones for the extraction of $(g-2)_{\tau}$ through comparison of the
experimental data to the predictions, it is important to have a good
theoretical control of their shape accounting for the NLO corrections
presented here.

Figure 9: Predicted $\gamma\gamma\to\tau^{+}\tau^{-}$ differential cross
sections at NLO QED accuracy as a function of ditau mass
$\mathrm{d}\sigma/\mathrm{d}m_{\tau\tau}$ (left) and
$\mathrm{d}\sigma/\mathrm{d}p_{\text{T}}^{\tau}$ (right) in Pb-Pb UPCs at
$\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV (for $m_{\tau\tau}>4$ GeV,
$|y^{\tau}|<2.5$) using the ChFF (red) and EDFF (blue) photon fluxes. The
lower plots show the corresponding
$\mathrm{d}\sigma^{\mathrm{NLO}}/\mathrm{d}\sigma^{\mathrm{LO}}$ $K$ factors.

Figure 10: Predicted $\gamma\gamma\to\tau^{+}\tau^{-}$ differential cross
sections at NLO QED accuracy as a function of ditau mass
$\mathrm{d}\sigma/\mathrm{d}m_{\tau\tau}$ (left) and
$\mathrm{d}\sigma/\mathrm{d}p_{\text{T}}^{\tau}$ (right) in p-p UPCs at
$\sqrt{s}=13.6$ TeV (for $m_{\tau\tau}>50$ GeV, $|y^{\tau}|<2.5$) using the
ChFF (red) and EDFF (blue) photon fluxes. The lower plots show the
corresponding $\mathrm{d}\sigma^{\mathrm{NLO}}/\mathrm{d}\sigma^{\mathrm{LO}}$
$K$ factors.

##  V Summary

We have presented theoretical calculations for the production cross sections
of muon and tau pairs in photon-photon collisions,
$\gamma\gamma\to\mu^{+}\mu^{-},\tau^{+}\tau^{-}$, at next-to-leading (NLO)
accuracy in quantum electrodynamics (QED). The total and differential cross
sections have been computed in the equivalent photon approximation (EPA)
implemented in the gamma-UPC$\,+\,$MadGraph5_aMC@NLO codes. The theoretical
predictions have been confronted to the existing exclusive
$\gamma\gamma\to\mu^{+}\mu^{-},\tau^{+}\tau^{-}$ data from ultraperipheral
collisions (UPCs) of proton (p) and lead (Pb) ions at the LHC, as well as to
results from $\mathrm{e}^{+}\mathrm{e}^{-}$ collisions at LEP for the
$\gamma\gamma\to\tau^{+}\tau^{-}$ process. For UPCs at the LHC, two types of
EPA photon fluxes, as well as associated survival probabilities of the photon-
emitting hadrons, based on the charge (ChFF) and electric-dipole (EDFF) form
factors have been tested. Three main findings are worth highlighting.

First, the NLO terms are found to modify the total fiducial cross sections by
up to 5%, as well as to increase the tails of the dilepton acoplanarity
$A_{\phi}^{\ell\ell}$ and transverse momentum $p_{\text{T}}^{\ell\ell}$
distributions, with respect to the leading-order (LO) predictions including
the effects of the (very small) virtuality of the colliding photons. Since the
extraction of the $\gamma\gamma\to\ell^{+}\ell^{-}$ signals is based on
applying selection criteria on the $A_{\phi}^{\ell\ell}$ and
$p_{\text{T}}^{\ell\ell}$ distributions, and/or fitting them to Monte Carlo
templates, the improved NLO calculations have quantitative implications for
the experimental cross section determinations, which are being measured with
increasing experimental precision (few percent uncertainties) at the LHC. The
NLO corrections also modify the differential cross sections as a function of
dilepton masses $m_{\ell\ell}$, rapidity $y_{\ell\ell}$, and
$|\cos{\theta^{\star}_{\ell\ell}}|$. In general, the corrections are small
(and, in some cases, positive) at low dilepton masses and lepton
$p_{\text{T}}^{\ell}$, but become increasingly large (and negative) for higher
$m_{\ell\ell}$ and $p_{\text{T}}^{\ell}$ values, with up to a 15% depletions
observed in the tails of some kinematic distributions compared to the LO
predictions. This latter behaviour is due to the
$\alpha\ln{m_{\ell}^{2}/m_{\ell\ell}^{2}}$ dependence of the quasi-collinear
logarithm that controls the final-state $\gamma$ radiation contribution. Those
results will have implications for upcoming, more precise, experimental
measurements of the anomalous magnetic moment of the tau lepton,
$(g-2)_{\tau}$, via the $\gamma\gamma\to\tau^{+}\tau^{-}$ process at different
$m_{\tau^{+}\tau^{-}}$ regions in UPCs at the LHC.

Second, the detailed UPC data-theory comparisons clearly favour the ChFF
photon fluxes over the EDFF ones, and the former should be used in the
calculations when confronting the data to the predictions for accurate
phenomenological interpretation. With the ChFF flux, one can also compute
nonexclusive photon-photon collisions with proton or nucleus overlap at NLO
accuracy. For the $\gamma\gamma\to\mu^{+}\mu^{-}$ processes, the distributions
of dimuon masses and rapidities, $|\cos{\theta^{\star}_{\mu\mu}}|$, and
$E_{\gamma,\mathrm{min,max}}$ agree nicely with the NLO prediction obtained
with the ChFF $\gamma$ flux, whereas the calculations with the EDFF flux
systematically underestimate the experimental results in all regions of phase
space.

Third, for photon-photon collisions at LEP, the NLO corrections for the
$\gamma\gamma\to\tau^{+}\tau^{-}$ final state measured within the DELPHI phase
space are positive and augment the LO cross sections by a factor of about 1%.
Although the 4% LEP experimental uncertainty is larger than the size of the
NLO modifications, it would have been more theoretically accurate to consider
them in order to extract $(g-2)_{\tau}$ from the data. Such considerations are
of particular relevance for future precise determinations of $(g-2)_{\tau}$
via $\gamma\gamma\to\tau^{+}\tau^{-}$ measurements at
$\mathrm{e}^{+}\mathrm{e}^{-}$ colliders, such as at Belle-II and FCC-ee, with
expected subpercent or even subpermil experimental uncertainties.

In summary, the inclusion of NLO QED terms (and the use of ChFF photon fluxes
for hadron beams) in the calculations of $\gamma\gamma\to\ell^{+}\ell^{-}$
processes improves the data–theory agreement, and proves an indispensable
ingredient for the extraction of precision quantities, in particular
$(g-2)_{\tau}$, and/or to search for deviations from the SM predictions in the
absolute or differential cross section due to elusive new physics
contributions in future measurements of $\gamma\gamma$ processes in UPCs at
the LHC and at $\mathrm{e}^{+}\mathrm{e}^{-}$ colliders. With the developments
presented in this work, we have extended the number of photon-fusion physics
processes calculable at NLO QED (or QCD) theoretical accuracy with the gamma-
UPC code including, to date, the exclusive production of muon and tau leptons,
light-by-light scattering, and top-quark pairs.

#### Acknowledgments.—

Support from the European Union’s Horizon 2020 research and innovation program
(grant agreement No. 824093, STRONG-2020, EU Virtual Access “NLOAccess”), the
ERC (grant 101041109 ‘BOSON’), the French ANR (grant ANR-20-CE31-0015,
“PrecisOnium”), and the French LIA FCPPN, are acknowledged. Views and opinions
expressed are however those of the authors only and do not necessarily reflect
those of the European Union or the European Research Council Executive Agency.
Neither the European Union nor the granting authority can be held responsible
for them. For the purpose of Open Access, a CC-BY public copyright license has
been applied by the authors to the present document and will be applied to all
subsequent versions up to the Author Accepted Manuscript arising from this
submission.

## References

  * (1) D. Morgan, M. R. Pennington, and M. R. Whalley, “A Compilation of data on two photon reactions leading to hadron final states,” J. Phys. G 20 Suppl. 8A (1994) A1–A147.
  * (2) M. R. Whalley, “A Compilation of data on two photon reactions,” J. Phys. G 27 (2001) A1–A121.
  * (3) M. Przybycien, “Two-photon physics at LEP,” Nucl. Phys. B Proc. Suppl. 179-180 (2008) 54–61.
  * (4) C. A. Bertulani, S. R. Klein, and J. Nystrand, “Physics of ultra-peripheral nuclear collisions,” Ann. Rev. Nucl. Part. Sci. 55 (2005) 271, arXiv:nucl-ex/0502005.
  * (5) A. J. Baltz, “The physics of ultraperipheral collisions at the LHC,” Phys. Rept. 458 (2008) 1, arXiv:0706.3356 [nucl-ex].
  * (6) J. de Favereau de Jeneret, V. Lemaitre, Y. Liu, S. Ovyn, T. Pierzchala, K. Piotrzkowski, X. Rouby, N. Schul, and M. Vander Donckt, “High energy photon interactions at the LHC,” arXiv:0908.2020 [hep-ph].
  * (7) C. F. von Weizsäcker, “Radiation emitted in collisions of very fast electrons,” Z. Phys. 88 (1934) 612.
  * (8) E. J. Williams, “Nature of the high-energy particles of penetrating radiation and status of ionization and radiation formulae,” Phys. Rev. 45 (1934) 729.
  * (9) S. J. Brodsky, T. Kinoshita, and H. Terazawa, “Two photon mechanism of particle production by high-energy colliding beams,” Phys. Rev. D 4 (1971) 1532.
  * (10) V. M. Budnev, I. F. Ginzburg, G. V. Meledin, and V. G. Serbo, “The two photon particle production mechanism. Physical problems. Applications. Equivalent photon approximation,” Phys. Rept. 15 (1975) 181.
  * (11) H.-S. Shao and D. d’Enterria, “gamma-UPC: automated generation of exclusive photon-photon processes in ultraperipheral proton and nuclear collisions with varying form factors,” JHEP 09 (2022) 248, arXiv:2207.03012 [hep-ph].
  * (12) STAR Collaboration, J. Adam et al., “Low-$p_{T}$ $e^{+}e^{-}$ pair production in Au$+$Au collisions at $\sqrt{s_{NN}}$ = 200 GeV and U$+$U collisions at $\sqrt{s_{NN}}$ = 193 GeV at STAR,” Phys. Rev. Lett. 121 (2018) 132301, arXiv:1806.02295 [hep-ex].
  * (13) ATLAS Collaboration, M. Aaboud et al., “Observation of centrality-dependent acoplanarity for muon pairs produced via two-photon scattering in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with the ATLAS detector,” Phys. Rev. Lett. 121 (2018) 212301, arXiv:1806.08708 [nucl-ex].
  * (14) ALICE Collaboration, S. Acharya et al., “Dielectron production at midrapidity at low transverse momentum in peripheral and semi-peripheral Pb–Pb collisions at ${\sqrt{s}}_{\textrm{NN}}$ = 5.02 TeV,” JHEP 06 (2023) 024, arXiv:2204.11732 [nucl-ex].
  * (15) CMS Collaboration, A. Hayrapetyan et al., “Observation of $\gamma\gamma\to\tau\tau$ in proton-proton collisions and limits on the anomalous electromagnetic moments of the $\tau$ lepton,” arXiv:2406.03975 [hep-ex].
  * (16) FP420 R&D Collaboration, M. G. Albrow et al., “The FP420 R&D Project: Higgs and New Physics with forward protons at the LHC,” JINST 4 (2009) T10001, arXiv:0806.0302 [hep-ex].
  * (17) CMS, TOTEM Collaboration, M. Albrow et al., “CMS-TOTEM Precision Proton Spectrometer.” CERN-LHCC-2014-021, http://cds.cern.ch/record/1753795, 2014.
  * (18) ATLAS Collaboration, M. Taševský, “Status of the AFP project in the ATLAS experiment,” AIP Conf. Proc. 1654 (2015) 090001.
  * (19) CMS, TOTEM Collaboration, A. Tumasyan et al., “Proton reconstruction with the CMS-TOTEM Precision Proton Spectrometer,” JINST 18 (2023) P09009, arXiv:2210.05854 [hep-ex].
  * (20) R. Bruce et al., “New physics searches with heavy-ion collisions at the CERN Large Hadron Collider,” J. Phys. G 47 (2020) 060501, arXiv:1812.07688 [hep-ph].
  * (21) CMS, TOTEM Collaboration, A. M. Sirunyan et al., “Observation of proton-tagged, central (semi)exclusive production of high-mass lepton pairs in pp collisions at 13 TeV with the CMS-TOTEM precision proton spectrometer,” JHEP 07 (2018) 153, arXiv:1803.04496 [hep-ex].
  * (22) ATLAS Collaboration, G. Aad et al., “Observation and measurement of forward proton scattering in association with lepton pairs produced via the photon fusion mechanism at ATLAS,” Phys. Rev. Lett. 125 (2020) 261801, arXiv:2009.14537 [hep-ex].
  * (23) G. Breit and J. A. Wheeler, “Collision of two light quanta,” Phys. Rev. 46 (1934) 1087.
  * (24) JADE Collaboration, W. Bartel et al., “Lepton pair production in double tagged two photon interactions,” Z. Phys. C 30 (1986) 545.
  * (25) L3 Collaboration, M. Acciarri et al., “Production of e, $\mu$ and $\tau$ pairs in untagged two photon collisions at LEP,” Phys. Lett. B 407 (1997) 341–350.
  * (26) C. R. Vane et al., “Electron positron pair production in Coulomb collisions of ultrarelativistic sulphur ions with fixed targets,” Phys. Rev. Lett. 69 (1992) 1911.
  * (27) CERES/NA45 Collaboration, R. Bauer et al., “Measurement of electromagnetically produced $\mathrm{e}^{+}\mathrm{e}^{-}$ pairs in distant S-Pt collisions,” Phys. Lett. B 332 (1994) 471.
  * (28) STAR Collaboration, J. Adams et al., “Production of $\mathrm{e}^{+}\mathrm{e}^{-}$ pairs accompanied by nuclear dissociation in ultra-peripheral heavy ion collision,” Phys. Rev. C 70 (2004) 031902, arXiv:nucl-ex/0404012.
  * (29) STAR Collaboration, J. Adam et al., “Measurement of $\mathrm{e}^{+}\mathrm{e}^{-}$ momentum and angular distributions from linearly polarized photon collisions,” Phys. Rev. Lett. 127 (2021) 052302, arXiv:1910.12400 [nucl-ex].
  * (30) PHENIX Collaboration, S. Afanasiev et al., “Photoproduction of $\mathrm{J}/\psi$ and of high mass $\mathrm{e}^{+}\mathrm{e}^{-}$ in ultra-peripheral Au+Au collisions at $\sqrt{s}=200$ GeV,” Phys. Lett. B 679 (2009) 321, arXiv:0903.2041 [nucl-ex].
  * (31) CDF Collaboration, A. Abulencia et al., “Observation of exclusive electron-positron production in hadron-hadron collisions,” Phys. Rev. Lett. 98 (2007) 112001, arXiv:hep-ex/0611040.
  * (32) CDF Collaboration, T. Aaltonen et al., “Search for exclusive Z boson production and observation of high mass $\rm p\bar{p}\to\gamma\gamma\to p+\ell\ell+\bar{p}$ events in $\rm p\bar{p}$ collisions at $\sqrt{s}$ = 1.96 TeV,” Phys. Rev. Lett. 102 (2009) 222002, arXiv:0902.2816 [hep-ex].
  * (33) ALICE Collaboration, E. Abbas et al., “Charmonium and $\mathrm{e}^{+}\mathrm{e}^{-}$ pair photoproduction at mid-rapidity in ultra-peripheral Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV,” Eur. Phys. J. C 73 (2013) 2617, arXiv:1305.1467 [nucl-ex].
  * (34) CMS Collaboration, S. Chatrchyan et al., “Search for exclusive or semi-exclusive photon pair production and observation of exclusive and semi-exclusive electron pair production in pp collisions at $\sqrt{s}=7$ TeV,” JHEP 11 (2012) 080, arXiv:1209.1666 [hep-ex].
  * (35) CMS Collaboration, A. M. Sirunyan et al., “Evidence for light-by-light scattering and searches for axion-like particles in ultraperipheral PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV,” Phys. Lett. B 797 (2019) 134826, arXiv:1810.04602 [hep-ex].
  * (36) CMS Collaboration, A. Hayrapetyan et al., “Measurements of the light-by-light scattering and the Breit–Wheeler processes, and searches for axion-like particles in ultraperipheral PbPb collisions at 5.02 TeV.” CMS-PAS-HIN-21-015, 2024.
  * (37) ATLAS Collaboration, G. Aad et al., “Measurement of exclusive $\gamma\gamma\to\ell^{+}\ell^{-}$ production in proton-proton collisions at $\sqrt{s}=7$ TeV with the ATLAS detector,” Phys. Lett. B 749 (2015) 242, arXiv:1506.07098 [hep-ex].
  * (38) ATLAS Collaboration, M. Aaboud et al., “Evidence for light-by-light scattering in heavy-ion collisions with the ATLAS detector at the LHC,” Nature Phys. 13 (2017) 852, arXiv:1702.01625 [hep-ex].
  * (39) ATLAS Collaboration, G. Aad et al., “Exclusive dielectron production in ultraperipheral Pb+Pb collisions at $\sqrt{s_{\textrm{NN}}}$ = 5.02 TeV with ATLAS,” JHEP 2306 (2023) 182, arXiv:2207.12781 [nucl-ex].
  * (40) S. Klein, A. H. Mueller, B.-W. Xiao, and F. Yuan, “Lepton Pair Production Through Two Photon Process in Heavy Ion Collisions,” Phys. Rev. D (2020) 094013, arXiv:2003.02947 [hep-ph].
  * (41) J. D. Brandenburg, J. Seger, Z. Xu, and W. Zha, “Report on progress in physics: observation of the Breit–Wheeler process and vacuum birefringence in heavy-ion collisions,” Rept. Prog. Phys. 86 (2023) 083901, arXiv:2208.14943 [hep-ph].
  * (42) D. d’Enterria and G. G. da Silveira, “Observing light-by-light scattering at the Large Hadron Collider,” Phys. Rev. Lett. 111 (2013) 080405, arXiv:1305.7142 [hep-ph]. [Erratum: Phys.Rev.Lett. 116, 129901 (2016)].
  * (43) ATLAS Collaboration, G. Aad et al., “Observation of light-by-light scattering in ultraperipheral Pb+Pb collisions with the ATLAS detector,” Phys. Rev. Lett. 123 (2019) 052001, arXiv:1904.03536 [hep-ex].
  * (44) V. A. Khoze, A. D. Martin, R. Orava, and M. G. Ryskin, “Luminosity monitors at the LHC,” Eur. Phys. J. C 19 (2001) 313–322, arXiv:hep-ph/0010163.
  * (45) A. G. Shamov and V. I. Telnov, “Precision luminosity measurement at LHC using two photon production of $\mu^{+}\mu^{-}$ pairs,” Nucl. Instrum. Meth. A 494 (2002) 51–56, arXiv:hep-ex/0207095.
  * (46) J. A. M. Vermaseren, “Two Photon Processes at Very High-Energies,” Nucl. Phys. B 229 (1983) 347.
  * (47) F. A. Berends, P. H. Daverveldt, and R. Kleiss, “Monte Carlo simulation of two photon processes. 1. Radiative corrections to multiperipheral $\mathrm{e}^{+}\mathrm{e}^{-}\to\mu^{+}\mu^{-}$ production,” Comput. Phys. Commun. 40 (1986) 271–284.
  * (48) PLUTO Collaboration, C. Berger et al., “Tagged two photon production of muon pairs,” Z. Phys. C 27 (1985) 249.
  * (49) TOPAZ Collaboration, H. Hayashii et al., “An experimental study of muon pair production in tagged two photon interactions,” Phys. Lett. B 279 (1992) 422–428.
  * (50) OPAL Collaboration, R. Akers et al., “A Study of muon pair production and evidence for tau pair production in photon-photon collisions at LEP,” Z. Phys. C 60 (1993) 593–600.
  * (51) L3 Collaboration, P. Achard et al., “Muon pair and tau pair production in two photon collisions at LEP,” Phys. Lett. B 585 (2004) 53–62, arXiv:hep-ex/0402037.
  * (52) H1 Collaboration, A. Aktas et al., “Muon pair production in ep collisions at HERA,” Phys. Lett. B 583 (2004) 28–40, arXiv:hep-ex/0311015.
  * (53) CDF Collaboration, T. Aaltonen et al., “Observation of exclusive charmonium production and $\gamma+\gamma\to\mu^{+}\mu^{-}$ in $\rm p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV,” Phys. Rev. Lett. 102 (2009) 242001, arXiv:0902.1271 [hep-ex].
  * (54) CMS Collaboration, S. Chatrchyan et al., “Exclusive photon-photon production of muon pairs in proton-proton collisions at $\sqrt{s}=7$ TeV,” JHEP 01 (2012) 052, arXiv:1111.5536 [hep-ex].
  * (55) ATLAS Collaboration, M. Aaboud et al., “Measurement of the exclusive $\gamma\gamma\to\mu^{+}\mu^{-}$ process in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” Phys. Lett. B 777 (2018) 303, arXiv:1708.04053 [hep-ex].
  * (56) ALICE Collaboration, S. Acharya et al., “Exclusive and dissociative J/$\psi$ photoproduction, and exclusive dimuon production, in p-Pb collisions at sNN=8.16 TeV,” Phys. Rev. D 108 (2023) 112004, arXiv:2304.12403 [nucl-ex].
  * (57) CMS Collaboration, A. M. Sirunyan et al., “Observation of forward neutron multiplicity dependence of dimuon acoplanarity in ultraperipheral Pb-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV,” Phys. Rev. Lett. 127 (2021) 122001, arXiv:2011.05239 [hep-ex].
  * (58) ATLAS Collaboration, G. Aad et al., “Exclusive dimuon production in ultraperipheral Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with ATLAS,” Phys. Rev. C 104 (2021) 024906, arXiv:2011.12211 [nucl-ex].
  * (59) S. R. Klein, J. Nystrand, J. Seger, Y. Gorbunov, and J. Butterworth, “STARlight: A Monte Carlo simulation program for ultra-peripheral collisions of relativistic ions,” Comput. Phys. Commun. 212 (2017) 258, arXiv:1607.03838 [hep-ph].
  * (60) L. A. Harland-Lang, M. Tasevsky, V. A. Khoze, and M. G. Ryskin, “A new approach to modelling elastic and inelastic photon-initiated production at the LHC: SuperChic 4,” Eur. Phys. J. C 80 (2020) 925, arXiv:2007.12704 [hep-ph].
  * (61) F. Cornet and J. I. Illana, “Tau pair production via photon-photon collisions at LEP,” Phys. Rev. D 53 (1996) 1181–1184, arXiv:hep-ph/9503466.
  * (62) F. del Aguila, F. Cornet, and J. I. Illana, “The Possibility of using a large heavy ion collider for measuring the electromagnetic properties of the tau-lepton,” Phys. Lett. B 271 (1991) 256.
  * (63) S. Atag and A. A. Billur, “Possibility of determining $\tau$ lepton electromagnetic moments in $\gamma\gamma\to\tau^{+}\tau^{-}$ process at the CERN-LHC,” JHEP 11 (2010) 060, arXiv:1005.2841 [hep-ph].
  * (64) L. Beresford and J. Liu, “New physics and tau $g-2$ using LHC heavy ion collisions,” Phys. Rev. D 102 (2020) 113008, arXiv:1908.05180 [hep-ph].
  * (65) M. Dyndal, M. Klusek-Gawenda, M. Schott, and A. Szczurek, “Anomalous electromagnetic moments of $\tau$ lepton in $\gamma\gamma\to\tau^{+}\tau^{-}$ reaction in Pb+Pb collisions at the LHC,” Phys. Lett. B 809 (2020) 135682, arXiv:2002.05503 [hep-ph].
  * (66) DELPHI Collaboration, J. Abdallah et al., “Study of tau-pair production in photon-photon collisions at LEP and limits on the anomalous electromagnetic moments of the tau lepton,” Eur. Phys. J. C 35 (2004) 159, arXiv:hep-ex/0406010.
  * (67) CMS Collaboration, A. Tumasyan et al., “Observation of $\tau$ lepton pair production in ultraperipheral lead-lead collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV,” Phys. Rev. Lett. 131 (2023) 151803, arXiv:2206.05192 [nucl-ex].
  * (68) ATLAS Collaboration, G. Aad et al., “Observation of the $\gamma\gamma\to\tau^{+}\tau^{-}$ process in Pb+Pb collisions and constraints on the $\tau$-lepton anomalous magnetic moment with the ATLAS detector,” Phys. Rev. Lett. 131 (2023) 151802, arXiv:2204.13478 [hep-ex].
  * (69) J. S. Schwinger, “On Quantum Electrodynamics and the magnetic moment of the electron,” Phys. Rev. 73 (1948) 416–417.
  * (70) S. Eidelman and M. Passera, “Theory of the tau lepton anomalous magnetic moment,” Mod. Phys. Lett. A 22 (2007) 159–179, arXiv:hep-ph/0701260.
  * (71) T. Aoyama et al., “The anomalous magnetic moment of the muon in the Standard Model,” Phys. Rept. 887 (2020) 1–166, arXiv:2006.04822 [hep-ph].
  * (72) A. A H, E. Chaubey, M. Fraaije, V. Hirschi, and H.-S. Shao, “Light-by-light scattering at next-to-leading order in QCD and QED,” Phys. Lett. B 851 (2024) 138555, arXiv:2312.16956 [hep-ph].
  * (73) A. A H, E. Chaubey, and H.-S. Shao, “Two-loop massive QCD and QED helicity amplitudes for light-by-light scattering,” JHEP 03 (2024) 121, arXiv:2312.16966 [hep-ph].
  * (74) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 07 (2014) 079, arXiv:1405.0301 [hep-ph].
  * (75) R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H. S. Shao, and M. Zaro, “The automation of next-to-leading order electroweak calculations,” JHEP 07 (2018) 185, arXiv:1804.10017 [hep-ph]. [Erratum: JHEP 11, 085 (2021)].
  * (76) H.-S. Shao, “HELAC-Onia: An automatic matrix element generator for heavy quarkonium physics,” Comput. Phys. Commun. 184 (2013) 2562–2570, arXiv:1212.5293 [hep-ph].
  * (77) H.-S. Shao, “HELAC-Onia 2.0: an upgraded matrix-element and event generator for heavy quarkonium physics,” Comput. Phys. Commun. 198 (2016) 238–259, arXiv:1507.03435 [hep-ph].
  * (78) S. Frixione, B. Fuks, V. Hirschi, K. Mawatari, H.-S. Shao, P. A. Sunder, and M. Zaro, “Automated simulations beyond the Standard Model: supersymmetry,” JHEP 12 (2019) 008, arXiv:1907.04898 [hep-ph].
  * (79) D. Pagani, H.-S. Shao, I. Tsinikos, and M. Zaro, “Automated EW corrections with isolated photons: $t\overline{t}\gamma$, $t\overline{t}\gamma\gamma$ and t$\gamma$j as case studies,” JHEP 09 (2021) 155, arXiv:2106.02059 [hep-ph].
  * (80) R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H.-S. Shao, and M. Zaro, “The complete NLO corrections to dijet hadroproduction,” JHEP 04 (2017) 076, arXiv:1612.06548 [hep-ph].
  * (81) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer, and T. Reiter, “UFO - The Universal FeynRules Output,” Comput. Phys. Commun. 183 (2012) 1201–1214, arXiv:1108.2040 [hep-ph].
  * (82) L. Darmé et al., “UFO 2.0: the ‘Universal Feynman Output’ format,” Eur. Phys. J. C 83 (2023) 631, arXiv:2304.09883 [hep-ph].
  * (83) Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * (84) N. Crepet, D. d’Enterria, and H.-S. Shao, “Improved modeling of high-energy photon-photon collisions at hadron colliders.” to be submitted.
  * (85) D. d’Enterria and C. Loizides, “Progress in the Glauber model at collider energies,” Ann. Rev. Nucl. Part. Sci. 71 (2021) 315, arXiv:2011.14909 [hep-ph].
  * (86) C. Loizides, J. Kamin, and D. d’Enterria, “Improved Monte Carlo Glauber predictions at present and future nuclear colliders,” Phys. Rev. C 97 (2018) 054910, arXiv:1710.07098 [nucl-ex]. [Erratum: Phys. Rev. C 99 (2019) 019901].
  * (87) L. A. Harland-Lang, V. A. Khoze, and M. G. Ryskin, “Elastic photon-initiated production at the LHC: the role of hadron-hadron interactions,” SciPost Phys. 11 (2021) 064, arXiv:2104.13392 [hep-ph].
  * (88) Belle-II Collaboration, W. Altmannshofer et al., “The Belle II Physics Book,” PTEP 2019 (2019) 123C01, arXiv:1808.10567 [hep-ex]. [Erratum: PTEP 2020, 029201 (2020)].
  * (89) FCC Collaboration, A. Abada et al., “FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2,” Eur. Phys. J. ST 228 no. 2, (2019) 261–623.
  * (90) M. Dam, “Tau-lepton Physics at the FCC-ee circular e+e- Collider,” SciPost Phys. Proc. 1 (2019) 041, arXiv:1811.09408 [hep-ex].
  * (91) S. I. Godunov, E. K. Karkaryan, V. A. Novikov, A. N. Rozanov, M. I. Vysotsky, and E. V. Zhemchugov, “Weak interaction corrections to muon pair production via the photon fusion at the LHC,” Phys. Rev. D 108 (2023) 093006, arXiv:2308.01169 [hep-ph].

##  Appendix A Initial photon virtualities in the gamma-UPC code

The virtuality dependence of the EPA photon fluxes emitted by a charge $Z$ is
given by Eq. (1). Such an equation can be more simply written as a function of
the photon virtuality $Q^{2}$ and longitudinal energy $E_{\gamma}$ (or,
equivalently, fraction of the beam momentum $x_{\gamma}$ carried by the
photon) alone, $N_{\gamma/Z}^{\mathrm{ChFF}}(x_{\gamma},Q^{2})$, using two
similar expressions. The first expression is obtained by taking $Q^{2}\equiv
k_{\perp}^{2}+\tfrac{E_{\gamma}^{2}}{\gamma_{\mathrm{L}}^{2}}=k_{\perp}^{2}+m_{\mathrm{N}}^{2}x_{\gamma}^{2}$
and $Q_{\mathrm{min}}^{2}=m_{\mathrm{N}}^{2}x_{\gamma}^{2}$, and reads

$\displaystyle\mathrm{d}N_{\gamma/Z}^{\mathrm{ChFF}}(x_{\gamma},Q^{2})$
$\displaystyle=$
$\displaystyle\frac{Z^{2}\alpha}{\pi}\frac{\mathrm{d}x_{\gamma}}{x_{\gamma}}\frac{\mathrm{d}Q^{2}}{Q^{2}}\left(1-\frac{Q^{2}_{\mathrm{min}}}{Q^{2}}\right)\left[F_{\mathrm{ch}}(Q)\right]^{2}.$
(2)

The second form uses the exact expression for the virtuality,
$Q^{2}=\tfrac{k_{\perp}^{2}+E_{\gamma}^{2}/\gamma_{\mathrm{L}}^{2}}{1-x_{\gamma}}=\tfrac{k_{\perp}^{2}+m_{\mathrm{N}}^{2}x_{\gamma}^{2}}{1-x_{\gamma}}$,
and is obtained by replacing
$\tfrac{E_{\gamma}^{2}}{\gamma_{\mathrm{L}}^{2}}=\left(x_{\gamma}m_{\mathrm{N}}\right)^{2}\overset{x_{\gamma}\ll
1}{\simeq}Q_{\mathrm{min}}^{2}\equiv\tfrac{m_{\mathrm{N}}^{2}x_{\gamma}^{2}}{1-x_{\gamma}}$,
and defining
$(1-x_{\gamma})Q^{2}=k_{\perp}^{2}+\tfrac{E_{\gamma}^{2}}{\gamma_{\mathrm{L}}^{2}}$,
thus $k_{\perp}^{2}=(1-x_{\gamma})(Q^{2}-Q^{2}_{\mathrm{min}})$, which leads
to

$\displaystyle\mathrm{d}N_{\gamma/Z}^{\mathrm{ChFF}}(x_{\gamma},Q^{2})$
$\displaystyle=$
$\displaystyle\frac{Z^{2}\alpha}{\pi}\frac{\mathrm{d}x_{\gamma}}{x_{\gamma}}\frac{\mathrm{d}Q^{2}}{Q^{2}}\left(1-\frac{Q^{2}_{\mathrm{min}}}{Q^{2}}\right)\left[F_{\mathrm{ch}}\left(\sqrt{1-x_{\gamma}}Q\right)\right]^{2}.$
(3)

Both expressions are identical in the case $x_{\gamma}\ll 1$, which is
fulfilled for most of the phenomenologically relevant phase space in
$\gamma\gamma$ collisions at colliders. Indeed, a typical $\gamma\gamma\to X$
collision process with $m_{X}\approx 5$ GeV at midrapidity in PbPb collisions
at the LHC probes $x_{\gamma}\approx 10^{-3}$ values, whereas the same
production of a final state with $m_{X}\approx 100$ GeV (corresponding to
$x_{\gamma}\approx 0.1$) has a cross section two orders-of-magnitude smaller
Shao:2022cly .

As explained in Section II, the matrix element calculations in MG5_aMC expect
collinear (i.e., $k_{\perp}$-independent) input photon distributions, and
therefore when running gamma-UPC \+ MG5_aMC together to compute a given
$\gamma\gamma$ process, the photon fluxes need to be fully integrated over $Q$
and $k_{\perp}$. In order to restore the full unintegrated dependencies of the
photon fluxes, our gamma-UPC setup incorporates a small extra $k_{\perp}$
implemented by directly modifying the 4-momenta of the incoming photons and,
correspondingly, outgoing produced particles in the Les Houches (lhe) file
output of the generated MC events. Such a “$k_{\perp}$ smearing” of the
initial and final states is performed by running a python script on the lhe
file that modifies the kinematics of each event as explained next.

Let us consider a generic photon-photon collision process
$\gamma_{1}\gamma_{2}\to X_{3}\cdots X_{n+2}$, with $n$ particles produced in
the final state. First, the incoming photons transverse momentum
$k_{\perp}\approx Q$ is sampled from the distribution given by Eq. (2) or (3)
depending on the choice of the `ion_Form` setting of the script, with uniform
azimuthal plane dependence. Then, the kinematics of the final particles is
appropriately reshuffled to respect momentum conservation and the on-mass-
shell conditions. Technically, this is implemented as in Ref. Frixione:2019fxg
where the 4-momenta of the two initial photons, in the c.m. frame of the two
initial colliding hadrons, are written as a function of the fractional
momentum $x_{\gamma}\equiv x$ (for simplicity, we omit the $\gamma$ subindex
in the following) carried out by each photon,

$\displaystyle k_{1}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{s_{{}_{\text{NN}}}}}{2}x_{1}\left(1,0,0,1\right),\quad
k_{2}=\frac{\sqrt{s_{{}_{\text{NN}}}}}{2}x_{2}\left(1,0,0,-1\right).$ (4)

The 4-momenta of the final particles, $k_{3},\ldots,k_{n+2}$, are then
required to satisfy the momentum conservation
$k_{1}+k_{2}=\sum_{j=3}^{n+2}{k_{j}}$ and onshell conditions
$k_{i}^{2}=m_{i}^{2}$, where $m_{1}=m_{2}=0$. After the momentum reshuffling,
we have

$\displaystyle k_{i}$ $\displaystyle\to$
$\displaystyle\bar{k}_{i},\quad\bar{k}_{1}^{2}=\bar{m}_{1}^{2}=-Q_{1}^{2},\quad\bar{k}_{2}^{2}=\bar{m}_{2}^{2}=-Q_{2}^{2},\quad\bar{k}_{j}^{2}=m_{j}^{2},j\geq
3.$ (5)

The new momentum conservation is
$\bar{k}_{1}+\bar{k}_{2}=\sum_{j=3}^{n+2}{\bar{k}_{j}}$. There are two
possible ways to proceed with the $k_{\perp}$ smearing of the collision
process: either starting reshuffling the momenta of the initial photons, or
modifying the kinematics of the final state:

  * •

Initial momentum reshuffling: This momentum reshuffling scheme applies for the
$n\geq 1$ case, i.e., for a photon-fusion process leading even to a single
one-particle final state. The invariant mass of the final-state system is
$\left(\sum_{j=3}^{n+2}{k_{j}}\right)^{2}=x_{1}x_{2}s_{{}_{\text{NN}}}$. Let
us use

$\displaystyle\bar{k}_{1}$ $\displaystyle=$
$\displaystyle\left(\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}\bar{x}_{1},\vec{\bar{k}}_{\perp,1},\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}\bar{x}_{1}\right),\quad\bar{k}_{2}=\left(\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}\bar{x}_{2},\vec{\bar{k}}_{\perp,2},-\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}\bar{x}_{2}\right),$
(6)

where $|\vec{\bar{k}}_{\perp,1}|=Q_{1}$ and $|\vec{\bar{k}}_{\perp,2}|=Q_{2}$.
For simplicity, we denote
$\vec{K}_{\perp}\equiv\vec{\bar{k}}_{\perp,1}+\vec{\bar{k}}_{\perp,2}$, where
we know that $(Q_{1}-Q_{2})^{2}\leq K_{\perp}^{2}\leq(Q_{1}+Q_{2})^{2}$. Now,
we must require
$(\bar{k}_{1}+\bar{k}_{2})^{2}=s_{{}_{\text{NN}}}\bar{x}_{1}\bar{x}_{2}-K_{\perp}^{2}=s_{{}_{\text{NN}}}x_{1}x_{2}$.
Thus,
$\bar{x}_{1}\bar{x}_{2}=x_{1}x_{2}+\tfrac{K_{\perp}^{2}}{s_{{}_{\text{NN}}}}$.
If $K_{\perp}^{2}\leq s_{{}_{\text{NN}}}(1-x_{1}x_{2})$, we can solve
$\bar{x}_{1}$ and $\bar{x}_{2}$ by additionally requiring
$\log{\tfrac{\bar{x}_{1}}{\bar{x}_{2}}}=\log{\tfrac{x_{1}}{x_{2}}}$. If there
is a physical solution, this solution is
$\bar{x}_{1}=\sqrt{\smash[b]{\tfrac{x_{1}}{x_{2}}\tfrac{K_{\perp}^{2}}{s_{{}_{\text{NN}}}}+x_{1}^{2}}}$
and
$\bar{x}_{2}=\sqrt{\smash[b]{\tfrac{x_{2}}{x_{1}}\tfrac{K_{\perp}^{2}}{s_{{}_{\text{NN}}}}+x_{2}^{2}}}$.
This gives us the final constraint on $K_{\perp}^{2}\leq
s_{{}_{\text{NN}}}\left(\mathrm{min}(1,\tfrac{x_{1}}{x_{2}},\tfrac{x_{2}}{x_{1}})-x_{1}x_{2}\right)$.
Otherwise, the values of $Q_{1}$ and $Q_{2}$ are regenerated. Lastly, the
whole final system needs to be boosted from the frame $k_{1}+k_{2}$ to
$\bar{k}_{1}+\bar{k}_{2}$. Since the fractional momentum $x$ changes, one can
decide whether to apply the additional rescaling weight,
$\frac{f_{1}(\bar{x}_{1})f_{2}(\bar{x}_{2})}{f_{1}(x_{1})f_{2}(x_{2})}$, or
not.

  * •

Final momentum reshuffling: We can only use this scheme for the $n\geq 2$
case, because we need to partition the final particles into two subsets. In
this scheme, let us keep $x_{1}$ and $x_{2}$, i.e.,

$\displaystyle\bar{k}_{1}$ $\displaystyle=$
$\displaystyle\left(\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}x_{1},\vec{\bar{k}}_{\perp,1},\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}x_{1}\right),\quad\bar{k}_{2}=\left(\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}x_{2},\vec{\bar{k}}_{\perp,2},-\tfrac{\sqrt{s_{{}_{\text{NN}}}}}{2}x_{2}\right)$
(7)

with $|\vec{\bar{k}}_{\perp,1}|=Q_{1}$ and $|\vec{\bar{k}}_{\perp,2}|=Q_{2}$.
Without loss of generality, we assume that the two subsets are
$\left\\{3,\ldots,{m+2}\right\\}$ and $\left\\{m+3,\ldots,n+2\right\\}$, where
$1\leq m\leq n-1$. Let us define $M_{1}^{2}\equiv
K_{1}^{2}\equiv\left(\sum_{j=3}^{m+2}{k_{j}}\right)^{2}$ and $M_{2}^{2}\equiv
K_{2}^{2}\equiv\left(\sum_{j=m+3}^{n+2}{k_{j}}\right)^{2}$, and keep the
invariant masses $M_{1}$ and $M_{2}$ of the two final systems unmodified. In
the case that $\bar{s}\equiv
x_{1}x_{2}s_{{}_{\text{NN}}}-K_{\perp}^{2}\leq(M_{1}+M_{2})^{2}$, we need to
regenerate $Q_{1}$ and $Q_{2}$ and
$\vec{K}_{\perp}=\vec{\bar{k}}_{\perp,1}+\vec{\bar{k}}_{\perp,2}$. Otherwise,
in the partonic rest frame, the energies of the two subsystems are
$\bar{K}_{1}^{0}=\varepsilon(\sqrt{\bar{s}},M_{1},M_{2})$ and
$\bar{K}_{2}^{0}=\varepsilon(\sqrt{\bar{s}},M_{2},M_{1})$, where
$\varepsilon(m,m_{1},m_{2})=\tfrac{m}{2}\left(1+\tfrac{m_{1}^{2}-m_{2}^{2}}{m^{2}}\right)$.
The mass on-shell conditions enforce
$\overrightarrow{\bar{K}}_{1}=\sqrt{\left(\bar{K}_{1}^{0}\right)^{2}-M_{1}^{2}}\tfrac{\overrightarrow{K}_{1}}{|\overrightarrow{K}_{1}|}$
and
$\overrightarrow{\bar{K}}_{2}=\sqrt{\left(\bar{K}_{2}^{0}\right)^{2}-M_{2}^{2}}\tfrac{\overrightarrow{K}_{2}}{|\overrightarrow{K}_{2}|}$.
Note that because
$\sqrt{\left(\bar{K}_{1}^{0}\right)^{2}-M_{1}^{2}}=\sqrt{\left(\bar{K}_{2}^{0}\right)^{2}-M_{2}^{2}}$
and in the partonic rest frame
$\tfrac{\overrightarrow{K}_{1}}{|\overrightarrow{K}_{1}|}=-\tfrac{\overrightarrow{K}_{2}}{|\overrightarrow{K}_{2}|}$,
it follows that
$\overrightarrow{\bar{K}}_{1}+\overrightarrow{\bar{K}}_{2}=\overrightarrow{0}$.
For $k_{j},3\leq j\leq m+2$, the $\bar{k}_{j}$ is obtained by first boosting
$k_{j}$ into the rest frame of $K_{1}$, and then boosting it back to the frame
of $\bar{K}_{1}$. A similar method is applied to obtain $\bar{k}_{j}$ (with
$m+3\leq j\leq n+2$) from $k_{j}$. Finally, the momenta of the final particles
are boosted from the $\left(\sqrt{\bar{s}},0,0,0\right)$ to
$\bar{k}_{1}+\bar{k}_{2}$ frames.

Since the initial-momentum reshuffling works out for the simplest production
of a single resonance in two-photon collisions, $\gamma\gamma\to X_{1}$, this
is the method currently implemented in our $k_{\perp}$ smearing python script.

The procedure to combine $\gamma\gamma\to\mathrm{X}$ kinematic distributions
including LO$\,+\,k_{\perp}$ smearing and NLO corrections has some arbitrary
freedom. The combined spectrum should follow closely the LO$\,+\,k_{\perp}$
result in the lowest values of the distributions, whereas the tails should
match exactly the NLO result. We adopt a rather simple approach here to
combine the intermediate region of the distributions where both contributions
overlap. For a given kinematic distribution differential in $O$ generated in a
given AA UPC, a weighted distribution is defined as follows

$\displaystyle\frac{\mathrm{d}\sigma^{\mathrm{AA}}_{\mathrm{merged}}}{\mathrm{d}O}$
$\displaystyle=$
$\displaystyle(1-w_{\mathrm{AA}}(O))\frac{\mathrm{d}\sigma^{\mathrm{AA}}_{\mathrm{NLO}}}{\mathrm{d}O}+w_{\mathrm{AA}}(O)\frac{\mathrm{d}\sigma^{\mathrm{AA}}_{\mathrm{LO+sm}}}{\mathrm{d}O},$
(8)

where the weights for the dilepton variables
$O=A_{\phi},\,p_{\text{T}}^{\ell\ell}$ are, respectively,

$\displaystyle w_{\mathrm{PbPb}}(A_{\phi})$
$\displaystyle=\frac{1}{1+(A_{\phi}/0.005)^{4}},\quad\quad\;\;\;w_{\mathrm{pp}}(A_{\phi})=\frac{1}{1+(A_{\phi}/0.03)^{4}},$
(9) $\displaystyle w_{\mathrm{PbPb}}(p_{\text{T}}^{\ell\ell})$
$\displaystyle=\frac{1}{1+(p_{\text{T}}^{\ell\ell}/0.01~{}\mathrm{GeV})^{4}},\quad
w_{\mathrm{pp}}(p_{\text{T}}^{\ell\ell})=\frac{1}{1+(p_{\text{T}}^{\ell\ell}/0.7~{}\mathrm{GeV})^{4}}.$

Figures 2 and 4 show distributions for dimuon acoplanarity $A_{\phi}^{\mu\mu}$
(left) and $p_{\text{T}}^{\mu\mu}$ (right) in p-p UPCs at $\sqrt{s}=13$ TeV
and Pb-Pb UPCs at $\sqrt{s_{{}_{\text{NN}}}}=5.02$ TeV, respectively, in the
fiducial phase space of Refs. ATLAS:2017sfe ; ATLAS:2020epq , computed with
gamma-UPC (ChFF photon flux) with LO$\,+\,k_{\perp}$-smearing (blue
histograms) and NLO QED (red histograms) accuracy, as well as their
corresponding combined distributions (filled grey histograms) following Eqs.
(8) and (9). The agreement with the acoplanarity distribution measured in Pb-
Pb UPCs is good as shown in the bottom left panel of Fig. 4.