1 Introduction

The diffractive exclusive electroproduction of vector mesons has been a very active field of study at the HERA collider at DESY. Especially it has served as a testbed for the perturbative QCD description of the Pomeron exchange mechanism which drives diffractive and elastic processes [1, 2]. New precise data, albeit in a different kinematic regime, are expected to be taken at the Brookhaven electron-ion collider (EIC) in a not too distant future [3, 4]. In this work we focus on the electroproduction of \(\rho \)-mesons, i.e. on the process \(\gamma ^* p \rightarrow \rho p\) at large virtuality \(Q^2\) of the photon and large \(\gamma ^* p\)-system center-of-mass energy W, with \( {x \simeq Q^2/W^2 \ll 1}\).

The large photon virtuality \(Q^2\) justifies the use of perturbation theory. A diagrammatic representation of the perturbative QCD factorization structure of the forward amplitude is shown in Fig. 1. At large \(Q^2\) the transverse internal motion of quark and antiquark in the \(\rho \)-meson can be integrated out, and all information is contained in the distribution amplitude (DA). The transverse momenta of gluons in the proton must be fully taken into account by utilizing the unintegrated gluon distribution (UGD). In the BFKL-approach the UGD is obtained by convoluting the BFKL two-gluon Green’s function with the proton impact factor (IF).

The IF for the transition \(\gamma ^*(\lambda _\gamma ) \rightarrow \rho (\lambda _\rho )\) depends on polarizations of photon and vector mesons. For the longitudinal polarizations \(\lambda _\gamma = \lambda _\rho = 0\) the dominance of small dipoles is evident, and the standard leading-twist distribution amplitude appears. In the case of the transverse polarizations the perturbative QCD factorization remains valid, but higher-twist DAs are needed [5, 6]. The difference between the impact factors makes the polarization dependence of diffractive production a sensitive probe of the UGD [7, 8]. Previously, for the case of \(\rho \)-meson production ratios of the forward amplitudes have been calculated in Ref. [8] and compared to data from HERA. For the case of \(\phi \) mesons the effect of a finite strange quark mass on cross sections has been discussed in Ref. [9], where also relations of the so-called genuine higher-twist DAs to weighted integrals over light-front wave functions have been given.

In this work, we wish to extend the efforts of [7,8,9] to the calculation of polarized \(\rho \)-meson production cross sections (as opposed to ratios of amplitudes) for a variety of unintegrated gluon distributions, and compare the results to HERA data. We will also give predictions for the kinematics relevant for the Electron-Ion Collider EIC.

Fig. 1
figure 1

Diagrammatic representation of the amplitude for the exclusive emission of a \(\rho \) meson in high-energy factorization. The off-shell impact factor (upper part) is built up as a collinear convolution between the hard factor for the photon splitting to a dipole and the non-perturbative \(\rho \)-meson DA (sea-green blob). The UGD is given by the convolution of the BFKL gluon Green’s function (yellow blob) and the non-perturbative proton impact factor (red blob)

Full size image

2 Polarized \(\rho \)-meson leptoproduction

The H1 and ZEUS collaborations have provided extended analyses of the helicity structure in the hard exclusive production of the \(\rho \) meson in ep collisions through the subprocess

$$\begin{aligned} \gamma ^*(\lambda _\gamma )p\rightarrow \rho (\lambda _\rho )p. \end{aligned}$$
(1)

Here \(\lambda _\rho \) and \(\lambda _\gamma \) represent the meson and photon helicities, respectively, and can take the values 0 (longitudinal polarization) and \(\pm 1\) (transverse polarizations). The helicity amplitudes \(T_{\lambda _\rho \lambda _\gamma }\) extracted at HERA [10] exhibit the hierarchy

$$\begin{aligned} T_{00} \gg T_{11} \gg T_{10} \gg T_{01} \gg T_{-11} , \end{aligned}$$
(2)

that follows from the dominance of a small-size dipole scattering mechanism, as discussed first in Ref. [11], see also [12]. As we previously mentioned, the H1 and ZEUS collaborations have analyzed data in different ranges of \(Q^2\) and W. In what follows we will refer only to the H1 ranges (3),

$$\begin{aligned} 2.5\,{\hbox {GeV}^2}< & {} Q^2<60 \,{ \hbox {GeV}^2},\nonumber \\ 35\, \text {GeV}< & {} W < 180\,\text {GeV}, \end{aligned}$$
(3)

illustrating the framework to probe the \(\rho \)-meson leptoproduction as a way to test the UGDs [7, 13, 14] where polarized cross sections will be the key observables.

2.1 Theoretical setup

In the high-energy regime, \(s\equiv W^2\gg Q^2\gg \Lambda _{\text {QCD}}^2\), which implies small \(x = Q^2/W^2\), the forward helicity amplitude for the \(\rho \)-meson electroproduction can be written, in high-energy factorization (also known as \(k_T\)-factorization), as the convolution of the \(\gamma ^*\rightarrow \rho \) IF, \(\Phi ^{\gamma ^*(\lambda _\gamma )\rightarrow \rho (\lambda _\rho )}(\kappa ^2,Q^2)\), with the UGD, \({{\mathcal {F}}}(x,\kappa ^2)\). Its expression reads

