Diffractive semi-hard production of a $J/\psi$ or a $\Upsilon$ from single-parton fragmentation plus a jet in hybrid factorization

Abstract

We investigate the inclusive hadroproduction of a heavy quarkonium ($J/\psi$ or $\Upsilon$), in association with a light-flavored jet, as a testing ground for the semi-hard regime of QCD. Our theoretical setup is the hybrid high-energy and collinear factorization, where the standard collinear approach is supplemented by the t-channel resummation of leading and next-to-leading energy logarithms à la BFKL. We present predictions for rapidity and azimuthal-angle differential distributions, hunting for stabilizing effects of the high-energy series under higher-order corrections. Our reaction represents an additional channel to test the production mechanisms of quarkonia at high energies and large transverse momenta and to possibly shed light on the transition region from heavy-quark pair production to single-parton fragmentation of $J/\psi$ and $\Upsilon$ states.

1 Hors d’œuvre

The study of production mechanisms of heavy-quark bound states in high-energy collisions has always played a key role in widening our knowledge of the dynamics of strong interactions. A relevant class of heavy-flavored hadrons is represented by mesons whose lowest Fock state is composed of a heavy quark and the corresponding antiquark. These mesons are called heavy quarkonia, or simply quarkonia. The era of quarkonium physics began with the November Revolution, when on 11 November 1974 the discovery of a new particle with a mass of about 3.1 GeV and the same quantum numbers of the photon was announced. That particle is the $J/\psi$ vector meson and owes its name to the two groups that simultaneously observed it, namely Collaboration leaded by B. Richter [1] at the Stanford Linear Accelerator Center (SLAC) and one headed by S. Ting [2] at the Brookhaven National Laboratory (BNL). Seven days after the SLAC and BNL announcements, the $J/\psi$ discovery was confirmed by the ADONE experiment in Frascati, under the direction of G. Bellettini [3]. The hadronic nature of the $J/\psi$ was proven by the cross-section ratio between hadrons and $\mu ^+\mu ^-$ in inclusive $e^+e^-$ inclusive annihilations that turned out to be much higher at the resonance, this meaning that $J/\psi$-like particles directly decay into hadronic states. This provided us with a proof of the existence of the charm-quark flavor and, more in general, with a clear evidence that quarks are real constituents and not just mathematical artifacts. The existence of the charm quark was first proposed in 1964 by Bjørken and Glashow [4], but there were few experimental evidences. In 1970, with the proposal of the Glashow–Iliopoulos–Maiani (GIM) mechanism [5], the existence of a new quark species, the charm one, became a necessary ingredient to explain why the flavor-changing neutral currents (FCNCs) are suppressed in loop diagrams.

Quarkonium studies represented one of the first testing grounds for the manifestation of fundamental properties of Quantum ChromoDynamics (QCD), such as the asymptotic freedom. Indeed, with the discovery of the $\psi (2S)$, the first radial excited state of the $J/\psi$, it became clear that strong interactions decrease at short distances, resulting in a Coulomb-like binding potential with a confinement part. The $J/\psi$ was the first charmonium state to be discovered. The name charmonium was given in analogy with a $e^+e^-$ bound pair, that is a positronium, and represents a family of hadrons which comprehends all the $c {\bar{c}}$ meson states. In the years following the November Revolution other charmonia were observed (as the P-wave $\chi _c$ and the pseudoscalar $\eta _c$), as well as open-charm D mesons [6]. In 1977 the first bottomonium state was observed. It was called $\Upsilon$, and it is a vector meson equivalent to the $J/\psi$, but with bottom quarks [7]. It was followed by the discovery of excited $b {\bar{b}}$ states (as the $\Upsilon (2S)$) and of open-bottom B mesons [8].

Despite the ease of performing experimental studies on quarkonia, especially on the vector mesons which can be easily observed in lepton–lepton reactions, the phenomenological description of their production rates still remains a challenge and different models for the quarkonium production mechanism have been proposed so far (see Refs. [9, 10] for further details). An early-day model, known as Color Evaporation Model (CEM) [11, 12], postulates that the quarkonium-constituent $(Q {\bar{Q}})$-pair undergoes a large number of soft interactions that make it a color singlet hadronizing into the observed quarkonium. Soft interactions lead to a complete color decorrelation between the produced partons and the pair, so that even a single gluon can evolve into the final hadron. The main limit of the CEM is the inability to afford a different description for production rates of distinct quarkonia, as the ones of $J/\psi$ and $\chi _c$ [13, 14].

A second attempt at describing quarkonium production relies on the so-called Color Singlet Mechanism (CSM) [15,16,17,18]. It is based on the assumption, opposite to the CEM one, that the $(Q {\bar{Q}})$-pair does not evolve due to the suppression of gluon emissions with powers of $\alpha _s(m_Q)$, where $m_Q$ is the heavy-quark mass. Since the quark spin and color remain unaltered during the hadronization, and since the quarkonium must be color-neutral, the $(Q {\bar{Q}})$-pair needs to be generated in a color-singlet state. Furthermore, since the quarkonium mass, $m_{{\mathcal {Q}}}$, turns out to be slightly larger than $2 m_Q$, the two constituent quarks can be assumed to be almost at rest in the hadron frame, namely with zero relative velocity, $v_{{\mathcal {Q}}}$, and the only non-perturbative contribution to the production mechanism is a non-relativistic Schrödinger wave function at the origin, ${\mathcal {R}}_{{\mathcal {Q}}}(0)$. The strength of the CSM consists in its relatively high predictive power. Indeed, the value of the radial wave function at the origin is the only free parameter which can be extracted from measurements of quarkonium leptonic-decay widths. However, the CSM suffers from infrared singularities arising at the next-to-leading order (NLO) in P-wave decay channels [19, 20].

To go beyond the CSM approximation of $(Q {\bar{Q}})$ static quantum numbers, it was proposed that higher Fock states (e.g., $|Q {\bar{Q}} g \rangle$), contribute to the quarkonium production and might remove the NLO singularity.In a $|Q {\bar{Q}} g \rangle$ state the quark-pair is in a color-octet (CO) configuration. More generally, the physical quarkonium is given as a linear combination of all possible Fock states. All contributions are organized as a double expansion in powers of $\alpha _s$ and $v_{{\mathcal {Q}}}$, of which the CSM represents the leading term in $v_{{\mathcal {Q}}}$ for S-wave channels (for higher waves this contribution may be zero). This double expansion is the building block of an effective field theory (EFT) known as Non-Relativistic QCD (NRQCD) [21,22,23]. Within this approach cross sections for quarkonium productions take the form of a sum of partonic cross sections creating a given Fock state, each of them being multiplied by a Long-Distance Matrix Element (LDME) encoding the non-perturbative hadronization. Although affording a very appealing solution to the quarkonium production puzzle, the NRQCD introduces many new non-perturbative parameters that makes non-trivial the comparison with data.

Another relevant issue comes from the fact that a quarkonium state can be produced via the decay or the de-excitation of a heavier hadron. As an example, a $J/\psi$ can be produced from the decay of B meson [24]. This channel is non-prompt, since there is a visible time delay due to the decay. Some experimental studies have managed to isolate the non-prompt contribution by detecting the distance between the creation vertex of the B particle and its decay or de-excitation to the quarkonium state. This procedure permits us to get the prompt component. Within the prompt channel, a charmed quarkonium, say the $J/\psi$, can be still generated from the decay of a $\chi _c$, or from the de-excitation of a $\psi (2S)$. This is known as indirect production component. Conversely, the fraction of hard-scattering produced quarkonia is known as direct component and is generally not accessible alone, but can be accessed by subtraction only.

The models described above rely on the assumption that the dominant production mechanisms is the short-distance production of a $(Q {\bar{Q}})$ pair in the hard scattering which then hadronizes into the detected quarkonium. Since the quark and the antiquark are created with a relative transverse separation of order $1/p_T$, this channel is expected to be suppressed at large $p_T$. According to the short-distance mechanism, the $(Q {\bar{Q}})$ pair is created with a transverse separation of order 1/Q, with Q the characteristic energy scale of the process [25]. In the large transverse-momentum range one has $Q \sim p_T \gg m_Q$. Thus, exchanges of particles having virtualities of the same order of $p_T$ quenches the probability amplitude, since they are related to scattering of particles in a quite small volume, $1/p_T^3$ [26, 27].

In Ref. [28] it was highlighted that in this kinematic regime an additional mechanism comes into play, namely the fragmentation of a large-$p_T$ parton (gluon or quark) which afterwards decays into the quarkonium state plus other partons. The hadronization process is described by a set of collinear Fragmentation Functions (FFs) which evolve according to the Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) equation [29,30,31,32,33]. NRQCD provides us a way to perturbatively calculate these FFs. Here, the $(Q {\bar{Q}})$ pair is produced with a separation of order $1/m_Q$. Although the fragmentation mechanism is often of higher perturbative order with respect to the short-distance one, it is enhanced by $(Q/m_Q)^2$. Therefore, it prevails at high energies, $Q \gg m_Q$. Figure 1 shows a selection of leading Feynman diagrams contributing to vector-quarkonium fragmentation. The NRQCD calculation of the gluon FF to a S-wave quarkonium in the CSM was completed in Ref. [28] for $J/\psi$ and $\eta _c$ at leading order (LO) and in Ref. [34] for $\eta _c$ at NLO. The corresponding results for the FF depicting the $c \rightarrow J/\psi +c+X$ were obtained in Ref. [35] at LO and in Ref. [36] at NLO. LO phenomenological studies on unveiling the transition region between the short-distance and the fragmentation mechanisms were performed in Refs. [37,38,39,40,41]. FFs from a heavy-quark pair in the S- and P-wave channels were investigated in Refs. [42, 43], respectively. In those works next-to-leading power studies in $p_T$ which have shown that the leading-power single-parton fragmentation is relevant at larger $p_T$-values, rather than what initially found in the 1990s. A next-to-NLO factorization analysis on quarkonium fragmentation functions can be found in Refs. [44, 45]. FFs for polarized quarkonia were studied in Refs. [46,47,48,49,50,51] (see Ref. [52] for a discussion). Reference [53] contains analyses on quarkonium factorization at evolution at large transverse momenta.

