#  Heavy quarkonium production
in ultraperipheral nuclear collisions

B. Z. Kopeliovich1<EMAIL_ADDRESS>M. Krelina2<EMAIL_ADDRESS>J. Nemchik2,3<EMAIL_ADDRESS>I. K. Potashnikova1<EMAIL_ADDRESS>1Departamento de Física, Universidad Técnica Federico Santa María, Avenida
España 1680, Valparaíso, Chile  2 Czech Technical University in Prague, FNSPE,
Břehová 7, 11519 Prague, Czech Republic  3 Institute of Experimental Physics
SAS, Watsonova 47, 04001 Košice, Slovakia

###### Abstract

Heavy quarkonium production in ultraperipheral nuclear collisions is described
within the QCD dipole formalism. Realistic quarkonium wave functions in the
rest frame are calculated solving the Schrödinger equation with a subsequent
Lorentz boost to high energy. We rely on several selected $Q\bar{Q}$
potentials, which provide the best description of quarkonium spectra and decay
widths, as well as data on diffractive electroproduction of quarkonia on
protons. Nuclear effects are calculated with the phenomenological dipole cross
sections fitted to DIS data. Higher twist effect related to the lowest
$Q\bar{Q}$ Fock component of the photon, as well as the leading twist effects,
related to higher components containing gluons, are included. The results for
coherent and incoherent photoproduction of charmonia and bottomonia on nuclei
are in a good accord with available data from the recent UPC measurements at
the LHC. They can also be verified in future experiments at the planned
electron-ion colliders.

###### pacs:

14.40.Pq,13.60.Le,13.60.-r

##  I Introduction

In comparison to hadro-production of heavy quarkonium states
[$V=J\\!/\\!\psi(1S)$, $\psi^{\,\prime}(2S)$, $\Upsilon(1S)$,
$\Upsilon^{\,\prime}(2S)$, $\Upsilon^{\,\prime\prime}(3S)$, …], where the
production mechanism is not well interpreted up to now, the main motivation
for their investigation in heavy-ion ultraperipheral collisions (UPC) at
Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC)
arises mainly from the following reasons:
– i) Here the mechanism of almost real photoproduction of heavy quarkonia is
well phenomenologically described Kopeliovich:1991pu ; Kopeliovich:1993pw ;
Nemchik:1996pp ; Nemchik:1996cw ; Nemchik:2002ug  and corresponding
theoretical uncertainties are rather well understood Krelina:2018hmt ;
Cepila:2019skb ; Krelina:2019egg  (a comprehensive review on quarkonium
phenomenology can be found in e.g. Refs. ivanov-04 ; brambilla-04 ;
brambilla-10 );
– ii) Relatively small size of heavy quarkonia should lead to minimization of
uncertainties inherent from the nonperturbative region, what allows to probe
and verify mainly a (semi)perturbative QCD (pQCD) origin of calculated nuclear
production amplitudes;
– iii) Besides the first data on charmonium production in UPC at RHIC reported
by the PHENIX collaboration Afanasiev:2009hy , during last few years a lot of
new data have been obtained at the LHC by the CMS Khachatryan:2016qhq , ALICE
Abelev:2012ba ; Abbas:2013oua ; Adam:2015sia  collaborations at c.m. collision
energy $\sqrt{s_{N}}=2.76\,\,\mbox{TeV}$, as well as by the LHCb LHCb:2018ofh
and ALICE Acharya:2019vlb  collaborations at even higher energy
$\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$. This can provide us with an additional
information about the dynamics of the corresponding photoproduction processes
and gives a further possibility to verify our ideas about the onset of various
nuclear effects occurring in interactions with nuclei in a wider energy range.

The UPC are associated with sufficiently large relative impact parameter of
colliding nuclei exceeding the sum of their radii. Then the electromagnetic
field generated by large nuclear charges induces an interaction of quasi–real
photons (where the photon virtuality $Q^{2}\to 0$) with nuclei followed by a
subsequent production of vector mesons (quarkonia). Such a quasi-real
diffractive electroproduction process is accompanied by two main phenomena
affecting the production rate.

The first phenomenon, representing the final state absorption of produced
quarkonia, is related to the effect known as the color transparency (CT).
Within the popular light-front (LF) color dipole approach Kopeliovich:1991pu ;
Kopeliovich:1993gk ; Kopeliovich:1993pw ; Nemchik:1994fp ; Nemchik:1994fq ;
Nemchik:1996cw ; Hufner:2000jb ; jan-00a ; jan-00b ; Kopeliovich:2001xj ;
Kopeliovich:2007wx , frequently used in the literature, such a diffractive
process can be described as a production of the $Q\bar{Q}$ pair by the real or
virtual photon. Treating the UPC, the transverse separation of a $Q\bar{Q}$
fluctuation is proportional to $1/m_{Q}$, where $m_{Q}$ is the quark mass. In
the photoproduction of heavy quarkonia off nuclei the large quark mass related
to the small sized $Q\bar{Q}$ pair should leads to a small absorption during
its propagation through the nucleus. This represents the basic idea of CT.
Here the large-sized photon configurations with larger absorption cross
sections are filtered out by the medium, what is known as the color filtering
phenomenon. The estimation of the corresponding length scale controlling the
evolution of the $Q\bar{Q}$ wave packet can be obtained in the rest frame of
the nucleus from the condition that the relative phase shift between the two
lowest levels, $V$ and $V^{\prime}$ becomes of the order of unity
Kopeliovich:1991pu ; Kopeliovich:2001xj

$l_{f}=\frac{2\,k}{\left.M_{V^{\prime}}\right.^{2}-\left.M_{V}\right.^{2}}\ ,$
(1)

where $k$ is the photon energy and $M_{V}$ and $M_{V^{\prime}}$ are the
quarkonium masses in $1S$ and $2S$ states, respectively. The scale $l_{f}$ in
Eq. (1) is usually called the formation length. In the present paper we assume
a very slow expansion of the $Q\bar{Q}$ transverse size relying on large
photon energies $k$, in particular in the LHC kinematic region. A more
sophisticated description of the $Q\bar{Q}$ evolution based on the rigorous LF
Green function formalism will be presented elsewhere.

Another important source of nuclear suppression is shadowing, which is caused
by a phase shift between amplitudes of photoproduction of $V$ on different
bound nucleons. When this phase shift reaches unity the interference becomes
destructive, i.e. the amplitudes tend to cancel each other. The corresponding
distance is called the coherence length (CL), Kopeliovich:1991pu ;
Kopeliovich:2001xj ,

$l_{c}=\frac{2\,k}{M_{V}^{2}}=\frac{s-M^{2}}{M\,M_{V}^{2}}\ ,$ (2)

where $M$ is the nucleon mass. It can be also interpreted as a lifetime of the
$Q\bar{Q}$ Fock component of the photon.

In our calculations of nuclear cross sections we have adopted for simplicity
an approximation of long-lived $Q\bar{Q}$ photon fluctuations widely used in
the literature. However, in order to include CL effects properly, we have
incorporated additionally corrections for the finite $l_{c}$ as is described
in Sect. II.4. The corresponding correction factors have been calculated up to
now within a simplified Glauber approximation Ivanov:2002kc ; Ivanov:2002eq ;
Ivanov:2007ms  whose applicability to quarkonium production, especially for
radially excited $2S$ states, is a questionable. For this reason, in the
present paper, the values of the finite-$l_{c}$ correction factors have been
obtained for the first time within the rigorous Green function formalism with
no restrictions for $l_{c}$.

Another initial state interaction phenomenon known to cause the nuclear
suppression in production of quarkonia is the gluon shadowing (GS). In the
present paper, similarly as for quarks discussed above, the GS is treated in
terms of the LF QCD dipole approach (see also Sect. II.5). It is related to
higher Fock components, $|Q\bar{Q}\,nG\rangle$, of the photon containing
gluons Mueller:1989st . In comparison with the lowest $|Q\bar{Q}\rangle$ Fock
state the corresponding lifetime of multi-gluon states is shorter
Kopeliovich:2000ra  due to their larger effective masses. This requires a
higher photon energy to reach a sufficient onset of the corresponding
shadowing corrections. The higher is the Fock state the higher must be the
photon energy to manifest its relevance for a shadowing.

Following analyses in Refs. Krelina:2018hmt ; Cepila:2019skb ; Krelina:2019egg
, in the present paper we treat a simple non-photon-like $V\rightarrow
Q\bar{Q}$ structure in order to avoid $D$-wave admixture excluding so the
additional non-specified contribution to the $S$-wave quarkonium
electroproduction cross sections. We have found in our present work
Krelina:2020bxt  that the undesirable relative contribution of $D$-wave
effects to nuclear cross sections of $J\\!/\\!\psi(1S)$ production in UPC does
not exceed $\sim 10\%-12\%$. However, for the $\psi^{\,\prime}(2S)$ state, it
plays an important negative role and causes $\sim 35\%$ enhancement of
corresponding production cross sections.

The advantage of the $S$-wave heavy quarkonia, studied in the present paper,
is based on a simple factorization of radial and spin-dependent components of
their wave functions. Here the former component is well defined in the
$Q\bar{Q}$ rest frame and can be obtained by solving the Schrödinger equation
for various realistic interaction potentials between $Q$ and $\bar{Q}$
proposed in the literature. In our calculations we choose two of them, the
power-like (POW) Martin:1980jx ; Barik:1980ai  and Buchmüller-Tye (BT)
Buchmuller:1980su  potentials, which provide the best description of available
data on charmonium electroproduction off protons as was analyzed in Ref.
Cepila:2019skb . The same reason is related to our choice of phenomenological
models for the dipole cross section, which represents the important ingredient
inherent from the LF color dipole formalism. Here we adopt two popular
parametrizations - Kopeliovich-Schafer-Tarasov (KST) from Ref.
Kopeliovich:1999am  and Golec-Biernat–Wusthoff (GBW) from Refs.
GolecBiernat:1998js ; GolecBiernat:1999qd .

Performing the boosting of the radial component of quarkonium wave functions
to the LF frame we rely on the widely used procedure Terentev:1976jk  for a
generation of the corresponding LF wave functions. In the spin-dependent
component, two-dimensional spinors describing the heavy $Q$ and $\bar{Q}$ in
the LF frame are related to those in the $Q\bar{Q}$ rest frame by the
transformation known as the Melosh spin rotation Melosh:1974cu ; Hufner:2000jb
; Krelina:2018hmt ; Cepila:2019skb ; Krelina:2019egg . The onset of spin
rotation effects for different LF wave functions, corresponding to various
realistic potential models for the $Q$-$\bar{Q}$ interaction, have been
analyzed in electroproduction of heavy quarkonia off protons in Refs.
Krelina:2018hmt ; Cepila:2019skb . However, in this work the study of UPC
requires the further extended investigation of spin effects in quarkonium
production off nuclei.

In this paper we extend our previous studies Ivanov:2002kc ; Ivanov:2002eq ;
Ivanov:2007ms  performing a comprehensive analysis of nuclear effects
occurring in heavy quarkonium production in UPC. Here we present for the first
time how the onset of GS and CL effects with no restrictions for $l_{c}$ is
manifested at different rapidities and c.m. collision energies corresponding
to RHIC and the LHC experiments. Besides, treating a simple non-photon-like
structure of the quarkonium vertex without $D$-wave admixture, we derive new
expressions for the nucleon electroproduction cross sections including the
Melosh spin transformation and containing an explicit dependence only on the
dipole cross section. This allows a simple straightforward generalization of
spin effects to nuclear targets, as well as to perform predictions for the
cross sections of coherent and incoherent photoproduction of various heavy
quarkonium states generated by different $Q-\bar{Q}$ interaction potentials in
the $Q\bar{Q}$ rest frame.

The paper is organized as follows. In the next Section we present expressions
for calculations of differential cross sections $d\sigma/dy$ corresponding to
coherent (elastic) and incoherent (quasi-elastic) heavy quarkonium production
in UPC as is described in Sections II.1 and II.2, respectively. Incorporation
of a small real part of the production amplitude, as well as spin rotation
effects is depicted in Sec. II.3. Consequently, in Sections II.4 and II.5 we
discuss two main effects modifying the nuclear suppression in UPC, the
corrections for a finite coherence length and the gluon shadowing. The former
is calculated for the first time within the LF dipole approach based on the
Green function formalism giving results that are substantially different from
values based on the standard vector dominance model (VDM). The next Sec. III
is devoted to comparison of model predictions with available data and to
analysis of particular nuclear effects in coherent and incoherent quarkonium
production in UPC. Finally, the last Sec. IV contains a summary with the main
concluding remarks about manifestations and the onset rate of various nuclear
effects occurring in UPC.

##  II Quarkonium production cross section in ultraperipheral collisions

The large charge $Z$ of heavy nuclei gives rise to strong electromagnetic
fields: in a heavy-ion UPC the photon field of one nucleus can produce a
photo-nuclear reaction in the other. Then the cross section for the
photoproduction of a vector meson $V$ by the Weizsäcker-Williams photons can
be written in the rest frame of the target nucleus $A$ as follows:

