$J/\psi$ production at NLO with a scale-dependent color-evaporation model

Abstract

Nearly ten years ago, Kang, Ma, Qiu, and Sterman derived an evolution equation for a $Q\bar{Q}$ pair fragmenting into a quarkonium. In this study we explore the consequence of this evolution for the color-evaporation model, focusing on $J/\psi$ transverse-momentum ($p_{t}$) distributions in proton-proton collisions. We show that, as expected, it softens the spectrum obtained by fixed-order calculations. While next-to-leading-order calculations strongly overestimate data at large $p_{t}$, ours, including the (approximate) $Q\bar{Q}$ evolution and next-to-leading-order cross sections computed with Madgraph, are in good agreement with experiments. Since our study with the color-evaporation model shows a significant effect of the $Q\bar{Q}$ evolution at large $p_{t}$, a new determination of long-distance-matrix elements of nonrelativistic QCD could be necessary. To describe data at small and intermediate $p_{t}$, we use the $k_{t}$-factorization approach, and we argue that quarkonia data could help constrain this formalism.

1 Introduction

Quarkonia production is important in Quantum Chromodynamics (QCD) because it provides insights into the fundamental dynamics of QCD, the structure of the QCD vacuum, the effects of heavy quark masses, the properties of the quark-gluon plasma, and allows for experimental tests of theoretical predictions. Models for quarkonia production, generally based on a factorization hypothesis, differ by their treatment of hadronization. For instance, the nonrelativistic QCD (NRQCD) [1] expression for the quarkonium cross section reads

$$\sigma^{H}=\sum_{[Q\bar{Q}(n)]}\hat{\sigma}_{[Q\bar{Q}(n)]}(\Lambda)\langle 0|\mathcal{O}_{[Q\bar{Q}(n)]}^{H}(\Lambda)|0\rangle,$$ (1)

where $\hat{\sigma}_{[Q\bar{Q}(n)]}$ are the process-dependent short-distance coefficients to produce a heavy-quark pair in a quantum state $n$ (color, spin, and orbital angular momentum states). The second factor in the right side of Eq. (1) is the long-distance-matrix element (LDME) describing the non-perturbative hadronization of the $Q\bar{Q}(n)$ pair into the quarkonium $H$. Both factors depend on the ultraviolet cutoff $\Lambda\sim\mathcal{O}(m_{Q})$, with $m_{Q}$ the heavy-quark mass.

On the opposite, the CEM uses the same weight for all $Q\bar{Q}$ states and has only one non-perturbative parameter, the hadronization factor $F_{H}$. It is then natural to think that the CEM is not fully consistent with QCD, as discussed in Ref. [2]. However, its simplicity and reasonably good description of numerous observables make it a useful tool for the study of quarkonia production.

An issue of the CEM is the too-hard spectrum due to next-leading order (NLO) contributions scaling as $1/p_{t}^{4}$ [3], where $p_{t}$ is the $J/\psi$ transverse momentum. An example of such contribution is shown in Fig. 0(a), along with the result of fixed-order NLO calculations obtained with madgraph5_aMC@NLO (MG5) [4] in Fig. 0(b).

((a))

((b))

Figure 1: (a) Example of NLO contribution scaling as $1/p_{t}^{4}$. (b) Comparison of fixed-order calculations (red triangles) with ATLAS data (blue dots) for the transverse-momentum distribution of prompt $J/\psi$ [5]. For this result, we used madgraph5_aMC@NLO at NLO using a 3-flavor-number scheme and a charm mass $m_{c}=1.3$ GeV.

In this study, we demonstrate that the evolution of the heavy-quark pair, presented in Sec. 2, brings CEM calculations in agreement with data for the $p_{t}$-distribution of the $J/\psi$ particle. This evolution, related to the scale-dependence of fragmentation functions, is a central prediction of perturbative QCD.

Our first goal is not to support the CEM, but to underline the importance of the heavy-quark pair evolution. Of course, having the accuracy of this model improved is also of interest. The present study suggests that evolution could significantly affect the determination of NRQCD’s LMDEs. We could reasonably hope that a new extraction of these non-perturbative functions, including the evolution, will alleviate the tensions between different LDME sets and improve the NRQCD description of experimental data.

