Search for a critical point of strongly-interacting matter in central 40Ar +
45Sc collisions at 13A–75A GeV/c beam momentum CERN-EP-2023-305 The critical
point of strongly interacting matter is searched for at the CERN SPS by the
NA61/SHINE experiment in central 40Ar + 45Sc collisions at 13A, 19A, 30A, 40A,
and 75A GeV/c. The dependence of the second-order scaled factorial moments of
proton multiplicity distributions on the number of subdivisions in transverse
momentum space is measured. The intermittency analysis uses statistically
independent data sets for every subdivision in transverse and cumulative-
transverse momentum variables. The results obtained do not indicate the
searched intermittent pattern. An upper limit on the fraction of correlated
protons and the intermittency index is obtained based on a comparison with the
Power-law Model.

##  1 Introduction

The experimental results are presented on intermittency analysis using second-
order scaled factorial moments of mid-rapidity protons produced in central
40Ar + 45Sc collisions at 13A, 19A, 30A, 40A, and 75A GeV/c beam momentum
($\sqrt{s_{\mathrm{NN}}}$ 111collision energy per nucleon pair in the center-
of-mass system $\approx$ 5.1–11.9 GeV). The measurements were performed by the
multi-purpose NA61/SHINE [1] apparatus being operated at the CERN Super Proton
Synchrotron (SPS). They are part of the strong interactions program of
NA61/SHINE devoted to studying the properties of strongly interacting matter,
such as the onset of deconfinement and critical endpoint (CP). Within this
program, a two-dimensional scan of collision energy and colliding nuclei size
was conducted [2].

In the QCD phase diagram, CP is a hypothetical endpoint of the first-order
phase transition line with properties of second-order phase transition. In the
proximity of CP, self-similar dynamics [3] is expected to lead to the
corresponding fluctuations of the chiral order parameter, belonging to the
3D-Ising universality class [4] and can be detected in transverse momentum
space within the framework of intermittency analysis of proton density
fluctuations using scaled factorial moments.

The idea of "intermittency", which was first introduced to the study of
turbulent flow [5], later became important in the physics of particle
production, especially as a way to study fluctuations. In the pioneering
article of Bialas and Peschanski [6] introducing intermittency analysis to
high-energy physics, it was proposed to study the scaled factorial moments of
the multiplicity of particles produced in high-energy collisions as a function
of the resolution size of rapidity interval.

This paper follows the proton intermittency analysis of 40Ar + 45Sc collisions
at 150A GeV/c [7] by the NA61/SHINE Collaboration. This analysis was performed
in intervals of transverse momentum and cumulative transverse momentum
distributions. Using the approach of intermittency analysis introduced in Ref.
[7], statistically independent data sets were used to obtain results for the
different numbers of intervals.

The paper is organized as follows. Section 2 introduces quantities exploited
for the CP search using the intermittency analysis. In Sec. 3, the
characteristics of the NA61/SHINE detector, relevant for the current study,
are briefly presented. The details of data selection and the analysis
procedure are presented in Sec. 4. Results obtained are shown in Sec. 5 and
compared with models in Sec. 6. A summary in Sec. 7 closes the paper.

Throughout this paper, the rapidity, $y=\text{atanh}\left(\beta_{L}\right)$,
is calculated in the collision center-of-mass frame by shifting rapidity in
laboratory frame by rapidity of the center-of-mass, assuming proton mass. The
$\beta_{L}=p_{L}/E$ ($c\equiv 1$) is the longitudinal (z) component of the
velocity, while $p_{L}$ and $E$ are particle longitudinal momentum and energy
in the collision center-of-mass system. The transverse component of the
momentum is denoted as $p_{T}=\sqrt{p_{x}^{2}+p_{y}^{2}}$, where $p_{x}$ and
$p_{y}$ are its horizontal and vertical components. The azimuthal angle is the
angle between the transverse momentum vector and the horizontal (x) axis.
Total momentum in the laboratory system is denoted as $p_{\textit{lab}}$.

The 40Ar + 45Sc collisions are selected by requiring a low energy value
measured by the forward calorimeter, the Projectile Spectator Detector (PSD).
This is the energy emitted into the region populated mostly by projectile
spectators. These collisions are called PSD-central collisions, and a
selection of collisions based on the PSD energy is called a PSD-centrality
selection.

##  2 Scaled factorial moments

###  2.1 Critical point and intermittency in heavy-ion collisions

A second-order phase transition leads to the divergence of the correlation
length ($\xi$). The infinite system becomes scale-invariant with the particle
density-density correlation function exhibiting power-law scaling, which
induces intermittent behavior of particle multiplicity fluctuations [8].

The maximum critical signal is expected when the freeze-out occurs close to
the CP. On the other hand, the energy density at the freeze-out is lower than
at the early stage of the collision. Clearly, the critical point should be
experimentally searched for in nuclear collisions at energies higher than that
of the onset of deconfinement – when quark-gluon plasma creation sets in.

The intermittent multiplicity fluctuations [6] were discussed as the signal of
CP by Satz [9], Antoniou et al. [10] and Bialas, Hwa [11]. This initiated
experimental studies of the structure of the phase transition region via
analyses of particle multiplicity fluctuations using scaled factorial moments
[12]. Later, additional measures of fluctuations were also proposed as probes
of the critical behavior [13, 14]. The NA61/SHINE experiment has performed a
systematic scan of collision energy and system size. To date none of the
anticipated signals have been observed for the critical behavior [15, 16, 17].
The new measurements may answer the question about the nature of the
transition region and, in particular, whether or not the critical point of
strongly interacting matter exists.

The scaled factorial moments $F_{r}(M)$ [6] of order $r$ are defined as:

$F_{r}(M)=\frac{\bigg{<}{\displaystyle{\frac{1}{M^{D}}\sum_{i=1}^{M^{D}}}n_{i}\;...\;(n_{i}-r+1)}\bigg{>}}{\bigg{<}{\displaystyle{\frac{1}{M^{D}}\sum_{i=1}^{M^{D}}}n_{i}}\bigg{>}^{r}}\>,$
(1)

Figure 1: Two-dimensional transverse momentum space is sub-divided into
$M\times M$ number of equally sized bins. The $n_{i}$ is the particle
multiplicity in a given sub-interval, $\Delta$ is the momentum region, and
$\delta$ is bin width (in this example $n_{i}$ equals 5 in the (4,3) sub-
interval).

where $M$ is the number of subdivision intervals in each of the $D$ dimensions
of the selected range $\Delta$ (see Fig. 1), $n_{i}$ is the particle
multiplicity in a given sub-interval and angle brackets denote averaging over
the analyzed events. In the presented analysis, $\Delta$ is divided into two-
dimensional ($D=2$) cells in $p_{x}$ and $p_{y}$. In case the mean particle
multiplicity, $\left<n_{i}\right>$, is proportional to the subdivision
interval size and for a Poissonian multiplicity distribution, $F_{r}(M)$ is
equal to 1 for all values of $r$ and $M^{D}$. This condition is satisfied in
the configuration space when the particle density is uniform. The momentum
distribution is, in general, non-uniform, and thus in the momentum space, it
is more convenient to use the cumulative variables [18] which, for a small
cell size, leave a power-law behavior unaffected and at the same time lead to
a uniformly distributed particle density.

If the system freezes out near CP, its properties are expected to be different
from those of an ideal gas [19]. Such a system represents a simple fractal and
$F_{r}(M)$ follows a power-law dependence:

$F_{r}(M)=F_{r}(\Delta)\cdot(M^{D})^{\varphi_{r}}\leavevmode\nobreak\ .$ (2)

Moreover, the exponent (intermittency index) $\varphi_{r}$ obeys the relation:

$\varphi_{r}=(r-1)\cdot(d_{r}/D)\>,$ (3)

where $d_{r}$, the anomalous fractal dimension, is independent of $r$ [11].
Such behavior is the analogue of critical opalescence in electromagnetically
interacting matter [3]. Importantly the critical properties given by Eqs. 2
and 3 are approximately preserved for a small cell size (large $M$) under
transformation to the cumulative variables [18, 20].

The ideal CP signal, Eqs. 2 and 3, derived for the infinite system in
equilibrium, are generally deteriorated by numerous experimental effects
present in high-energy nuclear collisions. This includes the system’s finite
size and evolution time, other dynamical correlations between particles, and
limited acceptance and finite resolution of measurements. Moreover, to
experimentally search for CP in high-energy collisions, the dimension and size
of the momentum interval must be chosen. Note that unbiased results can be
obtained only by analyzing variables and dimensions in which the singular
behavior appears [21, 22, 23]. Any other procedure is likely to distort the
critical-fluctuation signal.

Search for the CP via the study of proton fluctuations was suggested in
several publications [14, 24, 25, 26, 27, 28, 29, 30, 31]. In the case of a
pure system exhibiting critical fluctuations, for proton, $d=\varphi_{2}=5/6$
is expected [3].

###  2.2 Cumulative transformation

Scaled factorial moments are sensitive to the shape of the single-particle
momentum distribution. This dependence biases the signal of critical
fluctuations. To remove this dependence, one has two possibilities. The first
possibility is to exploit the mixed events, where each event is constructed
using particles from different experimental events, thereby removing all
possible correlations. The objective is to measure the following quantity:

$\Delta
F_{2}(M)=F_{2}^{\text{data}}(M)-F_{2}^{\text{mixed}}(M)\leavevmode\nobreak\ .$
(4)

