QCD evolution based evidence for the onset of gluon saturation in exclusive photo-production of vector mesons

TL;DR

This paper provides evidence for gluon saturation in protons by analyzing vector meson photo-production data from HERA and LHC, comparing non-linear and linear gluon evolution models.

Contribution

It demonstrates that non-linear gluon evolution models better describe experimental data, indicating the onset of gluon saturation in high-energy proton interactions.

Findings

Linear NLO evolution fails unless unnaturally large corrections are used

Non-linear BK evolution accurately fits energy dependence of data

Evidence supports high gluon densities and gluon saturation in protons

Abstract

We investigate photo-production of vector mesons J/Psi and Upsilon measured both at HERA and LHC, using 2 particular fits of inclusive unintegrated gluon distributions, based on non-linear Balitsky-Kovchegov evolution (Kutak-Sapeta gluon; KS) and next-to-leading order Balitsky-Fadin-Kuraev-Lipatov evolution (Hentschinski-Sabio Vera-Salas gluon; HSS). We find that linear next-to-leading order evolution can only describe production at highest energies, if perturbative corrections are increased to unnaturally large values; rendering this corrections to a perturbative size, the growth with energy is too strong and the description fails. At the same time, the KS gluon, which we explore both with and without non-linear corrections, requires the latter to achieve an accurate description of the energy dependence of data. We interpret this observation as a clear signal for the presence of high gluon densities in the proton, characteristic for the onset of gluon saturation.

Tables

Table 1. Table 1: Parameters of the boosted Gaussian vector meson wave functions for J / ψ 𝐽 𝜓 J/\psi and Υ Υ \Upsilon [ 16 , 18 ] .

Meson	$m_f/\text{GeV}$	${𝒩}_T$	${ℛ}^2$/${\text{GeV}}^{-2}$	$M_V$/GeV
$J/\psi$	$m_c=1.4$	$0.596$	$2.45$	$3.097$
${\rm Y}$	$m_b=4.2$	$0.481$	$0.57$	$9.460$

Equations

$$\frac{d\sigma }{dt}\left(\gamma p\to Vp\right)$$

$$B_D(W)$$

$${\sigma }^{\gamma p\to Vp}(W^2)$$

$$\frac{\Re \text{e}\mathcal{A}(W^2,t)}{\Im \text{m}\mathcal{A}(W^2,t)}$$

$$\Im \text{m}{\mathcal{A}}_T^{\gamma p\to Vp}(W,t=0)$$

$${\left({\Psi }_V^{\ast }\Psi \right)}_T(r,z)$$

$${\varphi }_{T,L}^{1s}(r,z)$$

$$2\int d^2\mathbf{b}\mathcal{N}\left(x,r,b\right)$$

$${\sigma }_{q\bar{q}}(x,r)$$

$${\sigma }_{q\bar{q}}^{(\text{HSS})}(x,r)$$

$${\hat{\sigma }}_{q\bar{q}}^{(\text{dom})}(x,r,M^2)$$

$${\hat{\sigma }}_{q\bar{q}}^{(\text{corr.})}(x,r,M^2)$$

$$\displaystyle\times\frac{\bar{\alpha}_{s}^{2}\beta_{0}\chi_{0}\left(\gamma\right)}{8N_{c}}\log{\left(\frac{1}{x}\right)}\Bigg{[}-\psi\left(\delta-\gamma\right)+\log\frac{M^{2}r^{2}}{4}-\frac{1}{1-\gamma}-\psi(2-\gamma)-\psi(\gamma)\Bigg{]}\;,$$

$$f(\gamma ,Q_0,\delta ,r)$$

$$\chi \left(\gamma ,M^2\right)$$

$$W(r)$$

Full text

1 Introduction

