Theoretical uncertainties in exclusive electroproduction of S-wave heavy quarkonia

Abstract

In this work, we revise the conventional description of $J/\psi (1S)$, $\Upsilon (1S)$, $\psi '(2S)$ and $\Upsilon '(2S)$ elastic photo- and electroproduction off a nucleon target within the color dipole picture and carefully study various sources of theoretical uncertainties in calculations of the corresponding electroproduction cross sections. For this purpose, we test the corresponding predictions using a bulk of available dipole cross section parametrisations obtained from deep inelastic scattering data at HERA. Specifically, we provide the detailed analysis of the energy and hard-scale dependencies of quarkonia yields employing the comprehensive treatment of the quarkonia wave functions in the Schrödinger equation based approach for a set of available $c-\bar{c}$ and $b-\bar{b}$ interquark interaction potentials. Besides, we quantify the effect of Melosh spin rotation, the $Q^2$-dependence of the diffractive slope and an uncertainty due to charm and bottom quark mass variations.

1 Introduction

One of the major widely used probes for interplay between hard and soft QCD physics is by means of bound states of heavy (charm and bottom) quarks known as quarkonia. Among these, the most well-studied are S-wave $J/\psi$, $\psi '$ and $\Upsilon$ states produced in high-energy particle collisions (for a detailed review on quarkonia physics, see e.g. Refs. [1,2,3]).

Despite a notable progress in theoretical description of heavy quarkonia production done over past few decades, the quarkonia production mechanism, as well as their propagation and dissociation in a hot medium, is considered to be an important probe for the medium created in heavy-ion collisions [4], and is still an actively developing research area. The problem concerns highly uncertain rates of $J/\psi$ and $\psi '$ mesons production in pp and pA collisions. These processes are also considered to be among the main tools for studying the soft QCD effects in hard processes.

The wealth of existing experimental data and theoretical studies show that the widely used simplifications in the analysis of exclusive quarkonia electroproduction observables may have significant impact on theoretical predictions and thus should be taken with care. One of the important ingredients of quarkonia production observables are the light-cone (LC) wave functions of heavy quarkonia. A popular simple model for the quarkonia wave functions is based upon an assumption that the potential between the bound heavy quarks is perfectly harmonic and no spin rotation in the quarkonia formation is considered (see e.g. Refs. [5, 6]). Such a treatment is usually performed in the conventional non-relativistic QCD (NRQCD) framework without an account for a non-trivial dependence on intrinsic transverse momenta and longitudinal momentum fractions of heavy quarks. In the case of charmonia production, a significance of non-perturbative and relativistic effects is often underestimated since the charm quark mass is not sufficiently large. Moreover, a spin rotation of heavy quark spinors from the $Q\bar{Q}$ rest frame to the infinite momentum frame known as Melosh rotation [7, 8], which influences mainly the angular part of the wave function, has a notable impact on (in particular, S-wave) charmonia differential observables [8, 9] while it is sometimes neglected in the existing calculations. A properly formulated radial part of the quarkonium wave function should be obtained by a numerical solution of the Schrödinger equation for a realistic $Q\bar{Q}$ potential and with an appropriate boosting and spin rotation. This work is aimed, in particular, at a proper accounting for and studying various sources of theoretical uncertainties in modeling the elastic $J/\psi$, $\psi '$ and $\Upsilon$ electroproduction processes in $\gamma p$ collisions in the color dipole picture [10].

One of the major problems of the QCD scattering theory is to correctly identify the underlying degrees of freedom which are the eigenstates of interaction. In last decades, it became clear that such eigenstates are color dipoles [5, 10,11,12], the universal elementary building blocks automatically accumulating both the hard and soft fluctuations [13, 14]. In particular, the light-cone (LC) color dipole framework has been developed and applied in treatment of both diffractive and inclusive quarkonia electro- and photoproduction processes in QCD in terms of certain superpositions of these eigenstates [5, 10]. The dipole picture has turned out to be rather successful in describing various hard processes beyond conventional QCD factorisation [15]. In particular, it is known to give as precise predictions e.g. for the Drell–Yan cross section as the Next-to-Leading-Order (NLO) collinear factorisation framework (see e.g. Ref. [16] and references therein). In deep inelastic scattering (DIS) or in vector meson production the virtual photon is expected to produce the quark–antiquark dipole with a transverse separation depending on the photon virtuality. The dipole formalism relies on a specific type of factorisation (or dipole factorisation) when a scattering cross section is written in impact parameter space as a convolution of the LC wave functions (for e.g. $\gamma ^{*}\rightarrow c\bar{c}$ and $c\bar{c}\rightarrow J/\psi$ fluctuations in the case of $J/\psi$ electroproduction) and the universal phenomenological dipole cross section fitted to the DIS data. One of the well-known coloured medium effects predicted by QCD is the so called colour transparency [17], an intrinsic feature of the dipole framework, when the medium becomes more transparent for smaller-size dipole configurations [11, 12, 18].

Traditionally, exclusive (or elastic) photo ($Q^2\simeq 0$) and electro ($Q^2\gg 0$) production of heavy quarkonia receives a lot of attention due to particularly clean signatures and precision measurements of the corresponding observables differential in hard scale, $Q^2$, energy, W, and momentum transfer squared, t. Such processes are highly relevant for e.g. a better understanding of the gluon density properties and their impact parameter profile in a target at very small x [19,20,21,22,23,24,25], as well as for probing the details of the quarkonia production and propagation processes. Indeed, the exclusive photo production cross sections are given by the gluon density in the second power, thus enabling us to probe it to a much better precision than in inclusive reactions whose differential cross sections are proportional to the first power of the (predominantly, gluon) parton density function (PDF). Our current study aims at a comprehensive analysis of these observables and related theory uncertainties in the dipole picture against the available data. The color dipole formalism [5, 10,11,12,13], well known for particularly successful description of various photo and hadro production reactions in both pp and pA collisions, enables to include systematically the QCD factorisation breaking and nuclear coherence effects as well as the initial-state and saturation phenomena [26, 27]. Consequently, there is a notable sensitivity of exclusive photoproduction observables to different gluon saturation models at low-x that is worth a careful study providing an efficient discriminating tool for various existing parametrisations for the low-x gluon density at a periphery of the target nucleon.

Fig. 1

A schematic illustration of the exclusive quarkonium electroproduction process, $\gamma ^{*}\,p\rightarrow V\,p$, in the dipole picture. On the left panel, the structure of the amplitude and kinematic variables in impact parameter space are depicted while its amplitude squared for the $J/\psi$ electroproduction is shown on the right panel

As a starting point, we would like to test various models of the LC wave functions with different $c-\bar{c}$ and $b-\bar{b}$ interquark interaction potentials and phenomenological dipole parametrisations against the recent data on $J/\psi$, $\psi '$ and $\Upsilon$ photo- and electroproduction as well as to study the impact of Melosh spin rotation in these observables. Of particular interest for us is the study of $\psi '$ production observables which are known to be highly sensitive to the shape of the wave function (in particular, to the position of its node) than $J/\psi$ observables [8]. Besides, we would like to estimate an overall theoretical uncertainty in the exclusive cross sections due to poorly known gluon density at low-x. This is accounted for by using several representative dipole parameterisations widely used in the literature. Finally, such effects as sensitivity to the heavy quark mass, to the skewness in the gluon density, and to the diffractive slope parameterisation are quantified.

The paper is organised as follows. In Sect. 2, we provide the basic details of the dipole approach to the exclusive quarkonia electroproduction with proper definitions of the kinematical variables, the elastic amplitude and the cross section. A thorough description of the normalised LC quarkonia wave functions, the boosting procedure and an overview of the interquark potential models used in our analysis is given in Sect. 3. Further on, in Sect. 4, we present a brief overview of the most frequently used dipole parameterisations that will be employed in our numerical analysis for estimation of the underlined uncertainties in QCD modeling of the target gluon density evolution at small-x. In Sect. 5 we compare our model predictions for the photo- and electroproduction cross sections of different quarkonia with available data and analyze the theoretical uncertainties caused by various sources and ingredients coming into the color dipole formalism. Final remarks and conclusions are summarized in Sect. 6.

2 Exclusive quarkonia electroproduction: dipole picture

In the framework of color dipole approach [5, 10,11,12,13], the projectile (real or virtual, with $q^2=-Q^2$) photon undergoes strong interactions via its Fock components containing quarks and gluons with the target proton in the frame where the target proton is at rest. In the dipole picture, such interactions are described by the universal dipole cross section, which is not derivable from the first principles but, instead, is fitted to e.g. HERA data (for more details, see below). In the case of exclusive vector meson electroproduction illustrated in Fig. 1 (left panel), such a lowest Fock state corresponds to the $Q\bar{Q}$ dipole whose transverse size r is nearly frozen in the high energy limit. Once the dipole scattering occurs, a coherent $Q\bar{Q}$ state forms a vector meson by means of a projection of the $Q\bar{Q}$ production amplitude on to a given LC quarkonium wave function. Let us now briefly describe the main ingredients of the dipole formulation of this process.

The forward amplitude for exclusive electroproduction of a vector meson V with mass $M_V$ in the target rest frame is given by (see e.g. Ref. [8] and references therein)

$$\begin{aligned}&\mathrm {Im}\mathcal {A}^{\gamma ^{*}p\rightarrow Vp}_{T,L}(x,Q^2)\nonumber \\&\quad =\int \mathrm {d}^{2}r\int \limits _{0}^{1}\mathrm {d}z \,\Psi ^{\dagger }_{V}(r,z)\, \Psi _{\gamma ^{*}_{T,L}}(r,z;Q^2)\sigma _{q\bar{q}}(x,r) , \;\; x\nonumber \\&\quad =\frac{M_V^2+Q^2}{s} , \end{aligned}$$

(2.1)

where x is the standard Bjorken variable [19], $s=Q^2+W^2$ is the square of the ep center-of-mass energy (with W being the $\gamma ^{*}p$ center-of-mass energy), $\Psi _{V}(r,z)$ is the vector meson V wave function, $\Psi _{\gamma ^{*}_{T,L}}(r,z;Q^2)$ is the LC distribution (or wave) function of a transversely (T) or longitudinally (L) polarized virtual photon for a $Q\bar{Q}$ fluctuation, $\vec {r}$ is the transverse size of the $Q\bar{Q}$ dipole, and $z=p_Q^+/p_{\gamma }^+$ is the boost-invariant fraction of the photon momentum $p_{\gamma }^+=E_{\gamma }+p_{\gamma }$ carried by a heavy quark (or anti-quark). The universal dipole cross section $\sigma _{q\bar{q}}(x,r)$ describes the dipole scattering off the target. It is typically fitted to the precision inclusive DIS data at HERA and then is used to describe a variety of other processes in ep and pp collisions (see below). In the NRQCD limit, one neglects relative motion of Q and $\bar{Q}$ such that $z=1/2$, and the LC wave function reduces to $\Psi _{V}(r,z)\propto \delta (z-1/2)$ [11]. In what follows, we go beyond this approximation.

The perturbative LC $\gamma ^{*} \rightarrow Q\bar{Q}$ wave function is given by [28, 29]

$$\begin{aligned} \Psi ^{(\mu ,\bar{\mu })}_{\gamma ^{*}_{T,L}}(r,z;Q^2)= & {} \frac{\sqrt{N_c\alpha _\mathrm{em}}}{2\pi }Z_Q\,\chi _Q^{\mu \dag } \hat{\mathcal {O}}_{T,L}\tilde{\chi }_{\bar{Q}}^{\mu }\,K_0(\varepsilon r) ,\nonumber \\ \varepsilon ^2= & {} z(1-z)Q^2+m_Q^2 , \end{aligned}$$

(2.2)

where $\varepsilon$ and $Z_Q$ are the energy and the electric charge of the heavy quark ($Z_c=2/3$, $Z_b=1/3$), $\chi ^\mu _Q$ and $\tilde{\chi }^{\bar{\mu }}_{\bar{Q}}\equiv i\sigma _y\chi ^{\bar{\mu } *}_{\bar{Q}}$ are the two-component spinors in the infinite momentum frame normalized as follows [30],

$$\begin{aligned} \sum \limits _{\mu ,\bar{\mu }}\left( \chi ^{\mu \dag }_Q\hat{A} \tilde{\chi }^{\bar{\mu }}_{\bar{Q}}\right) ^*\left( \chi ^{\mu \dag }_Q \hat{B}\tilde{\chi }^{\bar{\mu }}_{\bar{Q}}\right) = \mathrm {Tr}(\hat{A}^{\dagger } \hat{B}), \end{aligned}$$

(2.3)

and $K_0(\varepsilon r)$ is the modified Bessel function of the second kind. The operators $\hat{\mathcal {O}}_{T,L}$ are defined as follows,

$$\begin{aligned} \mathcal {O}_T&= m_Q\vec {\sigma }\cdot \vec {e}_\gamma + i(1-2z)(\vec {\sigma }\cdot \vec {n})(\vec {e}_\gamma \cdot \vec {\nabla }_r) +(\vec {n}\times \vec {e}_\gamma )\vec {\nabla }_r , \nonumber \\ \mathcal {O}_L&= 2Qz(1-z)\vec {\sigma }\cdot \vec {n} , \quad \vec {\sigma }=(\sigma _x,\sigma _y,\sigma _z) , \quad \vec {\nabla }_r \equiv \partial /\partial \vec {r} , \end{aligned}$$

(2.4)

where $\vec {e}_\gamma$ is the transverse photon polarisation vector, $\vec {n}=\vec {p}_{\gamma }/|\vec {p}_{\gamma }|$ is a unit vector along the photon momentum, and $\sigma _{x,y,z}$ are the Pauli matrices. In what follows, following Ref. [8] we neglect the effects associated with non-perturbative interactions within the heavy quark pair (including charm quarks) and that are not included into the corresponding interaction potential since $m_Q$ provides a sufficiently perturbative scale even in the photoproduction limit $Q^2 \rightarrow 0$.