$k\frac{d\sigma}{dk}=\int\,d^{2}\tau\int\,d^{2}b\,\,n(k,\vec{b}-\vec{\tau},y)\,\frac{d^{2}\sigma_{A}(s,b,y)}{d^{2}b}~{}~{}+~{}~{}\Bigl{\\{}y\rightarrow-y\Bigr{\\}}\,,$
(3)

where the rapidity variable
$y=\ln\bigl{[}\frac{s}{M_{V}\sqrt{s_{N}}}\Bigr{]}\approx\ln\bigl{[}\frac{2kM+M^{2}}{M_{V}\sqrt{s_{N}}}\bigr{]}$
with the c.m. collision energy $\sqrt{s_{N}}$.

The formula (3) is derived in the one-photon-exchange approximation. Here the
variable $\vec{\tau}$ is the relative impact parameter of a nuclear collision
and $\vec{b}$ is the impact parameter of the photon-nucleon collision relative
to the center of one of the nuclei. Treating the collision of identical nuclei
with nuclear radius $R_{A}$ in UPC, the nuclei interact with each other
through their electromagnetic fields, what leads to a condition that the
impact parameter $\tau$ should be larger than the sum of the corresponding
nuclear radii, i.e. $\tau>2R_{A}$.

The variable $n(k,\vec{b})$ in Eq. (3) represents the photon flux induced by
the projectile nucleus with Lorenz factor $\gamma$ and has the following form,

$n(k,\vec{b})=\frac{\alpha_{em}Z^{2}k^{2}}{\pi^{2}\gamma^{2}}\Biggl{[}K_{1}^{2}\left(\frac{bk}{\gamma}\right)+\frac{1}{\gamma^{2}}K_{0}^{2}\left(\frac{bk}{\gamma}\right)\Biggr{]}\,,$
(4)

where $Z$ is the ion charge, $\alpha_{em}=1/137.036$ is the fine-structure
constant, $K_{0,1}$ are the modified Bessel functions and the Lorentz factor
$\gamma=2\gamma_{col}^{2}$ with $\gamma_{col}=\sqrt{s_{N}}/2M$. The first and
the second term in Eq. (4) corresponds to the flux of photons transversely and
longitudinally polarized to the ion direction, respectively. The former photon
flux dominates in ultra-relativistic collisions with $\gamma\gg 1$.
Consequently, in heavy-ion UPC at RHIC and the LHC one can safely neglect the
second term in Eq. (4) treating the photons as almost real due to a very small
virtuality, $-q^{2}=Q^{2}<1/R_{A}^{2}$.

###  II.1 Coherent production of quarkonia

To calculate the cross sections for coherent ($coh$) quarkonium production,
$\gamma A\to V\\!A$, we use the light-front dipole approach Kopeliovich:1991pu
, which has been applied to describe $J\\!/\\!\psi$ photoproduction off
nucleons Hufner:2000jb ; Krelina:2018hmt ; Cepila:2019skb  and nuclei
Ivanov:2002kc ; Nemchik:2002ug . In this approach, assuming sufficiently large
photon energies, corresponding to most of kinematic regions studied in the
present paper when the CL (2) $l_{c}\gg R_{A}$, the nuclear cross section
takes a simple asymptotic form,

$\displaystyle\frac{d^{2}\sigma_{A}^{coh}(s,b)}{d^{2}b}\Biggr{|}_{l_{c}\gg
R_{A}}$ $\displaystyle=$ $\displaystyle\Biggl{|}\int
d^{2}r\int_{0}^{1}d\alpha\,\Psi_{V}^{*}(\vec{r},\alpha)\,\Biggl{(}1-\exp\left[-\frac{1}{2}\,\sigma_{Q\bar{Q}}(r,s)\,T_{A}(b)\right]\Biggr{)}\Psi_{\gamma}(\vec{r},\alpha)\Biggr{|}^{2}$
(5) $\displaystyle\equiv$ $\displaystyle\Biggl{|}\int
d^{2}r\int_{0}^{1}d\alpha\,\Sigma_{A}^{coh}(r,\alpha,s,b)\,\Biggr{|}^{2}.$

We rely on the optical approximation, which is rather accurate for heavy
nuclei. Here $T_{A}(b)=\int_{-\infty}^{\infty}dz\,\rho_{A}(b,z)$ is the
nuclear thickness function normalized as $\int d^{2}b\,T_{A}(b)=A$, where
$\rho_{A}(b,z)$ is the nuclear density function of realistic Wood-Saxon form
with parameters taken from DeJager:1987qc ; $\Psi_{V}(r,\alpha)$ is the LF
wave function for heavy quarkonium; $\Psi_{\gamma}(r,\alpha)$ is the LF
distribution or the wave function of the $Q\bar{Q}$ Fock component of the
quasi-real (transversely polarized) photon, where the $Q\bar{Q}$ fluctuation
(dipole) has the transverse size $\vec{r}$ and the variable
$\alpha=p_{Q}^{+}/p_{\gamma}^{+}$ is the boost-invariant fraction of the
photon momentum carried by a heavy quark (or antiquark).

The universal dipole-nucleon total cross section $\sigma_{Q\bar{Q}}(r,s)$
depends on transverse dipole separation $r$ and c.m. energy squared
$s=M_{V}\,\sqrt{s_{N}}\,exp[y]$. Energy dependence of the dipole cross section
can be alternatively included also via variable
$x=M_{V}^{2}/W^{2}=M_{V}\,exp[-y]/\sqrt{s_{N}}$.

Notice that the coherent cross section, Eq. (5), is different from the usual
Glauber expression Bauer:1977iq  due to presence of the dipole cross section
zkl . It effectively includes the Gribov inelastic shadowing corrections
Gribov:1968jf ; Kopeliovich:2016jjx  in all orders for the $Q\bar{Q}$ Fock
component of the photon.

###  II.2 Incoherent production of quarkonia

Besides ”elastic” coherent photoproduction $\gamma A\to VA$, where the nucleus
remains intact, the vector meson can be produced in a quasi-elastic process
$\gamma A\to VA^{*}$, where the nucleus is excited and decays to fragments.
Important is that additional meson production is excluded. In this case one
can sum over different products of nuclear excitation and employ the
conditions of completeness. Of course one channel of elastic photoproduction
must be subtracted. It is instructive to see the result within the Glauber
approximation Hufner:1996dr ,

$\displaystyle\frac{d^{2}\sigma_{A}^{inc}(s,b)}{d^{2}b}\Bigr{|}_{l_{c}\gg
R_{A}}^{Gl}$ $\displaystyle\propto$
$\displaystyle\exp\left[-\sigma^{VN}_{in}\,T_{A}(b)\right]-\exp\left[-\sigma^{VN}_{tot}\,T_{A}(b)\right]$
(6) $\displaystyle=$
$\displaystyle\exp\left[-\sigma^{VN}_{tot}\,T_{A}(b)\right]\Bigl{\\{}\exp\left[-\sigma^{VN}_{el}\,T_{A}(b)\right]-1\Bigr{\\}}\,.$

Here the inelastic $V$-$N$ cross section
$\sigma^{VN}_{in}=\sigma^{VN}_{tot}-\sigma^{VN}_{el}$, where the elastic cross
section

$\displaystyle\sigma^{VN}_{el}\approx\frac{(\sigma^{VN}_{tot})^{2}}{16\,\pi\,B^{VN}}$
(7)

and $B^{VN}$ is the slope of the differential elastic $V-N$ cross section.

The cross section of incoherent photoproduction has the form, analogous to
(6), but with additional integrations over the dipole size (see derivation in
Sect. VII of Ref. kps ),

$\displaystyle\frac{d^{2}\sigma_{A}^{inc}(s,b)}{d^{2}b}\Bigr{|}_{l_{c}\gg
R_{A}}$ $\displaystyle=$ $\displaystyle\int
d^{2}r_{1}\int_{0}^{1}d\alpha_{1}\,\Psi_{V}^{*}(\vec{r}_{1},\alpha_{1})\,\Psi_{\gamma}(\vec{r}_{1},\alpha_{1})\,\exp\left[-\frac{1}{2}\sigma_{Q\bar{Q}}(r_{1},s)T_{A}(b)\right]$
(8) $\displaystyle\times$ $\displaystyle\int
d^{2}r_{2}\int_{0}^{1}d\alpha_{2}\,\Psi_{V}^{*}(\vec{r}_{2},\alpha_{2})\,\Psi_{\gamma}(\vec{r}_{2},\alpha_{2})\,\exp\left[-\frac{1}{2}\sigma_{Q\bar{Q}}(r_{2},s)\,T_{A}(b)\right]$
$\displaystyle\times$
$\displaystyle\left\\{\exp\left[\frac{\sigma_{Q\bar{Q}}(r_{1},s)\,\sigma_{Q\bar{Q}}(r_{2},s)}{16\pi
B(s)}\,T_{A}(b)\right]-1\right\\}\,.$

The elastic cross section of a heavy quarkonium on a nucleon is rather small
and the exponential in the last row of Eq. (8) can be expanded. Then we arrive
at a simple result Ivanov:2002kc ,

$\displaystyle\frac{d^{2}\sigma_{A}^{inc}(s,b)}{d^{2}b}\Biggr{|}_{l_{c}\gg
R_{A}}$ $\displaystyle\approx$ $\displaystyle\frac{T_{A}(b)}{16\pi
B(s)}\,\Biggl{|}\int
d^{2}r\int_{0}^{1}d\alpha\,\Psi_{V}^{*}(\vec{r},\alpha)\,\Psi_{\gamma}(\vec{r},\alpha)\,\sigma_{Q\bar{Q}}(r,s)$
(9) $\displaystyle\times$
$\displaystyle\exp\left[-\frac{1}{2}\sigma_{Q\bar{Q}}(r,s)\,T_{A}(b)\right]\,\Biggr{|}^{2}$
$\displaystyle\equiv$ $\displaystyle\frac{T_{A}(b)}{16\pi B(s)}\,\Biggl{|}\int
d^{2}r\int_{0}^{1}d\alpha\,\Sigma_{A}^{inc}(r,\alpha,s,b)\,\Biggr{|}^{2}.$

###  II.3 Real part of the production amplitude and spin rotation effects

The both Eqs. (5) and (9) contain a small correction due to the real part of
the $\gamma N\to VN$ amplitude applying the following replacement bronzan-74 ;
Nemchik:1996cw ; forshaw-03 ,

$\displaystyle\sigma_{Q\bar{Q}}(r,s)\Rightarrow\sigma_{Q\bar{Q}}(r,s)\,\left(1-i\,\frac{\pi}{2}\,\frac{\partial\,\ln\,{\sigma_{Q\bar{Q}}(r,s)}}{\partial\,\ln
s}\right)\ .$ (10)

Treating the non-photon-like structure of the $V\to Q\bar{Q}$ vertex
Ivanov:2002kc ; Ivanov:2002eq ; Ivanov:2007ms ; Krelina:2018hmt ;
Cepila:2019skb  the Melosh spin transformation is incorporated performing
following substitutions in Eqs. (5) and (9),

$\displaystyle\Sigma_{A}^{coh}(r,\alpha,s,b)\Rightarrow\Sigma_{A}^{coh}(r,\alpha,s,b)\cdot\Bigl{[}\Sigma^{(1)}(r,\alpha)+\Sigma^{(2)}(r,\alpha)\Bigr{]}$
$\displaystyle\Sigma_{A}^{inc}(r,\alpha,s,b)\Rightarrow\Sigma_{A}^{inc}(r,\alpha,s,b)\cdot\Bigl{[}\Sigma^{(1)}(r,\alpha)+\Sigma^{(2)}(r,\alpha)\Bigr{]},$
(11)

where

$\displaystyle\Sigma^{(1)}(r,\alpha)=N\,K_{0}(m_{Q}r)\int\limits_{0}^{\infty}dp_{T}\,p_{T}\,J_{0}(p_{T}r)\Psi_{V}(\alpha,p_{T})\left[\frac{2\,m_{Q}^{2}(m_{L}+m_{T})+m_{L}\,p_{T}^{2}}{m_{T}(m_{L}+m_{T})}\right]$
(12)

and

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\Sigma^{(2)}(r,\alpha)=N\,K_{1}(m_{Q}r)\int\limits_{0}^{\infty}dp_{T}\,p_{T}^{2}\,J_{1}(p_{T}r)\Psi_{V}(\alpha,p_{T})\left[\frac{m_{Q}^{2}(m_{L}+2m_{T})-m_{T}\,m_{L}^{2}}{m_{Q}\,m_{T}(m_{L}+m_{T})}\right].$
(13)

