High-energy QCD dynamics from bottom flavor fragmentation at the Hi-Lumi LHC

We study the inclusive production of hadrons with bottom flavor at the LHC and its luminosity upgrade. We describe the collinear fragmentation of singly $b$-flavored hadrons, $B$ mesons and $\Lambda_b$ baryons, via the KKSS07 determination of fragmentation functions, while for charmed $B$ mesons, $B_c(^1S_0)$ and $B_c(^3S_1)$ particles, we employ the novel ZCFW22 set, built on the basis of state-of-the-art nonrelativistic QCD inputs. We use the JETHAD multimodular working environment to analyze rapidity and transverse-momentum distributions for observables sensitive to the associated emission of two hadrons or a hadron-plus-jet system. Our reference formalism is the NLL/NLO$^+$ hybrid collinear and high-energy factorization, where the standard collinear description is improved by the inclusions of energy logarithms resummed up to the next-to-leading approximation and beyond. We provide a corroborating evidence that $b$-flavor emissions act as fair stabilizers of the high-energy resummation, thus serving as valuable tools for precision studies of high-energy QCD. As a bonus, we highlight that the predicted production-rate hierarchy between noncharmed $b$-hadrons and charmed $B_c(^1S_0)$ mesons is in line with recent LHCb estimates. This serves as simultaneous benchmark both for the hybrid factorization and for the single-parton fragmentation mechanism based on initial-scale inputs calculated via the nonrelativistic QCD effective theory.