To describe data on the whole $p_{t}$ range, we use a transverse-momentum dependent formalism (see Sec. 3). Our second goal is to use quarkonium data to put constraints on transverse-momentum-dependent PDFs. Our main results are presented in Sec. 4, and we discuss the universality of the non-perturbative parameter, $F_{J/\psi}$, in Sec. 5. Finally, we show similar results obtained with the event generator EPOS4 in Sec. 6. The common features of EPOS4 with our calculations are the timelike cascade, giving the (approximate) scale evolution of the $Q\bar{Q}$ pair, and its hadronization based on the CEM. The similar results obtained with two different formalisms demonstrate the robustness of our main conclusion.

2 Evolution of the heavy-quark pair

Two groups published a series of papers discussing the scale evolution of quarkonia fragmentation functions. One of them worked with perturbative QCD [6, 7, 8, 9], while the other used the soft-collinear-effective theory [10, 11]. Contributions to quarkonia production are first separated into leading power, scaling as $1/p_{t}^{4}$, and next-to-leading power (NLP), scaling as $1/p_{t}^{6}$, and then factorized. It is the factorization of NLP contribution that involves the convolution of a short-distance partonic cross section for the production of a $Q\bar{Q}$ pair with the modified double-parton fragmentation function of a $Q\bar{Q}$ pair into a quarkonium $H$. This fragmentation function obeys the following evolution equation

$$\frac{\partial}{\partial\ln\mu^{2}}D_{[Q\bar{Q}(\kappa)]\rightarrow H}\left(z,\mu^{2}\right)=\frac{\alpha_{s}(\mu)}{2\pi}\sum_{n}\int_{z}^{1}\frac{dz^{\prime}}{z^{\prime}}P_{[Q\bar{Q}(n)]\rightarrow[Q\bar{Q}(\kappa)]}\left(\frac{z}{z^{\prime}}\right)\\ \times D_{[Q\bar{Q}(n)]\rightarrow H}\left(z^{\prime},\mu^{2}\right),$$ (2)

where $\kappa$ and $n$ denote the possible quantum numbers of the $Q\bar{Q}$ pair, and $P_{[Q\bar{Q}(n)]\rightarrow[Q\bar{Q}(\kappa)]}$ are the modified evolution kernels, see [8] for more details.

For the standard choice $\mu=p_{t}$, we see that at intermediate and large transverse momentum the color-octet (CO) and color-singlet (CS) states mix under Eq. (2). The different models for quarkonia production should be applied at $\mu_{0}\sim m_{Q}$, and the CEM uses a $\delta(z-1)$ as $z$ distribution [12, 13]. Naively, we expect the evolution to reduce the differences between the different production mechanisms. Indeed, a $Q\bar{Q}$ pair in a CS state at $\mu_{0}$ could have been produced in a CO state at $\mu=p_{t}$. The mixing of quantum states induced by the evolution step $Q\bar{Q}(n)\rightarrow Q\bar{Q}(\kappa)$ implies a mixing of the LO behavior $1/p_{t}^{a}$ ($a>0$) associated to these states. The naive (LO) $p_{t}$ dependence of the produced $Q\bar{Q}(n)$ is further modified by real emissions, partly responsible for the evolution equation, and resulting in energy loss. It is in fact this effect which will bring the CEM calculations in agreement with data.

In our study, we did not work with the solution of Eq. (2), but used the PYTHIA6 [14] timelike cascade to evolve the heavy-quark pair. A timelike cascade is an implementation of the DGLAP equation [15, 16, 17] for timelike partons, and a central tool for all realistic event generators. An important simplification made by these algorithms is that timelike partons do not interact. They simply split into two partons until they reach the lower cutoff of the cascade. By evolving the heavy quark and antiquark independently, we neglect several Feynman diagrams in Sec. IV C. of Ref. [8]. For instance, (virtual) diagrams with a gluonic line between the quark and antiquark in the amplitude (and nothing in the conjugate amplitude) are not included. Interference terms, where the gluon is emitted by the quark in the amplitude and by the antiquark in the conjugate amplitude are also not strictly included.