The quarkonium wave function is properly defined only in the $Q\bar{Q}$ rest frame where it can be directly found by solving the Schrödinger equation. Below, we discuss such solutions for several distinct interquark potentials. In order to obtain the production amplitude given by Eq. (2.1), the quarkonium wave function should be found in the infinite momentum frame. For a classical $Q\bar{Q}$ configuration, such a wave function could be computed from that in the $Q\bar{Q}$ rest frame by a applying a Lorentz boost. The quantum effects, however, are relevant such that a tower of Fock states emerges as a result of such a transformation [8], and the lowest Fock $|Q\bar{Q}\rangle$ components in these frames do not represent the same configuration. This issue is further discussed in the next Section.

In what follows, we assume a simple factorization of spatial and spin-dependent parts of the vector meson V wave function such as

$$\begin{aligned} \Psi _{V}^{(\mu ,\bar{\mu })}(z,\vec {p}_T)=U^{(\mu ,\bar{\mu })}(z, \vec {p}_T)\Psi _{V}(z,p_T) , \end{aligned}$$

(2.5)

where

$$\begin{aligned} U^{(\mu ,\bar{\mu })}(z,\vec {p}_T)=\frac{1}{\sqrt{2}}\xi _Q^{\mu \dag }\vec {\sigma } \vec {e}_{V}\tilde{\xi }_{\bar{Q}}^{\bar{\mu }}, \quad \tilde{\xi }_{\bar{Q}}^{\bar{\mu }}=i\sigma _y \xi _{\bar{Q}}^{\bar{\mu }*} , \end{aligned}$$

(2.6)

in terms of the vector meson polarisation vector $\vec {e}_{V}$ and quark spinors $\xi$ in the meson rest frame. The latter are related to spinors $\chi$ in the infinite momentum frame as follows

$$\begin{aligned} \xi ^\mu _Q = R(z,\vec {p}_T)\chi _Q^\mu , \quad \xi _{\bar{Q}}^{\bar{\mu }}=R(1-z,-\vec {p}_T)\chi _{\bar{Q}}^{\bar{\mu }} , \end{aligned}$$

(2.7)

known as the Melosh spin rotation [7, 8]. The R-matrix of such a rotation is given by

$$\begin{aligned} R(z,\vec {p}_T)=\frac{m_Q + zM_V - i(\vec {\sigma }\times \vec {n}) \vec {p}_T}{\sqrt{(m_Q+zM_V)^2+p_T^2}} . \end{aligned}$$

(2.8)

Using the quarkonium wave function given by Eq. (2.5) we assume that the vertex $\psi \rightarrow c\bar{c}$ differs from the photon-like $\gamma ^{*} \rightarrow c\bar{c}$ vertex with the structure $\psi _\mu \bar{u} \gamma ^\mu u$ used in Refs. [5, 6, 31, 32]. As was noticed in Ref. [8], the latter in the $c\bar{c}$ rest frame contains both S- and D-wave states whereas the D-wave weight is correlated strongly with the structure of the vertex and cannot be proved by any reasonable nonrelativistic $c-\bar{c}$ interaction potential.

Substituting

$$\begin{aligned}&\tilde{\xi }_{\bar{Q}}^{\bar{\mu }}=i\sigma _yR^*(1-z,-\vec {p}_T)(-i) \sigma _y^{-1}\tilde{\chi }_{\bar{Q}}^{\bar{\mu }} , \nonumber \\&\xi _Q^{\mu \dag }=\chi _Q^{\mu \dag }R^{\dag }(z,\vec {p}_T) , \end{aligned}$$

(2.9)

into Eq. (2.6) one gets finally

$$\begin{aligned}&U^{(\mu ,\bar{\mu })}(z,\vec {p}_T)\nonumber \\&\quad =\frac{1}{\sqrt{2}}\,\chi _Q^{\mu \dag }R^{\dagger }(z,\vec {p}_T)\, \vec {\sigma }\cdot \vec {e}_{V}\,\sigma _y R^*(1-z,-\vec {p}_T)\, \sigma _y^{-1}\tilde{\chi }_{\bar{Q}}^{\bar{\mu }} . \end{aligned}$$

(2.10)

The resulting dipole formula for the amplitude of photo and electroproduction of heavy quarkonia reads

$$\begin{aligned} \mathrm {Im}\mathcal {A}^{\gamma ^{*}p\rightarrow Vp}_{T,L}(x,Q^2)\! =\!\int \limits _0^1\mathrm {d}z \int \!\mathrm {d}^2r\, \Sigma _{T,L}(z,\vec {r};Q^2) \, \sigma _{q\bar{q}}(x,r) \end{aligned}$$

(2.11)

where

$$\begin{aligned} \Sigma _{T,L}(z,\vec {r};Q^2)= & {} \int \frac{\mathrm {d}^2 p_T}{2\pi } e^{-i\vec {p}_T\vec {r}}\,\Psi _{V}(z,p_T)\nonumber \\&\times \sum \limits _{\mu ,\bar{\mu }}U^{\dagger (\mu ,\bar{\mu })} (z,\vec {p}_T)\, \Psi ^{(\mu ,\bar{\mu })}_{\gamma ^{*}_{T,L}}(r,z;Q^2).\nonumber \\ \end{aligned}$$

(2.12)

The T and L amplitudes in a more explicit form are shown in Appendix C.

The total $\gamma ^{*}p\rightarrow Vp$ cross section is conventionally represented as sum of T and L contributions [8]

$$\begin{aligned}&\sigma ^{\gamma ^{*}p\rightarrow Vp}(x,Q^2)\nonumber \\&\quad =\sigma _T^{\gamma ^{*}p\rightarrow Vp}(x,Q^2) +\tilde{\varepsilon }\,\sigma _L^{\gamma ^{*}p \rightarrow Vp}(x,Q^2) = \frac{1}{16\pi B}\nonumber \\&\qquad \times \left( \Big \vert \mathrm {Im} \mathcal {A}^{\gamma ^{*}p\rightarrow Vp}_{T}(x,Q^2)\Big \vert ^{2}\! +\!\tilde{\varepsilon }\Big \vert \mathrm {Im}\mathcal {A}^{\gamma ^{*}p \rightarrow Vp}_{L}(x,Q^2)\Big \vert ^{2}\right) , \end{aligned}$$

(2.13)

with $\tilde{\varepsilon } = 0.99$. Here, B is the slope parameter fitted to the exclusive quarkonia electroproduction data. In the energy-independent approximation [33] it is taken to be $B=4.73\ \hbox {GeV}^{-2}$, while its possible energy and $Q^2$ dependence is discussed in more detail in Sects. 5.1.1 and 5.1.2.

Derivation of above formulas relies on the assumption that the S-matrix is purely real and so the amplitude $\mathcal {A}$ is purely imaginary. Following Refs. [5, 34, 35], the real part of the amplitude can be accounted for by multiplying the cross section $\sigma ^{\gamma ^{*}p\rightarrow Vp}_{T,L}(x,Q^2)$ by a factor $1+\tan ^2(\pi \lambda _{T,L}/2)$, where

$$\begin{aligned} \lambda _{T,L}&=\Bigg |\frac{\partial \ln \mathrm {Im} \mathcal {A}^{\gamma ^{*}p\rightarrow Vp}_{T,L}(x,Q^2)}{\partial \ln x}\Bigg |\nonumber \\&=\Bigg |\frac{1}{\mathrm {Im} \mathcal {A}^{\gamma ^{*}p\rightarrow Vp}_{T,L}} \int \limits _0^1\mathrm {d}z \int \mathrm {d}^2r \Sigma _{T,L}(z,\vec {r};Q^2) \frac{\partial \sigma _{q\bar{q}}(r,x)}{\partial \ln x}\Bigg |,\quad \end{aligned}$$

(2.14)

provided that only the dipole cross section depends on the Bjorken x.

3 Light-cone quarkonia wave function

One yet missing ingredient of the formula (2.1) is the LC quarkonium wave function $\Psi _{V}(r,z)$. Alike the LC photon-quark wave function $\Psi ^{T,L}_{\gamma ^{*}}(r,z;Q^2)$, it is defined in the infinite momentum frame. Let us discuss if and how this object can be obtained from the first principles.

In the $Q\bar{Q}$ rest frame and in impact parameter representation, the quarkonia wave function is typically found by solving the Schrödinger equation for a given choice of the heavy quark interaction potentials as discussed in Appendix B. In this work, we have employed five well-known parametrisations for the heavy quark interaction potentials illustrated in Fig. 2 for $c-\bar{c}$ (left panel) and $b-\bar{b}$ (right panel) cases and described in detail in Appendix A.

Fig. 2

An illustration of five distinct $c-\bar{c}$ (left panel) and $b-\bar{b}$ (right panel) interaction potentials, used in this work, as functions of 3D interquark distance ${\tilde{r}}$. For a detailed description and characteristic parameters of these potentials, see Appendix A

Since, in general, there is no direct relation between the rest-frame wave function of the lowest Fock $|Q\bar{Q}\rangle$ component and that in the infinite-momentum frame, the problem of building the latter is rather difficult and there is no generally acceptable solution yet. In the literature, there are recipes towards finding such a wave function, and in what follows we employ one particular widely used recipe of Ref. [36] known to give rather accurate predictions in the relevant kinematical regions (cf. Ref. [37]).

For practical purposes, it is convenient to turn to the momentum-space wave function,

$$\begin{aligned} \psi (p)=\frac{2}{\sqrt{2\pi }p}\int \limits _0^{\infty } \mathrm {d}{\tilde{r}}\,{\tilde{r}}\,\psi ({\tilde{r}})\sin (p\,{\tilde{r}}) , \quad \int |\psi (p)|^2\,\mathrm {d}^3 p=1 , \end{aligned}$$

(3.1)

in terms of the quark 3-momentum $p\equiv |\vec {p}\,|$ and the 3D interquark distance, ${\tilde{r}}\equiv |\vec {{\tilde{r}}}|$. Since the quarkonium production amplitude (2.1) is written in the infinite momentum frame, the corresponding wave function $\psi (p)$ should first be appropriately boosted to that frame. In terms of the LC variables, z and $p_T$, the invariant mass squared of the $Q\bar{Q}$ pair reads

$$\begin{aligned} M_{Q\bar{Q}}^2 = \frac{p_T^2+m_Q^2}{z(1-z)} , \end{aligned}$$

(3.2)

while the same quantity in the $Q\bar{Q}$ rest frame is given by

$$\begin{aligned} M_{Q\bar{Q}}^2 = 4(p^2+m_Q^2), \quad p^2 = p_L^2 + p_T^2 , \end{aligned}$$

(3.3)

where $p_L$ is the longitudinal component of the quark 3-momentum $\vec {p}$. These two relations, therefore, yield

$$\begin{aligned} p^2=\frac{p_T^2+(1-2z)^2m_Q^2}{4z(1-z)}, \quad p_L^2=\frac{(p_T^2+m_Q^2)(1-2z)^2}{4z(1-z)} , \end{aligned}$$

(3.4)

providing an appropriate conversion of kinematical variables between the infinite momentum and $Q\bar{Q}$ rest frames. Besides, following the recipe of Ref. [36] the conservation of probability density upon such a boosting

$$\begin{aligned} \mathrm {d}^3 p |\psi (p)|^2 = \mathrm {d}^2 p_T\mathrm {d}z|\psi (p_T,z)|^2 , \quad \mathrm {d}^3 p = \mathrm {d}p_L\mathrm {d}^2p_T \end{aligned}$$

(3.5)

results in the following Terent’ev relation [36] between the LC wave function $\psi (p_T,z)$ and its counterpart in the target rest frame $\psi (p)$

$$\begin{aligned} \psi (p_T,z)= & {} \left( \frac{p_T^2+m_Q^2}{16(z(1-z))^3}\right) ^{\frac{1}{4}}\,\psi (p) ,\nonumber \\&\times \,\int |\psi (p_T,z)|^2\mathrm {d}^2 p_T\mathrm {d}z = 1 , \end{aligned}$$

(3.6)

where $p=p(p_T,z)$ is given by Eq. (3.4). Note, the Terent’ev prescription for the Lorentz boosting presented above has been discussed and compared with exact calculations using the sophisticated Green function approach in Ref. [37]. It has been shown that for the $J/\psi$ wave function the Terent’ev prescription gives very accurate results for the averaged $\langle z \rangle \sim 0.5$. The LC wave function in the impact parameter representation is then given by

$$\begin{aligned} \Psi _{V}(r,z)=\int \limits _0^{\infty }\mathrm {d}p_T \, p_T J_0(p_Tr)\,\psi (p_T,z). \end{aligned}$$

(3.7)

Fig. 3

The LC wave function $\Psi _{V}(r,z)$ for $J/\psi (1S)$ (left panel) and $\psi '(2S)$ (right panel) mesons for different quark momentum fractions z. The distribution function $\Psi _{V}(r,z)$ is generated by two models for the $c-\bar{c}$ interaction potential: harmonic oscillator model denoted as HAR (green dotted lines) and Buchmüller-Tye parameterisation, or BT (black solid lines)

Fig. 4

The LC wave function $\Psi _{V}(r,z)$ for $\Upsilon (1S)$ (left panel), $\Upsilon '(2S)$ (middle panel) and $\Upsilon ''(3S)$ (right panel) mesons for different quark momentum fractions z. The distribution function $\Psi _{V}(r,z)$ is generated by two models for the $b-\bar{b}$ interaction potential: harmonic oscillator model denoted as HAR (green dotted lines) and Buchmüller-Tye parameterisation, or BT (black solid lines)

Fig. 5

Comparison of different parametrisations for the dipole cross section used in our calculations as described in the text

In Fig. 3 we show the numerical results for the boosted LC wave functions $\Psi _{V}(r,z)$ for two charmonia states, $J/\psi (1S)$ and $\psi '(2S)$, and which are obtained starting from numerical solutions of the Schrödinger equation for two models of the $c-\bar{c}$ interaction potential – the harmonic oscillator and Buchmüller-Tye parametrisation, discussed in Appendix A. While the overall shape of the wave functions appears to be consistent between the two models, they yield notable quantitative differences, especially for $\psi '(2S)$, where the positions of the node and the minimum are quite sensitive to the choice of the potential. The corresponding LC wave functions for $\Upsilon (1S)$, $\Upsilon '(2S)$ and $\Upsilon ''(3S)$ mesons are shown in Fig. 4.

