Towards stability of NLO corrections in high-energy factorization via modified multi-Regge kinematics approximation

The perturbatively-stable scheme of Next-to-Leading order (NLO) calculations of cross-sections for multi-scale hard-processes in DIS-like kinematics is developed in the framework of High-Energy Factorization. The evolution equation for unintegrated PDF, which resums log 1/z-corrections to the coefficient function in the Leading Logarithmic approximation together with a certain subset of Next-to-Leading Logarithmic and Next- to-Leading Power corrections, necessary for the perturbative stability of the formalism, is formulated and solved in the Doubly-Logarithmic approximation. An example of DIS-like process, induced by the operator tr [GμνGμν ], which is sensitive to gluon PDF already in the LO, is studied. Moderate (O(20%)) NLO corrections to the inclusive structure function are found at small xB< 10−4, while for the pT -spectrum of a leading jet in the considered process, NLO corrections are small (< O(20%)) and LO of kT -factorization is a good approximation. The approach can be straightforwardly extended to the case of multi-scale hard processes in pp-collisions at high energies.


Introduction
The High-Energy or k T -factorization formalism, first introduced in [1], later had been developed [2,3] as a tool to resum higher-order corrections to coefficient functions of Collinear Parton Model, enhanced by large logarithms log(1/z) of light-cone momentum fraction z, of a parton entering into a hard subprocess, relative to the characteristic lightcone momentum component of a final-state of interest. This kind of corrections become more and more important with increasing collision energy, since more phase-space for additional semi-hard emissions opens up. These emissions generate a transverse momentum recoil, which greatly affects kinematic distributions of the final-state of interest, e.g. dior multi-jet system [4][5][6], pair of heavy-flavoured mesons [7,8] or heavy quarkonia [9]. Therefore, k T -factorization calculation serves as an interesting alternative to fixed-order calculations of such observables in Collinear Parton Model (CPM) or with conventional Parton Showers (PS) (see ref. [10] for the review). As one can see, e.g. from references cited above, the k T -factorization canculation with a judicious choice of unintegrated-Parton Distribution Function (UPDF) quite often leads to a good description of various correlation observables already in the leading order (LO) in α s , as opposed to situation in CPM, where e.g. the description of ∆φ-spectrum can be quite poor even at Next-to-Leading order (NLO) and significant initial-state PS effects have to be taken into account.
In the present paper we continue development of the technique of NLO calculations in the gauge-invariant scheme of High-Energy Factorization, based on Lipatov's gaugeinvariant Effective Field Theory(EFT) for Multi-Regge processes in QCD [40,41]. Following refs. [4,7] we call this scheme -the Parton Reggeization Approach (PRA). The version of Modified Multi-Regge Kinematics (MMRK) approximation for QCD amplitudes with multiple real emissions was proposed in ref. [7] to justify the use of KMRW UPDFs together with tree-level amplitudes from High-Energy EFT [40,41]. Below we will write-down the evolution equation for UPDF, based on MMRK-approximation, which besides leading log(1/z)-terms allows one to resum a subset of subleading logarithmic and O(z) powersuppressed corrections to the UPDF. Then we perform an exploratory NLO calculation for the coefficient function of Deep-Inelastic-Scattering-like subprocess, driven by a gaugeinvariant operator tr[G µν G µν ] (where G µν is a Non-Abelian field-strength tensor), which couples to gluons already in the LO in α s . In CPM, the coefficients functions of this process are known up to O(α 3 s ) [42] and starting from NNLO they contain doubly-logarithmic terms ∝ α n s log 2n (1/z) origin of which had been explained in k T -factorization [43]. Besides inclusive DIS cross-section (or "structure-function") we also study the cross-section of production of a leading jet in this process in LO and NLO of PRA. The standard MRK approximation leads to large negative NLO corrections, which for x B > 10 −4 and Q 1 GeV turns NLO cross-section negative, signaling a severe perturbative instability. On the contrary, consistent implementation of our MMRK approximation results in a moderate NLO corrections, which shows, that it solves the major part of the problem of perturbative instability of BFKL formalism at NLO. We trace this results back to the improved treatment of region of initial-state collinear singularity (DGLAP region) in the MMRK approximation. In summary, we have come-up with a practical and manifestly perturbatively-stable recipe of NLO calculations in PRA, which allows one to improve accuracy of the predictions and establish the boundaries of applicability of the approach through the smallness of NLO correction.

JHEP08(2020)055
The paper is organized as follows: in section 2 we formulate the basic formalism of PRA for the particular process we have chosen to study and derive the evolution equation for UPDF in MRK approximation, then in section 3 we formulate our MMRK approximation, analyze it's performance in comparison to an exact QCD amplitude with one additional emission and write-down UPDF-evolution equation in MMRK approximation and corresponding NLO double-counting subtraction terms; in section 4 we describe our phase-space slicing strategy, compute corresponding analytic integrals and double-counting subtraction integral in the soft limit; in section 5 we recall the virtual part of NLO-correction under consideration, computed in ref. [44], and derive corresponding virtual subtraction terms; finally in section 6 we present and discuss some numerical results and formulate our conclusions. In the appendix A we derive an approximate doubly-logarithmic solution for our UPDF evolution equation, which we use for illustrative numerical calculations throughout this paper and in the appendix B we take a few iterations of evolution kernel (2.30) to demonstrate it's properties.