It was already verified [12] that this procedure, to a large extent, removes
the dependence of $\Delta F_{r}(M)$ on the shape of single-particle
distribution, at least for $r=2$.

The second possibility is to use cumulative transformation [18], which for a
one-dimensional single-particle distribution $f(x^{\prime})$ reads:

$Q_{x}=\int\limits_{a}^{x}f(x^{\prime})dx^{\prime}\Bigg{/}\int\limits_{a}^{b}f(x^{\prime})dx^{\prime}\leavevmode\nobreak\
,$ (5)

where $a$ and $b$ are lower and upper limits of the variable $x^{\prime}$. For
a two-dimensional distribution $f(x^{\prime},y^{\prime})$ and a given
$x^{\prime}=x$ the transformation reads

$Q_{y}(x)=\int\limits_{a}^{y}f(x,y^{\prime})dy^{\prime}\Bigg{/}\int\limits_{a}^{b}f(x,y^{\prime})dy^{\prime}.$
(6)

After the cumulative transformation, any single-particle distribution becomes
flat, ranging from 0 to 1, see example distributions in Ref. [7], and
therefore, it removes the dependence on the shape of the single-particle
distribution for uncorrelated particles [18]. It also distorts all non-scale-
invariant correlations. On the other hand, the transformation is proven to
preserve the critical behaviour [20] given by Eq. 2, at least for the second-
order scaled factorial moments.

Both methods are approximate. Subtracting moments for mixed data set may
introduce negative $\Delta F_{2}(M)$ values [12] and using cumulative
quantities mixes the scales of the momentum differences and therefore may
distort eventual power-law behavior due to the finite size of the momentum
interval, $\Delta$.

##  3 The NA61/SHINE detector

The NA61/SHINE detector (see Fig. 2) is a large-acceptance hadron spectrometer
situated in the North Area H2 beam-line of the CERN SPS [1]. The main
components of the detection system used in the analysis are four large-volume
Time Projection Chambers (TPC). Two Vertex TPCs (VTPC-1/2) are located
downstream of the target inside superconducting magnets with a maximum
combined bending power of 9 Tm. The magnetic field was scaled in proportion to
the beam momentum in order to obtain similar $y-p_{T}$ acceptance at all beam
momenta. The main TPCs (MTPC-L/R) and two walls of pixel Time-of-Flight
(ToF-L/R) detectors are placed symmetrically on either side of the beamline
downstream of the magnets. The TPCs were filled with Ar:CO2 gas mixture in
proportions 90:10 for the VTPCs and 95:5 for the MTPCs. The Projectile
Spectator Detector (PSD), a zero-degree hadronic calorimeter, is positioned
20.5 m (16.7 m) downstream of the MTPCs at beam momenta of 75A (13A, 19A, 30A,
and 40A), centered in the transverse plane on the deflected position of the
beam.

Figure 2:  The schematic layout of the NA61/SHINE experiment at the CERN SPS
[1] showing the components used for the Ar+Sc energy scan (horizontal cut, not
to scale). The detector configuration upstream of the target is shown in the
inset. The alignment of the chosen coordinate system is shown on the plot; its
origin (x=y=z=0) lies in the middle of VTPC-2, on the beam axis. The nominal
beam direction is along the z-axis. Target is placed at z = -580.00 cm. The
magnetic field bends charged particle trajectories in the x–z (horizontal)
plane. The drift direction in the TPCs is along the (vertical) y-axis.

The PSD consists of 44 modules that cover a transverse area of almost 2.5 m2.
The central part of the PSD consists of 16 small modules with transverse
dimensions of 10 x 10 cm2 and its outer part consists of 28 large 20 x 20 cm2
modules. Moreover, a brass cylinder of 10 cm (30A –75A GeV/c) or 5 cm (19A
GeV/c) length and 5 cm diameter (degrader) was placed in front of the center
of the PSD to reduce electronic saturation effects and shower leakage from the
downstream side caused by the Ar beam and its heavy fragments. No degrader was
used at 13A GeV/c.

Primary beams of fully ionized 40Ar nuclei were extracted from the SPS
accelerator at beam momenta of 13A, 19A, 30A, 40A, and 75A GeV/c. Two
scintillation counters, S1 and S2, provide beam signal, and a veto counter V1,
with a 1 cm diameter hole, defines the beam before the target. The S1 counter
also provides the timing reference (start time for all counters). Beam
particles are selected by the trigger system requiring the coincidence T1 =
$\textrm{S1}\wedge\textrm{S2}\wedge\overline{\textrm{V1}}$. The three beam
position detectors (BPDs) placed upstream of the target [1] precisely measure
individual beam particle trajectories. Collimators in the beam line were
adjusted to obtain beam rates of $\approx 10^{4}$/s during the 10.4 s spill
and a super-cycle time of 32.4 s.

The target was a stack of 2.5 x 2.5 cm2 area and 1 mm thick 45Sc plates of 6
mm total thickness placed $\approx$ 80 cm upstream of VTPC-1. Impurities due
to other isotopes and elements were measured to be 0.3% [32]. No correction
was applied for this negligible contamination.

Interactions in the target are selected with the trigger system by requiring
an incoming 40Ar ion and a signal below that of beam ions from S5, a small 2
cm diameter scintillation counter placed on the beam trajectory behind the
MTPCs. This minimum bias trigger is based on the breakup of the beam ion due
to interactions in and downstream of the target. In addition, central
collisions were selected by requiring an energy signal below a set threshold
from the 16 central modules of the PSD, which measure mainly the energy
carried by projectile spectators. The cut was set to retain only the events
with the $\approx$ 30% smallest energies in the PSD. The event trigger
condition thus was T2 =
T1$\wedge\overline{\textrm{S5}}\wedge\overline{\textrm{PSD}}$.

##  4 Analysis

The goal of the analysis was to search for the critical point of the strongly
interacting matter by measuring the second-order scaled factorial moments for
a selection of protons produced in central 40Ar + 45Sc collisions at 13A–75A
GeV/c, using statistically independent points and cumulative variables.

###  4.1 Event selection

The NA61/SHINE detector recorded events using 13A–75A GeV/c 40Ar beam
impinging on a stationary 45Sc target. However, not all of those events
contain well-reconstructed central Ar+Sc interactions. Therefore the following
criteria were used to select data for further analysis:

  1. (i)

T2 trigger set to select central and semi-central collisions,

  2. (ii)

beam particle detected in at least three planes out of four of BPD-1 and BPD-2
and in both planes of BPD-3,

  3. (iii)

no off-time beam particle detected within a time window of $\pm$ 4 $\mu$s
around the trigger particle,

  4. (iv)

no interaction-event trigger detected within a time window of $\pm$ 25 $\mu$s
around the trigger particle,

  5. (v)

a high-precision interaction vertex with z position (fitted using the beam
trajectory and TPC tracks) no further than 2 cm away from the center of the Sc
target.

###  4.2 Centrality selection

  * ■

T2 trigger,

  * ■ \+ ■ \+ ■

30A, 40A and 75A,

  * ■ \+ ■ \+ ■ \+ ■

13A and 19A GeV/$c$.

Figure 3: Event centrality selection using the PSD energy. Modules used in the
centrality class determination were chosen based on the anti-correlation
between the measured energy and the track multiplicity in a given event. All
modules were used for $13$A and $19$A GeV/$c$. For $30$A, $40$A, and $75$A
GeV/$c$ 28 central modules were chosen. For the T2 trigger, 16 central modules
were used.

The final results presented in this paper refer to the 0–10% of 40Ar + 45Sc
collisions with the lowest energy value measured by a subset of PSD modules
(see Fig. 3). The selection of the modules was to optimize the sensitivity to
projectile spectators. For more details see Ref. [33]. Online event selection
by the hardware trigger (T2) used a threshold on the sum of electronic signals
from the 16 central modules of the PSD set to accept $\approx$ 30%, 35%, 30%,
35%, and 20% of the inelastic 40Ar + 45Sc collisions at 13A–75A GeV/c. To
select 0–10% central collision events, an upper limit of energy values
measured by a subset of modules are 143, 264, 446, 666, and 1290.6 GeV
corresponding to 40Ar + 45Sc collisions at 13A–75A GeV/c.

The event statistics after applying the selection criteria are summarized in
Table 1.

Table 1: The statistics of selected events for 40Ar + 45Sc collisions at beam momentum ($p_{beam}$) of 13A–75A GeV/c. $p_{beam}$ | number of events after cuts (106)
---|---
(GeV/c) | T2 trigger | beam quality | beam off-time | vertex z position | 0–10% of most central
13A | 2.14 | 1.60 | 1.56 | 1.48 | 0.50
19A | 2.51 | 2.00 | 1.93 | 1.83 | 0.52
30A | 3.71 | 2.93 | 2.85 | 2.74 | 0.91
40A | 5.71 | 4.87 | 4.74 | 4.53 | 1.29
75A | 2.89 | 2.44 | 2.37 | 2.32 | 1.16

###  4.3 Single-track selection

To select tracks of primary charged hadrons and to reduce the contamination by
particles from secondary interactions, weak decays, and off-time interactions,
the following track selection criteria were applied:

  1. (i)

track momentum fit including the interaction vertex is required to have
converged,

  2. (ii)

total number of reconstructed points on the track is required to be greater
than 30,

  3. (iii)

the sum of the number of reconstructed points in VTPC-1 and VTPC-2 is required
to be greater than 15,

  4. (iv)

the ratio of the number of reconstructed points to the number of potential
(maximum possible) points is required to be greater than 0.5 and less than
1.1,

  5. (v)

