# Novel Cross Section Ratios as Possible Signals of Saturation in UPCs

Yuri V. Kovchegov<EMAIL_ADDRESS>Department of Physics, The Ohio State
University, Columbus, Ohio 43210, USA  Huachen Sun
<EMAIL_ADDRESS>Department of Physics, The Ohio State University,
Columbus, Ohio 43210, USA  Zhoudunming Tu<EMAIL_ADDRESS>Department of
Physics, Brookhaven National Laboratory, Upton, New York 11973, USA

###### Abstract

We propose new cross section ratios in ultra-peripheral A$+$A and $p+$A
collisions (UPCs) as a potential experimental signal of saturation physics. We
consider the ratio $R_{1}$ of elastic vector meson photoproduction cross
section to the inclusive hadron or jet photoproduction cross section. The
ratio can be measured in the $\gamma+$A and $\gamma+p$ collisions taking place
in the UPCs. We label the ratios $R_{1}({\rm A})$ and $R_{1}(p)$,
respectively. Constructing the double ratio $R_{\rm{UPC}}({\rm A})=R_{1}({\rm
A})/R_{1}(p)$, and performing a small-$x$ calculation both in the quasi-
classical approximation and by including small-$x$ evolution, we observe that
$R_{\rm{UPC}}({\rm A})$ exhibits a markedly different dependence on the
nuclear atomic number inside and outside the saturation region. This result
indicates that $R_{\rm{UPC}}({\rm A})$ and $R_{1}({\rm A})$ measurements in
UPCs may help in the experimental searches for the evidence of saturation
physics.

###### Contents

  1. I Introduction
    1. I.1 Saturation physics
    2. I.2 Double ratio in DIS
    3. I.3 A new double ratio in UPCs
    4. I.4 Outline
  2. II Ultra-Peripheral Collisions
  3. III Double ratio in UPCs coming from saturation physics
    1. III.1 Elastic vector meson production in UPCs
    2. III.2 Inclusive hadron or jet production
  4. IV Estimates of the double ratio in UPCs
    1. IV.1 Double ratio in the quasi-classical approximation
      1. IV.1.1 Elastic $J/\psi$ and $\rho$ production: heuristic estimates and numerical integration
      2. IV.1.2 Inclusive production cross section
      3. IV.1.3 $A$-scaling for the double ratio for $J/\psi$ and $\rho$ production
    2. IV.2 Double ratio with small-$x$ evolution
      1. IV.2.1 Elastic vector meson production
      2. IV.2.2 Inclusive production cross section
      3. IV.2.3 $A$-scaling for the double ratio for $J/\psi$ and $\rho$ production
  5. V Future opportunities at RHIC and LHC
  6. VI Summary and Outlook
  7. A Calculational details

##  I Introduction

###  I.1 Saturation physics

Understanding Quantum Chromodynamics (QCD) in high-energy scattering is one of
the open problems in physics. An interesting phenomenon one expects to find at
very high energy is the saturation of gluon and quark densities in the small
Bjorken $x$ region of the hadronic and nuclear wave functions [1]. (Small
values of Bjorken $x$ are probed in high-energy collisions.) The quark and
gluon parton distribution functions (PDFs) are known to grow with decreasing
$x$. However, this rise of PDFs cannot continue indefinitely: larger number of
partons leads to higher parton densities in the hadronic or nuclear wave
functions, which, in turn, lead to parton mergers becoming important, slowing
down the growth of the PDF with decreasing $x$. The regime where the parton
mergers compensate for the parton splittings is known as the parton saturation
(see Refs. [1, 2, 3, 4, 5, 6, 7, 8] for reviews). Unambiguously identifying
saturation signatures in the experimental data is one of the main challenges
of the modern hadronic physics, comprising a large part of the experimental
program for the future Electron-Ion Collider (EIC) [9, 10, 11, 12].

Theoretical description of parton saturation usually involves the linear
Balitsky–Fadin–Kuraev–Lipatov (BFKL) [13, 14] and non-linear
Balitsky–Kovchegov (BK) [15, 16, 17, 18] and Jalilian-
Marian–Iancu–McLerran–Weigert–Leonidov–Kovner (JIMWLK) [19, 20, 21, 22, 23,
24] evolution equations in $x$. While the linear BFKL evolution predicts the
rise of PDFs at small $x$, it does not contain saturation effects, and leads
to cross sections growing as a power of the center-of-mass energy, violating
unitarity at very high energies. Unitarity (in the sense of adhering to the
black disk limit) is restored when non-linear (parton merger) effects are
taken into account: this is accomplished by using the BK/JIMWLK evolution
equations. The initial conditions for the non-linear evolution equations are
given by the quasi-classical Glauber–Gribov–Mueller (GGM)/McLerran–Venugopalan
(MV) model [25, 26, 27, 28], which includes all-orders of multiple scatterings
of the projectile in the target. The above description is valid for the
scattering of a compact projectile on a large target, as seen in deep
inelastic scattering (DIS) experiments on nuclei (and, possibly, on protons)
and proton–nucleus ($p+$A) collisions at high energies. The saturation
formalism is also referred to as the color glass condensate (CGC).

The non-linear small-$x$ evolution combined with the GGM-MV initial conditions
dynamically generates a momentum scale intrinsic to the hadronic and nuclear
wave functions at small $x$ – the saturation scale $Q_{s}$. This scale
designates the transition from the low-density regime (with momentum scales
$Q>Q_{s}$) to the saturated high-density regime (for $Q\lesssim Q_{s}$). The
saturation scale grows with the atomic number $A$ of the nucleus and with
decreasing $x$, approximately as $Q_{s}^{2}\propto A^{1/3}\,(1/x)^{\lambda}$
with the power $\lambda\approx 0.3$. This means that, in high-energy
collisions (corresponding to low values of $x$) and/or for collisions
involving sufficiently large nuclei (with large $A$), the saturation scale
should get large. When it gets sufficiently above the QCD confinement scale,
$Q_{s}\gg\Lambda_{\textrm{QCD}}$, the corresponding strong coupling constant
becomes small, $\alpha_{s}(Q_{s}^{2})\ll 1$, allowing for a weak-coupling
description of the small-$x$ hadronic and nuclear wave functions and,
therefore, justifying the application of perturbative techniques to the
problem.

Tantalizing evidence for observation of the saturation regime has been
collected by accelerators around the world (see e.g. [9] for a summary of the
experimental results pertaining to saturation physics, along with the more
recent data [29, 30, 31, 32]). It is widely anticipated that the future EIC,
due to its reach in $x$ and $Q^{2}$, its unprecedented high luminosity, and a
variety of nuclear targets, would be able to seal the discovery of the
saturation phenomenon. However, with the EIC completion being projected for
the early 2030s, in the meantime, the hadronic physics community must make
every effort to study for potential saturation signals at the existing high-
energy accelerators, such as the Relativistic Heavy-Ion Collider (RHIC) and
the Large Hadron Collider (LHC). In this respect, a particularly promising
channel appears to be the ultra-peripheral collisions (UPCs) [33, 34, 35, 36,
37]. As we describe below, in the UPCs one usually collides two large nuclei
(A$+$A) or a proton and a nucleus ($p+$A) at very large impact parameters,
much larger than the sum of the radii for the scattering particles. The
interaction happens due to an almost-real photon coming from one of the nuclei
in A$+$A or from the nucleus in $p+$A collisions (the projectile) interacting
with the other nucleus or the proton (the target) by fluctuating into a quark-
antiquark ($q\bar{q}$) dipole, which then scatters on the target predominantly
via $t$-channel gluon exchanges (at high energy). Thus, the UPCs provide one
an opportunity to perform a small-$x$ DIS-type of experiment, either on the
proton (in $p+$A) or nuclear (in A$+$A) targets, except with an almost real
photon ($Q^{2}\approx 0$). The UPCs thus allow us to study the scattering
amplitude of a $q\bar{q}$ dipole on a target, searching for the signals of
saturation. The $q\bar{q}$ dipole scattering amplitude is an essential
component of saturation physics, directly entering the BK evolution equation.
The aim of this work is to propose a new saturation signal that can be studied
in the UPCs at RHIC and LHC.

###  I.2 Double ratio in DIS

One of the proposed key saturation signals at the EIC is the fraction of
diffractive DIS with respect to the total DIS cross section in electron-
nucleus ($e+$A) collisions [9]. As this ratio is found to be $\sim 13\%$ at
HERA in electron-proton ($e+p$) collisions [38], saturation calculations
predict that this ratio will be enhanced (larger than 13%) in $e+$A
collisions. To quantify this enhancement, a double ratio observable was
introduced between $e+$A and $e+p$ collisions and saturation models predict it
to be above one. The general idea behind this observable is that in the black
disk limit the elastic cross section comprises $50\%$ of the total cross
section. Indeed, for a $q\bar{q}$ dipole scattering on a target the total and
elastic scattering cross sections can be written as

$\displaystyle\sigma_{\rm tot}=2\,\int
d^{2}b_{\perp}\,N({\mathbf{r}},{\mathbf{b}},Y),$ (1a)
$\displaystyle\sigma_{\rm el}=\int
d^{2}b_{\perp}\,\left[N({\mathbf{r}},{\mathbf{b}},Y)\right]^{2},$ (1b)

where $N({\mathbf{r}},{\mathbf{b}},Y)$ is the (imaginary part of the) forward
scattering amplitude for the dipole of transverse separation $\mathbf{r}$
interacting with the target at the impact parameter $\mathbf{b}$ with the
dependence on the center-of-mass energy squared $s$ entering through the
rapidity variable $Y=\ln(s\,r_{\perp}^{2})\approx\ln(1/x)$ [17]. (In our
notation, two-dimensional transverse vectors are labeled
$\mathbf{v}=(v^{x},v^{y})$, their magnitude is $v_{\perp}=|{\mathbf{v}}|$,
while the two-dimensional integration measure is $d^{2}v_{\perp}$. For the
transverse momenta $p_{T}=|\mathbf{p}|$ and the integration measure is
$d^{2}p_{T}$.) In the black disk limit $N=1$, such that $\sigma_{\rm
el}/\sigma_{\rm tot}=1/2$ (see e.g. [7] for details). Away from the black disk
limit, where $N\ll 1$, the elastic to total cross sections ratio is small,
$\sigma_{\rm el}/\sigma_{\rm tot}\ll 1$. Saturation physics is responsible for
the approach to and the onset of the black disk limit, leading to the
anticipated increase of diffractive (elastic) cross section $\sigma_{\rm
diff}$ as a fraction of the total cross section $\sigma_{\rm tot}$, both with
increasing nuclear size ($A$) and with decreasing $x$.

The relevant double ratio is [9]

$\displaystyle R_{\rm{EIC}}=\frac{\left[\left(d\sigma_{\rm
diff}/dM_{X}^{2}\right)/\sigma_{\rm tot}\right]_{\rm
e+A}}{\left[\left(d\sigma_{\rm diff}/dM_{X}^{2}\right)/\sigma_{\rm
tot}\right]_{\rm e+p}},$ (2)

where $M_{X}^{2}$ is the invariant mass of the diffractively produced hadrons.
Saturation calculations generally predict $R_{\rm{EIC}}>1$.

Figure 1:  The ratio of diffractive over total cross-section for DIS on gold
normalized to DIS on the proton plotted for different values of $M^{2}_{X}$,
the mass squared of hadrons produced in the collisions, for models assuming
saturation and non-saturation. The statistical error bars are too small to
depict and the projected systematic uncertainty for the measurements is shown
by the orange bar. The theoretical uncertainty for the predictions of the LTS
model is shown by the grey band. Reprinted from [9].

Another reason this double ratio will likely be one of the key measurements at
the EIC is that it may provide maximum distinguishing power between the
saturation calculations [39, 40, 41, 42, 43, 44] and other theoretical
approaches, e.g., the leading-twist nuclear shadowing (LTS) model [45, 46].
Nuclear shadowing, as considered in [45, 46], is expected to cause the
fraction of diffractive over total DIS cross sections to decrease in $e+$A
with respect to the $e+p$ case, leading to a smaller than one double ratio,
$R_{\rm{EIC}}<1$. This is qualitatively different from the prediction of
saturation models. The two theoretical predictions for the double ratio to be
measured at EIC, coming from saturation physics and the leading-twist
shadowing model, are shown in Fig. 1 [9], where the double ratio
$R_{\rm{EIC}}$ is plotted versus the invariant mass of the diffractively
produced hadrons $M_{X}^{2}$ for fixed values of $x$ and $Q^{2}$, as indicated
in the plot legend.

This is a truly unique measurement that cannot be done at any facility other
than the EIC. It requires both $e+p$ and $e+$A DIS measurements at the same
center-of-mass energy (or in the same kinematics) and with a large rapidity
gap. In addition, most of the systematic uncertainty related to the detector
will cancel in this ratio.

###  I.3 A new double ratio in UPCs

In UPCs, although the same double ratio $R_{\rm{EIC}}$ as proposed above at
the EIC cannot be measured, a similar double ratio is possible. First, the
ratio between observables related to a proton and a heavy nucleus is possible,
because UPC measurements can be done in both $\rm{A+A}$ and $p+\rm{A}$
collisions. In $p+\rm{A}$ collisions, due to the large charge difference
between the proton and heavy ion, the dominant process is $\gamma+p$
scattering, where photons are emitted by the heavy ion. Second, the ratio
between diffractive DIS and total cross section cannot directly transform into
a UPC measurement. Experimentally, rapidity gap events in UPCs are very
difficult to measure and the event kinematics cannot be determined due to
limited acceptance. From a theory perspective, all UPC processes are in the
photoproduction limit ($Q^{2}\approx 0$) where no hard scale is available for
a reliable calculation, unless, at small $x$ and for large nuclei, the
saturation scale $Q_{s}$ becomes large enough.

While it is hard to replicate the ratio between diffractive DIS and total DIS
cross section in UPCs, we propose to replace it with a ratio between the
diffractive vector-meson (VM) photoproduction and inclusive particle or jet
photoproduction. The idea is similar to the above: as we will show below, the
elastic vector meson production cross section is quadratic in the dipole
amplitude $N$, while the inclusive particle production is linear in $N$ [47],
such that their ratio should be sensitive to the approach to the black disk
limit. Here the hard scale of the diffractive probe comes from the mass of the
VM, e.g., the $J/\psi$ meson, and from the transverse momentum of the jet or
charged particle in inclusive photoproduction, in addition to the saturation
scale $Q_{s}$. The new double ratio is defined as

$\displaystyle
R_{\rm{UPC}}=\rm{\frac{\left[\sigma^{VM}_{el}/\left(d\sigma^{jet}_{inclusive}/d^{2}p_{T}\right)\right]_{A+A}}{\left[\sigma^{VM}_{el}/\left(d\sigma^{jet}_{inclusive}/d^{2}p_{T}\right)\right]_{p+A}}}.$
(3)