Opening remarks
Exploring the heavy-flavor domain stands as a frontier research field whereby addressing unresolved quests of particle physics.Here, possible experimental evidences of the associated production of heavy quarks and Beyond-the-Standard-Model (BSM) particles are milestones in the search for long-awaited signals of New Physics.Still they can play a crucial role in precision studies of strong interactions, the values of charm-and bottom-quark masses making perturbative Quantum ChromoDynamics (QCD) applicable.
Special emphasis deserves the hadroproduction of the bottom quark, which represents the heaviest quark species that can fragment and generate hadrons.A rigorous description of ( b) pair production by the hands of collinear factorization and within the Next-to-Leading Order (NLO) perturbative accuracy was formalized at the end of the eighties [1][2][3][4], while detailed analyses on differential distributions at Next-to-NLO (NNLO) have been undertaken only recently [5,6].A key actor here is the bottom mass,   .It genuinely marks the transition point between two distinct factorization schemes.
For -flavored particles emitted with low transverse momentum, so that   ≳ | ⃗ |, a suitable description is found in the so-called Fixed Flavor Number Scheme (FFNS, see Ref. [7] for insight).It exclusively incorporates collinear Parton Distribution Functions (PDFs) and/or Fragmentation Functions (FFs) of light partons (light-flavored quarks and the gluon).To be precise, here we refer to a scheme with a number of active flavors   =   , with   being the number of light flavors. 1 Additionally, in the FFNS, heavy quarks manifest solely in the final state and their masses must not be set to zero.This is a required condition not to face perturbative-convergence issues due to the rise of power corrections of the form | ⃗ |∕  .Conversely, when the transverse momentum is large,   ≪ | ⃗ |, contributions scaling with ln(| ⃗ |∕  ) progressively grow and demand for an all-order logarithmic resummation of collinear logarithms of the factorization scale over the mass of the heavy quark(s) [8,9].In that case, the Zero-Mass Variable-Flavor Number Scheme (ZM-VFNS, or simply VFNS) emerges as the most suitable framework whereby combining fixed-order radiative corrections with the resummation [9][10][11][12][13][14]. Within the VFNS, all quark flavors are treated as massless particles and participate to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution.The flavor number increases by one every time a heavy-quark threshold is crossed.
An elegant matching between the FFNS scheme and the VFNS one is realized by means of the so-called General-Mass Variable-Flavor Number Scheme (GM-VFNS).It embodies calculations involving both massive quarks at lower scales and massless quarks at higher scales, while the heavy-quark masses serve as control parameters for the transition between the two descriptions.Several implementations of the GM-VFNS have been proposed so far [15][16][17][18][19].For an exhaustive exploration of these studies, we direct the reader to Refs.[10,12,20].
While the use of a GM-VFNS undoubtedly improves the description of both the kinematic limits, it still suffers from some ambiguities.On one hand, it systematically misses power corrections rising from initial-state quarks, which might be quenched, however, by the smallness of heavy-flavor PDFs.On the other hand, it does not behave as desired in the intermediate region,   ∼ | ⃗ |, since the difference between the pure ZM-VFNS contribution and its massless limit does not automatically vanish and it must be suppressed by hand (see Ref. [16] for a detailed discussion).
A simplest solution comes instead from the use of a ZM-VFNS with a suitable choice of the Heavy-Flavor Matching Point (HFMP), namely the transition scale starting from which the massive flavor is treated as if it were massless.As extensively shown in Ref. [21], numerically accurate predictions can be achieved with HFMP(s) ranging from five to 10 times the mass of the heavy quark(s).The main benefits of this choice are: the possibility of keeping the size of missing power corrections steadily below 4 ÷ 1%, a fully preserved resummation accuracy, and an improved continuity in predictions for observables.
Notably, bottom-flavored emissions are useful not only to study direct-production channels, but also to access the top sector.For this reason, mechanisms governing -quark fragmentation are believed to play an important role within the top-physics domain.The production of -particles at hadron colliders in NNLO QCD and its application to ( t) events with leptonic decays were extensively investigated in Refs.[22,23].The same approach finds application in precision studies of electroweak interactions, particularly in the analysis of photon radiation originating from massive charged fermions, like via ( b) decays from a Higgs boson [24].Refs.[25][26][27] investigate the role of the bottom quark in the associated production of lepton pairs.The semi-inclusive emission of a Higgs boson decaying into -quarks accompanied by a weak vector within NNLO QCD was considered in Ref. [28].Ref. [29] examines electroweak emissions from heavy fermions in  ± boson production processes with Monte Carlo techniques.
The total cross section for the production of a Higgs boson or a  boson in bottom fusion and via the FONLL [30,31] matching procedure 2 was computed in Refs.[32,33] and [34], respectively.The associated production of an electroweak boson with -quarks was investigated in Refs.[35,36].Effects of collinear logarithms and finite bottom-mass powers in ( b) tags were addressed in Refs.[37,38].Refs.[39,40] contain a detailed study on the inclusion of the resummation of soft logarithms at NNLO in the perturbative component of the bottom fragmentation function.Matching issues arising from the the noncommutativity of the light-mass limit and soft-emission one were recently highlighted [41].Analyses on heavy-quark fragmentation in lepton collisions at NNLO plus logarithmic corrections were presented in Refs.[42,43].QCD corrections to the hadroproduction of ( t b) systems were considered in Refs.[44][45][46][47][48].
Studies on bottom-flavored jets and hadrons at high energies can be found in Refs.[49][50][51][52].Ref. [53] investigates kinematic correlations of lepton pairs from semi-leptonic decays of  mesons.Multiple parton interactions in double  + productions at the LHC were considered in Ref [54].The VFNS collinear fragmentation to -mesons was first studied at NLO in Ref. [55].The obtained FFs came out from a fitting procedure on LEP1 data.A subsequent fit to data from LEP1 and SLAC-SLC was made [56].The resulted FF parametrization was then employed in the computation of NLO cross section for the inclusive hadroproduction of  mesons within the GM-VFNS.There, -quark finite-mass effects were also assessed [57].
Predictions for the semi-inclusive emission of singly -flavored particles, encompassing  mesons and Λ  baryons which we collectively denote as   hadrons, were recently compared with LHCb and CMS data [58].That analysis relied upon the assumption that a single set of functions could aptly describe the VFNS fragmentation all species of   particles.Consequently, FFs for any singly bottomed hadron species might be derived from the general ones by just multiplying it by an energy-independent branching fragmentation fraction.However, analyses conducted by the Heavy Flavor Averaging Group (HFAG) highlighted a departure from the universality assumption regarding the branching fraction, particularly in the case of Λ  tags, as evidenced by data from LEP and Tevatron [59].On the contrary, the universality seems to be valid of -meson emissions.A recent examination of semi-inclusive Λ  rates at LHCb and CMS emphasized the need to go beyond the branching-fraction description when the observed transverse momentum enlarges [60].Therefore, for the time being the use of the FF determination of Ref. [56,58] is safe provided that only  mesons or both  mesons and Λ  baryons are inclusively considered, but it is not valid for Λ  baryons only.Forthcoming data with diminished experimental uncertainties will contribute to a better understanding of this aspect.
Particular relevance in accessing the heavy-flavor sector deserves the case of mesons with both the charm and the bottom quark in their leading Fock level: | b⟩ for positive-charged particles, | c⟩ for negative ones.Because of the presence of two heavy quarks, these charmed  mesons can be thought of as generalized quarkonium states.However, contrarily to charmonia and bottomonia, they cannot annihilate into gluons.Thus, they are quite stable and exhibit narrow decay widths (for novel determinations of the   lifetime, see Refs.[61][62][63]).Being top quarks extremely short-lived and not able to hadronize, charmed  mesons actually represent the final frontier of meson spectroscopy (see, e.g.Ref. [64]).The CDF Collaboration at Tevatron observed the   ≡   ( 1  0 ) particle for the first time in 1998 [65].The  *  ≡   ( 3  1 ) resonance was detected by ATLAS only in 2014 [66].The simultaneous presence of a charm and a bottom quark in their leading Fock level makes charmed  mesons excellent testing grounds to shed light on the core nature of strong and weak interactions.Analyses on direct emissions as well as on indirect productions from electroweak decays stand as gold-plated channels whereby testing  ( * )  formation mechanisms as well as unraveling peculiar features of their decay products.As an example, they are useful to investigate rare Higgs decays [67][68][69].
Data analyzed by LHCb have shown that   meson production rates are three orders of magnitudes lower than   states, at most [70,71].This explains why charmed -meson channels are generally excluded from fits of -hadron fragmentation parameters, thus making our knowledge of their formation mechanism quite scarce.For this reason, any attempt at catching the core fragmentation dynamics of  ( * )   particles still needs to rely upon model-dependent studies and/or effective approaches.
An intriguing opportunity to access the production mechanism of charmed  mesons comes from an effective theory, known as NonRelativistic QCD (NRQCD) [72][73][74].Originally developed to explain issues in the description of charmonium and bottomonium production, the NRQCD framework builds on the assumption that, to the observed bound state, all the possible Fock states contribute through a linear combination.All these contributions are ordered in a double expansion in powers of the strong coupling,   , and the nonrelativistic relative velocity between two constituent heavy quarks, .
The question whether a NRQCD-inspired description is valid for the perturbative component of the collinear fragmentation of a single parton to the detected  ( * ) meson, or it should be rather adopted for the short-distance production of a charm-plus-bottom system directly in the hard subprocess, strictly depends on the transverse momenta at play.Indeed, while the short-distance mechanism suitably describes the production of the two constituent heavy quarks at low transverse momentum, namely when they feature a relative transverse separation of order 1∕| ⃗ |, when | ⃗ | grows the fragmentation of a single parton, followed by its inclusive hadronization into the physical state, starts to prevail.Phenomenological efforts to unravel the transition region between the two regimes were originally made in the context of charmonia [75][76][77][78][79].The main outcome was that the (gluon) fragmentation contribution starts to dominate when | ⃗ | ≳ 10 ÷ 15 GeV.An analogous threshold was then set also for charmed -mesons [80], whereas more recent analyses [81] pointed out how this lower bound might be larger.
Although NRQCD can model the initial-scale input of the  ( * )  fragmentation, the other key ingredient for a proper collinear description is the DGLAP evolution.Thus, by starting from NLO NRQCD calculations of heavy quarks [82] and gluon [83] fragmentation channels, first sets of VFNS-compatible, DGLAP-evolving FFs for   ( 1  0 ) and   ( 3  1 ) mesons were derived in Ref. [84].They were named ZCFW22 NLO functions.A similar strategy was applied to obtain the ZCW19 + NLO FFs for vector quarkonia [85,86].
In this article we will study the inclusive emission, at the LHC ad its luminosity upgrade, of a   hadron or a  ( * )  meson accompanied by another -particle or by a light-flavored jet, as a testing ground for the manifestation of imprints of high-energy QCD dynamics.We will make use of the hybrid collinear and high-energy factorization framework [85,[87][88][89][90][91], as implemented in the JETHAD code [86,88,91] to analyze rapidity-interval and double transverse-momentum differential rates for observables sensitive to the QCD resummation of high-energy logarithms within the next-to-leading approximation and beyond.
We will hunt for stabilizing effects of our distributions under radiative corrections and energy-scale variations, discovering that they are present and allow for a consistent description of these observables at LHC energies and at the natural scales provided by process kinematics.We will come out with a bonus result, namely that the aforementioned production-rate hierarchy between charmed   ( 1  0 ) mesons and noncharmed -hadrons is in fair agreement with recent LHCb findings [70,71].This will serve as simultaneous benchmark both for our resummed calculations and for the single-parton fragmentation mechanism from NRQCD.
The structure of this article reads as follows.Section 2 gives introductory remarks and technical details on NLL∕NLO (+) hybrid collinear and high-energy factorization applied to -flavor production.Section 3 contains phenomenological analyses on our high-energy distributions.Section 4 summarizes results and provides outlooks.

Bottom flavor production within NLL∕NLO + hybrid factorization
The first part of this Section (2.1) consists in a brief summary of recent phenomenological progresses concerning the QCD semi-hard regime.Then (2.2) the NLL∕NLO + hybrid collinear and high-energy is explained and adapted to the study of inclusive double -hadron and -hadron plus light-jet final states at the Hi-Lumi LHC.The third part (2.3) provides us with details on collinear ingredient.Finally in (2.4) the natural stability [112] of -flavor sensitive distributions is discussed.

High-energy resummation at work
The theoretical description of hadron collisions at high energies relies upon the well-defined collinear factorization, where long-distance and short-distance dynamics are decoupled.This remarkable features permits to factorize the nonperturbartive hadronic content from the perturbatively calculable information on parton scatterings.Nevertheless, there exist particular kinematic limits where the purely collinear description is spoiled by the emergence of large logarithmic corrections.Their impact on cross sections can be enough sizable to offset the narrowness of the QCD running coupling,   , thus harming the perturbative series convergence.In such cases the collinear factorization must be improved via the inclusion to all orders of one or more logarithmic resummations.
The kinematic sector matter of our analysis is the semi-hard regime (see Refs. [113,114] for novel applications), where the energy scale hierarchy √  ≫ {  } ≫ Λ QCD is strongly valid.Here, √  stands for the center-of-mass energy, {  } represents one or a set of reaction-related hard scales, and Λ QCD is the QCD hadronization scale.Whereas the second inequality simply states that perturbative QCD can be employed in our calculations, the first one heralds the fact that we have accessed the QCD Regge limit.Here, large ln(∕ 2  ) type logarithms enter the perturbative series with a power increasing with the order of   .Thus, any pure fixed-order perturbative techniques are not anymore applicable.As anticipated, they need to be enhanced by considering the all-order resummation of terms proportional to energy logarithms.
The Balitsky-Fadin-Kuraev-Lipatov (BFKL) approach [115][116][117] is the well-established tool to resum terms proportional to (  ln )  , namely the ones constituting the Leading Logarithmic (LL) approximation, and also factors proportional to   (  ln )  , namely the ones specific of the Next-to-Leading Logarithmic (NLL) approximation.Following the BFKL factorization method, any high-energy amplitude can be rewritten as a transverse-momentum convolution of a universal Green's function with two singly off-shell emission functions (also known as impact factors) describing the subprocess that leads, in the corresponding fragmentation region, to the emission of a forward or backward final-state particle from the stemming hadron.The Green's function evolves according to the BFKL integral equation, whose kernel has been calculated within the NLO fixed-order accuracy.Recent studies have been done to calculate its NNLO correction [118][119][120][121][122] (see also Ref. [123]).
A crucial problem rising from the BFKL description of Mueller-Navelet rapidity distributions and azimuthalangle correlations is the weight of the NLL part, which is of the same order, but of opposite sign than the pure LL case.This is at the origin of strong instabilities observed in the resummed series that become manifest when factorization and renormalization scales are varied around their natural values.Therefore, Muller-Navelet observables assume unphysical values when the jet rapidity separation becomes sufficiently large.In the same way, azimuthal Figure 1: Hybrid collinear and high-energy factorization for double-hadron detection (left) and hadron-plus-jet inclusive hadroproduction (right).Collinear parton densities are represented by red ovals.The green (blue) blob depicts the off-shell hard factor embodied in the hadron (jet) emission function, while indigo arrows depicts the emission of -flavored hadrons.The Green's function is described by the orange blob.The diagram was created with JaxoDraw 2.0 [211].
correlations exhibit an unphysical behavior both at small and at large rapidity distances.Different solutions have been tried to cure this problem.In particular, the Brodsky-Lepage-Mackenzie (BLM) prescription [201,202] as specifically built for semi-hard reactions [203] slightly dampens those instabilities on azimuthal correlations and thus averagely ameliorates the agreement with data.However, BLM is almost ineffectual on light di-hadron or hadron plus jet semihard distributions.The main reason is that the optimal values for renormalization scales prescribed by BLM are much higher than the process natural scales [88,113,157].Thus, total cross sections suffer from a large loss of statistics.
Fair signals of a stabilization of the high-energy resummation under NLL corrections and energy-scale variations recently emerged from semi-hard reactions whose final states show the signature of Higgs bosons [89,[204][205][206][207][208][209].Then, a fair stabilizing trend came out from the study of Λ  hyperons at the LHC [164], and also from analogous observables sensitive to  mesons and Λ  baryons [165].In those works it was proven that the stabilizing effect is encoded in the peculiar pattern of VFNS collinear FFs portraying the production of these singly heavy-flavored particles at high transverse momentum.Similar results were observed when vector quarkonia [85,86], charmed  mesons [84], or heavy-light tetraquarks [210] are emitted in semi-hard configurations.All these results gave us corroborating evidence that the found natural stabilization of BFKL emerges as an intrinsic property genuinely correlated with final states sensitive to heavy flavor [112].

NLL-resummed differential cross section
We will study the two kinds of processes depicted in Fig. 1 p where p( , ) is an initial-state proton possessing four-momentum  , , from which a parton with longitudinalmomentum fraction  , is struck.Then,  (1,2) ( 1,2 ,  1,2 ) stands for a singly bottom-flavored hadron   ≡ { + Λ  }, or a charmed-bottomed meson   ( 1  0 ) ≡   , or its resonance   ( 3  1 ) ≡  *  , emitted with four-momentum  1,2 , rapidity  1,2 and azimuthal angle  1,2 , and a light-favored jet is possibly tagged with four-momentum  2 , rapidity  2 and azimuthal angle  2 .Furthermore,  represents an inclusively radiated gluon system.The large final state transverse momenta, | ⃗  1,2 |, together with the high rapidity interval, Δ ≡  1 −  2 , are required ingredients to make the final state be both semi-hard and diffractive.In addition, transverse momenta of the observed particles need to be large enough to make the VFNS fragmentation approach be valid.
Stemming protons' four-momenta can be written as Sudakov vectors with the conditions  2  =  2  = 0 and 2(  ⋅   ) = , so that final-state four-momenta are cast via the following decomposition The outgoing-object longitudinal-momentum fractions are labeled as  1,2 .They are connected with the respective rapidities and to Δ via According to the pure collinear factorization, one would write the LO cross section for our processes (Eq.( 1)) as a one-dimensional convolution among the partonic-subprocess factor, proton PDFs, and -hadron FFs.Considering the double hadron case as depicted in the left diagram of Fig. 1, we have with the  and  indices running over parton species,  , ( 1,2 ,   ) being the proton PDFs, , ( 1,2 ∕ 1,2 ≡  1,2 ,   ) denoting -hadron FFs,  1,2 and  1,2 respectively standing for incoming-and fragmenting-partons' longitudinal fractions, and d σ representing the partonic cross section.Analogously, we write the following collinear-convolution formula for the -hadron plus light-jet case depicted in the right diagram of Fig. 1 d LO At variance with the collinear framework, to write the expression for the high-energy resummed cross section in the hybrid factorization, one first makes use of BFKL and then completes the picture by adding collinear PDFs and FFs.We suitably rewrite the differential cross section as a Fourier expansion in terms of azimuthal coefficients coefficients are calculated within the BFKL formalism at a full NLL accuracy.Making use of the MS renormalization scheme [212], we obtain [213] where ᾱ (  ) ≡   (  )  ∕ with   the color number,  0 = 11  ∕3 − 2  ∕3 is the first coefficient of the QCD -function with   the flavor number.A two-loop running-coupling choice with   (   ) = 0.118 and with dynamic thresholds for the flavor number   is made.The kernel entering the exponent of Eq. ( 7) encompasses the resummation of LL and NLL terms.Its analytic expression is  NLO (, ) = (, ) + ᾱ χ(, ) , (8) where stand for the LO kernel eigenvalues,  E is the Euler-Mascheroni constant, and () ≡ Γ ′ ()∕Γ() the Gamma-function logarithmic derivative.The χ(, ) function in Eq. ( 8) represents the NLO kernel correction in the Mellin space where  1,2⟂ are the transverse masses of the two outgoing states.For a -hadron one has 2  cosh 2 () ]} .
The Ξ(, ) function reads + ln and The two quantities stand for the emission functions of the forward -hadron and the light jet.At LO one has and with   ≡ ( 2  − 1)∕(2  ) and   ≡   being the Casimir factors for a gluon emission from a quark and a gluon, respectively.The  () function at the end of the last row of Eq. ( 7) embodies the logarithmic derivative of LO emission functions What is left in Eq. ( 7) are the NLO corrections to  , functions, 1,2 .The NLO correction to the  hadron at large transverse momentum was computed in Ref. [126] and is provided in the Appendix A. Our choice for the jet NLO emission function builds on analyses described in Refs.[126,128].It relies upon a jet algorithm obtained in the "smallcone" limitand with a cone-type distance [129,215,216], where the radius is set to  = 0.5 (see the Appendix B).
A proper "resummation versus fixed-order" study should rely upon comparing NLL predictions from our highenergy approach with NLO calculations in pure collinear factorization.Unfortunately, a numerical technology aimed at computing NLO distributions for two-particle reactions in hadron collisions is still missing.Thus, for the sake of comparison with reference fixed-order results, one can truncate the expansion of the azimuthal coefficients of Eq. ( 7) up to the ( 3  ) order.In this way, one obtains an effective high-energy fixed-order (HE-NLO + ) expression that collects the leading-power asymptotic signal present in a pure NLO calculation, whereas terms proportional by inverse powers of the partonic center-of-mass energy are neglected.The MS expressions for HE-NLO + azimuthal coefficients reads with the exponentiated kernel expanded and truncated at (  ).
We will also compare our NLL results with the pure LL background.The latter is obtained by discarding NLO corrections of both the kernel and the emission functions Eqs. (7) to (20) shed light on the structure of our hybrid factorization.According to BFKL, the cross section is high-energy factorized as a convolution between the Green's function and two singly off-shell emission functions.The latters embody collinear inputs, namely collinear convolutions between incoming protons' PDFs and outgoing hadrons' FFs.The NLL∕NLO + label tells us that the resummation of energy logarithms is performed within a full NLL accuracy and within NLO perturbative order.The '+' superscript in Eqs.(7) and (19) highlights that some Nextto-NLL (NNLL) extra terms come from the cross product between NLO corrections of the two emission functions.Finally, factorization (  ) and renormalization (  ) scales are set to the natural scales suggested by kinematics.In particular, we fix   =   =   ≡  1⟂ +  2⟂ .
Our guiding line here is to consider, as a natural scale, any energy scale with the same order of magnitude of the typical energy of the emitted particle.Since, in our case, two objects are emitted, there are in principle to scales (one for each fragmentation region).Natural values for them are  1⟂ and  2⟂ , respectively.For the sake of comparing our observables with future, possible predictions coming from other approaches, it is convenient to make a simple choice, namely to have just one scale.As a bonus, choosing the sum of  1⟂ and  2⟂ turns out to match the settings of many other numeric codes aimed at precision QCD phenomenology (see, e.g., Refs.[217][218][219]).