We have also performed calculations of the wave functions and total elastic electroproduction cross sections for a number of different $c\bar{c}$ (1S and 2S) and $b\bar{b}$ (1S, 2S and 3S) vector meson states for five distinct parameterisations of the interquark potentials and five different parameterisations for the dipole cross sections, $\sigma _{q\bar{q}}(r,x)$. Since the number of possible combinations of the parameterisations and states can be rather extensive, as a part of this project, we created a webpage on https://hep.fjfi.cvut.cz/vm.php, where one can find numerical datasets (grids) for each vector meson state (Figs. 3 and 4), interquark potential (Fig. 2) and the dipole parameterisation (Fig. 5).

The datasets are available for vector meson wave functions in the forms of a 3D radial solution of the Schrödinger equation in the $Q\bar{Q}$ rest frame, $\psi (\tilde{r})$ [shown in Eq. (B6)], the boosted LC wave function in momentum space, $\psi (p_T,z)$, given by Eq. (3.6), and the boosted LC wave function in impact parameter space, $\psi (r,z)$, given by Eq. (3.7). An interpolating routine written in C++ (including also an example for calculations) is also available on the same webpage. The web interface enables to generate plots for the electroproduction cross sections for a selected combination of the quarkonium wave function generated by the explicit $Q-\bar{Q}$ potential with the explicit dipole model for $\sigma _{q\bar{q}}(r,x)$. Calculations can be performed including or neglecting the Melosh spin rotation effects. Numerical results can be presented also in the form of a table, which can be used for practical purposes.

4 Dipole cross section

The essential ingredient of the color dipole approach, the universal dipole cross section $\sigma _{q\bar{q}}(r,x)$, has been first introduced a long ago in Ref. [13]. During past three decades, its kinematic (and energy) dependence has undergone remarkable development largely promoted by precision experimental information from the HERA collider. While an exact theoretical modelling of the dipole cross section (and the corresponding partial dipole amplitude) is not nearly close to its complete understanding, a number of phenomenological parametrisations accounting for the saturation phenomenon and the QCD-inspired Bjorken x- and hard-scale evolution have been proposed in the literature [26, 38,39,40,41,42,43,44,45,46,47,48,49]) and rely on the fits to the HERA DIS data.

Introducing the partial dipole amplitude ${\mathcal {N}} (\vec {r}, x, \vec {b})$, one conventionally determines the universal dipole cross section $\sigma _{q\bar{q}}(r,x)$ as an integral over the impact parameter $\vec {b}$:

$$\begin{aligned} \sigma _{q\bar{q}}(r,x) = 2\,\int d^2b\, {\mathcal {N}} (\vec {r}, x, \vec {b}), \quad r=|\vec {r}| . \end{aligned}$$

(4.1)

The evolution in x- or in rapidity $Y = \ln (1/x)$ in the high-energy ($x\ll 1$) regime is treated e.g. by an infinite hierarchy of the so-called Balitsky-JIMWLK equations [50,51,52,53,54,55] in the framework of the Color Glass Condensate (CGC) formalism [56, 57]. These equations reduce to the Balitsky–Kovchegov (BK) equation [54, 58] in the mean-field approximation. As it is rather difficult to obtain the $\vec {b}$-dependent solutions of the BK equation [59] while the impact-parameter profile is determined by essentially non-perturbative QCD phenomena, one usually imposes such approximations as the translational invariance of the amplitude disregarding the impact parameter dependence in numerical solutions. An alternative way usually admitted in the literature is to consider more phenomenological models for the $\vec {b}$ dependence, as well as accounting for the saturation phenomenon and the hard-scale evolution via DGLAP, that are fit to precision data e.g. from HERA. A naive comparison of the predictions of the dipole model calculations using several distinct parametrisations for the universal dipole cross section is a standard way to roughly estimate the associated theoretical uncertainties.

Since long ago, it was understood that at small dipole separations r the dipole cross section is essentially proportional to the gluon PDF in the target [60,61,62], i.e.

$$\begin{aligned} \sigma _{q\bar{q}}(r,x)\simeq \frac{\pi ^2}{3}\alpha _s \Big (\frac{\Lambda }{r^2}\Big )\,r^2\,xg\Big (x,\frac{\Lambda }{r^2}\Big ), \end{aligned}$$

(4.2)

with $\Lambda \approx 10$ [63]. Later on, in Ref. [45] it was suggested to merge this asymptotics with a naive saturated ansatz for the dipole cross section. Later on, in Ref. [42] it was proposed to introduce explicitly the $\vec {b}$-dependence into the corresponding ansatz for the partial dipole amplitude. The latter yields the following widely used parametrisation enabling for a description of exclusive observables at HERA and is known as the IP-Sat model,

$$\begin{aligned} \sigma _{q\bar{q}}(r, x)&= 2\,\int d^2b\,\nonumber \\&\quad \times \left[ 1-\exp \left( -\frac{\pi ^2}{2N_c}\,r^2\,\alpha _s(\mu ^2) \, xg(x,\mu ^2)T_G(b)\right) \right] \nonumber \\ \end{aligned}$$

(4.3)

given in terms of the number of colors in QCD, $N_c=3$, the strong coupling constant $\alpha _s(\mu ^2)$ determined at the hard scale $\mu$ connected to the size of the dipole r in a simple way as $\mu ^2=C/r^2 + \mu _0^2$. Here, the model parameters C, $\mu _0$ and $\sigma _0$ are extracted by fitting to the HERA data. Besides, the gluon PDF in the target nucleon $xg(x, \mu ^2)$ at small x is found as a solution of the conventional DGLAP equation [64,65,66] which takes into account the gluon splitting function $P_{gg}(z)$ only,

$$\begin{aligned} \frac{\partial xg(x,\mu ^2)}{\partial \ln \mu ^2 } = \frac{\alpha _s(\mu ^2)}{2\pi } \int _x^1 dz\, P_{gg}(z) \frac{x}{z} g\Big (\frac{x}{z}, \mu ^2\Big ). \end{aligned}$$

(4.4)

Here, the starting value of the target gluon density at $\mu ^2=\mu _0^2$ is given by

$$\begin{aligned} xg(x,\mu _0^2) = A_g x^{-\lambda _g} (1-x)^{5.6}. \end{aligned}$$

(4.5)

The $\vec {b}$-dependence is accounted for by means of a simple Gaussian profile

$$\begin{aligned} T_G(b)=\frac{1}{2\pi B_G}\,\exp \left( -\frac{b^2}{2 B_G}\right) , \end{aligned}$$

(4.6)

where $B_G$ is an additional free parameter in the model that can be extracted, in particular, from the measured t-dependent elastic electron–proton scattering. In general, $B_G$ in the IP-Sat model is taken to be different from the slope of the exclusive quarkonia electroproduction cross section defined in Eq. (2.13) (see e.g. a discussion in Ref. [40]). A comprehensive IP-Sat model fit of the complete (inelastic and elastic) set of HERA data has been performed in Ref. [43] leading to

$$\begin{aligned}&A_g = 2.373, \quad \lambda _g = 0.052,\nonumber \\&\mu _0^2 =1.428\,\mathrm{GeV}^2, \quad B_G = 4.0\,\mathrm{GeV}^2, \quad C = 4.0 , \end{aligned}$$

(4.7)

The Golec–Biernat–Wusthoff (GBW) model [26],

$$\begin{aligned} \sigma _{q\bar{q}}(r,x)= & {} \sigma _0 \,\left( 1 - e^{-\frac{r^2Q_s^2(x)}{4}}\right) ,\nonumber \\ Q_s^2(x)\equiv & {} R_0^{-2}(x) = Q_0^2 \left( \frac{x_0}{x} \right) ^\lambda , \end{aligned}$$

(4.8)

where $Q_s(x)$ is the x-dependent (and $\mu$-independent) saturation scale, gives rise to a fairly good description of a large variety of observables in high-energy ep and pp collisions. The same conclusion concerns also nuclear targets as well as for both inclusive and exclusive final states at not-so-large transverse momenta (or $Q^2$) and small $x\lesssim 0.01$. This model resembles the Glauber model of multiple interactions and can also be straightforwardly used to incorporate the saturation effects. The global of the DIS HERA data accounting for the charm quark contribution provides different sets of parameters. We use two sets - the one taken from [27] we label as GBWold and the one from [40] we label as GBWnew

$$\begin{aligned}&{\mathrm {GBWold{:}}}\ Q_0^2 = 1\, \mathrm{GeV}^2, \quad x_0 = 3.04 \times 10^{-4}, \nonumber \\&\quad \lambda = 0.288, \quad \sigma _0 = 23.03\, \mathrm{mb} \nonumber \\&{\mathrm {GBWnew{:}}}\ Q_0^2 = 1\, \mathrm{GeV}^2, \quad x_0 = 1.11 \times 10^{-4},\nonumber \\&\quad \lambda = 0.287, \quad \sigma _0 = 23.9\, \mathrm{mb} . \end{aligned}$$

(4.9)

Following the tradition and for the sake of completeness, we use this simple model as a reference in comparison with other more complicated parametrisations. Besides, we will also consider the solution to the running coupling BK equation calculated according to the procedure in Ref. [67]. The BK equation describes the evolution in rapidity Y of the scattering amplitude ${\mathcal {N}} (\vec {r}, x)$. This formulation is based on the work of [68,69,70].

$$\begin{aligned} \frac{\partial N(\vec {r},x)}{\partial Y}= & {} \int d \vec {r}_1 K(\vec {r},\vec {r}_1,\vec {r}_2)\Bigg (N(\vec {r}_1,x)+N(\vec {r}_2,x)\nonumber \\&-N(\vec {r},x)-N(\vec {r}_1,x)N(\vec {r}_2,x)\Bigg ) \end{aligned}$$

(4.10)

where $\vec {r}_2 = \vec {r}-\vec {r}_1$. The kernel incorporating the running of the QCD coupling [70] is given by

$$\begin{aligned} K(\vec {r},\vec {r}_1,\vec {r}_2)= & {} \frac{\alpha _s(r^2)N_c}{2\pi ^{2}} \left( \frac{r^2}{r_1^2 r_2^2} +\frac{1}{r_1^2} \left( \frac{\alpha _s(r_1^2)}{\alpha _s(r_2^2)}-1\right) \right. \nonumber \\&\left. +\frac{1}{r_2^2} \left( \frac{\alpha _s(r_2^2)}{\alpha _s(r_1^2)}-1\right) \right) , \end{aligned}$$

(4.11)

with

$$\begin{aligned} \alpha _{s}(r^2)=\frac{4\pi }{(11-\frac{2}{3}N_f) \ln \left( \frac{4C^2}{r^2\Lambda ^2_\mathrm{QCD}}\right) } \end{aligned}$$

(4.12)

where $N_f$ is the number of active flavours and C is a parameter to be fixed from data. We use the fixed number of flavours scheme with $\Lambda _\mathrm{QCD}=0.241$ MeV. For the initial form of the dipole scattering amplitude the McLerran–Venugopalan (MV) model [71] is used:

$$\begin{aligned} N(\vec {r},x=0.01)=1-\exp \Bigg (-\frac{\left( r^2Q^2_{s0}\right) ^\gamma }{4} \ln \left( \frac{1}{r\Lambda _\mathrm{QCD}}+e\right) \Bigg ) \end{aligned}$$

(4.13)

with the values of the parameters $Q^2_{s0}$, C and $\gamma$ taken from [70] yielding $\sigma _0=32.895$ mb, $Q^2_{s0} = 0.165\ \hbox {GeV}^2$, $\gamma =1.135$ and $C = 2.52$. The fit was performed under the assumption that $\alpha _s(r^2)$ freezes for values of r larger than $r_0$ defined by $\alpha _s(r^2_0)=0.7$. This model is denoted below as rcBK.

We have also included the collinearly improved kernel [72] given by

$$\begin{aligned} K(\vec {r},\vec {r}_1,\vec {r}_2)= & {} \frac{\bar{\alpha }_s N_c}{2\pi ^{2}}\Bigg (\frac{r^2}{r_1^2 r_2^2} \left( \frac{r^2}{\mathrm {min}(r_1^2,r_2^2)}\right) ^{\pm \bar{\alpha }_s A_1}\nonumber \\&\times \frac{J_1(2\sqrt{\bar{\alpha }_s | \ln (r_1^{2}/r^{2})\ln (r_2^{2}/r^{2})|})}{\sqrt{\bar{\alpha }_s |\ln (r_1^{2}/r^{2})\ln (r_2^{2}/r^{2})|}}\Bigg ), \end{aligned}$$

(4.14)

with $A_1=11/12$, the positive sign refers to the situation where $r < \mathrm {min}(r_1,r_2)$ and $\bar{\alpha }_s=\alpha _s (\mathrm {min}(r^{2},r_1^{2},r_2^{2})) \frac{N_c}{\pi }$. This kernel was used with variable number of flavours scheme [70], each flavour has its $\Lambda _\mathrm{QCD}$ calculated from the recurrent relation

$$\begin{aligned} \Lambda _{N_f-1}=m_f^{1-\frac{\beta _{N_f}}{\beta _{N_f-1}}} \Lambda _{N_f}^{\frac{\beta _{N_f}}{\beta _{N_f-1}}}, \end{aligned}$$

(4.15)

where $\beta _{N_f}=(11N_c-2N_f)/3$ and $m_f$ is the mass of the quark of flavour f. As a starting point one can take measured $\alpha _s(r^{2}=4C^{2}/M_Z^{2})=0.1189$ for $n_f=5$ at a scale of Z boson mass $M_Z=91.187$ GeV. This leads to the formula

$$\begin{aligned} \Lambda _{5}=M_Ze^{-\frac{2\pi }{\alpha _s(r^{2}=4C^{2}/M_Z^{2})\beta _{5}}}. \end{aligned}$$

(4.16)

A collinear version of MV initial conditions was published in [72]

$$\begin{aligned}&N(\vec {r},x=0.01)\nonumber \\&\quad =\left[ 1-\exp \left( -\left[ \frac{r^2Q^2_{s0}}{4} \bar{\alpha }_s(r^{2}) \left( 1+\ln \left( \frac{\alpha _{sat}}{\bar{\alpha }_s(r^{2})}\right) \right) \right] ^{p}\right) \right] ^{1/p}, \end{aligned}$$

(4.17)