However, virtual corrections are sometimes negative, cancellations can occur and the importance of the neglected Feynman diagrams is unclear. From phenomenological knowledge, one may argue that their net contribution is not significant. Indeed, the discussion for the $Q\bar{Q}$ pair is also true for any hard parton producing a jet; we neglect the Feynman diagrams associated with the interactions between timelike particles forming the jet. Since event generators and timelike cascades provide a good description of data, including exclusive observables, it seems to be a reasonable approximation. From this point of view, the heavy quark and antiquark have nothing special. In the real world, after a few collinear splittings, they will be (closely) surrounded by other partons and interact with them. But the heavy quarks are evolved as any other parton, by neglecting these interactions, and, in the end, will hadronize either into D mesons or charmonium.

Thanks to the simplicity of the CEM, the absence of evolution is easily observable, since, as demonstrated in Sec. 4, it is responsible for the overestimation of data by NLO calculations. This is another reason why we choose to work with the CEM. In the case of NRQCD, taking into account the evolution would certainly require a new determination of the LDMEs. The consequences of this could be modest for inclusive observables, and more visible for exclusive observables.

3 Heavy-quark pair production, evolution and hadronization

In order to describe $J/\psi$ production on the whole $p_{t}$ range, we use a transverse-momentum dependent formalism. A possible choice would be the TMD factorization, discussed for quarkonia production in proton-proton (pp) collisions using the soft-collinear and NRQCD effective field theories [19]. Another possibility is the $k_{t}$ factorization approach [20, 21, 22, 23, 24]. This formalism has been extensively used for the study of $J/\psi$ production in pp collisions [25, 26, 27, 28, 29, 30, 31, 32, 33]. In our case, we used a hybrid formalism, already used for Drell-Yan [34], where the off-shell cross section is replaced by on-shell cross section. For proton-proton collisions, the differential cross section reads

$$\frac{d\sigma(pp\to J/\psi+X)}{dx_{1}dx_{2}d^{2}p_{t}}(s,x_{1},x_{2},p_{t})=\sum_{a,b}\int^{k_{t,\text{max}}^{2}\sim s}_{0}d^{2}k_{1t}d^{2}k_{2t}F_{a/p}(x_{1},k_{1t};\mu)\\ \times F_{b/p}(x_{2},k_{2t};\mu)d\hat{\sigma}_{ab\to Q\bar{Q}}(x_{1}x_{2}s,q_{t};\mu)\otimes d_{Q\bar{Q}}(\mu,\mu_{0})\otimes D^{\text{CEM}}_{Q\bar{Q}\to J/\psi}(\mu_{0}),$$ (3)

where $\bm{q}_{t}=\bm{p}_{t}-\bm{k}_{1t}-\bm{k}_{2t}$, with $\bm{k}_{1t}$ and $\bm{k}_{2t}$ the transverse momentum of the two incoming partons. The variables $x_{1,2}$ are the longitudinal momentum fractions carried by the partons, $\mu$ is the factorization scale, and $s$ is the usual Mandelstam variable for the proton-proton system. The functions $F_{i/p}$ are the unintegrated PDFs (UPDFs) for a parton $i$ inside a proton. The transverse-momentum dependence of the initial partons is managed with CASCADE3, and we used the parton-branching (PB) UPDFs, set 2 [35], extracted from HERA data. They obey an evolution equation [36, 37], which can be combined with finite-order matrix elements either by matching [38, 34], or merging [39, 40]. The transverse-momentum distribution of the $J/\psi$ particle has already been studied with the event generator CASCADE using more traditional $k_{t}$-factorization calculations, e.g., with LO off-shell matrix elements and the color-singlet model, see, for instance, Ref. [26].