Indeed, the ratio $R_{\rm{UPC}}$ is not quite the same as $R_{\rm{EIC}}$. In
particular, one may argue that $\sigma_{\rm inclusive}=\langle
n\rangle\,\sigma_{\rm inel}$ with $\langle n\rangle$ the average multiplicity
of the produced hadrons and $\sigma_{\rm inel}$ the total inelastic cross
section. In a (perhaps hypothetical) case when both the nucleus and the proton
are scattering near the black disk limit, where $\sigma_{\rm
inel}\approx(1/2)\,\sigma_{\rm tot}$ and $R_{\rm{EIC}}\approx 1$, one can
approximate $R_{\rm{UPC}}\approx R_{\rm{EIC}}\,\langle n\rangle_{\gamma+{\rm
p}}/\langle n\rangle_{\gamma+{\rm A}}\approx\langle n\rangle_{\gamma+{\rm
p}}/\langle n\rangle_{\gamma+{\rm A}}$. Since the particle multiplicities are
different in $\gamma+p$ and $\gamma+$A collisions, $\langle
n\rangle_{\gamma+{\rm p}}/\langle n\rangle_{\gamma+{\rm A}}<1$, we see that
$R_{\rm{UPC}}<1$ in this example. Therefore, we cannot expect that
$R_{\rm{UPC}}$ is necessarily always greater than one in the saturation
picture, making this ratio different from $R_{\rm{EIC}}$.

However, what we will show below is that the $A$-scaling of $R_{\rm{UPC}}$ is
different inside and outside the saturation region. That is, we will
demonstrate that the $A$-dependence of $R_{\rm{UPC}}$ changes depending on
whether the produced vector meson is relatively small ($J/\psi$) or large
($\rho$), and also depending on the transverse momentum $p_{T}$ of the
produced hadron or jet in the inclusive cross section being larger than or
smaller than the saturation scale $Q_{s}$. We, therefore, expect that the
$A$-dependence of $R_{\rm{UPC}}$ would be of great interest in UPCs,
especially since the center-of-mass energy reach can be much higher in the
UPCs than that at the EIC. In fact, $R_{\rm{UPC}}$ measurements may lay the
groundwork for the future $R_{\rm{EIC}}$ at the EIC, possibly providing
necessary complimentary information needed to seal the saturation discovery at
the EIC.

###  I.4 Outline

The paper is structured as follows. After giving a brief introduction to UPCs
in Sec. II, we will re-derive the leading-order (LO) expressions for the
elastic vector meson production and the total hadron/jet production cross
section in the saturation/small-$x$ formalism in Sections III.1 and III.2,
respectively. While a next-to-leading order (NLO) calculation for elastic
vector meson production exists [48], and will need to be utilized in the
future more detailed phenomenological studies, our goal is to explore the
double ratio $R_{\rm{UPC}}$ at the more qualitative level, trying to identify
the difference between its $A$-dependence inside and outside the saturation
region. For such a goal, the LO expression for the cross sections should be
sufficient.

In Sec. IV.1 we evaluate the double ratio $R_{\rm{UPC}}$ using the GGM-MV
approximation for the dipole amplitude $N$, obtaining different $A$-dependence
for $R_{\rm{UPC}}$ inside and outside the saturation region for the $\rho$ and
$J/\psi$ meson production (in the numerator of the ratio). We repeat the
calculation in Sec. IV.2 for the dipole amplitude given by the geometric
scaling solution of the LO BK equation [49, 50, 1, 51, 52], reaching similar
conclusions as in the quasi-classical case. Based on our calculations, we
discuss the future opportunities to measure $R_{\rm{UPC}}$ in Sec. V, before
summarizing our results in Sec. VI.

##  II Ultra-Peripheral Collisions

Ultra-peripheral collisions (UPCs) represent a fascinating experimental
opportunity at heavy-ion collider facilities, offering a captivating avenue
for scientific exploration of the fundamental structure of nucleon and nuclei.
In these collisions, heavy ions such as gold (Au) or lead (Pb) approach each
other with impact parameters far exceeding the size of the ions themselves,
e.g., $b>2R_{A}$. See Fig. 2 for an illustration. This results in an
electromagnetic interaction dominated by the exchange of virtual photons,
where the photon probe interacts with the other nucleus. Therefore, the UPCs
are similar to the electron-ion scattering experiments with quasi-real photon
beams. Here we list a few of the important features of UPCs:

Figure 2:  Illustration of UPC Vector Meson (VM) photoproduction (left) and
jet photoproduction (right) in heavy-ion A$+$A collisions.

  * •

Virtuality. Due to the heavy nuclei, photons are emitted with small transverse
momentum that has a typical value of $\sim 20-40~{}\rm{MeV/c}$; the
virtuality, $Q^{2}$, is around $10^{-4}-10^{-3}~{}\rm{GeV^{2}}$. The hard
scale of the process is usually determined by the final-state particles, e.g.,
their mass or transverse momentum.

  * •

Rates. For QED processes, e.g., two-photon scatterings cross section, the rate
of the A+A reaction scales with the nuclear electric charge as $\sim Z^{4}$.
(Here $Z$ is the charge of the nucleus, in the units of electron charge
$|e|$.) For QCD processes, e.g., photon-Pomeron interactions, the rate scales
as $\sim Z^{2}$.

  * •

Energy reach. The center-of-mass energy between the photon and per nucleon in
the nucleus, $W_{\gamma*N}$, ranges from less than 10 GeV (RHIC) to $\sim
1~{}\rm{TeV}$ (LHC). This corresponds to momentum fraction $x$ in the range
from $10^{-1}$ to $10^{-5}$. Note that, since the saturation scale grows as
$Q^{2}_{s}\sim(A/x)^{1/3}$, the UPCs at LHC can probe larger $Q^{2}_{s}$
values than those to be probed in the collisions at the future EIC.

  * •

Observables. Most UPCs measurements are for the exclusive and diffractive
Vector-Meson (VM) photoproduction. Event kinematics and topology are easy to
reconstruct. Recent measurements show that semi-inclusive processes, e.g.,
inclusive and diffractive jets [53], and inclusive processes of measuring
charged particles only [54], are feasible.

There are two major heavy-ion collider facilities currently studying UPCs,
RHIC and the LHC. With the heavy-ion runs scheduled between 2023 and 2025 both
at RHIC and the LHC, there are a couple of experimental opportunities to the
UPCs. For example, forward detectors have been added to the STAR experiment,
where particle acceptance between 2.5 to 4 in pseudorapidity $\eta$ became
available. Typically, the forward rapidity corresponds to high-$x$ kinematic
region. In order to unambiguously discover saturation, the high-$x$ region
where we do not expect saturation should also be studied on the same footing
as the low-$x$ region. In addition, exploration of the large-$x$ region may
shine new light to the anti-shadowing region of nuclear PDFs. Currently,
resolving the photon energy ambiguity using the forward neutron emissions has
been shown possible, extending the kinematic reach at low-$x$ by two orders of
magnitude at the LHC.

##  III Double ratio in UPCs coming from saturation physics

###  III.1 Elastic vector meson production in UPCs

As described above, UPCs (in A$+$A or $p+$A) involve the interaction of a
photon of a very small virtuality originating from one nucleus with the other
nucleus or proton. At high energy (small $x$), the dominating process involves
the photon fluctuating into a $q\bar{q}$ pair that interacts with the nucleus.
One of the possible outcomes of this interaction is the exclusive production
of vector meson. The scattering amplitude for exclusive vector meson
production in UPCs at leading order is given by the diagram shown in Fig. 3,

Figure 3: Exclusive vector meson production amplitude in DIS and UPCs at high
energy.

where the oval on the upper right represents the vector meson light-cone wave
function [55], the ovals on the lower left and lower right represent the
target proton or nucleus which remains intact in the interaction, and the
vertical rectangle denotes the (mainly gluon) shock wave responsible for the
interaction. For the target to remain intact, the shock wave should transfer
no color away from the target, i.e., the exchange needs to be color-singlet.
To obtain the cross section, we need to square the scattering amplitude and
integrate over the phase space of the final state particles. For the
convenience of our later calculation, we should Fourier transform the
amplitude into transverse coordinate space. The result for the elastic vector
meson production cross section at leading order is [56, 57] (see [7] for a
pedagogical derivation in the transverse coordinate space)

$\sigma^{\gamma^{*}{\rm A}\to V{\rm
A}}_{\textrm{el}}=\frac{1}{2}\sum_{\lambda}\sum_{\lambda^{\prime}}\int
d^{2}b_{\perp}\left|\int\frac{d^{2}r_{\perp}}{4\pi}\int\limits_{0}^{1}\frac{dz}{z(1-z)}\Psi_{T}^{\gamma^{*}\to
q\bar{q}}(\mathbf{r},z)\,N(\mathbf{r},\mathbf{b},Y)\,\Psi^{V}_{T}(\mathbf{r},z)^{*}\right|^{2}.$
(4)

Here $\Psi_{T}^{\gamma^{*}\to q\bar{q}}(\mathbf{r},z)$ is the light-cone wave
function for a transversely polarized photon splitting into a $q\bar{q}$ pair
of the transverse separation $\mathbf{r}$ with the quark carrying a fraction
$z$ of the photon’s longitudinal momentum, $\Psi_{T}^{V}$ is the transversely
polarized vector meson wave function, and $N(\mathbf{r},\mathbf{b},Y)$ is the
dipole amplitude for a $q\bar{q}$ pair interacting with the nucleus at impact
parameter $\mathbf{b}$ and rapidity $Y=\ln(s\,r_{\perp}^{2})\approx\ln(1/x)$
introduced above. We use the definition of the light-cone wave functions
normalized as in [7]. As shown in Fig. 3, the photon comes in with
polarization $\lambda$, while the vector meson is produced with the
polarization $\lambda^{\prime}$. The former is averaged over, assuming that
the photon in UPCs is real, due to its very low $Q^{2}$: therefore, we keep
only the transverse polarizations for the photon. This, in turn, implies that
the outgoing vector meson has to be transversely polarized as well. We sum
over the vector meson’s transverse polarizations in Eq. (4). In addition, a
sum over quark polarizations $\sigma,\sigma^{\prime}$ (as denoted in Fig. 3)
and quark colors is implied inside the absolute value brackets in Eq. (4).

In the saturation/CGC framework, the dipole amplitude $N$ is found by using
the BK/JIMWLK nonlinear evolution equations [15, 16, 17, 18, 19, 20, 21, 22,
23, 24] with the GGM/MV quasi-classical initial conditions [25, 26, 27, 28].
Below, we will consider the quasi-classical and small-$x$ evolution regimes
separately, analysing $R_{\rm UPC}$ from Eq. (3) in each case. To proceed with
the calculation, we need to delve into a discussion about the photon and
vector meson wave functions.

The virtual photon wave function at the leading order is well known [58, 59].
The corresponding diagram for a photon splitting into a quark-antiquark pair
is shown below in Fig. 4. In the case of UPCs, the photon has a very small
$Q^{2}$, so the longitudinally polarized photon wave function is small
compared to the transversely polarized one: as already mentioned above, we
will ignore the former in our calculation. The transverse coordinate space
transversely-polarized photon wave function that appears in Eq. (4) is

$\Psi_{T}^{\gamma^{*}\to
q\bar{q}}(\mathbf{r},z)=\frac{eZ_{f}}{2\pi}\sqrt{z(1-z)}\,\delta_{ij}\left[-i\,\delta_{\sigma,-\sigma^{\prime}}\,(1-2z-\sigma\lambda)\,\frac{\mathbf{\epsilon}^{\lambda}\cdot\mathbf{r}}{r_{\perp}}\,\partial_{r_{\perp}}+\delta_{\sigma\sigma^{\prime}}\,\frac{m_{f}}{\sqrt{2}}\,(1+\sigma\lambda)\right]K_{0}(r_{\perp}a_{f}),$
(5)

where $\sigma,\sigma^{\prime}$ are the polarizations and $i,j$ are the colors
of the quark and antiquark, respectively. We have defined $a_{f}^{2}\equiv
z(1-z)\,Q^{2}+m_{f}^{2}$ with $m_{f}$ the mass of quarks with flavor $f$. The
magnitude of the electron charge is labelled by $e$, while $Z_{f}$ is the
fractional charge of the quark. In addition,
$\partial_{r_{\perp}}=\partial/\partial r_{\perp}$. As above, the photon has
polarization $\lambda$ and $z$ is the longitudinal fraction of photon’s
momentum carried by the quark.

Figure 4: Virtual photon wave function at leading order.

The vector meson wave function, on the other hand, is more complicated because
of its nonperturbative nature. We could, however, obtain a reasonable form of
the vector meson function if we are willing to make a few assumptions. For our
purposes we will follow the approach put forth in [60, 61, 62, 63, 64, 65].
The starting point of this approach is the assumption that the vector meson is
predominantly a quark–antiquark state that has the same polarization structure
as a virtual photon. By making the replacement

$\frac{eZ_{f}}{2\pi}z(1-z)K_{0}(r_{\perp}a_{f})\rightarrow\phi_{T}(r_{\perp},z),$
(6)

for the scalar part of the wave function, we can write down the transversely
polarized vector meson wave function as

$\Psi_{T}^{V}(\mathbf{r},z)=\frac{1}{\sqrt{z(1-z)}}\delta_{ij}\left[-i\,\delta_{\sigma,-\sigma^{\prime}}\,(1-2z-\sigma\lambda)\frac{\mathbf{\epsilon}^{\lambda}\cdot\mathbf{r}}{r_{\perp}}\partial_{r_{\perp}}+\delta_{\sigma\sigma^{\prime}}\frac{m_{f}}{\sqrt{2}}(1+\sigma\lambda)\right]\phi_{T}(r_{\perp},z).$
(7)

Note again that transversely polarized photon produces only transversely
polarized vector meson; therefore, in UPCs transverse vector mesons production
is dominant (with the longitudinal vector meson production being suppressed by
a factor of $Q^{2}/Q_{s}^{2}\ll 1$). We will employ only the transversely
polarized vector meson wave function.

The scalar part of the vector meson wave function can be reasonably described
using a “boosted Gaussian” form [66, 61, 62, 64]. We start by assuming that in
the rest frame of the vector meson, the wave function is Gaussian in the
transverse momentum of the quark and antiquark, so that

$\tilde{\phi}_{T}(p_{T},z)\propto\exp(-\frac{R^{2}}{2}p_{T}^{2}),$ (8)

where $p_{T}=|\mathbf{p}|$ is the transverse momentum of the quark (or
antiquark) and $R$ is the typical radius of a vector meson. The momentum
$\mathbf{p}$ is related to the invariant mass through

$M^{2}=4E^{2}=4(p_{T}^{2}+m_{f}^{2}).$ (9)

In the boosted frame, we have a fast moving vector meson with its quark and
antiquark having longitudinal momentum fractions $z$ and $1-z$, respectively.
The four momenta of the on-mass-shell quark and antiquark are

$\displaystyle p_{1}^{\prime}$
$\displaystyle=\left(zP^{+},\frac{k_{T}^{2}+m_{f}^{2}}{zP^{+}},\mathbf{k}\right),$
$\displaystyle p_{2}^{\prime}$
$\displaystyle=\left((1-z)P^{+},\frac{k_{T}^{2}+m_{f}^{2}}{(1-z)P^{+}},-\mathbf{k}\right),$
(10)