Here $N=Z_{Q}\,\sqrt{2N_{c}\,\alpha_{em}}/4\,\pi$, where the factor $N_{c}=3$
represents the number of colors in QCD, $Z_{Q}=\frac{2}{3}$ and $\frac{1}{3}$
are the charge-isospin factors for production of charmonia and bottomonia,
respectively, $J_{0,1}$ and $K_{0,1}$ are the Bessel functions of the first
kind and the modified Bessel functions of the second kind, respectively. The
variables $m_{T,L}$ in above formulas have the following form,

$m_{T}=\sqrt{m_{Q}^{2}+p_{T}^{2}}\,,\qquad
m_{L}=2\,m_{Q}\,\sqrt{\alpha(1-\alpha)}\,.$ (14)

In comparison with results from Ref. Ivanov:2002kc , the new form of Eqs.
(11)-(13) does not require to perform the so called “resummation procedure”
proposed in Ivanov:2002kc  in order to include properly the spin rotation
effects in nuclear photoproduction cross sections (5) and (9). However, it is
unclear whether these effects have been incorporated properly in Ref.
Henkels:2020kju  since the nucleon amplitude presented there depends also on a
derivative of $\sigma_{Q\bar{Q}}(r,s)$ and the corresponding explanation of
applicability of spin effects to nuclear targets is missing.

In our numerical calculations, following the results from Refs. Ivanov:2002kc
; Ivanov:2002eq ; Ivanov:2007ms ; Krelina:2018hmt ; Cepila:2019skb , we have
taken the charm and the bottom quark masses corresponding to values from
realistic phenomenological models for the $Q$-$\bar{Q}$ interaction potential,
such as POW and BT, used in our analysis. Similarly, the LF quarkonium wave
functions $\Psi_{V}(\alpha,r)$, after Lorentz boosting procedure from Ref.
Terentev:1976jk , are the counterparts of their nonrelativistic form generated
by these potentials in the $Q\bar{Q}$ rest frame.

###  II.4 Corrections for a finite coherence length

As was already mentioned above, the Green function approach allows to include
directly the effects of quantum coherence without any restrictions for the
coherence length, Eq. (2). However, as an alternative and a more simple way
instead of such a complicated method, one can use expressions (5) and (9)) for
nuclear cross sections in the limit of long CL, $l_{c}\gg R_{A}$, and then
provide additional corrections for a finite CL when $l_{c}\mathrel{\hbox
to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$<$}}R_{A}$. Such an incorporation of finite-$l_{c}$ effects via
the effective correction factors (form factors) $F(s,l_{c})$ has been
suggested in Ref. Hufner:1996jw  and employed in Ref. Ivanov:2002kc  for
calculations of charmonium photoproduction off nuclei. The corresponding
factors $F(s,l_{c})$ for the coherent and incoherent scattering have the
following form:

$\displaystyle F^{coh}(s,l_{c})$ $\displaystyle=$
$\displaystyle\int\\!d^{2}b\left|\,\int\limits_{-\infty}^{\infty}\\!dz\,\rho_{A}(b,z)\,F_{1}(s,b,z)\,e^{iz/l_{c}}\right|^{2}\bigl{/}\,\Bigl{(}...\Bigr{)}\Bigl{|}_{l_{c}\to\infty}\,,$
(15) $\displaystyle F^{inc}(s,l_{c})$ $\displaystyle=$
$\displaystyle\int\\!d^{2}b\,\int\limits_{-\infty}^{\infty}\\!dz\,\rho_{A}(b,z)\,\Bigl{|}F_{1}(s,b,z)-F_{2}(s,b,z,l_{c})\Bigr{|}^{2}\bigl{/}\,\Bigl{(}...\Bigr{)}\Bigl{|}_{l_{c}\to\infty}\,,$
(16) $\displaystyle F_{1}(s,b,z)$ $\displaystyle=$
$\displaystyle\,\exp\biggl{(}-\frac{1}{2}\,\sigma_{V\\!N}(s)\\!\\!\int\limits_{z}^{\infty}\\!dz^{\prime}\,\rho_{A}(b,z^{\prime})\biggr{)}\,,$
(17) $\displaystyle F_{2}(s,b,z,l_{c})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\sigma_{V\\!N}(s)\\!\\!\int\limits_{-\infty}^{z}\\!dz^{\prime}\,\rho_{A}(b,z^{\prime})\,F_{1}(s,b,z^{\prime})\,e^{-i(z-z^{\prime})/l_{c}}\,,$
(18)

where the variable $\sigma_{VN}(s)\equiv\sigma_{tot}^{VN}(s)$ represents the
quarkonium-nucleon total cross section introduced in the previous Sect. II.2.
Then the correction to the finite coherence length is usually incorporated
multiplying the nuclear cross sections in the limit of infinite coherence
length, Eqs. (5) and (9), by corresponding $l_{c}$-correction factors
$F^{coh}$ and $F^{inc}$ Ivanov:2002kc ,

$\displaystyle\frac{d^{2}\sigma_{A}^{coh}(s,b)}{d^{2}b}$ $\displaystyle=$
$\displaystyle\frac{d^{2}\sigma_{A}^{coh}(s,b)}{d^{2}b}\Biggl{|}_{l_{c}\gg
R_{A}}\cdot F^{coh}\\!\left(s,l_{c}(s)\right)\,,$ (19)
$\displaystyle\frac{d^{2}\sigma_{A}^{inc}(s,b)}{d^{2}b}$ $\displaystyle=$
$\displaystyle\frac{d^{2}\sigma_{A}^{inc}(s,b)}{d^{2}b}\Biggr{|}_{l_{c}\gg
R_{A}}\cdot F^{inc}\\!\left(s,l_{c}(s)\right)\,.$ (20)

The quarkonium-nucleon cross sections $\sigma_{VN}(s)$ in Eqs. (17) and (18)
have been calculated in Ref. Hufner:2000jb  for the $1S$ and $2S$ charmonium
states $J\\!/\\!\psi$ and $\psi^{\,\prime}$ including also the Melosh spin
effects. However, here we present the alternative and new expressions for
$\sigma_{VN}(s)$ based on Eqs. (12) and (13),

$\displaystyle\sigma_{VN}(s)=\frac{\int\limits_{0}^{1}d\alpha\,\int
d^{2}r\,\sigma_{Q\bar{Q}}(r,s)\Bigl{[}\Xi^{(1)}(r,\alpha)+\Xi^{(2)}(r,\alpha)\Bigr{]}^{2}}{\int\limits_{0}^{1}d\alpha\,\int
d^{2}r\,\Bigl{[}\Xi^{(1)}(r,\alpha)+\Xi^{(2)}(r,\alpha)\Bigr{]}^{2}}\,,$ (21)

where

$\displaystyle\Xi^{(1)}(r,\alpha)=\int\limits_{0}^{\infty}dp_{T}\,p_{T}\,J_{0}(p_{T}r)\Psi_{V}(\alpha,p_{T})\left[\frac{2\,m_{Q}^{2}(m_{L}+m_{T})+m_{L}\,p_{T}^{2}}{m_{Q}\,m_{T}(m_{L}+m_{T})}\right]\,,$
(22)
$\displaystyle\Xi^{(2)}(r,\alpha)=\int\limits_{0}^{\infty}dp_{T}\,p_{T}\,J_{1}(p_{T}r)\Psi_{V}(\alpha,p_{T})\left[\,p_{T}\,\frac{m_{Q}^{2}(m_{L}+2m_{T})-m_{T}\,m_{L}^{2}}{m_{Q}^{2}\,m_{T}(m_{L}+m_{T})}\right]\,.$
(23)

Figure 1:  (left panel) - Comparison of predictions for the quarkonium-nucleon
total cross sections from the exact expression (21) (solid lines) with results
where spin rotation effects are neglected (dashed lines). Calculations are
performed using quarkonium wave functions determined by the BT potential, as
well as the KST model for the dipole cross section. (right panel) - the same
as the left panel but with the GBW model for the dipole cross section.

As an example, using quarkonium wave functions generated by the BT potential
and adopting KST and GBW models for the dipole cross section, we plot in Fig.
1 by the solid lines the quarkonium-nucleon total cross sections for various
quarkonium states. Our calculations are based on Eq. (21), which includes the
spin rotation effects. Here the onset of such effects is tested by comparing
with results without any spin rotation depicted by the dashed lines. One can
see that the Melosh spin transformation manifests itself only for the
$J\\!/\\!\psi$ and $\Upsilon$ production, enhancing the corresponding
quarkonium-nucleon cross section by about $10\%-20\%$ and $8\%-10\%$,
respectively. For other quarkonium states the onset of spin effects is very
weak and can be neglected.

|  $\sigma_{0}$ (mb): KST |  $\Delta$: KST |  $\sigma_{0}$ (mb): GBW |  $\Delta$: GBW
---|---|---|---|---
$~{}~{}~{}J\\!/\\!\psi(1S)$ |  $6.04$ |  $0.206$ |  $4.78$ | $0.195$
$~{}~{}~{}\psi^{\,\prime}(2S)$ | $16.60$ |  $0.122$ | $11.28$ | $0.114$
$~{}~{}~{}\Upsilon(1S)$ |  $1.76$ |  $0.265$ |  $0.87$ | $0.250$
$~{}~{}~{}\Upsilon^{\,\prime}(2S)$ |  $7.75$ |  $0.188$ |  $3.79$ | $0.187$
Table 1:  Values of parameters $\sigma_{0}$ and $\Delta$ for the
parametrization Eq. (24). The quarkonium-nucleon cross sections for various
quarkonium states are calculated according to Eqs. (21)-(23) using wave
functions determined from the BT potential, as well as the KST and GBW
parametrizations for the dipole cross section.

The quarkonium-nucleon total cross sections $\sigma_{VN}(s)$ calculated from
Eq. (21) can be also parametrized as suggested in Ref. Hufner:2000jb . In Tab.
1 we present the values of parameters $\sigma_{0}$ and $\Delta$ inherent from
the following parametrization of $\sigma_{VN}(s)$,

$\sigma_{VN}(s)=\sigma_{0}\,\cdot\Bigl{(}\frac{s}{s_{0}}\Bigr{)}^{\Delta}\,,\qquad
s_{0}=1000\,\,\mbox{GeV}^{2}\,,$ (24)

where the parameters are averaged over the energy interval
$10^{2}\,\,\mbox{GeV}^{2}<s<10^{5}\,\,\mbox{GeV}^{2}$ and
$10^{3}\,\,\mbox{GeV}^{2}<s<10^{6}\,\,\mbox{GeV}^{2}$ for charmonium and
bottomonium states, respectively. Differences between values of these
parameters related to KST and GBW dipole models can be treated as a measure of
the theoretical uncertainty in determination of their magnitudes.

Figure 2:  (left panels) - comparison of $l_{c}$-correction factors based on
VDM and obtained from Eqs. (15)-(18) (dashed lines) with more proper
calculations within the color dipole approach based on the Green function
technique, Eqs. (25) and (26) (solid lines). The quarkonium wave functions are
generated by the BT potential. These factors correspond to coherent
respectively incoherent production of $J\\!/\\!\psi$ (thick lines) and
$\psi^{\,\prime}$ (thin lines) in $Au-Au$ UPC. (right panels) - the same as
left panels but for production of $\Upsilon$ and $\Upsilon^{\,\prime}$. Top
and bottom panels correspond to calculations using KST and GBW models for the
dipole cross section, respectively.

From known values of $\sigma_{VN}(s)$ one can calculate the $l_{c}$-correction
factors $F(s)$ using Eqs. (15)-(18). The corresponding values of $F(s)$ for
coherent and incoherent $J\\!/\\!\psi$ (thick lines) and $\psi^{\,\prime}$
(thin lines) photoproduction on the gold target are depicted by dashed lines
on left panels of Fig. 2 as function of the square of c.m. energy $s$. Here
the top and bottom panel corresponds to calculations using KST and GBW
parametrization for the dipole cross section, respectively. The analogous
results for $F(s)$ but for photoproduction of $\Upsilon$ and
$\Upsilon^{\,\prime}$ are depicted on right panels of the same Fig. 2.

Here we would like to emphasize that the above simple estimation of
corrections for a finite coherence length within VDM can be applied with a
reasonable accuracy only for photoproduction of $\rho$ mesons off nuclei (for
the review of photo- and electroproduction of vector mesons and VDM see Ref.
Bauer:1977iq , for example). However, for the heavy quarkonium production
(especially for production of radially excited $2S$ states) it cannot provide
reliable results (see estimations in Ref. Hufner:1997jg ). For this reason, it
is worth switching to a more rigorous quantum-mechanical description of
factors $F^{coh}$ and $F^{inc}$ within the sophisticated Green function
formalism. Then the corresponding functions $F_{1}$ and $F_{2}$ in Eqs. (15)
and (16) have the following form,