In Eq. (3), the differential cross section $\hat{\sigma}_{ab\to Q\bar{Q}}$ for the production of an unpolarized heavy-quark pair is computed at NLO with MG5, option -p, which includes all necessary subtraction terms. Option -p means the parton showers are not performed within MG5, but with other tools, in our case, CASCADE3 [41, 42] and pythia 6. We worked in a 3-flavor scheme with $m_{c}=1.3$ GeV, and used the CT14nlo_NF3 [43] PDFs. This choice is imposed by the fact that, in MG5, the variable-flavor-number scheme (VFNS) works with $m_{c}=0$ (i.e., it is a Zero-Mass VFNS). In [44], the authors argue that, in practice, it is generally better to use VFNS PDFs, even with cross sections computed in a 3-flavor scheme. We checked that changing the CT14nlo_NF3 for the CT14nlo gives similar results. With the different PDFs sets tested in our study, we observed a negligible or small variation of the final result compared to scale variation. The function $d_{Q\bar{Q}}(\mu,\mu_{0})$ implements the evolution of the heavy-quark pair, with a convolution on the variable $z$ representing the momentum fraction of the initial $Q\bar{Q}$ state carried by the final quarkonium. To be clear, we do not have an analytical expression for the evolution function $d_{Q\bar{Q}}$. All our calculations are performed by event generators using distribution probabilities, with the three steps being: 1. fixed-order production with MG5, 2. evolution with a timelike cascade, and 3. hadronization with the CEM. The latter is applied at $\mu_{0}\sim m_{J/\psi}$ and sums all pairs with an invariant mass between $2m_{c}$ and $2m_{D^{0}}$ [12, 13]

$$\frac{d\sigma}{dp_{t}}=F_{J/\psi}\int_{2m_{c}}^{2m_{D^{0}}}\frac{d\sigma^{Q\bar{Q}}}{dm_{Q\bar{Q}}dp_{t}}dm_{Q\bar{Q}}.$$ (4)

We fixed the parameter $F_{J/\psi}$ to the value found in Ref. [3] at NLO:

$$F_{J/\psi}=0.014.$$ (5)

In that sense, our calculations are parameter free. As visible from Eq. (3), our calculations include only the direct contribution, and will be compared to prompt-$J/\psi$ data. Finally, note that we did not use the improved CEM [45]

$$\frac{d\sigma}{d^{3}p}=F_{J/\psi}\int_{m_{J/\psi}}^{2m_{D^{0}}}\frac{d\sigma^{Q\bar{Q}}}{dm_{Q\bar{Q}}dp^{3}_{Q\bar{Q}}}\delta^{3}\left(p-\frac{m_{J/\psi}}{m_{Q\bar{Q}}}p_{Q\bar{Q}}\right)dm_{Q\bar{Q}},$$ (6)

because the reduced phase space ($m_{J/\psi}>2m_{c}$) requires more statistics. Note that the difference between the CEM and ICEM is visible for $p_{t}\lesssim 15$ GeV [46]. Several studies using LO off-shell matrix elements and the CEM are also available [29, 31, 33].

Contrary to the factorization formalism, the CEM does not organize the calculation into LP and NLP contributions. NLO diagrams with a gluon splitting into a $Q\bar{Q}$ pair (see Fig. 0(a)) contribute to both LP and NLP [7, 47]. In the factorization formalism, the LP contribution is removed by a subtraction term, avoiding mass singularities and double counting with the LP contribution $\sigma_{f}\otimes D_{f\to H}$. The latter describes the short-distance production of a parton $f$ followed by its fragmentation into the quarkonium $H$. In our calculations with MG5, no double counting occur because the contribution $\sigma_{f}\otimes D_{f\to H}$ is not included (then, several LP contributions are missing). Moreover, we keep the mass of heavy quark in the cross section, so the gluon cannot get on-shell and does not produce mass singularities. Then, CEM calculations do not include the subtraction term and retain some LP contributions, which are responsible for the overestimation of the cross section by fixed-order calculations. The missing LP contributions have certainly an effect on the determination of $F_{J/\psi}$, but not on the shape of the distribution at intermediate and large transverse momentum.

In the subsequent section, we compare our calculations with high-energy data, with the expectation that the evolution should improve the results obtained by fixed-order calculations. Indeed, a consequence of the evolution is energy loss by the heavy-quark pair, shifting fixed-order calculations of Fig. 0(b) to the left. The cascade usually starts at $\mu_{c}\simeq\mu\simeq p_{t}$, and we chose $\mu_{c}=\mu$ by simplicity. To larger $p_{t}$ corresponds a longer evolution and a larger shift in $p_{t}$, making the theoretical prediction softer. From a more theoretical point of view, we will obtain a softer spectrum because the evolution modifies the FO behavior of $1/p_{t}^{4}$ to [8]