Energies available at the LHC allow for a detailed study of dynamical effects of Quantum Chromodynamics (QCD). A prominent example is the Heavy Ion program which focuses on the study of the Quark Gluon Plasma (QGP). While many of the QGP properties are understood, the initial state leading to its formation poses still open questions [1, 2]. A closely related question is the formation of an over occupied system of gluons, which eventually leads to saturation of gluon densities [3]; finding convincing phenomenological evidence for gluon saturation is still one of the open problems of QCD. On microscopic level, gluon saturation is expected to arise as a consequence of recombination of so-called wee gluons. The net effect of such recombination is to slow down the growth of gluon number density with energy, commonly referred to as gluon saturation. The evolution from the low to large gluon densities is described by a set of nonlinear evolution equations, known as Balitsky-Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner evolution [4, 5, 6]; its frequent used mean field version is given by the Balitsky Kovchegov (BK) [4, 7] evolution equation. First hints for gluon saturation could be already found in data collected from HERA experiments [8], (for a review see: [9] and references therein). More recent studies focus on resolving hadronic final states at RHIC and LHC to explore the transverse momentum distribution in the potentially saturated proton and nucleus, e.g. hadron production [10, 13, 11, 14, 12] or azimuthal de-correlations of dijets [15]. Such studies allow to study signals of saturation linked indirectly to the energy dependence of gluon distributions.

In contrast to these studies, exclusive production of vector mesons provides a direct opportunity to explore the energy dependence of the cross section [16, 17, 18]. Since for this observable large amount of data on production of $J/\Psi,\Upsilon(1S)$ has been collected over wide range of energies, it allows to potentially observe the predicted slow down of the growth with energy of cross-sections, which is one of the core predictions associated with the presence of high and saturated gluon densities. Studies in the literature for this process, which take into account effects due to gluon saturation, exist both on the level of dipole models [16, 17, 18, 19, 20] and complete solutions to non-linear BK equation [21, 22]. The description of the cross sections is in general satisfactory. Here we would like to re-investigate the problem from a different angle: While the ability of non-linear evolution equations to describe collider data is a necessary requirement to establish the presence of gluon saturation, it is not sufficient: albeit a successful description of data is provided, one might still remain in the so-called dilute regime, where gluon densities are perturbative. That this could be actually the case, is at least suggested by the successful description of the same data set by frameworks which rely on collinear factorization [23, 24, 25, 26] and – maybe even more striking – linear NLO BFKL evolution [27]. The latter relies on a fit of an initial (low energy) transverse momentum distribution to combined HERA data by Hentschinski-Salas-Sabio Vera (HSS)[28], and has been recently explored in a number of phenomenological studies, see e.g. [29].

In the following we will argue that even though a description of exclusive photo-production data by linear QCD evolution can be achieved for the entire range of available center-of-mass energies, such a description requires unnaturally large higher order corrections, which in some region of phase space are larger than the dominant leading order term. To assess the importance of non-linear terms in low $x$ evolution equations, we use a particular solution to BK-evolution, with initial conditions fitted to combined HERA data by Kutak-Sapeta (KS) [30]. While the HSS-gluon relies on NLO-BKFL evolution and KS-evolution on LO-BK evolution, both evolution schemes supplement the original low $x$ evolution with a supplementary resummation of collinear logarithms [31]. Our strategy is then as follows: While exclusive photo-production of the $\Upsilon$ serves as the necessary cross-check in the perturbative, dilute region (provided by the relatively large bottom quark mass), saturation effects will be searched for in photo-production of the $J/\Psi$ mesons (where the hard scale is provided by the charm quark mass). In particular we will demonstrate that the seemingly flawless description of $J/\Psi$-data by NLO BFKL evolution is only possible due to the presence of a very large perturbative correction, which slows down the growth with energy. Choosing on the other hand a hard scale which renders this perturbative corrections small, the data are no longer described. The KS gluon (which is subject to non-linear QCD evolution) is on the other hand able to describe the energy dependence of $J/\Psi$-data. Turning off the non-linear effects in the evolution of the KS-gluon, we find that a description of the energy dependence is no longer possible. We are convinced that this observation provides a very strong phenomenological evidence for the presence of saturation effects in the high $W$ region of $J/\Psi$ photo-production observed at the LHC.

The outline of this letter is as follows: In Sec. 2 we provide details of our theoretical description, Sec. 3 is dedicated to a discussion of the large perturbative corrections of the NLO BFKL gluon in the large $W$ region while in Sec. 4 we present our conclusions.

2 Energy dependence of the photo-production cross-section