$\displaystyle F_{1}(s,b,z)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{1}d\alpha\int
d^{2}r_{1}\,d^{2}r_{2}\,\Psi^{*}_{V}(\vec{r}_{2},\alpha)\,G_{q\bar{q}}(z^{\prime},\vec{r}_{2};z,\vec{r}_{1})\,\sigma_{Q\bar{Q}}(r_{1},s)\,\Psi_{\gamma}(\vec{r}_{1},\alpha)\Bigl{|}_{z^{\prime}\to\infty}$
(25) $\displaystyle F_{2}(s,b,z)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\int\limits_{-\infty}^{z}dz_{1}\,\rho_{A}(b,z_{1})\,\int\limits_{0}^{1}d\alpha\int
d^{2}r_{1}\,d^{2}r_{2}\,d^{2}r\,\Psi^{*}_{V}(\vec{r}_{2},\alpha)$ (26)
$\displaystyle\times$ $\displaystyle
G_{q\bar{q}}(z^{\prime}\to\infty,\vec{r}_{2};z,\vec{r})\,\sigma_{Q\bar{Q}}(\vec{r},s)\,G_{q\bar{q}}(z,\vec{r};z_{1},\vec{r}_{1})\,\sigma_{Q\bar{Q}}(\vec{r}_{1},s)\,\Psi_{\gamma}(\vec{r}_{1},\alpha)\,,$

where $\Psi_{\gamma}(\vec{r}_{1},\alpha)=K_{0}(m_{Q}r)$ and the Green function
$G_{q\bar{q}}(z^{\prime},\vec{r}_{2};z,\vec{r}_{1})$ describes the propagation
of an interacting $Q\bar{Q}$ pair in a nuclear medium between points with
longitudinal coordinates $z$ and $z^{\prime}$ and with initial and final
separations $\vec{r}_{1}$ and $\vec{r}_{2}$. Here, for simplicity, we took the
quadratic form for the dipole cross section
$\sigma_{Q\bar{Q}}(r,s)=C(s)\,r^{2}$ and for the LF $Q\bar{Q}$ interaction
potential in the evolution equation for the Green function (see Ref.
Kopeliovich:2001xj , for example). Another simplification is related to a
constant nuclear density
$\rho_{A}(b,z)=\rho_{0}\Theta(R_{A}^{2}-b^{2}-z^{2})$, which can be safely
adopted for heavy nuclei used in our analysis. For the LF quarkonium wave
functions we used the following Gaussian shape Kopeliovich:2001xj  for the
$1S$ and $2S$ states,

$\displaystyle\Psi_{V}(r,\alpha)$ $\displaystyle=$ $\displaystyle
C_{V}\,a^{2}(\alpha)\,f(\alpha)\,\exp\Bigl{[}-\frac{1}{2}a^{2}(\alpha)r^{2}\Bigr{]}$
(27) $\displaystyle\Psi_{V^{\prime}}(r,\alpha)$ $\displaystyle=$
$\displaystyle
C_{V^{\prime}}\,a^{2}(\alpha)\,f(\alpha)\,\exp\Bigl{[}-\frac{1}{2}a^{2}(\alpha)r^{2}\Bigr{]}\Bigl{\\{}1+4h(\alpha)-\beta\,2a^{2}(\alpha)r^{2}\Bigr{\\}}\,,$
(28)

where

$f(\alpha)=\exp\Bigl{[}-h(\alpha)\Bigr{]}=\exp\Bigl{[}-\frac{m_{Q}^{2}}{2a^{2}(\alpha)}+\frac{4\alpha(1-\alpha)m_{Q}^{2}}{2a^{2}(\alpha)}\Bigr{]}\,,$
(29)

and the parameter $\beta$, controlling the position of the node, has been
determined from the orthogonality condition $\int
d^{2}r\,d\alpha\Psi_{V}(r,\alpha)\,\Psi_{V^{\prime}}(r,\alpha)=\delta_{VV^{\prime}}$.
We have found $\beta=0.908$ and $0.963$ for production of $\psi^{\,\prime}$
and $\Upsilon^{\,\prime}$, respectively. The function
$a^{2}(\alpha)=2\alpha(1-\alpha)\,m_{Q}\,\omega$, where the oscillatory
frequency $\omega=(\left.M_{V^{\prime}}\right.-\left.M_{V}\right.)/2\approx
0.3\,\,\mbox{GeV}$. This allowed to obtain an explicit analytical harmonic
oscillatory form for the Green function fg , what substantially simplified the
calculations of the $l_{c}$-correction factors $F^{coh}$ and $F^{inc}$.

The Fig. 2 clearly shows that more accurate calculations of factors $F^{coh}$
and $F^{inc}$ within a rigorous Green function formalism (solid lines) are
different from the simplified standard estimations based on the VDM (dashed
line). For the $1S$-quarkonium states the variation between both results is
not very large. Not so for radially excited quarkonia where the latter
approximation cannot be applied anymore for calculations of the
$l_{c}$-correction factors. The more detailed comparison of both approaches in
heavy quarkonium production will be presented elsewhere prepare .

One can see from Fig. 2 that effects of a finite CL are effective for energy
regions $s\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}10^{3}\,\,\mbox{GeV}^{2}$ and
$s\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$<$}}10^{4}\,\,\mbox{GeV}^{2}$ in production of charmonia and
bottomonia, respectively. This is a direct consequence of the CL dependence on
a quarkonium mass as given by Eq. (2). Figure 2 also demonstrates that a
contraction of the CL toward smaller values of $s$ leads to a significant
reduction of coherent cross sections for the $1S$ quarkonium states. However,
the corresponding incoherent cross sections are enhanced only slightly, by
$\sim 20\%-30\,\%$ and $\sim 7\%-12\,\%$ for production of $J\\!/\\!\psi(1S)$
and $\Upsilon(1S)$, respectively.

For radially excited $2S$ quarkonia the onset of the finite-$l_{c}$ effects
generated by the Green function formalism is significantly different from VDM.
This is most visible in the incoherent case as a direct consequence of the
nodal structure of quarkonium wave functions. Such a node effect is stronger
for the $\psi^{\,\prime}(2S)$ compared to $\Upsilon^{\,\prime}(2S)$ state
leading thus to much complicated non-monotonic behaviour of the factor
$F^{inc}(s)$ towards small values of $s$. On the other hand, a stronger energy
dependence from the region of small $Q\bar{Q}$ transverse separations below
the node position compared to large $Q\bar{Q}$ dipole sizes above the node
position causes a weakening of the node effect with the energy resulting thus
in a gradual convergence of factors $F^{inc}_{V^{\prime}}(s)$ and
$F^{inc}_{V}(s)$ towards large $s$.

Note that authors in the recent paper Henkels:2020kju  used the Glauber
approximation (15)-(18) for calculations of $F^{coh}$ and $F^{inc}$ repeating
only our previous studies Hufner:1996jw ; Ivanov:2002kc ; Ivanov:2002eq ;
Ivanov:2007ms . As we have already mentioned above, such an approximation
cannot be applied properly for the photo- and electroproduction of heavy
quarkonia off nuclei. Moreover, considering the same paper, the photon energy
given by Eq. (3.1) in Henkels:2020kju  cannot provide the $l_{c}$-values in
the target rest frame according to subsequent Eq. (3.4) in Henkels:2020kju .
This fact generates the question about a precise inclusion of the
finite-$l_{c}$ effects in calculations and, consequently, about a correctness
of predictions for nuclear production cross sections in regions of large
positive resp. negative rapidities. Besides, the corresponding analysis itself
of the onset of coherence effects as function of rapidity and collision energy
is missing. Especially, their strong onset is expected in the RHIC kinematic
region, which is also excluded in Ref. Henkels:2020kju  from a theoretical
study of quarkonium production in UPC. Here, apart from the question about a
validity of the Glauber approximation for estimation of factors $F^{coh}$ and
$F^{inc}$, the corresponding analysis of finite-$l_{c}$ effects is based on
the lack of calculational details for the total quarkonium-nucleon cross
section, which is an important ingredient of $l_{c}$-correction form factors
(see Eqs. (17) and (18)).

###  II.5 Gluon shadowing

Leading twist gluon shadowing (GS) was introduced within the dipole
representation in Kopeliovich:2001xj  and applied to photoproduction of vector
mesons on nuclei in Refs. Ivanov:2002kc ; Nemchik:2002ug ; Kopeliovich:2007wx
. Such effects are generated by higher Fock components of the photon
containing besides the $Q\bar{Q}$ pair also gluons, i.e. $|Q\bar{Q}G\rangle$,
$|Q\bar{Q}2G\rangle$ … $|Q\bar{Q}nG\rangle$. Since at small dipole sizes
$\vec{r}$ the dipole cross section $\sigma_{Q\bar{Q}}(r,s)$ depends on the
gluon distribution in the target, nuclear shadowing of the gluon distribution
is expected to reduce $\sigma_{Q\bar{Q}}(r,s)$ in nuclear reactions with
respect to processes on the nucleon,

$\sigma_{Q\bar{Q}}(r,s)\Rightarrow\sigma_{Q\bar{Q}}(r,s)\cdot
R_{G}(s,\mu^{2},b)\,.$ (30)

Here the gluon shadowing factor $R_{G}$, defined as the nucleus-to-nucleon
ratio of gluon densities,
$R_{G}(s,\mu^{2},b)=G_{A}(s,\mu^{2},b)/(A\,G_{N}(s,\mu^{2}))$, is probed at
the factorisation scale $\mu^{2}\sim\tau M_{V}^{2}=B/r_{S}^{2}$, where
$\tau=B/Y^{2}$. The large scale factor $B\approx 9-10$ has been estimated in
Ref. Nikolaev:1994ce  and the factor $Y$ is related to the scanning radius
$r_{S}$ Kopeliovich:1993gk ; Kopeliovich:1993pw ; Nemchik:1994fp ;
Nemchik:1996cw ; Kopeliovich:2001xj ; Kopeliovich:2007wx ; Krelina:2018hmt  as
$Y=r_{s}\,M_{V}$. The factors $Y$ have been obtained in Ref. Krelina:2018hmt
with corresponding values as presented in Tab. 2 for various quarkonium states
using wave functions generated by the BT and POW potentials.

|  $J\\!/\\!\psi(1S)$ |  $\psi^{\,\prime}(2S)$ |  $\Upsilon(1S)$ |  $\Upsilon^{\,\prime}(2S)$
---|---|---|---|---
$~{}~{}~{}~{}Y:BT$ | $5.3$ | $4.7$ | $5.7$ | $5.3$
$~{}~{}~{}~{}Y:POW$ | $5.8$ | $5.0$ | $6.0$ | $5.7$
Table 2:  Values of the factor $Y$ Krelina:2018hmt  for various quarkonium
states with wave functions generated by the BT and POW potentials.

Within the LF dipole formalism based on the path integral technique we
calculate the gluon shadowing correction corresponding to the lowest Fock
component $|q\bar{q}G\rangle$ containing only one gluon Kopeliovich:1999am ;
Kopeliovich:2001xj ; Ivanov:2002kc ; Nemchik:2002ug ; Kopeliovich:2008ek ;
Krelina:2020ipn . Consequently, eikonalization of the $R_{G}$ factor
Kopeliovich:2001ee  by Eq. (30) includes effectively the effect of higher
multi-gluon Fock states of the photon simulating thus the interference effects
in a nucleus (for an alternative estimation of gluon shadowing effects in
charmonium production in UPC see guzey-1 ; guzey-2 , for example).

Figure 3:  Gluon shadowing correction for photoproduction of $J\\!/\\!\psi$
(left panel) and $\Upsilon$ (right panel) off lead target as function of the
square of c.m. energy $s$ for the photon-lead system taking several fixed
values of the nuclear path length $L=1,3,5,10$ and $14\,\,\mbox{fm}$.

As an example, using the KST dipole model for $\sigma_{Q\bar{Q}}$, we plot in
Fig. 3 the magnitude of the gluon shadowing factor $R_{G}$ for $J\\!/\\!\psi$
and $\Upsilon$ photoproduction off lead target. Here we show that the GS rises
with c.m energy $s$ at several fixed values of the nuclear path length
$L=T_{Pb}(b)/\rho_{0}$, where $\rho_{0}$ is determined from the normalization
condition $\int d^{2}b\,dz\,\rho_{Pb}(b,z)=A$. One can see that such a
correction becomes effective only at relatively high energies starting from
$s\approx 10^{3}\,\,\mbox{GeV}^{2}$ and $s\approx 10^{4}\,\,\mbox{GeV}^{2}$
for production of charmonia and bottomonia, respectively. In comparison to
$J\\!/\\!\psi$ photoproduction, the onset of GS effects is much weaker for
$\Upsilon$ due to much larger factorization scale related to the bottomonium
mass.

##  III Model predictions vs available data

In order to calculate the nuclear cross sections (9) for the incoherent
(quasi-elastic) process one should know the slope parameter for the elastic
process $\gamma N\to VN$. Here we rely on the standard Regge form,

$B_{J\\!/\\!\psi}(s)=B_{0}+2\,\alpha^{\prime}(0)\ln\Big{(}\frac{s}{s_{0}}\Big{)}\,,$
(31)