$$\frac{m_{Q}^{2}}{p_{t}^{4}\mu^{2}}.$$ (7)

4 $J/\psi$ production at small, intermediate and large transverse momentum

Figure 2: Comparison of the prompt-$J/\psi$ differential cross section measured by ATLAS [5] as a function of transverse momentum $p_{T}$ (blue dots) with our calculations Eqs. (3) and (4) (red triangles). The mg5+cascade3 error bars includes only the statistical uncertainty. The ratio theory/exp is displayed in the bottom panel.

We start with ATLAS data on prompt $J/\psi$ at 8 TeV [5], already compared to the CEM at NLO (fixed order) in Fig. 1, see also Fig. 33 of Ref. [3]. While fixed-order calculations fail, we observe in Fig. 2 a good agreement between our full calculations (mg5+cascade3) and data. The results obtained for the other rapidity ranges are displayed at the end of this manuscript. At large $p_{t}$, the main difference between Figs. 0(b) and 2 is the timelike cascade. As expected, including the $Q\bar{Q}$ evolution improved the result. It was less clear, however, that this sole effect would bring CEM calculations in agreement with data, this model being quite simple. But we should be careful to not overinterpret these results. Our calculations do not include the feed down contributions, and the timelike cascade is only an approximation of the true evolution equation. Still, improving these points will certainly not change our conclusion. Note that, additionally of the energy lost by the heavy-quark pair, another effect of the timelike cascade is to change the invariant mass, sometimes destroying quarkonium candidates with $2m_{c}<m_{Q\bar{Q}}<2m_{D^{0}}$. The opposite situation is also possible, promoting $Q\bar{Q}$ pairs produced at fixed order with the incorrect invariant mass to quarkonium candidates. An advantage of event generators is to automatically include such kinematical effects.

Figure 3: Differential cross section, Eq. (3) (without the last factor), at $\sqrt{s}=200$ GeV as a function of the $Q\bar{Q}$ transverse momentum. The red and blue lines correspond to the UPDFs PB set 1 and 2. The green line corresponds to the case with no initial $k_{t}$.

While taking into account the initial transverse momentum is not primordial at large $p_{t}$ (with our formalism), this effect is essential at small $p_{t}$. This is illustrated in Fig. 3, where we plot the result obtained at $\sqrt{s}=200$ GeV, with the initial transverse momentum obtained from the PB set1 and set2, and with no initial transverse momentum in green. While both PB sets give a similar result, the green line corresponding to no initial transverse momentum is more than one order of magnitud below. Since the LO contribution is exactly zero in the case of no initial $k_{t}$, this line corresponds to the pure NLO contribution. We observe that this contribution starts to dominate at $p_{t}\sim 10$ GeV.

In Figs. 4 and 5, we compare our calculations at 7 TeV with LHCb [48], ALICE [49], and CMS [50] data.

Figure 4: Comparison of the $J/\psi$ differential cross section measured by the LHCb collaboration at 7 TeV [48] (blue dots) with our calculations Eqs. (3) and (4) (red triangles). The mg5+cascade3 error bands include statistical and scale uncertainties, see the text below Fig. 3 for more details. The ratio theory/exp is displayed in the bottom panel.

Figure 5: Comparison of the 7 TeV $J/\psi$ differential cross section measured by ALICE [49] (blue dots) and CMS [50] (green dots) with our calculations Eqs. (3) and (4) (red triangles). There is no theoretical prediction in the last bin because of statistics. The mg5+cascade3 error bands include statistical and scale uncertainties, see the text below Fig. 3 for more details.

As discussed in more detail in section 5, we worked with the same value of the non-perturbative parameter $F_{J/\psi}$. Our central values underestimate data at $p_{t}\lesssim 5$ GeV but are still in agreement within uncertainties. The theoretical uncertainties, shown as red and beige color bands, include the variation of factorization and/or renormalization scales, as well as statistical uncertainties. We obtained the red color band by varying the factorization scale between $m_{t}/2$ and $2m_{t}$, with $m_{t}$ the $J/\psi$ transverse mass, while the beige color band is built from the conventional 7-point variation. Both techniques give the same upper limit, but the 7-point variation gives a wider error band.