Basic formalism and UPDF evolution in MRK approximation
To simplify our presentation, we will always refer to a particular example of hard process -the DIS-like process (momenta of particles are given in parentheses): initiated by the gauge-invariant local QCD operator where λ is a coupling to an external source and G µν = −i [D µ , D ν ] /g s is a field-strength tensor of QCD with covariant derivative expressed as D µ = ∂ µ + ig s A µ and Hermitian gluon field A µ = A a µ T a , where T a are generators of SU(N c ). In the present paper we will concentrate on the case of pure gluodynamics, i.e. the theory with n F = 0. Within Standard Model, the operator (2.2) can be understood as an effective coupling of gluons to a Higgs boson through a loop of very heavy quark, and therefore process (2.1) can be visualized as a Higgs-exchange contribution to the usual electron-proton DIS. Of course phenomenologically such contribution is negligible, but since the operator (2.2) couples to gluons through the two-gluon vertex: where gluon momenta k 1,2 are incoming and q + k 1 + k 2 = 0, as well as through three and four-gluon vertices, proportional to corresponding vertices of QCD with q +k 1 +. . .+k 3 = 0 and q + k 1 + . . . + k 4 = 0, the process (2.1) have proven to be a useful tool for the formal studies in QCD, probing various aspects of evolution of gluon PDF, see e.g. refs. [42,45].
In the present paper we will consider the dimensionless inclusive "structure-function" of the process (2.1), which depends on usual DIS kinematic variables Q 2 = −q 2 and JHEP08(2020)055 x B = Q 2 /2(qP ) and in the LO of CPM is simply equal to where f g is a usual collinear PDF. In our numerical calculations below, we will put the factor πλ 2 /4 = 1. For simplicity, throughout this paper we choose to work in a center-ofmomentum frame of P and q, where the light-cone components 1 of these momenta can be expressed as: i.e. at x B 1 momentum q has large positive (forward) rapidity, while proton flies in the negative direction.
The general expression for inclusive DIS structure-function in CPM is well-known: where the coefficient function C is computed perturbatively as a power-series in a s = α s (µ 2 )/(2π) and the first (leading-twist) term of eq. (2.6) is valid up to corrections suppressed as (Q 2 ) −ν with ν > 0. The k T -factorization hypothesis [1][2][3] states, that higher-order corrections to the coefficient function, enhanced by log(1/z), can be further factorized-out: 7) where new coefficient function H is free from potentially large log(1/x 1 )-corrections, by µ Y we have denoted additional scale which arises due to factorization (2.7) and spurious dependence on x B is introduced to make sure, that variable x 1 has the same kinematical meaning in eq. (2.7) and eq. (2.8) below. The statement (2.7) is proven in QCD in leadingpower approximation in z 1 (i.e. up to O(z)-terms) for the series of Leading-Logarithmic (LL, ∝ [a s log k (1/z)] n with k = 1 or 2 depeding on a processes) [3,[14][15][16] and Next-to-Leading Logarithmic (NLL, ∝ a s [a s log k (1/z)] n ) [17][18][19] corrections. Evolution factor C is always single-logarithmic w.r.t. log(1/z) at leading power in z, however additional power of log(1/z) per a s can be generated by transverse-momentum integration in eq. (2.7) [43].
The situation at subleading power is significantly more complicated, see e.g. refs. [46,47], with doubly-logarithmic corrections arising for some quantities in LLA [48][49][50], however this corrections still can be organized into a sum of terms of a form (2.7) with different 1 For any four-momentum k we define Sudakov decomposition as k µ = n µ + k− + n µ − k+ /2 + k µ T with k± = n±k, n 2 ± = 0, n+n− = 2 and n±kT = 0, so that k 2 = k+k− − k 2 T and we do not distinguish covariant and contravariant light-cone components: k± = k ± . JHEP08(2020)055 coefficient functions H and evolution factors C. In the present paper, we work under assumption, that there exist a series of subleading-power (O(z)) corrections which can be written in a form (2.7) with the leading-power coefficient function H, so that all corrections are absorbed into C, and we assume that this series is numerically dominant for most of inclusive quantities. We make such an assumption, instead of systematically going orderby-order in z-expansion, because phenomenologically such an expansion can hardly be expected to be quickly-convergent, since higher-order corrections in QCD typically contain functions like log(1−z) or 1/(1−z) + with rather slowly-convergent expansion around z = 0.
Substituting eq. (2.7) to eq. (2.6), changing the order of integrals in x 1 and z and making the substitution z → zx B /x 1 one arrives at a standard High-Energy Factorization formula [1][2][3][4]7]: and we have rewritten the coefficient-function H in terms of "squared matrix element" (ME) |M PRA | 2 , which at leading power in z is computed as: where A EFT is an amputated Green's function of Lipatov's EFT for Multi-Regge processes in QCD [40] with one incoming Reggeized gluon R − carrying four-momentum q µ 1 = q − 1 n µ + /2 + q µ T 1 = x 1 P − n µ + /2 + q µ T 1 and some partonic final-state with total four-momentum p M . The convention (2.10) is introduced (see e.g. ref. [51]) to ensure, that in the onshell limit |q T 1 | → 0 the PRA squared ME reproduces the corresponding squared ME of CPM with Reggeized gluon substituted by an on-shell gluon with four-momentum n µ + q − 1 /2, averaged over it's color and helicity of this gluon: where φ 1 is an azimuthal angle of q T 1 . Also in eq. (2.8) we have introduced a flux-factor 2Sx 1 which is just a matter of convention, and one have to integrate over the final-state of M with the usual Lorentz-invariant phase-space volume element dΠ M . The LO PRA subprocess for the process (2.1) is: with the following squared ME, derived from vertex (2.3) and LO R − → g-mixing vertex ∆ ab −µ (q) = (−iq 2 )n + µ δ ab of EFT [40] (see e.g. eq. (13) in ref. [44]): (2.13)
To set our notation, below we will derive the real-emission term of evolution equation for C in the standard MRK approximation, which eventually will coincide with the LO BFKL equation with real emissions ordered in physical rapidity, but rewritten in terms of light-cone momentum fractions z and transverse momenta. To this end let's consider the k T -factorized expression for the contribution to the coefficient function of CPM (2.7) with n − 1 real emissions already factorized into evolution factor C n−1 and emission of one additional gluon with four-momentum k n in the PRA matrix element (see the figure 1): whereq µ 1 = (x 1 P − )n µ + /2 +q µ T 1 and we have introduced an intermediate "t-channel" fourmomentum q 1 . To kinematically factorize-out additional emission one performs the following approximation: (2.17) i.e. neglects the "small" light-cone component in the hard process, thus the light-cone components of q 1 are: and x 1 = q − 1 /P − 1 , in terms of which, the longitudinal measure of integration becomes: figure 1. With approximation (2.17) the t-channel momentum transfer in the hard process is equal to: The approximation (2.19) becomes accurate in the limit Q 2 → 0. In the Regge limit z n 1 one can put t −q 2 T 1 , but the latter approximation quickly degrades with an JHEP08(2020)055 increase of z n . In the kinematic constraint approach [52,53] to approximately take into account large collinearly-enhanced corrections to BFKL-evolution, one cuts-off the region of real-emission phase-space, where t −q 2 T 1 -approximation is no longer valid, i.e. one rejects the emissions with: Retaining some realistic approximation for t-channel momentum transfer, analogous to eq. (2.19), will provide the smooth cutoff in the same region of phase-space, thus leading to resummation of essentially the same class of higher-order corrections to BFKL-kernel as in the kinematic constraint approach.
The squared matrix element with one additional gluon emission, in the Regge limit z n 1 factorizes as follows: where factors in square brackets are introduced due to normalization prescription (2.10) and the factor 1/(2t) 2 is the squared tree-level propagator of Reggeized gluon. In the standard MRK-approximation, the factor in curly brackets is put to one as it was discussed above.
) amputated Green's function in the EFT [40], i.e. the square of Lipatov's vertex [14]: Collecting all pieces together, one can rewrite eq. (2.16) as: is the real-emission term of MRK evolution equation, which we have been looking for.