where the parameters $\alpha^{\prime}=0.171\,\mbox{GeV}^{-2}$, the slope of
the Pomeron trajectory, and $B_{0}=1.54\,\,\mbox{GeV}^{-2}$ were fitted in
Cepila:2019skb  to data on $J\\!/\\!\psi$ photoproduction with
$s_{0}=1\,\mbox{GeV}^{2}$. The slope for $1S$ bottomonium photoproduction was
fitted to data in Cepila:2019skb  and found to have a smaller value than for
$J\\!/\\!\psi$, $B_{\Upsilon}(s)\approx
B_{J\\!/\\!\psi}(s)-1\,\mbox{GeV}^{-2}$. For production of radially excited
$2S$ state of bottomonia, the node effect is negligibly small and one can
safely use the same magnitudes of the slope parameter for both $1S$ and $2S$
states, i.e. $B_{\Upsilon^{\,\prime}}(s)\sim B_{\Upsilon}(s)$. Not so for
production of radially excited charmonia where the difference of diffraction
slopes $\Delta_{B}(s)=B_{J\\!/\\!\psi}(s)-B_{\psi^{\,\prime}}(s)$ cannot be
neglected. Here we adopt a parametrization of the factor $\Delta_{B}(s)$ from
Ref. Cepila:2019skb  (see also Ref. jan-98 ).

Figure 4:  Rapidity distributions of coherent (left panels) and incoherent
(right panels) charmonium photoproduction in UPC at RHIC collision energy
$\sqrt{s_{N}}=200\,\,\mbox{GeV}$ (top panels) and at LHC energies
$\sqrt{s_{N}}=2.76\,\,\mbox{TeV}$ (middle panels) and
$\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$ (bottom panels). The nuclear cross sections
are calculated with charmonium wave functions generated by the POW (dashed
lines) and BT (solid lines) potentials and with KST (thick lines) and GBW
(thin lines) models for the dipole cross section. The data are taken from CMS
Khachatryan:2016qhq , ALICE Abelev:2012ba ; Abbas:2013oua ; Adam:2015sia ;
Acharya:2019vlb  and LHCb LHCb:2018ofh  collaborations.

We have calculated the coherent and incoherent heavy quarkonium
photoproduction according to Eq. (3). The corresponding nuclear cross sections
are calculated in the limit of large photon energy, when the CL (2) exceeds
substantially the nuclear radius $R_{A}$, using Eqs.(5) and (9).

For the dipole cross section inherent in nuclear cross sections we adopt KST
and GBW phenomenological parametrizations of a saturated form. Here,
analogously as in Ref. Kutak:2003bd , at small c.m. energies squared $s$
corresponding to large values of $x=M_{V}^{2}/s$ we included a modification of
the dipole model through an additional factor $(1-x)^{7}$, which follows from
the see quark counting rules.

The quarkonium wave functions in the $Q\bar{Q}$ rest frame are obtained as a
solution of the Schrödinger equation for two distinct realistic $Q-\bar{Q}$
interaction potentials, POW and BT. This requires to perform subsequently the
Lorentz boost of their radial parts to the LF frame. Here we used the standard
prescription from Ref. Terentev:1976jk . The spin-dependent part is a subject
to the Melosh spin transformation as presented in Refs. Hufner:2000jb ;
Krelina:2018hmt ; Cepila:2019skb , which contain also details about the photon
and quarkonium wave functions. Such spin transformation has been incorporated
in calculations of nuclear cross section via Eqs. (11)-(13).

We include in our predictions two main phenomena affecting the nuclear cross
sections: the gluon shadowing and the finite-$l_{c}$ corrections. Whereas the
former is described in Sect. II.5 and dominates at large photon energies, the
latter prefers smaller energies (when the CL of the short-lived $Q\bar{Q}$
fluctuations is comparable with the nuclear radius) and is included via
$l_{c}$-correction factors calculated within a rigorous Green function
formalism as described in Sect. II.4. Our predictions for differential cross
sections $d\sigma/dy$ of quarkonium production in UPC as function of the
rapidity $y$ have been performed for one RHIC $\sqrt{s_{N}}=200\,\,\mbox{GeV}$
and two different LHC c.m. collision energies $\sqrt{s_{N}}=2.76$ and
$5.02\,\,\mbox{TeV}$.

In Fig. 4 we present our results for the rapidity distributions of coherent
(left panels) and incoherent (right panels) charmonium photoproduction in UPC
obtained for $\sqrt{s_{N}}=200\,\,\mbox{GeV}$ (top panels),
$\sqrt{s_{N}}=2.76\,\,\mbox{TeV}$ (middle panels) and
$\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$ (bottom panels). Calculations have been
performed for charmonium wave functions generated by two distinct $Q\bar{Q}$
potentials, POW (dashed lines) and BT (solid lines). For the dipole cross
sections $\sigma_{Q\bar{Q}}$ we adopted two different parameterizations, KST
(thick lines) and GBW (thin lines). Here the model predictions are tested by
the LHC data from the CMS Khachatryan:2016qhq  and ALICE Abelev:2012ba ;
Abbas:2013oua ; Adam:2015sia  collaborations at c.m. collision energy
$\sqrt{s_{N}}=2.76\,\,\mbox{TeV}$, as well as by the LHCb LHCb:2018ofh  and
ALICE Acharya:2019vlb  data at $\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$.

One can see from Fig. 4 that values of $d\sigma/dy$ are strongly correlated
with the shape of quarkonium wave functions determined from various models for
$Q-\bar{Q}$ interaction potentials. Whereas in charmonium production in UPC
the POW (dashed lines) and BT (solid lines) models lead to rather different
predictions for $d\sigma/dy$, in the bottomonium case the both models give
very similar results (see Fig. 6), what is in correspondence with our previous
studies Cepila:2019skb  of quarkonium electroproduction off protons. However,
it is in contrast with results from Ref. Henkels:2020kju , where a rather weak
or strong sensitivity of $d\sigma/dy$ to the shape of charmonium or
bottomonium wave functions has been predicted, respectively. Such results have
been achieved only at the cost of universal heavy quark masses
$m_{c}=1.4\,\,\mbox{GeV}$ and $m_{b}=4.75\,\,\mbox{GeV}$ used in a given
analysis for all $Q-\bar{Q}$ interaction potentials without any justification.
Because of a strong sensitivity of quarkonium cross sections
$\sigma(\gamma^{*}N\to VN)$ to values of the quark mass as analysed in Ref.
Cepila:2019skb ), the above mentioned artificial and universal values for
$m_{c}$ and $m_{b}$, that are different from those extracted from various
distinct potential models, cannot provides us with proper and consistent
predictions for nuclear phenomena occurring in quarkonium production in UPC.

Figure 5:  The same as Fig. 4 but for the $\psi^{\,\prime}(2S)$ production in
UPC. The experimental value is taken from the ALICE Adam:2015sia
collaboration.

Figure 6:  The same as Fig. 4 but for the bottomonium production in UPC.

The experimental data on production of radially excited heavy quarkonia in UPC
are very scanty. The ALICE collaboration Adam:2015sia  has been measured
$d\sigma/dy$ for coherent production of $\psi^{\,\prime}(2S)$ at $y=0$ as is
depicted in the left panel of Fig. 5 together with our predictions. One can
see a reasonable agreement for both $Q-\bar{Q}$ interaction potentials, as
well as for both models of the dipole cross section. For completeness, the
right panel of Fig. 5 shows predicted values of $d\sigma/dy$ for the
incoherent production of $\psi^{\,\prime}(2S)$ in UPC that can be verified by
future measurements.

Figure 7:  Manifestation of particular nuclear effects in coherent charmonium
(left panels) and bottomonium (right panels) photoproduction in UPC at RHIC
collision energy $\sqrt{s_{N}}=200\,\,\mbox{GeV}$ (top panels) and at the LHC
energy $\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$ (bottom panels). The nuclear cross
sections are calculated with charmonium wave functions generated by the BT
potential adopting the KST model for the dipole cross section. The dotted
lines represent predictions in the high energy limit, Eqs. (5) and (9),
containing only the maximal shadowing from the long-lived lowest $Q\bar{Q}$
Fock states of the photon. The dashed and solid lines additionally include the
spin rotation effects and corrections to a finite CL with the gluon shadowing,
respectively. The data are taken from the PHENIX Afanasiev:2009hy , CMS
Khachatryan:2016qhq , ALICE Abelev:2012ba ; Abbas:2013oua ; Adam:2015sia ;
Acharya:2019vlb  and LHCb LHCb:2018ofh  collaborations.

The next Fig. 6 represents analogous predictions as Fig. 4 but for production
of bottomonia in UPC. Here the KST and GBW parametrizations give rather
different results in a correspondence with our analysis of the process $\gamma
N\to\Upsilon N$ in Ref. Cepila:2019skb . Differences in predictions using the
KST and GBW models for $\sigma_{Q\bar{Q}}(r)$ can be treated as a measure of
the underlined theoretical uncertainty.

Figure 7 demonstrates how the onset of particular nuclear effects is
manifested in distributions $d\sigma/dy$ at RHIC (top panels) and the LHC
(bottom panels) energies in coherent production of charmonia (left panels) and
bottomonia (right panels) in UPC. Here the dotted lines represent our
calculations in the standard high energy limit, including only the maximal
shadowing generated by the lowest long-lived $Q\bar{Q}$ Fock components of the
photon, without any other corrections and effects occurring in interactions
with nuclei. The dashed and solid lines incorporate additionally and
subsequently the spin rotation effects and corrections to a finite CL with the
gluon shadowing, respectively. Here our predictions correspond to the KST
model for the dipole cross section.

Figure 8:  The same as Fig. 7 but for coherent production of radially excited
$2S$ quarkonia in UPC.

Since in the RHIC kinematic region the gluon shadowing effects can be
neglected, the differences between solid and dashed lines in top panels of
Fig. 7 represent the onset of the finite-$l_{c}$ corrections which is more
pronounced in coherent production of bottomonia compared to $J\\!/\\!\psi$.
Whereas at the LHC collision energy $\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$ such
differences exhibit a dominance of former effects mainly in the midrapidity
($y=0$) region, the latter effects manifest themselves in regions of large
negative and positive rapidities. The both effects substantially reduce the
nuclear cross sections $d\sigma/dy$. In comparison to charmonium production at
the LHC, a much weaker onset of gluon shadowing effects in production of
bottomonia is caused by a larger corresponding scale $\propto
M_{\Upsilon}^{2}\gg M_{J\\!/\\!\psi}^{2}$ (see also Fig. 3).

Analogously, the Fig. 8 shows contributions of above nuclear effects to a
magnitude of rapidity distributions $d\sigma/dy$ for coherent production of
radially excited $2S$ quarkonia. Here in contrast to our predictions for
coherent $J\\!/\\!\psi$ production in Fig. 7, where the spin rotation effects
considerably enhance nuclear cross sections, the nodal structure of
$\psi^{\,\prime}(2S)$ wave function leads to a counter-intuitive reduction of
$d\sigma/dy$. This can be seen in the left panels of Fig. 8 as a difference
between dashed and dotted lines. In comparison to the RHIC energy range, such
a reduction is smaller at the larger LHC energy due to a weakening of the node
effect with energy. In coherent production of $\Upsilon^{\,\prime}(2S)$ at
$\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$ a very weak node effect leads to a standard
enhancement of $d\sigma/dy$ as is expected for $1S$ quarkonia (see the bottom-
right panel of Fig. 8).

Figure 9:  The same as Fig. 7 but for incoherent production of the ground
state $1S$ (top panels) and radially excited $2S$ (bottom panels) quarkonia in
UPC at $\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$.

The last Fig. 9 illustrates a manifestation of particular nuclear effects at
the LHC collision energy $\sqrt{s_{N}}=5.02\,\,\mbox{TeV}$ in incoherent
production of the ground state $1S$ (top panels) and radially excited $2S$
(bottom panels) quarkonia. Since we expect very small corrections for a finite
CL and weak gluon shadowing effects at small energies (see Figs. 2 and 3),
here we do not present our predictions for the RHIC kinematic region.
Differences between the solid and dashed lines represent the onset of net
gluon shadowing effects at the LHC which is maximal at $y=0$.

Finally, for completeness, we present in Tab. 3 the values of cross sections
$\sigma_{coh}(V)$ integrated over rapidity $y$ for coherent $AA\to VAA$
production of various quarkonium states
[$V=J\\!/\\!\psi,\psi^{\,\prime},\Upsilon,\Upsilon^{\,\prime}$] in UPC using
wave functions generated by the BT and POW potentials and adopting KST and GBW
models for the dipole cross section. The corresponding values of
$\sigma_{coh}$ clearly demonstrates that the onset of the node effect in
coherent production of $\psi^{\,\prime}(2S)$ and $\Upsilon^{\,\prime}(2S)$ is
weaker at larger collision energies, what is manifested as a rise of ratios
$\sigma_{coh}(\psi^{\,\prime})/\sigma_{coh}(J\\!/\\!\psi)$ and
$\sigma_{coh}(\Upsilon^{\,\prime})/\sigma_{coh}(\Upsilon)$ with
$\sqrt{s_{N}}$.

