Unwinding the rare $\rm\Omega$ sector: Fragmentation of
fully charmed baryons from HL-LHC to FCC

The fragmentation of heavy-flavored hadrons is inherently more intricate than that of light hadrons. This complexity arises because heavy-quark masses in their lowest Fock states fall within the perturbative QCD domain. As a result, while light-hadron FFs encode only nonperturbative dynamics at the initial energy scale, heavy-hadron FFs require a combination of perturbative and nonperturbative elements.

For singly heavy hadrons such as $D$ mesons, $B$ mesons, and ${\rm\Lambda}_{Q}$ baryons, the fragmentation process unfolds in two distinct stages Cacciari et al. (1997); Cacciari and Greco (1994); Jaffe and Randall (1994); Kniehl et al. (2005); Helenius and Paukkunen (2018, 2023). The first stage involves a parton $i$, produced in a hard-scattering event with large transverse momentum, fragmenting into a heavy quark $Q$. Since $\alpha_{s}(m_{Q})<1$, this step can be computed within perturbative QCD.

This contribution, known as the short-distance coefficient (SDC) for the $[i\to Q]$ fragmentation process, occurs on a much shorter timescale than hadronization. The first NLO determination of SDCs for singly heavy hadrons was presented in Ref. Mele and Nason (1990, 1991), with subsequent next-to-NLO refinements provided in Refs. Rijken and van Neerven (1996); Mitov and Moch (2006); Blumlein and Ravindran (2006); Melnikov and Mitov (2004); Mitov (2005); Biello and Bonino (2024).

At later timescales, the heavy quark $Q$ transitions into a bound hadronic state. This second stage of fragmentation is entirely nonperturbative and is commonly modeled through phenomenological approaches Kartvelishvili et al. (1978); Bowler (1981); Peterson et al. (1983); Andersson et al. (1983); Collins and Spiller (1985); Colangelo and Nason (1992) or effective field theories Georgi (1990); Eichten and Hill (1990); Grinstein (1992); Neubert (1994); Jaffe and Randall (1994).

In mathematical terms, for a parton $i$ fragmenting into a singly heavy hadron ${\cal H}_{Q}$ at the initial scale $\mu_{F,0}$ of the order of the mass of the heavy-quark, $m_{Q}$, one has Cacciari et al. (1997); Cacciari and Greco (1997)

$$D_{i}^{{\cal H}_{Q}}(z,\mu_{F,0})=\int_{z}^{1}\frac{{\rm d}\xi}{\xi}D_{i}^{Q}(\xi,\mu_{F,0})\,D_{\rm[np]}^{{\cal H}_{Q}}\left(\frac{z}{\xi}\right)\;.$$

(1)

Here, $D_{i}^{Q}$ is the perturbative part of the initial scale FF, namely the SDC for an outgoing (massless) parton to fragment, via a perturbative QCD cascade, into the constituent (massive) heavy quark. Conversely, $D_{\rm[np]}^{{\cal H}_{Q}}$ depicts the nonperturbative part of the fragmentation, which is taken to be universal, independent of the generating parton $i$. It also does not depend on $\mu_{F,0}$.

To construct a fully evolved set of FFs within the VFNS, it is crucial to account for energy evolution effects. Starting from initial-scale inputs, assumed to be free from scaling violations, the functions are evolved numerically using DGLAP timelike evolution equations at the required level of perturbative accuracy.

The two-step initial-scale fragmentation framework, originally devised for singly heavy hadrons, can be extended to quarkonia. In this case, the simultaneous presence of a heavy quark $Q$ and its antiquark $\bar{Q}$ within the lowest Fock state $|Q\bar{Q}\rangle$ adds an additional layer of complexity to the fragmentation description. The modern theoretical treatment of quarkonium production is based on nonrelativistic QCD (NRQCD) Caswell and Lepage (1986); Thacker and Lepage (1991); Bodwin et al. (1995); Cho and Leibovich (1996a, b); Leibovich (1997); Bodwin et al. (2005) (for a pedagogical introduction, see Refs. Grinstein (2000); Krämer (2001); Brambilla et al. (2004); Lansberg (2005, 2020)).

NRQCD treats heavy-quark and antiquark fields as nonrelativistic degrees of freedom, enabling a systematic factorization between SDCs, which govern the perturbative production of the $[Q\bar{Q}]$ pair, and long-distance matrix elements (LDMEs), which encapsulate the nonperturbative hadronization dynamics. A physical quarkonium state is then expressed as a linear superposition of all possible Fock states, organized through a systematic expansion in both the strong coupling $\alpha_{s}$ and the relative velocity $v_{\cal Q}$ of the heavy quark-antiquark pair.

A key feature of NRQCD is its ability to describe quarkonium production mechanisms across a wide range of transverse momenta. At low transverse momentum, $q_{T}$, the dominant process involves the short-distance formation of the $[Q\bar{Q}]$ pair in the hard scattering, followed by its nonperturbative evolution into a physical bound state. At higher $q_{T}$, however, an alternative mechanism—where a single parton fragments into the quarkonium state plus additional radiation—begins to compete with, and eventually surpasses, the short-distance process.

From a theoretical perspective, short-distance quarkonium production can be interpreted as a fixed-flavor number scheme (FFNS) Alekhin et al. (2010) two-parton fragmentation, capturing higher-twist power corrections (see Refs. Fleming et al. (2012); Kang et al. (2014); Echevarria (2019); Boer et al. (2023); Celiberto (2024a, 2025e, 2025f) for more details). By contrast, the single-parton mechanism at large $q_{T}$ is a collinear VFNS fragmentation, its energy evolution being governed by the DGLAP equations.

The first LO calculations of the initial-scale inputs for both gluon and heavy-quark FFs of $S$-wave color-singlet vector quarkonia were performed in the early 1990s Braaten and Yuan (1993); Braaten et al. (1993a). However, corresponding NLO refinements became available only recently Zheng et al. (2019a, 2022).

Building upon these inputs, a pioneering set of VFNS, DGLAP-evolving FFs for vector quarkonia, ZCW19⁺, was introduced in Refs. Celiberto and Fucilla (2022a); Celiberto (2023a). Shortly thereafter, the ZCFW22 extension was developed to incorporate $B_{c}$ mesons Celiberto (2022b, 2024d). Here, high-energy resummed production rates of charmed $B$ mesons, derived using the ZCFW22 FF framework, reinforced the experimental observation by the LHCb Collaboration Aaij et al. (2015b, 2017); Celiberto (2024d) that the relative production rate of $B_{c}(^{1}S_{0})$ mesons compared to singly bottomed $B$ mesons remains below 0.1% Celiberto (2024d). This result simultaneously benchmarked for both the high-energy resummation framework and the NRQCD fragmentation applied to $B_{c}$ mesons.

Remarkably, recent studies indicate that NRQCD factorization can be suitably extended to explore the production mechanism of di-$J/\psi$ excitations Aaij et al. (2020); Aad et al. (2023); Hayrapetyan et al. (2024b), interpreting these resonances as fully charmed, exotic tetraquarks Zhang and Ma (2020); Zhu (2021). From a NRQCD viewpoint, the production of a fully heavy $T_{4Q}$ tetraquark originates from the generation of a heavy-quark pair and its corresponding antiquark pair at short distances, with a characteristic energy scale set by the inverse of the heavy-quark mass. As with singly heavy hadrons and quarkonia, asymptotic freedom allows for the usual heavy-flavor fragmentation mechanism in two steps: an initial perturbative stage followed by a nonperturbative hadronization phase.

The first NRQCD calculation of the initial-scale FF for the $[g\to T_{4Q}]$ $S$-wave channel in a color-singlet configuration appeared in Ref. Feng et al. (2022a). Building on this foundation, the recent TQ4Q1.0 VFNS sets combined this NRQCD-based gluon channel with an initial-scale input for the $[Q\to T_{4Q}]$ function, derived by adapting a Suzuki-model approach Suzuki (1977, 1986); Amiri and Ji (1987); Moosavi Nejad and Amiri (2022).

Expanding this methodology, the first VFNS FFs for heavy-light tetraquarks, named TQHL1.0 determinations, were derived in Refs. Celiberto and Papa (2024); Celiberto (2024b). The latest release, the TQ4Q1.1 and TQHL1.1 families Celiberto and Gatto (2025); Celiberto (2025a), further refined the description by incorporating NRQCD-based modeling of also the $[Q\to T_{4Q}]$ initial-scale FF Bai et al. (2024), and improving the treatment of $X_{Qq\bar{Q}\bar{q}}$ fragmentation.

Finally, Ref. Celiberto (2025d) accompanies and presents the first release of collinear FFs for fully charmed $|ccc\bar{c}c\rangle$ pentaquarks. The resulting set PQ5Q1.0 embodies a consistent heavy-quark threshold DGLAP evolution from two possible initial-scale inputs for the heavy-quark channel: a compact direct multicharm state Farashaeian and Moosavi Nejad (2024a) or a dicharm-charm-dicharm configuration Farashaeian and Moosavi Nejad (2024b).

