Vector Quarkonia at the LHC with Jethad: A High-Energy Viewpoint

Abstract

In this review, we discuss and extend the study of the inclusive production of vector quarkonia, $J/\psi$ and $\mathrm{{\rm Y}}$, emitted with large transverse momenta and rapidities at the LHC. We adopt the novel ZCW19${}^{+}$ determination of fragmentation functions to depict the quarkonium production mechanism at the next-to-leading level of perturbative QCD. This approach is based on the nonrelativistic QCD formalism well adapted to describe the formation of a quarkonium state from the collinear fragmentation of a gluon or a constituent heavy quark at the lowest energy scale. We rely upon the $\mathrm{NLL}/{\mathrm{NLO}}^{+}$ hybrid high-energy and collinear factorization for differential cross-sections, where the collinear formalism is enhanced by the BFKL resummation of next-to-leading energy logarithms arising in the t-channel. We employ the method to analyze the behavior of the rapidity distributions for double-inclusive vector quarkonium and inclusive vector quarkonium plus jet emissions. We discover that the natural stability of the high-energy series, previously seen in observables sensitive to the emission of hadrons with heavy flavor detected in the rapidity acceptance of LHC barrel calorimeters, becomes even more manifest when these particles are tagged in forward regions covered by endcaps. Our findings present the important message that vector quarkonia at the LHC via hybrid factorization offer a unique chance to perform precision studies of high-energy QCD, as well as an intriguing opportunity to shed light on the quarkonium production puzzle.

Contents
1	Introduction........................................................................................................................	2
2	Vector Quarkonia at the LHC with Jethad....................................................................	6
	2.1      A Glance at High-Energy QCD Phenomenology................................................	6
	2.2      NLL/NLO+ Hybrid Factorization..........................................................................	8
3	NRQCD Collinear Fragmentation for Vector Quarkonia..........................................	11
	3.1      DGLAP-Evolved NRQCD FFs..............................................................................	12
	3.2      Natural Stability.....................................................................................................	14
4	Phenomenology...............................................................................................................	17
	4.1       Jethad v0.5.0: A Basic Overview.........................................................................	17
	4.2      Uncertainty Estimation.........................................................................................	18
	4.3      Final-State Ranges.................................................................................................	19
	4.4      ΔY-Distribution.....................................................................................................	19
5	Conclusions and Outlook..............................................................................................	23
A	Appendix A.......................................................................................................................	25
B	Appendix B.......................................................................................................................	27
C	Reference...........................................................................................................................	28

1. Introduction

It is widely recognized that the emissions of heavy-quark flavored objects in high-energy hadron collisions represent gold-standard channels to access the core dynamics of fundamental interactions. Quarks with heavy masses are sentinels of long-awaited imprints of New Physics, since they might couple with Beyond the Standard Model (BSM) particles. At the same time, since the values of their masses lie in a region where Quantum Chromo Dynamics (QCD) perturbative calculations are possible, they give a unique opportunity to perform precision studies of strong interactions.

QCD is the well-established quantum field theory that describes the strong interaction sector. It relies upon the non-Abelian $SU\left(N_c\right)$ gauge group, with $N_c\equiv 3$ being the number of colors [1,2,3,4]. Quarks are the spin-1/2 fermions that make up the matter part of hadrons. There exist six species of quarks and of their anti-matter counterparts, grouped into three families. In the QCD Lagrangian, they are described by fermionic fields belonging to the fundamental triplet representation of $SU\left(3\right)$. Conversely, gluons are the spin-1 massless bosons that mediate the strong force. They bind quarks together to form hadrons. In the QCD Lagrangian, they are described by bosonic fields belonging to the adjoint octet representation of $SU\left(3\right)$.

While QCD represents one of the founding pillars of the Standard Model (SM) of fundamental interactions, it can also serve as a “support basis” in relevant searches for New Physics (see Ref. [5] for a recent and comprehensive resource letter). An incomplete list of possible QCD portals beyond SM is as follows: axions, namely hypothetical elementary particles, originally proposed by R. Peccei and H. Quinn to resolve the strong Charge-Parity (CP) problem [6,7,8,9]; Non-Abelian dark gauge forces [10,11]; quarkyonic matter [12,13,14]; and higher-dimension QCD operators, effectively included in the Lagrangian (see, e.g., Refs. [15,16,17,18] and references therein).

Unraveling the production mechanisms of heavy-quark flavored hadrons plays a key role in understanding the true dynamics of strong interactions. Mesons whose lowest Fock state contains a heavy quark (Q) and its antiquark ($\overline{Q}$) are known as (heavy) quarkonia. The origins of quarkonium studies date back to the so-called November Revolution. On 11 November 1974, a new vector meson, whose mass was around 3.1 GeV, carrying photon quantum numbers, was discovered. It was named $J/\psi$ after the two groups that contemporaneously observed it: the collaboration headed by B. Richter [19] at the Stanford Linear Accelerator Center (SLAC) and one led by S. Ting [20] at the Brookhaven National Laboratory (BNL). A few days after the announcements by SLAC and BNL, the discovery of the $J/\psi$ was confirmed by the Frascati ADONE experiment, directed by G. Bellettini [21]. The fact that the $J/\psi$ was a hadron became clear from the analysis of the ratios of the cross-sections between hadrons and ${\mu }^{+}{\mu }^{-}$ in $e^{+}{e}^{-}$ inclusive annihilations, which eventually were much larger at the resonance, thus pointing out that $J/\psi$-like objects directly decay into lighter hadrons. At the same time, evidence was found that the charm quark exists and that quarks are real particles of which hadrons are made, and not simply mathematical artifacts. The concept of charm flavor was proposed for the first time by Bjorken and Glashow in 1964 [22], but its experimental evidence was scarce. Then, in 1970, the Glashow–Iliopoulos–Maiani (GIM) mechanism [23] proposed the existence of a new quark species, the charm flavor, which became necessary to explain the suppression of Flavor-Changing Neutral Currents (FCNCs) in Feynman diagrams embodying a loop.

Studies on quarkonium properties and formation mechanisms have been relevant to validate our knowledge of the fundamental properties of QCD, such as asymptotic freedom. In particular, the discovery of the first excited radial state of the $J/\psi$, called $\psi \left(2S\right)$, brought clear evidence that the strong force lowers at short distances. The $J/\psi$ represents the first discovered charmonium. This name reflects an affinity with the $e^{+}{e}^{-}$ bound pair, called the positronium, and it stands for an entire species of hadrons comprehending all the $c\overline{c}$ meson states. A few years after the first $J/\psi$ observation, other charmonia were discovered (such as the P-wave ${\chi }_c$ and the pseudoscalar ${\eta }_c$), as well as hadrons with open charm, such as D mesons [24]. The first bottomonium particle was observed in 1977 and was named $\mathrm{{\rm Y}}$. It represents a vector meson equivalent to the $J/\psi$, but with bottom quarks instead of charm quarks [25]. Excited $b\overline{b}$ states, such as the $\mathrm{{\rm Y}}\left(2S\right)$, and open-bottom B mesons were subsequently observed [26].

Toponium systems are hypothetical $t\overline{t}$ mesons whose lowest Fock state consists of a top and an anti-top quark. As for charmonia and bottomonia, toponia should also appear as vectors, called $\theta$ mesons, or scalars, called ${\eta }_t$ mesons (for a short overview, see Ref. [27] and references therein). At variance with lighter, observed quarkonia, $t\overline{t}$ mesons are thought to decay instantly since top quarks are short-lived [28,29,30,31,32]. Hence, they possess a very large decay width [33,34]. Estimates for toponium emission rates at the LHC were presented in Refs. [35,36,37,38]. Possible mechanisms for spin-1 and spin-0 toponium production at future lepton colliding machines were proposed in Refs. [39,40,41] and Ref. [42], respectively. Ref. [43] reports deviations between predictions and observations related to the production of top pairs at the LHC, followed by double-lepton decays. It also proposes a novel method to unveil toponium formation by reconstructing both top and anti-top constituents in dilepton hadroproduction. Ref. [44] focuses on probing CP violation via toponium leptoproduction.

A clear advantage in studying quarkonia comes from the fact that they can be easily detected, particularly vector mesons, which can be easily identified in lepton-initiated processes. At the same time, however, describing their production rates is a challenge. Difficulties encountered in attempting to describe quarkonium formation in distinct kinematic ranges give rise to the so-called quarkonium production puzzle. Several models attempting to capture the core features of the quarkonium formation mechanism have been proposed so far (see Refs. [45,46] for a review). The Color Evaporation Model (CEM) [47,48] was perhaps the first mechanism proposed. It assumes that the original $\left(Q\overline{Q}\right)$ pair, from which the final quarkonium stems, goes through a very large number of soft interactions, which eventually permits it to be a color singlet, after hadronization. These soft subprocesses completely decorrelate the initial $\left(Q\overline{Q}\right)$ color state with the final meson one. The predictive power of CEM is limited by the inability to provide us with a distinct description of the production rates of different quarkonia, such as $J/\psi$ and ${\chi }_c$ states [49,50].

The second model proposed was the Color Singlet Mechanism (CSM) [51,52,53,54]. In contrast with CEM, it postulates that the $\left(Q\overline{Q}\right)$ pair does not undergo any soft evolution. This happens because gluon emissions are suppressed with powers of ${\alpha }_s\left(m_Q\right)$, with ${\alpha }_s\left(\mu \right)$ the QCD running coupling and $m_Q$ the mass of the constituent heavy quark. Both the quark polarization and its color state do not change during the hadronization phase. Thus, the originating pair must be generated in a color-singlet configuration. Moreover, since the measured quarkonium mass, $m_{\mathcal{Q}}$, is slightly larger than $2m_Q$, the two constituent heavy quarks are almost at rest in the hadron frame, namely with almost zero relative velocity, $v_{\mathcal{Q}}$. The only nonperturbative contribution to the formation mechanism is a nonrelativistic Schrödinger wave function at the origin, labeled as ${\mathcal{R}}_{\mathcal{Q}}\left(0\right)$. While the CSM genuinely possesses high predictive power (its only free parameter is ${\mathcal{R}}_{\mathcal{Q}}\left(0\right)$, which can be extracted from experimental data on leptonic decay widths), it presents infrared divergences emerging at higher perturbative order in P-wave decay channels  [55,56].

