Hunting BFKL in semi-hard reactions at the LHC

The agreement between calculations inspired by the resummation of energy logarithms, known as BFKL approach, and experimental data in the semi-hard sector of QCD has become manifest after a wealthy series of phenomenological analyses. However, the contingency that the same data could be concurrently portrayed at the hand of fixed-order, DGLAP-based calculations, has been pointed out recently, but not yet punctually addressed. Taking advantage of the richness of configurations gained by combining the acceptances of CMS and CASTOR detectors, we give results in the full next-to-leading logarithmic approximation of cross sections, azimuthal correlations and azimuthal distributions for three distinct semi-hard processes, each of them featuring a peculiar final-state exclusiveness. Then, making use of disjoint intervals for the transverse momenta of the emitted objects, i.e. $\kappa$-windows, we clearly highlight how high-energy resummed and fixed-order driven predictions for semi-hard sensitive observables can be decisively discriminated in the kinematic ranges typical of current and forthcoming analyses at the LHC. The scale-optimization issue is also introduced and explored, while the uncertainty coming from the use of different PDF and FF set is deservedly handled. Finally, a brief overview of JETHAD, a numerical tool recently developed, suited for the computation of inclusive semi-hard reactions is provided.

The search for evidence of New Physics is in the viewfinder of current and forthcoming analyses at the Large Hadron Collider (LHC). This is the best time to shore up our knowledge of strong interactions though, the high luminosity and the record energies reachable widening the horizons of kinematic sectors uninvestigated so far. A broad class of processes, called diffractive semi-hard reactions process and Λ QCD the QCD scale), is stringently preserved, gives us a faultless chance to test perturbative QCD (pQCD) in new and quite original ways 1 . In the kinematical regime, known as Regge limit, where s is much larger than the Mandelstam variable t, fixedorder calculations in pQCD based on collinear factorization fail to catch the effect of large-energy logarithmic contributions, entering the perturbative series with a power increasing with the order and thus balancing the slightness of the strong coupling, α s . The Balitsky-Fadin-Kuraev-Lipatov (BFKL) approach [2-5] serves as the most adequate tool to perform the all-order resummation of these large-energy logarithms both in the leading approximation (LLA), namely of terms proportional to (α s ln s) n and the next-to-leading approximation (NLA), namely of terms proportional to α s (α s ln s) n . In the BFKL formalism, the imaginary part the amplitude of a hadronic process (thence its cross section, in the case of inclusive final states, via the optical theorem [6]) is elegantly expressed via a suitable convolution, pictorically depicted in Fig. 1, between two impact factors, portraying the transition from each parent particle to the final-state object(s) belonging to its fragmentation region, and a process-independent Green's function. The BFKL Green's function depends on energy and its evolution is regulated by an integral equation, whose kernel is known at the next-to-leading order (NLO) both for forward scattering [7,8], i.e. for t = 0 and color singlet in the t-channel, and for any fixed (not growing with energy) t and any possible two-gluon color configuration in the t-channel [9][10][11][12]. Conversely, impact factors strictly depend on the final state and on the hard scale, but not on the energy. Thus, the number of processes which can be investigated via the BFKL resummation is shortened by the narrow list of available impact factors, just few of them being calculated at the NLO: 1) colliding-parton (quarks and gluons) impact factors [13][14][15][16][17], which represent the landmark for the calculation of the 2) forward-jet impact factor [18][19][20][21][22] and of the 3) forward charged-light hadron one [23], 4) the impact factor specific of the (γ * → LVM) transition [24] at leading twist, where LVM stands for light vector meson, and 5) and the one detailing the (γ * → γ * ) subprocess [25][26][27][28][29][30][31].
With the aim of deepening our understanding of the high-energy regime, a notable variety of semi-hard reactions (see Ref. [32] for applications) has been proposed so far.
The theoretical description of the (γ * -γ * ) total cross section relies on the (γ * → γ * ) impact factor. Although being recognized as a privileged channel for the manifestation of the BFKL dynamics, comparisons between the large number of BFKL predictions [67][68][69][70][71][72][73] with the only available data from the LEP2 Cern experiment were ineffective due to the relatively low center-of-mass energy and the limited accuracy of the detector. Quite recently [74,75], the amplitude for the light-by-light elastic scattering via a single-quark loop was calculated in the so-called double-logarithmic approximation (DLA). This approach 2 accounts for the resummation of those contributions proportional to (α s ln s ln s) n that, at variance with gluon interactions in the BFKL ladder, appear in the quark-antiquark exchange [86][87][88][89], giving rise to a different asymptotic behavior in the power of s with respect to the single-logarithmic trend of the BFKL Pomeron (for more details, see Refs. [90,91]).
The exclusive leptoproduction of two LVMs was investigated by encoding the (γ * → LVM) impact factors in the definition of the BFKL amplitude. In particular, first studies at Born level were done in Ref. [92] with longitudinally polarized virtual photons, and in Ref. [93] by accounting also for the photon transverse polarization. Then, a LLA treatment supplemented by improvements of the LO BFKL kernel was proposed in Ref. [94], while a full NLA analysis was performed in Refs. [95,96].
LO impact factors for the production of forward heavy-quark pairs [97,98] represent the key ingredients for phenomenological analyses with partial NLA resummation effects, by including next-to-leading logarithmic corrections to the BFKL Green's function, of the heavy-quark di-jet photo- [99,100] and hadroproduction [98,100], tracing the path towards a prospective NLA BFKL treatment of heavy-flavored-meson production channels. With the same accuracy, J/Ψ-jet [101], Higgs-jet [102] and heavy-light di-jet [103] correlations have been proposed, while a full LLA description, together with consistency-inspired collinear improvements [104] to the BFKL kernel, of the forward Drell-Yan dilepton production plus a backward jet has been recently conducted [105].
All the possibilities presented above fall into a distinctive family of reactions, where at least two final-state particles are always emitted with large mutual separation in rapidity, together (in the inclusive channel) with a secondary, undetected gluon system. Another interesting subclass of semi-hard processes is represented by those final states featuring the tag of a single forward object in lepton-proton or proton-proton scatterings.
2 As a meaningful remark, we note that sub-eikonal corrections, inevitably neglected by BFKL, are instead genuinely included in DLA-based calculations. Thus, DLA resummation comes out as a powerful tool to investigate polarization effects in the high-energy/small-x regime. In this context, the Bartels-Ermolaev-Ryskin (BER) approach permits us to calculate flavor non-singlet [76] and singlet contributions [77] to the g 1 helicity structure function in the Deep Inelastic Scattering (DIS) at small-x. Another way to resum double-logarithmic powers in the definition of small-x polarized parton densities is represented by the Kovchegov-Pitonyak-Sievert (KPS) formalism [78][79][80][81][82], whose evolution equations are built in terms of (polarized) Wilson lines [83] and account for saturation effects. Originally aimed at the study of quark and gluon helicity distributions, the KPS scheme was then extended to valence-quark transversity [84] and orbital angular momentum [85] in the small-x limit. reduction in rapidity with respect to jets. It gives the possibility, however, to probe and constrain not only PDFs, but also fragmentation functions (FFs) entering the expression of the hadronic cross section. The limiting acceptances imposed by the hadron identification can be compensated by considering the concurrent emission of a forward (backward) hadron and a backward (forward) jet, as recently proposed in a first work on hadron-jet correlations [63]. Thus, in our study we take advantage of the window of opportunites offered by all the possible final-state combinations given in Fig. 2, shaped on CMS and CASTOR acceptances. Imposing disjoint intervals for final-state emissions in the transverse-momentum plane (κ-windows) to shed light on the transition region between the two resummations, we provide indisputable evidence that the onset of BFKL dynamics can be disentangled from DGLAP-sensitive contaminations at current energies and kinematic configurations of prospective experimental studies at the LHC.
At the same time, we address the issue related to the "optimal" value for the renormalization scale (together with a supplementary analysis on single hadron-species detection) and gauge the effect of the uncertainty coming from the use of different PDF and FF parametrizations.
Finally, we present the main features of JETHAD, an object-based interface we have been developing, pursuing the goal to provide the scientific community with a common reference software for the phenomenological study of inclusive semi-hard processes.
The paper reads as follows: in Section 2 a general framework for the three considered reactions is set; in Section 3 results for observables of interest are presented and discussed; Section 4 brings our conclusions.
We can present a general formula for the processes under consideration (Eq. (3)), making use of collinear factorization: where the (α, β) indices run over the parton types (quarks q = u, d, s, c, b; antiquarksq =ū,d,s,c,b; or gluon g), f α,β (x, µ F 1,2 ) are the initial proton PDFs; x 1,2 stand for the longitudinal fractions of the partons involved in the hard subprocess, whereas µ F 1 (µ F 2 ) is the factorization scale characteristic of the fragmentation region of the upper (lower) parent proton in panels of Fig. 2; dσ α,β (ŝ, µ F 1,2 ) denotes the partonic cross section and s ≡ x 1 x 2 s is the squared center-of-mass energy of the partonic collision.