The quark-diquark model provides a simplified, yet powerful framework for describing the internal structure of baryons, wherein two quarks are treated as a tightly bound subsystem, referred to as a diquark, that interacts with the third quark Gell-Mann (1964). This approach has been widely used in hadron spectroscopy and in the analysis of production and decay mechanisms of baryons, particularly those containing heavy quarks (see Refs. Maiani et al. (2005); Jaffe and Wilczek (2003); Guo et al. (2013); De Sanctis et al. (2016) for related discussions). Within this picture, diquarks are modeled either as scalar objects with spin-0 or as axial-vector states with spin-1. Their internal structure, inherently nonperturbative, is typically encoded through phenomenological form factors.

Scalar diquarks involve a single form factor and are generally associated with simpler spin configurations, while axial-vector diquarks require multiple form factors to describe their richer internal dynamics. Both types have been utilized in the modeling of spin-dependent processes in baryon formation, and also in the parameterization of polarized quark and gluon distributions inside the nucleon, especially within spectator models Bacchetta et al. (2008, 2010, 2020, 2024a); Chakrabarti et al. (2023); Banu et al. (2024).

Early implementations of the quark-diquark framework to describe the fragmentation of light and heavy baryons can be found in Refs. Nzar and Hoodbhoy (1995); Ma et al. (2002); Yang (2002) and Falk et al. (1994); Adamov and Goldstein (1997); Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2019), respectively. The same formalism has also been extended to the analysis of exotic hadrons, such as pentaquarks Maiani et al. (2015), and to the mass spectroscopy of doubly and fully heavy tetraquarks using relativistic quasipotential approaches based on diquark-antidiquark interactions Faustov et al. (2020, 2021, 2022).

In the present work, we adopt a modeling strategy that has become standard in the literature Adamov and Goldstein (1997); Martynenko and Saleev (1996); Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2019), in which the diquark is considered to be a scalar constituent. This choice is primarily motivated by computational convenience. It simplifies the spin algebra, reduces the number of independent form factors, and results in a more manageable analytical structure. It has been widely used in the modeling of unpolarized $J=1/2$ baryons in different approaches. For example, in Adamov and Goldstein (1997), the scalar diquark model is used to calculate collinear FFs within a perturbative QCD framework, emphasizing its tractability. Similarly, the works in Refs. Martynenko and Saleev (1996); Gomshi Nobary and Sepahvand (2007) show that adopting a scalar diquark allows one to construct models that are consistent with available data, while keeping the formalism analytically viable.

In contrast, the use of axial-vector diquarks introduces greater theoretical complexity, including the need for additional form factors and a more elaborate treatment of spin correlations, which limits their use in practical calculations. Nevertheless, axial-vector diquarks are essential for describing baryons with spin-$3/2$ and for incorporating polarization effects in both spin-$1/2$ and spin-$3/2$ states. In particular, any realistic modeling of the polarized fragmentation or production of the $\rm\Omega_{3c}^{*}$ resonance, with $J=3/2$, requires a vector diquark configuration.

The scalar diquark approximation employed in this study should therefore be interpreted as an effective model for initiating a perturbative QCD-based description of fragmentation into unpolarized triply heavy baryons such as ${\rm\Omega}_{3c}$ states. More sophisticated treatments incorporating spin-1 diquarks could be pursued in future work to provide a more complete account of polarization and spin structure.

Interestingly, our adoption of the scalar diquark picture in the description of unpolarized ${\rm\Omega}_{3c}$ and ${\rm\Omega}_{3b}$ production is fully consistent with the color and spin symmetry requirements of the baryonic wavefunction, provided that the process is described within the two-step fragmentation mechanism. In this framework, the heavy quark first fragments into a color-antitriplet heavy diquark, which subsequently hadronizes into the triply heavy baryon through nonperturbative QCD effects.

Although scalar diquarks alone do not manifest the full symmetrization structure expected for a $J=1/2$ baryon composed of three identical fermions, this structure can be effectively recovered through hadronization effects encoded in the purely nonperturbative component of the FF. This interpretation, while not explicitly stated in Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2020), is in line with the modeling strategies adopted therein, where scalar diquarks are used in the indirect channel to compute FFs into ${\rm\Omega}_{3c}$ and ${\rm\Omega}_{3b}$, yielding consistent quantum number assignments and phenomenologically viable results.

The FF initial-scale inputs used in this study are extracted from existing calculations in the literature, which focus on the production of triply heavy baryons in a color-singlet configuration. Specifically, we rely on perturbative results for the quark-induced channel computed at LO in Ref. Moosavi Nejad and Delpasand (2017) and at NLO in Ref. Moosavi Nejad (2017), while for the gluon-induced channel we use the NLO result from Ref. Delpasand and Moosavi Nejad (2019). In all cases, the baryon is modeled as a color-singlet state formed from three charms, and the diquark is treated as a spin-0 scalar constituent.

Two main production scenarios are typically considered for phenomenological explorations. The first is the direct fragmentation mechanism, in which the heavy quark fragments directly into the three-quark baryon without invoking any internal substructure. This approach involves three-body diagrams and includes the direct interaction among the three constituent quarks. The second scenario, which we adopt here, is the diquark fragmentation model. In this framework, two of the three quarks are assumed to form a tightly bound scalar diquark, which then combines with the remaining quark to form the baryon. As mentioned before, this approximation significantly reduces the complexity of the calculation and has been widely used to describe the unpolarized production of $J=1/2$ heavy baryons.

The so-called indirect mechanism, as referred to in Refs. Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2019), corresponds precisely to the two-stage fragmentation framework discussed in Section 2.1, adapted to the diquark picture. From a theoretical standpoint, this structure offers a robust and systematically extensible framework for modeling the fragmentation of heavy-flavored hadrons. First, it separates the perturbative production of an intermediate colored object (the diquark) from its nonperturbative hadronization into the physical baryon. Second, both the quark and gluon FFs have been calculated at NLO within this setup, making it a consistent and reliable foundation for a VFNS analysis.

The Suzuki model Suzuki (1986) is among the ingredients of the framework adopted in Refs. Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2019). In that context, it plays a role analogous to that of NRQCD in quarkonium production. Originally developed to describe the fragmentation of heavy quarks into baryons or mesons, the Suzuki picture provides both a structural and a computational foundation for the two-stage fragmentation process. It defines a perturbative FF transition amplitude, where the constituent heavy quark fragments into a constituent diquark and a spectator, followed by the hadronization of the diquark into a baryonic bound state. The model incorporates kinematic simplifications, such as the collinearity of final-state momenta, and introduces non-perturbative elements through the wave function at the origin and an effective diquark form factor.

The Suzuki model is used in both direct and diquark fragmentation scenarios studied in the literature, although its role is far more central in the latter. In the diquark-based case, it provides the full structure of the two-body fragmentation mechanism and enables the derivation of compact analytic expressions for the SDCs. This is particularly useful for modeling ${\rm\Omega}_{3c}$ baryons, where handling full QCD three-body fragmentation becomes analytically and symbolically prohibitive.

Adopting Suzuki to compute the perturbative part of the FF is motivated by both analytical practicality and physical consistency. While rooted in perturbative QCD, the approach forgoes operator-based expansions in favor of a more “diagrammatic” strategy. This yields a simplified yet effective tool for evaluating SDCs, especially in the case of complex multi-quark final states such as triply heavy baryons.

First, the Suzuki picture offers an analytically accessible framework that enables the closure of phase-space integrals and yields manageable expressions for amplitudes. Spin and color algebra are simplified, and the use of the wave function at the origin effectively absorbs the bound-state dynamics without invoking full nonperturbative QCD machinery. This is especially valuable in contexts where rigorous QCD calculations requiring renormalization, operator regularization, or twist expansions are unfeasible or unavailable.

Second, the Suzuki model retains an explicit connection to the internal structure of the hadron. Rather than treating the baryon as a pointlike final state, it describes the fragmentation process using constituent-level vertices, including spin-dependent terms, propagators for internal degrees of freedom, and color couplings. This construction ensures that the hard subprocess remains sensitive to the underlying baryonic wave function, preserving the link between perturbative production and hadron-level observables.

Third, the model is particularly well suited for baryons like ${\rm\Omega}_{3c}$, where direct short-distance transitions in full perturbative QCD involve the algebraic disentangling of intertwined Dirac, Lorentz, and color structures. By rephrasing the process as a two-body fragmentation, where the diquark is treated as an effective, quasi-elementary object, the Suzuki model circumvents these challenges while preserving the overall dynamical content.