number of points used to calculate energy loss
($\textrm{d}E\\!/\\!\textrm{d}x$) is required to be greater than 30,

  6. (vi)

the distance between the track extrapolated to the interaction plane and the
vertex (track impact parameter) is required to be smaller than 4 cm in the
horizontal (bending) plane and 2 cm in the vertical (drift) plane.

As the analysis concerns mid-rapidity protons, only particles with center-of-
mass rapidity (assuming proton mass) between 0 and 0.75 were considered. Only
particles with transverse momentum components, $p_{x}$ and $p_{y}$, values
lower than 1.5 GeV/c were accepted for the analysis.

####  4.3.1 Proton selection

To identify proton candidates, positively charged particles were selected.
Their ionization energy loss (d$E$/d$x$) in TPCs is taken to be greater than
0.5 and less than the proton Bethe-Bloch value increased by the 15% difference
between the values for kaons and protons while the momentum is in the
relativistic-rise region (from 4 to 125 GeV/c). The energy loss versus the
logarithm of the total momentum of the selected positive particles for 40Ar +
45Sc collisions at 13A–75A GeV/c is shown in Fig. 4. The selected region is
marked with a magenta line.

Figure 4:  Energy loss versus the logarithm of the total momentum of
positively charged particles measured with the NA61/SHINE Time Projection
Chambers in the selected 40Ar + 45Sc events at 13A–75A GeV/c. The picture’s
dashed blue, black, and green lines are nominal Bethe-Bloch lines for protons,
kaons, and pions, respectively. The graphical cut selecting proton candidates
is marked with a magenta line.

This procedure was found to select, on average, about 60% of protons. The
remaining average kaon contamination is of the order of a few percent,
depending on collision energy, see Ref. [34].

The corresponding random proton losses do not bias the final results of the
independent production of protons in the transverse momentum space. The
results for correlated protons will be biased by the selection (see Sec. 5).
Thus, the random proton selection should be considered when calculating model
predictions (see Sec. 6).

###  4.4 Acceptance maps

####  4.4.1 Single-particle acceptance map

A three-dimensional ($p_{x}$, $p_{y}$ and center-of-mass rapidity) acceptance
maps [35] for central 40Ar + 45Sc collisions at 13A–75A GeV/c were created to
describe the momentum region selected for this analysis. These maps were
calculated by comparing the number of Monte Carlo-generated mid-rapidity
protons before and after detector simulation and reconstruction. Only
particles from the regions with at least 70% reconstructed particles are
analyzed. The single-particle acceptance maps have to be used for calculating
model predictions and they are given in Ref. [35].

####  4.4.2 Two-particle acceptance map

The Time Projection Chambers (the main tracking devices of NA61/SHINE) cannot
distinguish tracks that are too close to each other in space. At a small
distance, their clusters overlap, and signals are merged. Consequently, the
TPC cluster finder frequently rejects overlapping clusters, and the tracks can
be lost. Moreover, the TPC track reconstruction may fail to merge two track
fragments. This can generate split tracks out of a single track.

The mixed data set is constructed by randomly swapping particles from
different events so that each particle in each mixed event comes from
different recorded events.

Figure 5:  The two-track distance distribution of the ratio of the number of
track pairs in data and mixed events. The example for central 40Ar + 45Sc
collisions at 75A GeV/c is given. The left plot shows the distribution before
the mTTD cut and the right one after. See text for details.

For each pair of particles in both recorded and mixed events, a Two-Track
Distance (TTD) is calculated. It is the average distance of their tracks in
x–y plane at eight different z planes (-506, -255, -201, -171, -125, 125, 352,
and 742 cm). The TPC’s limitation in recognizing close tracks is illustrated
in Fig. 5 (left) for 40Ar + 45Sc collisions at 75A GeV/c. The data to mixed
ratio is significantly below one for TTD less than $\approx$ 2 cm. The TTD cut
values for each data set are listed in Table 2.

Calculating TTD requires knowledge of the NA61/SHINE detector geometry and
magnetic field. Hence, it is restricted to the Collaboration members. Here, a
momentum-based Two-Track Distance (mTTD) cut is introduced.

The mTTD cut removes the remaining split tracks from the data after the
potential point ratio cut (see Sec. 4.3), and it provides the precise
definition of the biased region in which we do not have good efficiency for
measuring two-tracks. Due to its momentum-based definition of the biased
region, the mTTD cut can be used for model comparison of the experimental
results without having access to the internal NA61/SHINE resources.

The magnetic field bends the trajectory of charged particles in the x–z plane.
Thus, it is most convenient to express the momentum of each particle in the
following momentum coordinates:

$\begin{split}s_{x}&=p_{x}/p_{xz}=\cos(\Psi)\leavevmode\nobreak\ ,\\\
s_{y}&=p_{y}/p_{xz}=\sin(\lambda)\leavevmode\nobreak\ ,\\\
\rho&=1/p_{xz}\leavevmode\nobreak\ ,\end{split}$ (7)

where $p_{xz}=\sqrt{p_{x}^{2}+p_{z}^{2}}$.

For each pair of particles, a difference in these coordinates is calculated
as:

$\begin{split}\Delta s_{x}&=s_{x,2}-s_{x,1}\leavevmode\nobreak\ ,\\\ \Delta
s_{y}&=s_{y,2}-s_{y,1}\leavevmode\nobreak\ ,\\\
\Delta\rho&=\rho_{2}-\rho_{1}\leavevmode\nobreak\ .\end{split}$ (8)

The distributions of particle pairs’ momentum difference for pairs with TTD
less than $\approx$ 2 cm (as an example for Ar+Sc collisions at 75A GeV/c) are
parameterized with ellipses in the new momentum coordinates. Such
parameterized elliptical cuts are defined as:

$\begin{split}\left(\frac{\Delta\rho}{r_{\rho}}\right)^{2}+\left(\frac{\Delta
s_{y}}{r_{s_{y}}}\right)^{2}&\leq 1\leavevmode\nobreak\ ,\\\
\left(\frac{\Delta s_{x}}{r_{s_{x}}}\right)^{2}+\left(\frac{\Delta
s_{y}}{r_{s_{y}}}\right)^{2}&\leq 1\leavevmode\nobreak\ ,\\\
\left(\frac{\Delta\rho\cos\theta-\Delta s_{x}\sin\theta}{r_{\rho
s_{x}}}\right)^{2}+\left(\frac{\Delta\rho\sin\theta+\Delta
s_{x}\cos\theta}{r_{s_{x}\rho}}\right)^{2}&\leq 1\leavevmode\nobreak\
,\end{split}$ (9)

where $r_{\rho s_{x}}$ and $r_{s_{x}\rho}$ is the semi-major and semi-minor
axis of an ellipse formed by $\Delta\rho$ and $\Delta s_{x}$, and $\theta$ is
the angle from the positive horizontal axis to the ellipse’s major axis.

Proton pairs with momenta difference inside the ellipses (Eqs. 9) are
rejected. The mTTD cut provides the precise definition of the biased region
(see Fig. 5 right). The mTTD cut can replace the TTD cut. The parameters of
the mTTD cut (see Eqs. 9) are given in Table 2, and the mTTD cut is used as a
two-particle acceptance map for the data analysis (see Sec. 5) and comparison
with the models (see Sec. 6).

Table 2: Numerical values of the mTTD cut parameters (see Eqs. 9) used to analyze central 40Ar + 45Sc collisions at 13A–75A GeV/c. Particle pairs with momenta inside the ellipses are rejected.  $p_{beam}$ (GeV/c) | TTD cut (cm) | $r_{\rho}(\mbox{Ge\kern-1.00006ptV\kern-1.49994pt/\kern-0.50003pt{c}})^{-1}$ | $r_{s_{y}}$ | $r_{s_{x}}$ | $r_{\rho s_{x}}(\mbox{Ge\kern-1.00006ptV\kern-1.49994pt/\kern-0.50003pt{c}})^{-1}$ | $r_{s_{x}\rho}(\mbox{Ge\kern-1.00006ptV\kern-1.49994pt/\kern-0.50003pt{c}})^{-1}$ | $\theta$
---|---|---|---|---|---|---|---
13A | 3.5 | 0.470 | 0.004 | 0.047 | 0.470 | 0.004 | 5
19A | 2.8 | 0.121 | 0.003 | 0.010 | 0.121 | 0.003 | 8
30A | 2.8 | 0.123 | 0.002 | 0.013 | 0.123 | 0.002 | 13
40A | 2.2 | 0.043 | 0.002 | 0.010 | 0.043 | 0.002 | 15
75A | 2.2 | 0.080 | 0.002 | 0.011 | 0.020 | 0.002 | 31

The effect of the mTTD cut for 40Ar + 45Sc collisions at 75A GeV/c mixed
events, and the Power-law Model (see Sec. 6.2) is shown in Fig. 6. The
dependence of $F_{2}(M)$ on $M^{2}$ in cumulative transverse momentum space
for $M^{2}>1$ is systematically below $F_{2}(M=1)$ when TTD or mTTD cut is
applied to fully uncorrelated mixed events (left) and the Power-law Model with
uncorrelated particles only (right).

Figure 6: Example of the TTD or mTTD cut impact on mixed events for central
40Ar + 45Sc collisions at 75A GeV/c (left) and on the Power-law Model [36]
with uncorrelated particles only (right) are shown. The blue circles
correspond to the dependence of $F_{2}(M)$ on $M^{2}$ where neither TTD nor
mTTD cut is applied to the mixed events or the Power-law Model. Green and red
points correspond to either TTD or mTTD cut applied.