JHEP08(2020)055
The limits of z n -integration have to be defined in eq. (2.21). The lower limit is x due to conservation of (−)-component of momentum, while integrating up to z n = 1 will lead to non-regularized rapidity divergence due to denominator 1/(1 − z n ). Technically, the divergence arises because approximation (2.17) violates conservation of (+)-momentum component and hence additional emission is allowed to go arbitrarily forward (in the qdirection) in rapidity. Demanding, that rapidity of this emission is cut-off at some value Y µ one obtains the condition: where function ∆ is familiar from the definition of KMRW UPDF [36][37][38] while the rapidityscale µ Y is defined by the relation q − 1 = µ Y e −Yµ and evolution factor will depend on rapidity scale from now on. In DIS kinematics the good choice for Y µ is the rapidity of a parton emitted in the LO PRA subprocess (2.12): which removes large-logarithmic terms ∝ Y H form coefficient function H at NLO, as we will show in section 4. The rapidity of j-th gluon in the evolution cascade is given in terms of light-cone momentum fraction of an adjacent t-channel parton z j by: Applying the latter result to the rapidity-ordering condition for the next emission -y n > y n−1 one obtains the following choice of rapidity scale for the evolution-factor C n−1 in the integrand of eq. (2.21): µ Thus we have completely obtained all the details of real-emission part of our MRK evolution kernel.
Also, eq. (2.24) allows one to interpret the z n -integration measure in the eq. (2.21) as integration over rapidity y n : dz n z n (1 − z n ) = dy n .
To precisely write-down the D-dimensional virtual part of the evolution kernel, one can use the one-loop correction to Reggeized gluon propagator (with Born propagator 1/(2t 1 ) factorized-away) of refs. [54][55][56] or, equivalently, the eq. (53) in ref. [44] which we reproduce here for the later reference in section 5: where 1 is the parameter of regularization for rapidity divergences in loop integrals,

JHEP08(2020)055
which has been used in ref. [44]. The logarithm log(1/r) can be identified with the rapidity difference between two adjacent real emissions (see the discussion in the end of section 1 of ref. [44]), and therefore, the virtual part of the evolution equation will be proportional to the coefficient in front of this logarithm -the one-loop Regge-trajectory of a gluon: In rapidity-space, one iteration of the real-emission and virtual parts of the evolution kernel has the form [14][15][16] (see also review [57] or a textbook [58]): (2.27) where the factor of two takes into account one-loop contributions from amplitude and complex-conjugate amplitude.
The starting iteration of the evolution (the LO partonic "impact-factor" of the target, in BFKL terminology) is given in rapidity and x-space by: (2.28) and the evolution factor to all orders in α s in x-space is: is the perturbative initial condition. As we already have found in eq. (2.21) for the case of real-emission contribution, the eq. (2.27) can be equivalently rewritten in terms of light-cone fraction z: where the particualr rapidity-scale choice in the virtual part: µ ) is uniquely determined by requirements of exponentiation of the virtual part of the evolution and/or cancellation of infra-red divergences between real and virtual parts, as it is shown in the appendix B. eq. (2.30) is the final form of our MRK evolution equation for the UPDF evolution factor C.
It is well-known, that infra-red divergences cancel to all orders in α s in eq. (2.27) and hence they should also cancel in eq. (2.30) as it is discussed in more detail in the appendix B. But when one takes the q T 1 -convolution of the evolution factor C with the JHEP08(2020)055 coefficient-function H in factorization formula (2.7), the collinear divergences ∝ (α s / ) n are generated to all orders in α s . The latter should be absorbed by usual renormalization of collinear PDF in eq. (2.9) as it was first systematically done in ref. [3]. It is most convenient to perform this procedure in Fourier-conjugate x T -space, because in this space all collinear divergences are contained in the evolution factor C(x T ). We do this in appendix A for the simplified version of eq. (2.30) which strictly neglects all O(z)-corrections in the kernel and thus does not depend on the scale µ Y . The UPDF obtained from doubly-logarithmic solution of this simplified equation is used for illustrative numerical calculations throughout this paper.

Modified MRK approximation: subtraction terms and UPDF evolution
To demonstrate the necessity of improving MRK-approximation to reach the stability of NLO corrections in k T -factorization let us consider again, the real-emission NLO correction to the cross-section of the process (2.1), which is given by PRA subprocess: Let us introduce the convenient variable: which together with k T 1 , k T 2 , Q 2 and x B completely parametrizes exact 2 → 2 kinematics of this subprocess. In terms of the latter variable, the contribution of subprocess (3.1) to the SF is given by: where we have introduced integrand-function -w and reduced ME -f which are related with UPDF and squared ME of the subprocess (3.1) as follows: and From eq. (3.6) follows the simple expression for rapidity difference between gluons:

JHEP08(2020)055
which tells us, that for fixed transverse momenta, limitẑ → 0 corresponds to y 1 > y 2 (t-channel Regge limit, −t/s 1), while in the limitẑ → 1 one has y 2 > y 1 (u-channel Regge limit, −u/s 1). In general, the substitution: corresponds to permutation of final-state gluons. The squared ME of the subprocess (3.1) can be straightforwardly obtained using Feynman rules of EFT [40] (see e.g. refs. [4,7,44,59] for the detailed presentation) and is rather long and non-instructive expression, so we refrain from presenting it here. Some relevant limits of it are given below and in section 4.
As it was shown in section 2 the UPDF-evolution is obtained by factorizing-out an additional gluon emission from subprocess (3.1). Hence, to remove double-counting of eq. (3.3) with the evolution, one has to subtract the corresponding approximation for the squared ME from the exact integrand w. For standard MRK approximation, the t-channel subtraction term is given by eq. (3.4) with the following reduced ME: and UPDF evaluated at , in accordance with approximation (2.17). The θ-function in eq. (3.10) enforces the rapidity-ordering condition (2.22) and with µ Y = (Q 2 + k 2 T 1 )/|k T 1 | it is equivalent to the condition y 1 > y 2 , see eq. (3.8). The u-channel subtraction term is obtained from (3.10) via substitution (3.9).
Clearly, the NLO correction will be smaller if the subtraction term provides a better approximation to an exact squared ME. Since evolution equation is constructed by iterating the same approximation, improvement of the subtraction term will also make the evolution to capture more physics. In fact, as we will see in section 6, to obtain meaningful physical results it is crucial to come-up with better approximations to an exact ME, than eq. (3.10) can provide. The most important phase-space region, where improvements are necessary, is the DGLAP-region: q 2 integration over which at fixed q 2 T 1 generates the contribution enhanced by log(Q 2 /q 2 T 1 ). The latter large logarithm, when integrated over q 2 T 1 with UPDF, leads to sizeable numerical effects. For the squared ME, the DGLAP limit q 2 Q 2 asympthotics. Due to eq. (2.11), general (initial-state) collinear factorization theorem for squared MEs in QCD is applicable in this case (see e.g. eq. (4.9) of ref. [60]) and reduced ME is given by: . eq. (3.11) is a non-trivial function ofẑ, but in eq. (3.10) this function is approximated by a constant. Our goal is to improve this situation leaving the cancellation of infra-red divergences in eq. (2.30) intact. To this end we restore the (q 2 T 1 /t) 2 ("propagator-factor") in eq. with the following approximation for the t-channel momentum transfer (compare it with eq. (2.19)): so that reduced ME for the subtraction term in the modified MRK-approximation takes the form: (3.13) In the on-shell limit eq. (3.13) reproduces eq. (3.11) up to O(ẑ 2 )-terms, so we have partially achieved our goal. In the figure 2 we compare our subtraction terms with an exact integrand function (3.4) numerically. From the plot 2(a) one can see, that in the limit |k T 2 | |k T 1 | (or |k T 1 | |k T 2 |), both MRK and MMRK subtraction terms approximate an exact integrand very well, except from the region of final-state collinear singularity, which is located at ∆φ 2 12 + ∆y 2 12 1 ( figure 2(b)). When all transverse momenta are of the same order, the MRK subtraction term significantly overshoots an exact integrand outside the Regge limitsẑ → 0 andẑ → 1 (plots (d) and (e) in the figure 2), while MMRK expression (3.13) gives a more reasonable approximation in a whole range ofẑ. The same behavior is observed in the DGLAP (figure 2(c)) and on-shell (figure 2(f)) limits. In general, the MMRK subtraction term is smaller than standard MRK subtraction in whole phase-JHEP08(2020)055 space, which solves the problem of large negative NLO corrections, typical for BFKL-type calculations, as we will see in section 6.
Problems which we encounter here are familiar to the practitioners of NLO calculations in High-Energy QCD. For example, the similar severe over-subtraction problem was observed in the calculation of p T -spectrum of leading forward hadron in proton-nucleus collisions within Color-Glass-Condensate formalism at NLO [61][62][63] and was solved in ref. [64] by improvement of the kinematics of the subtraction term. The improvement of MRK approximation for tree-level MEs by propagator factors, as in eq. (3.13) is also not new. Such factors where first introduced in the High-Energy Jets (HEJ) approach [65,66]. Also, the new factor, introduced in eq. (3.13) strongly suppresses squared ME in the region which is completely removed in the kinematic constraint approach [52,53]. In ref. [53] it was shown that kinematic constraint approach correctly reproduces certain leading large-logarithmic terms of collinear origin which can be found in the NLO and NNLO expressions for BFKL kernel (the latter is not known in QCD but has been conjectured in N = 4 Supersymmetric Yang-Mils theory [67]). But e.g. the approach of ref. [25] resums the same series of collinear corrections by matching DGLAP and BFKL evolutions. So it seems, that kinematic constraint, MMRK-HEJ and direct resummation approaches are solving the same physical problems of BFKL evolution in a compatible way, but further investigations are needed to confirm this hypothesis.
Adding the same propagator-factor to real-emission term of evolution equation (2.30) one ends-up with: The MMRK evolution factor depends on an additional scale µ S , which characterizes the "hardness" of the next splitting in the evolution cascade and is an analog of factorization scale in CPM. Technically the scale µ S is needed to express kinematical difference between variablesẑ (3.2) and z n (2.18) and for DIS kinematics the optimal choice is µ 2 S = Q 2 + q 2 T 1 , while the choice (µ (n−1) S ) 2 = (q T + k T ) 2 in the integrand of eq. (3.15) is due to eq. (2.19). This equation have to be solved to obtain the UPDF in MMRK-approximation. It might be instructive first to study the "kinematic constraint" version of MMRK evolution equation, which is obtained from eq. (3.15) by the replacement: In the end of this section let us make a comment concerning the consistency between real and virtual parts of eq. (3.15). Our MMRK-approximation has altered only it's realemission part, leaving the cancellation of IR-divergences intact. The situation here is the JHEP08(2020)055 same as in the HEJ approach [65,66], where the standard LO gluon Regge trajectory is used together with modified real-emission amplitudes. However the well-known "bootstrap" property of BFKL-equation is lost in eq.  To make our NLO calculations more methodologically transparent, we decided to use a simple phase-space slicing method, similar to one proposed in ref. [68]. Our matrix element has non-integrable singularities in two non-overlapping phase-space regions: soft region, which we define by following cuts on dimensionless energies of gluons: where phase-space slicing parameter 0 < δ s 1. And (final-state) hard-collinear region, where: with 0 < δ c δ s and gluons 1 and 2 not satisfying condition (4.1). In terms of variables k T 1,2 andẑ, the soft condition for the first gluon has the form: where k T 2 q T 1 , 0 <ẑ <ẑ max withẑ max = δ s e 2Y H (where Y H is defined in eq. (2.23)), and for the second gluon it can be obtained using the substitution (3.9).
To facilitate the integration over hard-collinear region, we parametrize k T 1,2 in terms of q T 1 = k T 1 + k T 2 and new transverse vector ∆ as follows: which in particular allows one to conveniently express the invariant mass of the pair as: In terms of new variable, collinear condition (4.2) has the form: and requirement of both gluons to be non-soft translates into limits onẑ: JHEP08(2020)055 The soft limit of squared PRA amplitude can be computed using the usual eikonal Feynman rule for the emission of a soft gluon with four-momentum k from the hard gluon leg with momentum p: g s f abc p µ i /(kp i ). To take into account the presence of incoming Reggeized gluon R ± , an additional contribution, proportional to (−g s )f abc n µ ∓ /(k ∓ ) should be added to eikonal amplitude. Hence, for the case at hands, the soft limit is: which leads to the following reduced squared amplitude in the k 0 1 → 0 soft limit: (4.8) Eq. (4.8) has been cross-checked with an exact D−dimansional amplitude of the process (3.1) in the soft limit. The hard-collinear limit again can be obtained using the standard collinear factorization theorem for squared MEs [60], but this time for final-state singularity: with s given by eq. (4.5). Eq. (4.9) also has been verified against an exact squared PRA squared amplitude in D-dimensions. Substituting eq. (4.9) to eq. (3.3) with the parametrization (4.4) and cuts (4.6) and (4.7) one finds, that up to effects suppressed as O(δ c , δ s ) the UPDF can be taken out of ∆ andẑ integrals: with x 1 computed by eq. (2.15) and the following contribution to the NLO coefficient function: which can be straightforwardly integrated to give: where dimensionless couplingᾱ s has been introduced after eq. (2.25).
Integrating over the soft region (4.1) we can either add up contributions of k 0 1 → 0 and k 0 2 → 0 limits and divide the cross-section by 2!, or just integrate eq. (3.3) with reduced ME (4.8) over the region (4.3) and omit the Bose-symmetry factor 1/(2!):