The analogy with NRQCD is particularly instructive. Just as NRQCD allows for a systematic separation of short-distance and long-distance contributions in quarkonium production, the Suzuki model supports a clean factorization between perturbative dynamics, described by the SDCs, and nonperturbative hadronization. Both frameworks rely on the use of effective composite degrees of freedom (diquark or heavy-quark nonrelativistic fields), which serve to simplify the treatment of fragmentation processes while maintaining fidelity to the underlying physics.

A crucial component in the construction of our initial-scale FFs is the modeling of the nonperturbative transition from the constituent diquark ${\cal D}$ to the physical baryon ${\rm\Omega}_{3c}$. Following Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2019), in this work we adopt the phenomenological Peterson-Schlatter-Schmitt-Zerwas (PSSZ) parametrization Peterson et al. (1983), which has been widely used in heavy-flavor fragmentation studies and is particularly suited for scenarios involving hadrons with multiple heavy quarks. We have

$$D_{\cal D}^{{\rm\Omega}_{3c}}(z)\equiv D_{\rm[np]}^{\rm[PSSZ]}(z)={\cal N}_{\cal P}\frac{z(1-z)^{2}}{\left[(1-z)^{2}+b_{\cal P}z\right]^{2}}\;,$$

(2)

with the normalization factor

$$\begin{split}{\cal N}_{\cal P}&\,=\,\left\{\frac{b_{\cal P}^{2}-6b_{\cal P}+4}{(4-b_{\cal P})\sqrt{4b_{\cal P}-b_{\cal P}^{2}}}\right.\\ &\times\,\left[\arctan\frac{b_{\cal P}}{\sqrt{4b_{\cal P}-b_{\cal P}^{2}}}+\arctan\frac{2-b_{\cal P}}{\sqrt{4b_{\cal P}-b_{\cal P}^{2}}}\right]\\ &+\,\left.\frac{1}{2}\ln b_{\cal P}+\frac{1}{4-b_{\cal P}}\right\}\;.\end{split}$$

(3)

being obtained by summing over all hadrons containing the charm quark Cacciari and Greco (1997) (see also Section IV of Ref. Peterson et al. (1983)). From the original derivation of the model, we note that the $b_{\cal P}$ cannot be calculated exactly, since it is intrinsically connected to nonperturbative effects.

Indeed, for a singly charmed hadron, say a $D$ meson with a leading Fock state $|c\bar{q}\rangle$, the value of $b_{\cal P}$ is related to the charm-quark mass by $b_{\cal P}\approx{\rm\Lambda}/m_{c}$, where ${\rm\Lambda}$ stands for a hadronic scale of the order of the mass of the constituent light antiquark, $m_{\bar{q}}$. In such a case, one finds that $\langle z\rangle=1-\sqrt{b_{\cal P}}$, which is in line with the prediction of heavy-quark FFs scaling linearly with the heavy-quark mass Suzuki (1977); Bjorken (1978); Kinoshita (1986); Cacciari et al. (1997). On the contrary, for a fully charmed baryon ${\rm\Omega}_{3c}$, the $b_{\cal P}$ is connected to the charm and anticharm masses, which are above the perturbative threshold. Following Refs. Adamov and Goldstein (1997); Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2019), we set $b_{\cal P}=[m_{c}/(m_{c}+m_{\bar{c}})]^{2}=1/4$.

This choice is supported by previous studies, where the Peterson form is explicitly used to model the nonperturbative stage of the diquark hadronization in fragmentation to triply heavy baryons Moosavi Nejad (2017); Delpasand and Moosavi Nejad (2019). Work in Ref. Moosavi Nejad (2017) clearly emphasize that the Peterson function effectively captures the essential nonperturbative dynamics of the fragmentation process and is compatible with the overall two-step fragmentation scheme defined through the Suzuki model. Similarly, Ref. Delpasand and Moosavi Nejad (2019) supports the use of the Peterson parametrization on the basis of its widespread adoption and its ability to encode realistic hadronization patterns for heavy baryons. Moreover, the Peterson form complements the output of the perturbative calculation by providing a universal and physically interpretable nonperturbative component that fits naturally into collinear factorization frameworks.

The Peterson model, by contrast, was originally formulated to describe the fragmentation of heavy quarks in which the initiating heavy quark carries most of the hadron momentum. This structure is naturally suited to hadrons containing heavy constituents, since the resulting hadron tends to follow the direction and kinematics of the fragmenting parton. In the case of triply heavy baryons, where all valence quarks are heavy, the same reasoning applies: the baryon is formed predominantly from the massive fragments of the initial parton, and thus inherits a large fraction of its momentum. In this respect, the Peterson function provides a realistic and physically consistent parametrization of the nonperturbative transition, capturing the expected peak at large z and the suppression at small z, in line with the kinematics of systems such as triply heavy baryons.

This behavior is particularly well suited for modeling the hadronization of a heavy diquark into a triply heavy baryon, where the recombination process is dominated by massive constituents and the baryon inherits a substantial fraction of the momentum of the initiating quark or gluon. The simple analytic structure of the model, governed by a single shape parameter $b_{\cal P}$, allows for the adjustment to match the internal mass hierarchy and spatial configuration of the final-state baryon, while ensuring smooth integration within collinear convolution schemes.

In this context, the Peterson form meets several important criteria: it is grounded in heavy-flavor dynamics, it offers analytical tractability, it aligns with the factorized structure emerging from the Suzuki picture, and it preserves consistency with the DGLAP evolution in the VFNS. For these reasons, we consider it a natural and robust choice for describing the nonperturbative input to the ${\rm\Omega}_{3c}$ initial-scale FF in our two-step formalism.

The compact expression for the FF of a parton $i$ generating a ${\rm\Omega}_{3Q}$ baryon at the initial scale $\mu_{F,0}$ reads

$$D_{i}^{{\rm\Omega}_{3c}}(z,\mu_{F,0})=\int_{z}^{1}\frac{{\rm d}\xi}{\xi}D_{i}^{\cal D}(\xi,\mu_{F,0})\,D_{\cal D}^{{\rm\Omega}_{3c}}\left(\frac{z}{\xi}\right)\;.$$

(4)

Equation (4) essentially matches the general structure of the heavy-hadron initial-scale FF as in Eq. (1), with two suitable adjustments. We notice that

	$\displaystyle D_{i}^{Q}(\xi,\mu_{F,0})\quad$	$\displaystyle\longrightarrow\quad D_{i}^{\cal D}(\xi,\mu_{F,0})\;,$		(5a)
	$\displaystyle D_{\rm[np]}^{{\cal H}_{Q}}\left(\frac{z}{\xi}\right)\quad$	$\displaystyle\longrightarrow\quad D_{\cal D}^{{\rm\Omega}_{3c}}\left(\frac{z}{\xi}\right)\;.$		(5b)

Here we have SDCs describing the perturbative splitting of a massless parton $i$ into a diquark state ${\cal D}$ (Eq. (5a)), followed by the long-distance, nonperturbative hadronization of the diquark into the observed triply heavy baryon ${\rm\Omega}_{3c}$ (Eq. (5b)). In other words, the diquark ${\cal D}$ replaces the $Q$ entering the lowest Fock state of the singly heavy hadron ${\cal H}_{Q}$. As mentioned, the nonperturbative $D_{\cal D}^{{\rm\Omega}_{3c}}$ will be modeled according to PSSZ as in Eq. (2), while gluon and charm SDCs will be presented and discussed in the next section.

The initial-scale FF for the $[c\to{\rm\Omega}_{3c}]$ transition in the diquark-like scenario is diagrammatically shown in the left panel of Fig. 1. In this representation, the double lines denote dicharm states, ${\cal D}\equiv|cc\rangle$, or their corresponding antidiquark partners, $\bar{{\cal D}}\equiv|\bar{c}\bar{c}\rangle$. The orange blob illustrates the nonperturbative hadronization component of the FF, labeled as $D_{\rm[np]}^{{\cal H}_{Q}}\equiv D_{\cal D}^{{\rm\Omega}_{3c}}$. The black bullet vertex represents the scalar diquark form factor coupling.¹¹1While not explicitly shown in the formulæ of this Section, our treatment mirrors that of Ref. Moosavi Nejad (2017), where the form factor is assumed to exhibit a simple pole in transverse-momentum space at 5 GeV.

We note that the diagram on the left panel of Fig. 1 corresponds to the $[c\to(c,{\cal D})+\bar{\cal D}]$ splitting. By swapping all charm quarks with anticharms, one obtains the $[\bar{c}\to(\bar{c},\bar{\cal D})+{\cal D}]$ process. A detailed analysis of potential asymmetries in the production of triply charmed $\rm\Omega$ baryons and antibaryons is left for future dedicated studies. In this work, we assume full symmetry in the formation mechanisms of ${\rm\Omega}_{3c}$ and $\bar{{\rm\Omega}}_{3c}$ states. Accordingly, our predictions are based on observables sensitive to the inclusive, charge-averaged emission of baryons and antibaryons. Assuming this symmetry holds, the $c$ and $\bar{c}$ fragmentation channels are treated on equal footing (for comparison with the light-hadron case, see Ref. Bertone et al. (2018)).