$$\begin{aligned} T_{\lambda _\rho \lambda _\gamma }(s,Q^2)= & {} \frac{is}{(2\pi )^2}\int \dfrac{d^2\kappa }{(\kappa ^2)^2}\Phi ^{\gamma ^*(\lambda _\gamma )\rightarrow \rho (\lambda _\rho )} (\kappa ^2,Q^2) {{\mathcal {F}}}(x,\kappa ^2),\nonumber \\ x= & {} \frac{Q^2}{s}. \end{aligned}$$
(4)

Defining \(\alpha = \frac{\kappa ^2}{Q^2}\) and \(B =~2\pi \alpha _s \frac{e}{\sqrt{2}}f_\rho \), the expression for the IFs takes the following forms (see Ref. [5] for the derivation):

  • longitudinal case

    $$\begin{aligned} \Phi ^{\gamma _L\rightarrow \rho _L}(\kappa ,Q;\mu ^2) = 2 B\frac{\sqrt{N_c^2-1}}{Q\,N_c} \int ^{1}_{0}dy\, \varphi _1(y;\mu ^2)\left( \frac{\alpha }{\alpha + y{\bar{y}}}\right) ,\nonumber \\ \end{aligned}$$
    (5)

    where \(N_c\) denotes the number of colors and \(\varphi _1(y;\mu ^2)\) is the twist-2 DA which, up to the second order in the expansion in Gegenbauer polynomials, reads [15]

    $$\begin{aligned} \varphi _1(y; \mu ^2) = 6y{\bar{y}}\left( 1+a_2(\mu ^2)\frac{3}{2} \left( 5 (y-{\bar{y}})^2-1\right) \right) ;\nonumber \\ \end{aligned}$$
    (6)

  • transverse case

    $$\begin{aligned}&\Phi ^{\gamma _T\rightarrow \rho _T}(\alpha ,Q;\mu ^2)= \dfrac{(\epsilon _{\gamma }\cdot \epsilon ^{*}_{\rho }) \, 2 B m_{\rho } \sqrt{N_c^2-1}}{2 N_c Q^2}\nonumber \\&\quad \times \left\{ -\int ^{1}_{0} dy \frac{\alpha (\alpha +2 y {\bar{y}})}{y{\bar{y}} (\alpha +y{\bar{y}})^2}\right. \nonumber \\&\quad \times \left[ (y-{\bar{y}})\varphi _1^T(y;\mu ^2) + \varphi _A^T(y;\mu ^2)\right] \nonumber \\&\quad +\int ^{1}_{0}dy_2\int ^{y_2}_{0}dy_1 \frac{y_1{\bar{y}}_1\alpha }{\alpha +y_1{\bar{y}}_1}\nonumber \\&\quad \times \left[ \frac{2-N_c/C_F}{\alpha (y_1+{\bar{y}}_2) +y_1{\bar{y}}_2} -\frac{N_c}{C_F}\frac{1}{y_2 \alpha +y_1(y_2-y_1)}\right] \nonumber \\&\quad \left. \times M(y_1,y_2;\mu ^2)-\!\int ^{1}_{0}dy_2\!\int ^{y_2}_{0}dy_1\! \left[ \frac{2+N_c/C_F}{{\bar{y}}_1} +\frac{y_1}{\alpha +y_1{\bar{y}}_1}\right. \right. \nonumber \\&\quad \!\left. \left. \times \left( \frac{(2-N_c/C_F) y_1\alpha }{\alpha (y_1+{\bar{y}}_2)+y_1{\bar{y}}_2}-2\right) -\frac{N_c}{C_F}\frac{(y_2-y_1){\bar{y}}_2}{{\bar{y}}_1}\right. \right. \nonumber \\&\quad \left. \left. \times \frac{1}{\alpha {\bar{y}}_1+(y_2-y_1){\bar{y}}_2} \right] \!S(y_1,y_2;\mu ^2)\right\} , \end{aligned}$$
    (7)

    where \(C_F=\frac{N_c^2-1}{2N_c}\), while the functions \(M(y_1,y_2;\mu ^2)\) and \(S(y_1,y_2;\mu ^2)\) are defined in Eqs. (12)–(13) of Ref. [6] and are combinations of the twist-3 DAs \(B(y_1,y_2;\mu ^2)\) and \(D(y_1,y_2;\mu ^2)\) (see Ref. [15]), given by

    $$\begin{aligned} B(y_1,y_2;\mu ^2)&=-5040 y_1 {\bar{y}}_2 (y_1-{\bar{y}}_2) (y_2-y_1) \nonumber ,\\ D(y_1,y_2;\mu ^2)&=-360 y_1{\bar{y}}_2(y_2-y_1)\nonumber \\&\quad \times \left( 1+\frac{\omega ^{A}_{\{1,0\}}(\mu ^2)}{2} \left( 7\left( y_2-y_1\right) -3\right) \right) . \end{aligned}$$
    (8)

In Eqs. (6) and (8) the functional dependence of \(a_2\), \(\omega ^{A}_{\{1,0\}}\), \(\zeta ^{A}_{3\rho }\), and \(\zeta ^{V}_{3\rho }\) on the factorization scale \(\mu ^2\) can be determined from the corresponding known evolution equations [15], using some suitable initial condition at a scale \(\mu _0\).