Despite the ambiguities in the description of their production mechanisms and in the correct interpretation of experimental data, quarkonium emissions are considered excellent tools to investigate properties of strong interactions. First studies of quarkonium production in collinear factorization and with NLO accuracy can be found, e.g., in Refs. [54,55,56,57,58,59,60,61,62,63]. References [64, 65] contain NLO analyses on the transverse-momentum spectrum within CEM. Investigations on quarkonium polarization in NRQCD were done in Refs. [66,67,68]. Multi-Parton Interactions (MPI) effects in $J/\psi$ production were studied in Refs. [69, 70]. The impact of the nuclear modification of the gluon density on quarkonium was recently assessed [71,72,73]. Advances in the automation of NRQCD calculations for quarkonium emissions are reported in Refs. [74,75,76,77]. Issues in the large-$p_T$ description of $J/\psi$ production at NLO were attributed to an over-subtraction in the factorization the collinear divergences inside collinear Parton Distribution Functions (PDFs) [78,79,80], and a possible solution was found in the inclusion of high-energy effects on top of NLO calculations [81]. The nature of hidden-charm and bottom tetra- and penta-quarks was investigated [82, 83] via the pioneering hadro-quarkonium picture [84, 85]. Efforts in obtaining tetra-quark FFs from NRQCD-inspired calculations have been recently made [86, 87].

Semi-inclusive quarkonium-detection processes are widely recognized as golden channels to access the three-dimensional structure of hadrons via transverse-momentum-dependent gluon distributions [88,89,90,91]. Recent phenomenological analyses of gluon TMDs through observables sensitive to quarkonium production in electron–proton and proton–proton collisions were can be found in Refs. [92,93,94,95,96,97,98,99,100], respectively.¹ Studies on TMD factorization for heavy-quarkonium emissions were recently conducted in a Soft and Collinear Effective Theory (SCET) [102, 103], and then extended to di-jet and heavy-meson pair emissions in lepton-proton scatterings [104]. A SCET approach was employed to access the jet substructure and, more in particular, to define fragmenting jet functions, needed to describe the production of quarkonia inside jets [105,106,107,108,109]. The $p_T$-spectrum of quarkonium production from single light-parton fragmentation was recently investigated at the hands of SCET and NRQCD [110]. Phenomenological analyses of quarkonium emissions at new-generation colliders via a novel determination of polarized gluon TMDs [111,112,113,114,115,116,117] are underway.

Quarkonium studies at low-x are relevant to observe the growth with energy of cross sections [118] predicted by the linear Balitsky–Fadin–Kuraev–Lipatov (BFKL) evolution [119,120,121,122], as well as the onset of possible the non-linear gluon-recombination effects [123] leading to the saturation pattern. In Ref. [124] the emission of quarkonia in high-energy proton–nucleus scatterings was studied by making use of NRQCD and by accounting for small-x evolution and multiple-scattering effects within the Color Glass Condensate (CGC) formalism, a semi-classical EFT for saturation (see Refs. [125, 126] for a comprehensive overview). A selection of recent results on forward and central quarkonium production in (ultra-peripheral) hadron and lepton-hadron collisions can be found in Refs. [127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143]. Combining the information coming from single-forward emissions of quarkonia with the one encoded in single-forward detections of light vector mesons [144,145,146,147,148,149,150,151,152,153] and in the forward Drell–Yan di-lepton reaction [154,155,156] will permits us to perform an extensive scan of the intersection kinematic range between TMD and high-energy dynamics. In particular, the main outcome of Refs. [146, 151] is that the forward polarized $\rho$-meson leptoproduction at HERA and the EIC act as an excellent probe channel for the proton content at small-x, the corresponding $p_T$-factorized cross sections being sensitive to a small transverse-size dipole-scattering mechanism. On the other side, forward tags of excited states, such as $\psi (2S)$ and $\Upsilon (2S)$ mesons, could lead to quite larger sizes of dipoles [157,158,159], this calling for a TMD-inspired enhancement of the pure small-x description.

In this work we turn our attention to the inclusive hadroproduction of heavy vector mesons at high-energies, as the ones reachable at the Large Hadron Collider (LHC) and at new-generation hadron accelerators. This permits us to access kinematic regions where $\sqrt{s} \gg m_{{\mathcal {Q}}}$. The limit in which $\sqrt{s} \gg \{Q \} \gg \Lambda _\mathrm{QCD}$² is known as Regge–Gribov or semi-hard regime [160, 161]. As it is well known, in this regime large logarithms of the ratio $s/Q^2$ appear to all orders the perturbative series with a power increasing with the $\alpha _s$ order, thus spoiling its convergence. This simply means that the standard fixed-order approach for the computation of hard scattering cross sections is no longer a reliable framework to build theoretical predictions up. The need for a resummation to all orders that takes into account the effect of these logarithms, led authors of Refs. [119,120,121,122] to develop the previously mentioned BFKL approach, which is now established as the most powerful formalism for the description of the semi-hard QCD sector. Its validity holds within the leading approximation (LLA), which means resumming all terms proportional to $\alpha _s^n \ln (s/Q^2)^n$, and within the next-to-leading approximation (NLA), which means including all terms proportional to $\alpha _s^{n+1} \ln (s/Q^2)^n$.

When the observed particles are widely separated in rapidity, cross sections for hadronic processes take a peculiar factorized form, given as the convolution of two impact factors, related to the fragmentation of each colliding particle to an identified final-state object, and a process-independent Green’s function, which encodes the resummation of energy logarithms. The BFKL Green’s function satisfies an integral evolution equation, whose kernel is known up to NLO for any fixed, not growing with s, momentum transfer t and for any possible two-gluon color state in the t-channel [162,163,164,165,166,167,168]. Impact factors instead depend on processes, so they representing the most challenging part of the calculation. So far, the number of impact factors calculated within NLO accuracy is quite small.

Evidence of the onset of the BFKL dynamics was provided in the analysis of differential distributions for the hadroproduction of two objects widely separated in rapidity. An incomplete list of them includes: the inclusive detection of two light jets featuring large transverse momenta and well separated in rapidity (Mueller–Navelet channel [169]), for which several phenomenological studies have appeared so far (see, e.g., Refs. [170,171,172,173,174,175,176,177,178,179,180,181,182,183,184]), the inclusive emission of a light di-hadron system [185,186,187,188,189], multi-jet tags [190,191,192,193,194,195,196,197,198], hadron-plus-jet [199,200,201,202], Higgs-plus-jet [203,204,205,206,207,208,209], heavy-light di-jet system [210, 211], Drell–Yan-plus-jet [212], heavy-flavored hadrons [213,214,215,216,217,218], and notably $J/\psi$-plus-jet production [219].

A major problem rising in phenomenological studies of semi-hard processes via the high-energy resummation is connected to the fact that NLA corrections both to the BFKL Green’s function and impact factors have almost the same size and opposite sign of pure LLA terms. This generates instabilities in the BFKL series unstable which become strongly manifest when renormalization and factorization scales are varied from their natural values, namely the ones dictated by kinematics. Although employing some scale-optimization method, such as the Brodsky–Lepage–Mackenzie (BLM) procedure [220,221,222,223] effectively helped to reduce these instabilities in reactions featuring the emissions of light jets and/or hadrons [173, 175, 185, 187], it turned out that the given optimal values for scales were by far higher than the natural ones [224]. This led to a lowering of cross sections of one or more orders of magnitude, thus hampering any possibility of doing precision studies.

A first clue of a reached stability of semi-hard observables at natural scales was observed quite recently in reactions featuring the detection of particles with a large transverse mass, as Higgs bosons [203, 225] (see Refs. [226, 227] for novel NLO calculations) and heavy-quark jets [210]. Due to the lack of NLO calculations for the corresponding impact factors, these analyses were performed with a partial NLA accuracy. A strong evidence of stabilization effects in full NLA calculation emerged in a recent study on $\Lambda _c$-baryon production [216] In particular, it was pointed out that the peculiar pattern of VFNS FFs depicting the production of the charmed baryon acts as a fair stabilizer of the high-energy series. This result was subsequently confirmed also for bottomed hadrons [217].

Fig. 1

Left: One of the leading diagrams contributing to the heavy-quark fragmentation to a $^3S_1^{(1)}$ vector quarkonium ${\mathcal {Q}}$ at order $\alpha _s^2$. Right: One of the leading diagrams contributing to the gluon fragmentation to $^3S_1^{(1)}$ vector quarkonium ${\mathcal {Q}}$ at order $\alpha _s^3$. The green blob denotes the corresponding non-perturbative NRQCD LDME