Cross section in the NLA BFKL
We can decompose the cross section (see Ref. [36] for further details) into a suitable Fourier expansion of the azimuthal coefficients, C n , getting so: where ϕ ≡ ϕ 1 − ϕ 2 − π, with ϕ 1,2 the azimuthal angles of the tagged objects, O 1,2 , and the Jacobian for the change of variables has been taken into account. The zeroth coefficient, C n=0 , gives the ϕ-averaged contribution to the cross section, while the other ones, C n =0 , are simply called azimuthal-correlation coefficients. It is known that various expressions for C n , equivalent in the NLA BFKL approximation, exist (for an extensive study of them and of their peculiarities, see Ref. [41]). For the purposes of the analysis proposed in this work, it is sufficient to consider just one representation, the so-called exponentiated one, where the following global formula for the C n ≡ C NLA n coefficients holds with NLA BFKL accuracy: with n is an integer larger than one andᾱ s (µ R ) ≡ α s (µ R )N c /π, where N c is the color number and β 0 the first coefficient of the QCD β-function, responsible for running-coupling effects (Eq. (19)). Details on the BFKL kernel and on the impact factors can be found in Sections 2.1.1 and 2.1.2, respectively. It is worth to remark that the formula given in Eq. (7) is obtained in the so-called (n, ν)-representation, i.e. taking the projection onto the eigenfunctions of the LO BFKL kernel (for more details, see Refs. [95,96]). For the sake of completeness, we give the expression for the azimuthal coefficients in the LLA BFKL approximation: All these formulae for the azimuthal coefficients are obtained in the MS renormalization scheme. Corresponding expressions in the MOM scheme and with scale-optimization setting will be given in Section 3.4.1.

LO and NLO BFKL kernel
The expression in Eq. (7) for the LO BFKL characteristic function or LO BFKL kernel is with ψ(z) ≡ Γ (z)/Γ(z) the logarithmic derivative of the Gamma function, whileχ(n, ν), calculated in Ref. [166] (see also Ref. [167]), stands for the NLO correction to the BFKL kernel, handily represented in the form: where n f is the flavor number, and

LO impact factors in the (n, ν)-representation
Two kinds of LO impact factors, generically indicated in Eq. (7) as c i (n, ν, κ i , x i ) ≡ {c (H) (n, ν), c (J) (n, ν)}, where the labels (H) and (J) refer to the jet and the hadron case, respectively, will be considered.
The LO forward-hadron impact factor, c (H) (n, ν), consists in the convolution over the parton longitudinal fraction x between the quark/antiquark/gluon PDFs and the FFs portraying the detected hadron: with C F = (N 2 c −1)/(2N c ) and C A ≡ N c the Casimir factors connected to gluon emission from a quark and from a gluon, respectively. Correspondingly, c (J) (n, ν) is the LO forward-jet impact factor in the (n, ν)-representation: and the f (ν) function in the last term of Eq. (7) reads after fixing the final state, i.e. having selected one of the processes in Fig. 2. The remaining objects are NLO corrections to impact factors,ĉ i (n, ν, κ i , x i ) ≡ {ĉ (H) (n, ν),ĉ (J) (n, ν)}, their expressions being given in Appendix A (hadron) and in Appendix B (jet).

High-energy DGLAP
With the aim of providing a systematic comparison between BFKL-inspired predictions and fixed-order calculations, we derive a DGLAP-like general formula, valid for all the considered processes ( Fig. 2), where the C DGLAP n azimuthal coefficients are introduced as truncation to the O(α 3 s ) order of the corresponding NLA BFKL ones, C NLA n , up to the inclusion of terms beyond the LO accuracy. This allows us to catch the leading-power asymptotic features of a genuine NLO DGLAP description, neglecting at the same time those factors which are quenched by inverse powers of the energy of the partonic subprocess. Although being an alternative (and approximated) way to the standard, fixed-order DGLAP analysis, such procedure results to be adequate in the region of high-rapidity distance, Y , investigated in this work, as well as easy and flexible to implement.
Starting from Eq. (7), one gets the DGLAP limit by expanding the exponentiated BFKL kernel up to the first order in α s , and keeping impact factors at NLO. Our DGLAP master formula reads where the BFKL exponentiated kernel has been replaced by its expansion up to terms proportional to α s (µ R ). We stress that, although the ν-integral in Eq. (17) leads to distributions in the transverse-momentum space rather than ordinary functions, final expressions for the DGLAP azimuthal coefficients become regular when integrated over the final-state phase space, as in Eq. (22). This allows us to safely calculate them via a multidimensional strategy, where ν and phase-space integration are simultaneously performed (for more details, see Section 3.6).