Using symJETHAD Celiberto (2021a, 2022a, 2023a, 2024b, 2024c) and FeynCalc Mertig et al. (1991); Shtabovenko et al. (2016, 2020), we symbolically calculated and obtained the explicit form of the LO and NLO initial-scale charm SDC, finding agreement with Refs. Moosavi Nejad and Delpasand (2017) and Moosavi Nejad (2017), respectively. As for the LO case, we write

$$\begin{split}D_{c}^{\cal D}(z,\mu_{F,0})\,&\equiv\alpha_{s}^{2}{\rm d}_{c}^{\rm[LO]}(z;\mu_{F,0})\\[2.84544pt] \,&+\,\alpha_{s}^{3}\,{\rm d}_{c}^{\rm[NLO]}(z;\mu_{F,0})\;.\end{split}$$

(6)

with $\alpha_{s}\equiv\alpha_{s}(5m_{c})$ and

$${\rm d}_{c}^{\rm[LO]}(z;\mu_{F,0})={\cal N}_{c}(z)\frac{{\cal S}_{c}^{\rm[LO]}(z;{\cal R}_{q_{T}/c})}{{\cal T}_{c}^{\rm[LO]}(z;{\cal R}_{q_{T}/c})}\;,$$

(7)

where

$${\cal R}_{q_{T}/c}^{2}\equiv\langle{\vec{q}}_{T}^{\;2}\rangle/m_{c}^{2}\;,$$

(8)

$${\cal N}_{c}(z)=\frac{\pi^{2}}{\sqrt{6}}\,f_{\cal B}^{2}C_{F}^{2}\,z^{3}(1-z)^{3}\;,$$

(9)

$$\begin{split}{\cal S}_{c}^{\rm[LO]}&(z;{\cal R}_{q_{T}/c})\,=\,4(41z^{2}-80z+96+256/z^{2})\\ \,&+\,{\cal R}_{q_{T}/c}^{2}(8z^{3}+5z^{2}-16z+16)\\ \,&+\,{\cal R}_{q_{T}/c}^{4}\,4z^{2}\;,\end{split}$$

(10)

and

$$\begin{split}{\cal T}_{c}^{\rm[LO]}&(z;{\cal R}_{q_{T}/c})\,=\,({\cal R}_{q_{T}/c}^{2}+(4-3z)^{2})^{2}\\ \,&\times\,{\cal R}_{q_{T}/c}^{2}\,z^{2}+z^{2}-16z+16)^{2}\;.\end{split}$$

(11)

The NLO correction to the charm SDC was originally derived in Ref. Moosavi Nejad and Sartipi Yarahmadi (2016) for singly heavy mesons, and then employed for triply-heavy baryons (see Ref. Moosavi Nejad (2017)). It reads

$${\rm d}_{c}^{\rm[NLO]}(z;\mu_{F,0})={\cal N}_{c}(z)\frac{{\cal S}_{c}^{\rm[NLO]}(z;{\cal R}_{q_{T}/c})}{{\cal T}_{c}^{\rm[NLO]}(z;{\cal R}_{q_{T}/c})}\;,$$

(12)

where

$$\begin{split}{\cal S}_{c}^{\rm[NLO]}&(z;{\cal R}_{q_{T}/c})\,=\,96\pi z\left[\right.{\cal R}_{q_{T}/c}^{14}\,(17-20z)z^{14}\\ \,&+\,{\cal R}_{q_{T}/c}^{12}\,(240z^{4}-1664z^{3}\\ \,&+\,3813z^{2}-4130z+1711)z^{12}\\ \,&+\,{\cal R}_{q_{T}/c}^{10}\,(4368z^{6}-39072z^{5}+140635z^{4}\\ \,&-\,269144z^{3}+300490z^{2}\\ \,&-\,189616z+52243)z^{10}\\ \,&+\,{\cal R}_{q_{T}/c}^{8}\,(20976z^{8}-208040z^{7}\\ \,&+\,870151z^{6}-2070634z^{5}+3255115z^{4}\\ \,&-\,3533740z^{3}+2600201z^{2}\\ \,&-\,1202330z+268205)z^{8}\\ \,&+\,{\cal R}_{q_{T}/c}^{6}\,4(1-z)^{2}(11256z^{8}-75035z^{7}\\ \,&-\,151779z^{6}+2535886z^{5}\\ \,&-\,8913460z^{4}+15864355z^{3}\\ \,&-\,15844619z^{2}+8395110z-1823490)z^{6}\\ \,&+\,{\cal R}_{q_{T}/c}^{4}\,4(1-z)^{4}(12324z^{8}-22830z^{7}\\ \,&-\,2045973z^{6}+14641206z^{5}-42277675z^{4}\\ \,&+\,65338590z^{3}-57387186z^{2}\\ \,&+\,26330832z-4933872)z^{4}\\ \,&+\,{\cal R}_{q_{T}/c}^{2}\,48(1-z)^{6}(563z^{8}+2272z^{7}\\ \,&-\,252314z^{6}+1752016z^{5}-5068605z^{4}\\ \,&+\,7770780z^{3}-6800544z^{2}\\ \,&+\,3505032z-849528)z^{2}\\ \,&+\,144(1-z)^{8}(41z^{8}+438z^{7}-36063z^{6}\\ \,&+\,258240z^{5}-833328z^{4}+1614600z^{3}\\ \,&-\,1793016z^{2}+676512z-198288)\left.\right]\;,\end{split}$$

(13)

and

$$\begin{split}{\cal T}_{c}^{\rm[NLO]}(z;{\cal R}_{q_{T}/c})\,&=\,[{\cal R}_{q_{T}/c}^{2}\,z^{2}+4(3-2z)^{2}]^{2}\\ \,&\times\,[{\cal R}_{q_{T}/c}^{2}\,z^{2}+(6-z)^{2}]^{2}\\ \,&\times\,[{\cal R}_{q_{T}/c}^{2}\,z^{2}+(1-z)^{2}]\\ \,&\times\,[{\cal R}_{q_{T}/c}^{2}\,z^{2}+36(1-z)^{2}]^{2}\\ \,&\times\,[{\cal R}_{q_{T}/c}^{2}\,z^{2}+z^{2}-35z+36]\;.\end{split}$$

(14)

Based on LO kinematics, we set the initial scale of charm fragmentation to $\mu_{F,0}\equiv 5m_{c}$ (see Fig. 1, left panel).

Our charm FF differs from the one originally introduced in Ref. Moosavi Nejad (2017) in two main aspects. First, in that study, the constant ${\cal N}_{c}(z)$ in Eq. (9) was not explicitly computed but was fixed by means of a normalization condition. Second, the choice of the $\langle{\vec{q}}_{T}^{\;2}\rangle$ parameter, directly connected with the use of the Suzuki picture, requires further refinement.

The Suzuki model effectively accounts for spin correlations and serves as a proxy for transverse-momentum-dependent (TMD) FFs. In the collinear limit, instead of integrating over the squared transverse momentum of the outgoing charm quark, one may replace it with its average value, $\langle{\vec{q}}_{T}^{\;2}\rangle$. This substitution renders $\langle{\vec{q}}_{T}^{\;2}\rangle$ a free parameter that must be fixed phenomenologically. Increasing $\langle{\vec{q}}_{T}^{\;2}\rangle$ shifts the FF peak to lower values of $z$ while reducing the overall normalization Gomshi Nobary (1994).

The initial-scale quark FF presented in Ref. Moosavi Nejad (2017) was obtained by adopting $\langle{\vec{q}}_{T}^{\;2}\rangle=1\,\text{GeV}^{2}$, regarded as an upper-limit estimate for the average squared transverse momentum. In this work, we introduce a more refined selection of $\langle{\vec{q}}_{T}^{\;2}\rangle$, based on a balanced and phenomenologically motivated approach consistent with the exploratory goals of our analysis. This improvement is based on a preliminary refinement introduced in our previous investigation of collinear fragmentation into fully charmed tetraquarks Celiberto et al. (2024a). There, we drew on phenomenological insights from fragmentation studies of various hadrons produced in proton collisions. In particular, heavy-quark FFs for both light hadrons Celiberto et al. (2016a, 2017); Bolognino et al. (2018a); Celiberto (2021a) and heavy-flavored ones Celiberto et al. (2021a, b); Celiberto and Fucilla (2022a); Celiberto (2022b) were observed to be typically probed at average longitudinal momentum fractions $\langle z\rangle>0.4$.