JHEP08(2020)055
This integral can be calculated using the well-known Mellin-space representation of the θ-fucntion: x y γ 1 γ + δ , (4.12) together with the formula for two-dimensional Euclidean "bubble" integral with general indices: Also one can notice, that sinceẑ <ẑ max 1, the factor 1/(1 −ẑ) in eq. (4.11) can be omitted up to terms O(δ s ). Hence after expansion over δ s and we get: where ξ = e −2Y H . Finally, we have to take into account the double-counting subtraction with the evolution. It doesn't influence the collinear limit, since subtraction terms (3.10) or (3.13) are not singular in the region (4.2) and hence lead to O(δ c )-suppressed contributions. However, subtraction terms (3.10) or (3.13) have non-trivial soft limit. Theû-channel subtraction term, which is obtained from eq. (3.10) or (3.13) by the substitution (3.9), in the region (4.3) reduces to: where ξ µ = (µ 2 Y q 2 T 1 )/(Q 2 + q 2 T 1 ) 2 = e −2(Y H −Yµ) , see eq. (2.23). One should integrate this expression over region (4.3) and subtract the result from eq. (4.14). We have checked by explicit calculation, that the "propagator-factor" in eq. (3.13) makes no difference, up to O(δ s )-terms, so the double-counting subtraction in the soft limit turns out to be the same for MRK and MMRK approximations.
Due to a rapidity-ordering θ-function in eq. (4.15), we have to split the integration overẑ at a pointẑ m = δ s /(ξ µ + ξ), so that the subtraction term for the coefficient function takes the form:

JHEP08(2020)055
where Ω 2−2 = 2π 1− /Γ(1 − ). Calculating this integral, one obtains: Taking all the results of this section together we get: This expression has several important features. First, logarithms of δ s in the coefficient in front of 1/ , which are present in eqs. (4.10), (4.14) and (4.17) have canceled, giving IR-divergence a chance to cancel against the loop correction. The term log ξ µ / will also do so, as we will show in the next section. Second, subtraction (4.17) removed all terms proportional to log ξ = −2Y H , since this logarithm is resummed in the evolution, and only terms decreasing as e −2Y H are left. And third, if one takes the choice of rapidity scale (2.23), corresponding to ξ µ = 1, then all terms ∝ log ξ µ go away and one is left with: Let us discuss a bit the physical meaning of singularities, arising in the subtraction term (4.17). The integration region of eq. (4.16) is sketched in the figure 3 with the lines of constant rapidity of the first gluon and it's constant energy overlaid.
Going along the line of constant energy in the direction of decreasing rapidity y 1 (i.e. in a direction of a proton), one ends-up in the region of small k 2 T 1 . Hence, the 1/k 2 T 1 -singularity in eq. (4.15) is actually a rapidity divergence, corresponding to a fact, that probability of emitting a soft gluon is flat in rapidity. The actual soft divergence is located in a corner (ẑ → 0, k 2 T 1 → 0), where one arrives going along the line of constant rapidity in a direction of decreasing energy. These two divergences overlap in a corner (ẑ ẑ max , k 2 T 1 q 2 T 1 ), producing an 1/ 2 -pole in eq. (4.17). In the soft integral (4.14), the 2/ 2 -term has two sources. First mechanism is the same overlap of rapidity and soft divergences as in subtraction term, and second -the overlap of soft and final-state collinear divergences. The 1/ 2 -pole contribution from the second source, surviving after double-counting subtraction, will cancel against the loop correction.