To overcome these divergence issues associated with the CSM, it was argued that higher Fock states, such as $|Q\overline{Q}g\rangle$ and so on, could play a role in the production mechanism. Contextually, they could solve Next-to-Leading Order (NLO) divergences. The $\left(Q\overline{Q}\right)$ subsystem contained in a $|Q\overline{Q}g\rangle$ state is clearly in a Color Octet (CO) state. More generally, the physical quarkonium is built up as a linear combination of all the possible Fock states. All these terms are ordered within a double expansion in powers of both ${\alpha }_s$ and $v_{\mathcal{Q}}$. In this picture, the CSM represents the leading contribution in $v_{\mathcal{Q}}$ for S-wave channels (for higher waves, this term might vanish). The $({\alpha }_s\otimes v_{\mathcal{Q}})$ expansion constitutes the building block of NonRelativistic QCD (NRQCD) [57,58,59,60,61,62,63]. According to this effective field theory, high-energy distributions sensitive to quarkonium emissions are cast as sums of partonic subprocess hard factors producing a given Fock state in the final state, each of them being multiplied by a Long-Distance Matrix Element (LDME) embodying the nonperturbative hadronization mechanism. Despite giving an elegant solution to the quarkonium production puzzle, the NRQCD introduces many new nonperturbative parameters, which can make the comparison with data cumbersome.

What makes the overall picture even more intricate is the possibility of having a quarkonium produced through the decay or de-excitation of heavier particles. Thus, $J/\psi$ production can follow the decay of a B meson [64]. This channel is labeled as nonprompt, since it features a visible time delay. In some circumstances, the nonprompt contribution can be experimentally isolated by measuring the distance between the creation vertex of the B and its subsequent decay or de-excitation. In this way, one obtains the prompt component. Within the prompt channel, a charmonium, e.g., the $J/\psi$, can still be emitted as the decay product of a ${\chi }_c$, or via the de-excitation of a $\psi \left(2S\right)$. This represents the indirect production component. The total number of quarkonia produced in hard subprocesses is called the direct fraction and it is generally not measurable alone but can be accessed by subtraction methods.

All the models introduced above are based on the assumption that the dominant formation mechanism is the production of a $\left(Q\overline{Q}\right)$ pair in the hard scattering, namely at short distances, followed by its hadronization into the physical quarkonium state. The two constituent heavy quarks are generated with a relative transverse-momentum separation of order $1/|{\overrightarrow{p}}_T|$. Therefore, short-distance mechanisms are expected to be suppressed at large $|{\overrightarrow{p}}_T|$. Furthermore, the two quarks have a separation of order $1/\mu$, with $\mu$ any process-typical energy scale [65]. At large transverse momenta, one has $\mu \sim |{\overrightarrow{p}}_T|\gg m_Q$. In this way, since the exchanges of particles with virtualities of the order of $|{\overrightarrow{p}}_T|$ are related to the scatterings of particles in a small volume, $1/|{\overrightarrow{p}}_T{|}^3$ [66,67], they act as suppression factors for the probability amplitude.

The authors of Ref. [68] pointed out that, in moderate or large transverse-momentum regimes, an additional mechanism is at work. It relies upon the fragmentation of a single parton, which then inclusively hadronizes into the observed quarkonium. According to collinear factorization, which is the most appropriate formalism to describe these regimes, the full hadronization must be encoded in collinear Fragmentation Functions (FFs), whose time-like evolution [69,70] is naturally governed by the Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) system of equations [71,72,73,74,75]. Here, NRQCD turns out to be a powerful tool to perturbatively calculate nonevolved, initial scale inputs for these FFs. Now, the heavy-quark pair has a separation of order $1/m_Q$. On the one hand, fragmentation production is often of a higher perturbative order than the short-distance mechanism. On the other hand, fragmentation is magnified by powers of ${(\mu /m_Q)}^2$. Thus, it becomes dominant at large energy scales, $\mu \gg m_Q$.

NRQCD calculations of the gluon FF to an S-wave particle in the CSM were presented in Ref. [68] for $J/\psi$ and ${\eta }_c$ at Leading Order (LO), and in Ref. [76] for ${\eta }_c$ with higher accuracy (NLO). The corresponding calculation for the $(c\to J/\psi +c+\mathcal{X})$ FF was given in Ref. [77] at LO and in Ref. [78] at NLO. Refs. [79,80,81,82,83] contain phenomenological analyses aimed at shedding light, at the LO, on the intersection range between the short-distance and the fragmentation approximations. Functions portraying the fragmentation of constituent heavy quarks into S- and P-wave quarkonia were studied in Refs. [84,85]. We refer to Refs. [86,87] for some FF advancements at next-to-NLO. FFs for polarized quarkonia were investigated in Refs. [88,89,90,91,92,93] (see Ref. [94] for a digression). The authors of Ref. [95] analyze quarkonium factorization and evolution at large energies.

All the issues found on the path towards capturing the correct production mechanism(s) do not prevent quarkonium studies from providing gold-standard tools to access the inner dynamics of QCD. Pioneering analyses of quarkonium emissions within collinear factorization at NLO were performed in Refs. [32,96,97,98,99,100,101,102,103,104]. Refs. [105,106] study the NLO transverse-momentum spectrum in the CEM approximation. NRQCD advancements in quarkonium polarization can be found in Refs. [107,108,109]. The possible emergence of Multi-Parton Interaction (MPI) dynamics in quarkonium observables was quantified in Refs. [110,111]. The weight of nuclear modifications of the gluon density on quarkonium was gauged for the first time in Refs. [112,113,114]. Automated NRQCD calculations were recently implemented [115,116,117,118]. Problems in the NLO description of $J/\psi$ distributions at large transverse momenta can spread from an over-subtraction in the factorization of collinear singularities inside collinear Parton Distribution Functions (PDFs) [119,120,121]. Matching high-energy (low-x) resummation effects on top of NLO calculations recently emerged as a promising means to solve these issues [122,123]. Color reconnection effects in $J/\psi$ hadroproduction have been recently investigated [124]. Some recent results on forward and central quarkonium production in (ultra-peripheral) hadron and lepton–hadron scatterings can be found in Refs. [125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142]. The formation mechanism of hidden heavy-flavored tetra- and penta-quarks was studied [143,144] by means of the hadro-quarkonium approach [145,146]. NRQCD-based studies have been used to model tetra-quark fragmentation [147,148,149,150].

Observables sensitive to the single forward production of quarkonium states, as well as those for forward light vector mesons [151,152,153,154,155,156] and forward Drell–Yan pairs [157,158,159], are useful tools to extensively scan the intersection kinematic corners between the Transverse-Momentum-Depende (TMD) and the high-energy dynamics. Studies given in Refs. [154,156] unambiguously pointed out that spin-dependent distributions for the single forward electroproduction of a $\rho$ meson at HERA and the Electron-Ion Collider (EIC) give us a unique opportunity to access the content and properties of low-x gluons in the proton. There, transverse-momentum factorized cross-sections are sensitive to the scattering of color dipoles with a small transverse distance. Contrariwise, forward detections of excited states, such as $\psi \left(2S\right)$ and $\mathrm{{\rm Y}}\left(2S\right)$ quarkonia, might probe dipoles with larger transverse sizes [160,161,162], thus calling for a low-$|{\overrightarrow{p}}_T|$ improvement in the pure small-x picture.

Quarkonium emissions at low transverse momenta in lepton–hadron and hadron–hadron collisions offer a unique opportunity to perform 3D reconstructions of the nucleon content by means of transverse-momentum-dependent gluon densities [163,164,165,166]. We refer to Refs. [167,168,169,170,171,172,173,174,175,176,177,178] for recent phenomenological advancements. The validity of TMD factorization in the presence of quarkonium detections was recently questioned via a Soft and Collinear Effective Theory (SCET) [179,180] and then evaluated for di-jet and heavy-meson pair tags in lepton–proton reactions [181]. An SCET approach turns out to be useful to access jet fragmenting functions, a basic component to study quarkonia inside jets [182,183,184,185,186]. Quarkonium TMD distributions were recently investigated in Ref. [187].

In this review, we consider the inclusive hadroproduction, at LHC energies and kinematic configurations, of a forward vector quarkonium ($J/\psi$ or $\mathrm{{\rm Y}}$ in the direct channel) in association with a backward particle, which can be another vector quarkonium or a light-flavored jet. The two final-state objects have transverse momenta and are largely separated in rapidity. On the one hand, moderate values of partons’ longitudinal-momentum fractions make a description in terms of collinear PDFs valid. On the other hand, large rapidity intervals translate into significant exchanges of transverse momenta in the t-channel, which call for a $|{\overrightarrow{p}}_T|$-factorization treatment aimed at resumming large logarithmic contributions in energy, associated with rapidity-ordered gluon emissions. As a last step, we must select the most adequate mechanism for the quarkonium formation in the kinematic sector of our interest. From a collinear factorization perspective, large transverse momenta permit us to adopt a Variable-Flavor Number Scheme (VFNS) approach [188,189]. Here, the hard factor for the production of a light parton is convoluted with an FF that describes the hadronization into the vector quarkonium.

Following the recent scheme proposed in Ref. [190], we first make use of NRQCD to build initial-scale FF inputs, and then we encode the resummation of DGLAP-type logarithms. We take the NRQCD input from a recent NLO computation of the constituent heavy-quark FF channel [78], together with the gluon one [68]. We stress that these NRQCD FFs describe a vector quarkonium in a color-singlet state. In other words, they embody ${}^3{S}_1^{\left(1\right)}$ LDMEs (we refer to Section 3.1 for more details). Putting together all the components, we define a hybrid high-energy and collinear factorization, already set as the reference formalism for the study of semi-inclusive forward–backward hadroproduction [191,192,193,194,195], now enriched with a novel element: the single-parton quarkonium fragmentation [190].

Another formalism, close in spirit to our hybrid factorization and suitable for single forward detections, was proposed in Refs. [196,197,198,199,200]. We also refer to Refs. [201,202] for analyses of small-x resummed inclusive or differential distributions for Higgs and heavy-flavor hadroproduction via the Hell method [203,204,205], which relies upon the Altarelli–Ball–Forte (ABF) approach [206,207,208,209,210,211,212], aimed at combining collinear factorization with small-x resummation, based upon high-energy factorization theorems [213,214,215,216,217,218,219,220]. We believe that a strong formal connection between our hybrid factorization and the ABF formalism exists. Studies along this direction are relevant, but they go beyond the scope of this review. We postpone them for a future work.

Our study is complementary to the one presented in Ref. [221], where the inclusive semi-hard $J/\psi$ plus jet process was investigated within the short-distance approximation and with partial NLO accuracy. In the future, matching between the two approaches will certainly improve the description of quarkonium detections at high energies. For the sake of completeness, we mention the computation of doubly off-shell LO coefficient functions for inclusive quarkonium production in central rapidity regions [222,223].

This review is structured as follows. In Section 2.2, we introduce the $\mathrm{NLL}/{\mathrm{NLO}}^{(+)}$ hybrid high-energy and collinear factorization, while Section 3 contains technical details on the vector quarkonium collinear fragmentation mechanism, based on constituent heavy-quark and gluon initial-scale inputs on top of which standard DGLAP evolution is switched on. In Section 4, we present our phenomenological analysis of rapidity-differential distributions. Then, in Section 5, we draw our conclusions and highlight future perspectives.