###  4.5 Statistically-independent data points

The intermittency analysis gives the dependence of scaled factorial moments on
the number of subdivisions of transverse momentum or cumulative transverse
momentum intervals. In the past intermittency analyses, the same data set was
used to obtain results for each number of subdivisions. The results for
different $M^{2}$ were statistically correlated. Therefore, the full
covariance matrix was required for proper statistical treatment of the
results. This is numerically not trivial [37]. Here, the statistically-
independent data subsets were used to obtain results for different subdivision
numbers. In this case, the results for different subdivision numbers are
statistically independent. Only diagonal elements of the covariance matrix are
non-zero, and thus the complete relevant information needed to interpret the
results is easy to present graphically and use in the statistical tests. The
procedure decreases the number of events used to calculate each data point
increasing statistical uncertainties, and therefore forcing to reduce the
number of data points.

The number of events used in each subset was selected to obtain similar
statistical uncertainties for different subsets. Table 3 presents the
fractions of all selected events used to calculate each of the 10 points.

Table 3:  Fraction of the total number of selected events used to calculate second-order scaled factorial moments for the chosen number of cumulative momentum cells.  number of cells $M^{2}$ | $1^{2}$ | $50^{2}$ | $70^{2}$ | $86^{2}$ | $100^{2}$ | $111^{2}$ | $122^{2}$ | $132^{2}$ | $141^{2}$ | $150^{2}$
---|---|---|---|---|---|---|---|---|---|---
fraction of all events (%) | 0.5 | 3.0 | 5.0 | 7.0 | 9.0 | 11.0 | 13.0 | 15.5 | 17.0 | 19.0

###  4.6 Uncertainties and biases

The standard expression for the scaled factorial moments, Eq. 1 can be
rewritten as

$F_{2}(M)=2M^{2}\frac{\langle N_{2}(M)\rangle}{\langle N\rangle^{2}},$ (10)

where $N_{2}(M)$ denotes a total number of particle pairs in $M^{2}$ bins in
an event. Then the statistical uncertainties can be calculated using the
standard error propagation:

$\frac{\sigma_{F_{2}}}{|F_{2}|}=\sqrt{\frac{(\sigma_{N_{2}})^{2}}{\langle
N_{2}\rangle^{2}}+4\frac{(\sigma_{N})^{2}}{\langle
N\rangle^{2}}-4\frac{(\sigma_{N_{2}N})^{2}}{\langle N\rangle\langle
N_{2}\rangle}}.$ (11)

The final results presented in Sec. 5 are not corrected for possible biases.
Their magnitude was estimated by comparing results for pure Epos1.99 [38] and
Epos1.99 subjected to the detector simulation, reconstruction, and data-like
analysis. Figure 14 shows the comparison for 40Ar + 45Sc collisions at 13A–75A
GeV/c. Their differences are significantly smaller than the statistical
uncertainties of the experimental data and increase with $M^{2}$ up to about
0.1 at large $M^{2}$ values. Note that protons generated by Epos1.99 do not
show a significant correlation in the transverse momentum space, see Sec. 6.
In this case, the momentum resolution does not affect the results
significantly.

In the case of the critical correlations, the impact of the momentum
resolution may be significant, see Ref. [39] and Sec. 6 for detail. Thus a
comparison with models including short-range correlations in the transverse
momentum space requires smearing of momenta according to the experimental
resolution, which can be approximately parameterized as:

$\begin{split}p_{x}^{\text{smeared}}&=p_{x}^{\text{original}}+\delta
p_{x}\quad\text{and}\\\ p_{y}^{\text{smeared}}&=p_{y}^{\text{original}}+\delta
p_{y},\end{split}$ (12)

where $\delta p_{x}$ and $\delta p_{y}$ are randomly drawn from a Gaussian
distribution with $\sigma_{x}$ and $\sigma_{y}$ values for central 40Ar + 45Sc
collisions at 13A–75A GeV/c from the Table. 4.

Table 4: The $\sigma_{x}$ and $\sigma_{y}$ parameters, see text for details, used for smearing of proton transverse momentum in simulated 40Ar + 45Sc collisions at 13A–75A GeV/c.  $p_{beam}$ (GeV/c) |  $\sigma_{x}$ (MeV/c) |  $\sigma_{y}$ (MeV/c)
---|---|---
13A | 4.8 | 3.5
19A | 4.6 | 3.4
30A | 4.0 | 3.2
40A | 3.5 | 3.1
75A | 3.1 | 3.1

Uncertainties on the final results presented in Sec. 5 correspond to
statistical uncertainties.

##  5 Results

This section presents results on second-order scaled factorial moments (Eq. 1)
of $\approx 60\%$ randomly selected protons (losses due to proton
misidentification) with momentum smeared due to reconstruction resolution (Eq.
12) produced within the acceptance maps defined in Sec. 4.4 by strong and
electromagnetic processes in 0–10% central 40Ar + 45Sc interactions at 13A–75A
GeV/c. The results are shown as a function of the number of subdivisions in
transverse momentum space – the so-called intermittency analysis. The analysis
was performed for cumulative and original transverse momentum components.
Independent data sub-sets were used to calculate the results for each
subdivision.

Uncertainties correspond to statistical ones. Biases estimated using the
Epos1.99 [40] model (see Sec. 6) are significantly smaller than the
statistical uncertainties of the experimental data.

###  5.1 Two-particle correlation function

The two-particle correlation function, $\Delta p_{T}$, of selected proton
candidates within the analysis acceptance for central 40Ar + 45Sc collisions
at 13A–75A GeV/c is shown in Fig. 7. The correlation function is the ratio of
normalized $\Delta p_{T}$ distributions for data and mixed events. The data
distribution includes the mTTD cut, whereas the mixed one does not. The
decrease of the correlation function at $\Delta p_{T}\approx 0$ is due to
anti-correlation introduced by the mTTD cut.

Figure 7: Two-particle correlation functions in $\Delta p_{T}$ for selected
proton candidates within analysis acceptance for central 40Ar + 45Sc
collisions at 13A–75A GeV/c are shown. The data distribution includes the mTTD
cut, whereas the mixed one does not.

###  5.2 Subdivisions in cumulative transverse momentum space

Figures 8 and 9 present the dependence of the factorial moment on the number
of subdivisions in cumulative-transverse momentum space for the maximum
subdivision number of $M=150$ and $M=32$, respectively. The latter, coarse
subdivision, was introduced to limit the effect of experimental momentum
resolution; see Ref. [39] and below for details. The experimental results are
shown for 0–10% central 40Ar + 45Sc collisions at 13A–75A GeV/c. As a
reference, the corresponding results for mixed events are also shown.

By construction, the multiplicity distribution of protons in mixed events for
$M=1$ equals the corresponding distribution for the data. In mixed events, the
only correlation of particles in the transverse momentum space is due to the
mTTD cut. Both the data and mixed events include the mTTD cut. The mTTD cut is
necessary to properly account for the detector resolution – losses of close-
in-space tracks.

The results for subdivisions in cumulative transverse momentum space,
$F_{2}(M)$ for $M>1$ are systematically below $F_{2}(M=1)$. It is likely due
to the anti-correlation generated by the mTTD cut to the data (see Sec. 4.4.2
for details).

The experimental results show no increase of $F_{2}(M)$ with $M^{2}$. There is
no indication of the critical fluctuations for selected proton candidates.

Figure 8:  Results on the dependence of the scaled factorial moment of proton
multiplicity distribution on the number of subdivisions in cumulative
transverse momentum space $M^{2}$ for $1^{2}\leq M^{2}\leq 150^{2}$. Results
are shown for 0–10% central 40Ar + 45Sc collisions at 13A–75A GeV/c. Closed
red circles indicate the experimental data. Corresponding results for mixed
events (open triangles) are also shown. Both the data and mixed events include
the mTTD cut. Only statistical uncertainties are indicated.

Figure 9:  Results on the dependence of the scaled factorial moment of proton
multiplicity distribution on the number of subdivisions in cumulative
transverse momentum space $M^{2}$ for $1^{2}\leq M^{2}\leq 32^{2}$. Results
are shown for 0–10% central 40Ar + 45Sc collisions at 13A–75A GeV/c. Closed
red circles indicate the experimental data. Corresponding results for mixed
events (open triangles) are also shown. Both the data and mixed events include
the mTTD cut. Only statistical uncertainties are indicated.

###  5.3 Subdivisions in transverse momentum space

Figure 10 and 12 present the results which correspond to the results shown in
Figs. 8, and 9 but subdivisions are done in the transverse momentum space. By
construction, $F_{2}(1)$ values are equal for subdivisions in cumulative
transverse-momentum space and transverse-momentum space. But for the latter,
$F_{2}(M)$ rapidly increases from the value $F_{2}(M=1)$ to approximate
plateau at $M\approx 20$. This dependence is primarily due to the non-uniform
shape of the single-particle transverse momentum distributions, see Sec. 2.2.
It can be accounted for by comparing the experimental data results with those
obtained for the mixed events using $\Delta F_{2}(M)$. The dependence of
$\Delta F_{2}(M)$ on the number of sub-divisions, $M^{2}$ are shown in Figs.
11 and 13 for fine and coarse binning.

The experimental results presented in Figs. 10–13 do not show any significant
difference to the results for mixed events with the mTTD cut on $M^{2}$
($\Delta F_{2}(M)\approx 0$). There is no indication of the critical
fluctuations for selected protons.