Virtual correction and subtractions
We have computed the one-loop correction to an amplitude of the process (2.12) in ref. [44]. Apart from IR and UV divergences, which we regularize dimensionally, it contains rapidity divergence, physically corresponding to rapidity of a gluon in a loop going far negative. This divergence required additional regularization, which we preform by tilting the Wilson lines in the definition of EFT [40] from the light-cone, as was first proposed in [54][55][56]: n µ ± →ñ µ ± = n µ ± + r · n µ ∓ , where 0 < r 1 is the regularization parameter. As it was already noted in section 2 after eq. (2.25), such regularization roughly corresponds to a smooth cutoff for gluon rapidity at (− log r −1 )/2 and rapidity divergence manifests itself as log r-term, arising before one expands loop integrals in .

JHEP08(2020)055
In the ref. [44] we have also shown, that in the full amplitude, which includes oneloop corrections to both scattering vertices and t-channel Reggeized gluon propagator (see the right panel of figure 5 in ref. [44]), the log r-terms cancel and the (one-Reggeon exchange) EFT result precisely reproduces the (negative-signature part of) the dimensionallyregularized one-loop (O(g 3 s )) QCD amplitude of the process O(q) + g(P ) → g(k 2 , Y H ) + g(k 1 , y 1 ), (5.2) in the Regge limit.
Here we use this fact to derive the universal subtraction prescription for the virtual rapidity divergence in the one-loop coefficient function (5.1), consistent with MRK (2.30) and MMRK (3.15) evolution equations. To this end we start with the interference of one-loop and tree-level corrections to the subprocess (5.2) which contributes to the CPM coefficient-function of the process (2.1) in the NNLO(O(α 2 s )). The EFT predicts leading power Regge (Y H − y 1 1 or z → 0) limit of this interference to be proportional to the squared tree-level matrix element of the subprocess (5.2) with the one-loop coefficient: where H (NLO), g virt. unsubtr. is the interference of one-loop corrected and tree-level scattering vertices g(P ) → R + (q 1 ) + g(k 2 ) in the EFT, and Π (1) is the one-loop correction to the Reggeon propagator (2.25). We stress again, that all log r-divergences cancel in eq. (5.3).
On the other hand, eq. (5.3) already contains large-logarithmic contribution which is reproduced by one iteration of virtual part of evolution equation (2.30) or (3.15) (see e.g. eq. (B.1) in our appendix B): where Y (1,2) µ are rapidity scales for OR − g and gR + g scattering vertices, with the optimal choice Y (1) µ → y 1 , Y (2) µ → Y H and ω g is a gluon Regge trajectory (2.26). Subtracting the latter evolution contribution from eq. (5.3) and rearranging the terms, one obtains the following expressions for subtracted one-loop corrections to both scattering vertices: which are also free from log r-divergences.
By similar reasoning, one can obtain the subtracted one-loop correction to the scattering vertex with any regularization for rapidity divergences, including one proposed in JHEP08(2020)055 ref. [69], which opens-up a possibility to automatize the NLO calculations in a variety of small-x physics frameworks. The non rapidity-divergent-part of eq. (2.25) depends on a chosen rapidity regulator, but the subtracted results (5.4) and (5.5) should be regularizationindependent.
Using eq. (5.4) and subtracting the known (see ref. [42] and references therein) counterterm for the UV-renormalization of the operator (2.2) in the M S-scheme: 2δZ / we obtain the following subtracted one-loop coefficient function: which we have rewritten in terms of log ξ µ = 2(Y µ − Y H ). The divergence structure in eq. (5.6) precisely matches that of eq. (4.18), so the final answer for analytic part of our NLO correction is: The remaining µ-dependence in the coefficient cancels against the running of a coupling λ of O(x) to an external source. The ξ µ -dependence should cancel with the UPDF evolution, and with the rapidity-scale choice (2.23) one obtains: For further reference we also write-down the δ s,c -independent part of q 2 T 1 Q 2asymptotics of eq. (5.8): which contains very strong negative doubly-logarithmic contribution at small q T 1 . We will discuss numerical implications of this in the next section.