With the aim of hunting for the aforementioned stabilizing effects in semi-hard emissions of heavy-quarkonium states, in this article we investigate the high-energy behavior of the inclusive hadroproduction of a $J/\psi$ or a $\Upsilon$ accompanied by a light-quark jet at the LHC. The two detected objects are widely separated in rapidity. We remark that only the direct $J/\psi$ and $\Upsilon$ production components are considered in our study. Both the meson and the jet feature large transverse momenta, and they are well separated in rapidity. Kinematic typical of current LHC analyses feature moderate values of partons’ longitudinal-momentum fractions.³ This justify a description in terms of collinear PDFs. From a purely collinear-factorization viewpoint, the large required transverse momenta makes valid using a variable-flavor number-scheme (VFNS) approach [228, 229], in which the cross section for the production of a light parton is constructed and then convoluted with a FF that describes the transition from light quarks to heavy bound states ($J/\psi$ or $\Upsilon$ in our case). Then, making use of NRQCD, the quarkonium FF is constructed as a product between short-distance coefficient functions, which encode the resummation of DGLAP-type logarithms, and the corresponding LDMEs. We will employ a recent NLO determination for heavy-quark to $J/\psi$ and $\Upsilon$ FFs [36], together with the corresponding gluon to $J/\psi$ and $\Upsilon$ functions [28]. From a high-energy perspective, large rapidity intervals in the final state lead to non-negligible exchanges of transverse momenta in the t-channel which call for a $p_T$-factorization procedure, naturally provided by the BFKL resummation within NLA accuracy. Combining all these ingredient makes a full NLO treatment feasible. Our hybrid high-energy and collinear factorization, already set as the reference formalism for the description of inclusive two-particle semi-hard emissions, is thus enriched with a new element, the NRQCD, which provides us with a useful model for quarkonium fragmentation. Another formalism, close in spirit with our hybrid factorization and suited for single forward detections, was proposed in Refs. [230, 231].

As already mentioned, the fragmentation mechanism becomes increasingly competitive as well as the transverse momentum of tagged particle grows. We remind, however, that since we are neglecting the heavy-quark mass at the level of the hard part, our formalism is not suited to properly describe the intermediate region, where $|{\vec p}_T| \sim m_Q$. Here powers of the $\sqrt{m_Q^2 + |{\vec p}_T|^2}/|{\vec p}_T|$ ratio need to be taken into account. In this sense, our study is complementary to the one proposed in Ref. [219] (see also preliminary results in Refs. [232, 233]), where the inclusive semi-hard $J/\psi$-plus-jet production was investigated by making use of the short-distance $(Q {\bar{Q}})$ mechanism for the quarkonium production. In that work massive impact factors for the transitions $(g \rightarrow Q \bar{Q} \rightarrow J/\psi )$ and $(g \rightarrow Q \bar{Q} g \rightarrow J/\psi + g)$, encoding the corresponding LDMEs, where computed and then used to consider both color-singlet and color-octet configurations of the meson.⁴ In such a calculation all powers of $\sqrt{m_Q^2 + |{\vec p}_T|^2}/|{\vec p}_T|$ ratio are present, but since the final state features no DGLAP evolution, logarithms are not resummed. In the future, a possible matching between the two approaches could help us to enhance the description of heavy-flavor in the context of the high-energy resummation.

Coming back to our large-$p_T$ fragmentation treatment, possible ambiguities could emerge. On the one side, a pure collinear-factorization description relies on the fragmentation of a single parton to the identified hadron. At large transverse momenta the parton is always treated as a light one and the hadronization mechanism is portrayed by zero-mass (ZM) VFNS FFs that evolve with DGLAP. On the other side, a NRQCD FF always contains the perturbative splitting of the parton produced via the hard scattering to a massive $(Q {\bar{Q}})$ pair. Synthesizing these two viewpoints could hide some complications which are currently matter of investigation [236]. We leave these insights to future studies, and take the NRQCD functions of Ref. [36] as proxy models for $J/\psi$ and $\Upsilon$ collinear FFs.

Furthermore, one might argue that our analysis is not complete due to the absence of the CO contribution. Indeed, Ref. [36] contains the FF of the constituent heavy quark, Q, to a color-singlet quarkonium state, ${\mathcal {Q}}$, namely the contribution proportional to the $^3S_1^{(1)}$ LDME (see Sect. 2.2 for more details), and the same holds for the gluon to quarkonium FF [28]. We stress, however, that our primary goal is hunting for stabilizing effects of high-energy resummed distributions for $J/\psi$ and $\Upsilon$ production under higher-order corrections and energy-scale variations. Therefore we rely on the CSM contribution only. The inclusion of higher Fock-state contributions in our single-parton fragmentation description is postponed to a forthcoming analysis.

The structure of this article reads as follows. In Sect. 2 we introduce our theoretical setup, presenting the highlights of our hybrid high-energy and collinear factorization (Sect. 2.1), as well as details on the heavy-meson fragmentation mechanism and the light-jet selection function (Sect. 2.2). In Sect. 3, we discuss our phenomenological study of differential distributions. Finally (Sect. 4) we come out with a summary and future challenges.

2 Theoretical setup

In this section we present theoretical ingredients to build our observables. First we give analytic expressions of azimuthal-angle coefficients for the considered reaction (Fig. 2), calculated by means of our hybrid high-energy and collinear factorization (Sect. 2.1). Then we discuss our choice for quarkonium FFs and for the jet algorithm (Sect. 2.2).

Fig. 2

Hybrid high-energy/collinear factorization for the inclusive quarkonium + jet production. Red blobs refer to proton collinear PDFs, green blobs stand for quarkonium single-parton FFs, while blue arrows denote a light-flavored jet emission. The BFKL ladder, given by the yellow blob, is connected to impact factors via Reggeon (zigzag) lines. Diagrams were done by making use of the JaxoDraw 2.0 interface [237]

2.1 High-energy resummed cross section

The process under investigation is

$$\begin{aligned} p(P_a) + p(P_b) \rightarrow {{\mathcal {Q}}}(p_{{\mathcal {Q}}}, y_{{\mathcal {Q}}}) + X + {\mathrm{jet}}(p_J, y_J), \end{aligned}$$

(1)

where $p(P_{a,b})$ stands for an initial proton with momentum $P_{a,b}$, ${{\mathcal {Q}}}(p_{{\mathcal {Q}}}, y_{{\mathcal {Q}}})$ is a $J/\psi$ or a $\Upsilon$ emitted with momentum $p_{{\mathcal {Q}}}$ and rapidity $y_{{\mathcal {Q}}}$, the light jet is tagged with momentum $p_J$ and rapidity $y_J$, and X denotes all the undetected products of the reaction. High observed transverse momenta, $|\vec p_{{{\mathcal {Q}}},J}|$, together with a large rapidity separation, $\Delta Y= y_{{\mathcal {Q}}}- y_J$, are required conditions to get a diffractive semi-hard configuration in the final state. Furthermore the transverse-momentum ranges need to be enough large to ensure the validity of description of the quarkonium production mechanism in terms of single-parton VFNS collinear fragmentation.

The momenta of the two colliding protons are taken as Sudakov vectors satisfying $P_a^2= P_b^2=0$ and $2 (P_a\cdot P_b) = s$, and the four-momenta of final-state particles can be decomposed as

$$\begin{aligned} p_{{{\mathcal {Q}}},J}= & {} x_{{{\mathcal {Q}}},J} P_{a,b} + \frac{\vec p_{{{\mathcal {Q}}},J}^{\,2}}{x_{{{\mathcal {Q}}},J} s}P_{b,a} + p_{{{\mathcal {Q}}},J\perp },\nonumber \\&\quad p_{{{\mathcal {Q}}},J\perp }^2=-\vec p_{{{\mathcal {Q}}},J}^{\,2}, \end{aligned}$$

(2)

where the outgoing particle longitudinal momentum fractions, $x_{{{\mathcal {Q}}},J}$, are connected to the corresponding rapidities via the relation $y_{{{\mathcal {Q}}},J}=\pm \frac{1}{2}\ln \frac{x_{{{\mathcal {Q}}},J}^2 s}{\vec p_{{{\mathcal {Q}}},J}^2}$. One has $\mathrm{d}y_{{{\mathcal {Q}}},J} = \pm \frac{\mathrm{d}x_{{{\mathcal {Q}}},J}}{x_{{{\mathcal {Q}}},J}}$, and $\Delta Y\equiv y_{{\mathcal {Q}}}- y_J = \ln \frac{x_{{\mathcal {Q}}}x_J s}{|\vec p_{{\mathcal {Q}}}||\vec p_J|}$.

In a genuine QCD collinear treatment, the LO cross section of our process (Eq. (1)) is cast as a convolution of the partonic hard subprocess with the parent-proton PDFs and the quarkonium FFs

$$\begin{aligned}&\frac{\mathrm{d}\sigma ^\mathrm{LO}_\mathrm{coll.}}{\mathrm{d}x_{{\mathcal {Q}}}\mathrm{d}x_J\mathrm{d}^2\vec p_{{\mathcal {Q}}}\mathrm{d}^2\vec p_J}=\sum _{\alpha ,\beta =q,{{\bar{q}}},g}\int _0^1 \mathrm{d}x_a \int _0^1 \mathrm{d}x_b\ f_\alpha \left( x_a\right) f_\beta \left( x_b\right) \nonumber \\&\quad \times \int _{x_{{\mathcal {Q}}}}^1 \frac{\mathrm{d}z}{z}D^{{{\mathcal {Q}}}}_{\alpha }\left( \frac{x_{{\mathcal {Q}}}}{z}\right) \frac{\mathrm{d}{\hat{\sigma }}_{\alpha ,\beta }\left( {\hat{s}}\right) }{\mathrm{d}x_{{\mathcal {Q}}}\mathrm{d}x_J\mathrm{d}^2\vec p_{{\mathcal {Q}}}\mathrm{d}^2\vec p_J}. \end{aligned}$$

(3)