where $\bar{\alpha }_{sat}=\frac{N_c}{\pi }\alpha _{sat}, \bar{\alpha }_{s}(r^{2})=\frac{N_c}{\pi }\alpha _{s}(r^{2})$ and $\alpha _{sat}$ is fixed to 1. Parameters were fitted in [72] with $\sigma _0=31.4055$ mb, $Q^2_{s0}=0.4 \ \hbox {GeV}^{2}$, $C=2.586$ and $p=0.807$. This model is denoted as colBK. The dipole cross section is obtained from the scattering amplitude as $\sigma _{q\bar{q}}(r,x) =\sigma _0\,{\mathcal {N}}(r,x)$, where the normalisation $\sigma _0$ is fitted to the HERA data.

When decreasing the hard scale $Q^2\rightarrow \Lambda _\mathrm{QCD}^2$ relevant for e.g. photoproduction, one may reach small x values even for moderate and low energies. As was argued e.g. in Ref. [15], the Bjorken variable x becomes inappropriate in the soft limit. For such processes as e.g. the pion–proton scattering as well as the diffractive Drell–Yan and gluon radiation processes the saturation scale $Q_s^2 > rsim Q^2$ becomes a function of the dipole–target collision c.m. energy squared $\hat{s}$, and not Bjorken x. The corresponding parametrisation of the dipole cross section based upon the same saturated ansatz as in Eq. (4.8) is found by a replacement $\sigma _0 \rightarrow \overline{\sigma }_0(\hat{s})$ and $R_0 \rightarrow \overline{R}_0(\hat{s})$ where [15]

$$\begin{aligned}&\overline{R}_0(\hat{s})=0.88\,\mathrm {fm}\,(s_0/\hat{s})^{0.14},\\&\overline{\sigma }_0(\hat{s})=\sigma _\mathrm{tot}^{\pi p}(\hat{s}) \Big (1+\frac{3\overline{R}_0^2(\hat{s})}{8\langle r_\mathrm{ch}^2 \rangle _{\pi }}\Big ). \end{aligned}$$

Here, $\sigma _\mathrm{tot}^{\pi p}(\hat{s})=23.6(\hat{s}/s_0)^{0.08}$ mb is the pion–proton total cross section [73], $\langle r_\mathrm{ch}^2 \rangle _{\pi }= 0.44$ fm$^2$ is the mean pion radius squared [74], and $s_0=1000\,\mathrm{GeV}^2$. Interestingly enough, this parametrisation known as the KST model has been shown to give the correct description of the pion–proton cross section at scales up to $Q^2\sim 20\ \hbox {GeV}^2$. This parametrisation, together with the ones above, will be used in our analysis of exclusive real and virtual photoproduction of quarkonia.

Another parametrization was proposed by Iancu, Itakura and Munier [39]

$$\begin{aligned} \sigma _{q\bar{q}}(r,x)= & {} \sigma _0 N_0 \left( \frac{rQ_s(x)}{2}\right) ^{2\gamma _{eff}(r,x)} \quad rQ_s(x)\le 2\nonumber \\= & {} \sigma _0 \left( 1-e^{-A\ln ^{2}(BrQ_s(x))}\right) \quad \, rQ_s(x) > 2 \nonumber \\ \gamma _{eff}(r,x)= & {} \gamma +\frac{1}{\kappa \lambda Y}\ln \left( \frac{2}{rQ_s(x)}\right) \quad Y=\ln \left( \frac{1}{x}\right) ,\nonumber \\ \end{aligned}$$

(4.18)

where $\gamma _{eff}(r,x)$ is an effective anomalous dimension, $\kappa =9.9$, $N_0=0.7$, $Q_s^{2}(x)=\left( \frac{x_0}{x}\right) ^\lambda \ \hbox {GeV}^{2}$ and $\sigma _0=2\pi R_p^{2}$. Parameters A and B are chosen to ensure continuity between both parts of the parametrization at $rQ_s(x)= 2$ as

$$\begin{aligned} A= & {} -\frac{N_0^{2}\gamma ^{2}}{(1-N_0)^{2}\ln (1-N_0)}\nonumber \\ B= & {} \frac{1}{2}(1-N_0)^{-\frac{1-N_0}{N_0\gamma }}. \end{aligned}$$

(4.19)

Parameters $\sigma _0, R_p, \gamma , x_0$ and $\lambda$ have to be fitted to data. We use the fit performed in [41] with $\gamma =0.7376$, $\lambda =0.2197$, $x_0=0.1632\times 10^{-4}$ and $\sigma _0=27.33$ mb. This model will be denoted as IIM.

The last parametrization used in our analysis was proposed in Refs. [40, 75] as a modification of the IIM parametrization to include the explicit impact parameter dependence by introducing the modified saturation scale

$$\begin{aligned} Q_s^{2}(x)\rightarrow Q_s^{2}(x,b)=\left( \frac{x_0}{x}\right) ^\lambda \left( e^{-\frac{b^{2}}{2B_{G}}}\right) ^{\frac{1}{\gamma }}. \end{aligned}$$

(4.20)

Parameters $B_{G},N_0, \gamma , x_0$ and $\lambda$ has to be fitted to data. The most recent set of parameters comes from [44] and sets $\gamma =0.6492$, $N_0=0.3658$, $\lambda =0.2023$, $x_0=0.00069$ and $B_{G}=5.5 \ \hbox {GeV}^{-2}.$ This model is denoted as b-CGC.

In order to get the impact parameter independent dipole cross section from IPSat and b-CGC parametrizations an integral over the impact parameter was performed. As an illustration, in Fig. 5 we show the shape of different parametrisations for $\sigma _{q\bar{q}}(r,x)$ as a function of the dipole transverse separations r as two fixed values of the Bjorken variable $x=10^{-2}$ (left panels) and $x=10^{-4}$ (right panels).

Fig. 6

The dependence of the diffractive slope on c.m. energy B(W) given by Eq. (5.1) for the process $\gamma ^{*}\,p\rightarrow J/\psi \,p$ with two characteristic parameters $B_0$ and $\alpha '(0)$ determined by the fit to combined H1 [76, 77] and ZEUS [78, 79] data for two distinct $Q^2$ regions – low-$Q^2$ (photoproduction domain, left) and high-$Q^2$ (electroproduction domain, right)

At large dipole separations, we observe a substantial variation between both the shapes and magnitudes of $\sigma _{q\bar{q}}(r,x)$, and the differences in dipole models tend to rise with decreasing x. Quite interestingly, such differences become also large at very small dipole sizes $r\lesssim 0.05\div 0.06\,\,\text{ fm }$, i.e. in the perturbative region. Thus, the measurements of exclusive electroproduction of quarkonia at very large scales $Q^2 > rsim 300\div 400\,\,\text{ GeV }^2$ may provide additional constraints on the dipole parametrizations and means to further reduce theoretical uncertainties in the small-x treatment of the gluon density. Using the precision data in the hard and soft limits, one could ultimately start ruling out the models.

5 Numerical results vs data

5.1 Theoretical uncertainties caused by determination of the diffraction slope

Let us turn to discussion of numerical results on the $\gamma ^{*}\,p \rightarrow V\,p$ process in comparison with the data available from HERA. In order to calculate the total photo- and electroproduction cross section Eq. (2.13) with amplitudes given by Eqs. (C1) and (C2) one should know the magnitudes of the slope parameter as a function of the photon energy W and virtuality $Q^2$.

5.1.1 Diffraction slope for the process $\gamma ^{*}\,p\rightarrow J/\psi (\Upsilon )\,p$

For the c.m. energy behavior of the diffraction slope B(W) we use the standard form based on the Regge phenomenology,

$$\begin{aligned} B(W)=B_0 + 4\,\alpha '(0) \ln \Big (\frac{W}{W_0}\Big ) , \quad W_0 = 90 \, \,\text{ GeV }, \end{aligned}$$

(5.1)

where $\alpha '(0)$ represents the slope of the Pomeron trajectory.

Table 1 Parameters $B_0$ and $\alpha '(0)$ of the diffraction slope B obtained by a fit to different data sets at HERA in photo- ($Q^2 < 1\,\,\text{ GeV }^2$) and electroproduction ($Q^2 > 1\,\,\text{ GeV }^2$) of ground-state 1S-charmonium

Both parameters $B_0$ and $\alpha '(0)$ for the process $\gamma ^{*}\,p\rightarrow J/\psi \,p$ have been obtained by a fit to data from H1 [76, 77] and ZEUS [78, 79] collaborations at HERA as well as by our overall fit to the combined data from both collaborations as is shown in Fig. 6. Our fit resulted with $\chi ^2/ndf = 0.6$ for photoproduction and $\chi ^2/ndf = 3.75$ for electroproduction. The corresponding values are presented in Table 1.

Fig. 7

The exclusive $J/\psi$ electroproduction cross section as a function of c.m. energy W at fixed $Q^2=0.05\ \hbox {GeV}^2$ (left panel) and the scaling variable $Q^2+M_{J/\psi }^2$ at fixed $W=90\ \hbox {GeV}$ (right panel). The model calculations were performed with the $J/\psi$ wave function, obtained by using the BT potential [82], and with the phenomenological KST dipole cross section [15]. Here and below, the $Q^2$-dependent slope parameter labeled as “our fit” is determined from Eq. (5.2) with the energy behavior at $Q^2\rightarrow 0$ found in Table 1 and indicated there as “this work”. The results also incorporate the Melosh spin rotation effects. The data are provided by H1 [76, 77], ZEUS [78, 79], ALICE [83], E401 [84] and E516 [85] Collaborations

The values of $\alpha '(0)$ extracted from the available HERA data are in accordance with theoretical predictions in Ref. [80] based on the color dipole formalism and presented already in 1994. It was shown that the slope of the Pomeron trajectory is strongly correlated with the magnitude of the gluon correlation radius.

Since the data for the diffraction slope at $Q^2\gg 0$ are scarce, for the $Q^2$ dependence of the slope parameter $B(Q^2)$ we use the empirical parametrization from Ref. [81] based on the color dipole model calculations and valid for production of $J/\psi$ and $\Upsilon$ within the range of $Q^2\lesssim 500\,\,\text{ GeV }^2$,

$$\begin{aligned} B(W,Q^2)\approx B(W,Q^2=0) - B_1\,\ln \, \Big (\frac{Q^2 + M_V^2}{M_{J/\psi }^2}\Big ) , \end{aligned}$$

(5.2)

where the energy dependence of $B(W,Q^2=0)$ is determined using Eq. (5.1) with parameters found in Table 1, and $B_1 = 0.45\,\,\text{ GeV }^{-2}$. We tested that such a parametrization gives values of the slope parameter in a reasonable agreement with the existing data [76, 78] on electroproduction of $J/\psi$ at HERA.

Here we would like to emphasize that for the photo- and electroproduction of 1S bottomonium the corresponding diffraction slope can be estimated also from Eq. (5.2) as $B_{\Upsilon }(W,Q^2)\approx B_{J/\psi }(W,Q^2+M_{\Upsilon }^2)$.

5.1.2 Diffraction slope for the process $\gamma ^{*}\,p\rightarrow \psi '(\Upsilon ')\,p$

Detailed analysis of the diffraction slope in photo- and electroproduction of 2S-radially excited heavy quarkonia $\psi '(2S)$ and $\Upsilon '(2S)$ is presented in Ref. [81]. It was shown within the color dipole formalism that the inequality $B(2S) < B(1S)$ comes from the nodal structure of corresponding quarkonium wave functions. This is a direct consequence of the destructive interference of the contributions to the production amplitude coming from regions of small and large dipole separations. For production of bottomonia states, 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 '}(2S)\sim B_{\Upsilon }(1S)$.

However, for production of 2S-radially excited charmonium, the difference of diffraction slopes $\Delta _B=B(1S)-B(2S)$ can not be neglected. Model calculations within the color dipole formalism [81] at $W=90\,\,\text{ GeV }$ give the values $\Delta _B^T\sim 0.25\,\,\text{ GeV }^{-2}$ and $\Delta _B^L\sim 0.45\,\,\text{ GeV }^{-2}$ for photoproduction of T and L polarized $\psi '(2S)$, respectively, as a clear manifestation of the node effect. The quantity $\Delta _B$ gradually vanishes with $Q^2$ and can be neglected at $Q^2 > rsim 20\,\,\text{ GeV }^2$ as a result of a weak node effect at small dipole sizes. However, $\Delta _B$ rises towards small energies and at $W=15\,\,\text{ GeV }$ reaches much higher values, i.e. $\Delta _B^T\sim 0.38$ and $\Delta _B^L\sim 0.9\,\,\text{ GeV }^{-2}$ [81] for photoproduction of T and L polarized $\psi '(2S)$, respectively.

Fig. 8

The same as Fig. 7 but for the photo- and electroproduction of $\Upsilon (1S)$ state. The model predictions for $\sigma ^{\gamma ^{*}\,p\rightarrow \Upsilon \,p}(W,Q^2)$ are compared with the existing data from H1 [33], ZEUS [86, 87], CMS [88] and LHCb [89] experiments

In our calculations, we employ the following parametrization of the color dipole model predictions of the positive-valued part of $\Delta _B$ [81],

$$\begin{aligned} \Delta _B^{T,L} (W,Q^2)= c^{T,L}(W) \Bigg [1 - d(W)\, \ln \bigg (\frac{Q^2+M_{\psi '}^2}{M_{\psi '}^2}\bigg )\Bigg ]\ge 0 . \end{aligned}$$

(5.3)

Otherwise, $B(1S) = B(2S)$ for $\Delta _B^{T,L} (W,Q^2) \lesssim 0$ is adopted. Here, the energy-dependent coefficients are $c^T(W)=0.24-0.08\,\ln (W/W_0)\,\,\text{ GeV }^{-2}$ and $c^L(W)=0.45-0.24\,\ln (W/W_0)\,\,\text{ GeV }^{-2}$ for production of T and L polarized $\psi '(2S)$ state, respectively, and the factor $d(W) = 1.65 + 0.3\,\ln (W/W_0)$.

In what follows, in all figures we denote by “our fit” the model calculations that use the parametrization of the slope parameter given by Eq. (5.2), where the first term $B(W,Q^2=0)$ is determined from Table 1.