Strong-coupling setting
We use a two-loop running-coupling setup with α s (M Z ) = 0.11707 and five quark flavors, n f , active. Its expression in the MS scheme reads: where It is possible to obtain the corresponding espression for the strong coupling in the MOM scheme, whose definition is related to the 3-gluon vertex (an essential ingredient in the BFKL framework), by performing the following finite renormalization: with where I = −2 1 0 du ln(u) 2.3439 and ξ is a gauge parameter, fixed at zero in the following. The infrared improvement of the running coupling via well-known procedures, as the Webber parametrization [168], here is not needed, since energy scales are always in the perturbative region. They are strictly related to transverse-momentum ranges (see Section 3.2), their large values protecting us from a region where the diffusion pattern [169] (see also Refs. [170,171]) becomes relevant.

PDF and FF selection
Potential sources of uncertainty are expected to arise from the choice of particular PDF and FF sets, rather then other ones. While preliminary tests done using the three most popular NLO PDF parametrizations (MMHT 2014 [172], CT 2014 [173] and NNPDF3.0 [174]) have shown no significant discrepancy in the kinematic regions of our interest, this does not hold anymore for FFs, where the selection of distinct sets leads to non-negligible effects. This feature was pointed out in a first phenomenological analysis [60] on ϕ-averaged cross sections for the inclusive di-hadron production, with the hadrons' ranges tailored on the acceptances of the CMS detector, and then confirmed also for azimuthal correlations (i.e. ratios of phase-space integrated azimuthal coefficients, as explained in Section 3.2) in the case of inclusive hadron-jet emission [63], when the jet is tagged inside the CASTOR backward detector. The simplest strategy to address this issue is in line with the recommendations suggested by the PDF4LHC community [176] for Standard Model phenomenology, and consists in selecting an individual PDF set and taking its convolution with a given FF. This method will be stringently applied in our analysis on the scale optimization (Section 3.4) and on our results for azimuthal distributions (Section 3.5.2), where we will give predictions for the MMHT 2014 [172] NLO PDF together with the NNFF1.0 [175] NLO FF. Conversely, we will adopt a more extensive solution when gauging the effect discussed above for the comparison of our BFKL results with fixed-order DGLAP calculations, which is actually the main outcome of this paper (Section 3.5.1). On the one hand, as for the PDFs, we will choose the NLO PDF4LHC15 100 parametrization [176], which represents a statistical combination of the three sets mentioned before [172][173][174]. On the other hand, we will build our own combination for the FFs by merging predictions obtained with the four following NLO parametrizations: AKK 2008 [177], DSS 2007 [178,179], HKNS 2007 [180] and NNFF1.0 [175]. It is worth to remark that our prescription refers to the average of final results for our observables, while no new FF set has been created. A higher-level analysis, based on the so-called bootstrap replica method (see Ref. [181] for a combined study including fit procedures inspired by neural-network techniques), is employed just in one case, namely our theory-versus-esperiment study for azimuthal correlations in the Mueller-Navelet channel at √ s = 7 TeV (Section 3.3). We make use of a sample of 100 replicas of the central value of the PDF4LHC15 100 set. This collection of replicas is obtained by randomly altering the central value with a Gaussian background featuring the original variance. The extensive application of the replica method to all the considered observables goes beyond our scope and is left to further investigations which will include all the systematic effects.

Jet algorithm and event selection
Several types of reconstruction algorithms can be implemented in the definition of the (NLO) jet impact factor, although the choice of a particular one instead of the others seems not to produce a crucial effect on semi-hard typical observables, as pointed out in Section 2.3 of Ref [37]. The most popular jet-reconstruction functions essentially belong to two distinct classes (for further details, see e.g., Refs. [182,183] and Refs. therein): cone-type and sequential-clustering algorithms (the well known (anti-)κ ⊥ kind [184,185] falls into this last family). We will adopt a simpler version, infrared-safe up to NLO perturbative accuracy and suited to numerical calculations, given in Ref. [21] in terms of the so-called "small-cone" approximation (SCA) [186,187], i.e. for small-jet cone aperture in the rapidity-azimuthal angle plane (see Ref. [22] for a detailed study on different jet algorithms and their small-cone, approximated versions).
A potential issue, highlighted first [41] in the case of Mueller-Navelet jet production (left panel of Fig. 2) but present in any process featuring jet emission in the final state, is related with the experimental event selection in a situation when more than one jet is tagged in a single event. As an example, let us consider events with three objects in the final state, two of them being jets emitted in the forward region (with large positive rapidities, y J,1 > y J,2 > 0), and the third one being a backward object (hadron or jet, in our analyses), with negative rapidity, y O < 0. Traditionally, as in the current CMS analysis, such event would be classified as a unique one, with a single forward jet plus another, backward particle detected, the largest rapidity interval being, Y = y J,1 − y O . This selection criterion is suited to experimental analyses, but it does not match the theoretical definition in NLA BFKL calculations. Examining the derivation of the NLO jet impact factor [18,20], it follows that our calculation describes an inclusive, forward-in-rapidity jet production in the fragmentation region of the parent hadron, whereas possible additional parton radiation is attributed to the inclusive hadron system, X (Eq. (3)). The issue may be clarified either from the experimental side, by adapting the jet-selection paradigm, or from the theoretical side, by improving the NLO calculation, as proposed in Ref. [188].  7) and (17), we consider the integrated coefficients over the phase space for the two emitted objects, O 1,2 (κ 1,2 , y 1,2 ), while their rapidity distance, Y = y 1 − y 2 , is kept fixed: where the superscript [resum] stands indistinctly for LLA, NLA or DGLAP. This allows us to impose and study different ranges in transverse momenta and rapidities, based on realistic kinematic configurations used in experimental analyses at the LHC. Let us introduce the following representation for final-state particles, which will facilitate us to distinguish, process by process, the considered ranges: a) Mueller-Navelet: O 1 (κ 1 , y 1 ) + X + O 2 (κ 2 , y 2 ) ≡ jet(κ J,1 , y J,1 ) + X + jet(κ J,2 , y J,2 ) ; b) hadron-jet: Here labels a), b) and c) refer to the respective panels in Fig. 2, while a slight, process-suitable change of notation in the final-state variables has been made. For all the considered reactions, the maximum values for the transverse momenta are determined by general requirements, based on kinematics. On the one side, in the case of CMS (di-)hadron emission, κ max H is constrained by the lower cutoff of the used FF set and is always fixed at κ max H,CMS 21.5 GeV. On the other side, when (at least) a jet is tagged inside CMS or CASTOR, the value for the upper cutoff, κ max J , follows from the given bounds in rapidity, at fixed center-of-mass energy, √ s and will be individually discussed for each kinematic configuration introduced below. The rapidity interval, Y , is everywhere taken positive: 0 < Y ≤ y max 1 − y min 2 .
We will perform our analysis for √ s = 13 TeV, κ max H ≡ κ max H,CMS 21.5 GeV and κ max J ≡ κ max J,CMS = 60 GeV. This last value has been kept equal with respect to the symmetric case at √ s = 7 TeV (Section 3.1.1).
Results will be shown for √ s = 13 TeV, κ max J,CST 17.68 GeV and κ max J,CMS = 60 GeV.