|  AuAu ($\sqrt{s_{N}}=200\,\mbox{GeV}$) |  PbPb ($\sqrt{s_{N}}=2.76\,\mbox{GeV}$) |  PbPb ($\sqrt{s_{N}}=5.5\,\mbox{TeV}$)
---|---|---|---
$~{}J\\!/\\!\psi$: GBW, POW | $237.0~{}~{}\,\mu\mbox{b}$ | $20.9~{}~{}\,\mbox{mb}$ | $36.0~{}~{}\,\mbox{mb}$
$~{}J\\!/\\!\psi$: GBW, BT | $162.8~{}~{}\,\mu\mbox{b}$ | $15.4~{}~{}\,\mbox{mb}$ | $26.3~{}~{}\,\mbox{mb}$
$~{}\psi^{\,\prime}$ : GBW, BT | $10.12~{}~{}\,\mu\mbox{b}$ | $2.61~{}~{}\,\mbox{mb}$ | $4.91~{}~{}\,\mbox{mb}$
$~{}J\\!/\\!\psi$: KST, POW | $235.1~{}~{}\,\mu\mbox{b}$ | $19.7~{}~{}\,\mbox{mb}$ | $33.8~{}~{}\,\mbox{mb}$
$~{}J\\!/\\!\psi$: KST, BT | $159.9~{}~{}\,\mu\mbox{b}$ | $13.7~{}~{}\,\mbox{mb}$ | $24.4~{}~{}\,\mbox{mb}$
$~{}\psi^{\,\prime}$ : KST, BT | $15.9~{}~{}\,\mu\mbox{b}$ | $2.62~{}~{}\,\mbox{mb}$ | $5.01~{}~{}\,\mbox{mb}$
$~{}\Upsilon~{}~{}$: GBW, POW |  $0.97~{}~{}\,\mbox{nb}$ | $8.18~{}~{}\,\mu\mbox{b}$ | $18.3~{}~{}\,\mu\mbox{b}$
$~{}\Upsilon~{}~{}$: GBW, BT |  $0.99~{}~{}\,\mbox{nb}$ | $8.39~{}~{}\,\mu\mbox{b}$ | $19.3~{}~{}\,\mu\mbox{b}$
$~{}\Upsilon^{\,\prime}$ : GBW, BT |  $0.17~{}~{}\,\mbox{nb}$ | $2.24~{}~{}\,\mu\mbox{b}$ | $5.44~{}~{}\,\mu\mbox{b}$
$~{}\Upsilon~{}~{}$: KST, POW |  $2.83~{}~{}\,\mbox{nb}$ | $21.4~{}~{}\,\mu\mbox{b}$ | $50.6~{}~{}\,\mu\mbox{b}$
$~{}\Upsilon~{}~{}$: KST, BT |  $2.88~{}~{}\,\mbox{nb}$ | $22.1~{}~{}\,\mu\mbox{b}$ | $52.5~{}~{}\,\mu\mbox{b}$
$~{}\Upsilon^{\,\prime}$ : KST, BT |  $0.53~{}~{}\,\mbox{nb}$ | $6.74~{}~{}\,\mu\mbox{b}$ | $16.1~{}~{}\,\mu\mbox{b}$
Table 3:  Total (integrated over rapidity $y$) cross sections
$\sigma_{coh}(V)$ for coherent $AA\to V\\!AA$ vector meson
($V=J\\!/\\!\psi,\psi^{\,\prime},\Upsilon,\Upsilon^{\,\prime}$) production in
UPC using quarkonium wave functions obtained from the POW and BT potentials
and adopting GBW and KST models for the dipole cross section.

##  IV Conclusions

In this paper we treat the heavy quarkonium
[$J\\!/\\!\psi(1S),\psi^{\,\prime}(2S),\Upsilon(1S),\Upsilon^{\,\prime}(2S)$,…]
production in heavy-ion UPC in the energy range accessible by experiments at
RHIC and the LHC. Here the main observations are the following:

  * •

Calculations of the corresponding photoproduction cross sections for the
coherent, as well as incoherent process have been performed within the LF
color dipole approach. Here we used the standard formulas for nuclear cross
sections in the high energy limit where the CL is much longer than the nuclear
radius, $l_{c}\gg R_{A}$. However, we included also corrections for a finite
coherence length via form factors, which have been treated for the first time
within a rigorous Green function formalism. The standard Glauber approximation
based on the VDM is rather vague for the $1S$ quarkonium states. However, it
cannot be applied anymore for production of radially excited $2S$ quarkonia
(see differences between VDM and Green function results in Fig. 2). Such the
finite-$l_{c}$ corrections for a more proper inclusion of the quark shadowing
effects are absent in the most of recent papers.

  * •

We have found that the finite-$l_{c}$ corrections are important especially in
the RHIC energy range modifying significantly rapidity distributions
$d\sigma/dy$ at small photon energies when $l_{c}\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}R_{A}$, i.e. in
regions of large negative and positive rapidities. In the LHC kinematic region
such a modification of $d\sigma/dy$ is rather small. Whereas in coherent
production of quarkonia these corrections reduce substantially nuclear cross
sections, in the incoherent process they lead to a weak enhancement of
$d\sigma/dy$ (see Fig. 2). In comparison to $1S$ quarkonium states we predict
much weaker onset of CL effects in incoherent production of radially excited
$2S$ quarkonia due to nodal structure of their wave functions. For incoherent
production of $\Upsilon^{\,\prime}(2S)$ these effects are practically
invisible and can be neglected (see right panels of Fig. 2).

  * •

At larger photon energies we included also the gluon shadowing corrections
related to higher multi-gluon Fock components of the photon. Here we did not
adopt any phenomenological parametrization of the gluon shadowing factor
$R_{G}$ from the literature but calculated this factor directly within the LF
color dipole formalism based on the Green function technique. This allows to
obtain the magnitude of GS as function of the nuclear impact parameter, as
well as to include naturally the correct energy-dependent lifetime of such
higher Fock states of the photon. Note that using high-energy expressions,
Eqs. (5) and (9), an alternative calculations of GS effects based on a very
popular Balitsky-Kovchegov equation Balitsky:1995ub ; Kovchegov:1999yj , which
sums up the all but only long-lived multi-gluon photon fluctuations, cannot
provide us with reliable results and can be applied only at asymptotically
high energies when transverse sizes of all Fock states are “frozen” during
propagation through the nucleus.

  * •

In the LHC kinematic region we have found a significant onset of GS effects
causing $\sim 10\%-20\%$ reduction of nuclear cross sections for the coherent
and incoherent charmonium photoproduction. They gradually diminish toward
smaller energies and can be neglected in kinematic regions accessible by RHIC
experiments (see Fig. 3).

  * •

The quarkonium wave functions are well defined in the $Q\bar{Q}$ rest frame.
They have been included in our calculations and obtained by solving the
Schrödinger equation for two distinct realistic $Q-\bar{Q}$ interaction
potentials, POW and BT, which give the best agreement with available data on
charmonium real and virtual photoproduction off protons. Consequently, the
corresponding LF wave functions have been generated performing the boosting to
the LF frame using the so called Terent’ev prescription from Ref.
Terentev:1976jk , whose validity has been verified in Ref. Kopeliovich:2015qna
for symmetric $c\bar{c}$ photon fluctuations.

  * •

The spin-dependent part of the wave function for $S$-wave quarkonia can be
safely factorized from the radial component. In order to avoid an undesirable
$D$-wave admixture in the $Q\bar{Q}$ rest frame, here we treat a simple non-
photon-like structure of the $V\to Q\bar{Q}$ transition. Consequently, we
perform explicitly the transformation of 2-dim. heavy (anti)quark spinors from
the rest to the LF frame known as the Melosh spin rotation. We derived new
formulas for coherent and incoherent nuclear cross sections incorporating such
a transformation. We have found that spin effects have a significant impact on
the magnitude of nuclear cross sections. Standardly, they cause a substantial
enhancement of nuclear production rates. Such a scenario of enlargement
corresponds to results on quarkonium photoproduction off protons Hufner:2000jb
; Krelina:2018hmt ; Cepila:2019skb . However, in production of radially
excited $2S$ charmonia, the Melosh spin transformation leads to a counter-
intuitive reduction of nuclear rapidity distributions $d\sigma/dy$ as a
consequence of the nodal structure of corresponding radial wave functions. In
production of the $\Upsilon^{\,\prime}(2S)$ state, the node effect becomes
much weaker and leads again to the standard scenario of
$d\sigma/dy$-enhancement.

  * •

Besides the quarkonium wave functions generated by POW and BT potentials, we
have studied also differences in our predictions employing KST and GBW
phenomenological parametrizations for the dipole cross section in order to
estimate a corresponding measure of the theoretical uncertainty in our current
analysis. We have found that whereas in charmonium production in UPC the main
source of theoretical uncertainties is related mainly to our choice of
quarkonium wave functions, in the bottomonium case the main role of
variability in predictions is played by various models for
$\sigma_{Q\bar{Q}}(r,s)$.

  * •

Our calculations are in a rather good agreement with available data on
coherent production of $J\\!/\\!\psi(1S)$ and $\psi^{\,\prime}(2S)$ in UPC at
the LHC (see Figs. 4 and 5). They can be tested not only by the future
measurements of UPC at RHIC and the LHC but also by future experiments at
planned electron-ion colliders.

Acknowledgment:

This work was supported in part by grants ANID - Chile FONDECYT 1170319, by
ANID PIA/APOYO AFB180002, by USM internal project PI_LI_19_13, and by the
project of the European Regional Development Fund No.
CZ.02.1.01/0.0/0.0/16_019/0000778. The work of J.N. was partially supported by
Grant No. LTT18002 of the Ministry of Education, Youth and Sports of the Czech
Republic and by the Slovak Funding Agency, Grant No. 2/0007/18.