2. Vector Quarkonia at the LHC with Jethad

This part of the manuscript contains a brief review of recent key advancements in the phenomenology in the semi-hard sector of QCD (Section 2.1), as well as technical details of the hybrid high-energy and collinear factorization for inclusive vector quarkonium and vector quarkonium + jet hadroproduction at $\mathrm{NLL}/{\mathrm{NLO}}^{+}$ (Section 2.2).

2.1. A Glance at High-Energy QCD Phenomenology

One of the greatest results in high-energy particle physics is the possibility to decouple, in the context of hadron scatterings, the long-distance from the short-distance dynamics. This decoupling allows us to factorize nonperturbartive components from perturbative computations through the well-established collinear picture. At the same time, there exist kinematic regions where large logarithms arise. They can become so large as to compensate for the smallness of the QCD running coupling, thus spoiling the convergence of the perturbative series. Therefore, the standard collinear factorization needs to be enhanced by the inclusion of one or more all-order resummations.

We turn our attention to the so-called semi-hard regime [224] (see Refs. [225,226,227] for relevant applications), where the scale hierarchy $\sqrt{s}\gg \left\{Q\right\}\gg {\Lambda }_{\mathrm{QCD}}$ stringently holds. Here, s is the center-of-mass energy squared, $\left\{Q\right\}$ a set of process-characteristic hard scales, and ${\Lambda }_{\mathrm{QCD}}$ the QCD scale. While the second inequality simply tells us that the use of perturbation techniques is allowed, the first highlights that we have entered the Regge limit of QCD. In this case, large logarithms of the form $lns/Q^2$ enter the perturbative series with a power growing with the order. When ${\alpha }_s\left(Q^2\right)ln(s/Q^2)\sim 1$, pure fixed-order perturbative calculations are no longer reliable, but it must be improved by resumming these large-energy logarithms to all orders. The most powerful mechanism for this high-energy resummation is represented by the Balitsky–Fadin–Kuraev–Lipatov (BFKL) formalism [228,229,230,231], which tells us how to resum all the contributions proportional to ${({\alpha }_{s}lns)}^n$, the so-called Leading Logarithmic (LL) approximation, as well as the ones proportional to ${\alpha }_s{({\alpha }_{s}lns)}^n$, the so-called Next-to-Leading Logarithmic (NLL) approximation. In the BFKL formalism, a given scattering amplitude can be written as a convolution of a process-independent Green’s function with two singly-off-shell coefficient functions depicting the transition from each colliding particle to the respective final-state object. Using the BFKL jargon, these coefficients are known as forward-production impact factors. The BFKL Green’s function is controlled an integral evolution equation, whose kernel is known at the NLO for any fixed, not growing with energy, momentum transfer, t, and for any possible two-gluon color exchange in the t-channel [232,233,234,235,236,237,238]. Recent, partial advancements of the BFKL next-to-NLL level can be found in Refs. [239,240,241,242,243].

Despite the full NLO accuracy of the kernel, the predictive power of the BFKL at NLL is limited by the number of impact factors computed within NLO: (1) colliding parton (quarks and gluons) impact factors [244,245], namely the core component for (2) forward jet [246,247,248,249,250,251] and (3) forward light hadron [252] impact factors, then (4) ${\gamma }^{*}$ to Light Vector Meson (LVM) transition [253] and (5) ${\gamma }^{*}$ to ${\gamma }^{*}$ impact factors [254,255,256,257,258,259], and finally (6) the impact factor depicting the emission of a forward Higgs boson in the infinite top-mass limit [260,261,262] (see also Refs. [263,264]). The forward production impact factors known at LO only are those for (7) Drell–Yan pairs [157,265], (8) heavy-quark pairs hadro- and photo-produced  [266,267,268], and (9) forward $J/\psi$ emitted low transverse momentum in the short-distance approximation [221] (see also Refs. [269,270,271]).

On the one hand, these off-shell coefficient functions have been employed for phenomenological studies of several semi-inclusive reactions featuring forward plus backward two-particle final states. An incomplete list is as follows: the diffractive exclusive electroproduction of two light vector mesons [272,273,274,275,276], the inclusive hadroproduction of two jets with large transverse momenta and well separated in rapidity (Mueller–Navelet channel [277]), for which several analyses have been conducted so far (see, e.g., Refs. [191,249,251,278,279,280,281,282,283,284,285,286,287,288,289,290]), the inclusive emission of two light-charged hadrons [291,292,293,294,295], three- and four-jet hadroproduction [296,297,298,299,300,301,302,303,304,305,306,307,308], $J/\psi$ plus jet [221], hadron plus jet [192,309,310,311,312,313], Higgs plus jet [314,315,316,317,318], heavy-light dijet systems [193,319], forward Drell–Yan pairs with a possible backward-jet detection [320], and heavy-flavored hadrons’ hadroproduction [190,194,221,266,267,268,321,322,323,324,325,326,327].

On the other hand, single forward emissions permit us to unveil the proton content at low x via the BFKL Unintegrated Gluon Distribution (UGD), whose evolution is driven by Green’s function. Gold-standard channels for the testing of the UGD are the exclusive light vector–meson leptoproduction at HERA [151,152,153,154,155,328,329,330,331] and the EIC [156,332,333,334,335], the exclusive quarkonium photoemission [162,336,337], and the inclusive detection of Drell–Yan dilepton systems [157,158,159,338]. Then, the information on the gluon content at small x carried by the BFKL UGD has allowed us to enhance the collinear description via the first determination of small-x resummed PDFs [339,340,341], as well as to model low-x improved gluon TMD densities at leading twist, namely within the lowest level of the operator product expansion of the gluon correlator [342,343,344,345,346,347,348,349,350]. We refer to Refs. [351,352] for recent investigations of the interplay between the BFKL dynamics and the TMD factorization, and to Ref. [353] for a study of the connection between the deep-inelastic-scattering dipole cross-section and the UGD.

A major issue emerging from the study of Mueller–Navelet distributions is the size of the NLL corrections, which are of the same order but generally with an opposite sign with respect to the LL background. This leads to heavy instabilities in the high-energy series arising when renormalization and factorization scale variations are moved from their natural values, suggested by process kinematics. Thus, the differential cross-sections can easily become negative when the rapidity interval between the two jets increases. Furthermore, azimuthal-angle correlations turn out to be unphysical, both in small and large rapidity interval ranges. Several strategies have been proposed to solve these issues. Among them, the Brodsky–Lepage–Mackenzie (BLM) method [354,355,356,357], as specifically designed for semi-hard processes [282], has become quite popular, since it permits us to lightly suppress these instabilities on azimuthal correlations and to moderately improve the agreement with experimental data. However, the use of BLM is almost ineffective on cross-sections for light di-hadron and light hadron jet observables. This happens because the optimal renormalization scales are substantially higher than the natural scales [192,225,309]. As a result, one observes a significant loss of statistics for total cross-sections. Therefore, any effort at reaching the precision level is unsuccessful.

First, unambiguous imprints of fair stability of the high-energy resummation under higher-order corrections and energy-scale variations were observed only recently in Large Hadron Collider (LHC) semi-inclusive final states featuring the tags of objects possessing large transverse masses, such as Higgs bosons [314,358,359,360,361]. Subsequently, a very remarkable result at NLL was obtained by studying semi-hard correlations for a singly heavy-flavored hadron species, namely ${\Lambda }_c$ baryons. A strong stabilization pattern emerged from an analysis of double ${\Lambda }_c$ and ${\Lambda }_c$ plus jet productions at the LHC [194], and also from similar distributions sensitive to single-bottomed hadrons [321]. Indisputable evidence was provided that the characteristic trend of VFNS collinear FFs describing the formation of those heavy hadrons at large transverse momenta translates into the natural stabilization of the high-energy series, with the fair dampening of those instabilities associated with higher-order corrections. Analogous patterns were then observed for vector quarkonium [190] and charmed B-meson [325] semi-hard observables built up by combining the high-energy resummation with collinear PDFs and DGLAP-evolved FFs with an NRQCD initial-scale input. This suggests that the found natural stability is an intrinsic feature globally shared by heavy-flavor observables.

2.2. NLL/NLO${}^{+}$ Hybrid Factorization

We consider the reactions as in Figure 1:

$$\begin{array}{cc}\hfill \mathrm{p}\left(P_a\right)+\mathrm{p}\left(P_b\right)& \to \mathcal{Q}(p_{{\mathcal{Q}}_1},y_{{\mathcal{Q}}_1})+\mathcal{X}+\mathcal{Q}(p_{{\mathcal{Q}}_2},y_{{\mathcal{Q}}_2}),\hfill \\ \hfill \mathrm{p}\left(P_a\right)+\mathrm{p}\left(P_b\right)& \to \mathcal{Q}(p_{\mathcal{Q}},y_{\mathcal{Q}})+\mathcal{X}+\mathrm{jet}(p_J,y_J),\hfill \end{array}$$

(1)

with $\mathrm{p}\left(P_{a,b}\right)$ being an incoming proton with four-momentum $P_{a,b}$, $\mathcal{Q}(p_{{\mathcal{Q}}_{(1,2)}},y_{{\mathcal{Q}}_{(1,2)}})$ denoting a vector quarkonium ($J/\psi$ or $\mathrm{{\rm Y}}$) produced with four-momentum $p_{{\mathcal{Q}}_{(1,2)}}\equiv p_{1,2}$ and rapidity $y_{{\mathcal{Q}}_{(1,2)}}\equiv y_{1,2}$, and the light-favored jet being detected with four-momentum $p_J\equiv p_2$ and rapidity $y_J\equiv y_2$. The inclusive, undetected remnant is $\mathcal{X}$. The large final-state transverse momenta, $|{\overrightarrow{p}}_{1,2}|$, together with the high rapidity interval, $\Delta Y\equiv y_1-y_2$, are necessary conditions to make the final state diffractive and semi-hard. Moreover, the transverse-momentum ranges must to be large enough to preserve the validity of the single-parton collinear fragmentation as a leading mechanism for quarkonium production.

Incoming protons’ four momenta can be taken as Sudakov vectors obeying $P_a^2=P_b^2=0$ and $2(P_a·P_b)=s$, so that the four momenta of the observed object can be decomposed in the following way:

$$p_{1,2}=x_{1,2}{P}_{a,b}+\frac{{\overrightarrow{p}}_{1,2}^2}{x_{1,2}s}{P}_{b,a}+p_{1,2\perp },p_{1,2\perp }^2=-{\overrightarrow{p}}_{1,2}^2.$$

(2)