Based on a numerical analysis, we determined that setting $\langle{\vec{q}}_{T}^{\;2}\rangle_{T_{4c}}=70\,\text{GeV}^{2}$ for charm-to-$T_{4c}$ FFs results in an average $z$ greater than 0.4 and preserves a comparable order of magnitude between the quark and gluon channels. Extending this strategy to fully charmed $|ccc\bar{c}c\rangle$ pentaquarks, we adopted $\langle{\vec{q}}_{T}^{\;2}\rangle_{P_{5c}}=90\,\text{GeV}^{2}$ in Ref. Celiberto (2025d). Here, we fix $\langle{\vec{q}}_{T}^{\;2}\rangle_{{\rm\Omega}_{3c}}$ by correlating it with the peak position of the initial charm FFs. Following the same methodology used for $\langle{\vec{q}}_{T}^{\;2}\rangle_{T_{4c}}$ and $\langle{\vec{q}}_{T}^{\;2}\rangle_{P_{5c}}$, we performed a numerical scan and selected $\langle{\vec{q}}_{T}^{\;2}\rangle_{{\rm\Omega}_{3c}}=60\,\text{GeV}^{2}$.

A deeper justification supports our choice, going beyond heuristic arguments. Foundational studies on heavy-flavor fragmentation Suzuki (1977); Bjorken (1978); Kinoshita (1986); Peterson et al. (1983) have shown that heavy-quark FFs typically peak at large values of $z$, with binding effects scaling in proportion to the heavy-quark mass. To illustrate this behavior, consider the fragmentation process leading to a $D$ meson, whose lowest Fock state is $|c{\bar{q}}\rangle$, with total momentum $k$ and mass $M_{D}$.

In this context, the constituent heavy quark and the light antiquark are assumed to move with approximately the same velocity, $v\equiv v_{c}\simeq v_{q}$. Their momenta can thus be expressed as $k_{c}\equiv zk=m_{c}v$ for the charm quark and $k_{q}={\rm\Lambda}_{q}v$ for the light antiquark, where ${\rm\Lambda}_{q}$ is a hadronic mass scale of order ${\rm\Lambda_{\rm QCD}}$. Since the mass of the $D$ meson satisfies $M_{D}\approx m_{c}$, one finds $m_{c}v\approx k=k_{c}+k_{q}=zm_{c}v+{\rm\Lambda}_{q}v$. This implies the approximate relation $\langle z\rangle_{c}\approx 1-{\rm\Lambda}_{q}/m_{c}$, where the subscript ‘$c$’ indicates the $[c\to D]$ FF.

As pointed out in Refs. Celiberto et al. (2024a); Celiberto (2025d), this behavior does not generally apply to fully heavy-flavored states such as quarkonia, triply charmed baryons, tetracharms, or pentacharms. In the valence configuration of these systems, no soft scale is present, as their lowest Fock components contain only heavy quarks. For fully heavy rare baryons or exotics like $T_{4c}$ or $P_{5c}$ states, the complex dynamics among the three, four, or five heavy constituents hinders a simple kinematic determination of the peak position of the FF.

As a preliminary example, we present the $z$ dependence of the $[c\to{{\rm\Omega}_{3c}}]$ initial-scale FFs used in the OMG3Q1.0 set. To estimate uncertainties at the starting scale, we apply a simplified, diagonal DGLAP evolution at NLO, similar to that in Ref. Celiberto et al. (2024a), but using only the $P_{qq}$ kernel for the charm FF. The left panel of Fig. 2 shows the $z$-weighted charm FF, with $\mu_{F,0}=5m_{c}$ and scale variations in the range $4m_{c}<\mu_{F}<6m_{c}$. The $z$-weighted FF peaks in the $0.45<z<0.5$ range and drops to zero at the endpoints, as expected.

The small negative dip visible near $z\gtrsim 0.8$ may appear anomalous at first glance. However, dedicated numerical tests show that this feature originates from the choice of a relatively large value of $\langle{\vec{q}}_{T}^{\;2}\rangle$, higher than the one adopted in Ref. Moosavi Nejad (2017). This has no practical impact, as our tests confirm that FFs contribute dominantly to hadroproduction observables in the region $0.4\lesssim z\lesssim 0.6$ (see Section 3.3). For $z>0.8$, we have explicitly verified that no instability arises in high-energy cross-section predictions.

The initial-scale FF for the $[g\to{\rm\Omega}_{3c}]$ transition in the diquark-like scenario is diagrammatically shown in the right panel of Fig. 1. Although both the charm- and gluon-initiated channels are considered leading within their respective topologies, they differ in their perturbative order. As shown in Fig. 1, the gluon-induced process involves one additional QCD vertex compared to the charm-induced one, thus contributing at a higher order in the strong coupling expansion. Specifically, the $[c\to{\rm\Omega}_{3c}]$ channel is known both at LO, corresponding to the left diagram, and at NLO, as discussed in Section 3.1. In contrast, the $[g\to{\rm\Omega}_{3c}]$ channel is currently known only at $\mathcal{O}(\alpha_{s}^{3})$, which corresponds to the diagram shown in the right panel.

This difference comes from the fact that a gluon, which does not carry a charm quantum number, must generate at least three charm quarks to produce ${\rm\Omega}_{3c}$. Therefore, the first perturbative order at which such a process can occur is $\mathcal{O}(\alpha_{s}^{3})$, since the gluon must emit a $[c\bar{c}]$ pair to supply the free charm quark, generate a $|cc\rangle$ system to form the diquark (either directly or through an equivalent intermediate process) and include at least one additional QCD vertex to complete the baryon formation diagram. For consistency with this perturbative hierarchy, we denote the initial-scale gluon FF as the NLO function, while reserving the LO label for the lowest-order charm FF.

With symJETHAD Celiberto (2021a, 2022a, 2023a, 2024b, 2024c) and through the interface to FeynCalc Mertig et al. (1991); Shtabovenko et al. (2016, 2020), we symbolically computed the explicit form of the NLO initial-scale gluon SDC, originally provided in Ref. Delpasand and Moosavi Nejad (2019). One has

$$D_{g}^{\cal D}(z,\mu_{F,0})\,\equiv\,\alpha_{s}^{3}\,{\rm d}_{g}^{\rm[NLO]}(z;\mu_{F,0})\;.$$

(15)

with $\alpha_{s}\equiv\alpha_{s}(6m_{c})$ and

$${\rm d}_{g}^{\rm[NLO]}(z;\mu_{F,0})={\cal N}_{g}(z)\frac{{\cal S}_{g}^{\rm[NLO]}(z;{\cal R}_{q_{T}/c})}{{\cal T}_{g}^{\rm[NLO]}(z;{\cal R}_{q_{T}/c})}\;,$$

(16)

where

$${\cal N}_{g}(z)=\frac{2\pi^{3}}{3}\,f_{\cal B}^{2}C_{F}^{2}\,z^{3}(1-z)^{2}\;,$$

(17)

$$\begin{split}{\cal S}_{g}^{\rm[NLO]}(z;{\cal R}_{q_{T}/c})\,&=\,8(16z^{2}-32z+15)\\ \,&+\,{\cal R}_{q_{T}/c}^{2}\,2z^{2}(4z^{2}-20z+17)\\ \,&+\,{\cal R}_{q_{T}/c}^{4}\,z^{4}\end{split}$$

(18)

and

$$\begin{split}{\cal T}_{g}^{\rm[NLO]}&(z;{\cal R}_{q_{T}/c})\,=\,(4+z^{2}{\cal R}_{q_{T}/c}^{2})^{5}\;.\end{split}$$

(19)

As in Section 3.1 (see Eq. (8)), here we set ${\cal R}_{q_{T}/c}^{2}\equiv\langle{\vec{q}}_{T}^{\;2}\rangle/m_{c}^{2}$, with $\langle{\vec{q}}_{T}^{\;2}\rangle\equiv\langle{\vec{q}}_{T}^{\;2}\rangle_{{\rm\Omega}_{3c}}=60\,\text{GeV}^{2}$. In contrast to the charm case, the LO kinematics of gluon fragmentation set the initial scale at $\mu_{F,0}\equiv 6m_{c}$ (see Fig. 1, right diagram).

Similarly to the charm-induced case, we apply an analogous strategy to examine the $[g\to{{\rm\Omega}_{3c}}]$ fragmentation channel. As a preliminary illustration, we show the $z$ dependence of the gluon initial-scale FF adopted in the OMG3Q1.0 set. To evaluate uncertainties at the starting scale, we perform a simplified, diagonal NLO DGLAP evolution using only the $P_{gg}$ kernel, following the same methodology employed in Ref. Celiberto et al. (2024a) for gluon-induced tetraquark fragmentation.

The right panel of Fig. 2 displays the $z$-weighted gluon FF, evaluated at $\mu_{F,0}=6m_{c}$, and evolved within the range $5m_{c}<\mu_{F}<7m_{c}$. As expected, the function vanishes at $z\to 0$ and $z\to 1$, and exhibits a pronounced peak around $z\simeq 0.25$. Compared to charm FF, the gluon curve is wider and slightly shifted to lower $z$, reflecting a more uniform sharing of energy among final-state partons in gluon-initiated splittings.