Here the $\alpha ,\beta$ indices indicate the parton species (quarks $q = u, d, s, c, b$; antiquarks ${\bar{q}} = {\bar{u}}, {\bar{d}}, {\bar{s}}, {\bar{c}}, {\bar{b}}$; or gluon g), $f_{\alpha ,\beta }\left( x, \mu _F \right)$ are the colliding-proton PDFs and $D^{{{\mathcal {Q}}}}_{\alpha }\left( x/z, \mu _F \right)$ denote the outgoing-quarkonium FFs; $x_{a,b}$ are the longitudinal-momentum fractions of the partons initiating the hard subprocess and z the longitudinal fraction of the single parton that fragments into ${{\mathcal {Q}}}$. Finally, $\mathrm{d}\hat{\sigma }_{\alpha ,\beta }\left( {\hat{s}} \right)$ stands for the partonic cross section, with ${\hat{s}} \equiv x_a x_b s$ the squared center-of-mass energy of the partonic collision. For the sake of brevity, the explicit dependence of PDFs, FFs and partonic cross section on the factorization scale, $\mu _F$, has been dropped.

At variance with collinear factorization, in order to obtain the expression for the high-energy resummed cross section in our hybrid framework, one needs first to consider the high-energy factorization which is naturally encoded in the BFKL formalism, and then enhance the description by embodying collinear PDFs and FFs. We start by writing the differential cross section as a Fourier sum of the so-called azimuthal-angle coefficients

$$\begin{aligned} \frac{\mathrm{d}\sigma }{\mathrm{d}y_{{\mathcal {Q}}}\mathrm{d}y_J \mathrm{d}\vec p_{{\mathcal {Q}}}\mathrm{d}\vec p_J \mathrm{d}\phi _{{\mathcal {Q}}}\mathrm{d}\phi _J} = \frac{1}{(2\pi )^2} \left[ {\mathcal {C}}_0 + 2 \sum _{n=1}^\infty \cos (n \varphi )\, {\mathcal {C}}_n \right] , \nonumber \\ \end{aligned}$$

(4)

with $\varphi _{{{\mathcal {Q}}},J}$ the azimuthal angles of the detected particles and $\varphi = \phi _{{\mathcal {Q}}}- \phi _J - \pi$. The azimuthal coefficients ${\mathcal {C}}_n \equiv {\mathcal {C}}_n^{\mathrm{NLA}}$ can be calculated in the BFKL approach and they encode the resummation of energy logarithms up to the NLA accuracy. A NLA formula obtained in the $\overline{\mathrm{MS}}$ renormalization scheme reads (for details on the derivation see, e.g., Ref. [171])

$$\begin{aligned} {\mathcal {C}}_n^{\mathrm{NLA}}= & {} \int _0^{2\pi } \mathrm{d}\phi _{{\mathcal {Q}}}\int _0^{2\pi } \mathrm{d}\phi _J\, \cos (n \varphi ) \, \frac{\mathrm{d}\sigma _\mathrm{NLA}}{\mathrm{d}y_{{\mathcal {Q}}}\mathrm{d}y_J\, \mathrm{d}|\vec p_{{\mathcal {Q}}}| \, \mathrm{d}|\vec p_J| \mathrm{d}\phi _{{\mathcal {Q}}}\mathrm{d}\phi _J}\; \nonumber \\= & {} \frac{e^{\Delta Y}}{s} \int _{-\infty }^{+\infty } \mathrm{d}\nu \, e^{{\Delta Y} \bar{\alpha }_s(\mu _R)\left\{ \chi (n,\nu )+\bar{\alpha }_s(\mu _R) \left[ \bar{\chi }(n,\nu )+\frac{\beta _0}{8 N_c}\chi (n,\nu )\left[ -\chi (n,\nu )+\frac{10}{3}+4\ln \left( \frac{\mu _R}{\sqrt{|\vec p_{{\mathcal {Q}}}| |\vec p_J|}}\right) \right] \right] \right\} } \nonumber \\&\times \, \alpha _s^2(\mu _R) \, \left[ c_{{\mathcal {Q}}}^\mathrm{NLO}(n,\nu ,|\vec p_{{\mathcal {Q}}}|, x_1)[c_J^\mathrm{NLO}(n,\nu ,|\vec p_J|,x_2)]^*\,+ \bar{\alpha }_s^2(\mu _R) \, \Delta Y\frac{\beta _0}{4 N_c}\chi (n,\nu )f(\nu ) \right] , \end{aligned}$$

(5)

with $\bar{\alpha }_s(\mu _R) \equiv \alpha _s(\mu _R) N_c/\pi$, $N_c$ the color number and $\beta _0 = 11N_c/3 - 2 n_f/3$ the first coefficient of the QCD $\beta$-function. A two-loop running-coupling setup with $\alpha _s\left( M_Z\right) =0.11707$ and a dynamic flavor number $n_f$ is employed. The expression for the LO BFKL characteristic function reads

$$\begin{aligned} \chi \left( n,\nu \right) =2\left\{ \psi \left( 1\right) -\mathrm{Re}\left[ \psi \left( i\nu +\frac{1}{2}+\frac{n}{2} \right) \right] \right\} \end{aligned}$$

(6)

where $\psi (z) = \Gamma ^\prime (z)/\Gamma (z)$ is the logarithmic derivative of the Gamma function. The $\bar{\chi }(n,\nu )$ function in Eq. (5) contains NLO corrections to the BFKL kernel and was calculated in Ref. [238] (see also Ref. [239]). The two functions

$$\begin{aligned}&c_{{{\mathcal {Q}}},J}^\mathrm{NLO}(n,\nu ,|\vec p\,|,x) \nonumber \\&\quad = c_{{{\mathcal {Q}}},J}(n,\nu ,|\vec p\,|,x) + \alpha _s(\mu _R) \, {\hat{c}}_{{{\mathcal {Q}}},J}(n,\nu ,|\vec p\,|,x) \end{aligned}$$

(7)

are the impact factors for the production of a heavy-quarkonium state and for the emission of a light jet. Their LO parts read

$$\begin{aligned}&c_{{\mathcal {Q}}}(n,\nu ,|\vec p\,|,x) = 2 \sqrt{\frac{C_F}{C_A}} (|\vec p\,|^2)^{i\nu -1/2}\,\int _{x}^1\frac{\mathrm{d}\zeta }{\zeta } \left( \frac{\zeta }{x} \right) ^{2 i\nu -1} \nonumber \\&\quad \times \left[ \frac{C_A}{C_F}f_g(\zeta )D_g^{{\mathcal {Q}}}\left( \frac{x}{\zeta }\right) +\sum _{\alpha =q,{\bar{q}}}f_\alpha (\zeta )D_\alpha ^{{\mathcal {Q}}}\left( \frac{x}{\zeta }\right) \right] \end{aligned}$$

(8)

and

$$\begin{aligned}&c_J(n,\nu ,|\vec p\,|,x) = 2 \sqrt{\frac{C_F}{C_A}} (|\vec p\,|^2)^{i\nu -1/2}\nonumber \\&\quad \times \,\left( \frac{C_A}{C_F}f_g(x) +\sum _{\beta =q,{\bar{q}}}f_\beta (x)\right) , \end{aligned}$$

(9)

respectively. The $f(\nu )$ function is defined in terms of the logarithmic derivative of LO impact factors

$$\begin{aligned} f(\nu ) = \frac{i}{2} \, \frac{\mathrm{d}}{\mathrm{d}\nu } \ln \left( \frac{c_{{\mathcal {Q}}}}{c_J^*}\right) + \ln \left( |\vec p_{{\mathcal {Q}}}| |\vec p_J|\right) . \end{aligned}$$

(10)

The remaining functions in Eq. (5) are the NLO impact-factor corrections, ${\hat{c}}_{{{\mathcal {Q}}},J}$. The NLO correction to the ${{\mathcal {Q}}}$ impact factor is calculated in the light-quark limit [240]. This option is fully consistent with our VFNS treatment, provided that the $p_{{\mathcal {Q}}}$ values at work are much larger than the heavy-quark mass(charm or bottom, see Sect. 3 for more details). Our choice for the jet NLO impact factor is discussed at the end of Sect. 2.2.

We present for completeness the expression for the azimuthal coefficients in the pure LLA approximation, namely by neglecting NLO terms in the BFKL kernel and impact factors

$$\begin{aligned} {\mathcal {C}}_n^{\mathrm{LLA}}= & {} \frac{e^{\Delta Y}}{s} \int _{-\infty }^{+\infty }\mathrm{d}\nu \, e^{{\Delta Y} \bar{\alpha }_s(\mu _R)\chi (n,\nu )} \, \alpha _s^2(\mu _R) \, c_{{\mathcal {Q}}}\nonumber \\&\times \,(n,\nu ,|\vec p_{{\mathcal {Q}}}|, x_1)[c_J(n,\nu ,|\vec p_J|,x_2)]^* . \end{aligned}$$

(11)

In Eqs. (5)–(11) it is highlighted the way how our hybrid factorization is realized. The cross section is factorized à la BFKL as a convolution between the Green’s function and impact factors, the latter ones encoding the collinear functions.

We use expressions given in this section at the natural scales provided by the process. In particular we set $\mu _F = \mu _R = \mu _N \equiv \sqrt{m_{{{\mathcal {Q}}}\perp } |\vec p_J|}$, with $m_{{{\mathcal {Q}}}\perp } = \sqrt{ m_{{\mathcal {Q}}}^2 + |\vec p_{{\mathcal {Q}}}|^2}$ the quarkonium transverse mass. We fix $m_{{\mathcal {Q}}}\equiv m_{J/\psi } = 3.0969$ GeV or $m_{{\mathcal {Q}}}\equiv m_{\Upsilon } = 9.4603$ GeV, according to the emitted meson. Collinear PDFs are obtained via the MMHT14 NLO parameterization [241] as provided by the LHAPDFv6.2.1 library [242]. All calculations of our azimuthal coefficients are done in the $\overline{\mathrm{MS}}$ renormalization scheme.