Let us come back for a moment to our description of low $p_{t}$ data. In Sec. 3, we mentioned that our calculations miss some LPs. It is true for any calculation which does not include the fragmentation contribution $\sigma_{f}\otimes D_{f\to H}$, in particular, for the NRQCD cross section (1). By definition, LP contributions produce the same distribution and are dominant at large $p_{t}$. Then, the missing LPs can be compensated by an appropriate (larger) choice of the parameter $F_{J/\psi}$, which explains our good description of ATLAS data. But then, one could expect an overestimation of low-$p_{t}$ data where the NLPs are dominant. It is clearly not the case. One hypothesis is that the missing LPs do not contribute significantly to the cross section, at least in the transverse momentum range of Fig. 2. Another explanation is that the CEM at NLO misses important contributions at low $p_{t}$, for instance, higher-order contributions or the fragmentation of a gluon emitted during the spacelike cascade into a quarkonium. It would be interesting to see a dedicated study on the role of these different contributions.

Meanwhile, we conclude that the phenomenological CEM, with the evolution of the heavy-quark pair included, can describe data on the whole $p_{t}$ range shown in this study. Again, our main goal is not to defend the CEM, but to show the effect of scale evolution in a simple case. It could be interesting test the CEM further, with less inclusive observables such as the $z$-distribution of a $J/\psi$ within a jet. Fragmentation contributions matched with NRQCD have been compared to this observable, see, for instance, Refs. [51, 52].

5 Universality of $F_{J/\psi}$

In Figs. 6 and 7, we compare our calculations to LHCb data at 13 TeV [53, 54] and PHENIX data at 200 GeV [55]. We used the same value of the non-perturbative parameter, i.e., $F_{J/\psi}=0.014$, and obtained results similar to those of Sec. 4.

Figure 6: Comparison of the $J/\psi$ differential cross section measured by the LHCb collaboration at 13 TeV [53, 54] (blue dots) with our calculations Eqs. (3) and (4) (red triangles). The mg5+cascade3 error bands include statistical and scale uncertainties, see the text below Fig. 3 for more details. The ratio theory/exp is displayed in the bottom panel.

Figure 7: Comparison of the $J/\psi$ differential cross section measured by the PHENIX collaboration at 200 GeV [55] (blue dots) with our calculations Eqs. (3) and (4) (red triangles). The mg5+cascade3 error band includes statistical and scale uncertainties, see the text below Fig. 3 for more details. The ratio theory/exp is displayed in the bottom panel.

We conclude on the universality of $F_{J/\psi}$, at least on this range of energies. On the opposite, using the $k_{t}$-factorization with LO off-shell matrix elements, Ref. [33] reaches a different conclusion, with $F_{J/\psi}(200\text{ GeV})$ about 10 times larger than at 7 TeV (except for the ALICE experiment, with a factor of 2.7). The formalisms used in [33] being quite different, it is hard to pinpoint the precise reason explaining this difference. However, concentrating on small and intermediate $p_{t}$, we can make the following comments. This kinematical region, being dominated by the LO contribution, it seems unlikely that the difference between Ref. [33] and the present work can be explained by the fact that we included the NLO contribution. Another difference is the use of an on-shell cross section instead of an off-shell one. However, these two quantities are supposed to be close at small $k_{t}$.

Another possible explanation is the $k_{t}$ dependence of the KMRW UPDFs [56, 57], already discussed in [58, 59, 60] in the context of D-meson production, see also [61, 62]. These UPDFs, used in [33], are constrained for $k_{t}\in[0,\mu]$ by their relation to collinear PDFs

$$f_{k/h}(x;\mu)=\int_{0}^{\mu^{2}}F_{k/h}(x,k_{t};\mu)dk_{t}^{2}.$$ (8)

Here $f_{k/h}$ gives the collinear distribution of a parton of flavor $k$ in the hadron $h$. On the opposite, the cross section, Eq. (3), is integrated up to $k_{t}^{2}\sim s$. At fixed $\mu$, or equivalently fixed $p_{t}$, and increasing $s$, the unconstrained part of the UPDFs, $k_{t}\in[\mu,\sqrt{s}]$, contributing to the cross section increases. To compensate that, we expect a decrease of $F_{J/\psi}$, which acts as a normalization factor. This decrease of $F_{J/\psi}$ with energy is indeed observed in Ref. [33]. On the opposite, the PB UPDFs used in our study are constrained on the full phase space