Final-state observables
The phenomenological observables of our interest fall into the following three classes: 1. ϕ-averaged cross section, C 0 , and the ratio, R NLA/LLA 0 , between C 0 calculated in the full NLA BFKL accuracy and its counterpart in the LLA BFKL approximation; 2. Azimuthal-correlation moments (or simply azimuthal correlations), R nm ≡ C n /C m , with the ratios of the form R n0 among them having the immediate physical interpretation of averaged values of the cosine of n-multiples of ϕ, cos(nϕ) ; 3. Azimuthal distribution of the two emitted objects, which actually represents one of the most directly accessible observables in experimental analyses: On the one side, azimuthal distributions are very favorable observables to be compared with experimental analyses, even more than the azimuthal-correlation moments. This stems from the fact that, since measured distributions cannot cover the whole (2π) range of azimuthal angle due to detector limitations, observables differential on ϕ permit us to dampen the descending accuracy loss. One the other side, numeric calculations of Eq. (24) are not easy to perform, due to the large number of azimuthal coefficients required and due to instabilities in the ν-integration that progressively increase with n. This requires sizeable effort from the numerical side (see Section 3.6.1).

Prelude: Mueller-Navelet jets at 7 TeV LHC, theory versus experiment
In this Section predictions for different azimuthal correlations, R nm ≡ C n /C m , are given for the Mueller-Navelet jet production channel (left panel of Fig. 2) in the symmetric CMS configuration (see Section 3.1.1) and compared with recent CMS data at √ s = 7 TeV [50]. We actually improve the analysis conducted in Ref. [41], by implementing the "exact", numeric BLM scale-optimization method instead of using the two approximated, analytic ones (we are advisedly anticipating details on the scale-optimization issue, discussed in detail in Section 3.4.1). At variance with the previous work [41], where the MSTW 2008 [190] NLO PDF set (precursor of the MMHT 2014 one [172]) was employed, we adopt the NLO PDF4LHC15 100 parametrization [176]  and we estimate the uncertainty of our predictions by calculating the standard deviation at the hand of the replica method (see Section 2.4 for further details). As in Ref. [41], final calculations are done in the MS renormalization scheme. We remind that it is possible to recover the analytic expression for the azimuthal coefficients, C n , in the MS scheme from the corresponding ones, given in the MOM one (Eq. (28)), by making the substitutions (note that the value of the renormalization scale is left unchanged) given in Section 2.3. The two factorization scales and the renormalization one are set equal to the one provided, as outcome for each value of Y , by the BLM procedure: Results for R 10 and R 20 , reported in Fig. 3, unambiguously highlight that the pure LLA calculations overestimate by far the decorrelation between the two tagged jets, whereas the agreement with experimental data definitely improves, when NLA BFKL corrections, together with the BLM prescription, are included. For this process and for the considered transverse-momentum kinematics, the uncertainty coming from the numerical multidimensional integration over the final-state phase space is very small. Therefore, error bands strictly refer to the uncertainty hailing from the use of different PDF members inside the PDF4LHC15 100 collection of 100 replicas. As expected, its effect is very small when ratios of azimuthal coefficients are taken. Ancillary plots below primary panels of Fig. 3 clearly indicate that the relative standard deviation of the replicas' results is smaller than 1% in all cases. Here, nodal points bring information on uncertainty oscillations inherited from the collinear-PDF parametrization.

Need for scale optimization
It is widely recognized that calculations purely based on the BFKL resummation suffer from large higherorder corrections both in the kernel of the gluon Green's function and in the non-universal impact factors. Inter alia, NLA BFKL contributions to the zeroth conformal spin come up with roughly the same order of the corresponding LLA predictions, but with opposite sign. The primary manifestation of this instability lies in the large uncertainties arising from the renormalization (and factorization) scale choice, which hamper any genuine attempt at high-precision calculations. All that calls for some "stabilizing" procedure of the perturbative series, which can consist in: (i ) the inclusion of next-to-NLA terms (a priori unknown), such as those ones prescribed by the renormalization group, as in the all-order [39,[191][192][193][194][195][196][197][198][199][200][201][202][203] or in the consistency-condition [104] collinear-improvement procedures, or by energy-momentum conservation [204]; (ii ) the restriction of undetected-radiation activity to only allow final-state gluons separated by a minimum distance in rapidity (rapidity-veto approach [48,205,206]); (iii ) an "optimal" choice of the values of the energy and renormalization scales, which, though arbitrary within the NLA, can have a sizeable numerical impact through subleading terms.
As for the case (iii ), the most common optimization methods are those ones inspired by: the principle of minimum sensitivity (PMS) [207,208], the fast apparent convergence (FAC) [209][210][211], and the Brodsky-Lepage-Mackenzie method (BLM) [212][213][214][215][216]. Even though the selection of one or the other optimization procedure should not consistently affect the prediction of physical observables, it produces, de facto, a non-negligible impact. Thus, as pointed out in Ref. [41], preference to a particular method should be assigned by grading the agreement with experimental data in a given setup and, consequently, assumed to apply also in other configurations. Predictions derived with BLM turned to be in fair agreement [38,41] with the only kinematic configuration for which experimental analyses are to date performed, i.e. symmetric CMS cuts in the Mueller-Navelet channel (for a detailed comparison of the three listed methods, se the Discussion Section of Ref. [41]). In view of this outcome, the use of the BLM scheme was extended in time to the study of other semi-hard reactions.
In Section 3.4.1 we briefly introduce the key features of BLM, which relies on the removal of the renormalization scale ambiguity by absorbing the non-conformal β 0 -terms into the running coupling. Then, we apply it to the cases of our interest and present results for BLM optimal scale values and cross sections.