As a supplementary analysis, we will present values of energy scales obtained by applying the BLM optimization method [220,221,222,223]. We will not calculate our observables at these scales, but we will compare these last with the ones typical of inclusive emissions of other heavy-flavored hadrons, thus corroborating the statement that heavy-flavor collinear FFs act as stabilizers of the high-energy resummation.

The BLM method prescribes that the optimal renormalization-scale value, $\mu _R^\mathrm{BLM}$, is the one that allows us to remove the $\beta _0$ dependence of a given observable. In Ref. [177] a general procedure to apply the BLM prescription on BFKL-resummed azimuthal coefficients was set up. The condition for the BLM scale setting on a given coefficient, $C_n$, is the solution of the following integral equation

$$\begin{aligned} C_n^{(\beta )}(s, \Delta Y) = \int \mathrm{d}\Phi (y_{{{\mathcal {Q}}},J}, |\vec p_{{{\mathcal {Q}}},J}|, \Delta Y) \, \, {\mathcal {C}}_n^{(\beta )} = 0, \end{aligned}$$

(12)

where $\mathrm{d}\Phi (y_{{{\mathcal {Q}}},J}, |\vec p_{{{\mathcal {Q}}},J}|, \Delta Y)$ is the final-state differential phase space (see Sect. 3),

(13)

and

$$\begin{aligned} {\omega }(\nu ) = \frac{f(\nu )}{2} - \frac{1}{6} - \frac{2}{3} I + \ln \left( \frac{\mu ^\mathrm{BLM}_R}{\sqrt{|\vec p_{{\mathcal {Q}}}| |\vec p_J|}} \right) . \end{aligned}$$

(14)

In Eq. (13) $\bar{\alpha }_s^\mathrm{{MOM}} = N_c/\pi \, \alpha _s^\mathrm{{MOM}}$ stands for the expression of the strong coupling in the momentum renormalization scheme (MOM), given by inverting the relation

$$\begin{aligned} \alpha _s^{\overline{\mathrm{MS}}} = \alpha _s^\mathrm{MOM} \left[ 1 + (\tau ^{\beta } + \tau ^\mathrm{conf}) \frac{\alpha _s^\mathrm{MOM}}{\pi } \right] , \end{aligned}$$

(15)

with

$$\begin{aligned} \tau ^\beta = - \left( \frac{1}{2} + \frac{I}{3} \right) \beta _0 \end{aligned}$$

(16)

and

$$\begin{aligned} \tau ^\mathrm{conf}= \frac{N_c}{8}\left[ \frac{17}{2}I +\frac{3}{2}\left( I-1\right) \xi +\left( 1-\frac{1}{3}I\right) \xi ^2-\frac{1}{6}\xi ^3 \right] , \nonumber \\ \end{aligned}$$

(17)

where $I=-2\int _0^1\mathrm{d}z \frac{\ln ~z}{z^2-z+1} \simeq 2.3439$ and $\xi$ is gauge parameter fixed at zero in the following. We write the optimal renormalization scale as $\mu _R^\mathrm{BLM} = C_\mu ^\mathrm{BLM} \mu _N$, with $\mu _N$ the process natural scale defined previously, and we look for values of the $C_\mu ^\mathrm{BLM}$ parameter which solve Eq. (12).

2.2 Quarkonium fragmentation and jet selection function

We build our NLO collinear FF sets for the direct $J/\psi$ or $\Upsilon$ meson production by taking, as a starting point, the recent work done in Ref. [36]. There, a NLO calculation was performed for the heavy-quark FF depicting the transition $c \rightarrow J/\psi$ or the $b \rightarrow \Upsilon$ one, where c (b) indistinctly refer to the charm (bottom) quark and its antiquark. It essentially relies on the NRQCD factorization formalism taken with NLO accuracy (see, e.g., Refs. [14, 22, 23, 243,244,245] and references therein), which allows us to write the FF function of a parton i fragmenting into a heavy quarkonium ${\mathcal {Q}}$ with longitudinal fraction z as⁵

$$\begin{aligned} D^{{\mathcal {Q}}}_i(z, \mu _F) = \sum _{[n]} {\mathcal D}^{{\mathcal {Q}}}_{i}(z, \mu _F, [n]) \langle {\mathcal {O}}^{{\mathcal {Q}}}([n]) \rangle . \end{aligned}$$

(18)

In Eq. (18), ${\mathcal D}_{i}(z, \mu _F, [n])$ denotes the perturbative short-distance coefficient containing terms proportional to $\ln (\mu _F/m_{{\mathcal {Q}}})$ (to be resummed via DGLAP evolution), $\langle {\mathcal {O}}^{{\mathcal {Q}}}([n]) \rangle$ stands for the non-perturbative NRQCD LDME, and $[n] \equiv \,^{2S+1}L_J^{(c)}$ represents the quarkonium quantum numbers in the spectroscopic notation (see, e.g., Ref. [251]), the (c) superscript identifying the color state, singlet (1) or octet (8). Limiting ourselves to a spin-triplet (vector) and color-singlet quarkonium state, $^3S_1^{(1)}$, the analytic form of the initial-scale FF depicting the constituent heavy-quark to quarkonium transition, $Q \rightarrow {{\mathcal {Q}}}$ (we inclusively refer to $c \rightarrow J/\psi$ or $b \rightarrow \Upsilon$), reads (for details on its derivation, see Sections II and III of Ref. [36])

$$\begin{aligned}&D^{{\mathcal {Q}}}_Q(z, \mu _F \equiv \mu _0) = D^{{\mathcal {Q}}, \mathrm{LO}}_Q(z) \nonumber \\&\quad + \frac{\alpha _s^3(3m_Q)}{m_Q^3} \, |{\mathcal {R}}_{{\mathcal {Q}}}(0)|^2 \, \Gamma _Q^{{\mathcal {Q}}, \mathrm{NLO}}(z) , \end{aligned}$$

(19)

with $m_c = 1.5$ GeV or $m_b = 4.9$ GeV, and the NRQCD radial wave-function at the origin of the quarkonium state set to $|{\mathcal {R}}_{J/\psi }(0)|^2 = 0.810$ GeV$^3$ or to $|{\mathcal {R}}_{\Upsilon }(0)|^2 = 6.477$ GeV$^3$, according to potential-model calculations (Ref. [252] and references therein). The expression for the LO initial-scale FF was originally calculated in Ref. [35] and reads

$$\begin{aligned}&D^{{\mathcal {Q}}, \mathrm{LO}}_Q(z) = \frac{\alpha _s^2(3m_Q)}{m_Q^3} \,\frac{8z(1-z)^2}{27\pi (2-z)^6} \, |{\mathcal {R}}_{{\mathcal {Q}}}(0)|^2\nonumber \\&\quad \times \, (5z^4 - 32z^3 + 72z^2 -32z + 16) , \end{aligned}$$

(20)

and the polynomial function $\Gamma _Q^{{\mathcal {Q}}, \mathrm{NLO}}(z)$ entering the expression for the NLO-FF correction is

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

(21)

or

$$\begin{aligned} \Gamma _Q^{\Upsilon , \mathrm{NLO}}(z)= & {} - 14.00334z^{10} + 46.94869z^9 \nonumber \\&- 55.23509z^8 + 16.69070z^7 \nonumber \\&+ 22.09895z^6 - 26.85003z^5 \nonumber \\&+ 13.41858z^4 - 3.50293z^3 \nonumber \\&+ 0.46758z^2 - 0.03099z - 0.00226. \end{aligned}$$

(22)

Coefficients of z-powers in Eqs. (21) and (22) are obtained via a polynomial fit to the numerically-calculated NLO FFs. Starting from $\mu _F \equiv \mu _0 = 3 m_Q$, in Ref. [36] a DGLAP-evolved formula for the $D^{{\mathcal {Q}}}_Q(z, \mu _F)$ function was derived and then applied to phenomenological studies of $J/\psi$ and $\Upsilon$ production via $e^+ e^-$ single inclusive annihilation (SIA).

As pointed out in Ref. [38], both $(c \rightarrow J/\psi )$ and $(g \rightarrow J/\psi )$ fragmentation channels are similar in size. The relative weight of the heavy-quark and gluon contributions is also driven by the size of the hard scattering producing these partons. Therefore, the number of large $p_T$-gluons emitted could be of the same order, if not larger, than the heavy-quark one. Moreover, in a hadroproduction process such as the one considered in our study, the gluon FF is enhanced by the collinear convolution at LO with the corresponding gluon PDF (see Eq. (8)). Thus we expect, in our case, a stronger sensitivity on the gluon-fragmentation channel with respect to the case of a SIA-like reaction. Therefore, we include in our analysis also the contribution coming from the gluon fragmentation. The gluon to vector-quarkonium LO fragmentation mechanism start at $\alpha _s^3$, namely at the same order of the heavy-quark FF in Eq. (19). The $(g \rightarrow {^3S_1^{(1)}} g g)$ fragmentation function was computed in Ref. [28] and reads

$$\begin{aligned}&D_g^{{\mathcal {Q}}} (z, 2 m_Q) = \frac{5}{36 (2\pi )^2} \alpha _s^3(2 m_Q) \frac{|{\mathcal {R}}_{{\mathcal {Q}}}(0)|^2}{m_{{\mathcal {Q}}}^3} \nonumber \\&\quad \times \int _0^z \mathrm{d}\xi \int _{(\xi +z^2)/2z}^{(1+\xi )/2} \mathrm{d}\tau \frac{1}{(1-\tau )^2 (\tau -\xi )^2 (\tau ^2-\xi )^2} \nonumber \\&\quad \times \sum _{i=1}^{2} z^i \left[ f^{(g)}_i (\xi , \tau ) + g^{(g)}_i (\xi , \tau ) \frac{1+\xi -2 \tau }{2 (\tau -\xi ) \sqrt{\tau ^2-\xi }}\nonumber \right. \\&\quad \left. \times \ln \left( \frac{\tau - \xi + \sqrt{\tau ^2-\xi }}{\tau - \xi - \sqrt{\tau ^2-\xi }} \right) \right] , \end{aligned}$$