$$f_{k/h}(x;\mu)=\int_{0}^{\infty}\frac{F_{k/h}(x,{\bf k}_{t};\mu)}{\pi}d^{2}{\bf k}_{t},$$ (9)

and the numerical value of these functions in the region $[\mu,\sqrt{s}]$ is negligible. Then, the $k_{t}$ integral of Eq. (3) is effectively cut off at $k_{t}\sim\mu$, resulting in a energy independent $F_{J/\psi}$.

Whatever the correct explanation is, it shows the potential of quarkonia data to constrain the (gluon) UPDFs.

6 Comparison with EPOS4

6.1 Short presentation of the formalism

The EPOS4 project is an attempt to construct a realistic model for describing relativistic collisions of different systems, from proton-proton (pp) to nucleus-nucleus (AA), at energies from several TeV per nucleon down to several GeV, see Ref. [63] and references therein. In this section, we discuss the relevant features of this event generator for $J/\psi$ production.

EPOS4 is based on parallel scatterings, with the initial nucleons (and their partonic constituents) being involved, happening instantaneously at very high energies. The theoretical tool is S-matrix theory, using a particular form of the proton-proton scattering S-matrix (Gribov-Regge (GR) approach [64, 65, 15, 66]), which can be straightforwardly generalized to be used for nucleus-nucleus collisions. This S-matrix approach has a modular structure based on so-called “cut Pomerons”, representing inelastic parton-parton scatterings. Each cut pomeron has the ladder structure shown in Fig. 8.

      (a)                                   (b)                                         
      
                       (c)                                                                  
      

      (d)                                   (e)                                         
      

Figure 8: Heavy flavor production in parton ladders (cut pomeron). SLC refers to spacelike cascade, and TLC to timelike cascade. Even if we show mainly gluonic lines, we use a GM-VFN scheme and similar diagrams with quark lines can also be drawn.

Two initial spacelike partons evolve according to the DGLAP equation [15, 16, 17], with the center of the ladder corresponding to a leading-order $2\to 2$ cross section. The timelike partons emitted during the spacelike evolution and those produced in the $2\to 2$ scattering will also evolve according to the DGLAP equation: this is the timelike cascade discussed earlier. In Fig. 8, the heavy quark and antiquark are represented by red lines. We should stop here to compare the different formalisms. Diagrams (a), (b), and (d) are not included by the MG5+cascade3 study. On the other hand, the $gg\to gQ\bar{Q}$ diagram in the center of cut pomeron (e) is taken into account by the NLO cross section provided by MG5. However, EPOS4 works with LO cross section, and the splitting of the gluon into the heavy-quark pair happens in the timelike cascade. It is worth noting that the $Q\bar{Q}$ pair can be produced at any step during the timelike cascade, and that all kind of timelike partons can be produced in the $2\to 2$ scattering. Consequently, and contrary to our MG5+cascade3 calculations, EPOS4 includes contributions similar to the LP contributions $\sigma_{f}\otimes D_{f\to H}$ of the factorization formalism. Finally, the factorization formalism and our MG5+cascade3 calculations do not include, for instance, diagram (a), responsible for the production of low-$p_{t}$ $J/\psi$.

Once the timelike cascade ends, we compute the $J/\psi$ yield according to the CEM: we sum over all $c-\bar{c}$ pairs, compute their masses $m=\sqrt{p(c)\cdot p(\bar{c})}$, and count them with probability $F_{J/\psi}=0.028$ , in case of $m_{J/\psi}<m<2m_{D}$. There is no particular reason for using $m_{J/\psi}$ rather than $2m_{c}$, and it affects only the value of $F_{J/\psi}$. In the next version of EPOS4, we will use $2m_{C}$.