Note that the \(\kappa \)-dependence of the IFs is different in the cases of longitudinal and transverse polarizations and this poses a strong constraint on the \(\kappa \)-dependence of the UGD in the HERA energy range. The main point will be to demonstrate, considering different models of UGD, that the uncertainties of the theoretical description do not prevent us from some, at least qualitative, conclusions about the shape of the UGD in \(\kappa \).

The DAs \(\varphi ^T_1(y;\mu ^2)\) and \(\varphi ^T_A(y;\mu ^2)\) in Eq. (7) encompass both genuine twist-3 and Wandzura–Wilczek (WW) contributions [6, 15]. The former are related to \(B(y_1,y_2;\mu ^2)\) and \(D(y_1,y_2;\mu ^2)\); the latter are those obtained in the approximation in which \(B(y_1,y_2;\mu ^2)=D(y_1,y_2;\mu ^2)=0\), and in this case readFootnote 1

$$\begin{aligned}&\varphi ^{T\;WW}_A(y;\mu ^2)= \frac{1}{2}\left[ -{{\bar{y}}} \int _0^{y}\,dv \frac{\varphi _1(v;\mu ^2)}{{{\bar{v}}}} - y \int _{y}^1\,dv \frac{\varphi _1(v;\mu ^2)}{v} \right] , \nonumber \\&\varphi ^{T\;WW}_1(y;\mu ^2)= \frac{1}{2}\left[ -{{\bar{y}}} \int _0^{y}\,dv \frac{\varphi _1(v;\mu ^2)}{{{\bar{v}}}} + y \int _{y}^1\,dv \frac{\varphi _1(v;\mu ^2)}{v} \right] .\nonumber \\ \end{aligned}$$
(9)

The other interesting point will be the extension of information, collected by the helicity-amplitude analysis, reachable from the calculation of the cross section. As a matter of fact, the expressions for the polarized cross sections \(\sigma _L\) and \(\sigma _T\), calculated using Eqs. (5), (6) in Eqs. (4) and (7), (8) in Eq. (4), respectively, are

$$\begin{aligned} \sigma _L\,(\gamma ^*\,p \rightarrow \rho \,p)= & {} \frac{1}{16 \pi b(Q^2)}\left| \frac{T_{00}(s, t = 0)}{W^2}\right| ^2\,, \end{aligned}$$
(10)
$$\begin{aligned} \sigma _T\,(\gamma ^*\,p \rightarrow \rho \,p)= & {} \frac{1}{16 \pi b(Q^2)}\left| \frac{T_{11}(s, t = 0)}{W^2}\right| ^2, \end{aligned}$$
(11)

where \(b(Q^2)\) in Eqs. (10) and (11) is the diffraction slope, for which we adopt the parametrization [16]:

$$\begin{aligned} b(Q^2) = \beta _0 - \beta _1\,\log \left[ \frac{Q^2+m_\rho ^2}{m^2_{J/\psi }}\right] +\frac{\beta _2}{Q^2+m_\rho ^2}, \end{aligned}$$
(12)