We study the process 111Besides HERA data we also use the LHC p-p and Pb-p data where highly boosted p and Pb respectively become a source of photons leading to Ultra Peripheral Collisions $\gamma(q)+p(p)\to V(q^{\prime})+p(p^{\prime})$ where $V=J/\Psi,\Upsilon(1S)$ and $\gamma$ denotes a quasi-real photon with virtuality $Q\to 0$; $W^{2}=(q+p)^{2}$ is the squared center-of-mass energy of the $\gamma(q)+p(p)$ collision. The $x$ value probed in such a collision is obtained as $x\simeq M_{V}^{2}/W^{2}$ with $M_{V}$ the mass of the vector meson. With the momentum transfer $t=(q-q^{\prime})^{2}$, the differential cross-section for the exclusive production of a vector meson can be written in the following form

[TABLE]

where $\mathcal{A}(W^{2},t)$ denotes the scattering amplitude for the reaction $\gamma p\to Vp$ for color singlet exchange in the $t$-channel, with an overall factor $W^{2}$ already extracted. For a more detailed discussion of the kinematics we refer to [27].

2.1 The theoretical setup of our study

In the following we determine the total photo-production cross-section, based on an inclusive gluon distribution. This is possible following a two step procedure, frequently employed in the literature: First one determines the differential cross-section at zero momentum transfer $t=0$ (which can be expressed in terms of the inclusive gluon distribution); in a second step the $t$-dependence is modeled which then allows us to relate the differential cross-section at $t=0$ to the integrated cross-section. Here we follow the prescription given in [23, 24], where an exponential drop-off with $|t|$, $\sigma\sim\exp\left[-|t|B_{D}(W)\right]$ is used with an energy dependent $t$ slope parameter $B_{D}$, as motivated by Regge theory,

[TABLE]

Following [23, 24], we use for the numerical values $\alpha^{\prime}=0.06$ GeV*-2*, $W_{0}=90$ GeV and $b_{0}^{J/\Psi}=4.9$ GeV*-2* in the case of the $J/\Psi$, while $b_{0}^{\Upsilon}=4.63$ GeV*-2* for $\Upsilon$ production. The total cross-section for vector meson production is therefore obtained as

[TABLE]

The uncertainty introduced by the modeling of the $t$-dependence mainly affects the overall normalization of the cross-section with a mild logarithmic dependence on the energy. To determine the scattering amplitude, we first note that the dominant contribution is provided by its imaginary part. Corrections due to the real part of the scattering amplitude can be estimated using dispersion relations, in particular

[TABLE]

As noted in [27, 32], the dependence of the slope parameter $\lambda$ on energy $W$ provides a sizable correction to the $W$ dependence of the complete cross-section. We therefore do not assume $\lambda=$ const., but instead determine the slope $\lambda$ directly from the $W$-dependent imaginary part of the scattering amplitude. The latter is obtained from [20, 19]

[TABLE]

where $\mathcal{N}\left(x,r,b\right)$ is the dipole amplitude and $T$ denotes transverse polarization of the quasi-real photon. Here

[TABLE]

with $\epsilon^{2}=m_{f}^{2}$ for real photons. Furthermore $r=\sqrt{{\bm{r}}^{2}}$, while $f=c,b$ denotes the flavor of the heavy quark and $\hat{e}_{f}=2/3$, $-1/3$. For the scalar parts of the wave functions $\phi_{T,L}(r,z)$, we employ the boosted Gaussian wave-functions with the Brodsky-Huang-Lepage prescription [33]. For the ground state vector meson ($1s$) the scalar function $\phi_{T}(r,z)$, has the following general form [20, 34],

[TABLE]

The free parameters $N_{T}$ and $\mathcal{R}_{1s}$ of this model have been determined in various studies from the wave function normalization and the decay width of the vector mesons. In the following we use the values found in [16] ( $J/\Psi$) and [18] ( $\Upsilon$). The parameters are summarized in Tab. 1.

In the forward limit $t=0$, one further has,

[TABLE]

where $\sigma_{q\bar{q}}$ denotes the inclusive dipole cross-section which is related to the unintegrated gluon density ${\cal F}$ through [35]

[TABLE]

In [27] this expression has been used to calculate the BFKL impact factor in transverse Mellin space from the light-front wave function overlap Eq. (6). In the following study we chose a slightly different route and calculate the dipole cross-section directly from the regarding KS and HSS unintegrated gluon densities. These gluon densities have been obtained in the following way:

• •