(23)

with the six $f^{(g)}_i$ and $g^{(g)}_i$ functions being given in Eqs. (4)–(9) of Ref. [253]. We stress that the mechanism considered here is the direct production of the quarkonium from the parent gluon. Another contribution to $J/\psi$-meson production in high energy processes, not considered in our analysis, is the production of a P-wave charmonium state $\chi _c$, followed by its radiative decay $\chi _c \rightarrow J/\psi + \gamma$ (see Refs. [35, 254]). The different initial energy scales at which quarks and gluons FF are taken is due to the production mechanism itself. The heavy-quark fragmentation involves at least three heavy quarks in the final state (Fig. 1, left panel), thus the running coupling in Eqs. (19) and (20) is calculated at $\mu _R = 3m_Q$. Conversely, the gluon fragmentation involves two heavy quarks only (Fig. 1, right panel), and this explains why the running coupling in Eq. (23) is taken at $\mu _R = 2m_Q$.

For the present study we follow twofold strategy. As a first step, starting from the initial-scale FF in Eq. (19), we generate the corresponding functions for all parton species. This is a required step in order to perform analyses by means of our high-energy VFNS treatment (see Sect. 2.1). For a given quarkonium, $J/\psi$ or $\Upsilon$, we set the corresponding constituent heavy (anti-)quark FF, $D^{J/\psi }_c \equiv D^{J/\psi }_{{\bar{c}}}$ or $D^{\Upsilon }_b \equiv D^{\Upsilon }_{{\bar{b}}}$, to be equal to the parameterization given in Eq. (19) at the initial scale $\mu _0 = 3 m_Q$. Then, we compute the DGLAP-evolved functions via the APFEL++ library [255,256,257], thus getting for each meson a LHAPDF set of FFs that embodies all parton flavors. From now, according to names of Authors of Ref. [36], we will refer to these sets as the $J/\psi$ and $\Upsilon$ ZCW19 NLO FF parameterizations. This first strategy of generating the remnant parton FFs only via DGLAP evolution could be justified by the fact that, by definition, the lowest Fock state of a quarkonium is a $(Q {\bar{Q}})$ pair, whose contribution heavily dominates among the other flavors. Moreover, this choice is in line with general assumptions made in the extraction of FFs of other heavy-flavored hadrons, as $D^{(*)}$ mesons [258,259,260,261], B mesons [262, 263] and $\Lambda _c$ [264] baryons.

As a second step, having in mind the outcome of the previously mentioned study on gluon fragmentation to vector quarkonium [38], we improve our FF treatment by starting the DGLAP evolution from the gluon input of Eq. (23) at the initial scale $\mu _F = 2 m_Q$. Then, when the energy grows up to reach $\mu _F = 3 m_Q$, the heavy-quark channel is opened and the corresponding FF (Eq. (19)) starts its evolution. In this way, for each meson we get another LHAPDF FF set, named ZCW19$^+$, which embodies both the gluon and the heavy-quark NRQCD inputs to fragmentation. For our phenomenology we will mainly employ the ZCW19$^+$ functions, and we will use the ZCW19 ones for cross-checks.

For comparison with previous studies on inclusive semi-hard emissions of heavy-flavored hadrons [216, 217], in Fig. 3 we plot $\mu _F$-dependence of our ZCW19 (upper panels) and ZCW19$^+$ (lower panels) FF sets for a value of the quarkonium longitudinal momentum fraction that roughly corresponds to the average value of z at which collinear FFs are typically probed in the consider kinematic configurations, namely $z = 5 \times 10^{-1} \simeq \langle z \rangle$. As expected, the FF for the constituent heavy (anti-)quark, $c({\bar{c}})$ for $J/\psi$ (left panels) and $b({\bar{b}})$ for $\Upsilon$ (right panels), strongly prevails over the other ones. Notably, as already observed in previous analyses on $\Lambda _c$ baryons [216] and bottomed hadrons [217], the gluon FF grows with $\mu _F$ up to reach a plateau. Although being much smaller than the constituent heavy-quark one, the gluon-FF contribution is strongly enhanced by the gluon PDF in the diagonal convolution embodied in the LO quarkonium impact factor (see Eq. (8)) and this feature is kept also at NLO, where (qg) and (gq) non-diagonal channels are active, but their weight cannot compensate the (gg) one. As pointed out in the aforementioned studies (see Section 3.4 of Ref. [216] and Appendix A of Ref. [217]), the $\mu _F$-behavior of the gluon FF plays a crucial role in the stabilization of the high-energy resummation under higher-order corrections. More in particular, a clear evidence was brought that a smooth-behaved, non-decreasing with $\mu _F$ gluon FF leads to a strong stabilizing effect on rapidity-differential cross sections, that also reflects in a partial stabilization final-state azimuthal distributions. Furthermore, this translates in a weaker sensitivity of semi-hard observables on renormalization- and factorization-scale variations. In view of these consideration, we expect that these stabilizing effects could fairly emerge in our quarkonium-plus-jet phenomenological analysis. We observe that the inclusion of the initial-scale gluon FF leads to a global lowering of the DGLAP-evolved heavy-quark one, which is smaller in the ZCW19$^+$ case. The evolved gluon FF is qualitatively similar in both case, the corresponding ZCW19$^+$ function being larger (smaller) than the ZCW19 one in the low (high) $\mu _F$ range.

We describe jet emissions at NLO perturbative accuracy in terms of a reconstruction algorithm calculated within the so-called “small-cone” approximation (SCA) [265, 266], i.e. for a small-jet cone aperture in the rapidity/azimuthal angle plane. More in particular, we adopt the version derived in Ref. [267], which is infrared-safe up to NLO perturbative and well-suited for numerical computations, with a cone-jet function selection [176] and for the jet-cone radius fixed at ${\mathcal {R}}_J = 0.5$, as usually done in recent experimental analyses at CMS [268]. The expression for the NLO jet vertex can be obtained, e.g., by combining Eq. (36) of Ref. [171] with Eqs. (4.19)-(4.20) of Ref. [176].

Fig. 3

Energy-scale dependence of ZCW19 (upper) and ZCW19$^+$(lower) NLO FFs depicting $J/\psi$ (left) and $\Upsilon$ (right) emission, for $z = 5 \times 10^{-1}$

Fig. 4

Upper panels: $\Delta Y$-distribution in the $J/\psi$ $+$ jet (left) and in the $\Upsilon$ $+$ jet (right) channel, for $\sqrt{s} = 13$ TeV. Quarkonium fragmentation is described in terms of ZCW19$^+$ functions. Central panels: the same observable in double quarkonium production channels. Lower panels: predictions for the quarkonium $+$ heavy-hadron detection. Text boxes inside panels show transverse-momentum and rapidity ranges. The combined effect of scale variation and phase-space numerical integration is embodied in uncertainty bands

3 Numerical analysis

The numerical elaboration of all the considered observables was done by making use of the JETHAD modular work package [224]. The sensitivity of our results on scale variation was assessed by letting $\mu _R$ and $\mu _F$ to be around their natural values or their BLM optimal ones, up to a factor ranging from 1/2 to two. The $C_{\mu }$ parameter entering plots represents the ratio $C_\mu = \mu _{R,F}/\mu _N$. Error bands in our figures embody the combined effect of scale variation and phase-space multi-dimensional integration, the latter being steadily kept below 1% by the JETHAD integrators. All calculations of our observables were done in the $\overline{\mathrm{MS}}$ scheme. BLM scales are calculated by solving the integral equation (12) in the MOM scheme.

3.1 $\Delta Y$-distribution

The first observable that we consider in our phenomenological analysis is the so-called $\Delta Y$-distribution. It essentially corresponds to the $\varphi$-summed cross section differential in the rapidity interval, whose expression is obtained by integrating the ${\mathcal {C}}_0$ azimuthal coefficient (see Eq. (5)) over rapidities and transverse momenta of the two emitted objects, and imposing the delta condition coming from fixing the value of $\Delta Y$. One has

$$\begin{aligned} C_0= & {} \int _{y_{{\mathcal {Q}}}^{\mathrm{min}}}^{y_{{\mathcal {Q}}}^{\mathrm{max}}} \mathrm{d}y_{{\mathcal {Q}}}\int _{y_J^{\mathrm{min}}}^{y_J^{\mathrm{max}}} \mathrm{d}y_J \int _{p_{{\mathcal {Q}}}^{\mathrm{min}}}^{p_{{\mathcal {Q}}}^{\mathrm{max}}} \mathrm{d}|\vec p_{{\mathcal {Q}}}| \int _{p_J^{\mathrm{min}}}^{p_J^{\mathrm{max}}} \mathrm{d}|\vec p_J| \, \,\nonumber \\&\times \delta (\Delta Y- (y_{{\mathcal {Q}}}- y_J)) \, \, {\mathcal {C}}_0\left( |\vec p_{{\mathcal {Q}}}|, |\vec p_J|, y_{{\mathcal {Q}}}, y_J \right) . \end{aligned}$$

(24)