Figure 10:  Results on the dependence of the scaled factorial moment of proton
multiplicity distribution on the number of subdivisions in transverse momentum
space $M^{2}$ for $1^{2}\leq M^{2}\leq 150^{2}$. Results are shown for 0–10%
central 40Ar + 45Sc collisions at 13A–75A GeV/c. Closed red circles indicate
the experimental data. Corresponding results for mixed events (open triangles)
are also shown. Both the data and mixed events include the mTTD cut. Only
statistical uncertainties are indicated.

Figure 11: Results on the dependence of $\Delta F_{2}(M)$ of proton
multiplicity distribution on the number of subdivisions in transverse momentum
space $M^{2}$ for $1^{2}\leq M^{2}\leq 150^{2}$. Results for 0–10% central
40Ar + 45Sc collisions at 13A–75A GeV/c are shown. Both the data and mixed
events include the mTTD cut. Only statistical uncertainties are indicated.

Figure 12:  Results on the dependence of the scaled factorial moment of proton
multiplicity distribution on the number of subdivisions in transverse momentum
space $M^{2}$ for $1^{2}\leq M^{2}\leq 32^{2}$. Results are shown for 0–10%
central 40Ar + 45Sc collisions at 13A–75A GeV/c. Closed red circles indicate
the experimental data. Corresponding results for mixed events (open triangles)
are also shown. Both the data and mixed events include the mTTD cut. Only
statistical uncertainties are indicated.

Figure 13: Results on the dependence of $\Delta F_{2}(M)$ of proton
multiplicity distribution on the number of subdivisions in transverse momentum
space $M^{2}$ for $1^{2}\leq M^{2}\leq 32^{2}$. Results for 0–10% central 40Ar
+ 45Sc collisions at 13A–75A GeV/c are shown. Both the data and mixed events
include the mTTD cut. Only statistical uncertainties are indicated.

##  6 Comparison with models

This section presents a comparison of the experimental results with two
models. The first one, Epos1.99 [40], takes into account numerous sources of
particle correlations, in particular, conservation laws and resonance decays,
but without critical fluctuations. The second one, the Power-law Model [41],
produces pairs of particles correlated by the power law inside the pair and
fully uncorrelated particles.

###  6.1 Epos

Minimum bias 40Ar + 45Sc events have been generated with Epos1.99 for
comparison. Signals from the NA61/SHINE detector were simulated with Geant3
software, and the generated events were reconstructed using the standard
NA61/SHINE procedure.

The number of forward spectators is used to select central Epos1.99 events of
40Ar + 45Sc collisions at 13A–75A GeV/c. To calculate model predictions (pure
Epos), 0–10% most central collisions were analyzed. Protons and proton pairs
within the single-particle and two-particle acceptance maps were selected.
Moreover, 60% of accepted protons were randomly selected for the analysis to
account for the limited proton identification in the experiment (see Fig. 4).

Results for the reconstructed Epos events were obtained as follows. The model
events were required to have the reconstructed primary vertex. Selected
protons and proton pairs (matching the generated particles used for
identification) were subjected to the same cuts as used for the experimental
data analysis (see Sec. 4). The results for the pure and the reconstructed
Epos events are compared in Fig. 14 for central 40Ar + 45Sc collisions at
13A–75A GeV/c. One concludes that for the Epos-like physics, the experimental
biases are smaller than the statistical uncertainties of the data.

Figure 14:  Results on the dependence of the scaled factorial moment of proton
multiplicity distribution on the number of subdivisions in cumulative
transverse momentum space, $M^{2}$ for $1^{2}\leq M^{2}\leq 150^{2}$, for
events generated with Epos1.99. Results are shown for 0–10% central 40Ar +
45Sc collisions at 13A–75A GeV/c. Closed black circles represent reconstructed
Epos, and open rectangles indicate pure (smeared) Epos. Both pure and
reconstructed Epos include the mTTD cut. Only statistical uncertainties are
indicated.

Finally, the experimental results are compared with the pure Epos predictions
for central 40Ar + 45Sc collisions at 13A–75A GeV/c and are shown in Fig. 15.
No significant differences are found.

Figure 15:  Results on the dependence of the scaled factorial moment of proton
multiplicity distribution on the number of subdivisions in cumulative
transverse momentum space $M^{2}$ for $1^{2}\leq M^{2}\leq 150^{2}$. Closed
red circles indicate the experimental data. Corresponding results for the
Epos1.99 model (open squares) were also shown for comparison. Results are
shown for 0–10% central 40Ar + 45Sc collisions at 13A–75A GeV/c. Results for
data and Epos were obtained using the mTTD cut. Only statistical uncertainties
are indicated.

###  6.2 Power-law Model

Inspired by expectations of the power-law correlations between particles near
the critical point, the Power-law Model [36] was developed to compare with the
experimental result. It generates momenta of uncorrelated and correlated
protons with a given single-particle transverse momentum distribution in
events with a given multiplicity distribution. The correlated protons are
produced in pairs, and their correlation inside the pair is given the power
law as described in Ref. [36]. The model has two controllable parameters:

  1. (i)

fraction of correlated particles, r,

  2. (ii)

strength of the correlation (the power-law exponent), $\phi$.

The transverse momentum of particles is drawn from the input transverse
momentum distribution. Correlated-particle pairs’ transverse momentum
difference follows a power-law distribution:

$\rho(|\Delta\overrightarrow{p_{T}}|)\sim|\Delta\overrightarrow{p_{T}}|^{-\phi},$
(13)

where the exponent $0\leq\phi<1$. For $r=0$, the Power-law Model results
correspond to the mixed event results. The exponent $\phi$ is related to the
intermittency index $\varphi_{2}$ as:

$\varphi_{2}=\frac{\phi+1}{2}.$ (14)

Azimuthal angle distribution is assumed to be uniform. The momentum component
along the beamline, $p_{z}$, is calculated assuming a uniform rapidity
distribution from $0$ to $0.75$ and proton mass.

Many high-statistics data sets have been produced using the Power-law Model
with multiplicity and inclusive transverse momentum distributions taken from
experimental data. Each data set has a different fraction of correlated
particles (varying from 0 to 2%) and a different power-law exponent (varying
from 0.00 to 0.95). The following effects have been included:

  1. (i)

Gaussian smearing of momentum components to mimic reconstruction resolution of
the momentum (see Eqs. 12),

  2. (ii)

single-particle acceptance map (see Sec. 4.4),

  3. (iii)

two-particle acceptance map (see Sec. 4.4.2),

  4. (iv)

random exchange of 40% of correlated particles with uncorrelated ones to
simulate 60% acceptance of protons (preserves the desired multiplicity
distribution but requires generating more correlated pairs at the beginning)
(see Sec. 4.3.1).

The influence of each of the above effects separately and all of them applied
together on $F_{2}(M)$ is shown in Fig. 16 for r = 0.02 and $\phi$ = 0.8, and
fine subdivisions.

Figure 16: Dependence of the scaled factorial moment on the number of
subdivisions in the cumulative transverse momentum for the Power-law Model
with the power-law exponent set to 0.80 and fraction of correlated particles
to 2%. Each line presents a result with a different effect included
separately, and the red circles all of them together for central 40Ar + 45Sc
collisions at 19A GeV/c (left) and 75A GeV/c (right).

Next, all generated data sets with all the above effects included have been
analyzed the same way as the experimental data. Obtained $F_{2}(M)$ results
have been compared with the corresponding experimental results and $\chi^{2}$
and a p-value were calculated. Statistical uncertainties from the model with
similar statistics to the data were used for the calculation. Examples of such
comparisons for central 40Ar + 45Sc collisions at 19A GeV/c (left) and 75A
GeV/c (right) are presented in Fig. 17.

Figure 17:  Examples of comparison of results for two Power-law Model data
sets with the experimental data for central 40Ar + 45Sc collisions at 19A
GeV/c (top) and 75A GeV/c (bottom). The left panel includes model predictions
assuming only uncorrelated protons, whereas the right one shows predictions
for 2.0% of correlated protons with power-law exponent $\phi=0.80$.

Figure 18 shows obtained exclusion plots as a function of the fraction of
correlated protons and the intermittency index (calculated from power-law
exponent) for central 40Ar + 45Sc collisions at 13A–75A GeV/c and 150A
GeV/c.222Note, that in the previously published paper on 40Ar + 45Sc at 150A
GeV/c [7], the same notation was used for intermittency index ($\varphi_{2}$)
and power-law exponent ($\phi$)

Figure 18:  Exclusion plots for the Power-law Model [36] parameters – the
fraction of correlated protons and the power-law exponent for central 40Ar +
45Sc collisions at 13A–75A GeV/c and 150A GeV/c. The white areas above the
line correspond to p-values less than 5%. The exclusion plots were obtained
using data (see Fig. 8) for the fine subdivisions ($1^{2}\leq M^{2}\leq
150^{2}$).

White areas above the line correspond to a p-value of less than 5% and may be
considered excluded (for this particular model). Fluctuations of the exclusion
lines in Fig. 18 are due to limited statistics of the experimental data.

Results for the coarse subdivision have low statistical uncertainties (see
Fig. 9), thus small deviations from the behavior expected for uncorrelated
particle production due to non-critical correlations (conservation laws,
resonance decays, quantum statistics, and others), as well as possible
experimental biases, may lead to a significant decrease of the p-values. Thus,
exclusion plots for the coarse subdivisions were not calculated.