We note that charm FF dominates by several orders of magnitude over its gluon counterpart. This suppression is not merely an artifact of the scalar diquark model, rather it reflects a structural feature of fragmentation into triply heavy baryons. While final states with an even number of heavy quarks, such as quarkonia or fully heavy tetraquarks, can naturally emerge from gluon splittings, the exclusive production of three heavy quarks in the valence state—with no accompanying light partons—poses a significant challenge for gluon fragmentation.

The scalar diquark framework, which presumes an $\mathcal{O}(\alpha_{s})$ suppression compared to the quark-induced channel, emphasizes this aspect without rigidly imposing it. Instead, it offers a phenomenologically motivated picture wherein the hierarchical suppression of gluon-induced baryon production is not imposed, but rather emerges as a natural outcome of the model.

To complete the construction of our OMG3Q1.0 collinear FFs for triply charmed ${\rm\Omega}$ baryons, we perform a proper DGLAP evolution of the initial-scale inputs presented before. A key distinction from light-hadron fragmentation is that in this case both the (anti-)charm and gluon channels are subject to evolution thresholds. This feature stems directly from the perturbative nature of partonic splittings into the lowest Fock state: gluon-initiated and charm-initiated processes correspond to the right and left panels of Fig. 1, respectively, and are governed by their associated SDCs.

As outlined earlier, kinematic constraints dictate that the minimum invariant mass required for quark fragmentation is $5m_{c}$, which we adopt as the threshold scale for the charm channel. In contrast, gluon-induced production requires a higher threshold $6m_{c}$, which we define as the starting scale for the (anti)charm FF. To handle the evolution while properly incorporating these thresholds, we adopt a dedicated strategy based on the recently developed HF-NRevo scheme Celiberto (2024a, 2025e, 2025f).

This framework is tailored to describe DGLAP evolution of heavy-hadron FFs starting from nonrelativistic inputs and relies on three core components: interpretation, evolution, and uncertainty quantification. The first component allows us to interpret the low-transverse-momentum production mechanism as a two-parton fragmentation process, as outlined in Section 2.1. This forms the basis for a consistent matching between FFNS and VFNS calculations. The third component provides a systematic approach to estimate the impact of missing higher-order uncertainties (MHOUs) by varying the evolution thresholds.

Although originally developed to bridge precision QCD predictions with a hadron-structure-oriented vision for quarkonium fragmentation, the HF-NRevo methodology has recently been extended to address exotic matter production, with promising results. Applications to scalar ($0^{++}$) and tensor ($2^{++}$) $T_{4c}$ states Celiberto et al. (2024a); Celiberto and Gatto (2025), as well as to rare baryons such as ${\rm\Omega}_{3c}$ (as explored in this work), have established HF-NRevo as a flexible framework for evolving FFs that involve both constituent heavy-quark and gluon initial-scale inputs. This dual-channel structure introduces the need for a dedicated treatment of evolution thresholds, tailored to each partonic species.

In this exploratory study on ${\rm\Omega}_{3c}$ fragmentation, we postpone the implementation of matching procedures and uncertainty quantification. Instead, we concentrate on developing a robust strategy for handling DGLAP evolution in the presence of quark and gluon thresholds.

According to HF-NRevo, the DGLAP evolution of heavy-hadron FFs is implemented in two successive stages. In the context of rare baryon production, we begin by evolving the channel with the lowest threshold, namely the charm-initiated one. The initial condition is provided by the (anti)charm FF at $\mu_{F,0}=5m_{c}$, as obtained in Section 3.1. This FF is then evolved up to $6m_{c}$, the gluon threshold, using only the $P_{qq}$ kernel. Since this evolution is both expanded in powers of $\alpha_{s}$ and decoupled from other partonic contributions, it can be performed analytically using the symJETHAD plugin Celiberto (2021a, 2022a, 2023a, 2024b, 2024c).

The second stage begins by combining the evolved (anti)charm FF at $6m_{c}$ with the gluon FF input introduced in Section 3.2. Starting from this common scale, we apply a full numerical DGLAP evolution to produce the NLO OMG3Q1.0 set, which we release in LHAPDF format. We define the starting point of this step as the evolution-ready scale, $Q_{0}$, corresponding to the highest threshold among active partonic species, $6m_{c}$ in this case. The $Q_{0}$ scale is the energy at which the numerical evolution is initialized. For this purpose, we use APFEL++ Bertone et al. (2014); Carrazza et al. (2015); Bertone (2018). In future developments, we also intend to interface with EKO Candido et al. (2022); Hekhorn and Magni (2023) to further expand our evolution technology.

One could argue that our treatment is incomplete because of the omission of light- and bottom-quark channels. To the best of our knowledge, no calculation currently exists for the collinear fragmentation of a nonconstituent quark into a triply heavy baryon. As a result, in our two-step evolution framework, light and bottom quarks are not assigned any initial-scale FFs and instead emerge dynamically through DGLAP evolution. However, drawing from analogies with NRQCD analyses of the color-singlet pseudoscalar Braaten and Yuan (1993); Braaten et al. (1993a); Artoisenet and Braaten (2015); Zhang et al. (2019); Zheng et al. (2021a, b) and vector Braaten and Yuan (1993); Braaten et al. (1993a); Zheng et al. (2019b) charmonia, these channels are expected to be strongly suppressed compared to those initiated by gluons and charm quarks.

The plot on the left of Fig. 3 displays the factorization scale dependence of the OMG3Q1.0 NLO functions that describe collinear VFNS fragmentation into ${\rm\Omega}_{3c}$ baryons. For comparison, the panel on the right displays the corresponding $\mu_{F}$ evolution of the TQ4Q1.1 NLO functions Celiberto and Gatto (2025), which describe the collinear VFNS fragmentation into the tensor $T_{4c}$ resonance ($2^{++}$), widely regarded as the leading candidate for the $X(6900)$ tetraquark Aaij et al. (2020). For conciseness, both plots are evaluated at a representative value of the momentum fraction, $z=0.5\simeq\langle z\rangle$, which approximates the typical region where FFs provide dominant contributions to semi-hard final states Celiberto et al. (2016a, 2017); Bolognino et al. (2018a); Celiberto (2021a); Celiberto et al. (2021a, b); Celiberto and Fucilla (2022a); Celiberto (2022b, 2023b).

A direct comparison of the results in Fig. 3 reveals that the (anti)charm-to-${\rm\Omega}_{3c}$ fragmentation channel emerges as the dominant contribution. It significantly outweighs the light-parton (not shown due to their negligible size) and bottom-quark channels across the entire $\mu_{F}$ range considered. Moreover, it exceeds the gluon contribution by approximately two orders of magnitude. This behavior contrasts with the case of the $T_{4c}(2^{++})$ state, where the gluon FF dominates over the charm FF, with ratios ranging from a factor of 5 up to an order of magnitude. The prevailing role of the charm-initiated FF in ${\rm\Omega}_{3c}$ production directly reflects the initial-scale hierarchy between the two channels, as discussed earlier (see Fig. 2). While the gluon FF increases with $\mu_{F}$ due to its feeding from timelike splittings, the charm FF exhibits a smoother trend and receives no comparable enhancement. This explains why the initial five-order hierarchy between the charm and gluon FFs is reduced to roughly two orders of magnitude after DGLAP evolution.

Despite the charm FF being substantially larger in magnitude than its gluon counterpart, the latter still plays a key role in predicting ${\rm\Omega}_{3c}$ (semi-)inclusive production rates at hadron colliders. This relevance stems from the overwhelming dominance of the gluon parton distribution function (PDF) over quark PDFs, which enhances the contribution of both the $[gg\to gg]$ partonic subprocess and the gluon-initiated fragmentation. As a result, the simultaneous inclusion of initial-scale inputs for both the charm and gluon channels in our OMG3Q1.0 determinations strengthens the robustness and reliability of the proposed methodology.

Finally, we observe that the gluon-to-${\rm\Omega}_{3c}$ FF exhibits a mild increase with $\mu_{F}$, while its gluon-to-$T{4c}$ counterpart displays a slow decrease. Both behaviors are smoothly varying with respect to the factorization scale, a property that carries significant phenomenological implications. Indeed, it has recently been shown that gluon FFs characterized by a smooth $\mu_{F}$-dependence act as effective “stabilizers” in high-energy resummed observables sensitive to the semi-inclusive production of singly Celiberto et al. (2021a, b) and multiply Celiberto and Fucilla (2022a); Celiberto (2022b); Celiberto and Papa (2024) heavy-flavored hadrons.

This remarkable feature has been referred to as the natural stability Celiberto (2023c) of high-energy resummation (see Section 4). The natural stability arising from the heavy flavor fragmentation will constitute the cornerstone of our forthcoming phenomenological analysis (see Section 5).