Following the choice made in previous works on forward emissions of heavy-flavored hadrons [216, 217], we set the rapidity range of the quarkonium so that its detection is done only by the CMS barrel detector and not by endcaps, i.e. $|y_{{\mathcal {Q}}}| < 2.4$. Then, for the sake of consistency with the VFNS treatment, where energy scales need to be much larger than thresholds for DGLAP evolution of heavy quarks, we allow the quarkonium transverse momentum to be in the range 20 GeV $< |\vec p_{{\mathcal {Q}}}|<$ 60 GeV. The light jet is tagged in kinematic configurations typical of current studies at the CMS detector [268], namely 35 GeV $< p_J<$ 60 GeV and $|y_J| < 4.7$.

In upper panels of Fig. 4 we compare the NLA $\Delta Y$-behavior of $C_0$ with the corresponding prediction at LLA. Predictions were obtained by making use of the ZCW19$^+$ set. We note that values of $C_0$ are everywhere larger than 0.5 pb in the $J/\psi$-plus-jet channel (left) and $5\times 10^{-2}$ pb in the $\Upsilon$-plus-jet one (right). This leads to a very promising statistics, although being substantially lower than the one for heavy-baryon and heavy-light meson emissions [216, 217]. The downtrend with $\Delta Y$ of our distributions both at LLA and NLA is a common feature of all the hadronic semi-hard reactions investigated so far. It emerges as the interplay of two competing effects. On the one hand, the pure high-energy evolution leads to the well-known growth with energy of partonic cross sections. On the other hand, collinear parton distributions and fragmentation functions quench hadronic cross sections when $\Delta Y$ increases. We observe a partial stabilization of the high-energy series, with NLA bands almost overlapped to LLA ones at lower values of $\Delta Y$, and their mutual distance becoming wider in the large-$\Delta Y$ range. This feature is in line with previous studies on semi-hard heavy-flavor production, where the impressive stability of $C_0$ on next-to-leading order corrections observed in di-hadron channels ($\Lambda _c + \Lambda _c$ [216] and double b-flavored hadron [217]) is partially lost when light-jet emissions are instead considered. When a vector quarkonium is produced, the uncertainty band of $C_0$ at NLA is almost always smaller than the LLA one. To further examine this aspect, we compare our quarkonium-plus-jet predictions with results obtained when the light jet is replaced by another heavy-flavored hadron produced in the same final-state kinematic window. In central panels of Fig. 4 we show the $\Delta Y$-behavior of $C_0$ for a double quarkonium production (double $J/\psi$ or $\Upsilon$). Here, the behavior of $C_0$ is similar to the one observed in quarkonium-plus-jet channel. In particular, LLA and NLA bands decouple when $\Delta Y$ grows. In lower panels of Fig. 4 we consider the production of a quarkonium plus a hadron sharing the same heavy flavor. More in particular, the $J/\psi$ is accompanied by a $\Lambda _c^\pm$ baryon, described in terms of the KKSS19 NLO FF parameterization [264], while the $\Upsilon$ is emitted together a b-flavored hadron (${\mathcal {H}}_b$), portrayed by the KKSS07 NLO FF set [263, 269]. We observe that, at large $\Delta Y$, the NLA band for the $J/\psi$-plus-$\Lambda _c$ production is slightly wider than the double $J/\psi$ one (left central panel) and has almost the same size of the $J/\psi$-plus-jet one (left upper panel). Conversely, the NLA band for the $\Upsilon$-plus-${\mathcal {H}}_b$ tag is almost always smaller than the double $\Upsilon$ one (right central panel) and the $\Upsilon$-plus-jet one (right upper panel). All these observations support the statement that the high-energy stabilizing power of $J/\psi$ collinear FFs is stronger than the $\Lambda _c$ FF one, while $\Upsilon$ FFs act as weaker stabilizers with respect to the ${\mathcal {H}}_b$ ones.

Fig. 5

$\Delta Y$-distribution in the $J/\psi$ $+$ jet (left) and in the $\Upsilon$ $+$ jet (right) channel, for $\sqrt{s} = 13$ TeV. Quarkonium fragmentation is described in terms of ZCW19 (upper) or ZCW19$^+$ (lower) functions. A study on progressive energy-scale variation in the range $1< C_\mu < 30$ is made. Ancillary panels below primary plots exhibit reduced distributions, namely divided by $C_0$ taken at $C_\mu = 1$

Fig. 6

BLM scales for the quarkonium-plus-jet (left) and the double quarkonium (right) production as functions of the rapidity separation, $\Delta Y$, for $C_0$, and for $\sqrt{s} = 13$ TeV. Predictions for $J/\psi$ and $\Upsilon$ emissions, described in terms of ZCW19$^+$ FFs, are compared with configurations where other heavy-flavored hadron species are detected: $\Lambda _c$ baryons and ${\mathcal {H}}_b$ hadrons. Text boxes inside panels show transverse-momentum and rapidity ranges

In Fig. 5 we study the $\Delta Y$-behavior of the NLA-resummed $\Delta Y$-distribution under a progressive variation of $\mu _R$ and $\mu _F$ scales in the wider range given by $1< C_\mu < 30$. Upper (lower) plots contain predictions obtained via ZCW19$^{(+)}$ FFs. Ancillary panels below primary plots show the reduced distributions, i.e. divided by $C_0$ taken at $C_\mu = 1$. As a first remark, we note that the magnitude of $C_0$ is almost independent on the choice of the quarkonium FF set. Notably, the sensitivity of ZCW19 predictions on the $C_\mu$ parameter is rather weak for the $J/\psi$-plus-jet channel, while it becomes more evident at larger values of $\Delta Y$ for the $\Upsilon$-plus-jet reaction. Conversely, when ZCW19$^{+}$ functions are used, $J/\psi$-plus-jet emissions exhibit almost the same patter, while $\Upsilon$-plus-jet detections become less sensitive to scale variation when $\Delta Y$ grows. This is an indication that employing a complete NRQCD input, built in terms of both initial-scale heavy-quark and gluon functions, leads to a stronger stabilizing pattern than the one achievable by relying on the quark input only. Combining this information with the arguments presented above, we can easily explain the size of error bands for NLA results in upper panels of Fig. 4. Indeed, the major contribution to scale uncertainty comes, at least for the $J/\psi$-plus-jet emission, from variation of $\mu _R$ and $\mu _F$ in the $1/2< C_\mu < 1$ lower sub-range. Here, the value of energy scales comes closer, even if still larger, to DGLAP-evolution thresholds connected to heavy-quark masses and this generates some instabilities which in turn lead to an increased scale-variation sensitivity of our predictions.

As a supplementary analysis, we present results for the scale parameter $C_\mu \equiv C_\mu ^\mathrm{BLM}$ obtained by applying the BLM optimization method on $C_0$. In Fig. 6 we compare the $\Delta Y$-dependence of $C_\mu ^\mathrm{BLM}$ for our quarkonium production channels, described by means of the ZCW19$^+$ set, with the corresponding one for the emission of other c- and b-flavored non-quarkonium bound states, again $\Lambda _c$ baryons and b-hadrons. In Refs. [216, 217] it was shown that the BLM-scale values obtained when those heavy-flavored hadrons are detected in the final state are sensibly lower than the ones typical of lighter-hadron tags (see, e.g., Refs. [199, 207, 224]). Since the use of the BLM prescription operationally leads to a growth of $C_\mu ^\mathrm{BLM}$ to dampen the weight of next-to-leading order corrections, processes where smaller BLM scales are found natively embody an intrinsic, partial stabilization of the high-energy series, before adopting optimization procedure. The inspection of results presented in Fig. 6 fairly confirms the statement that those stabilization effects are present also when $J/\psi$ and $\Upsilon$ mesons are considered. More in particular, left panel of Fig. 6 show how $C_0$-related BLM scales are clearly smaller in the double $J/\psi$ channel with respect to the double $\Lambda _c$ one, while they are slightly larger in the double $\Upsilon$ channel with respect to the double b-hadron one. Then, in line with studies done in Ref. [217], b-flavored quarkonium and non-quarkonium detections provide us with smaller scales than the ones obtained for c-flavored tags. The same hierarchy is found in the heavy-hadron plus light jet channel (Fig. 6, right panel), matter of investigation in this article.

The general outcome of results presented in this Sections is that studies on forward $J/\psi$ and $\Upsilon$ emissions via ZCW19$^{(+)}$ NLO collinear FFs lead to a favorable statistics for the $\Delta Y$-differential cross section. Effects of stabilization of the high-energy resummation under higher-order corrections and scale variation are presents and are even stronger, at least for the c-flavor channel, when quarkonia are detected instead of other heavy-flavored hadrons. At the same time, the simultaneous tag of a quarkonium together a light-flavored jet could partially shadow these effects, and further, dedicated studies on the improvement of the theoretical description of NLO jet emissions are needed.

Another relevant point that deserves a discussion is the weight of MPI contributions to cross sections at large $\Delta Y$. Let us focus on the double parton scattering (DPS), namely where at most two hard-scattering subprocesses are considered. In Ref. [178] it was shown that DPS contributions to Mueller–Navelet jets can be potentially relevant for $\Delta Y$-distributions at large s and low $p_T$, while, in the case of angular correlations between the two jets, they lead to results compatible with a NLA description based on single parton scattering only. Analyses on double $J/\psi$ [65, 70], $J/\psi$ $+$ $\Upsilon$ [65], $J/\psi$ $+$ Z [64], and $J/\psi$ $+$ W [270] production processes have provided a clear indication that DPS effects are relevant and they need to be accounted for in studying the low-$p_T$ spectrum. Possible MPI signatures are expected also in triple $J/\psi$ events [271, 272]. Therefore, the search for DPS imprints in our vector-quarkonium plus jet observables is an important task that needs to be addressed in future studies.

3.2 Azimuthal distribution