Here, outgoing particle longitudinal momentum fractions, $x_{1,2}$, are related to the corresponding rapidities through $y_{1,2}=\pm \frac{1}{2}ln\frac{x_{1,2}^{2}s}{{\overrightarrow{p}}_{1,2}^2}$, so that $\mathrm{d}{y}_{1,2}=\pm \frac{\mathrm{d}{x}_{1,2}}{x_{1,2}}$, and $\Delta Y\equiv y_1-y_2=ln\frac{x_1{x}_{2}s}{|{\overrightarrow{p}}_1\left|\right|{\overrightarrow{p}}_2|}$.

A purely collinear treatment in QCD would cause the LO cross-sections of our reactions (Equation (1)) to be written as a convolution among the partonic hard factor, the proton PDFs, and the quarkonium FFs. In the double quarkonium channel (panel a) of Figure 1), one has

$$\begin{array}{cc}\hfill & \frac{\mathrm{d}{\sigma }_{[p+p\to \mathcal{Q}+\mathcal{Q}]}^{\mathrm{LO}}}{\mathrm{d}{x}_1\mathrm{d}{x}_2{\mathrm{d}}^2{\overrightarrow{p}}_1{\mathrm{d}}^2{\overrightarrow{p}}_2}=\sum _{a,b}{\int }_0^1\mathrm{d}{x}_a{\int }_0^1\mathrm{d}{x}_b{f}_a(x_a,{\mu }_F)f_b(x_b,{\mu }_F)\hfill \\ \hfill & \times {\int }_{x_1}^1\frac{\mathrm{d}{z}_1}{z_1}{\int }_{x_2}^1\frac{\mathrm{d}{z}_2}{z_2}{D}_a^{\mathcal{Q}}\left(\frac{x_1}{z_1},{\mu }_F\right)D_b^{\mathcal{Q}}\left(\frac{x_2}{z_2},{\mu }_F\right)\frac{\mathrm{d}{\widehat{\sigma }}_{a,b}}{\mathrm{d}{x}_a\mathrm{d}{x}_b\mathrm{d}{z}_1\mathrm{d}{z}_2{\mathrm{d}}^2{\overrightarrow{p}}_1{\mathrm{d}}^2{\overrightarrow{p}}_2},\hfill \end{array}$$

(3)

where the ($a,b$) indices run over quarks, antiquarks, and the gluon; $f_{a,b}\left(x,{\mu }_F\right)$ denote parent protons’ PDFs; $D_{a,b}^{\mathcal{Q}}\left(\frac{x}{z},{\mu }_F\right)$ denote quarkonium FFs; $x_{a,b}$ are the longitudinal-momentum fractions of the process-initiating partons; $z_{1,2}$ are the longitudinal-momentum fractions of fragmenting partons; and $\mathrm{d}{\widehat{\sigma }}_{a,b}(x_a,x_b)$ is the partonic hard factor. In the same way, in the quarkonium plus jet channel (panel b) of Figure 1), one writes

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\sigma }_{[p+p\to \mathcal{Q}+\mathrm{jet}]}^{\mathrm{LO}}}{\mathrm{d}{x}_1\mathrm{d}{x}_2{\mathrm{d}}^2{\overrightarrow{p}}_1{\mathrm{d}}^2{\overrightarrow{p}}_2}& =\sum _{a,b}{\int }_0^1\mathrm{d}{x}_a{\int }_0^1\mathrm{d}{x}_b{f}_a(x_a,{\mu }_F)f_b(x_b,{\mu }_F)\hfill \\ \hfill & \times {\int }_{x_1}^1\frac{\mathrm{d}z}{z}{D}_a^{\mathcal{Q}}\left(\frac{x_1}{z}\right)\frac{\mathrm{d}{\widehat{\sigma }}_{\alpha ,\beta }}{\mathrm{d}{x}_1\mathrm{d}{x}_2\mathrm{d}z{\mathrm{d}}^2{\overrightarrow{p}}_1{\mathrm{d}}^2{\overrightarrow{p}}_2}.\hfill \end{array}$$

(4)

Conversely, to obtain the formula for the high-energy resummed cross-section in our hybrid factorization, we first apply the high-energy factorization that is genuinely embodied in the BFKL framework, and then we complement the description by plugging collinear inputs, PDFs, and FFs. It is convenient to rewrite the differential cross-section as a Fourier sum of azimuthal-angle coefficients

$$\frac{\mathrm{d}{\sigma }^{\mathrm{NLL}/{\mathrm{NLO}}^{+}}}{\mathrm{d}{y}_1\mathrm{d}{y}_2\mathrm{d}{\overrightarrow{p}}_1\mathrm{d}{\overrightarrow{p}}_2\mathrm{d}{\phi }_1\mathrm{d}{\phi }_2}=\frac{1}{{\left(2\pi \right)}^2}\left[{\mathcal{C}}_0^{\mathrm{NLL}/{\mathrm{NLO}}^{+}}+2\sum _{n=1}^{\infty }cos\left(n(\mathrm{\Phi }-\pi )\right){\mathcal{C}}_n^{\mathrm{NLL}/{\mathrm{NLO}}^{+}}\right],$$

(5)

where ${\phi }_{1,2}$ are final-state particle azimuthal angles and $\mathrm{\Phi }\equiv {\phi }_1-{\phi }_2$. The azimuthal coefficients are calculated in the BFKL approach and they embody the resummation of energy logarithms up to the NLL level. Working in the $\overline{\mathrm{MS}}$ renormalization scheme [363], one obtains (for the technical derivation see, e.g., Ref. [249])

$$\begin{array}{cc}\hfill {\mathcal{C}}_n^{\mathrm{NLL}/{\mathrm{NLO}}^{+}}& ={\int }_0^{2\pi }\mathrm{d}{\phi }_1{\int }_0^{2\pi }\mathrm{d}{\phi }_{2}cos\left(n(\mathrm{\Phi }-\pi )\right)\frac{\mathrm{d}{\sigma }^{\mathrm{NLL}/{\mathrm{NLO}}^{+}}}{\mathrm{d}{y}_1\mathrm{d}{y}_2\mathrm{d}|{\overrightarrow{p}}_1\left|\mathrm{d}\right|{\overrightarrow{p}}_2|\mathrm{d}{\phi }_1\mathrm{d}{\phi }_2}\hfill \\ \hfill & =\frac{e^{\Delta Y}}{s}{\int }_{-\infty }^{+\infty }\mathrm{d}\nu e^{\Delta Y{\overline{\alpha }}_s\left({\mu }_R\right){\chi }^{\mathrm{NLO}}(n,\nu )}\hfill \\ \hfill & \times {\alpha }_s^2\left({\mu }_R\right)\{c_1^{\mathrm{NLO}}(n,\nu ,|{\overrightarrow{p}}_1|,x_1\left)\right[c_2^{\mathrm{NLO}}(n,\nu ,|{\overrightarrow{p}}_2|,x_2{\left)\right]}^{*}\hfill \\ \hfill & +\left.{\overline{\alpha }}_s^2\left({\mu }_R\right)\Delta Y\frac{{\beta }_0}{4N_c}\chi (n,\nu )f\left(\nu \right)\right\},\hfill \end{array}$$

(6)

where ${\overline{\alpha }}_s\left({\mu }_R\right)\equiv {\alpha }_s\left({\mu }_R\right)N_c/\pi$ with $N_c$ the number of colors; then, ${\beta }_0=11N_c/3-2n_f/3$ as the first coefficient of the QCD $\beta$-function with $n_f$ the number of flavors. A two-loop running-coupling setup with ${\alpha }_s\left(M_Z\right)=0.11707$ and a dynamic $n_f$ is adopted. The BFKL kernel at the exponent of Equation (6) encodes the resummation of leading and next-to-leading energy logarithms. It reads

$$\begin{array}{c}\hfill {\chi }^{\mathrm{NLO}}(n,\nu )=\chi (n,\nu )+{\overline{\alpha }}_s\widehat{\chi }(n,\nu ),\end{array}$$

(7)

with

$$\begin{array}{c}\hfill \chi \left(n,\nu \right)=-2{\gamma }_{\mathrm{E}}-2\mathrm{Re}\left\{\psi \left(\frac{1+n}{2}+i\nu \right)\right\}\end{array}$$

(8)

being the LO BFKL eigenvalues, ${\gamma }_{\mathrm{E}}$ the Euler–Mascheroni constant, and $\psi \left(z\right)\equiv {\mathrm{\Gamma }}^{\prime }\left(z\right)/\mathrm{\Gamma }\left(z\right)$ the logarithmic derivative of the Gamma function. Then, the $\widehat{\chi }(m,\nu )$ function in Equation (7) denotes the NLO correction to the BFKL kernel

$$\widehat{\chi }\left(n,\nu \right)=\overline{\chi }(n,\nu )+\frac{{\beta }_0}{8N_c}\chi (n,\nu )\left(-\chi (n,\nu )+10/3+2ln\frac{{\mu }_R^2}{{\mu }_C^2}\right).$$

(9)

Here, ${\mu }_C\equiv \sqrt{m_{1\perp }{m}_{2\perp }}$, with $m_{(1,2)\perp }$ being the transverse masses of the two final-state objects. For the quarkonia, one has $m_{\mathcal{Q}\perp }=\sqrt{m_{\mathcal{Q}}^2+{\left|{\overrightarrow{p}}_{\mathcal{Q}}\right|}^2}$, with $m_{\mathcal{Q}}\equiv m_{J/\psi }=3.0969$ GeV or $m_{\mathcal{Q}}\equiv m_{\mathrm{{\rm Y}}}=9.4603$ GeV. The light-jet transverse mass coincides with its transverse momentum, $m_{\mathcal{Q}\perp }=\left|{\overrightarrow{p}}_J\right|$. The characteristic $\overline{\chi }(n,\nu )$ function was calculated in Ref. [364] and reads

$$\begin{gathered} \overline{\chi }(n,\nu )=-\frac{1}{4}\left\{\frac{{\pi }^2-4}{3}\chi (n,\nu )-6\zeta \left(3\right)-\frac{{\mathrm{d}}^2\chi }{\mathrm{d}{\nu }^2}+2\varphi (n,\nu )+2\varphi (n,-\nu )\right. \\ \left.+\frac{{\pi }^{2}sinh\left(\pi \nu \right)}{2\nu {cosh}^2\left(\pi \nu \right)}\left[\left(3+\left(1+\frac{n_f}{N_c^3}\right)\frac{11+12{\nu }^2}{16(1+{\nu }^2)}\right){\delta }_{n0}-\left(1+\frac{n_f}{N_c^3}\right)\frac{1+4{\nu }^2}{32(1+{\nu }^2)}{\delta }_{n2}\right]\right\}, \end{gathered}$$

(10)

with