When developing matrix elements in terms of multiple scattering diagrams, the large majority of the diagrams cancel when it comes to inclusive cross sections. This is the Abramovsky-Gribov-Kancheli cancellations [66]. The challenge for EPOS4 is to use the full parallel multiple scattering scenario, but in such a way that for inclusive cross section the cancellations actually work. This is the new part in EPOS4, strongly based on an interplay between parallel scatterings and saturation. It is discussed in a separate publication [67]. But EPOS4 keeps all contributions, and goes beyond inclusive cross sections. This is important, for instance, when studying high multiplicity events, or for $J/\psi$ pair production. We plan to study the latter in an separate work.

6.2 Results

The two main similarities between EPOS4 and our mg5+cascade3 calculations are the use of a timelike cascade, and the hadronization of the produced $Q\bar{Q}$ pairs based on the CEM. Fig. 9 shows that EPOS4 results for ATLAS data are qualitatively the same as mg5+cascade3 calculations.

Figure 9: Comparison of the prompt-$J/\psi$ differential cross section measured by ATLAS [5] as a function of transverse momentum (black dots) with EPOS4 simulations in red. The rapidity range is indicated in the top of each plot.

There are, however, several relevant differences between the two formalisms:

• EPOS4 uses LO cross sections, and the madgraph NLO contribution is taken into account by the timelike cascade.

• It includes all LP contributions. Indeed, EPOS4 produces any kind of partons which evolve with the timelike cascade, and can split into a $Q\bar{Q}$ pair at any time.

• EPOS4 works with a GM-variable-flavor scheme, with $m_{c}=1.27$ GeV. It implies that contributions similar to graph (b) of Fig. 8 are included, while absent for mg5+cascade3.

• The details of the implementation of the timelike cascade may be different.

Despite these differences, both formalisms show similar results, in relatively good agreement with data, and far from the result obtained with NLO fixed-order calculations, Fig. 0(b). The key ingredient is the $Q\bar{Q}$ evolution.

Still, one may be surprised by the fact that both calculations give similar results for ATLAS data, where the LP contributions dominate. While EPOS4 allows a produced gluon to emit, in principle, any number $n$ of partons before splitting into the heavy-quark pair, our calculations with MG5 include only the case $n=0$. The explanation is probably that the produced $Q\bar{Q}$ pair keeps emitting partons, mainly gluons. Then, the global picture is that a color object propagates, emitting partons, and we find a heavy-quark pair at the end of the process. Qualitatively, it does not really matter if the partons are emitted by a gluon or the heavy-quark pair. The final result is energy loss making the spectrum softer. However, as already discussed in Sec. 4, the missing LP contributions in our mg5+cascade3 calculations could affect the value of $F_{J/\psi}$.

Finally, we remark that the probability for the propagating gluon to emit a large number $n$ of gluons before producing the $Q\bar{Q}$ pair is suppressed because the virtuality $Q^{2}$ of timelike partons decreases after each emission. Creating a $Q\bar{Q}$ pair requires at least a virtuality of $Q^{2}_{\text{min}}=4m_{Q}^{2}$, limiting the number of emissions before its production.

7 Complementary figures: Comparison with ATLAS data

In this section, all figures are the same as Fig. 2, but for different rapidity ranges.

8 Conclusion

We have shown that the evolution of the $Q\bar{Q}$ pair predicted by pQCD, and in particular real emissions, solves the issue of NLO calculations with the CEM. However, one should keep in mind that our calculations used several approximations discussed in Sec. 4. Still, there is a little doubt on the significant effect of the evolution on the quarkonium transverse-momentum distribution for factorization scales roughly larger than 20 GeV. Of course, this is also true for the usual fragmentation functions obeying to the DGLAP equation. We expect that the evolution could also have a significant impact on the determination of NRQCD LDMEs, for instance, at large $p_{t}$ or large $Q^{2}=-q^{2}$ in semi-inclusive DIS, where $q$ is the virtual-photon momentum. The latter case is relevant since the EIC will reach $Q$ lager than 20 GeV [68].

In Sec. 5, we showed that with the value $F_{J/\psi}=0.014$, we could describe data from $\sqrt{s}=13$ TeV to $\sqrt{s}=200$ GeV with satisfying accuracy. In other words, $F_{J/\psi}$ seems to be universal, as expected. Other calculations based on $k_{t}$-factorization did not reach this conclusion, and we speculated on the possible explanations. One of these explanations involves the definition of UPDFs, and we believe that quarkonium data have the potential to constrain these functions.