where $\mathbf{k}$ is the transverse momentum of the quark (or antiquark) in
the boosted frame and we assume that the meson is moving in the light-cone
plus direction with a large momentum $P^{+}$. The invariant mass is

$M^{2}=(p^{\prime}_{1}+p^{\prime}_{2})^{2}=\frac{k_{T}^{2}+m_{f}^{2}}{z(1-z)}.$
(11)

Equating the invariant masses in two frames, Eqs. (9) and (11), we can relate
the transverse momenta $p_{T}$ and $k_{T}$,

$p_{T}^{2}=\frac{k_{T}^{2}+m_{f}^{2}}{4z(1-z)}-m_{f}^{2}.$ (12)

Further, if we assume that the wave function is boost-invariant (i.e., the
scalar part of the wave function takes on the same value after a Lorentz
transformation), then in the boosted frame the wave function has the form (cf.
Eq. (8))

$\tilde{\phi}_{T}(k_{T},z)\propto\exp\left[-\frac{R^{2}}{8}\left(\frac{k_{T}^{2}+m_{f}^{2}}{z(1-z)}-4m_{f}^{2}\right)\right].$
(13)

Fourier-transforming the wave function into the transverse coordinate space,
we have

$\phi_{T}(r_{\perp},z)\propto
z(1-z)\exp\left[-\frac{2z(1-z)\,r_{\perp}^{2}}{R^{2}}-\frac{m_{f}^{2}R^{2}}{8z(1-z)}\right],$
(14)

where we have absorbed the $\frac{m_{f}^{2}R^{2}}{2}$ term in the exponent
into the overall normalization, since it is independent of $r_{\perp}$ and
$z$. Substituting the result (14) into Eq. (7), along with the normalization
factor ${\cal N}_{T}$, we obtain an expression for the transversely polarized
vector meson wave function [60, 61, 62, 63, 64],

$\displaystyle\Psi_{T}^{V}(r_{\perp},z)={\cal
N}_{T}\sqrt{z(1-z)}\,\delta_{ij}\left[i\,\delta_{\sigma,-\sigma^{\prime}}\,(1-2z-\sigma\lambda)\frac{4z(1-z)}{R^{2}}\,\mathbf{\epsilon}^{\lambda}\cdot\mathbf{r}\right.$
(15)
$\displaystyle+\left.\delta_{\sigma\sigma^{\prime}}\frac{m_{f}}{\sqrt{2}}(1+\sigma\lambda)\right]\exp\left[-\frac{2z(1-z)\,r_{\perp}^{2}}{R^{2}}-\frac{m_{f}^{2}R^{2}}{8z(1-z)}\right].$

The wave function (15) is normalized such that

$\displaystyle\sum_{\sigma,\sigma^{\prime},i,j}\,\int\limits_{0}^{1}\frac{dz}{z(1-z)}\,\int\frac{d^{2}r_{\perp}}{4\pi}\,\left|\Psi_{T}^{V}(r_{\perp},z)\right|^{2}=1,$
(16)

resulting from the assumption that the vector meson is predominantly a
quark–antiquark pair (in all frames). The relation (16) can be used to find
${\cal N}_{T}$: we will simply use the values in Table I of [65] to read off
${\cal N}_{T}$ and the vector meson radius $R$. The resulting parameters for
the wave function Eq. (15) for the $\rho$ and $J/\psi$ mesons are listed in
Table 1 below.

Before we proceed, we note here that vector meson light-cone wave functions
are, indeed, model-dependent. It would therefore be hard to estimate the
accuracy of a specific numerical prediction based on those wave functions.
However, our goal in this work is more qualitative than quantitative. For our
purposes, we need the $J/\psi$ wave function to be dominated by smaller
$q\bar{q}$ pairs, while the $\rho$-meson wave function needs to favor larger
dipole sizes. The wave function (15) certainly satisfies these properties.

Figure 5: A diagrammatic representation of Eq. (17). Vertical rectangles
denote the shock waves in the amplitude and in the complex conjugate
amplitude, while the solid vertical line denotes the final state cut.

For future calculations, it is easier to simplify Eq. (4) a little.
Substituting Eqs. (5) and (15) into Eq. (4), summing over polarizations
$\lambda,\lambda^{\prime},\sigma,\sigma^{\prime}$ and quark colors, and
integrating over the angles of $\mathbf{r}$ while assuming the dipole
amplitude to be independent of the $\mathbf{r}$ direction, we arrive at

$\displaystyle\sigma^{\gamma^{*}{\rm A}\to V{\rm A}}_{\textrm{el}}=$
$\displaystyle\,\frac{\alpha_{EM}Z_{f}^{2}}{\pi}\,N_{c}^{2}\,{\cal
N}_{T}^{2}\int\limits_{0}^{1}dz\,dz^{\prime}\int\limits_{0}^{\infty}dr_{\perp}\,dr^{\prime}_{\perp}\,r_{\perp}\,r_{\perp}^{\prime}$
(17)
$\displaystyle\times\left\\{\frac{4}{R^{2}}z(1-z)\left[z^{2}+(1-z)^{2}\right]\,r_{\perp}\,a_{f}\,K_{1}(r_{\perp}a_{f})+m_{f}^{2}\,K_{0}(r_{\perp}a_{f})\right\\}\,\exp\left[-\frac{2z(1-z)\,r_{\perp}^{2}}{R^{2}}-\frac{m_{f}^{2}R^{2}}{8z(1-z)}\right]$
$\displaystyle\times\left\\{\frac{4}{R^{2}}z^{\prime}(1-z^{\prime})\left[z^{\prime\,2}+(1-z^{\prime})^{2}\right]\,r_{\perp}^{\prime}\,a_{f}^{\prime}\,K_{1}(r_{\perp}^{\prime}a_{f}^{\prime})+m_{f}^{2}K_{0}(r_{\perp}^{\prime}a_{f}^{\prime})\right\\}\exp\left[-\frac{2z^{\prime}(1-z^{\prime})r_{\perp}^{\prime
2}}{R^{2}}-\frac{m_{f}^{2}R^{2}}{8z^{\prime}(1-z^{\prime})}\right]$
$\displaystyle\times\int
d^{2}b_{\perp}\,N(r_{\perp},b_{\perp},Y)\,N(r^{\prime}_{\perp},b_{\perp},Y).$

Here $a_{f}^{\prime\,2}\equiv z^{\prime}(1-z^{\prime})\,Q^{2}+m_{f}^{2}$ and
$\alpha_{EM}=e^{2}/(4\pi)$. The assumption of
$N(\mathbf{r},\mathbf{b},Y)=N(r_{\perp},b_{\perp},Y)$ that we employ here is
valid in the quasi-classical approximation for a large nucleus [25, 26, 27,
28] and for the leading contribution to high-energy scattering (see, e.g.,
[7]). Diagrammatically, Eq. (17) can be represented by the graph in Fig. 5,
where we separately integrate over the transverse sizes
$r_{\perp},r_{\perp}^{\prime}$ of the dipole in the amplitude and in the
complex conjugate amplitude, and over the impact parameter $\mathbf{b}$,
common to both amplitudes.

###  III.2 Inclusive hadron or jet production

The ratio $R_{\rm{UPC}}$ in Eq. (3) also contains the inclusive hadron or jet
production cross section. Our goal here is to assess the behavior of
$R_{\rm{UPC}}$ with varying nuclear atomic number $A$. We will assume that at
small-$x$ fragmentation and jet showering happen mainly outside the cold
nuclear medium of the target nucleus and are, therefore, $A$-independent.
Hence, we will omit such final-state interactions from our analysis and
concentrate on calculating the quark production cross section, which we will
then use to derive the $A$-dependence of the double ratio $R_{\rm{UPC}}$.

Figure 6: A diagram contributing to the scattering amplitude for the inclusive
single quark production in the virtual photon–target proton (nucleus)
scattering. The measured quark is marked by a cross.

To calculate the inclusive quark production cross section, we again consider
the dipole picture of the photon–nucleus scattering, which is the dominant
contribution at high energy (small $x$). A typical contributing diagram is
shown in Fig. 6, constructed by analogy to Fig. 3. The difference now is that
the target proton or nucleus can break up, which means the interaction is no
longer color-singlet. Moreover, the $q\bar{q}$ pair resulting from the virtual
photon’s splitting does not recombine into a vector meson. Both the quark and
antiquark in the pair can fragments into independent hadrons or jets. We will
calculate the quark production cross section (with the produced quark denoted
by the cross in Fig. 6): in the eikonal approximation employed here, the
antiquark production cross section is exactly the same.

The scattering amplitude in Fig. 6 can be factorized into a convolution of the
photon light cone wave function with the dipole amplitude for the $q\bar{q}$
pair scattering on the nucleus. In a light-cone perturbation theory
calculation [55], to the diagram in Fig. 6 one has to add the non-interacting
contribution where the $\gamma\to q\bar{q}$ splitting takes place to the right
of the shock wave. The result for the differential cross section
$\frac{d\sigma}{d^{2}p_{T}}$ was first derived in [67] and is also given by
Eq. (85) of [68] in terms of the dipole $S$-matrix. Recalling that this
$S$-matrix is $S=1-N$ for the dipole amplitude $N$, we write

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}=\frac{1}{2(2\pi)^{3}}$
$\displaystyle\ \int\limits_{0}^{1}\frac{dz}{z(1-z)}\,\int
d^{2}x_{\perp}d^{2}x_{\perp}^{\prime}d^{2}y_{\perp}\,e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\,\frac{1}{2}\,\sum_{\lambda}\Psi^{\gamma\to
q\bar{q}}_{T}(\mathbf{x}-\mathbf{y},z)\,\left[\Psi_{T}^{\gamma\to
q\bar{q}}(\mathbf{x}^{\prime}-\mathbf{y},z)\right]^{*}$ (18)
$\displaystyle\times\,\left[N\left(\mathbf{x}-\mathbf{y},\frac{\mathbf{x}+\mathbf{y}}{2},Y\right)+N\left(\mathbf{x}^{\prime}-\mathbf{y},\frac{\mathbf{x}^{\prime}+\mathbf{y}}{2},Y\right)-N\left(\mathbf{x}-\mathbf{x}^{\prime},\frac{\mathbf{x}+\mathbf{x}^{\prime}}{2},Y\right)\right],$

where the sum over the final-state quark polarizations and colors is implied
in the product of light-cone wave functions. As we have mentioned earlier, in
UPCs the virtual photon has a very small virtuality ($Q^{2}\approx
400-900~{}\text{MeV}^{2}$). Therefore we only need to consider transversely
polarized photons, with Eq. (18) including the averaging over the virtual
photon polarizations. Also, note that we are continuing to use the
$N(\mathbf{r},\mathbf{b},Y)$ notation for the dipole amplitudes, resulting in
a slightly more complicated than usual arguments of the dipole amplitudes in
Eq. (18).

The expression (18) can also be further simplified. Substituting the virtual
photon light-cone wave function from Eq. (5) into Eq. (18), averaging over
photon polarizations and summing over final-state quark and antiquark
polarizations, we obtain

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}=\frac{2N_{c}\,\alpha_{EM}Z_{f}^{2}}{(2\pi)^{4}}\int
d^{2}x_{\perp}d^{2}x_{\perp}^{\prime}d^{2}y_{\perp}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\int\limits_{0}^{1}dz$
(19)
$\displaystyle\times\,\left\\{\left[z^{2}+(1-z)^{2}\right]a_{f}^{2}\frac{(\mathbf{x}-\mathbf{y})\cdot(\mathbf{x}^{\prime}-\mathbf{y})}{|\mathbf{x}-\mathbf{y}||\mathbf{x}^{\prime}-\mathbf{y}|}K_{1}(|\mathbf{x}-\mathbf{y}|a_{f})\,K_{1}(|\mathbf{x}^{\prime}-\mathbf{y}|a_{f})+m_{f}^{2}K_{0}(|\mathbf{x}-\mathbf{y}|a_{f})K_{0}(|\mathbf{x}^{\prime}-\mathbf{y}|a_{f})\right\\}$
$\displaystyle\times\,\left[N\left(\mathbf{x}-\mathbf{y},\frac{\mathbf{x}+\mathbf{y}}{2},Y\right)+N\left(\mathbf{x}^{\prime}-\mathbf{y},\frac{\mathbf{x}^{\prime}+\mathbf{y}}{2},Y\right)-N\left(\mathbf{x}-\mathbf{x}^{\prime},\frac{\mathbf{x}+\mathbf{x}^{\prime}}{2},Y\right)\right].$

Some of the integrals in Eq. (19) can be carried out explicitly. This is done
in Appendix A, resulting in

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}=\frac{N_{c}\,\alpha_{EM}Z_{f}^{2}}{(2\pi)^{3}}\int\limits_{0}^{1}dz\,\int
d^{2}r_{\perp}d^{2}b_{\perp}e^{-i\mathbf{p}\cdot\mathbf{r}}\left\\{[z^{2}+(1-z)^{2}]\left[4ia_{f}K_{1}(r_{\perp}a_{f})\frac{\mathbf{p}\cdot\mathbf{r}}{(p_{\perp}^{2}+a_{f}^{2})r_{\perp}}-2\,K_{0}(r_{\perp}a_{f})\right.\right.$
(20)
$\displaystyle\left.+\,r_{\perp}\,a_{f}\,K_{1}(r_{\perp}a_{f})\Bigg{]}+\frac{4m_{f}^{2}}{p_{\perp}^{2}+a_{f}^{2}}K_{0}(r_{\perp}a_{f})-\frac{m_{f}^{2}r_{\perp}}{a_{f}}K_{1}(r_{\perp}a_{f})\right\\}N(\mathbf{r},\mathbf{b},Y).$

In arriving at Eq. (20), we have employed the fact that

$\displaystyle\int d^{2}b_{\perp}\,N(\mathbf{r},\mathbf{b},Y)=\int
d^{2}b_{\perp}\,N(-\mathbf{r},\mathbf{b},Y),$ (21)

which is valid due to $N$ being a scalar quantity and our unpolarized proton
or nuclear target having no preferred transverse direction.

##  IV Estimates of the double ratio in UPCs

Our next goal is to study the dependence of the double ratio $R_{\rm{UPC}}$
from Eq. (3) on the atomic number $A$ for different vector mesons. We will
pick the $\rho$ meson as an example of a vector meson whose wave function is
dominated by “large” dipole sizes, and $J/\psi$ meson as an example of a
small-dipole-size meson. Below we evaluate the $A$-dependence of
$R_{\rm{UPC}}$ for these two mesons analytically and numerically. We will
first work in the quasi-classical GGM/MV approximation [25, 26, 27, 28], then
we will include LO small-$x$ BK/JIMWLK evolution [15, 16, 17, 18, 19, 20, 21,
22, 23, 24].

###  IV.1 Double ratio in the quasi-classical approximation