JHEP08(2020)055 6 Numerical results
The analytic part of NLO correction, obtained in the previous section, should be added to the numerical integral (3.3) evaluated in D = 4 space-time dimensions over the region of phase-space where neither condition (4.1) nor condition (4.2) is satisfied and with subtractions (3.10) (MRK) or (3.13) (MMRK) included at the integrand level. Then, for sufficiently small values of δ c δ s 1, the dependence on this parameters is guaranteed to cancel. In the present section we will show some results of exploratory numerical calculations, performed with UPDF described in the appendix A. Throughout this section we use the scale-choice µ F = µ = Q. Our numerical calculations have been performed with the help of parallel version of the well-known VEGAS adaptive Monte-Carlo integration algorithm, implemented in the CUBA library [70]. The main purpose of this section is to show, that MMRK subtraction term (3.13) indeed leads to NLO correction with more reasonable physical behavior than MRK subtraction term (3.10).
The cancellation of δ s and δ c dependence is demonstrated in the figure 4. One can see, that the result for NLO correction is independent on δ c within integration accuracy practically for all points and the plateau in δ s is reached rather quickly both for MRK and MMRK subtraction terms. To obtain reliable numerical results at smallest values of δ s one have to use quadruple-precision arithmetic in the squared PRA ME subroutine, however such calculations are rather computationally-costly. In the calculations below, we will use fixed values of δ s = 2 × 10 −4 and δ c = 10 −3 δ s with which double-precision calculation is sufficiently stable and accuracy of the NLO cross-section better than 1% is reached.
Next, let us examine the relative size of NLO correction to inclusive structure function with MRK and MMRK subtractions at various x B and Q (figure 5). The NLO CPM curve is also shown in the figure 5 and NLO CPM K factor increases towards small x B and reaches up to a factor of two for Q = 10 GeV. The NNLO correction in CPM is negative at small x B and subtracts almost entire NLO correction, while N 3 LO correction [42] is positive again and has the same order of magnitude. This kind of behaviour demonstrates lack of stability of CPM calculation at small x B . The LO PRA curve in figure 5 grows towards small x B with the same rate as NLO of CPM, calculated with the same collinear PDF, which was used to generate the doubly-logarithmic UPDF. However, this kind of rapid growth is probably an artifact of doubly-logarithmic approximation, which produces UPDF with nonphysically hard high-q T tail. With the solution of the full MMRK evolution equation (3.15) the small-x B growth of the structure-function will be significantly sloweddown. However the detailed q T -shape of UPDF does not really affect the relative size of NLO correction. We also do not show the scale-variation bands, since doubly-logarithmic UPDF does not depend on scale µ Y and corresponding logarithms in the NLO correction will have nothing to cancel against. The more detailed numerical study will be performed once the solution of eq. (3.15) will become available.
For both values of Q = 10 and 50 GeV, the NLO PRA correction with MRK subtraction is negative and for Q = 50 GeV it becomes larger than LO PRA term at x B > 10 −4 , demonstrating a severe perturbative instability of the calculation with MRK-subtraction. On the contrary, MMRK-subtracted NLO results look reasonable for Q = 10 GeV over JHEP08(2020)055 the whole range of x B < 0.1 and for Q = 50 GeV, the NLO correction becomes larger than 50% only at x > 10 −2 . The K-factor of the MMRK-caclculation flattens-out towards small values of x B as one would expect to see for log(1/z)-resummed calculation, unlike the K-factor of CPM calculation.
It is interesting to investigate the reason, why NLO correction to the inclusive structure function for moderate values of x B > 10 −2 becomes larger at higher scales. In the left panel of the figure 6 we plot the inclusive structure function differential w.r.t. transverse momentum of an incoming Reggeon q T 1 at LO and NLO of PRA. At NLO we find a large negative contribution to the cross-section at small values of |q T 1 | Q. This contribution originates from the doubly-logarithmic term ∼ − log 2 (Q 2 /q 2 T 1 ) in eqs. (5.8) and (5.9), as it is evident from comparison of the full NLO PRA results with approximate result obtained with an asymptotic coefficient function (5.9) (dash-dotted histogram in the left panel of the figure 6). These doubly-logarithmic effects arise due to exchange of collinear virtual gluons in the PRA coefficient fucntion and can be factorized from it into a universal Sudakovtype doubly-logarithmic formfactor which will become a part of UPDF. An interesting feature of our formalism is, that this kind of doubly-logarithmic effects is entirely due to virtual exchanges and all doubly-logarithmic effects due to real emissions already had been factorized-out. We will discuss this problem in more detail elsewhere.
In fact, most of the cross-section is accumulated at moderate values of 1 GeV < |q T 1 | < Q, where doubly-logarithmic effects are much less important than single-logarithmic effects which we discuss in the present paper. With x B -decreasing, contribution of moderate and large |q T 1 | grows (see figure 8 in the appendix A), which explains, why NLO PRA K-factor flattens-out towards small x B and becomes much less sensitive to the value of Q 2 .