BLM scale setting and BFKL cross section
The BLM-optimized renormalization scale, µ BLM R , is the value of µ R that makes the non-conformal, β 0dependent terms entering the expression of the observable of interest vanish. The inspection of the analytic structure of semi-hard observables shows that two groups of non-conformal contribution exist 6 . The first one originates from the NLA BFKL kernel, while the second one from the NLO impact-factor corrections. This makes µ BLM R non-universal, but dependent on the energy of the process.
In order to apply the BLM procedure, a finite renormalization from the MS to the physical MOM scheme needs to be executed (see Section 2.3). Then, the condition for the BLM scale setting for a given azimuthal coefficient, C n , is determined by solving the following integral equation: where It is convenient to introduce the ratio of the BLM to the natural scale suggested by the kinematic of the process, i.e. µ N ≡ √ κ 1 κ 2 , so that C µ R ≡ µ BLM R /µ N , and hunt for those values of C µ R which solve Eq. (25). As a last step, the value of the scale given by the procedure is plugged into expressions for the integrated coefficients, getting the following definition, valid in the NLA BFKL accuracy and in the MOM scheme: wherec i (n, ν, κ i , x i ) are the NLO impact-factor corrections after the subtraction of those terms, entering their expressions, which are proportional to β 0 and can be universally expressed through the LO impact factors: The analogous BLM-MOM expression for the DGLAP case in the high-energy limit reads: Finally, a valid formula for the azimuthal coefficients in the pure LLA BFKL accuracy can be obtained by making in Eq. (8) the substitution (note that the functional form of the strong coupling and the value of the renormalization scale are simultaneously varied) Eq. (25) requires a numerical solution, whose algorithm is implemented in JETHAD (see Section 3.6.1 for further details) and universally holds for any semi-hard final state. Before the advent of JETHAD, some analytic, approximate BLM methods were developed. In those cases (see Eqs. (42)-(46) of Ref. [43]), the BLM scale is treated as a function of the variable ν and is fixed in order to make vanish either the NLO impact-factor β 0 -dependent terms, the NLA BFKL kernel ones, or the whole integrand of Eq. (25). In all the calculations of this work the two factorization scales, µ F 1,2 , are set equal to the renormalization scale µ R , as assumed by most of the PDF parametrizations. All results are obtained in the MOM renormalization scheme (except for the ones presented in Section 3.3, where the MS scheme is selected).
Values of the BLM scales for C 0 and C 1 , for all the three considered production channels (Fig. 2), are presented in Fig. 4. At fixed rapidity distance, Y , the C µ R ratio definitely increases with the number of charged-light hadrons emitted in the final state, up to around 30 ÷ 40 units of natural scales. This phenomenon, already observed in previous analyses on the inclusive di-hadron production [60,62], could be controlled by the behavior of collinear FFs as functions of µ F . More in particular, when FFs are convoluted with PDFs in our LO impact factors (Eq. (14)), the dominant contribution to PDFs in the kinematic sector of our interest is given by the gluon, and thus only the behavior of the gluon FF plays a role. Preliminary studies on FF parametrizations depicting the emissions of different hadron species, including heavy-flavored ones [217], have highlighted how the µ Fdependence of (gluon) FFs can enhance or worsen the stability of our azimuthal coefficients. Smooth-behaved, non-decreasing gluon FF functions, characteristic of heavy-flavored bound states (see, e.g., parametrizations for Λ c baryons [218], charged D * mesons [219][220][221] and b-flavored hadrons [222,223]), have a stabilizing effect. This situation is in some aspects closer to jet emissions, where no FFs are employed. Conversely, light-flavored FF sets, typical of the study presented in this work, lead to an increased sensitivity of predictions on energy scales, as well as to a stronger discrepancy between LLA and NLA. Thus, in our case, the larger scales prescribed by the BLM method come out as a technical consequence of the stability requirement. In addition, the function  , obtained with BLM optimization to the corresponding ones calculated at natural scales, i.e. µ R = µ N ≡ √ κ 1 κ 2 and µ F 1,2 = κ 1,2 , after having checked that the alternative choice for µ F 1,2 , µ F 1 = µ F 2 = µ N , produces almost the same results with respect to the previous one. In particular, Fig. 5 shows the Y -dependence of the R NLA/LLA 0 ratio for the three considered reactions (Fig. 2) in the asymmetric CMS configuration. We notice that NLA corrections, generally with opposite sign with respect LLA results, become larger and larger in absolute value at increasing rapidity interval, thus making the R NLA/LLA 0 smaller and smaller (we stress that the label "NLA" in our plots always stands for LLA plus higher-order terms, and not just the latter ones). This is an expected phenomenon in the BFKL approach which, however, becomes milder when scales are optimized. The adoption of the BLM method leads to a scale choice that permits to mimic the most relevant subleading terms, thus stabilizing the perturbative series. From the operational point of view, this results in a reduction of the distance between the LLA and the NLA, namely it raises the R NLA/LLA 0 ratio in our plots, making it ideally close to one. This is exactly what happens in the Mueller-Navelet case (upper left panel of Fig. 5), while its value at BLM scales exceeds one in the case of hadron detection (remaining panels). The explanation for this, apparently odd behavior, has not to be hunted in the use of BLM, but rather in a combination of two distinct effects, already present in the expression for C 0 given in Eq. (28) and therefore independent from the scale choice. On the one side, going from the MS to the MOM renormalization scheme generates a non-exponentiated, positive extra factor proportional to T conf . On the other side, it is easy to prove that the C gg coefficient in Eq. (A2) gives a large and positive contribution to the hadron NLO impact-factor correction. All that makes C 0 at the NLA larger than the corresponding LLA one before switching any eventual scale-optimization procedure on. For the sake of completeness, the Y -dependence of the R NLA/LLA 0 ratio for the di-jet and the hadron jet case (first two panels of Fig. 2) in the CASTOR-jet configuration is presented in Fig. 6.
Ultimately, the Y -dependence of the NLA ϕ-averaged cross section, C 0 , for all the considered reactions is examined. Final-state configurations molded on the asymmetric CMS (CASTOR-jet) event selection are presented in the left (right) panel of Fig. 7. Here, a definite hierarchy among them is stringently respected, except for values of Y close to the upper bound given by the kinematics: hadron detections dominate in statistics with respect to jet tags.

Hadron-species analysis
In this Section we hunt for possible contributions to the optimal scale values intrinsically coming from the description of the parton-to-hadron fragmentation. Predictions for single hadron-species detections (π ± , K ± , p(p)) are matched to the standard case, where the sum over all species is intended. For the sake of simplicity, we take the inclusive hadron-jet production (central panel of Fig. 2) as reference reaction and consider both the asymmetric CMS and the CASTOR-jet final-state configurations.
The inspection of results of the BLM scales (Fig. 8) for values of the conformal spin, n = 0 and n = 1, indicates that π ± emissions lead to larger scales which are, however, of the same order of the other single-species cases and of the sum one. Therefore, no clues of FF intrinsic effects emerge from our analysis. This can be due to the fact that, while FFs depend on the mass, M H , of the corresponding light-flavored hadron, H, the hard subprocess is mass-independent. Thus, residual logarithmic contributions in M H can survive, generating a mild sensitivity of the scale-choice procedure to the given hadron species.
For completeness of presentation, the Y -dependence of the R NLA/LLA 0 ratio with BLM optimization is shown in Fig. 9, whereas the ϕ-averaged NLA BLM cross section is given in Fig. 10.