The production of heavy-flavored hadrons offers a valuable perspective on high-energy QCD, where energy logarithms become large enough to challenge perturbative expansions. The Balitsky-Fadin-Kuraev-Lipatov (BFKL) approach Fadin et al. (1975); Kuraev et al. (1977); Balitsky and Lipatov (1978) addresses this by systematically resumming energy logarithms to all orders, covering both leading logarithmic (LL) terms $\alpha_{s}^{n}\ln(s)^{n}$ and next-to-leading logarithmic (NLL) corrections $\alpha_{s}^{n+1}\ln(s)^{n}$.

In the BFKL framework, cross sections are computed as transverse-momentum convolutions involving a universal NLO Green’s function Fadin and Lipatov (1998); Ciafaloni and Camici (1998) and process-dependent, singly off-shell emission functions, also known as forward impact factors. These functions incorporate collinear elements such as PDFs and FFs, leading to a hybrid-factorization formalism that merges high-energy and collinear QCD dynamics.

Over the years, the BFKL resummation has been applied to a variety of processes, including Mueller-Navelet jets Mueller and Navelet (1987); Ducloué et al. (2013); Colferai and Niccoli (2015); Celiberto et al. (2015a, b, 2016b); Celiberto (2017); Caporale et al. (2018); de León et al. (2021); Celiberto and Papa (2022); Baldenegro et al. (2024), di-hadron Celiberto et al. (2016a, 2017); Celiberto (2017); Celiberto et al. (2020); Celiberto (2022a) and hadron-jet systems Bolognino et al. (2018a, 2019a, 2019b); Celiberto (2021a); Celiberto et al. (2020); Mohammed (2022); Celiberto (2023b), multijets Caporale et al. (2016, 2017a); Celiberto (2016); Caporale et al. (2017b); Celiberto (2017), forward-Higgs Hentschinski et al. (2021); Celiberto et al. (2022a, 2021c); Mohammed (2022); Celiberto and Papa (2023); Celiberto et al. (2023, 2024b, 2024c, 2024d, 2022b, 2024e), Drell-Yan Celiberto et al. (2018a); Golec-Biernat et al. (2018), and heavy-flavored emissions Celiberto et al. (2018b); Boussarie et al. (2018); Bolognino et al. (2019c, d); Celiberto et al. (2021a, b); Celiberto and Fucilla (2022a); Celiberto (2023a, c); Bolognino et al. (2023); Celiberto (2022b); Celiberto and Fucilla (2022b); Celiberto (2024d). Single forward production rates have provided insights into small-$x$ gluon dynamics via the UGD, studied at HERA Anikin et al. (2011); Besse et al. (2013); Bolognino et al. (2018b, c, 2019e, 2020); Celiberto (2019); Bolognino (2021); Łuszczak et al. (2022), the EIC Bolognino et al. (2021a, 2022a); Bolognino (2021); Bolognino et al. (2022b, c). This has subsequently enabled the development of resummed PDFs Ball et al. (2018); Abdolmaleki et al. (2018); Bonvini and Giuli (2019); Silvetti and Bonvini (2023); Silvetti (2023); Rinaudo (2024) and improved low-$x$ TMDs Bacchetta et al. (2020, 2024a); Celiberto (2021b); Bacchetta et al. (2022a, b, c, d, e, f); Celiberto (2022c); Bacchetta et al. (2024b).

In the context of heavy-hadron emissions, processes like ${\rm\Lambda}_{c}$ Celiberto et al. (2021a) and $b$-hadron production Celiberto et al. (2021b) have revealed mechanisms to address challenges in natural-scale descriptions of semi-hard processes. Unlike light-particle emissions, which suffer from large NLL corrections and nonresummed threshold effects Ducloué et al. (2014); Caporale et al. (2014); Bolognino et al. (2018a); Celiberto (2021a), heavy-flavored hadrons exhibit a natural stabilization trend Celiberto (2023c), driven by collinear VFNS fragmentation.

These findings have inspired further studies, including the development of VFNS DGLAP-evolving FFs based on NRQCD inputs Braaten et al. (1993a); Zheng et al. (2019b); Braaten and Yuan (1993); Chang and Chen (1992); Braaten et al. (1993b); Ma (1994); Zheng et al. (2019a, 2022); Feng et al. (2022b); Feng and Jia (2023), extending from vector quarkonia Celiberto and Fucilla (2022a); Celiberto (2023a) to $B_{c}(^{1}S_{0})$ and $B_{c}(^{3}S_{1})$ mesons Celiberto (2022b, 2024d). Natural stability has also opened a portal to exotic matter, serving as a phenomenological playground to study the fragmentation of the leading power doubly Celiberto and Papa (2024); Celiberto and Gatto (2025) or fully heavy tetraquarks Celiberto et al. (2024a); Celiberto and Gatto (2025); Celiberto (2025a) and pentacharms Celiberto (2025d).

As an application to hadron-collider phenomenology, we examine the following semi-inclusive reaction (see Fig. 4)

$${\rm p}(P_{a})\,{\rm p}(P_{b})\to{{\rm\Omega}_{3c}}(q_{1},y_{1},\phi_{1})+{\cal X}+{\rm jet}(q_{2},y_{2},\phi_{2})\;,$$

(20)

where a triply-charmed baryon ${\rm\Omega}_{3c}$ is detected with four-momentum $q_{1}$, rapidity $y_{1}$, and azimuthal angle $\phi_{1}$. In addition, a light jet is tagged with four-momentum $q_{2}$, rapidity $y_{2}$, and azimuthal angle $\phi_{2}$. Both objects feature high transverse momenta, such that $|\vec{q}_{1,2}|\gg{\rm\Lambda}_{\rm QCD}$, and a large rapidity distance, $\Delta Y\equiv y_{1}-y_{2}$. An undetected final-state gluon system, ${\cal X}$, is produced inclusively. We express the final-state transverse momenta in the Sudakov-vector basis, defined by the parent protons’ momenta, $P_{a,b}$, where $P_{a,b}^{2}=0$ and $({P_{a}}\cdot{P_{b}})=s/2$. This leads to the decomposition

$$q_{1,2}=x_{1,2}\,P_{a,b}-\frac{q_{1,2\perp}^{2}}{x_{1,2}s}\,P_{b,a}+q_{{1,2\perp}}\;,$$

(21)

with $q_{1,2\perp}^{2}\equiv-\vec{q}_{1,2}^{\,2}$. In the center-of-mass frame, the following relations hold between the final-state longitudinal momentum fractions, rapidities, and transverse momenta:

$$x_{1,2}=\frac{|\vec{q}_{1,2}|}{\sqrt{s}}e^{\pm y_{1,2}}\;,\qquad{\rm d}y_{1,2}=\pm\frac{dx_{1,2}}{x_{1,2}}\;,$$

(22)

and therefore

$$\qquad\Delta Y\equiv y_{1}-y_{2}=\ln\frac{x_{1}x_{2}s}{|\vec{q}_{1}||\vec{q}_{2}|}\;.$$

(23)

In a purely collinear factorization approach, the LO differential cross section for these reactions can be written as a one-dimensional convolution involving the on-shell hard factor, proton PDFs, and ${\rm\Omega}_{3c}$ FFs:

	$\displaystyle\frac{{\rm d}\sigma^{\rm LO}_{\rm[collinear]}}{{\rm d}x_{1}{\rm d}x_{2}{\rm d}^{2}\vec{q}_{1}{\rm d}^{2}\vec{q}_{2}}=\sum_{m,n=q,{\bar{q}},g}\int_{0}^{1}{\rm d}x_{a}\!\!\int_{0}^{1}{\rm d}x_{b}\ f_{m}\left(x_{a}\right)f_{n}\left(x_{b}\right)$
	$\displaystyle\quad\times\,\int_{x_{1}}^{1}\frac{{\rm d}\zeta}{\zeta}\,D^{{\rm\Omega}_{3c}}_{m}\left(\frac{x_{1}}{\zeta}\right)\frac{{\rm d}{\hat{\sigma}}_{m,n}\left(\hat{s}\right)}{{\rm d}x_{1}{\rm d}x_{2}{\rm d}^{2}\vec{q}_{1}{\rm d}^{2}\vec{q}_{2}}\;,$		(24)

where the indices $m,n$ sum over all partons except the (anti)top quark, which does not participate in hadronization. For simplicity, the explicit dependence on the factorization scale $\mu_{F}$ is omitted in Eq. (4.2). The functions $f_{m,n}\left(x_{a,b},\mu_{F}\right)$ represent the proton collinear PDFs, while the $D^{{\rm\Omega}_{3c}}_{m}\left(x_{1}/\zeta,\mu_{F}\right)$ functions are the ${\rm\Omega}_{3c}$ collinear FFs. The variables $x_{a,b}$ correspond to the longitudinal momentum fractions of the parent partons, and $\zeta$ is the momentum fraction of the outgoing parton fragmenting into the ${\rm\Omega}_{3c}$ particle. Lastly, ${\rm d}\hat{\sigma}_{m,n}\left(\hat{s}\right)$ is the partonic hard factor, where $\hat{s}\equiv x_{a}x_{b}s$ is the partonic center-of-mass energy squared.