####  IV.1.1 Elastic $J/\psi$ and $\rho$ production: heuristic estimates and
numerical integration

In Eq. (4) (or Eq. (17)), the dipole amplitude $N(\mathbf{r},b,Y)$ describes
the interaction of the quark-antiquark pair with the target nucleus. In
processes where saturation effects are taken into account, one has to include
multiple gluon exchanges between the $q\bar{q}$ pair and the nucleus.
Including $t$-channel gluon exchanges to all orders in the GGM/MV model leads
to the dipole amplitude [25]

$N(\mathbf{r},\mathbf{b},Y)=1-\exp\left\\{-\frac{r_{\perp}^{2}Q_{s}^{2}(\mathbf{b})}{4}\ln\frac{1}{r_{\perp}\Lambda}\right\\}$
(22)

with the saturation scale $Q_{s}$ given by

$\displaystyle
Q_{s}^{2}(\mathbf{b})=4\pi\alpha_{s}^{2}\,\frac{C_{F}}{N_{c}}\,T(\mathbf{b}).$
(23)

Here $T(\mathbf{b})$ is the nuclear profile (thickness) function,

$\displaystyle
T(\mathbf{b})\equiv\int\limits_{-\infty}^{\infty}dz\,\rho(\mathbf{b},z),$ (24)

where $\rho(\mathbf{b},z)$ is the nucleon number density in the nucleus,
$\Lambda$ is the infrared (IR) cutoff, and $C_{F}$ is the fundamental Casimir
operator of SU($N_{c}$). We should point out that the dipole amplitude given
in Eq. (22) does not include small-$x$ evolution: this is why
$Q_{s}^{2}(\mathbf{b})$ here is independent of energy/rapidity $Y$, leading to
similarly energy-independent dipole amplitude $N(\mathbf{r},\mathbf{b},Y)$ in
Eq. (22).

Our goal now is to determine the dependence of the elastic VM production cross
section on the atomic number $A$. After a closer inspection of Eq. (17), we
see that the $A$-dependence is contained entirely in the $\mathbf{b}$-integral

$\displaystyle\int
d^{2}b_{\perp}\,N(r_{\perp},b_{\perp},Y)\,N(r^{\prime}_{\perp},b_{\perp},Y)$
(25)

over the transverse area of the nucleus. This integral is hard to evaluate
exactly analytically. Therefore, we have to make approximations for the dipole
amplitude $N(\mathbf{r},\mathbf{b},Y)$ based on whether $r_{\perp}$ and
$r^{\prime}_{\perp}$ are larger or smaller than $1/Q_{s}(\mathbf{b})$, which
corresponds to the dipole $r_{\perp}$ and/or the dipole $r^{\prime}_{\perp}$
being inside or outside the saturation regime (see Fig. 5). Since the
integrations over $r_{\perp}$ and $r_{\perp}^{\prime}$ range over all positive
values between $0$ and $\infty$, we have three cases to consider: (i)
$r_{\perp},r_{\perp}^{\prime}\ll 1/Q_{s}$, (ii)
$r_{\perp},r_{\perp}^{\prime}\gtrsim 1/Q_{s}$, and (iii) $r_{\perp}\ll
1/Q_{s},r_{\perp}^{\prime}\gtrsim 1/Q_{s}$. The case when
$r^{\prime}_{\perp}\ll 1/Q_{s},r_{\perp}\gtrsim 1/Q_{s}$ gives the same
contribution as the case (iii), due to the
$\mathbf{r}\leftrightarrow\mathbf{r}^{\prime}$ symmetry of Eq. (17). As
follows from Eq. (17), the dipole sized $r_{\perp}$ and $r^{\prime}_{\perp}$
are controlled by the convolutions of the virtual photon and vector meson wave
functions with the dipole size dependence of the amplitude $N$.

In these three regions we obtain different $A$-scaling, using the following
arguments:

  1. (i)

$r_{\perp},r_{\perp}^{\prime}\ll 1/Q_{s}$: We approximate the dipole amplitude
(22) outside the saturation region by expanding it to the lowest order in
$r_{\perp}\,Q_{s}(\mathbf{b})$, such that

$\displaystyle
N(\mathbf{r},\mathbf{b},Y)\bigg{|}_{r_{\perp}\,Q_{s}(\mathbf{b})\ll
1}\approx\frac{r_{\perp}^{2}Q_{s}^{2}(\mathbf{b})}{4}\ln\frac{1}{r_{\perp}\Lambda}\propto
A^{1/3},$ (26a) $\displaystyle
N(\mathbf{r}^{\prime},\mathbf{b},Y)\bigg{|}_{r^{\prime}_{\perp}\,Q_{s}(\mathbf{b})\ll
1}\approx\frac{r_{\perp}^{\prime
2}Q_{s}^{2}(\mathbf{b})}{4}\ln\frac{1}{r_{\perp}^{\prime}\Lambda}\propto
A^{1/3},$ (26b)

where the last proportionality follows from $Q_{s}^{2}(\mathbf{b})\propto
T(\mathbf{b})\propto A^{1/3}$. Since the area integral scales as $\int
d^{2}b_{\perp}\sim A^{2/3}$, we conclude that

$\int
d^{2}b_{\perp}\,N(\mathbf{r},\mathbf{b},Y)\,N(\mathbf{r}^{\prime},\mathbf{b},Y)\bigg{|}_{r_{\perp},r_{\perp}^{\prime}\ll
1/Q_{s}}\propto A^{4/3}.$ (27)

  2. (ii)

$r_{\perp},r_{\perp}^{\prime}\gtrsim 1/Q_{s}$: Inside the saturation region we
approximate

$N(\mathbf{r},\mathbf{b},Y)\bigg{|}_{r_{\perp}\,Q_{s}(\mathbf{b})\gtrsim
1}\approx
N(\mathbf{r}^{\prime},\mathbf{b},Y)\bigg{|}_{r^{\prime}_{\perp}\,Q_{s}(\mathbf{b})\gtrsim
1}\approx 1,$ (28)

such that

$\int
d^{2}b_{\perp}\,N(\mathbf{r},\mathbf{b},Y)\,N(\mathbf{r}^{\prime},\mathbf{b},Y)\bigg{|}_{r_{\perp},r_{\perp}^{\prime}\gtrsim
1/Q_{s}}\propto A^{2/3}.$ (29)

  3. (iii)

$r_{\perp}\ll 1/Q_{s},r_{\perp}^{\prime}\gtrsim 1/Q_{s}$ (or
$r^{\prime}_{\perp}\ll 1/Q_{s},r_{\perp}\gtrsim 1/Q_{s}$): With one dipole
being outside the saturation region, and another one being inside, we have

$\displaystyle
N(\mathbf{r},\mathbf{b},Y)\bigg{|}_{r_{\perp}\,Q_{s}(\mathbf{b})\ll
1}\approx\frac{r_{\perp}^{2}Q_{s}^{2}}{4}\ln\frac{1}{r_{\perp}\Lambda},$ (30a)
$\displaystyle
N(\mathbf{r}^{\prime},\mathbf{b},Y)\bigg{|}_{r^{\prime}_{\perp}\,Q_{s}(\mathbf{b})\gtrsim
1}\approx 1.$ (30b)

This leads to

$\int
d^{2}b_{\perp}\,N(\mathbf{r},\mathbf{b},Y)\,N(\mathbf{r}^{\prime},\mathbf{b},Y)\bigg{|}_{r_{\perp}\ll
1/Q_{s},r_{\perp}^{\prime}\gtrsim 1/Q_{s}}\propto A.$ (31)

Hence, we conclude that the elastic vector meson production cross section
scales with $A$ as a power of $A$,

$\displaystyle\sigma_{\textrm{el}}^{\gamma^{*}{\rm A}\to V{\rm A}}\propto
A^{\alpha},$ (32)

with $\alpha$ between $2/3$ and $4/3$. The precise power of the scaling
depends on the size of the vector meson: if the size of the vector meson is
small (e.g., $J/\psi$), then the integral contribution would be dominated by
region (i), and $\sigma_{\textrm{el}}^{\gamma^{*}{\rm A}\to J/\psi\,{\rm
A}}\propto A^{4/3}$; if the size of the vector meson is large (e.g., $\rho$),
then the integral contribution would be dominated by region (iii), and
$\sigma_{\textrm{el}}^{\gamma^{*}{\rm A}\to\rho\,{\rm A}}\propto A^{2/3}$.
Therefore, a transition from outside the saturation region into the saturation
region should lead to the decrease of the (effective) power $\alpha$ defined
in Eq. (32).

Notice that in Eq. (17) the integrand as a function of the dipole sizes
$r_{\perp}$ and $r_{\perp}^{\prime}$ is dominated by the Gaussian and the
modified Bessel functions (which decrease exponentially at large $r_{\perp}$
and $r_{\perp}^{\prime}$), so that the main contribution comes from the
regions where $r_{\perp},r_{\perp}^{\prime}<\frac{1}{a_{f}},R$. For $J/\psi$
production in UPCs, where $Q^{2}\approx 0$ and $a_{f}\approx m_{c}\approx
1.27\ \text{GeV}$, this corresponds to
$r_{\perp},r_{\perp}^{\prime}<\frac{1}{m_{c}}\approx 0.79\text{ GeV}^{-1}$. At
relatively low $x$ ($x$ between $10^{-3}$ and $10^{-4}$), the typical
saturation scale for a gold nucleus ($A=197$) is about $Q_{s}\approx 1\
\text{GeV}$ (see, e.g., Fig. 3.14 in [9]). We see that the
$r_{\perp},r_{\perp}^{\prime}$-integrals in Eq. (17) are dominated by the non-
saturated region (i), so that $\sigma^{\gamma^{*}A\to J/\psi\,A}\propto
A^{4/3}$. However, these integrals do include contributions from larger
$r_{\perp},r^{\prime}_{\perp}$, coming from the saturation region. Therefore,
in an exact evaluation of Eq. (17), one may expect to see an $A$-scaling that
is slightly slower than $A^{4/3}$, especially at the largest $A$ when
$1/Q_{s}$ starts to become comparable to the size of $J/\psi$ and saturation
effects start to settle in.

Figure 7: A log-log plot of the elastic $J/\psi$ production cross section
versus the atomic number $A$ of the target nucleus. The quasi-classical dipole
amplitude used here is given in Eq.(33).

The results of our analysis above are supported by numerical calculations. To
evaluate Eq. (17) numerically for the case of UPCs, we put $Q^{2}=900$ MeV2.
The parameters of the vector meson wave functions given by Eq. (15) are listed
here in Table 1, taken from [65].

Meson | $m_{f}/\rm{GeV}$ | ${\cal N}_{T}$ | ${R}^{2}/\rm{GeV}^{-2}$
---|---|---|---
$J/\psi$ | 1.27 | 0.578 | 2.3
$\rho$ | 0.14 | 0.911 | 12.9
Table 1: Parameters for $J/\psi$ and $\rho$ mesons wave functions given in Eq.
(15). Values are taken from [65].

We use the quasi–classical dipole amplitude that is given by (see, e.g., [69])

$N(\mathbf{r},\mathbf{b},Y)=1-\exp{-\frac{r_{\perp}^{2}Q_{s}^{2}(\mathbf{b})}{4}\ln\left(\frac{1}{r_{\perp}\Lambda}+e\right)},$
(33)

to be able to integrate over $r_{\perp}>1/\Lambda$ without having an
unphysical negative $N$ in that region (cf. Eq. (22)). For forward kinematics
at RHIC we use the value of the saturation scale for Au nuclei at $x\approx
2\times 10^{-3}$ and $b_{\perp}=0$ from Fig. 3.9 of [9] to write

$\displaystyle
Q_{s}^{2}(\mathbf{b})=(1\,\text{GeV}^{2})\,\frac{T(\mathbf{b})}{2R_{A}\,\rho}\,\left(\frac{A}{197}\right)^{1/3}\,\theta(R_{A}-b_{\perp})=(1\,\text{GeV}^{2})\,\sqrt{1-\frac{b_{\perp}^{2}}{R_{A}^{2}}}\,\left(\frac{A}{197}\right)^{1/3}\,\theta(R_{A}-b_{\perp}),$
(34)

where, for simplicity, we consider the nucleus to be a sphere of radius
$R_{A}$ with a uniform nucleon number density $\rho=0.40$ fm${}^{-3}=0.0031$
GeV3, such that $T(\mathbf{b})=2\rho\,\sqrt{R_{A}^{2}-b_{\perp}^{2}}$. We
approximate $R_{A}\approx 1.2A^{1/3}$ fm$=6\,A^{1/3}$ GeV-1 and choose
$\Lambda=0.2$ GeV for the IR cutoff. We also use $\alpha_{EM}=1/137$,
$\alpha_{s}=1/3$, $N_{c}=3$, and $Z_{f}=2/3$ for $J/\psi$.

In Fig. 7 we present a log-log plot of the $J/\psi$ elastic production cross
section versus $A$. The slope between the two adjacent data points starts at
about $\alpha=1.26$ at $A=1$, then gradually decreases to about $\alpha=1.12$
at $A=208$ (Pb). (The power $\alpha$ is defined in Eq. (32) above.) We see
that the $A$-scaling power $\alpha$ of the elastic $J/\psi$ production in UPCs
is smaller than the (analytical) non-saturated value of 4/3, indicating that
there are some saturation effects present in the process. The decrease of the
slope with increasing $A$ is indeed what we expected, since saturation is more
manifest in a larger nucleus. Note also that the magnitude of the cross
section $\sigma_{\textrm{el}}^{\gamma^{*}{\rm A}\to J/\psi\,{\rm A}}$ for Au
nuclei, $A=197$, is only off by a factor of 3 to 4 from the recent
experimental data recently reported by the STAR collaboration, despite the
relative crudeness of the approximations we have made here.

Figure 8: A log-log plot of the elastic $\rho$ production cross section versus
the atomic number $A$ of the target nucleus. The quasi-classical dipole
amplitude used here is given in Eq.(33).

Similarly, we perform a numerical evaluation of Eq. (17) for the $\rho$-meson
using the same quasi-classical dipole amplitude (33). The $\rho$-meson wave
function parameters are also given in Table 1. For the $\rho$ meson one uses
an effective value of $Z_{f}=1/\sqrt{2}$ [65]. In Fig. 8 we show the log-log
plot of the elastic $\rho$ production cross section in UPCs versus $A$. The
slope between the adjacent data points starts at about $\alpha=0.90$ at $A=1$,
and gradually decreases to about $\alpha=0.75$ at $A=208$. We see that the
$A$-scaling power $\alpha$ is smaller than that for $J/\psi$, as we expected,
since $\rho$ is a larger meson which is more affected by saturation effects,
which tend to reduce $\alpha$. Moreover, extrapolating the plot in Fig. 8 to
the unrealistically high $A=300$, we see that the power $\alpha$ approaches
the “ideal” value of 2/3 for a fully saturated process in the regime (ii)
considered above, thus agreeing with our analytic estimate.