The KS gluon is obtained as a solution of the momentum space version of BK equation with modifications according to the Kwieciński-Martin-Satsto (KMS) prescription [31] which includes a kinematical constraint to impose energy momentum conservation, complete DGLAP splitting function and contribution of quarks as well as the 1-loop QCD running coupling. As a consequence this evolution equation reduces in the collinear limit to leading order DGLAP evolution. The intial conditions of the KS gluon distribution have been determined using a fit [30] combined HERA data for the proton structure function $F_{2}$ [63] in the region $x<10^{-2}$.

• •

Below we will further consider the linear version of the KS gluon which is obtained as a solution of momentum space version of the leading order BFKL equation with modifications according to the KMS prescription described above. The initial conditions of the linear equation have been obtained from a fit to HERA data, similar to the non-linear case, but with photon virtualities restricted to the region $Q^{2}>4.5$ GeV2.

• •

The HSS gluon is obtained as a solution to the NLO BFKL equation in transverse Mellin space applying both a resummation of collinear logarithms within the NLO BFKL kernel and a resummation of large running coupling corrections (optimal renormalization scale setting). The initial conditions have been fitted [28] to combined HERA data for proton structure function $F_{2}$ in the region $x<10^{-2}$ and photon virtualiyt $Q^{2}\geq 1.2$ GeV2.

For a detailed discussion of the framework underlying both gluon distributions we refer to the literature: [30, 31] (KS) and [28, 36] (HSS).

2.2 Numerical results using standard implementations

The main uncertainty left is the scale at which the strong coupling constant $\alpha_{s}$ is to be evaluated in Eq. (9). In the case of the $J/\Psi$, which is characterized by a relatively small hard scale $m_{c}\simeq 1.4$ GeV and therefore large value of the strong coupling constant $\alpha_{s}\simeq 0.31$, this leads to a sizable ambiguity in the normalization of the total cross-section, since the latter depends through Eq. (1) quadratically on $\alpha_{s}$. Using similar conventions as used in original fits of the KS and HSS gluon distributions, we fix this scale to a typical hard scale of the process. For the KS gluon we chose both for the photo-production of $\Upsilon$ and $J/\Psi$ vector mesons, the mass of the respective heavy quark as the hard scale; for the HS-gluon it was found in [27] that a scale related to the size of the $J/\Psi$ wave function is more suitable in the case of $J/\Psi$-production, $M_{J/\Psi}^{\text{HS}}=8/\mathcal{R}_{J/\Psi}^{2}=3.27$ GeV while we use the bottom quark mass for $\Upsilon$-production. The results of our study for the fixed scale case can be found in Fig. 1, where continuous, black lines correspond to the KS-gluon and dotted, green lines to the HS-gluon at fixed scales; the dashed green lines, corresponding to a special scale setting of the HSS gluon, and the linear KS-gluon will be discussed in the forthcoming section. We observe that both the KS-gluon distribution and the HSS-gluon distribution provide an excellent description of the energy dependence of the data. While the KS-gluon requires in the current study a K-factor of $1.4-1.5$, we note that the size of such a correction strongly depends on the precise scale choice of the strong coupling constant and the precise parametrization used for the $t$-slope parameter $B_{D}$, Eq. (2). Indeed, using the parametrization of the $B_{D}$ parameter suggested in [18], would bring the K-factor of the KS-gluon close to one in the case of $\Upsilon$-photo-production. We further note that our study does not include a frequently employed phenomenological corrective factor which can be determined through relating the inclusive collinear gluon distribution to generalized parton distribution through a Shuvaev transform. While the precise numerical value of this factor depends on the energy dependent slope parameter $\lambda$, it generally provides a correction which is of the order of the K-factor found for the KS-gluon. Our general conclusion is that the theoretical accuracy of the observable is currently not sufficient to fix unambiguously the normalization. At the same time, the energy dependence is described in an excellent way, both by KS-evolution (non-linear BK evolution combined with DGLAP corrections and kinematical constraint and collinear resummation) and the HSS-gluon (NLO BFKL evolution with collinear resummation).

3 The need for non-linear low $x$ evolution