Conversely, the differential cross section within our hybrid high-energy and collinear factorization approach is expressed as a transverse-momentum convolution involving the BFKL Green’s function and the two singly off-shell emission functions. We rewrite the cross section as a Fourier series of the azimuthal-angle coefficients, ${\cal C}_{k\,\geq\,0}$. We write

$$\frac{(2\pi)^{2}\,{\rm d}\sigma}{{\rm d}y_{1}{\rm d}y_{2}{\rm d}|\vec{q}_{1}|{\rm d}|\vec{q}_{2}|{\rm d}\phi_{1}{\rm d}\phi_{2}}\!=\!\left[{\cal C}_{0}+2\!\sum_{k=1}^{\infty}\!\cos(k\phi)\,{\cal C}_{k}\right]\,,$$

(25)

where $\phi\equiv\phi_{1}-\phi_{2}-\pi$.

In the $\overline{\rm MS}$ renormalization scheme and using BFKL, we obtain (see Ref. Caporale et al. (2013) for more details)

$${\cal C}_{k}^{\rm NLL/NLO^{+}}=\frac{e^{\Delta Y}}{s}\int_{-\infty}^{+\infty}\!\!\!{\rm d}\nu\,e^{{\Delta Y}\bar{\alpha}_{s}(\mu_{R})\chi^{\rm NLO}(k,\nu)}$$

(26)

$$\times\,\alpha_{s}^{2}(\mu_{R})\biggl{\{}{\cal E}_{{\rm\Omega}_{3c}}^{\rm NLO}(k,\nu,|\vec{q}_{1}|,x_{1})[{\cal E}_{J}^{\rm NLO}(k,\nu,|\vec{q}_{2}|,x_{2})]^{*}$$

$$+\left.\alpha_{s}^{2}(\mu_{R})\Delta Y\frac{\beta_{0}}{4\pi}\,\chi(k,\nu)\left[\ln\left(|\vec{q}_{1}||\vec{q}_{2}|\right)+\frac{i}{2}\,\frac{{\rm d}}{{\rm d}\nu}\ln\frac{{\cal E}_{{\rm\Omega}_{3c}}}{{\cal E}_{J}^{*}}\right]\right\}.$$

In this equation, $\bar{\alpha}_{s}(\mu_{R})\equiv\alpha_{s}(\mu_{R})N_{c}/\pi$ is the QCD running coupling, with $N_{c}=3$ representing the number of colors, and $\beta_{0}=11N_{c}/3-2n_{f}/3$ being the first coefficient of the QCD $\beta$-function. We choose a two-loop running coupling setup with $\alpha_{s}\left(m_{Z}\right)=0.118$, and the dynamic number of flavors, $n_{f}$.

The high-energy kernel in the exponent of Eq. (26) captures the resummation of energy logarithms at NLL accuracy:

$\displaystyle\chi^{\rm NLO}(k,\nu)=\chi(k,\nu)+\bar{\alpha}_{s}\hat{\chi}(k,\nu)\;,$

(27)

where $\chi(k,\nu)$ is the eigenvalue of the LO kernel:

$\displaystyle\chi\left(k,\nu\right)=-2\gamma_{\rm E}-2\,{\rm Re}\left\{\psi\left(\frac{1+k}{2}+i\nu\right)\right\}\,,$

(28)

with $\gamma_{\rm E}$ being the Euler–Mascheroni constant and $\psi(z)\equiv\Gamma^{\prime}(z)/\Gamma(z)$ the logarithmic derivative of the Gamma function. The function $\hat{\chi}(k,\nu)$ in Eq. (27) represents the NLO correction to the high-energy kernel:

		$\displaystyle\hat{\chi}$	$\displaystyle\left(k,\nu\right)=\bar{\chi}(k,\nu)+\frac{\beta_{0}}{8N_{c}}\chi(k,\nu)$
		$\displaystyle\times$	$\displaystyle\left\{-\chi(k,\nu)+10/3+2\ln\left[\left(\mu_{R}^{2}/(|\vec{q}_{1}||\vec{q}_{2}|\right)\right]\right\}\;,$

with the characteristic $\bar{\chi}(k,\nu)$ function calculated in Ref. Kotikov and Lipatov (2000).

The two quantities

$\displaystyle{\cal E}_{{\rm\Omega}_{3c},J}^{\rm NLO}(k,\nu,|\vec{q}_{1,2}|,x_{1,2})={\cal E}_{{\rm\Omega}_{3c},J}+\alpha_{s}(\mu_{R})\,\hat{\cal E}_{{\rm\Omega}_{3c},J}$

(30)

represent the NLO emission functions for the rare baryon (${\rm\Omega}_{3c}$) and the light jet ($J$), respectively. These objects are calculated in Mellin space and projected onto the eigenfunctions of the LO kernel. For the ${\rm\Omega}_{3c}$ emission function, we adopt the NLO calculation from Ref. Ivanov and Papa (2012a). Although originally tailored for light-flavored hadron production, this computation aligns with our VFNS approach for heavy baryons, provided that transverse momenta exceed the heavy-quark thresholds relevant for DGLAP evolution.

At the LO level, the ${\rm\Omega}_{3c}$ emission function reads

$$\begin{split}{\cal E}_{{\rm\Omega}_{3c}}(k,\nu,|\vec{q}_{1}|,x_{1})=\delta_{c}\,|\vec{q}_{1}|^{2i\nu-1}\int_{x}^{1}\frac{{\rm d}\zeta}{\zeta}\;\hat{x}^{1-2i\nu}\\ \times\,\Big{[}\rho_{c}f_{g}(\zeta)D_{g}^{{\rm\Omega}_{3c}}\left(\hat{x}\right)+\sum_{m=q,\bar{q}}f_{m}(\zeta)D_{m}^{{\rm\Omega}_{3c}}\left(\hat{x}\right)\Big{]}\;,\end{split}$$

(31)

where $\hat{x}=x/\zeta$, $\delta_{c}=2\sqrt{C_{F}/C_{A}}$, and $\rho_{c}=C_{A}/C_{F}$. Here, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ and $C_{A}=N_{c}$ are the Casimir factors associated with gluon emissions from quarks and gluons, respectively. The complete NLO formula for ${\cal E}_{{\rm\Omega}_{3c}}^{\rm NLO}$ is available in Ref. Ivanov and Papa (2012a).

The LO jet emission function is expressed as

$$\!\!c_{J}(k,\nu,|\vec{q}_{2}|,x)=\delta_{c}|\vec{q}_{2}|^{2i\nu-1}\!\Big{[}\rho_{c}f_{g}(x)+\!\!\!\sum_{n=q,\bar{q}}f_{n}(x)\Big{]}\,.$$

(32)

The corresponding NLO correction can be obtained by combining Eq. (36) of Ref. Caporale et al. (2013) with Eqs. (4.19) and (4.20) of Ref. Colferai and Niccoli (2015). These calculations are grounded in results from Refs. Ivanov and Papa (2012a, b), optimized for numerical analyses. They employ a small-cone jet selection function Furman (1982); Aversa et al. (1989) based on a cone-type algorithm Colferai and Niccoli (2015).

Equations (26), (31), and (32) clearly demonstrate the implementation of our hybrid high-energy and collinear factorization scheme. In this approach, the cross section is expressed through the BFKL formalism, with the Green’s function and emission functions as key components. The Green’s function handles the resummation of large logarithmic contributions in the high-energy limit, while the emission functions encapsulate the PDFs and FFs, effectively bridging collinear factorization with high-energy dynamics.

The ‘$+$’ superscript in the ${\cal C}_{k}^{\rm NLL/NLO^{+}}$ label signifies that the expression of azimuthal coefficients in Eq. (26) incorporates corrections beyond the NLL accuracy. These enhancements stem from two key sources: the exponentiated NLO corrections to the high-energy kernel and the interplay between the NLO corrections to the impact factors. As a result, the azimuthal coefficients provide a more refined representation, capturing intricate effects essential for precise predictions in processes where both high-energy and collinear logarithms are influential.

By discarding all NLO contributions in Eq. (26), one retrieves the pure LL limit of the angular coefficients, thus having

$$\begin{split}&{\cal C}_{k}^{\rm LL/LO}=\frac{e^{\Delta Y}}{s}\int_{-\infty}^{+\infty}{\rm d}\nu\,e^{{\Delta Y}\bar{\alpha}_{s}(\mu_{R})\chi(k,\nu)}\\[5.12128pt] &\hskip 4.26773pt\times\,\alpha_{s}^{2}(\mu_{R})\,{\cal E}_{{\rm\Omega}_{3c}}(k,\nu,|\vec{q}_{1}|,x_{1})[{\cal E}_{J}(k,\nu,|\vec{q}_{2}|,x_{2})]^{*}\;.\end{split}$$

(33)