$$\begin{gathered} \varphi (n,\nu )=-{\int }_0^1\mathrm{d}x\frac{x^{-1/2+i\nu +n/2}}{1+x}\left\{\frac{1}{2}\left({\psi }^{\prime }\left(\frac{n+1}{2}\right)-\zeta \left(2\right)\right)+{\mathrm{Li}}_2\left(x\right)+{\mathrm{Li}}_2(-x)\right. \\ \left.+lnx\left[\psi (n+1)-\psi \left(1\right)+ln(1+x)+\sum _{k=1}^{\infty }\frac{{(-x)}^k}{k+n}\right]+\sum _{k=1}^{\infty }\frac{x^k}{{(k+n)}^2}\left[1-{(-1)}^k\right]\right\} \\ =\sum _{k=0}^{\infty }\frac{{(-1)}^{k+1}}{k+(n+1)/2+i\nu }\left\{{\psi }^{\prime }(k+n+1)-{\psi }^{\prime }(k+1)\right. \\ \left.+{(-1)}^{k+1}\left[{\beta }_{\psi }(k+n+1)+{\beta }_{\psi }(k+1)\right]-\frac{\psi (k+n+1)-\psi (k+1)}{k+(n+1)/2+i\nu }\right\}, \end{gathered}$$

(11)

where

$${\beta }_{\psi }\left(z\right)=\frac{1}{4}\left[{\psi }^{\prime }\left(\frac{z+1}{2}\right)-{\psi }^{\prime }\left(\frac{z}{2}\right)\right],$$

(12)

and

$${Li}_2\left(x\right)=-{\int }_0^x\mathrm{d}y\frac{ln(1-y)}{y}.$$

(13)

The two quantities

$$c_{1,2}^{\mathrm{NLO}}(n,\nu ,|{\overrightarrow{p}}_{1,2}|,x_{1,2})=c_{1,2}(n,\nu ,|{\overrightarrow{p}}_{1,2}|,x_{1,2})+{\alpha }_s\left({\mu }_R\right){\widehat{c}}_{1,2}(n,\nu ,|{\overrightarrow{p}}_{1,2}|,x_{1,2})$$

(14)

represent the forward impact factors depicting the emission of a quarkonium and of a light jet. In the LO limit, we have

$$\begin{array}{cc}\hfill c_{\mathcal{Q}}(n,\nu ,|{\overrightarrow{p}}_{\mathcal{Q}}|,x_{\mathcal{Q}})& =2\sqrt{\frac{C_F}{C_A}}\left(\right|{\overrightarrow{p}}_{\mathcal{Q}}{{|}^2)}^{i\nu -1/2}{\int }_{x_{\mathcal{Q}}}^1\frac{\mathrm{d}z}{z}{\left(\frac{z}{x_{\mathcal{Q}}}\right)}^{2i\nu -1}\hfill \\ \hfill & \times \left[\frac{C_A}{C_F}{f}_g\left(z\right)D_g^{\mathcal{Q}}\left(\frac{x_{\mathcal{Q}}}{z}\right)+\sum _{a=q,\overline{q}}{f}_a\left(z\right)D_a^{\mathcal{Q}}\left(\frac{x_{\mathcal{Q}}}{z}\right)\right]\hfill \end{array}$$

(15)

and

$$c_J(n,\nu ,|{\overrightarrow{p}}_J|,x_J)=2\sqrt{\frac{C_F}{C_A}}\left(\right|{\overrightarrow{p}}_J{{|}^2)}^{i\nu -1/2}\left(\frac{C_A}{C_F}{f}_g\left(x_J\right)+\sum _{b=q,\overline{q}}{f}_b\left(x_J\right)\right),$$

(16)

respectively. Here, $C_F\equiv (N_c^2-1)/\left(2N_c\right)$ and $C_A\equiv N_c$ correspond to Casimir constants related to a gluon emission from a quark and a gluon, correspondingly. The $f\left(\nu \right)$ function contains the logarithmic derivative of LO impact factors

$$f\left(\nu \right)=\frac{i}{2}\frac{\mathrm{d}}{\mathrm{d}\nu }ln\left(\frac{c_1}{c_2^{*}}\right)+ln\left(|{\overrightarrow{p}}_1\left|\right|{\overrightarrow{p}}_2|\right).$$

(17)

What remains in Equation (6) is the NLO corrections to the forward impact factors, ${\widehat{c}}_{1,2}$. The NLO correction to the vector quarkonium impact factor is obtained in the light-quark limit [252] and its expression is given in Appendix A. This choice is fully consistent with our VFNS scheme, provided that the observed transverse momenta of the quarkonium are sufficiently larger than the mass of the constituent heavy quark. Our selection for the light-jet NLO impact factor is based on studies performed in Refs. [250,252], suited to numeric analyses, which embody a jet algorithm calculated in the “small-cone” limit (SCA) [365,366] with a cone-type selection [251]. Its analytic formula is reported in Appendix B.

We compare our NLL predictions with current LL predictions, which can be obtained by neglecting all the NLO pieces of the BFKL kernel and the impact factors. One obtains

$${\mathcal{C}}_n^{\mathrm{LL}/\mathrm{LO}}=\frac{e^{\Delta Y}}{s}{\int }_{-\infty }^{+\infty }\mathrm{d}\nu e^{\Delta Y{\overline{\alpha }}_s\left({\mu }_R\right)\chi (n,\nu )}{\alpha }_s^2\left({\mu }_R\right)c_1(n,\nu ,|{\overrightarrow{p}}_1|,x_1\left)\right[c_2(n,\nu ,|{\overrightarrow{p}}_2|,x_2{\left)\right]}^{*}.$$

(18)

Equations (6)–(18) highlight the way in which the hybrid factorization is constructed. The cross-section is high-energy factorized à la BFKL as a convolution between Green’s function and the two impact factors. Then, the latter contain the collinear inputs. The $\mathrm{NLL}/{\mathrm{NLO}}^{+}$ label highlights that the complete resummation of the energy logarithms is performed at NLL and within the NLO perturbative order. The ‘+’ superscript in Equation (6) indicates that some next-to-NLL contributions arise from the cross-product of the two NLO impact factor corrections.

We employ the formulæ presented in this Section at the natural scales provided by the process. In particular, we set ${\mu }_F={\mu }_R={\mu }_N\equiv m_{{\mathcal{Q}}_1\perp }+m_{{\mathcal{Q}}_2\perp }$, or ${\mu }_F={\mu }_R={\mu }_N\equiv m_{\mathcal{Q}\perp }+\left|{\overrightarrow{p}}_J\right|$, according to the considered final state (see Figure 1). Collinear PDFs are obtained from the NNPDF4.0 NLO parametrization [367,368] as implemented in LHAPDF v6.5.4 [369]. The NNPDF4.0 PDF determination is obtained from global fits and via the so-called replica method originally derived in Ref. [370] in the context of neural network techniques (for more details of the ambiguities arising from correlations between different PDF sets, see Ref. [371]). All results are obtained in the $\overline{\mathrm{MS}}$ renormalization scheme [363].

3. NRQCD Collinear Fragmentation for Vector Quarkonia

In Section 3.1, we give technical details regarding the way that we construct DGLAP-evolved, collinear FFs for vector $J/\psi$ and $\mathrm{{\rm Y}}$ states, starting with ${}^3{S}_1^{\left(1\right)}$ NRQCD inputs at the initial energy scale. In Section 3.2, we discuss the natural stabilization property emerging from the gluon fragmentation channel. We compare this pattern with corresponding patterns obtained for other hadron species.

3.1. DGLAP-Evolved NRQCD FFs

NLO collinear FFs describing the direct inclusive emission of a vector, $J/\psi$ or $\mathrm{{\rm Y}}$, can be constructed by starting from initial-scale NRQCD inputs. According to NRQCD factorization, the function describing the collinear fragmenting process of a parton a into a quarkonium state $\mathcal{Q}$ having longitudinal fraction z can be cast as

$$D_a^{\mathcal{Q}}(z,{\mu }_F)=\sum _{\left[n\right]}{\mathcal{D}}_a^{\mathcal{Q}}(z,{\mu }_F,\left[n\right])\langle {\mathcal{O}}^{\mathcal{Q}}\left(\left[n\right]\right)\rangle .$$

(19)

where ${\mathcal{D}}_a(z,{\mu }_F,\left[n\right])$ represents a short-distance, perturbative coefficient embodying DGLAP-type logarithms to be resummed through collinear evolution; $\langle {\mathcal{O}}^{\mathcal{Q}}\left(\left[n\right]\right)\rangle$ denotes a nonperturbative NRQCD LDME; and $\left[n\right]\equiv {}^{2S+1}{L}_J^{\left(\mathcal{C}\right)}$ denotes the comprehensive set of quarkonium quantum numbers in the spectroscopic notation [372], with the $\left(\mathcal{C}\right)$ superscript labeling the color state, singlet (1) or octet (8). In our study, we consider a spin-triplet state, namely a vector, taken in color singlet, so that ${}^{2S+1}{L}_J^{\left(\mathcal{C}\right)}\equiv {}^3{S}_1^{\left(1\right)}$.

The first building block is the heavy-quark fragmentation, $c\left(\overline{c}\right)\to J/\psi$ or $b\left(\overline{b}\right)\to \mathrm{{\rm Y}}$. The left plot of Figure 2 is one of the representative diagrams for this channel, with the green blob depicting the $\langle {\mathcal{O}}^{\mathcal{Q}}{(}^3{S}_1^{\left(1\right)})\rangle$ nonperturbative NRQCD LDME. Here, a standard color-triplet heavy quark, $Q\equiv c$ or $Q\equiv b$, is emitted from the hard subprocess. It radiates collinear gluons via DGLAP (not shown in the picture), thus losing energy. Then, when its energy is around $3m_Q$, it perturbatively fragments, according to NRQCD, into a $\left(Q\overline{Q}Q\right)$ state. At this stage, the $\left(Q\overline{Q}\right)$ subsystem, carrying color-singlet quantum numbers, nonperturbatively hadronizes into the observed vector quarkonium, $Q\equiv J/\psi$ or $Q\equiv \mathrm{{\rm Y}}$. The LO FF for the $Q\to \overline{Q}$ transition in color singlet was calculated in Ref. [77], whereas its NLO correction was recently obtained in Ref. [78]. We write 1

$$D_Q^{\mathcal{Q}}(z,{\mu }_F=3m_Q)=D_Q^{\mathcal{Q},\mathrm{LO}}\left(z\right)+\frac{{\alpha }_s^3({\mu }_R=3m_Q)}{m_Q^3}{\left|{\mathcal{R}}_{\mathcal{Q}}\left(0\right)\right|}^2{\mathrm{\Gamma }}_Q^{\mathcal{Q},\mathrm{NLO}}\left(z\right),$$

(20)