First, we test how uncertainties in determination of the diffraction slope for the process $\gamma ^{*}\,p\rightarrow J/\psi \,p$ lead to uncertainties in model predictions for the real and virtual photoproduction cross sections. For this purpose, we use the realistic BT potential [82] (see also Appendix A) in determination of the charmonium wave functions as well as the phenomenological KST dipole cross section [15], which provides a good description of hadronic processes, also at small scales corresponding to the nonperturbative region of large dipole sizes.

Figure 7 shows the color dipole model calculations versus the HERA data on the photo- and electroproduction cross sections as a function of the c.m. energy W at fixed $Q^2 = 0.05\,\,\text{ GeV }^2$ (left panel) and the scaling variable $Q^2+M_{J/\psi }^2$ at fixed $W=90\,\,\text{ GeV }$ (right panel) using the different parametrizations for the diffraction slope as depicted in Table 1. Corresponding amplitudes (C1) and (C2) for production of T and L polarized charmonia, entering the expression (2.13) for the electroproduction cross section, contain corrections for the Melosh spin rotation effects. Since the data on the $Q^2$ behavior of the diffraction slope are very scarce we took the results of model calculations [81], which can be simply parametrized by Eq. (5.2) and provide a reasonable description of the HERA data.

One can see from the left panel of Fig. 7 that model predictions using the constant value for the slope parameter $B=4.73\,\,\text{ GeV }^{-2}$ underestimate the data at lower c.m. energies $W\lesssim 100\,\,\text{ GeV }$. However, they lead to an overestimation of the ALICE experimental value [83] at higher $W\sim 700\,\,\text{ GeV }$. An agreement with the data can be improved by taking the energy-dependent diffraction slope with parameters from Table 1. All these parametrizations lead to very similar values for the diffraction slope at small energies but start to differentiate from each other at higher energies corresponding to the LHC energy range. Here, the best description of the data is achieved by the fit to only H1 data, as well as by our fit to the combined H1 and ZEUS data sets.

The right panel of Fig. 7 shows the model predictions for electroproduction cross section $\sigma ^{\gamma ^{*}\,p\rightarrow J/\psi \,p}(W,Q^2)$ as a function of the scaling variable $Q^2+M_{J/\psi }^2$ at fixed value of c.m. energy $W=90\,\,\text{ GeV }$. The $Q^2$ dependence of the slope parameter is given by the empirical formula Eq. (5.2), whereas for $B(W=90\,\,\text{ GeV },Q^2=0)$ we take different parametrizations from Table 1. As a result, the shape of the corresponding theoretical curves is almost identical describing the available data from H1 and ZEUS collaborations reasonably well.

Fig. 9

The same as Fig. 7 but for the real and virtual photoproduction of $\Upsilon '(2S)$ bottomonia

Fig. 10

The same as Fig. 7 but for the real and virtual photoproduction of $\psi '(2S)$ charmonium

Note that differences in model predictions using various parametrizations for the diffraction slopes can be treated as a measure of the underlined theoretical uncertainty.

Theoretical uncertainties in predictions of the real and virtual photoproduction cross sections of the elastic process $\gamma ^{*}\,p\rightarrow \Upsilon (1S)\,p$ inherent for determination of the slope parameter are depicted in Fig. 8. Here, similarly to the case of 1S-charmonium production, we compare the model predictions for different parametrizations of the slope parameter from Table 1. In the case of electroproduction of 1S bottomonium, the corresponding diffractive slope can be approximately estimated from Eq. (5.2) as follows: $B_{\Upsilon }(W,Q^2)\approx B_{J/\psi }(W,Q^2+M_{\Upsilon }^2)$. This is a consequence of the scaling properties in production of different vector mesons [81]. Here, we assume a similar value of the Pomeron trajectory slope $\alpha '(0)$ describing the energy dependence of the diffractive slope, see Eq. (5.1), for charmonium as well as for bottomonium production. This is supported by calculations of $\alpha '(0)$ performed in Ref. [80] within the color dipole formalism.

The left panel of Fig. 8 clearly demonstrates that inclusion of the energy dependent slope parameters brings our model predictions to a better agreement with the available data. As was already emphasized above, the differences in model predictions corresponding to different parametrizations of the diffractive slope can be considered as a good measure of the underlined theoretical uncertainty.

For the photo- and electroproduction of 2S radially-excited $\psi '(2S)$ and $\Upsilon '(2S)$ states the nodal structure of the corresponding wave functions (see Figs. 3, 4) causes an inequality $B(2S)\lesssim B(1S)$. The corresponding difference $B(1S)-B(2S)$ was calculated in Ref. [81] within the color dipole formalism and can be parametrized as is given by Eq. (5.3). For the photo- and electroproduction of $\Upsilon '(2S)$ the node effect can be neglected and we can safely take the same slope parameter as for the $\Upsilon (1S)$ state, namely, $B_{\Upsilon '}(2S)\sim B_{\Upsilon }(1S)$. The corresponding model predictions, taking four different parametrizations for the diffraction slope from Table 1, are presented in Fig. 9.

In comparison with $\Upsilon '(2S)$ eletroproduction, a stronger node effect in production of 2S radially-excited charmonium causes a larger difference of diffractive slopes given by Eq. (5.3) such that it can not be neglected. Consequently, one expects that $B_{J/\psi (1S)} > rsim B_{\psi '}(2S)$ [81]. The corresponding model predictions for $\sigma ^{\gamma ^{*}\,p\rightarrow \psi '(2S)\,p}(W,Q^2)$ including the different parametrizations for the slope parameter from Table 1, as well as the corrected diffraction slope in eletroproduction of 2S radially-excited charmonium, $B_{\psi '}(2S) = B_{J/\psi }(1S) - \Delta _B$, are shown in Fig. 10. One can see that the node effect leads to an enhancement of the $J/\psi (2S)$ photoproduction cross section, especially at small c.m. energies W, as well as at small values of $Q^2$ enabling a better agreement with the data (see also Fig. 11).

We would like to emphasize that one should distinguish between manifestations of the node effect in amplitude for production of 2S radially-excited quarkonia and in the magnitude of the corresponding diffraction slope. The nodal structure of the wave function for radially-excited states causes cancellations in the production amplitude from regions of large and small transverse sizes above and below the node position. Here, investigation of the ratio $R\equiv R(W,Q^2)$ of the $\psi '(2S)$-to-$J/\psi (1S)$ photo- and electroproduction cross sections allows to minimize the theoretical uncertainties connected to a determination of the corresponding slope parameters for $\gamma ^{*}\,p\rightarrow J/\psi (1S)\,p$ and $\gamma ^{*}\,p\rightarrow \psi '(2S)\,p$ processes.

Fig. 11

The color dipole model predictions for the $\psi '(2S)$-to-$J/\psi (1S)$ ratio of electroproduction cross sections as functions of c.m. energy W at fixed $Q^2=0.05\,\,\text{ GeV }^2$ (left panel), as well as functions of $Q^2$ at fixed $W=90\,\,\text{ GeV }$ (right panel) versus the existing data from H1 [90, 91], ZEUS [92] and fixed target experiments [93,94,95,96,97]. The solid and dashed lines correspond to calculations with and without the correction $\Delta _B$ given by Eq. (5.3) in determination of the slope parameter for the process $\gamma ^{*}\,p\rightarrow \psi '(2S)\,p$, respectively. The model calculations were performed with the charmonium wave functions obtained by using the BT potential [82] and with the phenomenological KST dipole cross section [15]. The Melosh spin rotation effects are included in this calculation

Fig. 12

The same as Fig. 7 but for different realistic $c-\bar{c}$ interaction potentials as described in Appendix A