BFKL versus DGLAP
In this Section we present and discuss features of the core results of our work. The examination of azimuthalcorrelation moments, given in Section 3.5.1, is extended and accompanied by the analysis of azimuthal distributions, portrayed in Section 3.5.2.

Azimuthal correlations
NLA BFKL predictions for the azimuthal ratios, R nm ≡ C n /C m , in the Mueller-Navelet channel (left panel of Fig. 2) are compared in Fig. 11 with the corresponding high-energy DGLAP ones in the (asymmetric CMS configuration) for √ s = 13 TeV. This extends and completes the study conducted in Refs. [44,45], where a sharp distance between BFKL-resummed and fixed-order calculations came out in the inclusive di-jet hadroproduction with partially asymmetric configurations in the κ-plane, namely 35 GeV < κ J,1 < 60 GeV and 45, 50 GeV < κ J,2 < 60 GeV, and for |y J,1,2 | < 4.7, at √ s = 7 TeV. Here, besides the use of fully asymmetric transverse-momentum cuts (disjoint κ-windows), the main improvement (panels of Fig. 11  proton(p 1 ) + proton(p 2 ) → hadron(κ H , y H ) + X + jet(κ J , y J ) sum π ± K ± p(p) proton(p 1 ) + proton(p 2 ) → hadron(κ H , y H ) + X + jet(κ J , y J ) sum π ± K ± p(p) respective ones in Figs. 2 and 3 of Ref. [44]) stands in the adoption of the "exact" BLM scale-optimization procedure (see Section 3.4.1) instead of approximated, semi-analytic BLM choices, and in a fair reduction of the numeric uncertainty. Analogous results for the di-jet correlations in the more exclusive, CASTOR-jet range, are displayed in Fig. 12. In both the cases, the pattern of the NLA BFKL R n0 series presents a plateau at large rapidity distance, Y , more emphasized in the CASTOR-jet configuration, which visibly evolves into a turn-up for R 10 (upper left panel of Fig. 11 and of Fig. 12). The main reason for this behavior is that the increase of Y values moves the parton longitudinal fractions towards the so-called threshold region, where the energy of the di-jet system approaches the value of the center-of-mass energy, √ s. Hence, PDFs are sounded in ranges close to the end-points of their definitions, where they exhibit large scaling violations and uncertainties 7 . In these configurations, our formalism misses the sizeable effect of threshold double logarithms which enter the perturbative series and have to be resummed to all orders. We note that error bands in the DGLAP case are larger with respect to the resummed ones. This originates from the highly-oscillatory behavior of the ν-integrand  proton(p 1 ) + proton(p 2 ) → hadron(κ H , y H ) + X + jet(κ J , y J ) sum π ± K ± p(p) FIG. 10: Y -dependence of the ϕ-averaged cross section, C0, for inclusive hadron-jet production (central panel in Fig. 2) in the NLA BFKL accuracy, in the asymmetric CMS (left) and in the CASTOR-jet (right) configuration, and for √ s = 13 TeV. Predictions for single hadron-species emissions (π ± , K ± , p(p)) are compared with the standard case, where the sum over all species is taken.
in Eq. (17), not anymore faded by the exponential factor as in the NLA BFKL case, in Eq. (7) (we refer the interested reader to Section 3.6.2 for more details on the numerical-uncertainty estimate of all the presented results).   proton(p 1 ) + proton(p 2 ) → jet(κ J,1 , y J,1 ) + X + jet(κ J,2 , y J,2 ) proton(p 1 ) + proton(p 2 ) → jet(κ J,1 , y J,1 ) + X + jet(κ J,2 , y J,2 ) We show predictions for the R nm ratios in the hadron-jet channel (central panel of Fig. 2) in Fig. 13 (asymmetric CMS ) and in Fig. 14 (CASTOR-jet), whereas analogous results for the di-hadron production (right panel of Fig. 2) are given in Fig. 15 (asymmetric CMS ). Here, presence of identified hadrons in the final state has a dual impact. On the one part, the effect of further splittings (of parton momenta) entering the definition of hadron FFs prevails over the azimuthal recorrelation of the emitted objects at large Y , thus shadowing the turn-up tail observed in the Mueller-Navelet case. On the other part, the convolution between PDFs and FFs in the LO/NLO hadron impact factor(s) (Eqs. (14) and (A1)) stabilizes the oscillations of the ν-integrand, thus reducing the numeric uncertainty in the calculation of the high-energy DGLAP series, whose error bands are now of the same magnitude of the NLA BFKL ones.
The unphysical effect that brings R 10 above one for small values of the rapidity interval, Y , is well known in the context of semi-hard reactions and has a straightforward explanation. In this kinematic limit, the hardsubprocess energy, √ŝ , is small and the employment of BFKL resummation, which systematically neglects terms damped by powers ofŝ, becomes inadequate. The fact that, in the DGLAP case, values of R 10 larger than one are still present for larger Y -values originates from the nature of our approach. Being it, de facto, an high-energy calculation where the NLA resummation is truncated to a given order in the running coupling (see Section 2.2), collinear contaminations, which are largely present at low values of the conformal spin and at small-Y , may survive and manifest themselves also in the remaining part of the Y -range.
The overall outcome of the analysis conducted in this Section is a clear separation between BFKL and DGLAP predictions on azimuthal correlations of the two detected objects. This effects holds for all the considered proton(p 1 ) + proton(p 2 ) → hadron(κ H , y H ) + X + jet(κ J , y J ) proton(p 1 ) + proton(p 2 ) → hadron(κ H , y H ) + X + jet(κ J , y J ) Results in the asymmetric CMS (CASTOR-jet) configuration are given in upper (lower) panels.
reactions and becomes more and more evident as the rapidity interval, Y , raises. The found pattern matches in toto the farsighted idea of Mueller and Navelet [33], namely that the wealth of undetected gluons radiated in the final state, theoretically described at the hand of the high-energy resummation, markedly heighten the decorrelation in the azimuthal plane between the emitted particles. This makes a substantial difference with respect to the DGLAP case, where only a limited number of gluon emissions, fixed by the truncation order of the perturbatuve series, is allowed. The adoption of asymmetric cuts for the transverse momenta in the final state fades the Born contribution, thus spotlighting the discrepancy between the two approaches.