where $m_c=1.5$ GeV or $m_b=4.9$ GeV and $|{\mathcal{R}}_{J/\psi }{\left(0\right)|}^2=0.810$ GeV${}^3$ or $|{\mathcal{R}}_{\mathrm{{\rm Y}}}{\left(0\right)|}^2=6.477$ GeV${}^3$ is the NRQCD radial wave function at the origin of the heavy quarkonium, according to potential model studies [380]. The LO initial-scale function in Equation (20) reads

$$D_Q^{\mathcal{Q},\mathrm{LO}}\left(z\right)=\frac{{\alpha }_s^2\left(3m_Q\right)}{m_Q^3}\frac{8z{(1-z)}^2}{27\pi {(2-z)}^6}{\left|{\mathcal{R}}_{\mathcal{Q}}\left(0\right)\right|}^2(5z^4-32z^3+72z^2-32z+16),$$

(21)

while

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_Q^{J/\psi ,\mathrm{NLO}}\left(z\right)& =-9.01726z^{10}+18.22777z^9+16.11858z^8-82.54936z^7\hfill \\ & +106.57565z^6-72.30107z^5+28.85798z^4-6.70607z^3\hfill \\ & +0.84950z^2-0.05376z-0.00205\hfill \end{array}$$

(22)

or

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_Q^{\mathrm{{\rm Y}},\mathrm{NLO}}\left(z\right)& =-14.00334z^{10}+46.94869z^9-55.23509z^8+16.69070z^7\hfill \\ & +22.09895z^6-26.85003z^5+13.41858z^4-3.50293z^3\hfill \\ & +0.46758z^2-0.03099z-0.00226\hfill \end{array}$$

(23)

are the NLO FF core parts suitably cast in the form of fitted polynomial functions. From a perturbative QCD perspective, the LO heavy-quark FF starts at $\mathcal{O}\left({\alpha }_s^2\right)$, while its NLO correction starts at $\mathcal{O}\left({\alpha }_s^3\right)$.

The second building block is the gluon fragmentation, $g\to J/\psi$ or $g\to \mathrm{{\rm Y}}$, whose leading channel is portrayed in the right plot of Figure 2. Here, from the hard subprocess, a standard gluon is extracted. As for the heavy-quark case, it emits collinear radiation through DGLAP (not shown in the diagram), thus losing energy. Then, when its energy is around $2m_Q$, it perturbatively fragments into a $\left(Q\overline{Q}gg\right)$ state, according to NRQCD. The $\left(Q\overline{Q}\right)$ subsystem nonperturbatively hadronizes into the final-state quarkonium, $\mathcal{Q}$. We note that two extra gluon emissions are encoded in the NRQCD mechanisms. Indeed, being the parent gluon a colored object, in order to produce a colorless $\left(Q\overline{Q}\right)$ pair, at least another gluon needs to be emitted. However, due to the Landau–Yang selection rule [381,382], a spin-1 massive particle, such as our quarkonium, cannot be produced via a fermion triangle with two other on-shell vectors, such as gluons. Therefore, another secondary gluon must be emitted. One immediately recognizes that, from a perturbative QCD point of view, the leading gluon to vector quarkonium fragmentation channel is opened at $\mathcal{O}\left({\alpha }_s^3\right)$. Thus, to be consistent with the perturbative level of the heavy-quark FF, we do not need any higher-order correction for the gluon function. The FF for the $g\to \mathcal{Q}$ transition was calculated in Ref. [68] and reads 2

$$\begin{gathered} D_g^{\mathcal{Q}}(z,{\mu }_F=2m_Q)=\frac{5{\alpha }_s^3({\mu }_R=2m_Q)}{36{\left(2\pi \right)}^2}\frac{|{\mathcal{R}}_{\mathcal{Q}}{\left(0\right)|}^2}{m_{\mathcal{Q}}^3}{\int }_0^z\mathrm{d}\upsilon {\int }_{(1+z^2)/2z}^{(1+\upsilon )/2}\frac{\mathrm{d}\omega }{{(1-\omega )}^2{(\omega -\upsilon )}^2{({\omega }^2-\upsilon )}^2} \\ \sum _{l=0}^2{z}^l\left[f_l^{\left(g\right)}(\upsilon ,\omega )+g_l^{\left(g\right)}(\upsilon ,\omega )\frac{1+\upsilon -2\omega }{2(\omega -\upsilon )\sqrt{{\omega }^2-\upsilon }}ln\left(\frac{\omega -\upsilon +\sqrt{{\omega }^2-\upsilon }}{\omega -\upsilon -\sqrt{{\omega }^2-\upsilon }}\right)\right], \end{gathered}$$

(24)

where $f_{l=0,1,2}^{\left(g\right)}$ and $g_{l=0,1,2}^{\left(g\right)}$ are analytic functions calculated in Ref. [385]. We remark that, here, we consider the quarkonium’s direct production from the parent gluon. Another channel for high-energy $J/\psi$ emissions, not considered in our study, is the production of a P-wave charmonium ${\chi }_c$, followed by its radiative decay, given by the ${\chi }_c\to J/\psi +\gamma$ subprocess [77,386].

A major outcome of the phenomenological study of Ref. [80] is that the $c\to J/\psi$ and $g\to J/\psi$ fragmentation channels are similar in weight. Their relative size also depends on the hard scattering producing these partons. In particular, the number of gluons emitted at large transverse momentum might be of the same order or larger than the number of heavy quarks. Including the gluon channel is also relevant for our hadroproduction reactions (Figure 1), since the gluon FF is magnified by the collinear convolution with the corresponding gluon PDF, which is diagonal at LO (see Equation (15)).

Starting from the NRQCD inputs in Equations  (20)–(24), we generate a first and novel DGLAP-evolved set of collinear FFs for ${}^3{S}_1^{\left(1\right)}$ vector quarkonia. Several tools, such as QCD-PEGASUS [387], HOPPET [388], QCDNUM [389], APFEL(++) [390,391,392], and EKO [393], have been released as publicly available interfaces for DGLAP evolution. In contrast with collinear PDFs, whose evolution is space-like, the one for collinear FFs is time-like [69,70]. We employ APFEL++, in which the time-like evolution has been already implemented. According to the diagrams in Figure 2, the threshold for switching on the DGLAP evolution of gluons is set to ${\mu }_F=2m_Q$, whereas constituent heavy-quark and anti-quark FFs start to evolve at ${\mu }_F=3m_Q$. Light-quark and nonconstituent heavy-quark species are generated by evolution only. In this way, for each quarkonium species, $J/\psi$ and $\mathrm{{\rm Y}}$, we generate a phenomenology-ready FF set in LHAPDF format, named ZCW19${}^{+}$ [190], which encodes both the gluon NRQCD input and the constituent heavy-quark one.

3.2. Natural Stability

In this section, we shed light on the stabilization pattern arising from the fragmentation mechanism depicting the production of heavy-flavored hadrons in the VFNS. The connection between the behavior of DGLAP-evolving VFNS FFs and the stability of the $\mathrm{NLL}/{\mathrm{NLO}}^{(+)}$ hybrid high-energy and collinear factorization was discovered for the first time [323] in the context of singly heavy-flavored hadrons, such as D mesons [394,395,396,397,398,399,400,401,402], ${\Lambda }_c$ hyperons [395,403,404], and b-flavored (${\mathcal{H}}_b$) hadrons [405,406,407,408,409,410,411,412]. This remarkable property was first seen in observables sensible to the forward inclusive emission of ${\Lambda }_c$ [194,195,322], $D^{*\pm }$ [324,326], and ${\mathcal{H}}_b$ [195,321] particles. Then, further corroborating evidence emerged from the study of vector quarkonia [190,327] and charmed B mesons [325] produced via the NRQCD single-parton fragmentation mechanism [68,77,78,373,374,375,376,383,384]. Milder stabilization effects were also observed when s-flavored, cascade ${\mathrm{\Xi }}^{-}/{\overline{\mathrm{\Xi }}}^{+}$ baryons were tagged [313].

In the upper plots of Figure 3, we show the ${\mu }_F$-dependence of ZCW19${}^{+}$ collinear FFs describing the production of a $J/\psi$ (left) or a $\mathrm{{\rm Y}}$ (right). We fix the quarkonium longitudinal-momentum fraction z to a value that roughly corresponds to the average value at which FFs are probed in the kinematic ranges typical of our analyses (see Section 4.3). We have $z=0.45\simeq \langle z\rangle$. As expected, the constituent heavy (anti-)quark, $c\left(\overline{c}\right)$ for $J/\psi$ and $b\left(\overline{b}\right)$ for $\mathrm{{\rm Y}}$, heavily dominates over the other flavors. Notably, the gluon FF smoothly increases as ${\mu }_F$ grows, to reach a plateau. A similar trend characterizes the lower plots in Figure 3, where collinear ZCFW22 FFs depicting the emission of a $B_c^{\pm }{(}^1{S}_0)$ meson (left) or a $B_c^{*\pm }{(}^3{S}_1)$ one (right) are presented. The VFNS ZCFW22 sets for charmed B mesons were recently obtained in Ref. [325] by performing a DGLAP evolution of NQRCD initial-scale heavy-quark [376] and gluon [383] inputs at NLO. The strategy adopted there is analogous to the one defined in Ref. [190] for vector quarkonia and presented in Section 3 of this review.

The upper plots of Figure 4 contain the ${\mu }_F$ shape of FFs describing the production of two singly heavy-flavored hadrons, ${\Lambda }_c$ baryons (left) and ${\mathcal{H}}_b$ particles (right). The first FFs are depicted by the KKSS19 determination [403], a VFNS set based on a Bowler-like [413] initial-scale parametrization for charm and bottom (anti-)quark flavors. The second functions are portrayed by KKSS07 FFs [406,409], a VFNS set based on a power-like [414] initial-scale Ansatz for the bottom flavor. Moreover, in these two cases, heavy-quark FFs heavily prevail (for b-hadrons, only bottom functions dominate, since charm fragmentation is generated by DGLAP evolution), and the gluon FFs increase with ${\mu }_F$ (for b-hadrons, the growth is stronger). Finally, the lower plots of Figure 4 are for two species of light-flavored hadrons: $\Lambda$ hyperons (left) and $\mathrm{\Xi }$ baryons (right). Hyperons are described by AKK08 FFs [415], while cascade baryons are described by the novel SHKS22 neural network determination [416]. In the $\Lambda$ case, the strange FF dominates but, more importantly, all parton FF channels decrease with ${\mu }_F$. Cascade particle FFs present an intermediate situation, with the gluon channel smoothly growing up to roughly $20÷40$ GeV, and then slightly falling off.