## References

  * (1) B.Z. Kopeliovich and B.G. Zakharov; Phys. Rev. D 44, 3466 (1991).
  * (2) B.Z. Kopeliovich, J. Nemchik, N.N. Nikolaev and B.G. Zakharov; Phys. Lett. B 324, 469 (1994).
  * (3) J. Nemchik, N. N. Nikolaev, E. Predazzi and B. Zakharov, Phys. Lett. B 374, 199-204 (1996).
  * (4) J. Nemchik, N. N. Nikolaev, E. Predazzi and B. Zakharov, Z. Phys. C 75, 71-87 (1997).
  * (5) J. Nemchik, Phys. Rev. C 66, 045204 (2002).
  * (6) M. Krelina, J. Nemchik, R. Pasechnik and J. Cepila, Eur. Phys. J. C 79, no.2, 154 (2019).
  * (7) J. Cepila, J. Nemchik, M. Krelina and R. Pasechnik, Eur. Phys. J. C 79, no.6, 495 (2019).
  * (8) M. Krelina, J. Nemchik and R. Pasechnik, Eur. Phys. J. C 80, no.2, 92 (2020).
  * (9) I.P. Ivanov, N.N. Nikolaev and A.A. Savin; Phys. Part. Nucl. 37, 1 (2006).
  * (10) N. Brambilla et al. [Quarkonium Working Group], arXiv:hep-ph/0412158.
  * (11) N. Brambilla et al.; Eur. Phys. J. C 71, 1534 (2011).
  * (12) S. Afanasiev et al. [PHENIX], Phys. Lett. B 679, 321-329 (2009).
  * (13) V. Khachatryan et al. [CMS], Phys. Lett. B 772, 489-511 (2017).
  * (14) B. Abelev et al. [ALICE], Phys. Lett. B 718, 1273-1283 (2013).
  * (15) E. Abbas et al. [ALICE], Eur. Phys. J. C 73, no.11, 2617 (2013).
  * (16) J. Adam et al. [ALICE], Phys. Lett. B 751, 358-370 (2015).
  * (17) A. Bursche [LHCb], Nucl. Phys. A 982, 247-250 (2019).
  * (18) S. Acharya et al. [ALICE], Phys. Lett. B 798, 134926 (2019).
  * (19) B. Kopeliovich, J. Nemchik, N. N. Nikolaev and B. Zakharov, Phys. Lett. B 309, 179 (1993).
  * (20) J. Nemchik, N. N. Nikolaev and B. Zakharov, Phys. Lett. B 341, 228 (1994).
  * (21) B. Kopeliovich, J. Nemchik, A. Schafer and A. Tarasov, Phys. Rev. C 65, 035201 (2002).
  * (22) B. Kopeliovich, J. Nemchik and I. Schmidt, Phys. Rev. C 76, 025210 (2007).
  * (23) J. Nemchik, N. N. Nikolaev and B. Zakharov, Phys. Lett. B 339, 194-200 (1994).
  * (24) J. Hufner, Y.P. Ivanov, B.Z. Kopeliovich and A.V. Tarasov; Phys. Rev. D 62, 094022 (2000).
  * (25) J. Nemchik; Phys. Rev. D 63, 074007 (2001).
  * (26) J. Nemchik; Eur. Phys. J. C 18, 711 (2001).
  * (27) Y. Ivanov, B. Kopeliovich, A. Tarasov and J. Hufner, Phys. Rev. C 66, 024903 (2002).
  * (28) Y. P. Ivanov, B. Kopeliovich, A. Tarasov and J. Hufner, AIP Conf. Proc. 660, 283 (2003).
  * (29) Y. Ivanov, B. Kopeliovich and I. Schmidt, arXiv:0706.1532 [hep-ph].
  * (30) A. H. Mueller, Nucl. Phys. B 335, 115 (1990).
  * (31) B.Z. Kopeliovich, J. Raufeisen and A.V. Tarasov, Phys. Rev. C 62, 035204 (2000).
  * (32) M. Krelina and J. Nemchik, arXiv:2010.00329 [hep-ph].
  * (33) A. Martin; Phys. Lett. B 93, 338 (1980).
  * (34) N. Barik and S.N. Jena; Phys. Lett. B 97, 265 (1980).
  * (35) W. Buchmuller and S.H.H. Tye; Phys. Rev. D 24, 132 (1981).
  * (36) B.Z. Kopeliovich, A. Schafer and A.V. Tarasov; Phys. Rev. D 62, 054022 (2000).
  * (37) K.J. Golec-Biernat and M. Wusthoff; Phys. Rev. D 59, 014017 (1998).
  * (38) K.J. Golec-Biernat and M. Wusthoff; Phys. Rev. D 60, 114023 (1999).
  * (39) M.V. Terentev; Sov. J. Nucl. Phys. 24, 106 (1976) [Yad. Fiz. 24, 207 (1976)].
  * (40) H.J. Melosh; Phys. Rev. D 9, 1095 (1974).
  * (41) H.De Vries, C.W.De Jager and C.De Vries, Atom. Data Nucl. Data Tabl. 36, 495 (1987).
  * (42) T. Bauer, R. Spital, D. Yennie and F. Pipkin, Rev. Mod. Phys. 50, 261 (1978).
  * (43) B. Z. Kopeliovich, L. I. Lapidus and A. B. Zamolodchikov, JETP Lett. 33, 595 (1981).
  * (44) V. N. Gribov, Sov. Phys. JETP 29, 483 (1969) [Zh. Eksp. Teor. Fiz. 56, 892 (1969)].
  * (45) B. Z. Kopeliovich, Int. J. Mod. Phys. A 31, no.28n29, 1645021 (2016).
  * (46) J. Hufner, B. Kopeliovich and J. Nemchik, Phys. Lett. B 383, 362-366 (1996).
  * (47) B. Z. Kopeliovich, I. K. Potashnikova and I. Schmidt, Phys. Rev. C 73, 034901 (2006).
  * (48) J.B. Bronzan, G.L. Kane and U.P. Sukhatme; Phys. Lett. B 49, 272 (1974).
  * (49) J.R. Forshaw, R. Sandapen and G. Shaw; Phys. Rev. D 69, 094013 (2004).
  * (50) C. Henkels, E. G. de Oliveira, R. Pasechnik and H. Trebien, Phys. Rev. D 102, 014024 (2020).
  * (51) J. Hufner, B. Kopeliovich and A. B. Zamolodchikov, Z. Phys. A 357, 113 (1997).
  * (52) J. Hufner and B. Z. Kopeliovich, Phys. Lett. B 426, 154-160 (1998).
  * (53) R.P. Feynman and A.R. Gibbs, Quantum Mechanics and Path Integrals, McGraw-Hill Book Company, NY 1965.
  * (54) M. Krelina, B. Z. Kopeliovich and J. Nemchik, paper in preparation.
  * (55) N. N. Nikolaev and B. Zakharov, Phys. Lett. B 332, 184-190 (1994).
  * (56) B. Kopeliovich, J. Nemchik, I. Potashnikova and I. Schmidt, J. Phys. G 35, 115010 (2008).
  * (57) M. Krelina and J. Nemchik, Eur. Phys. J. Plus 135, 444 (2020).
  * (58) B. Kopeliovich, A. Tarasov and J. Hufner; Nucl. Phys. A 696, 669 (2001).
  * (59) V. Guzey and M. Zhalov, JHEP 10, 207 (2013).
  * (60) V. Guzey, E. Kryshen, M. Strikman and M. Zhalov, Phys. Lett. B 726, 290-295 (2013).
  * (61) J. Nemchik, N.N. Nikolaev, E. Predazzi, B.G. Zakharov and V.R. Zoller; J. Exp. Theor. Phys. 86, 1054 (1998).
  * (62) K. Kutak and J. Kwiecinski, Eur. Phys. J. C 29, 521 (2003).
  * (63) I. Balitsky, Nucl. Phys. B 463, 99 (1996).
  * (64) Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999).
  * (65) B.Z. Kopeliovich, E. Levin, I. Schmidt and M. Siddikov; Phys. Rev. D 92, 034023 (2015).


# Logarithmic systolic growth for hyperbolic surfaces in every genus

Mikhail G. Katz  Department of Mathematics, Bar Ilan University, Ramat Gan
5290002 Israel<EMAIL_ADDRESS>and  Stéphane Sabourau  Univ Paris Est
Creteil, CNRS, LAMA, F-94010 Creteil, France Univ Gustave Eiffel, LAMA,
F-77447 Marne-la-Vallée, France<EMAIL_ADDRESS>
###### Abstract.

More than thirty years ago, Brooks and Buser–Sarnak constructed sequences of
closed hyperbolic surfaces with logarithmic systolic growth in the genus.
Recently, Liu and Petri showed that such logarithmic systolic lower bound
holds for every genus (not merely for genera in some infinite sequence) using
random surfaces. In this article, we show a similar result through a more
direct approach relying on the original Brooks/Buser–Sarnak surfaces.

###### Key words and phrases:

systole, hyperbolic surface, logarithmic systolic growth

###### 2020 Mathematics Subject Classification:

Primary 53C45; Secondary

##  1\. Introduction

The systole of a closed hyperbolic surface $S$, denoted by $\operatorname{{\rm
sys}}(S)$, is the length of a shortest (noncontractible) closed geodesic of
$S$. By Mumford’s compactness theorem, the systole function achieves a maximum
among all closed hyperbolic surfaces of a given genus. By considering the area
of a disk of radius $\frac{1}{2}\operatorname{{\rm sys}}(S)$, one easily sees
(see [6, Lemma 5.2.1]) that the systole of every closed hyperbolic surface $S$
of genus $g$ satisfies the upper bound

$\operatorname{{\rm sys}}(S)\leq 2\log(4g-2).$ (1.1)

As far as lower bounds are concerned, Brooks [4] and Buser–Sarnak [1]
constructed a sequence of closed hyperbolic surfaces $S_{g_{p}}$ of genus
$g_{p}$ with

$g_{p}=(p^{3}-p)\nu+1$ (1.2)

obtained as congruence coverings of a fixed arithmetic surface $S$ defined
from a division quaternion algebra $A$ such that

$\operatorname{{\rm sys}}(S_{g_{p}})\geq 4\log(p)$ (1.3)

where $p$ is any odd prime and $\nu$ is fixed and depends only on $A$; see [4,
Lemma 3.1] and [1, Eq. (4.3)]. This immediately leads to the following
logarithmic systolic growth

$\operatorname{{\rm sys}}(S_{g_{p}})\geq\frac{4}{3}\log(g_{p})-C$ (1.4)

where $C$ is a positive constant depending only on $A$, but not on $p$.

Combining the upper and lower bounds (1.1) and (1.4), we obtain

$\frac{4}{3}\leq\limsup_{g\to\infty}\max_{S\in\mathcal{M}_{g}}\frac{\operatorname{{\rm
sys}}(S)}{\log g}\leq 2$ (1.5)

where $\mathcal{M}_{g}$ is the moduli space of closed hyperbolic surfaces of
genus $g$. See also Brooks [5], Gromov [9] and Schmutz Schaller [19].

The Brooks/Buser–Sarnak construction was generalized in [12] to principal
congruence covers of any closed arithmetic surface. Makisumi [14] proved that
the multiplicative constant $\frac{4}{3}$ in (1.4) is optimal for these
congruence covers. Other sequences of closed hyperbolic surfaces with
arbitrarily large genus satisfying a logarithmic growth have been obtained in
[17] and [18] from graphs of large (logarithmic) girth/systole. Still, these
constructions provide families of surfaces of fairly sparse genera.

Recently, Liu and Petri [13] showed using a probabilistic method that for
every $c<\frac{2}{9}$, there exists a constant $C>0$ such that for _every_
genus $g\geq 2$ (not merely for genera in some infinite sequence), there
exists a closed hyperbolic surface $S_{g}$ of genus $g$ with

$\operatorname{{\rm sys}}(S_{g})\geq c\log(g)-C.$

It follows that

$\liminf_{g\to\infty}\max_{S\in\mathcal{M}_{g}}\frac{\operatorname{{\rm
sys}}(S)}{\log g}\geq\frac{2}{9}.$

Note that their argument relies on an appropriate model of random surfaces and
does not provide an explicit construction.

In the present article, we show a similar result, see Theorem 3.1, without
using probabilistic arguments (albeit with a worse constant). Instead, we use
a more direct approach relying on the original surfaces of
Brooks/Buser–Sarnak.

######  Theorem 1.1.

For every $c<\frac{19}{120}$, there exists a constant $C>0$ such that for
every genus $g\geq 2$, there exists a closed hyperbolic surface $S_{g}$ of
genus $g$ with

$\operatorname{{\rm sys}}(S_{g})\geq c\log(g)-C.$

Actually, Buser and Sarnak [1] considered a similar question for the systole
of the Jacobian of a Riemann surface of large genus after showing that there
exists a sequence of Riemann surfaces of arbitrarily large genus whose square
of the Jacobian systole has a logarithmic growth in the genus. To show this
holds in every genus, they mention that one could seek to vary the parameters
in the definition of the quaternion algebra $A$ and to vary the congruence
groups in the arithmetic construction, but as they point out this seems rather
complicated; see [1, p. 47]. Instead, they rely on the circle method from
additive analytic number theory and Fay’s gluing of Rieman surfaces [7, §III]
to construct Rieman surfaces whose Jacobian decomposes as the product of
Jacobians of lower dimension up to a small error term (see [1, p. 47]), which
is enough for their purposes. As mentioned in [16, p. 127], it seems likely
(though this would require some argument) that the homological systole of the
underlying hyperbolic surfaces thus-obtained have a logarithmic growth. It is
however unclear whether this construction yields an explicit homological
systolic lower bound (as we have not been able to fill in all the details) and
no such explicit bound for the Jacobian systole is given in [1] (only for the
surfaces constructed as congruence covers of an arithmetic surface). Note also
that, though the homological systole of these surfaces is likely to have a
logarithmic growth, their systole is equally likely to go to zero. However,
since for hyperbolic surfaces of a fixed genus, the maximum of the systole is
equal to the maximum of the homological systole by Parlier [15], this would
not be an obstacle to deriving a logarithmic systolic growth in every genus.

To prove Theorem 1.1, we rely on a more flexible construction using
$\textrm{CAT}(-1)$ metrics rather than the more rigid hyperbolic metrics. This
change of setting is not restrictive for our purpose since the maximum of the
systole on $\textrm{CAT}(-1)$ surfaces of fixed genus is achieved by a
hyperbolic metric according to [10]. Along the way, we apply the prime gap
theorem to control the growth of the genus in Brooks/Buser–Sarnak’s congruence
surfaces. As a result, we obtain an explicit asymptotic systolic lower bound.

Recent advances in systolic geometry include [8] and [11].

_Acknowledgement._ The authors would like to thank Hugo Parlier for bringing
the article [13] to their attention after receiving a first version of the
present article.

##  2\. Preliminaries

For future reference, let us recall some classical formulas in hyperbolic
geometry. The area and the perimeter of a disk $D(r)$ of radius $r$ in the
hyperbolic plane are given by

$\displaystyle\operatorname{{\rm area}}D(r)$
$\displaystyle=4\pi\sinh^{2}\left(\frac{r}{2}\right)$ (2.1)
$\displaystyle\operatorname{{\rm length}}\partial D(r)$
$\displaystyle=2\pi\sinh(r)$ (2.2)

and the area of a closed hyperbolic surface $S$ of genus $g$ is given by

$\operatorname{{\rm area}}(S)=4\pi(g-1).$ (2.3)

Due to the arithmetic nature of the construction of the Brooks/Buser–Sarnak
surfaces, we will also need the following prime gap theorem.

######  Theorem 2.1 (Prime gap theorem [3]).

There exist $\lambda\geq 1$ and $\theta\in(0,1)$ such that the gap between any
prime $p$ and the next prime $p^{\prime}$ satisfies