Neglecting the impact of the node effect on the magnitude of the slope parameter $B_{\psi '}(2S)$, one can safely use the approximate equality $B_{J/\psi }(1S)\sim B_{\psi '}(2S)$ with a rather good accuracy. Consequently, $B_{J/\psi }(1S)$ and $B_{\psi '}(2S)$ cancel in the ratio $R(W,Q^2)$. Then the rise of R with c.m. energy W and $Q^2$ depicted by dashed lines in Fig. 11 is a characteristic manifestation of the node effect. Since the size of $\psi '(2S)$ is larger than $J/\psi (1S)$, one should naturally expect a stronger energy dependence for the $J/\psi (1S)$ electroproduction cross section because dipoles with a smaller transverse size have a steeper rise with energy. As a result, the ratio R(W) should decrease with energy. However, despite of this expectation, the nodal structure of the wave function for 2S radially-excited states causes an opposite effect, i.e. the rise of R(W) with energy. The steeper energy dependence at smaller dipole sizes below the node position diminishes the node effect at higher energies. This is a result of reduction of a cancellation in the 2S production amplitude from regions below and above the node position. This then leads to a steeper energy dependence of $\psi '(2S)$ compared to $J/\psi (1S)$ production cross section (compare Fig. 7 with Fig. 10). The rise of the ratio R(W) with c.m. energy W is depicted in the left panel of Fig. 11 where the model predictions are in accordance with the data, especially at higher energies $W > rsim 50\,\,\text{ GeV }$. Similarly, the node effect becomes weaker at larger $Q^2$ causing a rise of the ratio $R(Q^2)$ in a reasonable agreement with the existing data as is demonstrated in the right panel of Fig. 11.

Fig. 13

The same as Fig. 8 but for different realistic $b-\bar{b}$ interaction potentials as described in Appendix A

The node effect has some impact also on the magnitude of the diffractive slope for electroproduction of 2S radially-excited charmonium as was presented in Ref. [81] and discussed above. This leads to the following inequality $B(2S)\lesssim B(1S)$. The corresponding difference $\Delta _B$ was computed within the color dipole model in Ref. [81] and can be parametrized by Eq. (5.3). This correction $\Delta _B$ rises towards small W and $Q^2$ since the onset of the node effect becomes stronger and leads to an enhancement of the ratio $R(W,Q^2)$ as shown in Fig. 11 by solid lines. Notably, such an effect brings our predictions to a better agreement with the data at smaller energies $W\lesssim 20\,\,\text{ GeV }$.

5.2 Theoretical uncertainties caused by a shape of the $c-\bar{c}$ ($b-\bar{b}$) interaction potential

Here we analyze how determination of the quarkonium wave functions generated by various interquark interaction potentials leads to a different behavior of the photo- and electroproduction cross sections. The results for $J/\psi$ and $\Upsilon$ are shown in Figs. 12 and 13 in comparison with the available data. Our calculations were performed using the phenomenological KST parametrization for the dipole cross section and for the 1S quarkonium wave functions determined from the COR, HAR, LOG, POW and BT potentials described in Appendix A. Our observations are as follows:

• (i) The potentials labeled as HAR, BT and LOG well describe the photoproduction $J/\psi$ data, whereas the potential POW slightly overestimates the data while the potential COR significantly underestimates them by a factor of about $2\div 2.5$.

• (ii) Such a different behavior originates from different charm quark masses used in various potentials. The potentials BT and LOG use $m_c=1.5\,\,\text{ GeV }$, while HAR adopts $m_c=1.4\,\,\text{ GeV }$, POW – $m_c=1.3\,\,\text{ GeV }$ and COR takes $m_c=1.84\,\,\text{ GeV }$. Different potentials have only a small impact on the shape of wave functions for 1S-state charmonium (see Fig. 28). However, the photon wave function, Eq. (2.2) is extremely sensitive to the value of $m_c$ that enters the argument of the Bessel function $K_0$.

• (iii) The model predictions for the photoproduction cross section, quite expectedly, exhibit the following hierarchy: the smaller c-quark mass used in the realistic potential leads to higher values of the cross sections (see also Figs. 17, 18).

• (iv) Dependence of the $J/\psi$ electroproduction cross section (see the right panel of Fig. 12 on the scaling variable $Q^2+M_{J/\psi }^2$ follows from the structure of $\varepsilon ^2$ in Eq. (2.2). As was analyzed in Ref. [5] the nonrelativistic approximation with $z=0.5$ can be safely used for production of charmonia and, especially, bottomonia. In this approximation, $\varepsilon ^2$ takes the value $\propto Q^2+(2\,m_c)^2\approx Q^2+M_{J/\psi }^2$.

• (v) Similarly to photoproduction of 1S charmonium, the right panel of Fig. 12 shows a reasonable agreement of the data with our calculations using the COR, HAR, LOG and BT potentials. Differences in model predictions gradually decrease with $Q^2$ since the variation between the corresponding realistic potentials is weaker at smaller dipole transverse separations r (see Fig. 2). Only the HAR potential leads to much smaller values of the cross sections at large $Q^2$ due to a lack of the Coulomb-like behavior at small r.

• (vi) Model predictions for the $\Upsilon (1S)$ photoproduction cross section depicted in the left panel of Fig. 13 exhibit a rather good description of the data with the use of all five realistic potentials considered in this work. Thus, we confirm the universality property of the quarkonia production cross sections as functions of the scaling variable $Q^2+M_V^2$. Here, due to such universality the theoretical uncertainty given by a spread between the results obtained with different interquark potentials (see the left panel of Fig. 13) directly corresponds to the results for 1S charmonia electroproduction at $Q^2\sim M_{\Upsilon }^2$ (compare with the right panel of Fig. 12).

• (vii) A small variance of the model predictions made also for 1S bottomonia electroproduction using different realistic potentials is demonstrated in the right panel of Fig. 13. However, this variance rises with $Q^2$ due to growing differences between the considered $b-\bar{b}$ interaction potentials at small $\tilde{r}\lesssim 0.1\,\,\text{ fm }$ (see the right panel of Fig. 2).

5.3 Theoretical uncertainties caused by different parametrizations of the color dipole cross section

Fig. 14

The same as Fig. 7 but for different phenomenological dipole cross sections described in Sect. 4

Fig. 15

The same as Fig. 8 but for different phenomenological dipole cross sections described in Sect. 4

Fig. 16

The exclusive $J/\psi$ electroproduction cross section as a function of energy W at several fixed values of $Q^2=3.1\,\,\text{ GeV }^2$ (top left panel), $Q^2=6.8\,\,\text{ GeV }^2$ (top right panel), $Q^2=16.0\,\,\text{ GeV }^2$ (bottom left panel), and $Q^2=22.4\,\,\text{ GeV }^2$ (bottom right panel). The model predictions, including the Melosh spin rotation effects, were performed with the $J/\psi$ wave function using the BT potential [82] and for different phenomenological dipole cross sections described in Sect. 4. The data are taken from H1 [76] and ZEUS [79] Collaborations at HERA

The calculations performed in the framework of color dipole approach are strongly correlated with the shape of the dipole cross section, $\sigma _{q\bar{q}}(r,x)$. In our predictions using the BT realistic potential for determination of the quarkonium wave functions, in Figs. 14, 15 and 16 we test the eight main phenomenological parametrizations for $\sigma _{q\bar{q}}(r,x)$ found in the literature and discussed in Sect. 4. Here, the main observations are the following:

• (i) In the case of 1S charmonium photoproduction, the KST and GBWold dipole models give almost the same cross section at c.m. energies $W\lesssim 200\,\,\text{ GeV }$ describing the available data reasonably well. The other phenomenological parametrizations denoted as GBWnew, rcBK, coIBK, IM, bCGC and IPSat strongly underestimate the data by a factor of $2\div 3$ (see the left panel of Fig. 14.

• (ii) In the electroproduction of 1S charmonia, the KST and GBWold dipole cross sections lead to the cross sections that differ from each other by a factor of $2\div 3$ at high $Q^2$ (see the right panel of Fig. 14). Here, the KST parametrization provides the best description of $Q^2$-dependent data. Other six parametrizations used in our study grossly underestimate the data within the whole considered $Q^2$ interval.

• (iii) A similar conclusion as above can be made also from Fig. 16 where we studied the energy dependence of 1S charmonium electroproduction cross section at different fixed values of $Q^2$.

• (iv) In analogy to electroproduction of 1S charmonia, the model calculations using the KST phenomenological parametrization for the dipole cross section provide the best description of the available data on photoproduction of 1S bottomonia as shown in Fig. 15. Except for the rcBK parametrization at large W, all other parametrizations lead to a significant underestimation of these data in the whole range of W.

• (v) A huge variance of the model predictions for the photoproduction ($Q^2\rightarrow 0$) cross sections of $\Upsilon (1S)$ using various parametrizations for the dipole cross section remains also in the case of electroproduction results shown in the right panel of Fig. 15. Here, the spread between the results rises with $Q^2$ as a direct consequence of the growing differences between the dipole parametrizations at small transverse separations $r\lesssim 0.1\,\,\text{ fm }$. The latter is demonstrated by bottom panels of Fig. 5.

5.4 Sensitivity of model predictions to quark mass

The quark mass has a strong impact on magnitudes of the model predictions as was presented and discussed earlier in Sect. 5.2. Different realistic potentials (see Appendix A) used in our analysis of the quarkonium wave functions contain distinct values of quark masses ranging within the interval $m_c\in (1.3 - 1.84)\,\,\text{ GeV }$, for the charm quark, and $m_b\in (4.2 - 5.17)\,\,\text{ GeV }$, for the bottom quark. Here we test, taking the BT potential as a reference point with $m_c=1.48\,\,\text{ GeV }$ and $m_b=4.87\,\,\text{ GeV }$, how much our model predictions are modified by changing the quark mass from the minimal to maximal values corresponding to these intervals.

Fig. 17

The same as Fig. 7 but for the test of sensitivity of the model predictions to the typical charm quark mass $m_c$ variations

Fig. 18

The same as the left panel of Fig. 8 but for the test of sensitivity of the model predictions to the bottom quark mass $m_b$

Figure 17 clearly demonstrates the sensitivity of model results, taking the realistic BT potential and KST parametrization of the dipole cross section, to different quark mass values. Whereas our calculations, using the BT potential with $m_c=1.48\,\,\text{ GeV }$, lead to a reasonable description of the data, a modification of the charm quark mass to the lower ($m_c=1.3\,\,\text{ GeV }$) and higher ($m_c=1.84\,\,\text{ GeV }$) value causes a gross overestimation and underestimation of these data, respectively. Such a strong sensitivity to the value of the charm quark mass comes from the photon wave function, Eq. (2.2), which contains $m_c$ in the argument of the Bessel function $K_0$.

The sensitivity of model predictions to quark mass values gradually decreases with $Q^2$ since, in comparison to photoproduction limit $Q^2\rightarrow 0$, the quark mass scale plays a weaker role and can be neglected at large $Q^2\gg m_c^2$. Then, the model calculations naturally give very similar values for the 1S charmonium electroproduction cross section as is demonstrated in the right panel of Fig. 17.

Fig. 19

The same as Fig. 7 but showing the effect of the Melosh spin rotation in the exclusive $J/\psi$ electroproduction cross section shown as a function of c.m. energy W (left panel) and the scaling variable $Q^2+M_{J/\psi }^2$ (right panel)

Fig. 20

The same as Figs. 8 and 19 but for exclusive electroproduction of $\Upsilon (1S)$ bottomonium

A variation of the model predictions with quark mass is presented in Fig. 18 for the case of photo- and electroproduction of $\Upsilon (1S)$. In comparison to charmonium production, here the sensitivity of the cross section to bottom quark mass variations is weaker due to a smaller relative change in $m_b$ and also gradually decreases with $Q^2$ at large $Q^2\gg m_b^2$ as expected.

5.5 Spin rotation effects in electroproduction of 1S quarkonia

The effects of the Melosh spin rotation [see Eqs. (C1) and (C2) in Appendix C] in diffractive electroproduction of S-wave heavy quarkonia have been studied in detail in the framework of color dipole formalism in Ref. [9]. For this reason, we present here only the main features of the spin rotation and demonstrate how much the spin effects can modify the corresponding photo- and electroproduction cross sections.

The role of spin effects in photo- and electroproduction of 1S charmonium is presented in Fig. 19. It leads to an enhancement of the photoproduction cross section by $\approx 20\div 30\%$ leading a better agreement with the data (see the left panel of Fig. 19). This fact clearly supports an importance of the Melosh spin transformation, which is obviously neglected in many present studies of diffractive photo- and electroproduction of heavy quarkonia.

The right panel of Fig. 19 demonstrates that the onset of spin effects gradually diminishes with the scaling variable $Q^2+M_{J/\psi }^2$ and leads to a better description of the data, especially at small and medium $Q^2\lesssim 20\div 30\,\,\text{ GeV }^2$.

As was analyzed recently in Ref. [9], the universal properties in production of different vector mesons cause a similar onset of spin rotation effects in production of charmonia and bottomonia at the same fixed values of the scaling variable $Q^2+M_V^2$. For this reason, we predict a weak onset of these effects also in the photoproduction of $\Upsilon (1S)$ state corresponding to electroproduction of $J/\psi (1S)$ at $Q^2\sim M_{\Upsilon }^2$ (compare the right panel of Fig. 19 with the left panel of Fig. 20). The weak onset of the Melosh spin transformation in $\Upsilon (1S)$ photoproduction decreases further with $Q^2$ as is demonstrated in the right panel of Fig. 20.

5.6 Theoretical uncertainties in predictions for the $\psi '(2S)$-to-$J/\psi (1S)$ and $\Upsilon '(2S)$-to-$\Upsilon (1S)$ ratios

The theoretical uncertainties presented above in Sects. 5.2, 5.3, 5.4 and 5.5 can be tested by investigating also the ratios $R_{2S/1S}$ for charmonia $\psi '(2S)$-to-$J/\psi (1S)$ and bottomonia $\Upsilon '(2S)$-to-$\Upsilon (1S)$ photo- and electroproduction cross sections. Such a study enables us to minimize the uncertainties providing with more stable and accurate predictions, which can be verified by the future measurements.

Fig. 21

The color dipole model predictions for the $\psi '(2S)$-to-$J/\psi (1S)$ (top panels) and $\Upsilon '(2S)$-to-$\Upsilon (1S)$ (bottom panels) ratios of electroproduction cross sections as functions of c.m. energy W at fixed $Q^2=0.05\,\,\text{ GeV }^2$ (left panels) and $Q^2$ at fixed $W=90\,\,\text{ GeV }$ (right panels) versus the data from H1 [90, 91], ZEUS [92] and fixed target experiments [93,94,95,96,97]. The calculations were performed for the quarkonium wave functions generated by different realistic $c-\bar{c}$ and $b-\bar{b}$ interaction potentials as depicted in Appendix A and with the phenomenological KST dipole cross section [15]. The results include the Melosh spin rotation

5.6.1 Dependence on $c-\bar{c}$ and $b-\bar{b}$ interaction potentials

In Fig. 21, using the KST phenomenological parametrization (see Sect. 4) for the dipole cross section, we test the sensitivity of model predictions for the $\psi '(2S)$-to-$J/\psi (1S)$ and $\Upsilon '(2S)$-to-$\Upsilon (1S)$ ratios with respect to the choice of interaction potentials which are employed in deriving the corresponding quarkonium LC wave functions.

One can notice in top panels of Fig. 21 a good agreement of our calculations with the experimental data for all realistic potentials (COR, LOW, POW and BT), except for the HAR potential, which grossly overestimates the data at large W and $Q^2$. It is caused by the lack of Coulomb-like behavior in the HAR potential, which amplifies the role of the node effect. This is based upon a stronger enhancement of the small-r domain of the $\psi '(2S)$ and $\Upsilon '(2S)$ wave functions below the node position, therefore leading to a stronger reduction of the cancellation between low-r and high-r domains in the 2S production amplitude. Since the role of the Coulomb-like behavior increases in production of bottomonia, the bottom panels of Fig. 21 clearly demonstrate a huge difference in predictions between the HAR potential and all the other potentials. The latter generate only a small variance in the corresponding results for the $\Upsilon '(2S)$-to-$\Upsilon (1S)$ ratio.

5.6.2 Dependence on the phenomenological dipole cross sections

In Sect. 5.3 we studied a correlation of the model predictions for photo- and electroproduction of 1S quarkonium with a shape of the color dipole cross section, $\sigma _{q\bar{q}}(r,x)$. We found a huge variance in the model predictions for the electroproduction cross section by using eight different popular parametrizations for $\sigma _{q\bar{q}}(r,x)$ discussed in Sect. 4. Here, we test how large is the theoretical uncertainty in the model predictions for the ratio $R_{2S/1S}(W,Q^2)$ caused by such a variety of different treatments of the target gluon density encoded in these parametrizations.

The results of our calculations are depicted in Fig. 22. One can see that, in comparison to the electroproduction cross section, the study of $R_{2S/1S}$ ratio (utilizing, for example, the realistic BT potential) allows to reduce substantially the uncertainty of our predictions stemming from different existing parametrizations for $\sigma _{q\bar{q}}(r,x)$ (compare Fig. 22 with the results of Sect. 5.3).

On the other hand, such a study makes it possible to analyze how the node effect manifests itself for different shapes of the color dipole cross section. The onset of the node effect is controlled by an increase of the ratio $R_{2S/1S}$ with energy W and photon virtuality $Q^2$. The stronger is the cancellation in 2S production amplitude, the steeper is the rise of $R_{2S/1S}(W,Q^2)$ with a rate, which is slightly different for various dipole parametrizations. For production of 2S bottomonia the node effect is much weaker as one can see in the bottom panels of Fig. 22.

Note, that the rise of such variations in the model predictions towards small energies can be influenced by a worse accuracy in dipole phenomenological parametrizations at the corresponding (large) values of Bjorken $x > rsim 0.1$. This is due to a natural limitation of the color dipole approach that is expected to fail at sufficiently large Bjorken x.

Fig. 22

The same as Fig. 21 but for the quakonium wave functions generated by the realistic BT potential and for different dipole cross section parametrizations described in Sect. 4

Fig. 23

The same as Fig. 21 but for the test of sensitivity of the dipole model model predictions to the charm $m_c$ and bottom $m_b$ quark mass variations

5.6.3 Dependence on the mass of charm and bottom quark

The study of $R_{2S/1S}(W,Q^2)$ ratios in production of quarkonia also allows to minimize the underlined theoretical uncertainties in our knowledge of the corresponding quark mass value. In Fig. 23, we test a variance in the model predictions taking values of $m_c$ and $m_b$ determined from the BT potential used in the calculations as well as the minimal and maximal $m_c$ and $m_b$ values occurring along all the other realistic potentials studied in this work as was described above in Sect. 5.4. One can see that the sensitivity of $R_{2S/1S}$ to different values of $m_c$ and $m_b$ is much weaker in comparison to the results for the photo- and electroproduction cross sections (compare with Figs. 17 and 18).

5.6.4 Importance of spin effects

In comparison to production of 1S quarkonia (see Sect. 5.5), as a consequence of the node effect leading to a cancellation in the production amplitude from regions below and above the node position, the onset of spin rotation effects is much stronger in proto- and electroproduction of radially-excited $\psi '(2S)$, $\Upsilon '(2S)$ and $\Upsilon ''(3S)$ as was recently discussed in detail in Ref. [9]. Here, we predict a dramatic effect of the Melosh spin transformation in charmonium electroproduction causing an increase of the $R_{2S/1S}(W,Q^2)$ ratio by a factor of $2\div 3$ as is demonstrated in the top panels of Fig. 24.

Fig. 24

The same as Fig. 21 but for demonstration of the importance of the Melosh spin rotation effects

Fig. 25

The ratio of integrated cross sections for elastic electroproduction of longitudinally (L) and transversely (T) polarized $J/\psi (1S)$ (left panel) and $\Upsilon (1S)$ (right panel) as a function of the scaling variable $Q^2+M_{J/\psi }^2$ (left panel) and $Q^2+M_{\Upsilon }^2$ (right panel) at fixed c.m. energy $W=90\,\,\text{ GeV }$. The data are taken from H1 [76] and ZEUS [79] collaborations. The results, including also the Melosh spin rotation effects, were obtained using the phenomenological KST dipole cross section [15] and five different $c-\bar{c}$ and $b-\bar{b}$ interaction potentials described in Appendix A

One can see that such a substantial enhancement of $R_{2S/1S}$ due to the spin effects brings our predictions, using the KST dipole parametrization and the realistic BT potential, to the values close to the experimental data. Here, the rise of $R_{2S/1S}(W,Q^2)$ with c.m. energy W and with $Q^2$ is yet another manifestation of the node effect as was discussed in Ref. [9].

Due to a weaker node effect at larger $Q^2$, we predict that the spin rotation effects gradually diminish with $Q^2$ as is demonstrated in the right panels of Fig. 24 for charmonium and bottomonium ratios $R_{2S/1S}$. Since the same values of the scaling variable $Q^2+M_V^2$ lead to a similar onset of various effects in production of different quarkonia, we predict a weak onset of spin effects also in photo- and electroproduction of bottomonia (see also Ref. [9]) at the corresponding photon virtuality $Q^2(\Upsilon )\approx Q^2(J/\psi ) + M_{J/\psi }^2$ as is shown in the bottom panels of Fig. 24.

5.7 Theoretical uncertainties in predictions for the ratio $\sigma _L^{\gamma ^{*}\,p\rightarrow J/\psi (\Upsilon )\,p}/\sigma _T^{\gamma ^{*}\,p \rightarrow J/\psi (\Upsilon )\,p}$

The theoretical uncertainties in predictions, presented above in Sects. 5.2, 5.3, 5.4 and 5.5, can be eliminated to a large extent by investigating the ratio of the elastic electroproduction cross sections of longitudinally and transversely polarized quarkonia. In Fig. 25 we present our results for such ratios $R_{LT}^{J/\psi } =\sigma _L^{\gamma ^{*} \,p\rightarrow J/\psi \,p} / \sigma _T^{\gamma ^{*}\,p\rightarrow J/\psi \,p}$ and $R_{LT}^{\Upsilon } = \sigma _L^{\gamma ^{*}\,p\rightarrow \Upsilon \,p} / \sigma _T^{\gamma ^{*}\,p\rightarrow \Upsilon \,p}$ as functions of the scaling variables $Q^2+M_{J/\psi }^2$ and $Q^2+M_{\Upsilon }^2$, respectively. One can see a rather good agreement of $R_{LT}^{J/\psi }$ with the existing data for all considered $c-\bar{c}$ potentials. Our predictions for the ratio $R_{LT}^{\Upsilon }(Q^2)$ can be tested by future measurements.

Here we would like to emphasize that the variation in model predictions for the ratio $R_{LT}$ using different quarkonium wave functions generated by distinct potentials is much less pronounced than that observed in Sect. 5.2 for the standard photo- and electroproduction cross sections (see Figs. 12, 13).

Fig. 26

The same as Fig. 7 but for illustration of the onset of the skewness effect as a function of c.m. energy W and the scaling variable $Q^2+M_{J/\psi }^2$. The model calculations were performed with the $J/\psi$ wave function generated by the BT potential [82] and with the phenomenological GBWnew dipole cross section [40] including the Melosh spin rotation

5.8 The skewness effect in electroproduction of quarkonia

The skewness correction is frequently interpreted in the literature as an effect when the gluons attached to the $Q\bar{Q}$ fluctuation of the photon carry (very) different light-front fractions x and $x'$ of the proton momentum [98, 99]. The corresponding expression for the correction factor $R_g$ has the following form [98],

$$\begin{aligned} R_g(\lambda _{T,L})=\frac{2^{2\lambda _{T,L} +3}}{\sqrt{\pi }} \frac{\Gamma (\lambda _{T,L} + 5/2)}{\Gamma (\lambda _{T,L} + 4)} , \end{aligned}$$

(5.4)

with $\lambda _{T,L}$ determined by Eq. (2.14) for both the photon polarizations T and L. Consequently, the skewness correction is then accounted for by multiplying the cross section $\sigma _{T,L}^{\gamma ^{*}\,p\rightarrow V\,p}(x,Q^2)$, Eq. (2.13), by a factor of $R_g^2(\lambda _{T,L})$.

However, the shape of the correction factor in Eq. (5.4) has been derived in Ref. [98] within the next-to-leading order approximation assuming the strong inequalities $x'\ll x\ll 1$ in the small-t region and for the specific power-law form of the diagonal gluon density of the target. Here, such a small-x shape of the gluon density is not fully probed within kinematic regions studied in the present paper, and consequently, may not be fully consistent with those extracted from different dipole models used in our calculations.

The statement from Ref. [44] that the skewness correction given by Eq. (5.4) can be incorporated into the bCGC dipole model is not fully consistent for the case of electroproduction of heavy quarkonia. The dipole amplitude is related to the gluon structure function of the target only at sufficiently large $Q^2\sim \Lambda /r^2$ [see also Eq. (4.2)], where the dipole sizes $r\lesssim r_0\sim 0.3\,\,\text{ fm }$ ($r_0$ is the gluon propagation radius [100, 101]) and the large numerical factor $\Lambda \approx 10$ have been estimated in Ref. [63]. In the case of quarkonium production, this condition requires rather large values of the saturation scale squared corresponding to the bCGC dipole model, $Q_s^2(x)$. This leads to rather small values of the Bjorken variable $x\lesssim 10^{-5}\div 10^{-6}$ necessary for justification of Eq. (5.4) for the skewness correction. Such small x-values correspond to way too large c.m. energies $W > rsim 10^3\,\,\text{ GeV }$, which are far beyond the energy range studied in the present paper. The same conclusion concerns also the other dipole models since the corresponding saturation scales are similar to that in the bCGC parametrization.

Fig. 27

The same as Fig. 26 but for exclusive electroproduction of 1S bottomonia. The data are taken from from H1 [33], ZEUS [86, 87], CMS [88] and LHCb [89] Collaborations

Since the exact analytical expression for $R_g$ is not available in the literature, we present here only a phenomenological estimation of the onset of the skewness effect in electroproduction of heavy quarkonia relying on the known approximate relation, Eq. (5.4). The results are depicted in Figs. 26 and 27 for the case of electroproduction of charmonia and bottomonia in the ground state, respectively.

The model calculations have been performed, as an example, with the phenomenological GBWnew dipole cross section [40] and with the quarkonium wave functions generated by the realistic BT potential [82]. One can see from Figs. 26 and 27 that the skewness correction increases the photo- and electroproduction cross section of quarkonia by a factor of $\sim 1.5\div 1.6$. As was analyzed in Sect. 5.3, neglecting the skewness correction, only the KST and GBWold dipole parametrizations lead to the best description of the available data on quarkonium electroprodution, whereas other phenomenological dipole cross sections grossly underestimate these data. Consequently, one can expect that the onset of the factor $R_g$ in our calculations should cause a slight overestimation of data for the KST but would lead to an improvement of the data description using not only GBWnew but also other dipole parametrizations. Namely, such an effort to obtain a better agreement with the data typically generates the main reason to include formally the skewness effects adopting only an approximate relation (5.4) based on assumptions, which can not be naturally adopted or justified for an arbitrary process. This is the basic motivation for us not to include the skewness factor in the rest of the calculations in the previous sections, instead, showing the more justified color dipole model predictions and estimates for the underlined theoretical uncertainties.

6 Conclusions

We have presented an exploratory and comprehensive study of elastic photo- and electroproduction of heavy quarkonia within the color dipole formalism. The main motivation is based on a growing interest in this topic, mainly in connection with an extensive ongoing investigation of quarkonium production processes in ultra-peripheral collisions at the RHIC and LHC facilities. Although the color dipole approach is well-known already of about thirty years and a wealth of research has been done, it is frequently used in the literature without a deeper understanding of the underlined theoretical uncertainties in predictions caused by various effects and properties of particular ingredients entering into the production amplitudes. Consequently, in order to obtain a better agreement with the data, this leads to an ongoing effort to include some additional new phenomena or additional ingredients instead of a better understanding the corresponding uncertainties or performing more accurate calculations. For this reason, in this paper we try to describe and analyze various sources of theoretical uncertainties and study their impact on the magnitude of the corresponding electroproduction cross sections for a large variety of quarkonia states and physics inputs.

In the color dipole formalism the production amplitude, given by the factorized light-cone expression (2.1), has the following ingredients: (i) the perturbative light-cone wave functions for the heavy $Q\bar{Q}$ fluctuation of the photon, (ii) the light-cone wave functions for the S-wave quarkonia states, and (iii) a phenomenological dipole cross section $\sigma _{q\bar{q}}(r,x)$ describing the interaction of the $Q\bar{Q}$ fluctuation with a proton target.

A description of the photon wave function is well-known and quite well understood, so it should not cause major uncertainties in calculations of the production amplitude. On the other hand, the determination of the quarkonium light-cone wave functions remains rather uncertain. Here, we adopted the frequently used prescription from Ref. [36] for the transition from the $Q\bar{Q}$ rest frame to the infinite momentum one. The corresponding quarkonium wave functions in the $Q\bar{Q}$ rest frame have been obtained by solving the Schrödinger equation for various $Q-\bar{Q}$ interaction potentials. Such an ambiguity in determination of quarkonium wave functions represents one of more relevant sources of theoretical uncertainties.

The essential ingredient in our calculations of the photo- and electroproduction cross sections of heavy quarkonia is the dipole cross section $\sigma _{q\bar{q}}(r,x)$. Here, we adopted the total of eight main phenomenological parametrizations for $\sigma _{q\bar{q}}(r,x)$ found in the literature that exhibit a saturated form at large transverse separations (dipole sizes) r as well as roughly satisfy the characteristic small-r behavior, $\sigma _{q\bar{q}}(r,x)\propto r^2$ for $r\rightarrow 0$ (color transparency). The differences in the corresponding parametrizations for $\sigma _{q\bar{q}}(r,x)$ represent another source of theoretical uncertainties in calculations of dipole amplitudes and, subsequently, of the corresponding electroproduction cross sections.

In order to avoid a double counting, the effect of higher Fock states, $Q\bar{Q}G$, $Q\bar{Q}GG$, ..., containing gluons in the photon wave function can be reabsorbed into the energy (Bjorken x) dependence of $\sigma _{q\bar{q}}(r,x)$. On the other hand, the dipole cross section has a steeper rise with energy at smaller dipole sizes due to more intensive gluon radiation. Here all cross sections at different dipole sizes are expected to follow the universal asymptotic properties at very large energies controlled by the Froissart bound.

The model predictions for the exclusive quarkonium electroproduction cross sections depend on the magnitude of the diffraction slope B [see Eq. (2.13)] for the corresponding elastic process $\gamma ^{*}\,p\rightarrow V\,p$, where the vector meson $V = J/\psi (1S), \psi '(2S), \Upsilon (1S), \Upsilon '(2S), \Upsilon ''(3S)$, etc. The energy dependence of B(W) has been obtained by the fit to the available data at HERA [see Eq. 5.1 and Table 1]. Since the data on the $Q^2$ behavior of the slope parameter are very scarce, we adopted a phenomenological model from Ref. [81] leading to an empirical parametrization (5.2), which gives the values of $B(Q^2)$ in a reasonable agreement with the data. Within the same model, we have included also the differences in slope parameters $B(1S)-B(2S)$, corresponding to production of the 1S-ground state and 2S-radially excited quarkonia. These differences come as a direct manifestation of the node effect in the quarkonium wave functions, in particular, leading to a cancellation in the production amplitude coming from regions in r below and above the node position. We have verified that different parametrizations of the energy evolution, with modelled $Q^2$ behavior of the slope parameter, cause only rather small uncertainties in the model predictions using various combinations of the quarkonium wave functions and phenomenological dipole parametrizations for $\sigma _{q\bar{q}}(r,x)$.

Another source of uncertainties studied in this work refers to the effect of the Melosh spin rotation, which is often neglected in the literature. We found that such spin effects are very important, especially in elastic photoproduction of quarkonia. They lead to a $\approx 20\div 30\%$ rise of the $J/\psi (1S)$ photoproduction cross section contributing to a better agreement of the model predictions with the data. However, they cause even more dramatic effect in $\psi '(2S)$ photoproduction substantially increasing the corresponding cross sections, as well as the $\psi '(2S)$-to-$J/\psi (1S)$ ratio, by a factor of $2\div 3$ (see also Ref. [9]).

We have also presented and discussed a large sensitivity of the model predictions to the value of heavy quark mass $m_Q$ which is caused by the photon wave function, Eq. (2.2) containing $m_Q$ in the argument of the Bessel function $K_0$.

Although the skewness correction is frequently used in calculations of the quarkonium photo- and electroproduction cross sections, only an approximate relation, Eq. (5.4), is known for the corresponding correction factor $R_g$. Since the exact analytical formula for $R_g$ is not available in the literature, we estimated a magnitude of this effect relying on the known expression (5.4) and found that the skewness correction increases the quarkonium electroproduction cross section by a factor of $\sim 1.5\div 1.6$. However, it is questionable to what extent and with what accuracy the approximate relation, Eq. (5.4), can be applied to quarkonium electroproduction within the kinematic ranges studied in the present paper.

Finally, we have found that all these sources of theoretical uncertainties can be reduced to a large extent when investigating the ratios of the cross sections such as $R_{2S/1S}(W,Q^2) = \sigma ^{\gamma ^{*}\,p\rightarrow \psi '(2S) (\Upsilon '(2S))\,p}(W,Q^2) / \sigma ^{\gamma ^{*}\,p\rightarrow J/\psi (1S) (\Upsilon (1S))\,p}(W,Q^2)$, as well as $R_{L/T}(W,Q^2) = \sigma _L^{\gamma ^{*}\,p\rightarrow J/\psi (\Upsilon )\,p}(W,Q^2) / \sigma _T^{\gamma ^{*}\,p\rightarrow J/\psi (\Upsilon )\,p}(W,Q^2)$. We have demonstrated that, in comparison to the standard quarkonium electroproduction cross sections, the ratios $R_{2S/1S}$ and $R_{L/T}$ exhibit much smaller variations generated by these uncertainties and thus produce more stable and accurate results, which can be tested by the future experiments.

To summarize, in our current analysis performed within the color dipole formalism we have used for the first time a combination of several new ingredients simultaneously, such as the proper light-cone wave functions of heavy quarkonia generated by realistic interquark interaction potentials, together with the Melosh spin rotation and the most recent models for the saturated dipole cross section. We have successfully described the existing $J/\psi$, $\psi '$ and $\Upsilon$ photo- and electroproduction data off the nucleon target. This encourages us to extend consequently such an analysis, going beyond the NRQCD approximation, also for nuclear targets and verify our predictions for vector meson photoproduction by comparing with the recent data obtained from ultra-peripheral heavy-ion collisions at RHIC and LHC. The corresponding new predictions can be tested then by the future (e.g. LHeC) measurements.

Finally, we would like to emphasize that the most of the results presented in the current paper can also be obtained interactively on our webpage https://hep.fjfi.cvut.cz/vm.php, where the model predictions for the photo- and electroproduction cross sections can be readily computed for various combinations of the quarkonium wave functions with particular dipole parametrizations for $\sigma _{q\bar{q}}(r,x)$ including or neglecting the Melosh spin rotation effects. Such an online tool is expected to be very useful for QCD practitioners and experimentalists working in the research areas connected to quarkonia physics.

Appendices

Appendix A: Quarkonia potentials

In order to compute the quarkonium wave function, one needs to specify an interaction potential between heavy quarks. Here, we provide the details of several distinct models for interquark potentials used in our numerical analysis.

1.1 1. Harmonic oscillator

The potential for harmonic oscillator (denoted as HAR)

$$\begin{aligned} V({\tilde{r}}) =\frac{1}{2}\,m_Q\,\omega ^{2}\,{\tilde{r}}^{2} , \quad \omega =\frac{1}{2}(M_{2S}-M_{1S}) , \end{aligned}$$

(A1)

is the simplest and the most common choice that leads to the Gaussian shape of the wave function. The masses of charm c and bottom b quarks are taken to be $m_c=1.4\ \hbox {GeV}$ and $m_b=4.2\ \hbox {GeV}$, respectively. The parameter $\omega$ is fixed to $0.3\ \hbox {GeV}$, for charmonia, and to $0.28\ \hbox {GeV}$, for bottomonia. The Schrödinger equation with this potential has an analytic solution

$$\begin{aligned} u({\tilde{r}})=\exp \biggl [ -\frac{1}{4}\,m_Q\,\omega \, {\tilde{r}}^2\biggr ] , \end{aligned}$$

(A2)

however, we obtain a solution of the Schrödinger equation for the harmonic oscillator numerically.

1.2 2. Cornell potential

The Cornell potential (COR) given by

$$\begin{aligned} V({\tilde{r}})=- \frac{k}{{\tilde{r}}}+\frac{{\tilde{r}}}{a^{2}}\quad k=0.52 \quad a=2.34\,\mathrm {GeV}^{-1} , \end{aligned}$$

(A3)

with $m_c=1.84\ \hbox {GeV}$ and $m_b=5.17\ \hbox {GeV}$, was initially proposed in Refs. [102, 103] and was also used in quarkonia photoproduction studies in Refs. [8, 42].

1.3 3. Logarithmic potential

The logarithmic potential (LOG) given by

$$\begin{aligned} V({\tilde{r}})=- 0.6635\,\,\text{ GeV }+ (0.733\,\,\text{ GeV })\, \log \bigl ({\tilde{r}}\cdot 1\,\,\text{ GeV }\bigr ) ,\nonumber \\ \end{aligned}$$

(A4)

with $m_c=1.5\,\,\text{ GeV }$ and $m_b=5.0\,\,\text{ GeV }$, is motivated by Ref. [104] and was also used in quarkonia photoproduction studies in Ref. [8].

1.4 4. Power-law potential

The effective power-law potential (POW) is given by

$$\begin{aligned} V({\tilde{r}})=- 6.41\,\,\text{ GeV }+ (6.08\,\,\text{ GeV })\,({\tilde{r}}\cdot 1\,\,\text{ GeV })^{0.106} , \end{aligned}$$

(A5)

with $m_c=1.334\,\,\text{ GeV }$ and $m_b=4.721\,\,\text{ GeV }$, is motivated by Ref. [105, 106] and the values were taken from Ref. [107].

1.5 5. Buchmüller-Tye potential

The Buchmüller-Tye potential (BT) [82] has a Coulomb-like behaviour at small ${\tilde{r}}$ and a string-like behaviour at large ${\tilde{r}}$. Its structure is similar to the Cornell potential but with additional corrections, particularly effective at small ${\tilde{r}}$. Namely,

$$\begin{aligned} V({\tilde{r}})=\frac{k}{{\tilde{r}}}-\frac{8\pi }{27} \frac{v(\lambda {\tilde{r}})}{{\tilde{r}}} , \end{aligned}$$

(A6)

for ${\tilde{r}}\ge 0.01\,\,\text{ fm }$, and

$$\begin{aligned} V({\tilde{r}}) =\,&-\frac{16\pi }{25}\frac{1}{{\tilde{r}} \,\ln \bigl (w({\tilde{r}})\bigr )}\nonumber \\&\times \!\left( 1 + 2\left( \gamma _E+\frac{53}{75}\right) \frac{1}{\ln \bigl (w({\tilde{r}})\bigr )} -\frac{462}{625}\, \frac{\ln \Big ( \ln \bigl ( w({\tilde{r}}) \bigr )\Big )}{\ln \bigl ( w({\tilde{r}})\bigr )}\right) , \end{aligned}$$

(A7)

for ${\tilde{r}} < 0.01\,\,\text{ fm }$. Here,

$$\begin{aligned}&w({\tilde{r}})=\frac{1}{\lambda ^2_\mathrm{MS}\,{\tilde{r}}^2}, \quad \lambda _\mathrm{MS}=0.509\,\,\text{ GeV },\nonumber \\&\quad k=0.153\,\,\text{ GeV }^2,\quad \lambda =0.406\,\,\text{ GeV }, \end{aligned}$$

(A8)

$\gamma _E=0.5772$ is the Euler constant, and the function v(x) is provided numerically in Ref. [82]. This potential uses the following quark mass values: $m_c=1.48\,\,\text{ GeV }$ and $m_b=4.87\,\,\text{ GeV }$.

Fig. 28

The radial part of the wave function $\psi ({\tilde{r}})$ for the $J/\psi (1S)$ (left panel) and $\psi '(2S)$ (right panel) mesons as a solution of the Schrödinger equation for five distinct $c-\bar{c}$ interaction potentials described in the Appendix A

Fig. 29

The radial part of the wave function $\psi ({\tilde{r}})$ for the $\Upsilon (1S)$ (left panel) and $\Upsilon '(2S)$ (right panel) mesons as a solution of the Schrödinger equation for five distinct $b-\bar{b}$ interaction potentials described in the Appendix A

Appendix B: Spatial quarkonium wave function in the $Q\bar{Q}$ rest frame

The spatial part of the quarkonium wave function satisfies the Schrödinger equation [8]

$$\begin{aligned} \left( -\frac{\Delta }{2\mu }+V({\tilde{r}})\right) \Psi _{nlm}(\vec {\tilde{r}}) =E_{nl}\Psi _{nlm}(\vec {\tilde{r}}) , \quad \mu =\frac{m_Q}{2} ,\nonumber \\ \end{aligned}$$

(B1)

where $\mu$ is the reduced mass of the $Q\bar{Q}$ pair, and the operator $\Delta$ acts on the coordinate ${\tilde{r}}$ and has the following form

$$\begin{aligned} \Delta= & {} \sum \limits _{i=1}^{3}\frac{\partial ^{2}}{\partial x_i^{2}} =\frac{1}{{\tilde{r}}^{2}}\frac{\partial }{\partial {\tilde{r}}} \left( {\tilde{r}}^{2} \frac{\partial }{\partial {\tilde{r}}}\right) \nonumber \\&+\frac{1}{{\tilde{r}}^{2}\sin \theta }\, \frac{\partial }{\partial \theta }\left( \sin \theta \, \frac{\partial }{\partial \theta }\right) +\frac{1}{{\tilde{r}}^{2}\sin ^{2}\theta } \,\frac{\partial ^{2}}{\partial \varphi ^{2}} . \end{aligned}$$

(B2)

Factorizing the spatial wave function into the radial and angular parts,

$$\begin{aligned} \Psi _{nlm}(\vec {\tilde{r}})=\psi _{nl}({\tilde{r}})\, Y_{lm}(\theta ,\varphi ) \end{aligned}$$

(B3)

the Schrödinger equation (B1) with (B2) can be expressed as the following two equations,

$$\begin{aligned}&\frac{1}{{\tilde{r}}}\frac{\partial ^{2}}{\partial {\tilde{r}}^{2}} \bigl ({\tilde{r}}\psi ({\tilde{r}})\bigr )\nonumber \\&\quad +m_Q\bigl (E-V({\tilde{r}})\bigr )\psi ({\tilde{r}}) =\frac{l(l+1)}{{\tilde{r}}^{2}}\psi ({\tilde{r}})\nonumber \\&\frac{1}{\sin \theta }\,\frac{\partial }{\partial \theta } \left( \sin \theta \,\frac{\partial Y(\theta ,\varphi )}{\partial \theta }\right) \nonumber \\&\quad +\frac{1}{\sin ^{2}\theta }\,\frac{\partial ^{2} Y(\theta ,\varphi )}{\partial \varphi ^{2}}=- l(l+1) Y(\theta ,\varphi ) \end{aligned}$$

(B4)

with $l=0$ for S-wave states, $l=1$ for P-waves, etc. The first differential equation for the radial wave function $\psi ({\tilde{r}})$ in the $Q\bar{Q}$ rest frame can be rewritten in a more convenient form

$$\begin{aligned}&\frac{\partial ^{2}u({\tilde{r}})}{\partial {\tilde{r}}^{2}} =(V_\mathrm{eff}\bigl ({\tilde{r}}) -\epsilon _Q\bigr )u({\tilde{r}}),\nonumber \\&V_\mathrm{eff}({\tilde{r}}) =m_QV({\tilde{r}})+\frac{l(l+1)}{{\tilde{r}}^{2}} ,\quad \epsilon _Q = m_Q\,E. \end{aligned}$$

(B5)

where the new radial wave function $u({\tilde{r}})$ is related to $\psi ({\tilde{r}})$ satisfying the following normalization,

$$\begin{aligned} u({\tilde{r}})= & {} \sqrt{4\pi }\,{\tilde{r}}\psi ({\tilde{r}}) ,\nonumber \\&\times \int \limits _0^{\infty }|u({\tilde{r}})|^2\mathrm {d}{\tilde{r}}=1, \quad \int |\psi ({\tilde{r}})|^2\mathrm {d}^3 {\tilde{r}}=1 . \end{aligned}$$

(B6)

The Schrödinger equation (B5) can be solved numerically, e.g. as a special case of the second-order differential equation by means of the Numerov method [108] or converting this equation into a set of the first-order differential equations by means of the Runge–Kutta method [109], for each of the five distinct $Q-\bar{Q}$ interaction potentials discussed in Appendix A. The numerical results for the radial wave function $\psi ({\tilde{r}})$ generated by various $c-\bar{c}$ interaction potentials are shown in Fig. 28 for the $J/\psi (1S)$ (left panel) and $\psi '(2S)$ (right panel) states. The corresponding results for the $\Upsilon (1S)$ and $\Upsilon '(2S)$ radial wave functions are depicted in Fig. 29.

One can see that the variation in the results for $\psi ({\tilde{r}})$, using various interaction potentials, increases towards small ${\tilde{r}}$ in the region where a Coulomb-like behavior of potentials becomes important. The enhanced sensitivity of numerical results to the choice of the $Q-\bar{Q}$ interaction potential appears especially for the 2S radially-excited charmonium state due to the nodal structure of the corresponding radial wave function.

Appendix C: Expressions for amplitudes $\mathcal {A}_{L,R}$

The resulting expressions for the amplitudes of quarkonia photo- and electroproduction in the polarised photon-nucleon scattering read [8]

$$\begin{aligned}&\mathrm {Im}\mathcal {A}_{L}(x,Q^2)=\int \limits _0^1\mathrm {d}z\int \mathrm {d}^2r\,\Sigma _{L}(z,r;Q^2)\,\sigma _{q\bar{q}}(x,r) , \nonumber \\ \Sigma _{L}= & {} Z_q\,\frac{\sqrt{N_c\alpha _{em}}}{2\pi \sqrt{2}} \,4Qz(1-z)K_0(\varepsilon r)\nonumber \\&\times \int p_T\mathrm {d}p_TJ_0(p_Tr) \Psi _{V}(z,p_T) \frac{m_Tm_L+m_Q^2}{m_Q(m_T+m_L)} ,\nonumber \\ \end{aligned}$$

(C1)

for a longitudinally polarised photon,¹ and

$$\begin{aligned} \mathrm {Im}\mathcal {A}_{T}(x,Q^2)= & {} \int \limits _0^1\mathrm {d}z\int \mathrm {d}^2r \left[ \Sigma ^{(1)}_{T}(z, r;Q^2)\sigma _{q\bar{q}}(x,r)\right. \nonumber \\&\left. + \Sigma ^{(2)}_{T}(z, r;Q^2)\frac{\mathrm {d}\sigma _{q\bar{q}}(x,r)}{\mathrm {d}r}\right] , \nonumber \\ \Sigma ^{(1)}_{T}= & {} Z_q\,\frac{\sqrt{N_c\alpha _{em}}}{2\pi \sqrt{2}} \,2K_0(\varepsilon r)\nonumber \\&\quad \times \int \mathrm {d}p_TJ_0(p_Tr)\Psi _{V}(z,p_T) p_T\nonumber \\&\quad \times \frac{m_T^2+m_Tm_L-2p_T^2z(1-z)}{m_T+m_L} , \nonumber \\ \Sigma ^{(2)}_{T}= & {} Z_q\,\frac{\sqrt{N_c\alpha _{em}}}{2\pi \sqrt{2}}\,2K_0(\varepsilon r)\nonumber \\&\quad \times \int \mathrm {d}p_T J_1(p_Tr) \Psi _{V}(z,p_T)\frac{p_T^2}{2}\nonumber \\&\quad \times \frac{m_T+m_L+m_T(1-2z)^2}{m_T(m_T+m_L)} , \end{aligned}$$

(C2)

for a transversely polarised photon. In the above formulas,

$$\begin{aligned} m_T^2 = m_Q^2 + p_T^2 , \quad m_L^2 = 4m_Q^2\,z(1-z) , \end{aligned}$$

(C3)

such that the meson mass squared reads

$$\begin{aligned} M_V^2 = \frac{m_T^2}{z(1-z)} . \end{aligned}$$

(C4)