To summarize our observations here, we see that saturation effects tend to
reduce the effective power $\alpha$ from Eq. (32). A transition to the
saturation regime can happen in two ways. We can increase the atomic number
$A$ of the nucleus, while considering the elastic production cross section for
the same vector meson. Saturation would lead to the decrease of the slope
$\alpha$ with $A$, as observed separately in Fig. 7 and in Fig. 8.
Alternatively, saturation can be reached by increasing the size of the vector
meson. Hence, in going from $J/\psi$ to $\rho$ meson, by comparing Fig. 7 to
Fig. 8, we see that $\alpha$ decreases as well between these two figures due
to saturation effects.

####  IV.1.2 Inclusive production cross section

Now let us use the above estimates to evaluate the $A$-scaling of the double
ratio $R_{\rm{UPC}}$ from Eq. (3). Since the elastic cross section has been
found above, we also need the $A$ scaling of the inelastic hadron/jet
production cross section. To obtain that we will use Eq. (20), in which the
$r_{\perp}$ integral is restricted to $r_{\perp}<\min\\{1/p_{T},1/a_{f}\\}$ by
the Fourier exponential and by the modified Bessel functions coming from the
virtual photon wave function. Since the typical jet or net hadron production
cross section is dominated by light flavors, we see that in UPCs $a_{f}$ for
Eq. (20) would be rather small, of the order of several tens of MeV, such that
$a_{f}\ll p_{T}$: hence, we only have the $r_{\perp}<1/p_{T}$ constraint on
the range of $r_{\perp}$, before we account for the dipole amplitude $N$.

For $p_{T}\gg Q_{s}$ we have $r_{\perp}\ll 1/Q_{s}$ in Eq. (20), placing us
outside the saturation region. Using the approximation from Eqs. (26) for the
dipole amplitude, we see that $N\propto A^{1/3}$, such that

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}\bigg{|}_{p_{T}\gg Q_{s}}\propto A$
(35)

after we integrate over the impact parameters $\mathbf{b}$ between the dipole
and the nucleus.

For low momenta, $p_{T}\lesssim Q_{s}$ (while still $p_{T}\gg\Lambda_{QCD}$
with $\Lambda_{QCD}$ the QCD confinement scale), we have the typical
$r_{\perp}\sim 1/p_{T}\gtrsim 1/Q_{s}$, such that the $r_{\perp}$ integral in
Eq. (20) is largely inside the saturation region. We then apply the
approximation $N\approx 1$ in Eq. (20) (see (28)) to obtain

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}\bigg{|}_{p_{T}\lesssim Q_{s}}\propto
A^{2/3}$ (36)

after integration over $\mathbf{b}$. The results in Eqs. (35) and (36) are
consistent with the expectations of saturation physics for particle production
in dilute-dense scattering, like the proton–nucleus collisions considered in
[47]. Similar to the elastic vector meson production process considered above,
saturation effects tend to reduce the power of $A$ for the inclusive
production cross section as well.

####  IV.1.3 $A$-scaling for the double ratio for $J/\psi$ and $\rho$
production

Let us first summarize our results for the elastic vector meson production and
the inelastic hadron/jet production, before constructing cross section ratios.
The elastic vector production cross section scales approximately as

$\displaystyle\sigma^{\gamma^{*}{\rm A}\to V{\rm
A}}_{\textrm{el}}\propto\begin{cases}A^{4/3},&\text{outside the saturation
region},\\\ A^{2/3},&\text{inside the saturation region}.\end{cases}$ (37)

Here transition in or out of the saturation region is obtained by varying $A$
and/or the typical dipole size $r_{\perp}$. The latter is varied in UPCs by
choosing different-size vector mesons, with $J/\psi$ being smaller than
$\rho$. The inelastic cross section has an $A$-scaling that is dependent on
the transverse momentum $p_{T}$ of the produced quark:

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}\propto\begin{cases}A,&p_{T}\gg
Q_{s},\\\ A^{2/3},&p_{T}\ll Q_{s}.\end{cases}$ (38)

The transition into the saturation regime depends on the value of $p_{T}$
compared to $Q_{s}$, with the latter depending on $A$ as well.

Taking the ratio of the elastic vector meson production cross section to the
inelastic hadron/jet production cross section, we have the following results.
For $J/\psi$ production, assuming that $J/\psi$ is small enough for the
process to be outside the saturation region, we first define and calculate the
single ratio

$R_{1}^{J/\psi}({\rm A})\equiv\frac{\sigma^{\gamma^{*}{\rm A}\to J/\psi\,{\rm
A}}_{\textrm{el}}}{\frac{d\sigma}{d^{2}p_{T}}}\propto\begin{cases}A^{\frac{1}{3}},&p_{T}\gg
Q_{s},\\\ A^{\frac{2}{3}},&p_{T}\ll Q_{s}.\end{cases}$ (39)

If we take the double ratio proposed in Eq. (3), the prefactors which depend
on $p_{T}$ only would mostly drop out (given that we evaluate the inelastic
cross sections at the same $p_{T}$) and we obtain the following $A$-scaling
for $R_{\rm{UPC}}$,

$R_{\rm{UPC}}^{J/\psi}\equiv\frac{R_{1}^{J/\psi}({\rm
A})}{R_{1}^{J/\psi}(p)}=\begin{cases}A^{\frac{1}{3}},&p_{T}\gg Q_{s},\\\
A^{\frac{2}{3}},&p_{T}\ll Q_{s}.\end{cases}$ (40)

We see that inside the saturation region for the inclusive cross section
($p_{T}\ll Q_{s}$), the $A$ scaling of $R_{\rm{UPC}}^{J/\psi}$ comes with a
higher power of $A$, than outside the saturation region ($p_{T}\gg Q_{s}$).

Figure 9: A plot of our approximate estimates for the double ratio
$R_{\rm{UPC}}$ as a function of $A$, with $R_{\rm{UPC}}$ given by Eqs. (40)
and (41) in the quasi-classical approximation.

For the $\rho$ meson production, assuming it is large enough for its elastic
production cross section to mainly probe the saturation region, we have the
following double ratio due to a different $A$-scaling in the elastic process,

$R_{\rm{UPC}}^{\rho}=\begin{cases}A^{-\frac{1}{3}},&p_{T}\gg Q_{s},\\\
A^{0}=1,&p_{T}\ll Q_{s}.\end{cases}$ (41)

Again, we see that the power of $A$ is lower outside the saturation region, by
the same amount of $1/3$ as in Eq. (40).

Note also, that in going from $J/\psi$ to the larger $\rho$ meson, the powers
of $A$ decreased in $R_{\rm{UPC}}$. So the double ratio appears to be
sensitive to how one approaches the saturation region: while lowering $p_{T}$
leads to the double ratio having a higher power of $A$, increasing the dipole
(meson) size leads to a decrease of the double ratio for the same $A$ and
$p_{T}$.

A plot of our result for $R_{\rm{UPC}}$ in Eqs. (40) and (41) is shown in Fig.
9, demonstrating different $A$ dependence of the double ratio for $J/\psi$ and
$\rho$ mesons. Note that a more careful numerical evaluation of the double
ratio, akin to that done in Figs. 7 and Fig. 8, is left for the future work:
for the realistic inelastic hadron production cross section one would need to
augment Eq. (20) by the fragmentation functions. In RHIC kinematics, large-$x$
effects may need to be included as well, for the higher $p_{T}$ values.
Additionally, small-$x$ evolution has to be included into the dipole amplitude
$N$ for very small $x$ values, as reached in the forward direction at RHIC and
LHC. This is what we will do next, to test the validity of our conclusions
here.

###  IV.2 Double ratio with small-$x$ evolution

####  IV.2.1 Elastic vector meson production

The results of Section IV.1, given by Eqs. (40) and (41), were obtained using
the quasi-classical dipole scattering amplitude (22). However, in practice,
radiative corrections need to be incorporated if one wants to obtain a more
realistic result, especially at small $x$. In the dipole picture, where the
color dipole collides with the target nucleus, the dipole radiates gluons
during the process if there is a sufficient scattering energy available. These
gluon emissions/absorptions result in corrections to the dipole amplitude
$N(\mathbf{r},\mathbf{b},Y)$ that are described by the BK/JIMWLK evolution
equations [15, 16, 17, 18, 19, 20, 21, 22, 23, 24], making
$N(\mathbf{r},\mathbf{b},Y)$ depend on the rapidity interval $Y$ between the
dipole and the target nucleus. One consequence of the BK evolution is the
extended geometric scaling behavior outside the saturation region (see [70,
49, 50, 71]) and the geometric scaling behavior inside the saturation region
(see [72]). Outside the saturation region in the extended geometric region the
dipole amplitude takes the form [51]

$N(\mathbf{r},\mathbf{b},Y)\propto\left[r_{\perp}^{2}Q_{s}^{2}(\mathbf{b},Y)\right]^{\gamma_{cr}}$
(42)

where $\gamma_{cr}\approx 0.6275$ is the solution of

$\displaystyle\chi(\gamma_{cr})=\gamma_{cr}\,\frac{d\chi(\gamma_{cr})}{d\gamma_{cr}}$
(43)

with $\chi(\gamma)$ being the eigenvalue of the LO BFKL kernel given by

$\chi(\gamma)=2\,\psi(1)-\psi\left(\gamma\right)-\psi\left(1-\gamma\right),\hskip
28.45274pt\text{where}\hskip
28.45274pt\psi(\gamma)=\frac{\Gamma^{\prime}(\gamma)}{\Gamma(\gamma)}.$ (44)

With the small-$x$ evolution included, the saturation scale $Q_{s}$ is now a
function of rapidity $Y$ given by

$Q_{s}(\mathbf{b},Y)=Q_{s}(\mathbf{b})\,\exp\left[\frac{\alpha_{s}N_{c}}{\pi}\frac{\chi(\gamma_{cr})}{2\gamma_{cr}}Y\right]\approx
Q_{s}(\mathbf{b})\,\exp\left[2.44\frac{\alpha_{s}N_{c}}{\pi}Y\right]$ (45)

where $Q_{s}(\mathbf{b})$ is the quasi–classical saturation scale from Eq.
(34) in Sec. IV.1. Note that we still have $Q_{s}(\mathbf{b},Y)^{2}\propto
A^{1/3}$. We see that $N(\mathbf{r},\mathbf{b},Y)$ in Eq. (42) depends only on
the combined quantity $r_{\perp}Q_{s}(\mathbf{b},Y)$, instead of depending on
two variables $r_{\perp}$ and $Y$ separately, a behavior known as geometric
scaling [73, 72, 70, 49, 50, 71, 51]. Note that our treatment of the impact
parameter $\mathbf{b}$ dependence of the saturation scale is approximate [74,
75, 76, 77, 78, 79], but should suffice in the spirit of our mainly
qualitative arguments here.

Deep inside the saturation region, the dipole amplitude is given by the Levin-
Tuchin (LT) formula [72],

$N(\mathbf{r},\mathbf{b},Y)=1-S_{0}\exp\left\\{-\frac{\gamma_{cr}}{2\,\chi(\gamma_{cr})}\,\ln^{2}[r_{\perp}^{2}Q_{s}^{2}(\mathbf{b},Y)]\right\\},$
(46)

where $S_{0}$ is a constant specified along the saturation boundary curve
$r_{\perp}=1/Q_{s}(\mathbf{b},Y)$. Again, we see the geometric scaling
behavior: $N(\mathbf{r},\mathbf{b},Y)$ depends only on the combined quantity
$r_{\perp}Q_{s}(\mathbf{b},Y)$. For a detailed discussion and derivation of
these scaling results see Chapter 4.5 in [7].

Let us now examine the $A$-dependence of the elastic vector meson production
in UPCs, if we use the dipole amplitude in the scaling region, given by either
Eq. (42) or Eq. (46), in Eq. (17). Going back to Eq. (17), with the dipole
amplitude now depending on $Y$, we apply the same arguments as in Sec. IV.1.
We divide the integration over $r_{\perp}$ and $r_{\perp}^{\prime}$ into three
regions: (i) $r_{\perp},r_{\perp}^{\prime}\ll 1/Q_{s}$, (ii)
$r_{\perp},r_{\perp}^{\prime}\gtrsim 1/Q_{s}$, (iii) $r_{\perp}\ll
1/Q_{s},r_{\perp}^{\prime}\gtrsim 1/Q_{s}$ or vice versa. For $J/\psi$
production, assuming it is small enough to be entirely outside the saturation
region, the integral in Eq. (17) is dominated by region (i), so we use Eq.
(42) to approximate the dipole amplitude $N(\mathbf{r},\mathbf{b},Y)$. Since
$Q_{s}^{2}(\mathbf{b},Y)\propto A^{1/3}$, we have $N\sim A^{\gamma_{cr}/3}$,
such that $\sigma^{\gamma^{*}{\rm A}\to J/\psi\,{\rm A}}_{\textrm{el}}\propto
A^{\frac{2}{3}(1+\gamma_{cr})}$. For $\rho$ production, assuming the dipole is
large enough to be inside the saturation region, the integral in Eq. (17) is
dominated by region (ii), so we use the approximation that
$N(\mathbf{r},\mathbf{b},Y)\approx 1$ (instead of the LT formula) for
convenience. We would then have $\sigma^{\gamma^{*}{\rm A}\to\rho\,{\rm
A}}_{\textrm{el}}\propto A^{\frac{2}{3}}$, same as in the quasi–classical
regime of Sec. IV.1.

Since the $r_{\perp},r_{\perp}^{\prime}$ integrals in Eq. (17) run over the
entire range of $r_{\perp},r_{\perp}^{\prime}\in[0,\infty]$, we would like to
perform a numerical integration to obtain a more precise estimate of the
$A$-scaling for the elastic cross sections. This time, we use a dipole
amplitude that is a solution of the BK evolution equation. Unfortunately, the
exact analytic solution to the BK equation is not known, so we have to either
solve the BK equation numerically, or resort to an approximation. In our
present exploratory study, we chose the latter, and approximate the dipole
scattering amplitude by piecing together the extended geometric scaling
behavior formula (42) outside the saturation region ($r_{\perp}\lesssim
1/Q_{s}$) with the LT formula (46) inside the saturated regime
($r_{\perp}\gtrsim 1/Q_{s}$). In addition, we require the dipole amplitude to
be continuous at $r_{\perp}=1/Q_{s}$. We thus write

$\displaystyle N(\mathbf{r},\mathbf{b},Y)=$ $\displaystyle\
(1-S_{0})\,[r_{\perp}^{2}Q_{s}^{2}(\mathbf{b},Y)]^{\gamma_{cr}}\,\theta\left(1-r_{\perp}Q_{s}(\mathbf{b},Y)\right)$
(47)
$\displaystyle+\left[1-S_{0}\exp\left\\{-\frac{\gamma_{cr}}{2\,\chi(\gamma_{cr})}\,\ln^{2}[r_{\perp}^{2}Q_{s}^{2}(\mathbf{b},Y)]\right\\}\right]\,\theta\left(r_{\perp}Q_{s}(\mathbf{b},Y)-1\right).$