The intermittency index $\varphi_{2}$ for an infinite system at QCD critical
point is expected to be $\varphi_{2}$ = 5/6 [3], assuming that the latter
belongs to the 3-D Ising universality class. If this value is set as the
power-law exponent of the Power-law Model with fine subdivisions (see Fig.
18), for the NA61/SHINE data on 0–10% central 40Ar + 45Sc collisions at
13A–75A GeV/c, an upper limit on the fraction of correlated protons is of the
order of 1%. It should be underlined that these numbers are specific to the
Power-law Model. Analyses made with other models, which could provide
different values depending on the assumed scenarios, remain beyond the scope
of this paper.

##  7 Summary

This paper reports on the search for the critical point of strongly
interacting matter in central 40Ar + 45Sc collisions at beam momenta of 13A,
19A, 30A, 40A, and 75A GeV/c. Results on second-order scaled factorial moments
of proton multiplicity distribution at mid-rapidity are presented. Protons
produced in strong and electromagnetic processes in 40Ar + 45Sc interactions
and selected by the single- and two-particle acceptance maps, as well as the
identification cuts, are used.

The scaled factorial moments are shown as a function of the number of
subdivisions of transverse momentum space – the so-called intermittency
analysis. The analysis was performed for cumulative and non-cumulative
transverse momentum components. Independent data sets were used to calculate
results for each subdivision. The influence of several experimental effects
was discussed and quantified. The results show no intermittency signal. A
summary of the proton intermittency from the Ar+Sc energy scan results is
shown in Fig. 19.

Figure 19: Summary of the proton intermittency results from the NA61/SHINE
Ar+Sc energy scan. Results on the dependence of the scaled factorial moment of
proton multiplicity distribution on the number of subdivisions in cumulative
transverse momentum space $M^{2}$ for $1^{2}\leq M^{2}\leq 150^{2}$ (left) and
$1^{2}\leq M^{2}\leq 32^{2}$ (right) are shown. The open circles represent
results on 0–20% central 40Ar + 45Sc collisions at 150A GeV/c [7]. Closed
circles indicate the experimental data results obtained within this work for
0–10% central 40Ar + 45Sc collisions at 13A, 19A, 30A, 40A, and 75A GeV/c.
Points for different energies are slightly shifted in horizontal axis to
increase readability.

The experimental data are consistent with the mixed events and the Epos model
predictions. An upper limit on the fraction of correlated protons of the order
of 1% was obtained based on a comparison with the Power-law Model.

Figure 20:  Diagram of chemical freeze-out temperature and baryon-chemical
potential. The dashed line indicates parameters in p+p interactions and the
dotted line in the central Pb+Pb collisions; points estimated and extrapolated
($\mu_{B}$ for p+p) based on Ref. [42]. The colored points mark reactions
(Ar+Sc and Pb+Pb [43]) in the $T-\mu_{B}$ phase diagram for which the search
for the critical point was conducted, and no evidence for the critical point
was found.

The ongoing critical point search studies via proton intermittency are
summarized on the diagram of chemical freeze-out temperature and chemical
potential (estimated based on Ref. [42]) and shown in Fig. 20. The
intermittency analysis of other reactions recorded within the NA61/SHINE
program on strong interactions is well advanced, and new results should be
expected soon.

## Acknowledgements

We would like to thank the CERN EP, BE, HSE and EN Departments for the strong
support of NA61/SHINE.

This work was supported by the Hungarian Scientific Research Fund (grant NKFIH
138136/138152), the Polish Ministry of Science and Higher Education
(DIR/WK/2016/2017/10-1, WUT ID-UB), the National Science Centre Poland (grants
2014/14/E/ST2/00018, 2016/21/D/ST2/01983, 2017/25/N/ST2/02575,
2018/29/N/ST2/02595, 2018/30/A/ST2/00226, 2018/31/G/ST2/03910,
2019/33/B/ST9/03059, 2020/39/O/ST2/00277), the Norwegian Financial Mechanism
2014–2021 (grant 2019/34/H/ST2/00585), the Polish Minister of Education and
Science (contract No. 2021/WK/10), the Russian Science Foundation (grant
17-72-20045), the Russian Academy of Science and the Russian Foundation for
Basic Research (grants 08-02-00018, 09-02-00664 and 12-02-91503-CERN), the
Russian Foundation for Basic Research (RFBR) funding within the research
project no. 18-02-40086, the Ministry of Science and Higher Education of the
Russian Federation, Project "Fundamental properties of elementary particles
and cosmology" No 0723-2020-0041, the European Union’s Horizon 2020 research
and innovation programme under grant agreement No. 871072, the Ministry of
Education, Culture, Sports, Science and Technology, Japan, Grant-in-Aid for
Scientific Research (grants 18071005, 19034011, 19740162, 20740160 and
20039012), the German Research Foundation DFG (grants GA 1480/8-1 and project
426579465), the Bulgarian Ministry of Education and Science within the
National Roadmap for Research Infrastructures 2020–2027, contract No.
D01-374/18.12.2020, Ministry of Education and Science of the Republic of
Serbia (grant OI171002), Swiss Nationalfonds Foundation (grant
200020117913/1), ETH Research Grant TH-01 07-3 and the Fermi National
Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of
Science, HEP User Facility managed by Fermi Research Alliance, LLC (FRA),
acting under Contract No. DE-AC02-07CH11359 and the IN2P3-CNRS (France).

The data used in this paper were collected before February 2022.

## References

  * [1] N. Abgrall et al., [NA61/SHINE Collab.] JINST 9 (2014) P06005, arXiv:1401.4699 [physics.ins-det].
  * [2] A. Aduszkiewicz et al., [NA61/SHINE Collab.], “Report from the NA61/SHINE experiment at the CERN SPS,” Tech. Rep. CERN-SPSC-2018-029. SPSC-SR-239, CERN, Geneva, Oct, 2018.  https://cds.cern.ch/record/2642286.
  * [3] N. G. Antoniou, F. K. Diakonos, A. S. Kapoyannis, and K. S. Kousouris Phys. Rev. Lett. 97 (2006) 032002, arXiv:hep-ph/0602051 [hep-ph].
  * [4] N. G. Antoniou, Y. F. Contoyiannis, and F. K. Diakonos Phys. Rev. E 62 (2000) 3125–3134, arXiv:hep-ph/0008047.
  * [5] B. B. Mandelbrot J. Fluid Mech. 62 no. 2, (1974) 331–358.
  * [6] A. Bialas and R. B. Peschanski Nucl. Phys. B273 (1986) 703–718.
  * [7] H. Adhikary et al., [NA61/SHINE Collab.] Eur. Phys. J. C 83 no. 9, (2023) 881, arXiv:2305.07557 [nucl-ex].
  * [8] J. Wosiek Acta Phys. Polon. B19 (1988) 863–866.
  * [9] H. Satz Nucl.Phys. B326 (1989) 613–618.
  * [10] N. G. Antoniou, E. N. Argyres, C. G. Papadopoulos, A. P. Contogouris, and S. D. P. Vlassopulos Phys. Lett. B245 (1990) 619–623.
  * [11] A. Bialas and R. Hwa Phys.Lett. B253 (1991) 436–438.
  * [12] T. Anticic et al., [NA49 Collab.] Eur. Phys. J. C 75 no. 12, (2015) 587, arXiv:1208.5292 [nucl-ex].
  * [13] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak Phys. Rev. Lett. 81 (1998) 4816–4819, arXiv:hep-ph/9806219 [hep-ph].
  * [14] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak Phys.Rev. D60 (1999) 114028, arXiv:hep-ph/9903292 [hep-ph].
  * [15] E. Andronov, [NA61/SHINE Collab.] Nucl. Phys. A 982 (2019) 835–838, arXiv:1807.10737 [nucl-ex].
  * [16] A. Aduszkiewicz et al., [NA61/SHINE Collab.] Eur. Phys. J. C 76 no. 11, (2016) 635, arXiv:1510.00163 [hep-ex].
  * [17] M. Maćkowiak-Pawłowska, [NA61/SHINE Collab.] EPJ Web Conf. 206 (2019) 03003.
  * [18] A. Bialas and M. Gazdzicki Phys. Lett. B 252 (1990) 483–486.
  * [19] F. Becattini, M. Gazdzicki, A. Keranen, J. Manninen, and R. Stock Phys. Rev. C 69 (2004) 024905, arXiv:hep-ph/0310049.
  * [20] S. Samanta, T. Czopowicz, and M. Gazdzicki Nucl. Phys. A 1015 (2021) 122299, arXiv:2105.00344 [nucl-th].
  * [21] A. Bialas and J. Seixas Phys. Lett. B250 (1990) 161–163.
  * [22] W. Ochs and J. Wosiek Phys. Lett. B214 (1988) 617–620.
  * [23] W. Ochs Phys. Lett. B247 (1990) 101–106.
  * [24] Y. Hatta and M. A. Stephanov Phys. Rev. Lett. 91 (2003) 102003, arXiv:hep-ph/0302002 [hep-ph]. [Erratum: Phys. Rev. Lett.91,129901(2003)].
  * [25] M. A. Stephanov Prog. Theor. Phys. Suppl. 153 (2004) 139–156, arXiv:hep-ph/0402115 [hep-ph]. [Int. J. Mod. Phys.A20,4387(2005)].
  * [26] K. Fukushima and T. Hatsuda Rept. Prog. Phys. 74 (2011) 014001, arXiv:1005.4814 [hep-ph].
  * [27] Y. Hatta and T. Ikeda Phys.Rev. D67 (2003) 014028, arXiv:hep-ph/0210284 [hep-ph].
  * [28] N. Antoniou, F. Diakonos, and A. Kapoyannis Phys.Rev. C81 (2010) 011901, arXiv:0809.0685 [hep-ph].
  * [29] F. Karsch and K. Redlich Phys. Lett. B695 (2011) 136–142, arXiv:1007.2581 [hep-ph].
  * [30] V. Skokov, B. Friman, and K. Redlich Phys. Rev. C83 (2011) 054904, arXiv:1008.4570 [hep-ph].
  * [31] K. Morita, V. Skokov, B. Friman, and K. Redlich Eur. Phys. J. C74 (2014) 2706, arXiv:1211.4703 [hep-ph].
  * [32] D. Banas, A. Kubala-Kukus, M. Rybczynski, I. Stabrawa, and G. Stefanek Eur. Phys. J. Plus 134 no. 1, (2019) 44, arXiv:1808.10377 [nucl-ex].
  * [33] A. Acharya et al., [NA61/SHINE Collab.] Eur. Phys. J. C 80 no. 10, (2020) 961, arXiv:2008.06277 [nucl-ex]. [Erratum: Eur.Phys.J.C 81, 144 (2021)].
  * [34] H. Adhikary et al. arXiv:2308.16683 [nucl-ex].
  * [35] NA61/SHINE acceptance map available at https://edms.cern.ch/document/2778197.
  * [36] T. Czopowicz arXiv:2309.13706 [hep-ex].
  * [37] B. Wosiek Acta Phys. Polon. B 21 (1990) 1021–1030.
  * [38] T. Pierog and K. Werner Phys. Rev. Lett. 101 (2008) 171101, arXiv:astro-ph/0611311 [astro-ph].
  * [39] S. Samanta, T. Czopowicz, and M. Gazdzicki J. Phys. G 48 no. 10, (2021) 105106, arXiv:2105.01763 [nucl-th].
  * [40] K. Werner Nucl. Phys. Proc. Suppl. 175-176 (2008) 81–87.
  * [41] T. Czopowicz, in preparation.
  * [42] F. Becattini, J. Manninen, and M. Gazdzicki Phys.Rev. C73 (2006) 044905, arXiv:hep-ph/0511092 [hep-ph].
  * [43] H. Adhikary, [NA61/SHINE Collab.] EPJ Web Conf. 274 (2022) 06008, arXiv:2211.10504 [hep-ex].