The study of the azimuthal-angle distribution in semi-hard two-particle final-state reactions is widely recognized as one of the most fertile grounds where to hunt for distinctive signals of the BFKL dynamics. Here, the weight of undetected gluons strongly ordered in rapidity, accounted for the resummation of high-energy logarithms, becomes more and more relevant when the rapidity interval $\Delta Y$ between the two objects grows. Thus, the decorrelation of the detected particles on the azimuthal plane increases at large $\Delta Y$ and the number of back-to-back events diminishes.

First observables where the azimuthal-angle decorrelation was observed are the so called azimuthal-correlation moments, namely the ratios of azimuthal coefficients $R_{n0} \equiv C_n/C_0$. Here, $C_0$ is the $\varphi$-summed $\Delta Y$-distribution introduced in Sect. 3.1, while $C_{n \ge 1}$ are the azimuthal coefficients integrated over the final-state phase space in the same way as $C_0$ (see Eq. (24)). The $R_{n0}$ moment has an immediate physical interpretation, being the average value of the cosine $\langle \cos \varphi \rangle$, and the $R_{n0}$ ones correspond to the higher moments $\langle \cos (n \varphi ) \rangle$. The study of ratios between azimuthal correlations, $R_{nm} \equiv C_n/C_m = \langle \cos (n \varphi ) \rangle / \langle \cos (m \varphi ) \rangle$, was proposed in Refs. [273, 274]. Comparison of NLA predictions for azimuthal ratios in the Mueller–Navelet light-di-jet channel have shown a nice agreement with LHC data at $\sqrt{s} = 7$ TeV and for symmetric $p_T$-ranges of the two emitted jets [172, 173, 175].

Fig. 7

LLA (upper panels) and NLA (lower panels) predictions for the $\varphi$-distribution in the $J/\psi$ $+$ jet (left) and in the $\Upsilon$ $+$ jet (right) channel, at $\sqrt{s} = 13$ TeV, and for three distinct values of $\Delta Y$. Quarkonium fragmentation is described in terms of ZCW19$^+$ functions. Text boxes inside panels exhibit final-state kinematic cuts. Uncertainty bands embody the combined effect of scale variation and phase-space multi-dimensional integration

As already mentioned, a major issue emerging in the NLA BFKL description of the Mueller–Navelet di-jet azimuthal ratios is the rise of instabilities of the high-energy series, which turned out to so be strong to prevent any possibility to perform realistic analyses at natural scales [172, 175, 177]. This phenomenon was later observed also in the inclusive semi-hard light di-hadron [187, 224] and hadron-jet production [224]. Recent studies on inclusive forward heavy-flavored hadrons have shown how the stabilization effect coming from the use of heavy-flavor VFNS FFs is quite strong where di-hadron final-state systems are considered, but it is less effective on hadron-plus-jet concurrent emissions [216, 217].

In this Section we focus on an alternative azimuthal-angle dependent observable, namely the $\varphi$-distribution, or simply azimuthal distribution

$$\begin{aligned} \frac{1}{\sigma } \frac{\mathrm{d}\sigma }{\mathrm{d}\varphi }= & {} \frac{1}{2 \pi } \left\{ 1 + 2 \sum _{n = 1}^{\infty } \cos (n \varphi ) \langle \cos (n \varphi ) \rangle \right\} \nonumber \\= & {} \frac{1}{2 \pi } \left\{ 1 + 2 \sum _{n =1 }^{\infty } \cos (n \varphi ) R_{n0} \right\} . \end{aligned}$$

(25)

Proposed for the first time in the context for Mueller–Navelet studies [172, 275], the $\varphi$-distribution brings with it two main advantages. On the theoretical side, it collects the high-energy signals coming from all the $C_n$ coefficients, thus representing one of the most solid observables where to hunt for resummation effects. On the experimental side, since due to detector limitations the measured distributions cannot cover the whole $2 \pi$ azimuthal-angle range, a $\varphi$-dependent observable turns out to be much easier to be compared with data, rather than a $R_{nm}$ ratio. At the same time, due to the large number of $C_n$ coefficients required, numerical computations of Eq. (25) can be time consuming. We checked the numerical stability of our calculation by progressively raising the effective upper limit of the n-sum in Eq. (25). An excellent numerical convergence was found at $n_{\mathrm{max}} = 20$.

In Fig. 7 we present predictions for the azimuthal distribution as a function of $\varphi$ and for three distinct values of the rapidity interval, $\Delta Y= 1, 3, 5$. Results for the $J/\psi$ $+$ jet hadroproduction are given in left panels, while the ones for the $\Upsilon$ $+$ jet hadroproduction are shown in right panels. Results were obtained by making use of the ZCW19$^+$ set. We adopted the same final-state kinematic cuts introduced in Sect. 3.1. A general feature of all the presented distributions is the presence of a clear peak at $\varphi = 0$. At this value the quarkonium and the jet are emitted back-to-back. In all cases the height of the peak visibly diminishes when $\Delta Y$ grows, while the distribution width slightly widens. This is distinctive signal of the onset of the high-energy dynamics. Indeed, at large values of $\Delta Y$ the weight of gluons strongly ordered in rapidity predicted by BFKL increases, thus bringing to a reduction of the azimuthal correlation between the quarkonium and the jet, so that the number of back-to-back events lowers.

Focusing on patterns for the $\varphi$-distribution in the $J/\psi$ $+$ jet channel, we observe that at $\Delta Y= 1$ the LLA peak is much more pronounced than the corresponding NLA one, at $\Delta Y= 3$ the LLA and NLA peaks are similar in height, and at $\Delta Y= 5$ the LLA peak is beyond the NLA one. This behavior as a straightforward explanation. The small-$\Delta Y$ range stays at the limit of applicability of the BFKL resummation, since the low values of the partonic center-of-mass energies reduces the phase space for secondary-gluon emissions. This leads to a stronger discrepancy between LLA and NLA results. Conversely, in the moderate-$\Delta Y$ regime the high-energy series shows a fair stability when NLA corrections are switched on. Finally, the large-$\Delta Y$ territory is very sensitive to the BFKL dynamics, this reflecting in a stronger weight of the re-correlation effects, predicted by the NLA resummation, over pure LLA results.

Considering $\Upsilon$ $+$ jet production, we observe that LLA predictions are very close to the $J/\psi$ $+$ jet corresponding ones, whereas NLA results exhibit a different pattern. In particular, peaks are globally smaller, width are larger and uncertainty bands around peaks are very close to each other as $\Delta Y$ increases. This is a further indication that $\Upsilon$ emissions have a less stabilizing power on the high-energy resummation than $J/\psi$ ones. At the same time, we remark that azimuthal distributions for both the considered channels feature a cleaner separation among results for different values of $\Delta Y$ with respect to what happens for other reactions. As an example, error bands for the $\varphi$-distribution for the inclusive heavy-light di-jet production are definitely more overlapped, both at LLA and NLA (see Fig. 9 of Ref. [210]).

The general outcome of results presented in this section is the possibility of studying the azimuthal distribution for quarkonium-plus-jet processes around natural values of $\mu _R$ and $\mu _F$ scales. This observable can be easily measured at the LHC, thus offering the possibility of doing stringent tests of the high-energy resummation.

4 Summary and outlook

We proposed the direct inclusive hadroproduction of a forward $^3S_1^{(1)}$ quarkonium, $J/\psi$ or $\Upsilon$, in association with a backward light jet in hybrid high-energy and collinear factorization. We described the vector-meson hadronization in terms of a novel model of heavy-quark to quarkonium VFNS collinear FFs [36], suited for the study of quarkonium production at large transverse momentum, where the single-parton fragmentation mechanism is expected to prevail on short-distance $(Q {\bar{Q}})$-pair production. We built two novel collinear FF sets, respectively named ZCW19 and ZCW19$^+$ functions. The first one was obtained trough the DGLAP evolution of a NRQCD input for the fragmentation of a heavy quark (charm or bottom) into a $^3S_1^{(1)}$ state at NLO [36]. The second one combines the heavy-quark input [36] with the gluon one, originally calculated in Ref. [28]. We believe that these functions can serve as useful tools for further phenomenological applications outside the high-energy domain, as well as for formal studies on the connection between the collinear factorization and the NRQCD effective formalism.

Motivated by recent studies on forward emissions of heavy-flavored hadrons in the same formalism [216, 217], we hunted for stabilizing effects of the high-energy resummation under higher-order corrections and energy-scale variation. These effects are present and are even stronger when a forward $J/\psi$ is tagged rather than another charmed hadron, as the $\Lambda _c$ baryon. Such clues are less evident when $\Upsilon$ emissions are compared with bottomed-hadron ones. This is due to larger values of the b-quark mass that push the threshold for DGLAP evolution of the FFs closer to the energies provided by the considered kinematic regime, thus slightly worsening the validity of a pure VFNS description. Future investigations are needed to gauge the dependence of our predictions on intrinsic features of the jet-emission description, as the jet selection function.

This work represents a further step in our ongoing program on heavy-flavored emissions, started from the analytic calculation of heavy-quark pair impact factors [213, 214, 276], and a first contact with the phenomenology of quarkonium production at high energies. We plan to extend our analyses by considering quarkonium detections in wider kinematic ranges, as the ones accessible at new-generation colliding facilities, as the Electron-Ion Collider (EIC) [277,278,279], NICA-SPD [280, 281], the High-Luminosity LHC [282,283,284], the Forward Physics Facility (FPF) [285,286,287], and the International Linear Collider (ILC) [288]. Here, our hybrid high-energy and collinear factorization could serve as an additional theoretical tool to perform quarkonium studies that cover a broad $p_T$-spectrum, with the aim of shedding light on the transition region from the short-distance $(Q {\bar{Q}})$-pair production to the single-parton fragmentation mechanism.