Collinear inputs
As for collinear PDFs, we make use of the NNPDF40_nlo_as_01180 NLO set [220,221] as implemented in LHAPDF v6.5.4 [222].It was obtained on the basis of global fits by means of the well-defined replica methodology originally proposed in Ref. [223] in the context of a neural-network approach (we refer to Ref. [224] for an in-depth analysis on ambiguities about correlations among different PDF sets).
Parton fragmentation to singly-bottomed   ≡ { + Λ  } hadrons is depicted via the KKSS07 NLO FF determination, which was originally obtained by fitting data on inclusive productions of noncharmed  mesons in lepton annihilation [56].Here, the threshold for bottom-and antibottom-quark evolution is set to 4.5 GeV and its initial-scale parametrization reads as a power-like function with three parameters [225].Light partons and the charm quarks are perturbatively generated with DGLAP, while their FFs are zero below the threshold for the bottom.The current KKSS07 set [58] for   particles was constructed from the original noncharmed -meson ones by switching off the branching fraction for the  →  ± subprocess, which was assumed to be   =   = 0.397 [56].Such a choice is justified by the assumption that a unique function can be employed in the description of noncharmed -meson fragmentation or of an inclusive sum of  and Λ  channels, but not for the Λ  -baryon case alone.
Our strategy for charmed -meson fragmentation is more involved.It was presented for the first time in our recent study focused on accessing the high-energy spectrum of QCD via the hadroproduction of these doubly heavy-flavored states [84].It takes inspiration from a first determination of VFNS, DGLAP-evolving FFs for vector quarkonia [85,86].The starting point is the NLO calculation, within the NRQCD effective theory, of initial-scale inputs for charm, bottom quark [82] and gluon [83] fragmentation channels to   ( 1  0 ) ≡   and   ( 3  1 ) ≡  *  mesons (see Refs. [226,227] for similar calculations).To generate DGLAP-evolved sets for our charmed  mesons we need to plug the NRQCD inputs into a numeric evolution code.Among them, we cite: QCD-PEGASUS [228], HOPPET [229], QCDNUM [230], APFEL(++) [231][232][233], and EKO [234,235].Contrariwise to PDFs, which obey a space-like evolution equation, FFs must be time-like evolved [236,237].In Ref. [84] we employed APFEL++ to generate novel and phenomenology-ready LHAPDF FF grids, which we named ZCFW22 NLO determinations.
For the sake of comparison, we employ the NNFF10_KAsum_nlo set [238] (see Ref. [239] for advancements) to describe charged-kaon fragmentation at NLO.