Azimuthal distribution
Here we give predictions for azimuthal distributions, which, as anticipated, are directly accessible observables in the experimental studies. The first examination of these quantities was performed few years ago [37] in the context of the Mueller-Navelet jet production, for symmetric ranges of transverse momenta, and with partial and full NLA BFKL accuracy.
Results for the azimuthal distribution in the di-jet channel are shown with NLA BFKL (high-energy DGLAP) accuracy in the left (right) panels of Fig. 16, for three distinct values of the rapidity interval, Y . Predictions in the asymmetric CMS configuration are given in upper panels, while the lower ones refer to the CASTOR-jet proton(p 1 ) + proton(p 2 ) → hadron(κ H,1 , y H,1 ) + X + hadron(κ H,2 , y H,2 ) selection. Fig. 17 shows, in the same way, the azimuthal distribution in the hadron-jet emission for two distinct values of Y , whareas plots of Fig. 18 correspond to the di-hadron distribution with asymmetric CMS cuts, for two distinct values of Y .
The peculiar shape of our series represent a further manifestation of the BFKL-versus-DGLAP "dicotomy". All distributions feature a clear peak in correspondence of the value of the azimuthal-angle difference, ϕ ≡ ϕ 1 − ϕ 2 − π, for which the two final-state objects are emitted back-to-back, i.e. ϕ = 0. Let us start by considering, for each panel, the first entry, which corresponds to the lower value of Y among the selected ones. Both in the BFKL and in the DGLAP cases, this peak dominates over the larger-Y ones, the DGLAP being narrower in comparison to the BFKL one. Then, when one moves towards the larger values of Y , peaks visibly shrink and widths moderately widen in the BFKL case, while that switch is much less evident for the DGLAP series. This means that, when the rapidity interval grows, the number of back-to-back events predicted by BFKL decreases, while it remains relatively unchanged according to DGLAP. Our outcome is in perfect agreement with the main statement raised in Section 3.5.1.
As a side consideration, we notice that trends obtained in the CASTOR-jet ranges always show a slightly larger decorrelation with respect to asymmetric CMS patterns (lower panels versus upper ones in Figs. 16 and 17). The added value of the more exclusive kinematic configurations offered by tagging a jet in the CASTOR ultra-backward detector translates in an intrinsic decorrelation gained by the final state.
All these features brace the message that kinematic configurations suitable to heighten high-energy effects in the context of the LHC phenomenology exist and have been effectively detected.
3.6. Numerical strategy 3.6.1. JETHAD: an object-based, process-independent interface All numerical calculations were performed using JETHAD, a hybrid Fortran2008/Python3 modular code we recently developed, suited for the computation of cross sections and related observables for inclusive semi-hard reactions. Natively equipped with performance acceleration, by making extensive use of parallel computing techniques, and interfaced with the most advanced (multi)dimensional integration routines, JETHAD allowed us to dynamically select the best integration algorithm, depending on the behavior of the considered integrand. Monte-Carlo inspired integrators with variance reduction via importance sampling, like the Vegas routine [235] as given in its concurrent version via the Cuba package 4.2 [236,237], were primarily selected to perform multidimensional integrations needed to calculate the C DGLAP n coefficients (Eq. (17)) in the Mueller-Navelet channel (left panel of Fig. 2). Conversely, adaptive-quadrature based functions, like DADMUL and WGauss as implemented in the last version of Cern program libraries [238], were mainly preferred for the computation of all observables related to (di-)hadron emission (central and right panels of Fig. 2) and for the one-dimensional integration over the longitudinal momentum fraction, ζ, entering the expression for the NLO hadron/jet impact factors (Eqs. (A1) and (B1)), respectively. All PDF parametrizations, as well as the NNFF1.0 FF set, were calculated through the Les Houches Accord PDF interpolator (LHAPDF) 6.2.1 [239], while native routines for the remaining FFs were direclty linked to the corresponding module in our code. In order to dynamically select the considered reaction through a common interface, a structure-based smart-managment system, where physical final-state particles are described in terms of object prototypes (i.e., Fortran structures), was incorporated inside JETHAD. Particle objects carry all information about basic and kinematic properties of their physical counterparts, from mass and charge, to kinematic ranges and rapidity tag. They are first loaded by the JETHAD user routine from the master database through a specific particle generator routine (custom-particle generation is allowed too). Then, they are cloned to the final-state object array and thus injected from the integrand routine, differential on the final-state variables, to the respective impact-factor module through the impactfactor controller. To the strong flexibility in the final-state generation, a wide choice for the initial-state selection corresponds. Thanks to a peculiar particle-ascendancy structure attribute, JETHAD is able indeed to recognize if an object is hadroproduced or emitted as subproduct of a leptonic interaction, and automatically determines which modules need to be initialized (PDFs, FFs, unintegrated densities, etc...), breaking down computingtime lags. Hence, JETHAD comes as an object-based interface, completely independent on the reaction being investigated. While inspired by the BFKL phenomenology, it is possible to perform calculations of the same observables in different approaches, by creating new, customized routines, which can be linked to the core structure of the code through a native point-to-routine system, thus making JETHAD a general, HEP-purposed tool. We pursue, as a medium-term goal, to release a first public version soon, providing the scientific community with a standard software for the analysis of inclusive semi-hard processes. Another code, LExA, based on the same framework and suited to the study of exclusive semi-hard reactions, is currently at an early-development stage.
A robust improvement of our technology would consist in successfully interfacing these codes with already existing software, suited to high-energy/small-x studies. An incomplete list includes: novel software for NLA studies [240][241][242][243] of Mueller-Tang (alias jet-gap-jet) configurations [244], the BFKLex Monte Carlo [245][246][247], designed to investigate the high-energy jet production in the multi-Regge limit, the KaTie generator [248] of off-shell matrix elements, and the TMDlib library [249], where different models of the small-x UGD are collected.