Here $S_{0}$ is some constant determined by the initial condition for the
forward scattering amplitude. Since the curve $Q_{s}(Y)$ which describes the
rapidity dependence of the saturation scale $Q_{s}$ is determined by requiring
the dipole amplitude $N(1/Q_{s}(Y),Y)$ to be constant along this curve, we
have the freedom to choose what this constant value is as long as it is of the
order of 1. This, in turn, leads to the freedom in $S_{0}$. A typical choice
of $S_{0}$ motivated by the GGM/MV model (cf. Eq. (22)) is to take
$S_{0}=e^{-\frac{1}{4}}.$ Using this $S_{0}$ in Eq. (47) and employing the
resulting dipole amplitude to perform the numerical integrations in Eq. (17)
leads to the log-log plot of $J/\psi$ elastic production cross section versus
$A$ depicted in Fig. 10 for a sample of small $x$ values,
$x=\\{10^{-2},10^{-3},10^{-4},10^{-5}\\}$. The saturation scale is given by
Eq. (45) with $Q_{s}(\mathbf{b})$ from Eq. (34). We put $\alpha_{s}=0.12$ in
Eq. (45) to account for the effects of higher-order corrections on the
$x$-dependence of the saturation scale in our leading-order formalism. (Note
that $Y\approx\ln(1/x)$.)

Figure 10: A log-log plot of the elastic $J/\psi$ production cross section
versus the atomic number $A$ of the target nucleus for four different
$x$-values, $x=\\{10^{-2},10^{-3},10^{-4},10^{-5}\\}$. The dipole amplitude
used here includes small-$x$ evolution and is given in Eq. (47).

In Fig. 10 we see the expected trend of the scattering cross section
increasing with increasing scattering energy (with decreasing $x$). The slopes
$\alpha$ from Eq. (32), calculated again using pairs of adjacent data points
in Fig. 10, is summarized in Table 2. We observe that the slope at larger
values of $x$, $x=10^{-2}$ or $x=10^{-3}$, is very close to the estimated
value of $\frac{2}{3}(1+\gamma_{cr})\approx 1.085$ in region (i) considered
above (outside of the saturation region), indicating that the elastic $J/\psi$
production process for larger $x$ is indeed non-saturated in our picture, but
is still inside the geometric scaling region. However, as we go to lower $x$,
the scaling power $\alpha$ decreases significantly, becoming much smaller than
in the quasi-classical case, while still remaining significantly above the
fully-saturated value of $2/3$. This implies that the increase of $Q_{s}$ at
small $x$ moves the $J/\psi$ production process more into the saturated
region.

$x$ | slope at $A=1$ | slope at $A=208$
---|---|---
$10^{-2}$ | 1.08 | 1.06
$10^{-3}$ | 1.06 | 1.04
$10^{-4}$ | 1.05 | 0.96
$10^{-5}$ | 0.98 | 0.85
Table 2: The effective scaling power $\alpha$ of $A$ for the $J/\psi$ elastic
production cross section (the slopes in Fig. 10) for different $x$ values:
$\sigma_{\textrm{el}}^{\gamma^{*}{\rm A}\to J/\psi\,{\rm A}}\propto
A^{\alpha}$. Figure 11: A log-log plot of the elastic $\rho$ meson production
cross section as a function of the atomic number $A$ of the target nucleus for
four different $x$-values, $x=\\{10^{-2},10^{-3},10^{-4},10^{-5}\\}$. The
dipole amplitude used here includes small-$x$ evolution and is given in Eq.
(47).

Performing a similar numerical integration by using the amplitude (47) in Eq.
(17), but now for the elastic $\rho$ meson production yields the log-log plot
shown in Fig. 11. The slope of the data points is summarized in Table 3. The
scaling powers $\alpha$ are smaller than those for $J/\psi$ production, which
is expected since saturation effects are stronger for the larger mesons like
$\rho$. For small $x$ values such as $x=10^{-5}$, the slope is very close to
the $2/3$ scaling power for a fully saturated process (region (ii) above).
However for larger $x$, even for $x=10^{-3}$, it deviates quite noticeably
from the above estimate of $\alpha=2/3$ for a perfectly saturated process, and
$\alpha$ is larger than that extracted from Fig. 8 for the quasi–classical
case. This is because $N(\mathbf{r},\mathbf{b},Y)$ approaches 1 more slowly in
the LT formula (cf. Eq. (46)) than in the quasi-classical case (cf. Eq. (33)):
in the latter case, $N$ approaches 1 exponentially in $r_{\perp}^{2}$ instead
of the exponential of $\ln^{2}[r_{\perp}^{2}\,Q_{s}^{2}]$ in the former case
(the LT formula). Hence, the saturation effects in terms of $A$-scaling
manifest themselves less when evolution corrections are included.

$x$ | slope at $A=1$ | slope at $A=208$
---|---|---
$10^{-2}$ | 1.03 | 0.90
$10^{-3}$ | 0.92 | 0.79
$10^{-4}$ | 0.80 | 0.72
$10^{-5}$ | 0.73 | 0.69
Table 3: The effective scaling power $\alpha$ of $A$ for the $\rho$ meson
elastic production cross section (the slopes in Fig. 11) for different $x$
values: $\sigma_{\textrm{el}}^{\gamma^{*}{\rm A}\to\rho\,{\rm A}}\propto
A^{\alpha}$.

Still, the overall pattern for the elastic vector meson production remains the
same after inclusion of small-$x$ evolution corrections: the effect of
saturation is to reduce the power of $A$ in the cross section. Again, this
reduction can be achieved by either increasing $A$ or by increasing the dipole
size (in going from $J/\psi$ to $\rho$). Additionally, with small-$x$
evolution included, the reduction of the power of $A$ can be achieved by going
to smaller $x$.

####  IV.2.2 Inclusive production cross section

Similar to what we did for the elastic vector meson production, we can include
small-$x$ evolution in the inclusive production cross section (20) by using
the dipole amplitude $N(\mathbf{r},\mathbf{b},Y)$ given by the solution of the
BK equation. Again, we separately consider the $p_{T}\gg Q_{s}$ and
$\Lambda_{QCD}\ll p_{T}\ll Q_{s}$ cases.

Outside the saturation region, when $p_{T}\gg Q_{s}$, the contribution to the
$r_{\perp}$-integral in Eq. (20) mainly comes from the region where
$r_{\perp}\ll 1/p_{T}\ll 1/Q_{s}$, so we are in the extended geometric scaling
region. Hence we approximate $N(\mathbf{r},\mathbf{b},Y)$ by Eq. (42), which
leads to $N\propto A^{\frac{1}{3}\gamma_{cr}}$, and consequently

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}\bigg{|}_{p_{T}\gg Q_{s}}\propto
A^{\frac{2}{3}+\frac{1}{3}\gamma_{cr}}.$ (48)

Inside the saturation region, when $\Lambda_{QCD}\ll p_{T}\ll Q_{s}$, the
contribution to the integral in Eq. (20) mainly comes from the region where
$1/Q_{s}\lesssim r_{\perp}\lesssim 1/p_{T}$, so we are probing the saturation
region. Hence we approximate $N(\mathbf{r},\mathbf{b},Y)\approx 1$, and
consequently

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}\bigg{|}_{p_{T}\lesssim Q_{s}}\propto
A^{2/3},$ (49)

same as in the quasi–classical case.

####  IV.2.3 $A$-scaling for the double ratio for $J/\psi$ and $\rho$
production

Let us combine the results above for the elastic vector meson production and
the inclusive hadron/jet production. The analytic results of Sec. IV.2.1 can
be summarized as

$\displaystyle\sigma^{\gamma^{*}{\rm A}\to V{\rm
A}}_{\textrm{el}}\propto\begin{cases}A^{\frac{2}{3}(1+\gamma_{cr})},&\text{outside
the saturation region},\\\ A^{\frac{2}{3}},&\text{inside the saturation
region}.\end{cases}$ (50)

Similarly, the results of Sec. IV.2.2 are

$\displaystyle\frac{d\sigma}{d^{2}p_{T}}\propto\begin{cases}A^{\frac{2}{3}+\frac{1}{3}\gamma_{cr}},&p_{T}\gg
Q_{s},\\\ A^{\frac{2}{3}},&p_{T}\ll Q_{s}.\end{cases}$ (51)

We see that when the small-$x$ evolution is taken into account, we have the
following single ratio of the elastic $J/\psi$ production cross section to the
inelastic jet/hadron production cross section:

$R_{1}^{J/\psi}({\rm A})=\frac{\sigma^{\gamma^{*}{\rm A}\to J/\psi\,{\rm
A}}_{\textrm{el}}}{\frac{d\sigma}{d^{2}p_{T}}}\propto\begin{cases}A^{\frac{1}{3}\gamma_{cr}},&p_{T}\gg
Q_{s},\\\ A^{\frac{2}{3}\gamma_{cr}},&p_{T}\ll Q_{s}.\end{cases}$ (52)

Taking the double ratio gives

$R_{\rm{UPC}}^{J/\psi}=\begin{cases}A^{\frac{1}{3}\gamma_{cr}},&p_{T}\gg
Q_{s},\\\ A^{\frac{2}{3}\gamma_{cr}},&p_{T}\ll Q_{s}.\end{cases}$ (53)

For the $\rho$ meson, the double ratio now is

$R_{\rm{UPC}}^{\rho}=\begin{cases}A^{-\frac{1}{3}\gamma_{cr}},&p_{T}\gg
Q_{s},\\\ A^{0}=1,&p_{T}\ll Q_{s}.\end{cases}$ (54)

Figure 12: The double ratio $R_{\rm{UPC}}$ plotted versus $A$ for the case
where we include LO small-$x$ evolution into the elastic and inclusive cross
sections. The figure illustrates the results in Eqs. (53) and (54).

Using $\gamma_{cr}=0.6275$, we plot our results from Eqs. (53) and (54) for
$R_{\rm{UPC}}$ versus $A$ in Fig. 12. Figure 12 is qualitatively similar to
Fig. 9, with the somewhat lower powers of $A$ for some of the curves in Fig.
12. Further phenomenological analysis of the double ratio, which may involve
solving the evolution equations at next-to-leading (NLO) order [80, 81, 82,
83], along with constructing NLO expressions for the corresponding elastic and
inclusive cross sections (see, e.g., [84]), is left for the future work
improving on our present exploratory study.

##  V Future opportunities at RHIC and LHC

In the near term at RHIC, there are two running experiments, sPHENIX and STAR,
taking high-luminosity Au$+$Au and p$+$p data at 200 GeV center-of-mass energy
for Run 2024 and Run 2025. There is a possibility of a $p+$Au run in 2024, if
condition allows, which will be valuable for the baseline measurement of
scattering on the proton. These two runs will provide the highest-statistics
UPC samples for the entire RHIC era.

Specifically, the sPHENIX experiment is a brand new detector with the
capability of performing UPC measurements, while the STAR experiment has a
well-established UPC program and the newly installed forward detector system
will provide a large extension in the kinematic phase space for UPCs. The
advantage of UPC measurements at RHIC is the unique kinematics that
transitions from high-$x$ to low-$x$, where saturation physics is expected to
start playing a role ($x<0.01$). Furthermore, with the extension of the STAR
forward rapidity coverage, measurement of the UPC jet or semi-inclusive hadron
photoproduction events may become possible with dedicated triggers. Together
with VM production, e.g., $\rho$ and $J/\psi$, the single and double ratio
observables we proposed here can be measured experimentally. Due to the
excellent tracking system in STAR, it is possible to measure the
photoproduction of $\phi$ meson with soft kaon decays for the first time
(average momentum for the kaon daughter is around 100 $\rm MeV/c$), which will
provide an extra handle on the mass or dipole size dependence of our proposed
observables.

At the LHC, the Run 3 has just started with three Heavy-Ion runs scheduled.
All four LHC experiments have been participating in UPC measurements [85, 86,
87, 88, 89, 90, 91, 92, 93, 31]. Due to the higher center-of-mass energy, the
kinematic phase space for most experiments is expected to be deep into the
saturation regime. There is even a hint of saturation phenomenon recently
reported by the CMS experiment [29], although the interpretation is not
entirely clear. The ATLAS experiment, although it has not yet reported any VM
photoproduction measurement, has pioneered UPC jets study in the past years
based on the data from Run 2. There is a good experimental opportunity for the
ATLAS experiment to simultaneously measure elastic VM and inclusive jet
productions in UPCs to study the single and double ratios proposed in this
paper. In the longer term, when the LHC Run-4 data will become available,
there will be still a few years before the operation of the EIC. At Run-4, the
ALICE experiment may have a Forward Calorimeter (FoCal) upgrade [94] that will
be able to probe $x$ down to $x\sim 10^{-6}$. This will be another great
opportunity for testing our proposed measurements quantitatively.

##  VI Summary and Outlook

In this paper, inspired by the double ratio of the diffractive to total cross
section in DIS to be measured at the EIC [9], we have proposed a novel double
ratio, $R_{\rm{UPC}}({\rm A})$ in Eq. (3), that is also sensitive to
saturation effects and can be measured in the UPCs at RHIC and LHC. More
precisely, we have performed a calculation of $R_{\rm{UPC}}({\rm A})$ for
$J/\psi$ and $\rho$ meson production in the saturation framework. Within this
framework, we have considered the quasi–classical approximation first, and
included the small-$x$ evolution later. Our results for the double ratios are
given in Eqs. (40) and (41) in the quasi–classical case and by Eqs. (53) and
(54) for the small-$x$ evolution case. In both calculations, we observe a
different dependence of $R_{\rm{UPC}}({\rm A})$ on the atomic number $A$
inside and outside the saturation region.

Whether we are probing the physics inside or outside the saturation region is
determined by the transverse momentum $p_{T}$ of the hadron or jet: for
$p_{T}\lesssim Q_{s}$ the inclusive production cross section in the
denominator of the single ratio $R_{1}$ (defined in Eq. (39)) is largely
probing the physics inside the saturation region, while for $p_{T}\gg Q_{s}$
it is probing the physics outside the saturation region. The size of the
produced meson in the elastic cross section also affects whether saturation
region is probed or not in the elastic cross section in the numerator of
$R_{1}$. In the case of a relatively compact meson production, for $J/\psi$,
the double ratio $R_{\rm{UPC}}({\rm A})$ grows faster with $A$ inside the
saturation region ($p_{T}\lesssim Q_{s}$) than outside ($p_{T}>Q_{s}$), as
follows from Eqs. (40) and (53). For a larger meson, like $\rho$, the double
ratio $R_{\rm{UPC}}({\rm A})$ decreases with $A$ outside the saturation region
($p_{T}>Q_{s}$), and remains relatively flat in $A$ inside the saturation
region ($p_{T}\lesssim Q_{s}$), as per Eqs. (41) and (54).

If our, admittedly preliminary, estimates are confirmed by the future more
detailed calculations, possibly involving NLO corrections to the cross
sections and evolution equations, then one may be able to conclude that the
$A$-dependence of $R_{\rm{UPC}}({\rm A})$ would be a novel signal of
saturation physics, to be explored in the current and future UPC measurements
at RHIC and LHC.