As extensively explained in Refs. [190,194,321,323,325], the gluon collinear FF is a key actor in our $\mathrm{NLL}/{\mathrm{NLO}}^{(+)}$ hybrid factorization framework. Its energy dependence directly regulates the stability of the high-energy logarithmic series of our observables. Indeed, in the kinematic sectors of interest, i.e., when ${10}^{-4}\lesssim x\lesssim {10}^{-2}$, the gluon PDF strongly prevails over all (anti-)quark channels. Since the gluon FF is diagonally convoluted with the gluon PDF in the LO hadron impact factor (see Equation (15)), its behavior is strongly heightened. As discussed in Ref. [194], this feature is valid also at NLO, namely where the $\left(qg\right)$ and $\left(gq\right)$ nondiagonal channels are switched on (see Appendix A).

On one hand, the QCD running coupling becomes smaller when ${\mu }_R$ increases, and this happens both in Green’s function and in the impact factors (see Section 2.2). On the other, it is well known that the gluon PDF grows with ${\mu }_F$. When the latter is convoluted in the impact factor with a gluon FF that also grows with ${\mu }_F$, as happens for heavy-flavored species, the two effects compensate each other. This gives rise to the stability of heavy-hadron distributions under scale variation. The stronger the growth of the gluon FF with ${\mu }_F$, the clearer the stabilization pattern. Thus, ${\mathcal{H}}_b$ sensitive distributions are even more stable than ${\Lambda }_c$ distributions [321]. Conversely, when the gluon FF falls off with ${\mu }_F$, as happens for $\Lambda$ hyperons, no stabilization under scale variation is found. This hampers any possibility of performing precision studies of the high-energy distributions at natural scales [192]. Cascade $\mathrm{\Xi }$ baryons present an intermediate situation, with the smooth-behaved ${\mu }_F$ gluon FF leading to a partial and milder stabilization pattern [313].

The rise of the natural stability in the presence of both singly heavy-flavored hadrons and quarkonia gives strong evidence that this remarkable property is an intrinsic feature of heavy-flavor emissions, which becomes manifest whenever a heavy-hadron species is detected, independently of the Ansätze made for the corresponding collinear FFs.

4. Phenomenology

All the presented results were obtained with Jethad, a Python/Fortran hybrid and multimodular interface aimed at the calculation, management, and processing of physical observables defined in the context of different formalisms [192,195]. In particular, numeric computations of rapidity distributions were performed through the Fortran 2008 modular routines implemented in Jethad, while the Jethad Python 3.0 analysis environment was employed for the elaboration of the results. Section 4.1 encompasses the core features of the current version of the Jethad technology (v0.5.0, not yet public). Details of our strategy to assess the weights of uncertainties for high-energy predictions are given in Section 4.2. The final-state kinematic configurations to which our distributions are tailored can be found in Section 4.3. Numeric analyses and a discussion of the results for the rapidity distributions are presented in Section 4.4.

4.1. Jethadv0.5.0: A Basic Overview

The origins of the Jethad project date back to the end of 2017, when the constantly increasing number of semi-hard hadron [291,293] or/and jet [284,286,309] sensitive LHC final states, proposed as a testing ground for high-energy resummation, called for the development of a standard and reference interface, suited both to the computation and the elaboration of high-energy distributions. The first (private) recognized version, v0.2.7, served as a useful tool to perform a pioneering comparison between $\mathrm{NLL}/\mathrm{NLO}$ hybrid predictions and fixed-order computations for hadron–jet correlations in the high-energy limit of QCD [192]. The possibility of considering forward heavy-quark pair emissions both in photo- and hadroproduction was implemented in v0.3.0 [268]. Studies of Higgs emissions as well as transverse-momentum distributions became accessible with v0.4.2 [314]. The first version of the Python analyzer was combined with the Fortran core in v0.4.3 [193]. The first analyses of heavy-flavored hadrons via VFNS FFs at NLO started with v0.4.4 [194]. The ΔΥναμις (DYnamis) work package, suitable to investigate forward Drell–Yan dilepton emissions [159], was incorporated in Jethad v0.4.5. The possibility of accessing the proton content at low x through the UGD and gluon TMDs become possible after interfacing the previously stand-alone Leptonic-Exclusive-Amplitudes (LExA) modular code [154] with Jethad v0.4.6 [156]. Support for quarkonium studies from NRQCD single-parton fragmentation started with v0.4.7 [190]. The newest features in v0.5.0 essentially rely on (i) an improved system to perform studies of energy-scale variations, (ii) the extension of the list of observables that can be analyzed, with particular attention to singly and doubly differential transverse momentum distributions [290,313], and (iii) the possibility of removing the NLO expanded contribution from the fully NLL resummed series to support matching procedures with collinear factorization (see Equation (1) of Ref. [360]).

From the main core to the service modules and routines, Jethad has been designed to dynamically reach a high level of computational performance. Jethad’s multidimensional integrators make use of extensive parallel computing to actively select the best integration algorithm, depending on the shape of the integrand. Any process implemented in Jethad can be actively selected via a user-friendly, structure-based smart management interface, where physical final-state particles are portrayed by object prototypes (code structures in the working environment). Particle objects encode all the information about the basic and kinematic properties of their physical Doppelgängers, from mass and charge to kinematic ranges and rapidity tags. As a first step, they are loaded from a master database via a dedicated particle generation routine (custom particle generation is also possible). Subsequently, they are cloned into a final-state object array and therefore injected from the integrand routine, differential regarding final-state kinematic variables, to the respective hard-factor module via a suitable controller. The strong flexibility characterizing the generation of the physical final states is complemented by a window of possibilities in selecting the initial state. Indeed, a peculiar particle ascendancy structure attribute allows Jethad to dynamically recognize whether a particle is hadroproduced, leptoproduced, photoemitted, and so on. In this way, only the relevant modules will be initialized, thus breaking down computing time lags. Jethad is structured as an object-based interface, entirely independent of the reaction being investigated. While inspired by high-energy resummation and TMD factorization phenomenology, different approaches can be easily encoded in it by implementing new and customized modules and supermodules, which can be linked to the core structure of the code via a native point-to-routine system, thus making Jethad a general, particle physics environment.

With the aim of providing the scientific community with a standard software program designed for the computation and management of different types of processes described by distinct formalisms, we plan, as a medium-term goal, to release the first public version of Jethad.

4.2. Uncertainty Estimation

Accurate predictions for our high-energy observables rely upon identifying the main potential uncertainty sources and quantifying their effects. A widely employed strategy to assess the size of higher-order corrections, not accounted for by our formalism, consists in gauging the sensitivity of our distribution on renormalization- and factorization-scale variation. From a recent analysis on semi-inclusive emissions of light and heavy bound states [195], it emerged that the effect of varying energy scales around their natural values represents one of the main contributions to the global uncertainty. Therefore, we study the effect of simultaneously varying ${\mu }_R$ and ${\mu }_F$ around their natural values, in the range of 1/2 to 2. The $C_{\mu }$ parameter in the plots in Section 4.4 refers to the ratio $C_{\mu }={\mu }_R/{\mu }_N$. We also set ${\mu }_F={\mu }_R$. Another potential source of relevant uncertainty might come from proton PDFs. However, recent studies on semi-hard observables pointed out that the choices of different PDF parametrizations, as well as of different members inside the same set, brings no sizeable impact [192,195,309,321]. Thus, we use the central member of only one PDF determination, the NNPDF4.0 one. Two further potential sources of uncertainty might be connected to $a)$ a collinear improvement in the NLO BFKL Green’s function [417,418,419,420,421,422,423], which is based on the inclusion of renormalization group (RG) contributions needed to impose a compatibility with the DGLAP equation in the collinear limit, and from $b)$ employing another renormalization scheme. The point $a)$ was studied in detail in Ref. [195]. It was found that the effect of the collinear improvement on rapidity distributions was completely contained by the error bands produced by energy-scale variations, and it was even less relevant when other other azimuthal-angle sensitive observables were considered. In the same work, an estimation of point $b)$, namely the upper limit of the effect of moving from $\overline{\mathrm{MS}}$ to the MOM renormalization scheme, was provided. MOM predictions for rapidity distributions are systematically larger than those for $\overline{\mathrm{MS}}$ predictions, but are still affected by error bands due to scale variations. We remark that a consistent MOM study would rely upon MOM-evolved PDFs, not available so far. In view of these considerations, we generate uncertainty bands for our predictions by considering the net and combined effects of varying energy scales, together with the numeric uncertainty generated by final-state multidimensional integration (Section 4.3). The latter is steadily kept below $1\%$ by the Jethad numeric integrators.

4.3. Final-State Ranges

In our phenomenological study, we provide $\mathrm{NLL}/{\mathrm{NLO}}^{+}$ predictions for the rapidity distribution, i.e., the cross-section differential in the rapidity distance, $\Delta Y$, between the two outgoing objects. Its expression follows from the integration of the ${\mathcal{C}}_0$ azimuthal angle coefficient, defined in Section 2.2, over the final-state rapidities and transverse momenta, with $\Delta Y$ being taken to be fixed. One has

$$\frac{\mathrm{d}\sigma (\Delta Y,s)}{\mathrm{d}\Delta Y}=\int \mathrm{d}|{\overrightarrow{p}}_1|\int \mathrm{d}|{\overrightarrow{p}}_2|{\int }_{y_1^{\min },max(\Delta Y+y_2^{\min })}^{y_1^{\max },min(\Delta Y+y_2^{\max })}\mathrm{d}{y}_1{\mathcal{C}}_0\left(|{\overrightarrow{p}}_1|,|{\overrightarrow{p}}_2|,y_1,y_2\right){|}_{y_2=y_1-\Delta Y}$$

(25)

Here, a $\delta ((y_1-y_2)-\Delta Y)$ function has been used to remove and to enforce the extremes of integration in $y_1$ accordingly. As proposed in Ref. [190], the transverse momenta of the detected quarkonia lie in the range 20 GeV to 60 GeV. This choice is perfectly consistent with our VFNS treatment, which requires the energy scales to be sufficiently larger than the thresholds controlling the DGLAP evolution of heavy quarks. Conversely, light-flavored jets are reconstructed with transverse-momentum cuts typical of current works at the CMS detector [424], 35 GeV $<|{\overrightarrow{p}}_J|\equiv |{\overrightarrow{p}}_2|<$ 60 GeV. Adopting asymmetric transverse-momentum windows allows us to better disengage the core high-energy signal from the DGLAP background [192,284,285]. It also suppresses Sudakov logarithmic contributions rising from almost back-to-back configurations, which would call for another appropriate resummation mechanism [425,426,427,428,429]. Moreover, it dampens possible instabilities emerging in next-to-leading computation [430,431], as well as energy–momentum conservation violations [432].

Concerning rapidity ranges, we select two relevant cases.