Uncertainty estimate
The most relevant uncertainty comes from the numerical 4-dimensional integration over the transverse momenta, κ 1,2 , of the two final-state objects, the rapidity, y 1/2 , of one of them (the other one is fixed by the condition Y = y 1 − y 2 enforced in the definition of the integrated azimuthal coefficients (22)), and over ν. Its effect was directly estimated and given as output by master integrator routines. Further, secondary sources of uncertainty, are respectively: the one-dimensional integration over the parton longitudinal fraction, x, needed to perform the convolution between PDFs and, eventually, FFs in the LO/NLO hadron impact factors (Eq. (14)) for LO, Eq. (A1)), the one-dimensional integration over the longitudinal momentum fraction, ζ, in the NLO impact factor corrections (Eq. (A1) for hadrons, Eq. (B1) for jets), and the upper cutoff in the numerical integral over ν. While the first two ones turn out to be negligible with respect to the multidimensional integration, the last one deserves particular attention. As pointed out in Section 3.3 of Ref. [44], the C DGLAP n coefficients are expected to exhibit a stronger sensitivity to the upper cutoff of the ν-integration, ν max , due to the fact that oscillations rising in the ν-integrand in Eq. (17) are not quenched by the exponential factor as in the NLA and LLA BFKL expressions (Eqs. (7) and (8), respectively). This turns to be particularly true in the Mueller-Navalet channel (left panel of Fig. 2), while the peculiar structure of the x-integrand in the LO hadron impact factor (Eq. (14)) generates, as the net result, a significant damping of the ν-associated oscillations both in the DGLAP and in the NLA/LLA BFKL cases, allowing us to use smaller values of ν max in hadron-jet and di-hadron production (central and right panels of Fig. 2) with respect to the di-jet case. Values for ν max , given below (Tab. I), have been taken after checking that the uncertainty on the R nm ratios, coming from raising them by a value of 10, is negligible if compared to the valued error of the multidimensional integration.
Error bands of all predictions presented in this work are given in terms of the uncertainty on the final-state integration, combined, in case of hadron emission (central and right panels of Fig. 2), with the one coming from averaging over different FF sets (see Section 2.4 for further details). An exception is representend by the theory-versus-experiment analysis in the Mueller-Navelet channel (Section 3.3), where error bands in the two panels of Fig. 3 show the standard deviation of predictions for the R nm ratios calculated via the replica method.  7) and (8)) and in the high-energy DGLAP limit (Eq. (17)), for all considered reactions (Fig. 2).

CLOSING STATEMENTS
We brought evidence that distinctive signals of the high-energy resummation emerge in the NLA description of different semi-hard reactions. Taking advantage of the record energies and of the exclusive kinematic configurations provided by the LHC, their effects can be effectively disengaged from the ones arising from (pure) fixed-order, DGLAP-inspired calculations. In this direction, the use of (completely) asymmetric intervals for the transverse momenta of the detected objects plays a crucial role and definitely needs to be taken into account in the forthcoming experimental analyses on forward/backward final-state emissions. Among them, stringent measurements of the azimuthal-angle averaged cross section, C 0 , turn out to be essential not only in the discrimination of BFKL from other theoretical approaches, but also in the assessment of the intrinsic ambiguities of the given approach.
With the goal of providing, in the medium-term future, comparisons with genuine fixed-order calculations, the combined effect of a selection of potential uncertainties was gauged. Our preference fell on those uncertainties which are typically included in collinear-physics phenomenology, but they turn out to be novel in the semi-hard one. More in particular, we studied the sensitivity of azimuthal-correlation moments on PDF uncertainty via the so-called replica method [181] (this was done for the Mueller-Navelet channel at √ s = 7 TeV (Section 3.3)) and on FF uncertainty, by averaging among four distinct FF sets. A detailed analysis of all the potential sources of uncertainty is postponed to future studies, when data for our disjoint κ-windows will be available. Then, two complementary paths should be traced and concurrently treated in the context of our studies.
Pursuing the goal of dealing with more exclusive final states, the investigation of heavy-flavored emissions has been taken into account in our program [98][99][100], with the medium-term target of including the case of quarkonia (for a recent work on the inclusive J/ψ-plus-jet hadroproduction, see Ref. [101]). Here, the ultimate task relies upon a profound examination of the theoretical production mechanism [250][251][252] for these quark-bound states, which has not yet been completely outfound, and on the grading of the most popular models proposed so far [253][254][255]. Various benefits can be gained from considering single forward emissions. First, the experimental statistics appreciably increases. Second, it allows us to probe different frameworks for the unintegrated gluon distribution (UGD), including the ones [256,257] whose definition suitably embodies non-perturbative inputs driven by the transverse-momentum-dependent (TMD) factorization (see, e.g., Refs. [129,258,259] and references therein for an overview on the general framework, Refs.  for quite recent applications) together with effective small-x effects.
The analysis of more differential distributions covering broader kinematic ranges requires an ineludible effort into the enhancement of our formalism to accommodate other resummation mechanisms. A major outcome of a recent work on the inclusive hadroproduction of a Higgs boson and a jet well separated in rapidity [102] is that an exhaustive study of the Higgs transverse-momentum distribution lies on the exploration of neighboring kinematic regions, each of them representing a preferred channel to probe a specific resummation. That calls for a net overhaul of our framework (first from the analytic and then from the numeric point of view) which would not anymore rely only on the pure BFKL approach, but should rather evolve into an underlying staging where distinct resummations are plugged-in and play their part.
One step forward is represented by the Altarelli-Ball-Forte (ABF) formalism [288][289][290][291][292][293][294], where DGLAP and BFKL inputs are combined to improve the perturbative accuracy of the resummed series, by imposing consistency conditions (duality aspects), symmetrizing the BFKL kernel in the (anti-)collinear phase-space regions, and incorporating those contributions to running coupling which affect the small-x divergences. Thence, the interplay between collinear and high-energy factorization allows us to perform the resummation of coefficient functions, whereas the resummation of splitting functions is endowed by enforcing the consistency between the two evolution equations. The first determination of proton PDFs where NLO and next-to-NLO fixed-order calculations are supplemented by the NLA small-x resummation has been recently realized [158,295,296], by making use of the HELL [297,298] numerical code interfaced to the APFEL [299-301] PDF-evolution library.
Another intriguing possibility is represented by the Catani-Ciafaloni-Fiorani-Marchesini (CCFM) branching scheme [302][303][304][305]. Enforcing angular ordering (coherence) of soft-parton emission from the large-x to the small-x kinematic ranges, the CCFM framework provides with a unified evolution pattern for unintegrated gluon densities. In the totally inclusive configuration, this evolution interpolates between the DGLAP equation at moderate-x and the BFKL one at small-x. For large values of x and for high virtualities, CCFM dynamics matches DGLAP evolution, whereas in the asymptotic-energy limit it almost corresponds to BFKL (see Refs. [306][307][308]). It essentialy builds on the sum of ladder diagrams with angular ordering along the chain. Parton transverse momentum is generated via the genuine recoil effect due to gluon radiation controlling the UGD evolution. The direct employment of CCFM to evolution equations for parton distributions was proposed by Jan Kwieciński [309] through the so-called CCFM-K equations. The CCFM formalism was then generalized [310] in order to account for non-linear effects in the gluon evolution, thus giving us a chance to gauge the impact of parton saturation on exclusive observables. On the phenomenological side, a description of the Drell-Yan production at the hand of CCFM-K-evolved distributions was recently proposed [311].