JHEP08(2020)055
Single-logarithmic effects ∼ log(Q 2 /q 2 T 1 ) which drastically improve the quality of MMRK-approximation in comparison with the MRK calculation, come from the DGLAP region of NLO real-emission phase-space q 2 In CPM this is a region of initial-state collinear divergence, which is subtracted from NLO correction and governs PDF evolution. In PRA, there is no initial-state collinear divergence in the coefficient-function calculation. All collinear divergences are subtracted at the level of UPDF (see appendix A). But the DGLAP region still generates large contribution, enhanced by log(Q 2 /q 2 T 1 ) and proportional to the DGLAP splitting function p gg (ẑ), see eq. (3.11). Double-counting subtraction term partially subtracts this contribution, but for MRK-approximation, the quality of this subtraction is very poor, resulting in a huge negative NLO correction. MMRK subtraction term approximates an exact DGLAP splitting function better, as it was discussed in section 3, but the remaining mismatch still generates some "collinear" logarithmic term. However, at smaller values of x B and Q this logarithmic term does not present such a big problem as for x B closer to one and at higher Q Λ, where Λ ∼ 1 GeV is a scale of nonperturbative transverse momentum, which is present in UPDF. This factor also contributes to better stability of NLO correction at Q = 10 GeV vs. 50 GeV in the figure 5.
An interesting feature of PRA is, that many observables related with transverse momentum are available already in the LO. For the process at hands, such an observable is a leading jet p T -spectrum, which at LO is just a structure function (2.14), differential in transverse momentum q T 1 and jet rapidity Y H (2.23). Therefore now we have an opportunity to study the relative size of NLO correction in PRA also for the jet-p T spectrum.
To meaningfully define the jet observable at NLO we need to take into account, that UPDF evolution is not ordered in p T , so it can produce jets with transverse momenta higher than p T of a jet originating from the hard process. But evolution is ordered in rapidity, so we can avoid the need of fully-exclusive Monte-Carlo simulation by reasonably defining the "most forward" high-p T jet. We do this as follows: 1. If both gluons in the NLO subprocess (3.1) have ∆y 2 1,2 + ∆φ 2 1,2 < R 2 with jet-radius parameter R = 0.4 in our numerical calculations below, their four-momenta are added to form a four-momentum of a jet.
2. Otherwise the four-momentum of a gluon leading in p T and lying within rapidityacceptance is taken as a jet four-momentum 3. If rapidity of the gluon subleading in p T is y subl. < y jet , then it is unconstrained 4. If y subl. > y jet , we reqire it's p In all other respects, our NLO calculation for jet-p T spectrum proceeds the same way as for inclusive structure function, with no need to re-calculate the analytic part, just the phase-space slicing parameters have to be taken sufficiently small to avoid interference with jet definition.
Numerical results for jet-p T spectrum are shown in the figure 7 for two different values of Q = 10 and 50 GeV and the same value of "center-of-mass energy" S = Q 2 /x B . For jet p T -spectrum we use the same factorization and renormalization scale-choice as for JHEP08(2020)055 dF/dp  inclusive SF. Here we find larger NLO corrections at smaller scales, which most likely reflects a steeper decrease of PDF with increasing values of x at smaller scales. The MMRK approximation again leads to smaller NLO correction and at Q = 50 GeV, the NLO correction to jet-p T spectrum is negligible in most bins, suggesting that LO PRA is a good approximation for this observable in this kinematic region. Finally, the right panel of the figure 6 allows us to examine, how NLO correction to the jet p T -spectrum is distributed w.r.t. transverse momentum of incoming Reggeon q 2 T 1 . For the jet-p T spectrum, the loop correction and IR-cancellation effects are present only in two bins in the right panel of figure 6 where the LO term is nonzero, and one can see, that NLO result in this bins is very close to LO, so NLO correction in those bins is negligible. It is the behavior of subtraction term, which is responsible for the bulk of NLO correction to the jet-p T spectrum.
We emphasise, that doubly-logarithmic effect ∼ log 2 (Q 2 /q 2 T 1 ), which was de-stabilizing the NLO correction to the inclusive structure function is absent for the jet p T -spectrum with p jet T ∼ Q so the NLO-correction to this observable is particularly stable and suitable for calculations in k T -factorization and PRA. This conclusion is likely to generalize on a wide class of observables which are sensitive to the transverse momentum of incoming partons and absent in the LO of CPM. The ∆φ-spectrum of a dijet or multi-jet system at small values of ∆φ (away from back-to-back configuration) is a good example of such an observable [4][5][6].
As a conclusion we emphasize, that in the present paper we have formulated the technique of NLO calculations in PRA for gluon-induced processes and the MMRKapproximation for squared matrix element in QCD with emission of additional partons. This approximation should be used consistently as the subtraction term in NLO correction and in the UPDF evolution, leading to improved perturbative stability of the calculation both for inclusive structure function and jet cross-section.

JHEP08(2020)055
be subtracted from it. The renormalization factor which subtracts collinear divergences from the hard process can be defined in N -space as (see eq. (2.28) in ref. [3]): where S = exp [ (log 4π − γ E )] is the usual factor defining the M S-scheme, µ F is the factorization scale and γ k (N ) are the coefficients of expansion of DGLAP anomalous dimension γ gg (N, α s ) in powers ofα s . In agreement with known results [3,71], we find that the following series of coefficients leads to subtraction of collinear divergences from C(N, x T , µ) up to O(α 9 s ): We have checked up to O(α 9 s ), that the finite part of the coefficient function can be represented as: wherex T = x T e γ E /2, the non-cancellation of γ E and 1/2 is a consequence of working in x T -space. The exponential factor in last equation resums double-logarithms of the form α s log(x 2 T µ 2 F ) log(1/x) and corrections to it are at best -single-logarithmic and start at O(α 3 s ). Therefore, the double-logarithmic approximation is the basic approximation for UPDF. Converting the exponential factor back to x-space one obtains: where I 1 (x) is the Bessel function of the first kind. Before Fourier-transforming this expression numerically back to q T -space, we multiply it by the non-perturbative shape-function, which we take in a Gaussian form, suppressing large values of x T : where parameter Λ, equal to 1 GeV in our numerical calculations, characterizes the spread of "intrinsic" transverse-momentum of a gluon in a proton. Finally, to obtain the UPDF we take a Mellin convolution of the evolution factor with the collinear PDF as in eq. (2.9).
In the numerical calculations of the present paper we have used the NLO set of CTEQ-14 PDFs [72] as a collinear input and the NLO running of α s corresponding to this PDF set with α s (M Z ) = 0.106, as provided by LHAPDF library [73].
In the figure 8 we compare Doubly-Logarithmic UPDF with several other phenomenological UPDFs known in the literature. All UPDFs in this figure are obtained from HERAPDF-20-NLO-EIG PDF set [73,74] as a collinear input. In general, one observes that all distributions are broadening in q 2 T with decreasing x. There is a serious disagreement between: different approaches at small values of |q T | ≤ 1 GeV, since no dedicated fits where performed there. At moderate 1 GeV < |q T | < µ all approaches basically agree in shape, which is a manifestation of universality of doubly-logarithmic approximation. For |q T | > µ Doubly-logarithmic and Collins-Ellis-Blümlein UPDFs continue the same trend, while KMRW and PB UPDFs demonstrate more physical softer behavior. In general, the latter UPDFs tend to agree with data better, while Collins-Ellis-Blümlein tends to over-estimate absolute values of cross-sections. Agreement in shape and normalization of high-q 2 T tail between Collins-Ellis-Blümlein UPDF and our doubly-logarithmic approach is not surprising, since both approaches resum the same tower of double logarithms log(µ 2 /q 2 T ) log(1/x), but Collins-Ellis-Blümlein approach is set-up directly in q Tspace and is less systematic from the point of view of keeping correct M S-scheme definition. Also, Collins-Ellis-Blümlein-approach does not contain any non-perturbative shapefunction, which explains different small-q T behavior. JHEP08(2020)055 emission, and |k T + q T | |q T |, where the first emission is collinear, while the second emission generates transverse momentum. The 1/ pole of the infra-red divergent contribution from the first region: is cancelled by the virtual contribution (B.1), while the divergence from collinear region -remains. It is natural to expect the collinear divergence to appear, since we have two real emissions and one of them can take whole transverse momentum q T , while another generates collinear divergence which should be absorbed into collinear PDF.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.