At first sight it appears that both non-linear $W$ evolution (KS) and linear $W$ evolution (HSS) describe the $W$-dependence of data; one might be therefore lead to conclude that non-linear effects, associated with the presence of high and saturated gluon densities, are absent and the observables merely probes the linear, dilute regime. In the following we argue that such a conclusion is pre-mature. Indeed there are strong hints which suggest that we are at least in the transition region towards high and saturated gluon densities.

To fully access this question, we first recapitulate which possible impact large gluon densities could have on the observable. First of all, the presence of high density effects cannot be seen directly at the level of the observable. The scattering amplitude Eq. (5) depends only on the dipole amplitude, which itself can be expressed as the correlator of two Wilson lines which resum the gluonic field of the proton, see e.g. [50]. Even though the dipole amplitude resums the interaction of the $q\bar{q}$-dipole with an infinite number of gluons, the gluons couple to the $q\bar{q}$-dipole like a single gluon; the “reggeize” in the language of [51] and therefore appear like a single gluon. At the level of our phenomenological study, this property reveals itself through Eq. (9), which relates the dipole cross-section to the unintegrated gluon density. To make multiple re-scattering of partons on the target field visible, it would be necessary to resolve the hadronic final state of the dissociated photon, see e.g. [52, 53]. This not the case for photo-production of vector mesons. The only place where one could expect a signal for the presence of saturation effects is therefore the $x$-dependence of the underlying gluon distributions. As an immediate consequence, any framework which is based on a direct fit of the $x$-dependence at the $J/\Psi$ scale (such as collinear parton distribution functions) does not exclude presence of saturation effects; it merely demonstrates the ability to fit the resulting $x$-dependence of the underlying gluon distribution. While this initial $x$-distribution can be evolved through DGLAP evolution to events with higher hard scales, such events are generally characterized by larger values of $x$ ($x_{\Upsilon}>2.28\cdot 10^{-5}$ vs. $x_{J/\Psi}>2.99\cdot 10^{-6}$ in the current case). Taking further into account that DGLAP evolution is known to shift large $x$ input to lower $x$, it is therefore safe to say that the mere ability of DGLAP fits to accommodate low $x$ $J/\Psi$ photo-production data, does not exclude the potential presence of sizable non-linear effects for the data points at highest $W$-values.

Instead of DGLAP evolution, a suitable benchmark to establish presence/absence of gluon saturation is provided by linear NLO BFKL evolution, such as the HSS gluon. While the HSS gluon provides a very good description of both $\Upsilon$ and $J/\Psi$ photo-production data, the following observation can be made: Recalling the particularly solution of NLO BFKL evolution used for the HSS-fit, one finds at the at level of the dipole cross-section two terms

[TABLE]

where

[TABLE]

and

[TABLE]

is a function which collects both factors resulting from the proton impact factor and the transformation of the unintegrated gluon density to the dipole cross-section, see [28, 27] for details. The parameters $Q_{0}=0.28$ GeV, $\mathcal{C}=2.29$ and $\delta=6.5$ have been determined from the HERA data fit. Furthermore $\bar{\alpha}_{s}=\alpha_{s}N_{c}/\pi$ with $N_{c}$ the number of colors, and $\chi(\gamma,M^{2})$ is the next-to-leading logarithmic (NLL) BFKL kernel after collinear improvements; in addition large terms proportional to the first coefficient of the QCD beta function, $\beta_{0}=11N_{c}/3-2n_{f}/3$ have been resumed through employing a Brodsky-Lepage-Mackenzie (BLM) optimal scale setting scheme [54]. The NLL kernel with collinear improvements reads

[TABLE]

where $\chi_{i}$, $i=0,1$ denotes the LO and NLO BFKL eigenvalue and $\chi_{\text{RG}}$ resums (anti-)collinear poles to all orders; for details about the individual kernels see [28, 27]. The scale $M$ is a characteristic hard scale of the process. The second contribution $\hat{\sigma}^{\text{corr.}}$ is proportional to $\beta_{0}$ and acts in $\gamma$-space as a differential operator on the impact factors of external particles. These terms do not exponentiate and have been therefore treated in [28] as a perturbative correction to the BFKL Green’s function. Even though $\hat{\sigma}^{\text{corr.}}$ is suppressed by a factor of $\alpha_{s}^{2}$, enhancement by $\ln(1/x)$ will eventually compensate for the smallness of the strong coupling constant and invalidate the perturbative expansion.