• Standard kinematic configurations, characteristic of the ongoing LHC phenomenology, where the quarkonium is identified in the detector barrel only [433], while the jet is reconstructed also in the endcap calorimeters. We have $|y_{{\mathcal{Q}}_{(1,2)}}|<2.4$ and $|y_J|<4.7$, so that $\Delta Y<4.8$ in the double quarkonium channel and $\Delta Y<7.1$ in the quarkonium plus jet channel.
• Extended kinematic configurations, suitable to test the validity of the natural stability at the edges of application of the hybrid factorization, with the quarkonium detected both in the barrel and the endcaps; we have $|y_{{\mathcal{Q}}_{(1,2)}}|<4.7$ and $|y_J|<4.7$, so that $\Delta Y<9.4$ in the two production channels.

4.4. $\Delta Y$-Distribution

Figure 5 shows the behavior of the $\Delta Y$-distributions for the double vector quarkonium (left) and the inclusive vector quarkonium plus jet production (right) at $\mathrm{NLL}/{\mathrm{NLO}}^{+}$, for $\sqrt{s}=14$ TeV, and in standard rapidity configurations. The values of our distributions are all larger than ${10}^{-2}$ pb in the vector quarkonium plus jet channel (right). This is a very favorable statistic, which turns out to be lower than the one for heavy-baryon and heavy-light meson emissions [194,321], but substantially higher than the one for the double vector quarkonium channel (left). The falloff of these high-energy resummed cross-sections as $\Delta Y$ increases is a common pattern of all the hadronic semi-hard processes analyzed so far. It is the net effect of the interplay of two competing features. On the one hand, the pure high-energy dynamics would lead to the well-known growth with the energy of partonic hard factors. On the other hand, collinear DGLAP-evolving PDFs and FFs quench hadronic distributions when $\Delta Y$ increases.

In particular, the plots of Figure 5 contain a study of the $\Delta Y$-distributions under the progressive variation of the ${\mu }_R$ and ${\mu }_F$ scales in a wide range, regulated by the $C_{\mu }$ parameter, which runs from 1 to 30. The ancillary panels below the primary plots exhibit the reduced cross-sections, i.e., divided by their central values taken at $C_{\mu }=1$. Notably, the sensitivity of the $\mathrm{NLL}/{\mathrm{NLO}}^{+}$ predictions on $C_{\mu }$ variation is rather weak for all four considered channels (Figure 1) and for all the values of the probed $\Delta Y$ spectrum. This indicates that our VFNS quarkonium FFs, built in terms of both initial-scale heavy-quark and gluon NRQCD functions, and then evolved via DGLAP, lead to a strong stabilization pattern, as with the one observed for singly heavy-flavored hadrons and for charmed B mesons (see results on the same distributions given in Refs. [190,194,321,325]).

Figure 6 shows the behavior of the $\Delta Y$-distributions for the double vector quarkonium (left) and the inclusive vector quarkonium plus jet production (right) at $\mathrm{NLL}/{\mathrm{NLO}}^{+}$, for $\sqrt{s}=14$ TeV, and in extended rapidity configurations. Here, as anticipated, we allow the vector quarkonium to be tagged not only in the barrel calorimeter but also in the endcaps, to match the detection range of the light jet. On the one hand, from an experimental point of view, such configurations could be difficult to realize, due to hardware limitations in particle identification at high rapidities. On the other hand, however, this extension offers us a unique opportunity to test, from a phenomenological perspective, the validity of the natural stability at the edges of the applicability of our hybrid factorization. Indeed, at a very large rapidity separation between the two detected objects, we enter the so-called threshold regions. In this case, forward final-state objects stem from a parton struck with a large longitudinal-momentum fraction x. As a result, the weight of undetected gluon radiation diminishes and, correspondingly, large Sudakov-type double logarithms (threshold double logarithms) appear in the perturbative series. In principle, they could spoil the convergence of the perturbative series, and they should be accounted for to all orders by a proper resummation mechanism [434,435,436,437,438,439,440,441,442,443,444,445,446,447,448,449,450,451], which is different from our high-energy one. At the same time, large $\Delta Y$ values probe our distribution in a high-energy regime, which might be “asymptotic”, namely where the BFKL dynamics might be the dominant mechanism or, at least, it could strongly modulate the pure fixed-order background.

The global pattern emerging from the inspection of the plots in Figure 6 unambiguously highlights that the $\Delta Y$-distributions are very stable under progressive energy-scale variation, in all the investigated channels. The natural stability is not spoiled by the extended rapidity kinematics, but rather magnified. Finally, we observe that the statistic is lower than that obtained in the standard rapidity range, but still favorable for the quarkonium plus jet channel.

We complete the analysis given in this section by providing predictions for $\mathrm{NLL}/{\mathrm{NLO}}^{+}$$\Delta Y$-differential cross-sections as compared with the corresponding ones taken in a pure $\mathrm{LL}/\mathrm{LO}$ limit. Results for the four channels under investigation (see Figure 1) after applying standard and extended final-state cuts are, respectively, presented in Figure 7 and Figure 8. The shaded uncertainty bands embody the combined effects of energy-scale variations and numeric multidimensional integrations (see Section 4.2 for details). Additionally, in this case, the lower ancillary panels refer to reduced distributions, namely divided by central values taken at $C_{\mu }=1$. The overall pattern is one of fair stability under $\mathrm{NLL}/\mathrm{NLO}$ corrections, with the $\mathrm{NLL}/{\mathrm{NLO}}^{+}$ bands generally narrower and partially nested inside the corresponding $\mathrm{LL}/\mathrm{LO}$ ones. Their overlap is weaker in the quarkonium plus jet channel. This is expected, since having only one heavy hadron in the final state, instead of two, clearly halves the stabilizing power coming from the VFNS collinear FFs.

We conclude this section by mentioning a relevant aspect that goes beyond the aim of this review, but still deserves a discussion. We refer to the impact of MPI effects on the differential distributions at increasing values of $\Delta Y$. In particular, we consider the Double-Parton Scattering (DPS), which accounts for two hard-scattering subreactions at most. Ref. [283] deals with a study of DPS corrections to Mueller–Navelet jets, which can be potentially relevant for $\Delta Y$-distributions at large center-of-mass energies and moderate transverse momenta. From investigations of double $J/\psi$ [106,111], $J/\psi$ plus $\mathrm{{\rm Y}}$ [106], $J/\psi$ plus Z [105], and $J/\psi$ plus W [452] final states, an important indication emerged that the DPS contribution might be sizeable and it should be considered to improve the description of the moderate transverse-momentum spectrum. Then, MPI imprints could arise also in triple $J/\psi$ events [453,454]. Therefore, the search for MPI signatures in our observables represents an important development that should be carried out in the medium-term future.

5. Conclusions and Outlook

We investigated the direct inclusive hadroproduction of a forward vector quarkonium, $J/\psi$ or $\mathrm{{\rm Y}}$, in association with another backward quarkonium or a light-flavored jet in a hybrid high-energy and collinear factorization formalism. We relied upon the single-parton fragmentation approximation, which represents the most adequate tool to describe quarkonium hadronization at large transverse momentum [68,80,81,455]. By starting from a NRQCD input at the initial energy scale for the constituent heavy-quark [78] and the gluon [68] fragmenting into the observed meson, we achieved the first and novel determination of VFNS-consistent, DGLAP-evolving quarkonium FFs, which we named ZCW19${}^{+}$. We believe that this FF set can serve as a useful tool for further phenomenological applications lying outside the high-energy domain, as well as a support basis for forthcoming analyses aimed at unveiling the connection between collinear factorization and NRQCD effective theory.

Inspired by recent findings on the inclusive forward emissions of singly heavy-flavored hadrons studied in the context of the $\mathrm{NLL}/{\mathrm{NLO}}^{(+)}$ hybrid factorization [194,321], we searched for a stabilization pattern in the high-energy series under higher-order corrections and energy-scale variation. Our key finding is that these effects are present also when vector quarkonia are detected at large transverse momenta and in forward rapidity directions. Notably, the natural stability of the $\Delta Y$-distributions turns out to be even more manifest when meson tags are tailored in larger rapidity windows, such the ones reachable when endcap detection is allowed. Thus, we provide strong evidence that this remarkable property is an intrinsic feature characterizing heavy-flavor emissions, which becomes manifest whenever a heavy-hadron species is detected, independently of the initial-scale input starting from which corresponding collinear FFs are obtained by DGLAP evolution. Following the results presented in this review, we can unambiguously state that the natural stability is so strong as to surmount and overpower any possible destabilizing effect coming from other resummation mechanisms, particularly from threshold contaminations highly affecting the forward emission domain [434,435,436,437,438,439,440,441,442,443,444,445,446,447,448,449,450,451].

Our analysis of vector quarkonia from NRQCD single-parton fragmentation constitutes an important step forward in the study of heavy-flavored emissions at $\mathrm{NLL}/{\mathrm{NLO}}^{(+)}$, starting from the analytic computation of heavy-quark pair impact factors [266,267,268]. It also represents the first and most relevant contact point with the high-energy quarkonium phenomenology. Future investigations are needed to address quarkonium production in broader domains, such as the regimes reachable at new-generation colliding machines [456,457,458,459,460,461,462,463,464,465,466,467,468,469,470,471,472,473,474,475,476,477,478,479,480,481,482,483].

Another important advancement comes from NLO studies of single-forward or almost back-to-back semi-inclusive emissions in the context of the gluon saturation framework (see, e.g., Refs. [484,485,486,487,488,489,490,491,492,493,494,495] and references therein). In this context, the impact of soft gluon radiation on angular asymmetries in di-jet production was addressed in Refs. [496,497,498,499,500]. A recent NLO computation [501] will make it possible to perform similar analyses also for diffractive di-hadron detections. Notably, NLO saturation is relevant to access the (un)polarized gluon content of protons and nucleons at small x [502,503,504,505,506,507]. The authors of Refs. [508,509,510,511,512] investigated quarkonium production in proton–proton and proton–nucleus collisions by combining both NRQCD low-x effects. Future studies will explore the common ground between our method to address quarkonium production via the $\mathrm{NLL}/{\mathrm{NLO}}^{(+)}$ hybrid factorization and higher-order calculations for NLO saturation in the exclusive emissions of heavy mesons [513,514].

Finally, in order to reach a higher level of precision in our knowledge of quarkonium fragmentation, we must compare NRQCD-inspired collinear FFs with the determinations of these quantities obtained by extraction via global fits. On this basis, neural network techniques already developed for the extraction of light-flavored hadron FFs [416,515,516,517,518,519,520,521] will represent an asset. Another important development will be linking the Jethad technology to fully automated codes devoted to the study of quarkonium production at NLO. We mainly refer to MadGraph5_aMC@NLO [522] and HELAC-Onia [116,117] as interfaced with the NLOAccess virtual platform [523,524].