Natural stabilization of the NLL series
Here we provide details on the stabilizing mechanism emerging from the VFNS collinear fragmentation of hadrons.The linking dynamics between heavy-flavor FFs and the observed stability of semi-inclusive distributions calculated within the hybrid collinear and high-energy factorization at NLL∕NLO (+) was first unveiled through a study of semi-hard emissions of heavy-light hadrons, namely:  particles [240][241][242][243][244][245][246][247][248], Λ  baryons [241,249,250], and   states [55-58, 60, 251-254].This striking property came out for the first time from an analysis of semi-inclusive detections of Λ  [91,164],  * ± [166,167], and   [91,165] hadrons in forward directions of rapidity at the LHC.Such a remarkable result was then corroborated by similar investigations on vector quarkonia [85,255] and  ( * ) mesons [84], whose formation mechanism was modeled on the basis of single-parton NRQCD fragmentation [82,83,226,227,[256][257][258][259][260][261].A further evidence came recently out in the context of heavy-light tetraquarks [210,262] described via a Suzuki-Nejad-Amiri-Ji initial-scale fragmentation input [263][264][265][266] (see Refs. [267,268] for further details).A less pronounced stabilization pattern was also observed in the case of strange-quark flavored, cascade Ξ − ∕ Ξ+ baryon tags [161].Panels in Fig. 2 contain the factorization-scale dependence of KKSS07 [56,58] and ZCFW22 [84] NLO functions respectively describing the collinear fragmentation of   hadrons (right upper plot) and  ( * )  mesons (lower plots).For comparison, the left upper plot carries analogous information of light-quark flavored species, namely charged kaons, described by the NNFF1.0[238] NLO collinear-FF determination.For the sake of simplicity, we show values of our FF sets just for one value of the final-state hadron longitudinal-momentum fraction, .In particular, we choose a value of  that roughly matches the mean value at which FFs are typically probed standard semi-hard configurations (see Section 3.2).Thus we have  = 0.475 ≃ ⟨⟩.As expected, the  channel of KKSS07 and ZCFW22 functions heavily prevails over the other flavors.Particular attention deserves, however, the behavior the gluon FF, which clearly decreases with   for kaons but smoothly increases for all the three -flavored species, up to reach a plateau.More quantitative information on the FFs of Fig. 2 would come from an analysis on the associated uncertainties.However, apart the kaon NNFF1.0 determination, the other ones do not carry such information.KKSS07 functions  describing noncharmed -hadron fragmentation were extracted without any uncertainty estimate, whereas the inclusion of uncertainties for the two ZCFW22 sets, possibly connected with perturbative scale-variation effects, is planned as a future extension of this work.
As extensively shown in Refs.[84,85,112,164,165,210], the gluon collinear FF plays a crucial role.De facto, it has a strong stabilization power on our NLL∕NLO + observables.Its behavior with   fairly controls the stability of the high-energy logarithmic series of our distributions.In the semi-hard regime proton PDFs are typically probed at 10 −4 ≲  ≲ 10 −2 , where the gluon PDF heavily dominates over quark ones.Given that the gluon FF diagonally multiplies with the gluon PDF in the LO hadron emission function (see Eq. ( 16)), its impact on the cross section is magnified.This effect is not quenched at NLO, where {} and {} nondiagonal channels are opened [164] (see the Appendix A).On one side,   decreases with   both in the Green's function and in the emission functions (see Section 2.2).At the same time, however, the gluon PDF increases as   enlarges.When that density multiplies a gluon FF which also increases with   , as it happens for -hadrons, the two effects balance each other.This leads to the stabilizing trend of heavy-flavor distributions under energy-scale variations.If instead the gluon FF falls down with   , as it happens for light-flavored species, such charged kaons, the stabilizing pattern is lost.This harms any chance of making precision tests of semi-hard observables at the natural energies given by kinematics [88].The emergence of the natural stability [112] in presence of both singly and doubly heavy flavored hadron species provides us with clear evidence that this striking feature is an intrinsic property of heavy-flavor detections, which becomes manifest whenever a heavy hadron emitted, independently of the choice made for the initial-scale input of the corresponding collinear fragmentation.

Hi-Lumi LHC distributions
We performed our numeric computations via the JETHAD v0.5.1 interface [86,88,91], a PYTHON+FORTRAN multimodular working environment aimed at the calculation, management, and processing of high-energy distributions obtained from different theoretical approaches.In particular, rapidity and transverse-momentum differential cross sections were calculated by means of JETHAD FORTRAN 2008 core modules, whereas its PYTHON 3.0 analyzer was used to elaborate results.Section 3.1 gives details on our procedure to gauge the size of uncertainties for our resummed observables.Final-state kinematic cuts are presented in Section 3.2.Then, Section 3.3 carries a discussion of our NLL∕NLO + predictions.

Uncertainty estimation
A widely-used methodology to quantify the weight of Missing Higher-Order Uncertainties (MHOUs) relies upon assessing how much our observables are sensitive to variations of factorization and renormalization scales around their natural values.It is well-known that the effect MHOUs gives perhaps the major contribution to the overall uncertainty [91].For this reason, we will gauge the size of simultaneously varying   and   around   ∕2 and 2  .The   parameter in panels of figures of Sections 3.3 and 3.4 is just the   ≡   ∕  =   ∕  ratio.This will represent our standard procedure to estimate the size of MHOUs.Then, for the sake of comparison with other approaches, where energy scales are allowed to assume values much larger than the natural ones (see, e.g., the BLM method [201][202][203]), we will also present an accessory, extended MHOU scan, with the   parameter ranging from 1 to 30.Another possibly relevant uncertainty is encoded in proton PDFs.Quite recent analyses on high-energy production rates have shown, however, that selecting distinct PDF parametrizations as well as distinct members inside the same set produces a very small effect [88,91,157,165].Therefore, we will calculate our observables by considering only the central member of the NNPDF4.0parametrization.Then, other potential uncertainties could come from a collinear improvement of the NLO kernel (see Refs. [269][270][271] and reference therein), which consists in including renormalization-group (RG) terms to make the BFKL equation compatible with the DGLAP one in the collinear limit, or from changing the renormalization scheme.The first case was deeply addressed in Ref. [91].As a result, the impact of collinear-improvement techniques on rapidity-differential rates is fairly contained inside error bands produced by MHOUs, and it is even smaller in the case of when other differential observables.The same article provides us with an overestimation of the effect of passing from MS [212] to MOM [272,273] renormalization scheme.MOM results for rapidity distributions are systematically higher than MS ones, but still contained into MHOU bands.We stress, however, that a proper MOM analysis should be based on MOM-evolved PDFs and FFs, not available so far.
We will produce uncertainty bands for our distributions by combining MHOUs with the numeric error generated by the multidimensional integration (see Section 3.2).The latter will be constantly kept below 1% by the JETHAD integration subinterfaces.

Kinematic constraints on final-state ranges
We will consider two main observables: () the rapidity-interval distribution, namely the cross section differential in the rapidity separation, Δ =  1 −  2 , between the two produced particles, and () the doubly-differential transversemomentum distribution.
The first observable genuinely corresponds the  0 azimuthal-angle coefficient, defined in Section 2.2, integrated over final-state transverse momenta and rapidities and taken at fixed values of Δ d(Δ ; ) where the '[order]' label for  0 inclusively refers to: NLL∕NLO + , HE-NLO + , or LL∕LO.Here we made use of a (Δ −  1 +  2 ) function to remove the integration over one of the two rapidities, say  2 .Transverse momenta of the forward -hadron and the backward one stay in the ranges 60 < | ⃗  1 |∕GeV < 120 and 30 < | ⃗  2 |∕GeV < 120, respectively.These cuts fairly meet the criteria of applicability of a VFNS-based fragmentation approach, in which energy scales must be sufficiently higher than thresholds for the DGLAP evolution of heavy-quark species.Then, light jets are tagged in transverse-momentum ranges lightly different but compatible with current and forthcoming studies at the LHC [274], say 50 < | ⃗  2 |∕GeV < 60.The use of such asymmetric transverse-momentum ranges eases disentangling the pure high-energy dynamics from the fixed-order signal [88,144,145].It also dampens large Sudakov logarithms due to almost back-to-back configurations which would require the employment of another resummation [275][276][277][278][279][280].Finally, it quenches instabilities connected to radiative corrections [281,282] as well as violations of the energy-momentum conservation relation [283].As for the rapidity ranges, we select typical cuts of current LHC studies, where hadrons are detected only in the barrel calorimeter [284], say from −2.4 to 2.4, whereas the jet can be also tagged in endcaps' regions [274], say from −4.7 to 4.7.
The second observable is given by the  0 azimuthal-angle coefficient, integrated over rapidities while Δ is kept fixed, but this time being differential in the transverse momenta of the two outgoing objects Here we consider the previously mentioned kinematic cuts for rapidities, whereas the two transverse momenta can run from 10 to 100 GeV.This observable was recently proposed, in the context of singly bottom-flavored hadrons [165] and cascade Ξ baryons [161], as a favorable testing ground whereby unraveling the interplay among distinct resummations.

Rapidity-interval distributions
In Fig. 3 we present rapidity-interval distributions for our reference processes, calculated within NLL∕NLO + hybrid factorization and taken at √  = 14 TeV.Left plots refer to double-hadron channels, while right plots are for the hadronjet ones.Final states with   hadrons,   ( 1  0 ) mesons, or   ( 3  1 ) resonances are considered in upper, central, and lower plots, respectively.More in particular, any forward -flavored hadron is always accompanied by a singly bottomed hadron   (left plots) or by a light-flavored jet (right plots).This is to avoid the presence of two identified charmed  mesons in the final state, which would lead to a substantial lowering in statistics.We will adopt this choice also in the next Section (3.4).The minimum values of our distributions are the ones obtained at the maximum values of Δ .They are uniformly larger than 10 −2 nb when only   hadrons (and light jets) are produced, and larger than 10 −5 nb in the two  ( * )   cases.This represents a very promising statistics that can be easily accessed via dedicated experimental studies at the forthcoming Hi-Lumi upgrade of the LHC.The decreasing trend with Δ of our NLL-resummed distributions is a general feature shared by all the semi-hard hadroproduction final states investigated so far.It comes out as the interplay between two competing effects.Indeed, although the high-energy dynamics encoded in the partonic hard factor would make the cross section grow with energy, collinear PDFs and FFs fall off very fast with Δ , since this kinematically translates to longitudinal-momentum fractions very close to one.The net result is a dampening effect with Δ of hadronic distributions.
The core analysis shown in Fig. 3 is a study of our rapidity-interval differential distributions under a progressive variation of factorization and renormalization scales in an broad window, controlled by the   parameter introduced in Section 3.1, which ranges from 1 to 30.This strategy actually represents an extended MHOU scan with respect to the traditional one, 1∕2 <   < 2. Ancillary panels drawn right below primary ones show reduced distributions, namely the ones divided by their central values calculated at   = 1.Remarkably, NLL∕NLO + results for all the six considered final states feature quite a weak sensitivity on   variation in the whole spectrum of Δ matter of our investigation.This supports the important message that a strong stabilization mechanism is encoded both in the KKSS07 functions extracted from global data on  mesons and Λ  baryons, as well as in the novel ZCFW22 FFs built on the basis of NRQCD initial-scale inputs for   ( 1  0 ) and   ( 3  1 ) generalized-quarkonium states and then evolved via DGLAP in a VFNS scheme.Results contained in plots of Fig. 3 bring to light the emergence of the natural stability [112] from semi-inclusive -flavor emissions at high energies.
As a complement, in Fig. 4 we show results for the same Δ -distributions, calculated within the NLL∕NLO + accuracy and compared with HE-NLO + and LL∕LO corresponding ones.Left (right) plots are for hadron-plus-hadron (hadron-plus-jet) observables.Upper, central, and lower plots respectively correspond to   ,   , and  *  sensitive final states.This time, shaded bands are built by combining the uncertainty of the multidimensional integration, numerically performed by the JETHAD integrators, with the one coming from a traditional MHOU scan, namely with the   parameter running from 1/2 to two.As done before, lower ancillary panels exhibit reduced distributions, i.e. divided by central values calculated when   = 1.The global outcome is a strong stability under higher-order corrections.This is confirmed by the fact that NLL∕NLO + uncertainty bands are generally narrower and partially nested inside corresponding LL∕LO ones, in particular in the large-Δ regime, where the impact of high-energy logarithms is expected to be large.Their overlap is slightly weaker for hadron-plus-jet cases (right plots).This is not unexpected, since the emission of just one -hadron rather then two reduces the stabilizing power by one half.
As discussed in previous works (see, e.g.Section 3 of Ref. [84]), shaded bands for fixed-order HE-NLO + results are almost completely nested in resummed NLL∕NLO + ones.Future analyses on more exclusive observables, such as angularity distributions will shed light on a clearer disentanglement between our resummation and the fixed-order signal.Finally, as a bonus result, we confirm the estimate, provided by the LHCb Collaboration [70,71], on the production-rate hierarchy between singly-bottomed hadrons and charmed  mesons.Production rates of the latters seems not to exceed 0.1% of the   ones, at most.To compare our predictions with these numbers, we can consider the hadron-plus-jet channel (right diagram of Fig. 1) at LL∕LO, whose differential cross section is linear in the hadron FFs (see Eq. ( 16)).The inspection of right plots of Fig. 4 in clearly in line with LHCb estimates for charmed -meson sates. 3The hierarchy remains de facto unaltered both at NLL∕NLO + and at HE-NLO + .Remarkably, this fact not only represents a reliability test for our hybrid factorization, but it also definitely corroborates the validity of the single-parton fragmentation mechanism from NRQCD, at least for the   ( 1  0 ) case.
Combining threshold and high-energy logarithms is a cumbersome task.A first double resummation was obtained in the context of the inclusive Higgs hadroproduction from gluon-gluon fusion [95].The key ingredient there was the possibility of decoupling the two dynamics in the Mellin space, so that the small- (high-energy) resummation and the large- (threshold) one respectively "control" the small- and large- tail of the the Mellin projection of the cross section [305,306].In this way, the two kinds of logarithms can be separately resummed.In our case, however, the semi-inclusive emission of two objects (see Fig. 1) genuinely leads to rapidity-differential observables, for which a theoretical framework to properly perform such a double resummation is still missing.
Inclusive emissions of two-particle final states in proton collisions, such as photon [318][319][320][321], Higgs [322] and  ± boson [323] pairs, as well as boson plus jet [324,325] and  plus photon [326] systems, were proposed as favorable channels whereby testing the TM resummation.TM third-order fiducial rates for Drell-Yan and Higgs tags were respectively investigated in Refs.[327][328][329][330] and [328,[331][332][333]. Ref. [324] presents a pioneering joint resummation of TM logarithms rising from the emission of two states.There, the simultaneous measurement of transverse momenta of the Higgs boson and the leading jet was analyzed up to the NNLL accuracy by means of the RADISH numeric technology [331].Ref. [334] deals with a study of TM distributions for the fully leptonic ( +  − ) tag at the LHC.
Heavy-flavor distributions bring to a further complication.Indeed, when the transverse momentum of a heavy hadron diminishes, its transverse mass comes closer and closer to heavy-quark masses serving as DGLAP thresholds, even crossing then when energy scales are moved below their natural values (as an example, this happens when the   parameter introduced in Section 3.1 is set to 1/2).In this case, the validity of a pure VFNS approach becomes questionable.
In this Section we consider cross sections at fixed values of Δ and differential in the transverse momenta of both the final-state objects.In the kinematic sectors matter of our investigations energy logarithms arising from the semi-hard scale hierarchy are large.At the same time, however, also TM logarithms are potentially sizable.Therefore, we propose this study without claiming to catch the core features of these observables thanks to our hybrid-factorization framework, but rather to set the stage for futures analyses aimed at shedding light on the formal and phenomenological distributions (right panels) instead generally oscillate around   = 1, which turns out to work as a critical point.This is a clear signal that our NLL∕NLO + doubly-differential cross sections are quite stable under MHOU analyses.Notably, no distribution comes with a peak, which might however be present in the low-TM range (excluded from our study), where TM-resummation effects would dominate.The information contained in the core of our 3D figures is supplemented by 2D slice projections highlighting the behavior of our observables at | ⃗  1 | = 0 and at | ⃗  2 | = 0. What emerges from the inspection of these 2D contour plots at low or intermediate transverse momentum is a smaller

statistics in the |
range than in the opposite case.This happens because cross sections are globally smaller when a bottomed hadron is produced rather than a light-flavored jet [84,165].
Finally, numeric checks on our transverse-momentum distributions have shown a noncharmed versus charmed -meson production-rate hierarchy analogous to the one found in the case of rapidity-interval distributions (see Section 3.3).  production rates are almost always three orders of magnitude smaller than   ones, whereas rates of the  *  resonance range from 10 −2 to 1∕60 times the ones of   hadrons.

Summary and Outlooks
We considered the semi-inclusive production of a forward noncharmed or charmed -hadron in association with a backward noncharmed -hadron or a jet in high-energy regimes that can be accessed at current LHC analyses and its forthcoming Hi-Lumi upgrade.We depicted the production of bottom-flavored hadrons by means of a VFNS collinearfragmentation approach, whose validity is justified by the large transverse-momentum regime matter of our analysis.
Singly heavy-flavored particles, namely inclusive states consisting of  mesons and Λ  baryons, were described via the KKSS07 [56,58] NLO FFs based on a fit to global data and assuming a three-parameter power-like initialscale parametrization for  and b fragmentation channels [225].As for charmed  mesons, i.e.   ≡   ( 1  0 ) states and  *  ≡   ( 3  1 ) resonances, we adopted a more sophisticated strategy.First, we modeled the initial-scale input for the fragmentation of heavy quarks (charm, bottom, and their antiparticles) and gluons to observed hadrons by means of corresponding NLO calculations within the NRQCD effective theory [226,227].Then, we plugged these inputs into the APFEL(++) DGLAP evolution code [231][232][233] to generate phenomenology-ready VFNS FF grids in LHAPDF format [222], which we named ZCFW22 NLO functions [84].These novel FF sets will play a role in future phenomenological applications lying outside the semi-hard regime.They also will serve as a support basis to shed light on the intersection corner between the collinear factorization and the NRQCD formalism.
By relying upon a traditional as well as an extended MHOU scan on our hybrid-factorization setup, we hunted for clues of stabilization of the high-energy resummation under next-to-leading logarithmic corrections.We provided strong evidence that these effects are present and they allow for a fair description of high-energy distributions, differential in the rapidity interval and in the final-state transverse momenta.As a bonus, we highlighted that the predicted production-rate hierarchy between noncharmed -hadrons and charmed   ( 1  0 ) mesons is in line with recent LHCb estimates.In particular, the production rate of   states is three orders of magnitude lower than the   one.This served as a simultaneous benchmark for both the hybrid factorization and the NRQCD fragmentation applied to charmed  mesons at our reference transverse masses (see Section 3.2).
As a prospect, we plan to investigate bottom-flavor emissions in the moderate transverse-momentum sector, namely at the intersection line between the VFNS and the FFNS scheme, and possibly set the stage for a matching between the two descriptions.Here, high-energy tags of charmed  mesons will bring novelty to the plan.In that case, indeed, the VFNS versus FFNS dichotomy translates into a short-distance versus fragmentation dichotomy.This calls for a further enhancement of our ZCFW22-based treatment, which currently relies upon NRQCD-modeled calculations for initial-scale inputs, but will serve in the medium-term future as a "launch pad" for extractions through fits to new data to be collected at the Hi-Lumi LHC.In this respect, neural-network strategies already employed to extract collinear FFs for lighter hadron species will be an asset [238,239,[362][363][364][365].As mentioned in Section 2.4, a proper determination of FF uncertainties with possible inclusion of MHOU effects, not yet encoded in our ZCFW22 sets, is needed.It might take inspiration from the MCscales approach [366] (for close-in-spirit studies on proton and pion PDFs, see Refs.[367][368][369][370] and Ref. [371]).
The fair stability exhibited by our NLL∕NLO + distributions motivates our interest in proposing the hybrid factorization as useful tool to enhance and extend the standard collinear description.To reach the precision level, our NLL∕NLO + hybrid approach needs to evolve into a multi-lateral formalism where different resummation mechanisms are at work.As a first steps, links with the soft-gluon [279,280,372,373] and the jet-radius [374][375][376][377][378] 4 resummations should be found.Possible common grounds with studies on jet angularities [379][380][381] also represent an intriguing perspective.
Then, the resummation of energy logarithms from BFKL could improve the description of DGLAP-evolved bottomflavor FFs, its effect becoming sizable already at moderate values of final-state hadron longitudinal fraction [382].Another engaging development would be exploring the validity of the use of a ZM-VFNS with higher HFMPs in the context of heavy-meson fragmentation [21].Finally, the striking evidence of intrinsic charm quarks in the proton [383], now corroborated by quite recent studies on its valence distribution [384], traces the path towards searches for imprints of the possible existence of an intrinsic bottom component.In this respect, any progress in widening the horizons of our knowledge of bottom physics is needed.We believe that studies aimed at shedding light on the bottom-flavor collinear fragmentation go along these directions.
d  ] , In our expressions  0 is an energy scale that genuinely appears within the BFKL approach.We fix  0 ≡   .Then we define τ = 1 −  and  =  − with 2  1 being the standard hypergeometric function.
The plus prescription entering Eqs. ( 24) and ( 27) comes out as where  () is a test function regularly behaved when  approaches one.

Appendix B: NLO light-jet emission function
We present here the NLO correction to the emission function for a forward light-flavored jet emission in the smallcone selection (see Refs. [129,213] 5.62  GeV or   ≡   ( * )  = 6.275GeV.Conversely, the transverse mass of a light-flavored jets can be safely set to transverse momentum, say  ⟂ ≡ | ⃗   |.The NLO χ(, ) function was computed in Ref.[214]