To prepare for the future experiments and to make sure the behavior predicted
here is unique to saturation physics, one has to compare the saturation
predictions for $R_{\rm UPC}$ to those coming from other calculations. In
particular, similar to [9], one may consider the LTS model [45, 46]. While our
$R_{\rm{UPC}}({\rm A})$ has not been calculated in the LTS model yet, a
similar (single) ratio of the diffractive to inclusive di-jet production cross
sections in UPCs was studied recently in [95]. The conclusion of [95] appears
to indicate that the diffractive to inclusive di-jet ratio in the LTS picture
decreases with the atomic number $A$ when comparing ultra-peripheral $p+p$ and
A$+$A collisions. This means that the corresponding double ratio, if
constructed, would be smaller than one and decreasing with $A$. This appears
to be different from the $A$-dependence we obtained above for $R_{\rm UPC}$ in
most (though not all) regimes considered, making us cautiously optimistic that
an LTS-based calculation for $R_{\rm UPC}$ would lead to a significantly
different $A$-dependence from that in the saturation picture, potentially
allowing one to distinguish between the two approaches using the future UPCs
data.

## Acknowledgments

The authors would like to thank Peter Steinberg for information about the UPC
program in the ATLAS experiment. We would like to thank the Brookhaven EIC
group for general discussion on UPCs and saturation physics.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Nuclear Physics under Award Number DE-SC0004286
and the work performed within the framework of the Saturated Glue (SURGE)
Topical Theory Collaboration (YK and HS).

The work of ZT is supported by the U.S. Department of Energy under Award DE-
SC0012704 and the Laboratory Directed Research LDRD-23-050 project.

##  Appendix A Calculational details

Using the following integral formulae,

$\displaystyle\int\frac{d^{2}q_{\perp}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q}\cdot\mathbf{x}}=2\pi\,K_{0}(a_{f}x_{\perp}),$
(55a) $\displaystyle\int
d^{2}q_{\perp}\frac{\mathbf{q}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q}\cdot\mathbf{x}}=2\pi\,i\,a_{f}\,K_{1}(a_{f}x_{\perp})\,\frac{\mathbf{x}}{x_{\perp}},$
(55b)

we can perform some of the integrals in Eq. (19) term by term:

$\displaystyle\int
d^{2}x_{\perp}^{\prime}d^{2}y_{\perp}e^{i\mathbf{p}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}\int
d^{2}x_{\perp}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{y})}\left\\{[z^{2}+(1-z)^{2}]a_{f}^{2}K_{1}(|\mathbf{x}^{\prime}-\mathbf{y}|a_{f})\frac{\mathbf{x}^{\prime}-\mathbf{y}}{|\mathbf{x}^{\prime}-\mathbf{y}|}\cdot\frac{1}{2\pi
ia_{f}}\right.$ (56a) $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\times\left.\int
d^{2}q_{\perp}\frac{\mathbf{q}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}+m_{f}^{2}K_{0}(|\mathbf{x}^{\prime}-\mathbf{y}|a_{f})\frac{1}{2\pi}\int\frac{d^{2}q_{\perp}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}\right\\}N\left(\mathbf{x}^{\prime}-\mathbf{y},\frac{\mathbf{x}^{\prime}+\mathbf{y}}{2},Y\right)$
$\displaystyle=$ $\displaystyle\int
d^{2}x_{\perp}^{\prime}d^{2}y_{\perp}e^{i\mathbf{p}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}\left\\{[z^{2}+(1-z)^{2}]a_{f}K_{1}(|\mathbf{x}^{\prime}-\mathbf{y}|a_{f})(-2\pi
i)\frac{\mathbf{p}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}{(p_{\perp}^{2}+a_{f}^{2})|\mathbf{x}^{\prime}-\mathbf{y}|}\right.$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left.\frac{2\pi
m_{f}^{2}}{p_{\perp}^{2}+a_{f}^{2}}K_{0}(|\mathbf{x}^{\prime}-\mathbf{y}|a_{f})\right\\}N\left(\mathbf{x}^{\prime}-\mathbf{y},\frac{\mathbf{x}^{\prime}+\mathbf{y}}{2},Y\right)$
$\displaystyle=2\pi\int
d^{2}x_{\perp}d^{2}x_{\perp}^{\prime}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\left\\{[z^{2}+(1-z)^{2}]\,a_{f}\,K_{1}(|\mathbf{x}-\mathbf{x}^{\prime}|a_{f})\,i\,\frac{\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}{(p_{\perp}^{2}+a_{f}^{2})|\mathbf{x}-\mathbf{x}^{\prime}|}+\frac{m_{f}^{2}}{p_{\perp}^{2}+a_{f}^{2}}K_{0}(|\mathbf{x}-\mathbf{x}^{\prime}|a_{f})\right\\}$
$\displaystyle\hskip
142.26378pt\times\,N\left(\mathbf{x}^{\prime}-\mathbf{x},\frac{\mathbf{x}+\mathbf{x}^{\prime}}{2},Y\right),$
$\displaystyle\int
d^{2}x_{\perp}d^{2}y_{\perp}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{y})}\int
d^{2}x_{\perp}^{\prime}e^{i\mathbf{p}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}\left\\{\left[z^{2}+(1-z)^{2}\right]a_{f}^{2}K_{1}(|\mathbf{x}-\mathbf{y}|a_{f})\frac{\mathbf{x}-\mathbf{y}}{|\mathbf{x}-\mathbf{y}|}\cdot\frac{1}{2\pi
ia_{f}}\right.$ (56b) $\displaystyle\hskip 85.35826pt\times\left.\int
d^{2}q_{\perp}\frac{\mathbf{q}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}+m_{f}^{2}K_{0}(|\mathbf{x}-\mathbf{y}|a_{f})\frac{1}{2\pi}\int\frac{d^{2}q_{\perp}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}\right\\}N\left(\mathbf{x}-\mathbf{y},\frac{\mathbf{x}+\mathbf{y}}{2},Y\right)$
$\displaystyle=\int
d^{2}x_{\perp}d^{2}y_{\perp}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{y})}\left\\{[z^{2}+(1-z)^{2}]a_{f}K_{1}(|\mathbf{x}-\mathbf{y}|a_{f})(2\pi
i)\frac{\mathbf{p}\cdot(\mathbf{x}-\mathbf{y})}{(p_{\perp}^{2}+a_{f}^{2})|\mathbf{x}-\mathbf{y}|}+\frac{2\pi
m_{f}^{2}}{p_{\perp}^{2}+a_{f}^{2}}K_{0}(|\mathbf{x}-\mathbf{y}|a_{f})\right\\}$
$\displaystyle\hskip
142.26378pt\times\,N\left(\mathbf{x}-\mathbf{y},\frac{\mathbf{x}+\mathbf{y}}{2},Y\right)$
$\displaystyle=2\pi\int
d^{2}x_{\prime}d^{2}x_{\perp}^{\prime}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\left\\{[z^{2}+(1-z)^{2}]a_{f}K_{1}(|\mathbf{x}-\mathbf{x}^{\prime}|a_{f})i\frac{\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}{(p_{\perp}^{2}+a_{f}^{2})|\mathbf{x}-\mathbf{x}^{\prime}|}+\frac{m_{f}^{2}}{p_{\perp}^{2}+a_{f}^{2}}K_{0}(|\mathbf{x}-\mathbf{x}^{\prime}|a_{f})\right\\}$
$\displaystyle\hskip
142.26378pt\times\,N\left(\mathbf{x}-\mathbf{x}^{\prime},\frac{\mathbf{x}+\mathbf{x}^{\prime}}{2},Y\right),$
$\displaystyle\int
d^{2}x_{\perp}d^{2}x_{\perp}^{\prime}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\int
d^{2}y_{\perp}\left\\{[z^{2}+(1-z)^{2}]a_{f}^{2}\left(\frac{1}{2\pi
ia_{f}}\right)^{2}\int
d^{2}q_{\perp}\frac{\mathbf{q}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q\cdot(\mathbf{x}-\mathbf{y})}}\right.$
(56c) $\displaystyle\times\,\int
d^{2}q_{\perp}^{\prime}\frac{\mathbf{q}^{\prime}}{q_{\perp}^{\prime
2}+a_{f}^{2}}e^{i\mathbf{q}^{\prime}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}\left.+m_{f}^{2}\left(\frac{1}{2\pi}\right)^{2}\int\frac{d^{2}q_{\perp}}{q_{\perp}^{2}+a_{f}^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}\int\frac{d^{2}q_{\perp}^{\prime}}{q_{\perp}^{\prime
2}+a_{f}^{2}}e^{i\mathbf{q}^{\prime}\cdot(\mathbf{x}^{\prime}-\mathbf{y})}\right\\}N\left(\mathbf{x}-\mathbf{x}^{\prime},\frac{\mathbf{x}+\mathbf{x}^{\prime}}{2},Y\right)$
$\displaystyle=\int
d^{2}x_{\perp}d^{2}x_{\perp}^{\prime}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\left\\{[z^{2}+(1-z)^{2}]\int
d^{2}q_{\perp}\frac{q_{\perp}^{2}}{(q_{\perp}^{2}+a_{f}^{2})^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}+m_{f}^{2}\int\frac{d^{2}q_{\perp}}{(q_{\perp}^{2}+a_{f}^{2})^{2}}e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\right\\}$
$\displaystyle\hskip
142.26378pt\times\,N\left(\mathbf{x}-\mathbf{x}^{\prime},\frac{\mathbf{x}+\mathbf{x}^{\prime}}{2},Y\right)$
$\displaystyle=\int
d^{2}x_{\perp}d^{2}x_{\perp}^{\prime}e^{-i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\left\\{[z^{2}+(1-z)^{2}]\pi\Big{[}2\,K_{0}(|\mathbf{x}-\mathbf{x}^{\prime}|a_{f})-|\mathbf{x}-\mathbf{x}^{\prime}|a_{f}\,K_{1}(|\mathbf{x}-\mathbf{x}^{\prime}|a_{f})\Big{]}\right.$
$\displaystyle\hskip
142.26378pt\left.+m_{f}^{2}\frac{\pi|\mathbf{x}-\mathbf{x}^{\prime}|}{a_{f}}K_{1}(|\mathbf{x}-\mathbf{x}^{\prime}|a_{f})\right\\}\,N\left(\mathbf{x}-\mathbf{x}^{\prime},\frac{\mathbf{x}+\mathbf{x}^{\prime}}{2},Y\right).$

Substituting Eqs. (56) into Eq. (19), relabeling variables appropriately, and
employing Eq. (21), leads to Eq. (20).

## References

  * [1] L. V. Gribov, E. M. Levin and M. G. Ryskin, _Semihard Processes in QCD_ , _Phys. Rept._ 100 (1983) 1–150.
  * [2] E. Iancu and R. Venugopalan, _The Color glass condensate and high-energy scattering in QCD_ , pp. 249–3363.  3, 2003.  hep-ph/0303204.  10.1142/9789812795533$\\_$0005\.
  * [3] H. Weigert, _Evolution at small $x_{bj}$: The Color Glass Condensate_, _Prog. Part. Nucl. Phys._ 55 (2005) 461–565, [hep-ph/0501087].
  * [4] J. Jalilian-Marian and Y. V. Kovchegov, _Saturation physics and deuteron-Gold collisions at RHIC_ , _Prog. Part. Nucl. Phys._ 56 (2006) 104–231, [hep-ph/0505052].
  * [5] F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan, _The Color Glass Condensate_ , _Ann.Rev.Nucl.Part.Sci._ 60 (2010) 463–489, [1002.0333].
  * [6] J. L. Albacete and C. Marquet, _Gluon saturation and initial conditions for relativistic heavy ion collisions_ , _Prog.Part.Nucl.Phys._ 76 (2014) 1–42, [1401.4866].
  * [7] Y. V. Kovchegov and E. Levin, _Quantum chromodynamics at high energy_ , vol. 33.  Cambridge University Press, 2012.
  * [8] A. Morreale and F. Salazar, _Mining for Gluon Saturation at Colliders_ , _Universe_ 7 (2021) 312, [2108.08254].
  * [9] A. Accardi et al., _Electron Ion Collider: The Next QCD Frontier_ , _Eur. Phys. J._ A52 (2016) 268, [1212.1701].
  * [10] D. Boer et al., _Gluons and the quark sea at high energies: Distributions, polarization, tomography_ , 1108.1713.
  * [11] A. Prokudin, Y. Hatta, Y. Kovchegov and C. Marquet, eds., _Proceedings, Probing Nucleons and Nuclei in High Energy Collisions: Dedicated to the Physics of the Electron Ion Collider: Seattle (WA), United States, October 1 - November 16, 2018_ , WSP, 2020.  10.1142/11684.
  * [12] R. Abdul Khalek et al., _Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report_ , 2103.05419.
  * [13] E. A. Kuraev, L. N. Lipatov and V. S. Fadin, _The Pomeranchuk singlularity in non-Abelian gauge theories_ , _Sov. Phys. JETP_ 45 (1977) 199–204.
  * [14] I. Balitsky and L. Lipatov, _The Pomeranchuk Singularity in Quantum Chromodynamics_ , _Sov.J.Nucl.Phys._ 28 (1978) 822–829.
  * [15] I. Balitsky, _Operator expansion for high-energy scattering_ , _Nucl. Phys._ B463 (1996) 99–160, [hep-ph/9509348].
  * [16] I. Balitsky, _Factorization and high-energy effective action_ , _Phys. Rev._ D60 (1999) 014020, [hep-ph/9812311].
  * [17] Y. V. Kovchegov, _Small-x $F_{2}$ structure function of a nucleus including multiple pomeron exchanges_, _Phys. Rev._ D60 (1999) 034008, [hep-ph/9901281].
  * [18] Y. V. Kovchegov, _Unitarization of the BFKL pomeron on a nucleus_ , _Phys. Rev._ D61 (2000) 074018, [hep-ph/9905214].
  * [19] J. Jalilian-Marian, A. Kovner and H. Weigert, _The Wilson renormalization group for low x physics: Gluon evolution at finite parton density_ , _Phys. Rev._ D59 (1998) 014015, [hep-ph/9709432].
  * [20] J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, _The Wilson renormalization group for low x physics: Towards the high density regime_ , _Phys. Rev._ D59 (1998) 014014, [hep-ph/9706377].
  * [21] H. Weigert, _Unitarity at small Bjorken x_ , _Nucl. Phys._ A703 (2002) 823–860, [hep-ph/0004044].
  * [22] E. Iancu, A. Leonidov and L. D. McLerran, _The renormalization group equation for the color glass condensate_ , _Phys. Lett._ B510 (2001) 133–144.
  * [23] E. Iancu, A. Leonidov and L. D. McLerran, _Nonlinear gluon evolution in the color glass condensate. I_ , _Nucl. Phys._ A692 (2001) 583–645, [hep-ph/0011241].
  * [24] E. Ferreiro, E. Iancu, A. Leonidov and L. McLerran, _Nonlinear gluon evolution in the color glass condensate. II_ , _Nucl. Phys._ A703 (2002) 489–538, [hep-ph/0109115].
  * [25] A. H. Mueller, _Small x Behavior and Parton Saturation: A QCD Model_ , _Nucl. Phys._ B335 (1990) 115.
  * [26] L. D. McLerran and R. Venugopalan, _Computing quark and gluon distribution functions for very large nuclei_ , _Phys. Rev._ D49 (1994) 2233–2241, [hep-ph/9309289].
  * [27] L. D. McLerran and R. Venugopalan, _Gluon distribution functions for very large nuclei at small transverse momentum_ , _Phys. Rev._ D49 (1994) 3352–3355, [hep-ph/9311205].
  * [28] L. D. McLerran and R. Venugopalan, _Green’s functions in the color field of a large nucleus_ , _Phys. Rev._ D50 (1994) 2225–2233, [hep-ph/9402335].
  * [29] CMS collaboration, _Photon-nucleus energy dependence of coherent J/ $\psi$ cross session in ultraperipheral PbPb collisions at $5.02~{}\mathrm{TeV}$ with CMS_, .
  * [30] ALICE collaboration, S. Acharya et al., _First measurement of the $|t|$-dependence of incoherent J/$\psi$ photonuclear production_, 2305.06169.
  * [31] ALICE collaboration, S. Acharya et al., _Energy dependence of coherent photonuclear production of J/ $\psi$ mesons in ultra-peripheral Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}$=5.02 TeV_, 2305.19060.
  * [32] STAR collaboration, M. S. Abdallah et al., _Evidence for Nonlinear Gluon Effects in QCD and Their Mass Number Dependence at STAR_ , _Phys. Rev. Lett._ 129 (2022) 092501, [2111.10396].
  * [33] C. A. Bertulani, S. R. Klein and J. Nystrand, _Physics of ultra-peripheral nuclear collisions_ , _Ann. Rev. Nucl. Part. Sci._ 55 (2005) 271–310, [nucl-ex/0502005].
  * [34] A. J. Baltz, _The Physics of Ultraperipheral Collisions at the LHC_ , _Phys. Rept._ 458 (2008) 1–171, [0706.3356].
  * [35] S. R. Klein and H. Mäntysaari, _Imaging the nucleus with high-energy photons_ , _Nature Rev. Phys._ 1 (2019) 662–674, [1910.10858].
  * [36] S. Klein and P. Steinberg, _Photonuclear and Two-photon Interactions at High-Energy Nuclear Colliders_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 323–354, [2005.01872].
  * [37] M. Arslandok et al., _Hot QCD White Paper_ , 2303.17254.
  * [38] H. Abramowicz and A. Caldwell, _HERA collider physics_ , _Rev. Mod. Phys._ 71 (1999) 1275–1410, [hep-ex/9903037].
  * [39] W. Buchmuller, M. F. McDermott and A. Hebecker, _Gluon radiation in diffractive electroproduction_ , _Nucl. Phys._ B487 (1997) 283–310, [hep-ph/9607290].
  * [40] Y. V. Kovchegov and L. D. McLerran, _Diffractive structure function in a quasi-classical approximation_ , _Phys. Rev._ D60 (1999) 054025, [hep-ph/9903246].
  * [41] H. Kowalski, T. Lappi and R. Venugopalan, _Nuclear enhancement of universal dynamics of high parton densities_ , _Phys. Rev. Lett._ 100 (2008) 022303, [0705.3047].
  * [42] H. Kowalski, T. Lappi, C. Marquet and R. Venugopalan, _Nuclear enhancement and suppression of diffractive structure functions at high energies_ , _Phys. Rev. C_ 78 (2008) 045201, [0805.4071].
  * [43] A. D. Le, A. H. Mueller and S. Munier, _Analytical asymptotics for hard diffraction_ , _Phys. Rev. D_ 104 (2021) 034026, [2103.10088].
  * [44] T. Lappi, A. D. Le and H. Mäntysaari, _Rapidity gap distribution of diffractive small- $x_{IP}$ events at HERA and at the EIC_, 2307.16486.
  * [45] L. Frankfurt, V. Guzey and M. Strikman, _Leading twist coherent diffraction on nuclei in deep inelastic scattering at small x and nuclear shadowing_ , _Phys. Lett. B_ 586 (2004) 41–52, [hep-ph/0308189].
  * [46] L. Frankfurt, V. Guzey and M. Strikman, _Leading Twist Nuclear Shadowing Phenomena in Hard Processes with Nuclei_ , _Phys. Rept._ 512 (2012) 255–393, [1106.2091].
  * [47] Y. V. Kovchegov and A. H. Mueller, _Gluon production in current nucleus and nucleon nucleus collisions in a quasi-classical approximation_ , _Nucl. Phys._ B529 (1998) 451–479, [hep-ph/9802440].
  * [48] H. Mäntysaari and J. Penttala, _Complete calculation of exclusive heavy vector meson production at next-to-leading order in the dipole picture_ , _JHEP_ 08 (2022) 247, [2204.14031].
  * [49] E. Iancu, K. Itakura and L. McLerran, _Geometric scaling above the saturation scale_ , _Nucl. Phys._ A708 (2002) 327–352, [hep-ph/0203137].
  * [50] A. H. Mueller and D. N. Triantafyllopoulos, _The energy dependence of the saturation momentum_ , _Nucl. Phys._ B640 (2002) 331–350, [hep-ph/0205167].
  * [51] S. Munier and R. Peschanski, _Traveling wave fronts and the transition to saturation_ , _Phys. Rev._ D69 (2004) 034008, [hep-ph/0310357].
  * [52] E. Levin and K. Tuchin, _Nonlinear evolution and saturation for heavy nuclei in DIS_ , _Nucl. Phys. A_ 693 (2001) 787–798, [hep-ph/0101275].
  * [53] ATLAS collaboration, _Photo-nuclear dijet production in ultra-peripheral Pb+Pb collisions_ , .
  * [54] ATLAS collaboration, G. Aad et al., _Two-particle azimuthal correlations in photonuclear ultraperipheral Pb+Pb collisions at 5.02 TeV with ATLAS_ , _Phys. Rev. C_ 104 (2021) 014903, [2101.10771].
  * [55] G. P. Lepage and S. J. Brodsky, _Exclusive processes in perturbative quantum chromodynamics_ , _Phys. Rev._ D22 (1980) 2157.
  * [56] M. G. Ryskin, _Diffractive J / psi electroproduction in LLA QCD_ , _Z. Phys. C_ 57 (1993) 89–92.
  * [57] S. J. Brodsky, L. Frankfurt, J. F. Gunion, A. H. Mueller and M. Strikman, _Diffractive leptoproduction of vector mesons in QCD_ , _Phys. Rev. D_ 50 (1994) 3134–3144, [hep-ph/9402283].
  * [58] J. D. Bjorken, J. B. Kogut and D. E. Soper, _Quantum Electrodynamics at Infinite Momentum: Scattering from an External Field_ , _Phys. Rev._ D3 (1971) 1382.
  * [59] N. N. Nikolaev and B. G. Zakharov, _Colour transparency and scaling properties of nuclear shadowing in deep inelastic scattering_ , _Z. Phys._ C49 (1991) 607–618.
  * [60] H. Kowalski and D. Teaney, _An impact parameter dipole saturation model_ , _Phys. Rev._ D68 (2003) 114005, [hep-ph/0304189].
  * [61] J. Nemchik, N. N. Nikolaev and B. G. Zakharov, _Scanning the BFKL pomeron in elastic production of vector mesons at HERA_ , _Phys. Lett. B_ 341 (1994) 228–237, [hep-ph/9405355].
  * [62] J. Nemchik, N. N. Nikolaev, E. Predazzi and B. G. Zakharov, _Color dipole phenomenology of diffractive electroproduction of light vector mesons at HERA_ , _Z. Phys. C_ 75 (1997) 71–87, [hep-ph/9605231].
  * [63] H. G. Dosch, T. Gousset, G. Kulzinger and H. J. Pirner, _Vector meson leptoproduction and nonperturbative gluon fluctuations in QCD_ , _Phys. Rev. D_ 55 (1997) 2602–2615, [hep-ph/9608203].
  * [64] J. R. Forshaw, R. Sandapen and G. Shaw, _Color dipoles and rho, phi electroproduction_ , _Phys. Rev. D_ 69 (2004) 094013, [hep-ph/0312172].
  * [65] H. Kowalski, L. Motyka and G. Watt, _Exclusive diffractive processes at HERA within the dipole picture_ , _Phys. Rev. D_ 74 (2006) 074016, [hep-ph/0606272].
  * [66] S. J. Brodsky, T. Huang and G. P. Lepage, _The Hadronic Wave Function in Quantum Chromodynamics_ , 6, 1980.
  * [67] A. H. Mueller, _Parton saturation at small x and in large nuclei_ , _Nucl. Phys._ B558 (1999) 285–303, [hep-ph/9904404].
  * [68] Y. V. Kovchegov and M. D. Sievert, _Calculating TMDs of a Large Nucleus: Quasi-Classical Approximation and Quantum Evolution_ , _Nucl. Phys._ B903 (2016) 164–203, [1505.01176].
  * [69] J. L. Albacete, N. Armesto, J. G. Milhano and C. A. Salgado, _Non-linear QCD meets data: A global analysis of lepton- proton scattering with running coupling BK evolution_ , _Phys. Rev._ D80 (2009) 034031, [0902.1112].
  * [70] A. M. Stasto, K. Golec-Biernat and J. Kwiecinski, _Geometric scaling for the total $\gamma^{*}p$ cross-section in the low x region_, _Phys. Rev. Lett._ 86 (2001) 596–599, [hep-ph/0007192].
  * [71] S. Munier and R. Peschanski, _Geometric scaling as traveling waves_ , _Phys. Rev. Lett._ 91 (2003) 232001, [hep-ph/0309177].
  * [72] E. Levin and K. Tuchin, _Solution to the evolution equation for high parton density qcd_ , _Nucl. Phys._ B573 (2000) 833–852, [hep-ph/9908317].
  * [73] L. V. Gribov, E. M. Levin and M. G. Ryskin, _Singlet structure function at small x: Unitarization of gluon ladders_ , _Nucl. Phys._ B188 (1981) 555–576.
  * [74] A. Kovner and U. A. Wiedemann, _Nonlinear qcd evolution: Saturation without unitarization_ , _Phys. Rev._ D66 (2002) 051502, [hep-ph/0112140].
  * [75] A. Kovner and U. A. Wiedemann, _Perturbative saturation and the soft pomeron_ , _Phys. Rev._ D66 (2002) 034031, [hep-ph/0204277].
  * [76] A. Kovner and U. A. Wiedemann, _No Froissart bound from gluon saturation_ , _Phys. Lett. B_ 551 (2003) 311–316, [hep-ph/0207335].
  * [77] K. J. Golec-Biernat and A. M. Stasto, _On solutions of the Balitsky-Kovchegov equation with impact parameter_ , _Nucl. Phys. B_ 668 (2003) 345–363, [hep-ph/0306279].
  * [78] J. Berger and A. Stasto, _Numerical solution of the nonlinear evolution equation at small x with impact parameter and beyond the LL approximation_ , _Phys. Rev. D_ 83 (2011) 034015, [1010.0671].
  * [79] J. Berger and A. M. Stasto, _Small x nonlinear evolution with impact parameter and the structure function data_ , _Phys. Rev. D_ 84 (2011) 094022, [1106.5740].
  * [80] I. Balitsky and G. A. Chirilli, _Next-to-leading order evolution of color dipoles_ , _Phys.Rev._ D77 (2008) 014019, [0710.4330].
  * [81] A. Kovner, M. Lublinsky and Y. Mulian, _Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov, Kovner evolution at next to leading order_ , _Phys. Rev._ D89 (2014) 061704, [1310.0378].
  * [82] A. V. Grabovsky, _Connected contribution to the kernel of the evolution equation for 3-quark Wilson loop operator_ , _JHEP_ 09 (2013) 141, [1307.5414].
  * [83] I. Balitsky and A. V. Grabovsky, _NLO evolution of 3-quark Wilson loop operator_ , _JHEP_ 01 (2015) 009, [1405.0443].
  * [84] R. Boussarie, A. V. Grabovsky, D. Y. Ivanov, L. Szymanowski and S. Wallon, _Next-to-Leading Order Computation of Exclusive Diffractive Light Vector Meson Production in a Saturation Framework_ , _Phys. Rev. Lett._ 119 (2017) 072002, [1612.08026].
  * [85] CMS collaboration, V. Khachatryan et al., _Coherent $J/\psi$ photoproduction in ultra-peripheral PbPb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV with the CMS experiment_, _Phys. Lett. B_ 772 (2017) 489–511, [1605.06966].
  * [86] ALICE collaboration, B. Abelev et al., _Coherent $J/\psi$ photoproduction in ultra-peripheral Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV_, _Phys. Lett. B_ 718 (2013) 1273–1283, [1209.3715].
  * [87] ALICE collaboration, S. Acharya et al., _Coherent photoproduction of $\rho^{0}$ vector mesons in ultra-peripheral Pb-Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}$ = 5.02 TeV_, _JHEP_ 06 (2020) 035 (2020), [2002.10897].
  * [88] ALICE collaboration, S. Acharya et al., _First measurement of coherent $\rho$0 photoproduction in ultra-peripheral Xe–Xe collisions at sNN=5.44 TeV_, _Phys. Lett. B_ 820 (2021) 136481, [2101.02581].
  * [89] ALICE collaboration, S. Acharya et al., _First measurement of the $|t|$-dependence of coherent $J/\psi$ photonuclear production_, _Phys. Lett. B_ 817 (2021) 136280, [2101.04623].
  * [90] ALICE collaboration, S. Acharya et al., _Coherent $J/\psi$ and $\psi^{\prime}$ photoproduction at midrapidity in ultra-peripheral Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}~{}=~{}5.02$ TeV_, _Eur. Phys. J. C_ 81 (2021) 712, [2101.04577].
  * [91] LHCb collaboration, R. Aaij et al., _$J/\psi$ photoproduction in Pb-Pb peripheral collisions at $\sqrt{s_{NN}}$= 5 TeV_, _Phys. Rev. C_ 105 (2022) L032201, [2108.02681].
  * [92] LHCb collaboration, _Study of coherent charmonium production in ultra-peripheral lead-lead collisions_ , 2206.08221.
  * [93] CMS collaboration, A. Tumasyan et al., _Probing small Bjorken- $x$ nuclear gluonic structure via coherent J/$\psi$ photoproduction in ultraperipheral PbPb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV_, 2303.16984.
  * [94] ALICE collaboration, _Physics of the ALICE Forward Calorimeter upgrade_ , .
  * [95] V. Guzey and M. Klasen, _How large is the diffractive contribution to inclusive dijet photoproduction in ultraperipheral collisions at the LHC?_ , _Phys. Rev. D_ 104 (2021) 114013, [2012.13277].