The NA61/SHINE Collaboration

H. Adhikary 14, P. Adrich 16, K.K. Allison 27, N. Amin 5, E.V. Andronov 23, T.
Antićić 3, I.-C. Arsene 13, M. Bajda 17, Y. Balkova 19, M. Baszczyk 18, D.
Battaglia 26, A. Bazgir 14, S. Bhosale 15, M. Bielewicz 16, A. Blondel 4, M.
Bogomilov 2, Y. Bondar 14, N. Bostan 26, A. Brandin 23, W. Bryliński 22, J.
Brzychczyk 17, M. Buryakov 23, A.F. Camino 29, P. Christakoglou 7, M. Ćirković
24, M. Csanád 9, J. Cybowska 22, T. Czopowicz 14, C. Dalmazzone 4, N. Davis
15, F. Diakonos 7, A. Dmitriev 23, P. von Doetinchem 28, W. Dominik 20, P.
Dorosz 18, J. Dumarchez 4, R. Engel 5, G.A. Feofilov 23, L. Fields 26, Z.
Fodor 8,21, M. Friend 10, M. Gaździcki 14,6, O. Golosov 23, V. Golovatyuk 23,
M. Golubeva 23, K. Grebieszkow 22, F. Guber 23, S.N. Igolkin 23, S. Ilieva 2,
A. Ivashkin 23, A. Izvestnyy 23, K. Kadija 3, A. Kapoyannis 7, N. Kargin 23,
N. Karpushkin 23, E. Kashirin 23, M. Kiełbowicz 15, V.A. Kireyeu 23, H.
Kitagawa 11, R. Kolesnikov 23, D. Kolev 2, Y. Koshio 11, V.N. Kovalenko 23, S.
Kowalski 19, B. Kozłowski 22, A. Krasnoperov 23, W. Kucewicz 18, M. Kuchowicz
21, M. Kuich 20, A. Kurepin 23, A. László 8, M. Lewicki 21, G. Lykasov 23,
V.V. Lyubushkin 23, M. Maćkowiak-Pawłowska 22, Z. Majka 17, A. Makhnev 23, B.
Maksiak 16, A.I. Malakhov 23, A. Marcinek 15, A.D. Marino 27, H.-J. Mathes 5,
T. Matulewicz 20, V. Matveev 23, G.L. Melkumov 23, A. Merzlaya 13, Ł. Mik 18,
A. Morawiec 17, S. Morozov 23, Y. Nagai 9, T. Nakadaira 10, M. Naskręt 21, S.
Nishimori 10, V. Ozvenchuk 15, A.D. Panagiotou 7, O. Panova 14, V. Paolone 29,
O. Petukhov 23, I. Pidhurskyi 14,6, R. Płaneta 17, P. Podlaski 20, B.A. Popov
23,4, B. Pórfy 8,9, M. Posiadała-Zezula 20, D.S. Prokhorova 23, D. Pszczel 16,
S. Puławski 19, J. Puzović 24†, R. Renfordt 19, L. Ren 27, V.Z. Reyna Ortiz
14, D. Röhrich 12, E. Rondio 16, M. Roth 5, Ł. Rozpłochowski 15, B.T.
Rumberger 27, M. Rumyantsev 23, A. Rustamov 1,6, M. Rybczynski 14, A. Rybicki
15, K. Sakashita 10, K. Schmidt 19, A.Yu. Seryakov 23, P. Seyboth 14, U.A.
Shah 14, Y. Shiraishi 11, A. Shukla 28, M. Słodkowski 22, P. Staszel 17, G.
Stefanek 14, J. Stepaniak 16, M. Strikhanov 23, H. Ströbele 6, T. Šuša 3, L.
Swiderski 16, J. Szewiński 16, R. Szukiewicz 21, A. Taranenko 23, A. Tefelska
22, D. Tefelski 22, V. Tereshchenko 23, A. Toia 6, R. Tsenov 2, L. Turko 21,
T.S. Tveter 13, M. Unger 5, M. Urbaniak 19, F.F. Valiev 23, M. Vassiliou 7, D.
Veberič 5, V.V. Vechernin 23, V. Volkov 23, A. Wickremasinghe 25, K. Wójcik
19, O. Wyszyński 14, A. Zaitsev 23, E.D. Zimmerman 27, A. Zviagina 23, and R.
Zwaska 25

† deceased

1 National Nuclear Research Center, Baku, Azerbaijan
2 Faculty of Physics, University of Sofia, Sofia, Bulgaria
3 Ruđer Bošković Institute, Zagreb, Croatia
4 LPNHE, University of Paris VI and VII, Paris, France
5 Karlsruhe Institute of Technology, Karlsruhe, Germany
6 University of Frankfurt, Frankfurt, Germany
7 University of Athens, Athens, Greece
8 Wigner Research Centre for Physics, Budapest, Hungary
9 Eötvös Loránd University, Budapest, Hungary
10 Institute for Particle and Nuclear Studies, Tsukuba, Japan
11 Okayama University, Japan
12 University of Bergen, Bergen, Norway
13 University of Oslo, Oslo, Norway
14 Jan Kochanowski University, Kielce, Poland
15 Institute of Nuclear Physics, Polish Academy of Sciences, Cracow, Poland
16 National Centre for Nuclear Research, Warsaw, Poland
17 Jagiellonian University, Cracow, Poland
18 AGH - University of Science and Technology, Cracow, Poland
19 University of Silesia, Katowice, Poland
20 University of Warsaw, Warsaw, Poland
21 University of Wrocław, Wrocław, Poland
22 Warsaw University of Technology, Warsaw, Poland
23 Affiliated with an institution covered by a cooperation agreement with CERN
24 University of Belgrade, Belgrade, Serbia
25 Fermilab, Batavia, USA
26 University of Notre Dame, Notre Dame, USA
27 University of Colorado, Boulder, USA
28 University of Hawaii at Manoa, Honolulu, USA
29 University of Pittsburgh, Pittsburgh, USA


# Quarkonia production in ultra-peripheral PbPb collisions at LHCb††thanks:
Presented at ”Diffraction and Low-$x$ 2022”, Corigliano Calabro (Italy),
September 24-30, 2022

Xiaolin Wang on behalf of the LHCb collaboration
Guangdong Provincial Key Laboratory of Nuclear Science, Guangdong-Hong Kong
Joint Laboratory of Quantum Matter, Institute of Quantum Matter, South China
Normal University, Guangzhou, China

###### Abstract

Measurements of coherent charmonium production cross sections together with
their ratio in ultra-peripheral PbPb collisions are studied at a nucleon-
nucleon centre-of-mass energy of $5.02\,\mathrm{TeV}$, the differential cross-
sections are measured as a function of rapidity and transverse momentum,
separately. The photo-production of ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$ mesons at low transverse momentum is studied in peripheral PbPb
collisions, which confirms coherent ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$ production in hadronic collisions. These latest results significantly
improve previous measurements and are compared with some theoretical
predictions.

##  1 Introduction