$p^{\prime}-p\leq\lambda p^{\theta}.$ (2.4)

For this bound, we can take $\theta=\frac{21}{40}=0.525$.

######  Remark 2.2.

For $p$ large enough, the multiplicative constant $\lambda\geq 1$ can be taken
arbitrarily close to $1$; see [3].

##  3\. Logarithmic systolic growth

We will use the notation $A_{g}\simeq B_{g}$ for $A_{g}=B_{g}(1+o(1))$ as $g$
goes to infinity (similarly with $A_{g}\lesssim B_{g}$ and $A_{g}\gtrsim
B_{g}$).

######  Theorem 3.1.

For every $g\geq 2$, there exists a closed hyperbolic surface $S_{g}$ of genus
$g$ with

$\operatorname{{\rm sys}}(S_{g})\gtrsim\frac{1-\theta}{3}\log(g)$ (3.1)

where $\theta\in(0,1)$ is the exponent in (2.4).

######  Remark 3.2.

With the value of $\theta$ given in Theorem 2.1, the coefficient
$\frac{1-\theta}{3}$ can be taken equal to
$\frac{19}{120}=0.18333...\geq\frac{1}{7}$.

###### Proof.

For every odd prime $p$, consider a closed hyperbolic surface $S_{g_{p}}$ of
genus $g_{p}$ satisfying (1.3). Let $r$ and $d$ be two positive reals (to be
determined later) such that $r+d<\frac{1}{4}\operatorname{{\rm
sys}}(S_{g_{p}})$. Consider a maximal collection $(D_{i})_{1\leq i\leq N}$ of
disjoint disks $D_{i}$ of radius $(r+d)$ in $S_{g_{p}}$. Denote by $x_{i}$ the
center of $D_{i}$. Since the collection of disks $D_{i}$ is maximal, the disks
$2D_{i}$ of radius $2(r+d)$ centered at $x_{i}$ cover $S_{g_{p}}$. Since their
radius is less than $\frac{1}{2}\operatorname{{\rm sys}}(S_{g_{p}})$, the
disks $2D_{i}$ can be isometrically embedded in the hyperbolic plane onto a
disk $2D$ of radius $2(r+d)$. It follows that the number $N$ of points $x_{i}$
satisfies

$N\cdot\operatorname{{\rm area}}(2D)\geq\operatorname{{\rm area}}(S_{g_{p}}).$

By the relations (2.1) and (2.3), we immediately obtain

$N\geq\frac{g_{p}-1}{\sinh^{2}(r+d)}.$ (3.2)

Now, remove $2k$ disks $D(x_{i},r)$ from $S_{g_{p}}$ with $k\leq\frac{N}{2}$,
and pairwise glue together the boundary components of the resulting surface.
This gives rise to a closed $\textrm{CAT}(-1)$ surface $\bar{S}_{g_{p}+k}$ of
genus $g_{p}+k$.

Before showing that the surface $\bar{S}_{g_{p}+k}$ satisfies the logarithmic
systolic growth of (3.1), we will first show that every genus large enough can
be attained as the genus of a surface $\bar{S}_{g_{p}+k}$. For that, for $p$
large enough, we want $g_{p}+\lfloor\frac{N}{2}\rfloor$ to be at least
$g_{p^{\prime}}$ to recover every genus between $g_{p}$ and $g_{p^{\prime}}$,
when $k$ runs between $0$ and $\frac{N}{2}$ (recall that $p^{\prime}$ is the
prime right after $p$). By (3.2), it is enough to have

$\frac{g_{p}-1}{2\sinh^{2}(r+d)}\geq g_{p^{\prime}}-g_{p}+1.$ (3.3)

Since $g_{p}=(p^{3}-p)\nu+1$, the genus gap $g_{p^{\prime}}-g_{p}$ satisfies

$\displaystyle g_{p^{\prime}}-g_{p}+1$ $\displaystyle=(p^{\prime
3}-p^{3})\nu+(p^{\prime}-p)\nu+1$ $\displaystyle\leq 3\nu\lambda
p^{2+\theta}+3\nu\lambda^{2}p^{1+2\theta}+\nu\lambda^{3}p^{3\theta}+\nu\lambda
p^{\theta}+1$ $\displaystyle\leq 8\nu\lambda^{3}p^{2+\theta}+1$ (3.4)

after making use of the bound $p^{\prime}\leq p+\lambda p^{\theta}$ (see
Theorem 2.1), expanding the power, and bounding each term $p^{1+2\theta}$,
$p^{3\theta}$ and $p^{\theta}$ by $p^{2+\theta}$. Thus, by (3.3), (1.2) and
(3.4), it is enough to have

$\sinh^{2}(r+d)\leq
p^{1-\theta}\,\frac{(1-p^{-2})\nu}{2(8\nu\lambda^{3}+p^{-2-\theta})}.$

Since $\sinh(r+d)\leq\frac{e^{r+d}}{2}$ and $p\geq 3$, it is enough to have

$r+d\leq\frac{1-\theta}{2}\log(p)-C(\lambda)$ (3.5)

where
$C(\lambda)=\frac{1}{2}\log\left(\frac{2(8\nu\lambda^{3}+3^{-2-\theta})}{(1-3^{-2})\nu}\right)$.
Note that this condition immediately implies that

$r+d\leq\log(p)\leq\frac{1}{4}\operatorname{{\rm sys}}(S_{g_{p}})$

for $p$ large enough, as previously required.

Now, in order to show that the surface $\bar{S}_{g_{p}+k}$ satisfies the
logarithmic systolic growth of (3.1), we will need the following lemma.

######  Lemma 3.3.

$\operatorname{{\rm sys}}(\bar{S}_{g_{p}+k})\geq\min\\{\operatorname{{\rm
sys}}(S_{g_{p}}),2\pi\sinh(r),2d\\}.$

###### Proof.

By the circle length formula (2.2), the closed geodesics of
$\bar{S}_{g_{p}+k}$ corresponding to a multiple of $\partial D(x_{i},r)$ are
of length at least

$\operatorname{{\rm length}}\partial D(x_{i},r)=2\pi\sinh(r).$

The closed geodesics of $\bar{S}_{g_{p}+k}$ intersecting $\partial
D(x_{i},r)$, but non homotopic to a multiple of $\partial D(x_{i},r)$, have to
leave $D(x_{i},r+d)$, and thus, are of length at least $2d$. Finally, the
closed geodesics of $\bar{S}_{g_{p}+k}$ which do not intersect $\partial
D(x_{i},r)$ lie in the original surface $S_{g_{p}}$, and thus, are of length
at least $\operatorname{{\rm sys}}(S_{g_{p}})$. This implies the desired
systolic lower bound. ∎

At this point, we can fix the values of $r$ and $d$ as follows. Take the real
$r$ so that $2\pi\sinh(r)=4\log(p)$. That is,

$r=\operatorname{arcsinh}\left(\tfrac{2}{\pi}\log p\right)\simeq\log\log p.$

Take also the real $d$ so that the equality case is attained in the equality
(3.5). That is,

$d=\frac{1-\theta}{2}\log(p)-r-C(\lambda)\simeq\frac{1-\theta}{2}\log(p).$

By Lemma 3.3, we deduce that

$\operatorname{{\rm sys}}(\bar{S}_{g_{p}+k})\gtrsim(1-\theta)\log(p).$

For $k\leq\lambda p^{\theta}$, we obtain that $g_{p}+k\simeq\nu p^{3}$. Thus,

$\operatorname{{\rm
sys}}(\bar{S}_{g_{p}+k})\gtrsim\frac{1-\theta}{3}\log(g_{p}+k).$

Since every genus $g$ large enough can be decomposed as $g=g_{p}+k$, where $p$
is an odd prime and $k\leq\lambda p^{\theta}$, the maximum of the systole over
all closed $\textrm{CAT}(-1)$ surfaces $\bar{S}_{g}$ of genus $g$ is roughly
at least $\frac{1-\theta}{3}\log g$. But the maximum of the systole over all
closed $\textrm{CAT}(-1)$ surfaces of genus $g$ is attained by a hyperbolic
metric; see [10]. Hence the result. ∎

######  Remark 3.4.

Our argument applies to surfaces with sparser genera than the ones of
Brooks/Buser–Sarnak (e.g., for polynomial genera in primes in an arithmetic
progression). Let $a$ and $q$ be two coprime numbers. Suppose that
$\\{S_{g_{p}}\mid p\mbox{ prime with }p\equiv a\,(\mbox{mod }q)\\}$ is a
family of closed hyperbolic surfaces of genus $g_{p}$ with

$\operatorname{{\rm sys}}(S_{g_{p}})\geq C_{1}\,\log(g_{1})-C_{2}$

for some positive constants $C_{1}$ and $C_{2}$, where $g_{p}=P(p)$ is
polynomial in $p$ as in (1.2). We can apply the same arguments as in the proof
of Theorem 3.1 to construct a family of closed hyperbolic surfaces $S_{g}$ for
every genus $g$ with a logarithmic systolic growth using the prime gap theorem
version of [2, Theorem 3] (see the MathSciNet review of the article for a
statement without misprints) instead of Theorem 2.1.

## Acknowledgments

Mikhail Katz was supported by the BSF grant 2020124 and the ISF grant 743/22.
Stéphane Sabourau was supported by the ANR project Min-Max (ANR-19-CE40-0014).

## References

  * [1] Buser, P.; Sarnak, P. On the period matrix of a Riemann surface of large genus (with appendix by J.H. Conway and N.J.A. Sloane). _Invent. Math._ 117 (1994), no. 4, 27–56.
  * [2] Baker, R. C.; Harman, G.; Pintz, J. The exceptional set for Goldbach’s problem in short intervals. In Sieve methods, exponential sums, and their applications in number theory (Cardiff, 1995), 1–54. _London Math. Soc. Lecture Note Ser._ 237, Cambridge University Press, 1997.
  * [3] Baker, R. C.; Harman, G.; Pintz, J. The difference between consecutive primes. II _Proc. London Math. Soc. (3)_ 83 (2001), no. 3, 532–562.
  * [4] Brooks, R. Injectivity radius and low eigenvalues of hyperbolic manifolds. _J. Reine Angew. Math._ 390 (1988), 117–129.
  * [5] Brooks, R. Platonic surfaces. _Comment. Math. Helv._ 74 (1999), no. 1, 156–170.
  * [6] Buser, P. Geometry and spectra of compact Riemann surfaces. _Modern Birkhäuser Classics_ , Birkhäuser, 2010.
  * [7] Fay, J.D. Theta functions on Riemann surfaces. _Lect. Notes Math._ vol. 352, Springer, 1973.
  * [8] Goodwillie, T.; Hebda, J.; Katz, M. Extending Gromov’s optimal systolic inequality. _Journal of Geometry_ 114 (2023), article 23. https://doi.org/10.1007/s00022-023-00685-3, https://mathscinet.ams.org/mathscinet/article?mr=4624541
  * [9] Gromov, M. Systoles and intersystolic inequalities. Actes de la Table Ronde de Géométrie Différentielle (Luminy, 1992), 291–362, _Sémin. Congr._ 1, Soc. Math. France, Paris, 1996.
  * [10] Jabbour, A.; Sabourau, S. Sharp systolic bounds on negatively curved surfaces. _Math. Annalen_ 389 (2024), 1131–1159.
  * [11] Katz, M.; Sabourau, S. Nonpositively curved surfaces are Loewner (2024). Journal of Geometric Analysis (2024). https://doi.org/10.1007/s12220-024-01732-4, https://arxiv.org/abs/2404.00757
  * [12] Katz, M.; Schaps, M.; Vishne, U. Logarithmic growth of systole of arithmetic Riemann surfaces along congruence subgroups. _J. Differ. Geom._ 76 (2007), no. 3, 399–422.
  * [13] Liu, M.; Petri, B. Random surfaces with large systoles (2023). https://arxiv.org/abs/2312.11428
  * [14] Makisumi, S. A note on Riemann surfaces of large systole. _J. Ramanujan Math. Soc._ 28 (2013), 359–377.
  * [15] Parlier, H. The homology systole of hyperbolic Riemann surfaces. _Geom. Dedicata_ 157 (2012), 331–338.
  * [16] Parlier, H. Simple closed geodesics and the study of Teichmüller spaces. Handbook of Teichmüller theory. Vol. IV, 113–134, IRMA Lect. Math. Theor. Phys., 19, Eur. Math. Soc., Zürich, 2014.
  * [17] Petri, B. Hyperbolic surfaces with long systoles that form a pants decomposition. _Proc. Am. Math. Soc._ 146 (2018), no. 3, 1069–1081.
  * [18] Petri, B.; Walker, A. Graphs of large girth and surfaces of large systole. _Math. Res. Lett._ 25(2018), no. 6, 1937–1956.
  * [19] Schmutz, P. Riemann surfaces with shortest geodesic of maximal length. _Geom. Funct. Anal_. 3 (1993), no.6̃, 564–631.