fixing the values of the constants as follows: \(\beta _0 = 6.5\) \(\hbox {GeV}^{-2}\), \(\beta _1 = 1.2\) \(\hbox {GeV}^{-2}\) and \(\beta _2 = 1.6\). In Fig. 2, we compare the parametrization of Eq. (12) with the data of the H1-collaboration from Ref. [10]. Above, we took into account only the helicity conserving amplitudes. These dominate the diffractive peak at small t, and hence also the integrated cross section of interest. While the impact factors for the helicity flip transitions \(\gamma ^*(T)\rightarrow \rho (L)\) with \(\Delta \lambda = \pm 1\) and \(\gamma ^*(T)\rightarrow \rho (T')\) with \(\Delta \lambda = \pm 2\), are available within the higher-twist factorization approach [5, 6], a discussion of the \(t-\)dependent cross section or the polarization density matrix of \(\rho \)-mesons would require also a generalization to off-forward UGDs and goes beyond the scope of this work.

The advantage of considering polarized cross sections relies on the possibility to constrain not only the shape and the behavior of results, but also the normalization, now essential and which is clearly irrelevant for the evaluation of the helicity-amplitude ratio \(T_{11}/T_{00}\). Although all UGDs described in Sect. 3 present fixed values for their parameters, the ABIPSW model, one of the first UGD parametrization adopted in phenomenological analysis [17], was defined up to the overall normalization. As regards the helicity-amplitude ratio \(T_{11}/T_{00}\), results for the ABIPSW UGD model show a fair agreement with data [6], this suggesting that the shape, controlled by the M parameter (see Eq. (13)), has been already guessed. The best value of the normalization parameter can be obtained via a simple global fit to experimental data of both polarized cross sections, \(\sigma _L\) and \(\sigma _T\).

Fig. 2
figure 2

\(Q^2\)-dependence of the diffractive slope for the exclusive \(\rho \)-meson leptoproduction (Eq. (12)) compared to the H1 data from Ref. [10]

Full size image

3 UGD models used in the present study

We have considered a selection of several models of UGD, without pretension to exhaustive coverage, but with the aim of comparing (sometimes radically) different approaches. We refer the reader to the original papers for details on the derivation of each model and limit ourselves to presenting here just the functional form \({{\mathcal {F}}}(x,\kappa ^2)\) of the UGD as we implemented it in the numerical analysis.

3.1 An x-independent model (ABIPSW)

The simplest UGD model is x-independent and merely coincides with the proton impact factor [6]:

$$\begin{aligned} {{\mathcal {F}}}(x,\kappa ^2)= \frac{A}{(2\pi )^2\,M^2} \left[ \frac{\kappa ^2}{M^2+\kappa ^2}\right] , \end{aligned}$$
(13)

where M corresponds to the non-perturbative hadronic scale, fixed as \(M = 1\) GeV. The constant A is irrelevant when we consider the ratio \(T_{11}/T_{00}\) for the \(\rho \)-meson leptoproduction, but it becomes essential to calculate cross sections. Therefore, we determined it by a global fit to all available data for polarized cross sections, getting \(A = 148.14 \text { GeV}^2\), with a relative uncertainty below \(0.05\%\).

3.2 Gluon momentum derivative

This UGD is given by

$$\begin{aligned} {{\mathcal {F}}}(x, \kappa ^2) = \frac{dxg(x, \kappa ^2)}{d\ln \kappa ^2} \end{aligned}$$
(14)

and encompasses the collinear gluon density \(g(x, \mu _F^2)\), taken at \(\mu _F^2=\kappa ^2\). It is based on the obvious requirement that, when integrated over \(\kappa ^2\) up to some factorization scale, the UGD must give the collinear gluon density. We have employed the CT14 parametrization [18], using the appropriate cutoff \(\kappa _{\text {min}} = 0.3\) GeV (see Section III A of Ref. [8] for further details).

Fig. 3
figure 3

\(Q^2\)-dependence of the longitudinally polarized cross section, \(\sigma _L\), for all the considered UGD models, at \(W = 75\) GeV together with the HERA data (left upper panel) and at \(W = 20, 30, 50\) GeV for EIC (the remaining panels). Uncertainty bands represent the effect of varying \(a_2(\mu _0 = 1\,\) GeV) between 0.0 and 0.6

Full size image
Fig. 4
figure 4

\(Q^2\)-dependence of the transversely polarized cross section, \(\sigma _T\), for all the considered UGD models, at \(W = 75\) GeV together with the HERA data (left upper panel) and at \(W = 20, 30, 50\) GeV for EIC (the remaining panels). Uncertainty bands represent the effect of varying \(a_2(\mu _0 = 1\,\)GeV) between 0.0 and 0.6

Full size image
Fig. 5
figure 5

\(Q^2\)-dependence of the polarized cross-section ratio, \(\sigma _L/\sigma _T\), for all the considered UGD models, at \(W = 75\) GeV together with the HERA data (left upper panel) and at \(W = 20, 30, 50\) GeV for EIC (the remaining panels). Uncertainty bands represent the effect of varying \(a_2(\mu _0 = 1\,\)GeV) between 0.0 and 0.6

Full size image
Fig. 6
figure 6

\(Q^2\)-dependence of the longitudinally polarized cross section, \(\sigma _L\), for GBW and IN UGD models, at \(W = 75\) GeV together with the HERA data (left upper panel) and at \(W = 20, 30, 50\) GeV EIC (the remaining panels). Uncertainty bands represent the effect of varying \(a_2(\mu _0 = 1\,\)GeV) between 0.0 and 0.6

Full size image
Fig. 7
figure 7

\(Q^2\)-dependence of the transversely polarized cross section, \(\sigma _T\), for GBW and IN UGD models, at \(W = 75\) GeV together with the HERA data (left upper panel) and at \(W = 20, 30, 50\) GeV relevant for EIC (the remaining panels). Uncertainty bands give the effect of varying \(a_2(\mu _0 = 1\,\)GeV) between 0.0 and 0.6. Full, WW and genuine contributions are shown separately

Full size image
Fig. 8
figure 8

\(Q^2\)-dependence of the polarized cross-section ratio, \(\sigma _L/\sigma _T\), for the GBW and IN UGD models, at \(W = 75\) GeV HERA (left upper panel) and at \(W = 20, 30, 50\) GeV relevant for EIC (the remaining panels). Uncertainty bands represent the effect of varying \(a_2(\mu _0 = 1\,\)GeV) between 0.0 and 0.6. Full and WW contributions are shown separately

Full size image

3.3 Ivanov–Nikolaev (IN) UGD: a soft-hard model

The UGD proposed in Ref. [19] is developed with the purpose of probing different regions of the transverse momentum. In the large-\(\kappa \) region, DGLAP parametrizations for \(g(x, \kappa ^2)\) are employed. Moreover, for the extrapolation of the hard gluon densities to small \(\kappa ^2\), an Ansatz is made [20], which describes the color gauge invariance constraints on the radiation of soft gluons by color singlet targets. The gluon density at small \(\kappa ^2\) is supplemented by a non-perturbative soft component, according to the color-dipole phenomenology.

This model of UGD has the following two-component form:

$$\begin{aligned} {{\mathcal {F}}}(x,\kappa ^2)= {{\mathcal {F}}}^{(B)}_\text {soft}(x,\kappa ^2) {\kappa _{s}^2 \over \kappa ^2 +\kappa _{s}^2} + {{\mathcal {F}}}_\text {hard}(x,\kappa ^2) {\kappa ^2 \over \kappa ^2 +\kappa _{h}^2},\nonumber \\ \end{aligned}$$
(15)

where \(\kappa _{s}^2 = 3\) \(\hbox {GeV}^2\) and \(\kappa _{h}^2 = [1 + 0.047\log ^2(1/x)]^\frac{1}{2}\) \(\hbox {GeV}^2\).

The soft term reads

$$\begin{aligned} {{\mathcal {F}}}^{(B)}_\text {soft}(x,\kappa ^2) = a_\text {soft} C_{F} N_{c} {\alpha _{s}(\kappa ^2) \over \pi } \left( {\kappa ^2 \over \kappa ^2 +\mu _\text {soft}^{2}}\right) ^2 V_{\text{ N }}(\kappa ),\nonumber \\ \end{aligned}$$
(16)

with \(\mu _\text {soft} = 0.1\) GeV. The parameter \(a_\text {soft} = 2\) gives a measure of how important is the soft part compared to the hard one. On the other hand, the hard component reads

$$\begin{aligned} {{\mathcal {F}}}_\text {hard}(x,\kappa ^2)= & {} {{\mathcal {F}}}^{(B)}_{\text {pt}}(\kappa ^2) {{{\mathcal {F}}}_{\text {pt}}(x,Q_{c}^{2}) \over {{\mathcal {F}}}_{\text {pt}}^{(B)}(Q_{c}^{2})} \theta (Q_{c}^{2}-\kappa ^{2})\nonumber \\&+{{\mathcal {F}}}_{\text {pt}}(x,\kappa ^2) \theta (\kappa ^{2}-Q_{c}^{2}), \end{aligned}$$
(17)

where \({{\mathcal {F}}}_{\text {pt}}(x, \kappa ^2)\) is related to the collinear gluon parton distribution function (PDF) as in Eq. (14) and \(Q_{c}^2 = 3.26\) \(\hbox {GeV}^2\) (see Section III A of Ref. [8] for further details). We refer to Ref. [19] for the expressions of the vertex function \(V_{\text{ N }}(\kappa )\) and of \({{\mathcal {F}}}^{(B)}_{\text {pt}}(\kappa ^2)\). Another relevant feature of this model is given by the choice of the coupling constant. In this regard, the infrared freezing of strong coupling at leading order (LO) is imposed by fixing \(\Lambda _\text {QCD} = 200\) MeV:

$$\begin{aligned} \alpha _s(\mu ^2) = \text {min} \left\{ 0.82, \, \frac{4 \pi }{\beta _0 \log \left( \frac{\mu ^2}{\Lambda ^2_\text {QCD}}\right) }\right\} . \end{aligned}$$
(18)

We wish to stress that this model was successfully tested in the unpolarized electroproduction of vector mesons at HERA.

3.4 Hentschinski–Sabio Vera–Salas (HSS) model

This model, originally used in the study of DIS structure functions [21], takes the form of a convolution between the BFKL gluon Green’s function and a LO proton impact factor. It has been employed in the description of single-bottom quark production at the LHC [22], to investigate the photoproduction of \(J/\Psi \) and \(\Upsilon \) in [23,24,25] and to study the forward Drell–Yan invariant-mass distribution [26]. We implemented the formula given in Ref. [22] (up to a \(\kappa ^2\) overall factor needed to match our definition), which reads

$$\begin{aligned} {{\mathcal {F}}}(x, \kappa ^2, M_h)= & {} \int _{-\infty }^{\infty } \frac{d\nu }{2\pi ^2}\ {{\mathcal {C}}} \ \frac{\Gamma (\delta - i\nu -\frac{1}{2})}{\Gamma (\delta )}\nonumber \\&\times \left( \frac{1}{x}\right) ^{\chi \left( \frac{1}{2}+i\nu \right) } \left( \frac{\kappa ^2}{Q^2_0}\right) ^{\frac{1}{2}+i\nu } \nonumber \\&\times \left\{ 1 +\frac{{\bar{\alpha }}^2_s \beta _0 \chi _0\left( \frac{1}{2} +i\nu \right) }{8 N_c}\log \left( \frac{1}{x}\right) \right. \nonumber \\&\quad \times \left. \left[ -\psi \left( \delta -\frac{1}{2} - i\nu \right) -\log \frac{\kappa ^2}{M_h^2}\right] \right\} ,\nonumber \\ \end{aligned}$$
(19)

where \(\beta _0=\frac{11 N_c-2 N_f}{3}\), with \(N_f\) the number of active quarks (put equal to four in the following), \({\bar{\alpha }}_s = \dfrac{\alpha _s\left( \mu ^2\right) N_c}{\pi }\), with \(\mu ^2 = Q_0 M_h\), and \(\chi _0(\frac{1}{2} + i\nu )\equiv \chi _0(\gamma ) = 2\psi (1) - \psi (\gamma ) - \psi (1-\gamma )\) is (up to the factor \({{\bar{\alpha }}}_s\)) the LO eigenvalue of the BFKL kernel, with \(\psi (\gamma )\) the logarithmic derivative of the Euler Gamma function. Here, \(M_h\) plays the role of the hard scale which in our case can be identified with the photon virtuality, \(\sqrt{Q^2}\). In Eq. (19), \(\chi (\gamma )\) (with \(\gamma = \frac{1}{2} + i\nu \)) is the NLO eigenvalue of the BFKL kernel, collinearly improved and BLM optimized. It reads

$$\begin{aligned} \chi (\gamma )= & {} {\bar{\alpha }}_s\chi _0(\gamma )+{\bar{\alpha }}^2_s\chi _1(\gamma )\nonumber \\&-\frac{1}{2}{\bar{\alpha }}^2_s\chi ^\prime _0(\gamma )\,\chi _0(\gamma ) + \chi _{RG}({\bar{\alpha }}_s, \gamma ), \end{aligned}$$
(20)

with \(\chi _1(\gamma )\) and \(\chi _{RG}({\bar{\alpha }}_s, \gamma )\) given in Section 2 of Ref. [22].

This UGD model is characterized by a peculiar parametrization for the proton impact factor, whose expression is

$$\begin{aligned} \Phi _p(q, Q^2_0) = \frac{{{\mathcal {C}}}}{2\pi \Gamma (\delta )} \left( \frac{q^2}{Q^2_0}\right) ^\delta e^{-\frac{q^2}{Q^2_0}}, \end{aligned}$$
(21)

which depends on three parameters \(Q_0\), \(\delta \) and \({{\mathcal {C}}}\) which were fitted to the combined HERA data for the \(F_2(x)\) proton structure function. We adopted here the so-called kinematically improved values (see Section III A of Ref. [8] for further details) and given by

$$\begin{aligned} Q_0 = 0.28\,\text {GeV}, \quad \delta = 6.5, \quad {{\mathcal {C}}} = 2.35. \end{aligned}$$
(22)

3.5 Golec–Biernat–Wüsthoff (GBW) UGD

This UGD parametrization derives from the effective dipole cross section \({\hat{\sigma }}(x,r)\) for the scattering of a \(q{\bar{q}}\) pair off a nucleon [27],

$$\begin{aligned} {\hat{\sigma }}(x, r^2) = \sigma _0 \left\{ 1-\exp \left( -\frac{r^2}{4R^2_0(x)} \right) \right\} , \end{aligned}$$
(23)

through a reverse Fourier transform of the expression

$$\begin{aligned}&\sigma _0 \left\{ 1-\exp \left( -\frac{r^2}{4R^2_0(x)} \right) \right\} ={2 \pi \over N_c} \int \frac{d^2\kappa }{\kappa ^4} \alpha _s {{\mathcal {F}}}(x,\kappa ^2) \nonumber \\&\quad \times \left( 1-\exp (i \mathbf {\kappa }\cdot \mathbf {r})\right) \left( 1-\exp (-i \mathbf {\kappa } \cdot \mathbf {r})\right) , \end{aligned}$$
(24)
$$\begin{aligned}&\alpha _s {{\mathcal {F}}}(x,\kappa ^2)= N_c \kappa ^4 \sigma _0 \frac{R^2_0(x)}{4 \pi ^2} e^{-\kappa ^2 R^2_0(x)}, \end{aligned}$$
(25)

with

$$\begin{aligned} R^2_0(x) = \frac{1}{{\text {GeV}}^2} \left( \frac{x}{x_0}\right) ^{\lambda _p} \end{aligned}$$
(26)

and the following values

$$\begin{aligned} \sigma _0 = 23.03\,\text {mb}, \quad \lambda _p = 0.288, \quad x_0 = 3.04 \cdot 10^{-4}. \end{aligned}$$
(27)

The normalization \(\sigma _0\) and the parameters \(x_0\) and \(\lambda _p > 0\) of \(R^2_0(x)\) have been determined by a global fit to \(F_2(x)\) in the region \(x < 0.01\).

3.6 Watt–Martin–Ryskin (WMR) model

This model is based on the idea that the \(\kappa \) dependence of the UGD comes from the last step of the evolution ladder.

The UGD introduced in Ref. [28] reads

$$\begin{aligned}&{{\mathcal {F}}}(x,\kappa ^2,\mu ^2) = T_g(\kappa ^2,\mu ^2)\,\frac{\alpha _s(\kappa ^2)}{2\pi }\,\int _x^1\!dz\;\nonumber \\&\quad \times \left[ \sum _q P_{gq}(z)\,\frac{x}{z}q\left( \frac{x}{z}, \kappa ^2\right) + P_{gg}(z)\,\frac{x}{z}g\left( \frac{x}{z},\kappa ^2\right) \,\Theta \left( \frac{\mu }{\mu +\kappa }-z\right) \,\right] \,, \nonumber \\ \end{aligned}$$
(28)

where the term

$$\begin{aligned} T_g(\kappa ^2,\mu ^2)= & {} \exp \!\left( -\int _{\kappa ^2}^{\mu ^2}\!d\kappa _t^2\, \frac{\alpha _s(\kappa _t^2)}{2\pi }\,\!\right. \nonumber \\&\times \left. \left( \int _{z^\prime _{{\text {min}}}}^{z^\prime _{{\text {max}}}} \!dz^\prime \;z^\prime \,P_{gg}(z^\prime ) + N_f\!\!\int _0^1\!dz^\prime \,P_{qg}(z^\prime ) \right) \right) ,\nonumber \\ \end{aligned}$$
(29)

gives the probability of evolving from the scale \(\kappa \) to the scale \(\mu \) without parton emission. Here \(z^\prime _{\mathrm {max}}\equiv 1-z^\prime _{\mathrm {min}}=\mu /(\mu +\kappa _t)\); \(N_f\) is the number of active quarks. This UGD model depends on an extra-scale \(\mu \), which we fixed at Q. The splitting functions \(P_{qg}(z)\) and \(P_{gg}(z)\) are given by

$$\begin{aligned} P_{qg}(z)= & {} T_R\,[z^2 + (1-z)^2], \\ P_{gg}(z)= & {} 2\,C_A \left[ \dfrac{1}{(1-z)_+} + \dfrac{1}{z}- 2 +z(1-z)\right] \\&+ \left( \frac{11}{6}C_A - \frac{N_f}{3}\right) \delta (1 -z), \end{aligned}$$

with the plus-prescription defined as

$$\begin{aligned} \int _{a}^{1} dz \frac{F(z)}{(1-z)_+} = \int _{a}^{1} dz \frac{F(z) - F(1)}{(1-z)} - \int _{0}^{a} dz \frac{F(1)}{(1 -z)}.\nonumber \\ \end{aligned}$$
(30)

3.7 Bacchetta–Celiberto–Radici–Taels (BCRT) distribution

We include in our analysis the unpolarized distribution part of the set of leading-twist T-even transverse-momentum-dependent (TMD) gluon PDFs calculated in Ref. [29] (see also Ref. [30]), suited for studies in wider kinematic ranges at new-generation collider machines, such as the EIC [4], HL-LHC [31], and NICA-SPD [32]. One has

$$\begin{aligned} {{\mathcal {F}}}(x,\kappa ^2) = \kappa ^2 \int _M^{\infty } d M_X \, \rho _X (M_X) \, {\hat{f}}_1^g(x, \kappa ^2; M_X), \end{aligned}$$
(31)

Here, \({\hat{f}}_1^g(x, \kappa ^2; M_X)\) stands for the distribution of unpolarized gluons in an unpolarized proton, calculated in the so-called spectator-model approximation, namely where the remainders of the proton after gluon emission are treated as a single spin-1/2 spectator particle of mass \(M_X\)

$$\begin{aligned} {\hat{f}}_1^g(x, \kappa ^2; M_X)= & {} \Big [ \big ( 2 M x g_1(\kappa ^2) - x (M+M_X) g_2(\kappa ^2) \big )^2 \, \nonumber \\&\times \big [ (M_X - M (1-x) )^2 + \kappa ^2 \big ] \nonumber \\&\quad + \, 2 \kappa ^2 \, (\kappa ^2 + x M_X^2) \, g_2(\kappa ^2)^2\nonumber \\&+ 2 \kappa ^2 M^2 \, (1-x) \, (4 g_1(\kappa ^2)^2 - x g_2(\kappa ^2)^2 ) \Big ] \nonumber \\&\times \Big [ (2 \pi )^3 \, 4 x M^2 \, (x M_X^2 - x (1 - x) M^2 + \kappa ^2)^2 \Big ]^{-1} .\nonumber \\ \end{aligned}$$
(32)

In Eq. (32) M is the proton mass, whereas \(g_{1,2}\) couplings depict the effective proton-gluon-spectator vertex interaction

$$\begin{aligned} g_{1,2}(\kappa ^2) = G_{1,2} \, \frac{(1 - x) \left[ \kappa ^2 + x M_X^2 - x (1 - x) M^2 \right] }{\left[ \kappa ^2 + x M_X^2 - x (1 - x) M^2 + (1 - x) \Lambda _X^2 \right] ^2}. \nonumber \\ \end{aligned}$$
(33)

The spectral function \(\rho _X\) in Eq. (31) weighs \({\hat{f}}_1\) over \(M_X\) in a continuous range and its expression reads

$$\begin{aligned} \rho _X (M_X) = \mu ^{2a} \left( \frac{A}{B + \mu ^{2b}} + \frac{C}{\pi \sigma } e^{-\frac{(M_X - D)^2}{\sigma ^2}} \right) , \end{aligned}$$
(34)

where \(\mu ^2 = M_X^2 -M^2\), and the model parameters in Eqs. (33) and (34) were obtained by a simultaneous fit of the \(\kappa \)-integral of spectator-model unpolarized and helicity functions on NNPDF3.1sx [33] and NNPDFpol1.1 [34] collinear PDFs, respectively, and they encode effective small-x effects coming from the BFKL resummation. In our study we employ values for these parameters obtained by averaging over the whole set of replicas provided (see Section 3 of Ref. [29] for details). They read: \(G_1 = 1.51\) \(\hbox {GeV}^2\), \(G_2 = 0.414\) \(\hbox {GeV}^2\), \(\Lambda _X = 0.472\) GeV, \(A = 6.1\), \(a = 0.82\), \(b = 1.43\), \(C = 371\), \(D = 0.548\) GeV, and \(\sigma = 0.52\) GeV.

We stress that exploratory studies by using this TMD model as a UGD are justified by the fact that we are considering observables for the exclusive \(\rho \)-meson leptoproduction in the forward limit. In a more general, off-forward case one should consider rather models for the gluon generalized parton distribution (GPD), instead of TMDs.

4 Results

All the results shown in this section were obtained by making use of the leptonic-exclusive-amplitudes (LExA) modular interface as implemented in the JETHAD code [35].

We present predictions for the polarized cross sections \(\sigma _L\) and \(\sigma _T\) and their ratio \(\sigma _L/\sigma _T\), as obtained with all the UGD models presented above, and compare them with HERA data. The behavior of the polarized cross sections \(\sigma _L\) and \(\sigma _T\) and the ratio \(\sigma _L/\sigma _T\) [36] in terms of \(Q^2\) for all UGDs, at \(W = 75\) GeV (comparison with HERA data) and at \(W = 20, 30, 50\) GeV (predictions for the EIC), is shown in Figs. 3, 4, and 5, taking into account the variation of the Gegenbauer coefficient \(a_2(\mu _0)\) in the same range used for the helicity-amplitude ratio analysis of Ref. [8]. As can be seen from the figures, the uncertainties due to the choice of UGDs are much bigger than those due to \(a_2\) uncertainty which opens a possibility to test UGD models for the reaction under consideration. The comparison exhibits a partial agreement with the experimental data, where once again none of the proposed models is able to describe the whole \(Q^2\) region. However, we can specify which UGD model is more suitable for the description of the polarized cross sections and which, indeed, for their ratio in the \(Q^2\) intermediate range. On one side, observing Figs. 6 and 7 , the GBW model, considered in its standard definition, i.e. without the evolution of saturation scale (for a detailed discussion about saturation effects see Ref. [37]), and the IN one appear to be the UGD models that allow us to match data for the single polarized cross sections, \(\sigma _L\) and \(\sigma _T\), in the most accurate way. Here the genuine twist-3 contribution is considered. On the other side, although the predictions for the cross section ratio \(\sigma _L/\sigma _T\) in Fig. 8 with the GBW model is quite resonable, if we regard increasing values of \(Q^2 \sim 10\) \(\hbox {GeV}^2\), the IN UGD modelFootnote 2 is able to slightly better catch also the low-\(Q^2\) region of data. We stress, however, that in this low-\(Q^2\) range the validity of our approach to IFs, based on the use of collinear DAs, could be questionable. Therefore our ability to discriminate among UGD models at small virtualities could be limited.

5 Conclusions

We have calculated cross sections for the diffractive electroproduction of \(\rho \) mesons in the energy range of HERA and EIC. We have used the high-energy factorization formalism for the forward amplitude [5, 6], utilizing an empirical parametrization of the diffractive slope to obtain the relevant integrated cross sections.

The impact factors for longitudinal and transverse mesons probe the transverse momentum dependence of the UGD in different ways, so that the polarization dependence of \(\rho \)-production has been proposed as a sensitive probe of the shape of the UGD [7, 8]. The impact factor for longitudinally polarized meson is obtained in terms of the leading twist DA. We find that, for the higher-twist DAs relevant for transverse mesons, the WW contributions dominate. We have investigated the effect of uncertainty related with the form of leading twist DA by varying the Gegenbauer coefficient \(a_{2}\).

We have performed calculations for a representative choices of UGDs available in the literature. These UGDs follow different strategies in their construction. Some explicitly connect to solutions of a specific evolution equation [21, 28], while others stress the presence of a substantial nonperturbative component [19, 27], while being adjusted to giving a good description of proton deep inelastic structure functions.

The latter two UGDs in fact do give the best description of the HERA cross sections. While our use of an empirical parametrization of the diffractive slope introduces an additional model element/uncertainty, the spread of predictions from various UGDs appears to be larger than the uncertainty in the slope. Some UGDs substantially underestimate the total cross section, which may mean that, in the relevant \(\kappa \)-range, their \(\kappa \)-shape and normalization miss some important insights on the nonperturbative proton dynamics.

The cross section ratio \(\sigma _L/\sigma _T\) indeed appears to have potential to discriminate further between UGDs, for the case at hand the IN UGD gives the best description of HERA data.

New data taken at the Electron-Ion Collider EIC also have a potential to further discriminate between the different UGDs. At the lower range of W an investigation of the matching/correspondence to TMD factorization approaches with on-shell partons would be interesting.

It has been recently pointed out how the inclusive diffractive tagging of particles with a heavy transverse mass, such as Higgs bosons [38], heavy-flavored jets [39,40,41] and charmed baryons [42], leads to a fair stabilization of the high-energy series under higher-order corrections. Future studies of reactions featuring the forward emission of those objects will provide us with additional and possibly clearer probe channels for the UGD. Furthermore, the detection of forward quarkonia (recently studied in the context of high-energy factorization in more central directions of rapidity [43,44,45,46,47,48]) in the low\(-p_T\) region will help us (i) to shed light on the quarkonium production mechanisms and (ii) to investigate kinematic regions at the frontier between the high-energy and the TMD dynamics.

The analysis presented in this work is at leading order and an obvious improvement would be its extension to the next-to-leading order. This calls for the setup of a theoretical scheme for the convolution of NLO IFs and UGD, which is not trivial, except for UGD models based on BFKL.