The behavior of the HSS-dipole cross-section is studied in Fig. 2. To identify the relevant region in dipole size $r$ for $J/\Psi$ photo-production we further define

[TABLE]

as the $z$-integrated wave function overlap. Working with a fixed hard scale $M^{2}=3.27$ GeV2 we find in Fig. 2, that the perturbative expansion is well under control for a typical HERA $x$ value of $x=3.55\cdot 10^{-4}$ (Fig. 2.a). Turning however to the lowest $x$ values probed at the LHC of $x=2.81\cdot 10^{-6}$ (Fig. 2.b) we observe that the correction term is generally large; for certain $r$ values, which are further enhanced through the $r$-dependence of $W(r)$, they even super-seed the dominant term, resulting into a negative dipole cross-section. While the dipole cross-section is not an observable, this clearly indicates a breakdown of the perturbative expansion for dipole sizes where the integrated wave function overlap $W(r)$ has its maximum value. The problem of unnaturally large higher order corrections can be fixed by choosing a hard scale related to the transverse size of the dipole. Following the scale setting used in fits [55] of the IP-sat model [56], we may therefore choose $M^{2}=\frac{4}{r^{2}}+\mu_{0}^{2}$ with $\mu_{0}^{2}=1.51~{}\text{GeV}^{2}$. With this scale, we find that the perturbative expansion indeed stabilizes: both for the HERA (Fig. 2.c) and LHC $x$-values (Fig. 2.d) the perturbative term is well under control. Turning with this choice for the hard scale however to data, we find that this scale setting (green dashed line in Fig. 1) describes very well the energy dependence of $\Upsilon$-photo-production as well as $J/\Psi$ photo-production in the HERA region $W<300$ GeV, but fails to describe $J/\Psi$ production at the LHC ($W>300$ GeV). The resulting growth with energy is too strong and the data are no longer described (Fig. 1, top). We therefore conclude that NLO BFKL evolution can only describe data in the region $W>300$ GeV if one accepts very large perturbative corrections, which super-seed for certain dipole sizes the dominant term and which slow down the growth of the cross-section. If the size of these perturbative corrections is reduced using a suitable hard scale, the growth of the HSS-gluon is too strong and cannot be accomodated by data.

The KS-gluon, which is subject to LO-BK evolution with collinear resummation provides a very good description of $J/\Psi$ data in the region $W>300$ GeV. To answer the question whether this description relies on the presence of non-linear terms in the evolution equation, we compare in addition to a linearized KS gluon (dashed black line in Fig. 1). We find that the growth in $W$ is in that case even stronger than for the HSS gluon with dipole size scale. The ability of the KS gluon to describe data is therefore strongly connected to non-linear evolution effects.

4 Conclusion and Discussion

We conclude that there are strong hints for the presence of the saturation effects in exclusive photo-production of $J/\Psi$ at small $x$. While both linear and non-linear low $x$ QCD evolution can describe the data, the former requires the presence of unnaturally large perturbative corrections. Rendering these corrections small, $\Upsilon$-data, characterized by a hard perturbative scale, are still well described while the growth with energy is too strong for $J/\Psi$ data in the LHC region with $W>300$ GeV. The successful description of data by the KS-gluon directly relies on including non-linear terms in the evolution; with those terms being absent, the description breaks down.

To observed slow-down of the growth with energy is one of the core predictions of gluon saturation. We therefore are convinced that the current study provides substantial evidence for the presence of the saturation effects. We wish to note that a related observation has already been made in [16] on the level of dipole models. The current study substantiates this observation through employing dipole cross-sections which are directly subject to linear and non-linear QCD evolution.

Nevertheless it must be noted that the current study is not without deficits: the description is based on LO wave function and LO BK evolution (although supplemented with collinear resummations). To establish the observation made in this letter it is therefore necessary to search for different observables which probe the low $x$ gluon in a similar kinematic regime and to increase further the theoretical accuracy of the underlying framework. Steps to address the latter are currently undertaken in [57, 58] (evolution equations) and [59, 60, 61, 62] (determination of higher order corrections).

5 Acknowledgments

This work is partly supported by the Polish National Science Centre grant no. DEC-2017/27/B/ST2/01985 and by COST Action CA16201 ”Unraveling new physics at the LHC through the precision frontier”.