The LHCb detector is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$ [1, 2], which has great performance with
precise vertex reconstruction, high particle momentum resolution, and
excellent particle identification capability. Measurements of quarkonia
production in (ultra-)peripheral heavy-ion collisions play an important role
in studying the photon-nucleus interaction, the mechanism of vector meson
production, and the partonic structure of nuclei. Meanwhile, coherent photo-
production would provide an excellent laboratory to study the nuclear
shadowing effects and the initial state of heavy ion collisions at small-$x$
at the LHC [3].

##  2 Charmonia production in ultra-peripheral PbPb collisions at LHCb [4]

Ultra-peripheral collisions(UPCs) occur when the impact parameter is larger
than the sum of the radii of two incoming nuclei [5]. The two ions interact
via their cloud of semi-real photons, and photon-nuclear interactions
dominate. In UPCs, ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ and
$\psi{(2S)}$ mesons are produced from the colorless reaction of a photon from
one nucleus and a Pomeron from the other. If the photon interacts with a
Pomeron emitted by the whole nucleus, there would be no nucleus break-up, and
this is called coherent production, which is the process we are going to
study.

The charmonia candidates are reconstructed through the
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow{\mu^{+}\mu^{-}}$
and ${\psi{(2S)}}\rightarrow{\mu^{+}\mu^{-}}$ decay channels using a PbPb data
sample corresponding to an integrated luminosity of $228\pm 10\,\mathrm{\mu
b}^{-1}$, which was collected by the LHCb experiment in 2018. According to the
characteristics of UPCs, only events with low activity could be kept in the
selection. After that, the signal extraction is done through two steps. A fit
on the dimuon invariant mass spectrum is needed to estimate the non-resonant
background yields within the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
and $\psi{(2S)}$ mass windows, then the fits to
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ and $\psi{(2S)}$ transverse
momentum distributions of selected candidates are performed to determine the
coherent production events from the inclusive charmonia yields, as shown in
Fig. 1.

${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$

$\psi{(2S)}$

Figure 1: Fit to the invariant mass distribution of dimuon candidates (up) and
the $\ln(p_{\mathrm{T}}^{*2})$ distribution fit of dimuon candidates within
the $2.0<y^{*}<4.5$ range for ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
candidates (left) and $\psi{(2S)}$ candidates (right), where the starred
notation indicates that the observable is defined in the nucleus-nucleus
centre-of-mass frame.

The measured differential cross sections as a function of $y^{*}$ and
$p_{\mathrm{T}}^{*}$ of coherent ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$ and $\psi{(2S)}$ are shown in Fig. 2, where the starred notation
indicates that the observable is defined in the nucleus-nucleus centre-of-mass
frame. The differential cross-section ratio of coherent $\psi{(2S)}$ to
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ production is calculated as a
function of rapidity for the first time and shown in Fig. 3. Compared with
theoretical predictions, which are grouped as perturbative-QCD calculations
[6, 7] and colour-glass-condensate (CGC) models [8, 9, 10, 11, 12, 13, 14,
15], it could be found that the measurements are in agreement with most of the
prediction curves in general.

Figure 2: Differential cross-section as a function $y^{*}$ (up) and
$p_{\mathrm{T}}^{*}$ (down) for coherent
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ (left) and $\psi{(2S)}$
(right) production. Compared to theoretical predictions, which are grouped as
perturbative-QCD calculations (red lines) and colour-glass-condensate models
(blue lines).

Figure 3: Differential cross-section ratio of $\psi{(2S)}$ to
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ as a function of $y^{*}$,
compared to theoretical predictions, which are separated into perturbative QCD
calculations (red lines) and colour-glass-condensate models (blue lines).

##  3 Study of ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ photo-production
in PbPb peripheral collisions at $\sqrt{s_{NN}}=5\text{\,Te\kern-1.00006ptV}$
[16]

In peripheral collisions, the impact parameter is less than the sum of radii
of the two colliding nuclei, so there is not only photo-nuclear interaction
but also hadronic interaction. For hadronic productions, the gluon-gluon
interaction produces a pair of charm quarks and lastly goes into the
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ meson. In this case,
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ meson will have transverse
momentum typically between 1 and 2 $\text{\,Ge\kern-1.00006ptV\\!/}c$, and for
coherent photo-production, ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
would have a very small transverse momentum, less than 300
$\text{\,Me\kern-1.00006ptV\\!/}c$. Thanks to the high precision of the LHCb
experiment, there are two visible clear peaks in the distribution of dimuon
$\ln(p_{\mathrm{T}}^{*2})$, which could be used to discriminate photo-produced
and hadronically produced ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$, as
shown in Fig. 4.

Figure 4: Fit to the invariant mass spectrum of
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ candidates (left) and the
$\ln(p_{\mathrm{T}}^{2})$ distribution fit for
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ candidates (right).

The differential photo-production yields of
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ versus the rapidity, the number
of participants in the collision, and the double-differential
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ photo-production yields versus
transverse momentum as shown in Fig. 5. For the phenomenological models, one
scenario does not consider the destructive effect due to the overlap between
the two nuclei, whereas the other takes it into account. In general, the
trends of ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ photo-production
measurements and theoretical predictions are consistent. Meanwhile, the two
theoretical curves do not show a significant difference because the collisions
are peripheral and the nuclear overlapping effect is expected to increase in
more central collisions. Currently, the measurement is limited to low values
of $N_{part}$ due to detector limitation, and the measurement of
$p_{\mathrm{T}}$ is consistent with the coherent
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ photo-production; the peak
value of transverse momentum is around 60 $\text{\,Me\kern-1.00006ptV\\!/}c$.

Figure 5: The results of the differential photo-production yields of
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ versus the rapidity and the
number of participants in collisions. Double-differential yields of the
coherent ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ produced in peripheral
PbPb collisions as a function of $p_{\mathrm{T}}$ (right).

##  4 Conclusion

We report the results of coherent ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$, $\psi$(2S) production cross sections in UPCs, as well as their ratio.
It is the first coherent $\psi{(2S)}$ and $\psi{(2S)}$ to
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ ratio measurement in forward
rapidity region for UPC at LHC, and the differential cross section of coherent
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ and $\psi{(2S)}$ production in
PbPb UPC is also measured as a function of $p_{\mathrm{T}}^{*}$ for the first
time. Also, it is the most precise measurement of
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ production in UPCs. The
production of photo-produced ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
mesons in peripheral PbPb collisions is measured and by far it is the most
precise coherent photo-produced ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
by LHCb.

## References

  * [1] LHCb Collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [2] LHCb Collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [3] S. P. Jones, A. D. Martin, M. G. Ryskin, and T. Teubner, _Exclusive $J/\psi$ and $\Upsilon$ photoproduction and the low $x$ gluon_, J. Phys. G43 (2016) 035002, arXiv:1507.06942
  * [4] LHCb collaboration, _Study of coherent charmonium production in ultra-peripheral lead-lead collisions_ , arXiv:2206.08221
  * [5] C. A. Bertulani, S. R. Klein, and J. Nystrand, _Physics of ultra-peripheral nuclear collisions_ , Ann. Rev. Nucl. Part. Sci. 55 (2005) 271, arXiv:nucl-ex/0502005
  * [6] V. Guzey, E. Kryshen, and M. Zhalov, _Coherent photoproduction of vector mesons in ultraperipheral heavy ion collisions: Update for run2 at the CERN Large Hadron Collider_ , Phys. Rev. C93 (2016) 055206, arXiv:1602.01456
  * [7] V. Guzey, M. Strikman, and M. Zhalov, _Accessing transverse nucleon and gluon distributions in heavy nuclei using coherent vector meson photoproduction at high energies in ion ultraperipheral collisions_ , Phys. Rev. C95 (2017) 025204, arXiv:1611.05471
  * [8] V. P. Gonçalves and M. V. T. Machado, _Vector meson production in coherent hadronic interactions: an update on predictions for RHIC and LHC_ , Phys. Rev. C84 (2011) 011902, arXiv:1106.3036
  * [9] J. Cepila, J. G. Contreras, and M. Krelina, _Coherent and incoherent $\mathrm{J/}\psi$ photonuclear production in an energy-dependent hot-spot model_, Phys. Rev. C97 (2018) 024901, arXiv:1711.01855
  * [10] B. Z. Kopeliovich, M. Krelina, J. Nemchik, and I. K. Potashnikova, _Heavy quarkonium production in ultraperipheral nuclear collisions_ , arXiv:2008.05116
  * [11] V. P. Gonçalves et al., _Color dipole predictions for the exclusive vector meson photoproduction in $pp$, $p\mathrm{Pb}$, and $\mathrm{Pb}\mathrm{Pb}$ collisions at run2 LHC energies_, Phys. Rev. D96 (2017) 094027, arXiv:1710.10070
  * [12] V. P. Gonçalves and M. V. T. Machado, _The QCD pomeron in ultraperipheral heavy ion collisions: IV. Photonuclear production of vector mesons_ , Eur. Phys. J. C40 (2005) 519, arXiv:hep-ph/0501099
  * [13] H. Kowalski, L. Motyka, and G. Watt, _Exclusive diffractive processes at HERA within the dipole picture_ , Phys. Rev. D74 (2006) 074016, arXiv:hep-ph/0606272
  * [14] H. Mäntysaari and B. Schenke, _Probing subnucleon scale fluctuations in ultraperipheral heavy ion collisions_ , Phys. Lett. B772 (2017) 832, arXiv:1703.09256
  * [15] T. Lappi and H. Mäntysaari, _Diffractive vector meson production in ultraperipheral heavy ion collisions from the color glass condensate_ , PoS DIS2014 (2014) 069, arXiv:1406.2877
  * [16] 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, arXiv:2108.02681