A new subtraction scheme at NLO exploiting the privilege of k T -factorization

We present a subtraction method for the calculation of real-radiation integrals at NLO in hybrid k T -factorization. The main difference with existing methods for collinear factorization is that we subtract the momentum recoil, occurring due to the mapping from an ( n + 1 ) -particle phase space to an n -particle phase space, from the initial-state momenta, instead of distributing it over the final-state momenta.


Introduction
In the calculation of real radiation contributions to scattering cross sections in perturbative quantum chromodynamics (QCD) one encounters divergent phase space integrals (cf.[1] for an accessible introduction).On one hand, these are preferably dealt with using dimensional regularization, in order to match the divergences from virtual contributions and have a manifest cancellation of those.On the other hand, one would like to calculate arbitrary differential distributions for arbitrary fiducial phase space regions, and preferably even generate event files, pointing at the Monte Carlo method of numerical integration as the obvious method of choice.One way to satisfy both requirements is the subtraction method, in which the divergent integrals are regularized by subtracting extra terms that point-wise cancel the singularities causing the divergences.These extra terms are then constructed such that they themselves can be integrated analytically within dimensional regularization, explicitly exposing the divergences ready to be cancelled against the virtual ones.While conceptually easy, explicit implementation faces various complications especially for increasing final-state multiplicity, and over the last decades several solutions have been devised both at next-to-leading order (NLO) [1][2][3][4][5][6][7] and next-to-nextto-leading order (NNLO) [8][9][10][11][12].
All the cited approaches were designed for collinear factorization, for which the partonic initial-state momenta are light-like.In this work, we consider k T -factorization [13,14], hybrid k T -factorization to be precise for which only one of the initial-state partons is space-like, and its accompanying PDF depends on transverse momentum (k T ) besides the longitudinal momentum fraction (x) [15][16][17].This factorization becomes phenomenologically relevant for scattering events in which the final-state products are boosted towards one initial-state direction, implying that one of the collinear momentum fractions was much smaller than the other.It has already been applied at NLO precision for selected processes [18][19][20][21][22][23][24].The matrix elements with spacelike initial-state partons cannot be obtained from a trivial kinematic continuation of completely on-shell matrix elements, and need to be defined properly.One approach to achieve this is the auxiliary parton method [25], which is particularly well-suited to be applied together with efficient and easily automatable methods for the calculation of matrix elements that are successfully applied in programs for collinear factorization.The approach is implemented as such at tree-level in the parton-level Monte Carlo generator KATIE [26].The auxiliary parton method embeds the process with space-like partons into a process in which the space-like parton is replaced by a pair of on-shell partons.One can extract the necessary extra Feynman rules to construct the desired space-like amplitudes directly as is done in KATIE, or one can use a complete expression of an embedding on-shell amplitude to extract the desired space-like amplitude.The latter method was applied in [27,28], to calculates space-like one-loop amplitudes for selected processes.
While most divergences appearing in collinear factorization eventually cancel at the partonlevel, some divergences that can be identified as being of the collinear kind remain, and are taken care of by renormalization of the PDFs.The fact that this works in a process-independent manner and to any order in the perturbative expansion is an essential ingredient in the predictive power of collinear factorization.In [29], a scheme was outlined to perform calculations in hybrid k T -factorization with the auxiliary parton method at NLO.It was shown that, while there are more divergences that do not cancel at the parton-level, their occurrence is process independent, opening the possibility of a renormalization procedure and perform predictions.
It is also shown in [29] that the process-dependent part of the real-radiation contribution to the cross section is obtained just as in collinear factorization, as the phase space integrals of the processes with one more final-state parton and a jet algorithm allowing for one final-state parton to be unresolved.It was argued that all the divergences coming from these integrals are exactly as in collinear factorization, except for the remnant collinear divergence related to the space-like initial-state parton.In this work, we will make this more explicit, by introducing a subtraction method to tackle the real-radiation phase space integrals, and providing the universal expressions of the divergent part of the integrated subtraction terms.
We follow the work of [1,4] as a guideline for the organization of the subtraction terms.The main difference is the treatment of the momentum recoil appearing in the necessary mapping of the (n + 1)-particle phase space to a n-particle phase space.Since we have 4 independent kinematic variables in the initial state (2 longitudinal components, and 2 transverse ones), we can simply subtract the recoil from those.This has the advantage that the final-state momenta are touched in a minimal way, which is helpful in the calculation of differential distributions, possibly minimizing the so-called mis-binning problem.As a drawback, the finite part of the integrated subtraction terms cannot be obtained completely analytically and requires the application of numerical integration.This, however, can easily be combined with the Monte Carlo procedure for the phase space integration.

Exposition of the subtraction scheme
We consider hadron scattering with light-like hadron momenta P, P. For arbitrary momentum K, we can then make the Sudakov decomposition with We use the labels χ, χ to indicate the partonic initial states, and the momenta in hybrid k Tfactorization are that is we omit the labels for the longitudinal and transverse components of those: x χ = x and xχ = x.Eventually, we will assume that calculations are performed in a frame for which P µ = (E, 0, 0, E) , Pµ = ( Ē, 0, 0, − Ē) , so S = 4E Ē .
The symbol E with an appropriate subscript will also be used for the energy of momenta in general.Also, we will employ unit vectors denoted with the letter n.So for a light-like momentum K µ we may write We will use a separate symbol for the sum of the parton-level initial-state momenta: Integration over the initial-state variables will be abbreviated as One thing we need to address already before we continue is our notation regarding dimensional regularization.We write and will often use ε in calculations because it can be imagined to be positive for infra-red divergent integrals.Furthermore, we define and The common overall factor 4πα s µ 2ϵ for singular limits of matrix elements divided by the phase space factor (2π) 3+ ε of the radiation can be written as The denominator factor π ϵ typically drops out in final integrated expressions.The Born-level formula for the cross section in hybrid k T -factorization is where dΦ is the differential phase space for the final-state momenta {p} n dΦ Q; and M 2 Q; {p} n is the tree-level matrix element of the considered parton-level process, summed over helicity-and color-degrees of freedom.It does not include symmetry factors nor factors to average over initial-state degrees of freedom -they are captured by S n .The function L includes both the collinear and the k T -dependent PDF, and also the flux factor where µ F is factorization scale.Like in [29], we will restrict ourselves to the situation in which the χ parton is a space-like gluon.The light-like χ parton can be gluon or quark or antiquark.The symbol J B denotes the jet function, which at Born-level eventually consists of a number of phase space cuts making sure the singularities of M 2 are avoided.

Limits of matrix elements
In the following, we assume that all final-state parton momenta are light-like.The real-radiation contribution at NLO involves scattering processes with one more final-state parton than at Born level.Furthermore, the jet function, now denoted J R , does not avoid all singularities of the treelevel squared matrix element anymore, but allows one pair of partons to become collinear, or one parton to become soft The jet function behaves in those limits such that where {p} r / n is obtained from {p} n+1 by removing momentum p r , with Of course, at the collinear limit p r ∥ p i , we have (1 + z ri )p i = p r + p i .We silently assume here again that final-state parton momenta are light-like.Let us consider the matrix element for collinear factorization, and imagine k ⊥ = 0.It is singular in these limits, but the singular behavior factorizes in a universal, and well-known, way.We want to write them in a factorized form and level of detail necessary for the further discussion.Also, we would like to write them in a functional sense, such that the symbols representing momenta are identical on the left-hand side and the right-hand side.We have The quantities R represent the singular behavior as function of the radiative momentum p r .The quantities Â involve squares of tree-level scattering amplitudes with n final-state momenta, but including color and spin correlations.The mappings T establish the aforementioned demand that the relations hold in the functional sense, and make sure that on-shellness and momentum conservation hold for the arguments of Â, while the symbols represent the same momenta on the left-hand side and the right-hand side of the relations.The sets of final-state momenta {p} r / n and {p} r /;i n do not satisfy momentum conservation, that is their sum is not equal momentum before the semicolon, and the mappings T fix that by deforming them a bit, keeping them however on-mass shell e.g.: They could for example be the phase space mappings employed in the well-known subtraction schemes to make sure that the subtraction terms are valid at least in a non-vanishing phase space region around the singularities.The limits for the matrix element are given in more detail by Here, M 2 color(a,b) refers to the color-correlated squared matrix element for the process with one final-state gluon fewer than the original process, see Appendix B. The labels 'ir' and 'χr' for the other matrix elements indicate that it can be for a process with something less trivial than one final-state gluon fewer, like a quark-antiquark pair replaced by a gluon, an initial-state gluon replaced by quark and a final-state antiquark removed, etc.The collinear splitting functions Q are defined in Appendix A as functions of the ratio of the momentum fractions of the momenta becoming collinear.It is understood that their subscripts refer to the flavors of the indicated particles in the process.The minus-sign for the initial-state case stems from the fact that we define the splitting functions for matrix elements with all momenta outgoing, implying that p i •p r is negative if p i is an initial-state momentum.The symbol ⊗ for the collinear cases indicates that there are spin-correlated matrix elements involved, see Appendix B. For the initial-state collinear limit, we prefer to write the expression in terms of the variables x r , xr .The χ-case has the same form with the bar-and non-bar quantities exchanged.The a-sum and b-sum for the soft limit are over all external partons, both initial-state and final-state.
The singular behavior causes the real radiation integral, with M 2 Q; {p} n+1 and J R {p} n+1 , to be divergent, and the task at hand is to split it into a finite part and a divergent part that can be explicitly expressed as a Laurent expansion in ϵ within dimensional regularization This decomposition is not unique, because σ div R (ϵ) will typically include finite pieces of O(ϵ 0 ).Only the 1/ϵ 2 and 1/ϵ terms in σ div R (ϵ) are universal.The formulas for the singular limits in the foregoing are for the case when both initial-state momenta are light-like.It turns out that if we assume the χ-parton to be space-like, then the formulas for the limits are just the same, only with space-like matrix elements.Despite nonvanishing k ⊥ , there is still a collinear singularity when a final-state gluon becomes collinear to the hadron momentum P associated with the χ-gluon.It was shown in [29] that the formulas for the singular limits have forms identical to the on-shell case, just with a particularly simple splitting function In fact, there are no spin correlations involved and the symbol ⊗ is just multiplication in this case.With this work we numerically confirm the correctness of this formula by observing that our subtraction terms indeed cure the singularity.The sum of the n final-state momenta must be equal to the total initial-state momentum, the quantities before the semicolon, in order for the matrix elements to be well-defined.The mappings T before involved only final-state momenta, as is the case in most subtraction schemes, but this is not apriori necessary, and the initial-state momenta could be involved in the re-distribution of the recoil too.For collinear factorization, the initial-state momenta can only accomodate longitudinal momentum components.For hybrid factorization, they can accomodate a full 4momentum, and we can simply re-distribute the recoil over the initial-state momenta as Because the expressions on the right-hand-side are well-defined all over phase space, and not only at the limit, they can be used to construct subtraction terms in order establish Eq. ( 31).
The formulas above cannot be the only ingredient for the subtraction terms, because of double subtraction of soft-collinear singularities.We also want to augment them with phase space restrictions limiting their support around the singularities they are supposed to subtract.Finally, it will prove to be beneficial for soft terms with final-state collinear singularities to apply the final-state collinear type of phase space mapping.Fact is, however, that subtraction terms can be cast into the two argument types Q − p r + z ri p i ; {p} r /;i n and Q − p r ; {p} r / n .We introduce new objects R and A, now without hat, that will be precisely defined later, but fit the aforementioned forms.

Organization of the subtraction terms
We define the finite radiative integral as where the r-sum is over all final-state partons, and with where also the i-sum is over all final-state partons with R F rr (p r ) ≡ 0. We make sure to avoid double counting in the precise definition of the terms.The terms with label 'F' refer to finalstate singularities and the ones with 'I' to initial-state singularities.Some of the terms carry the label r because they depend on the flavor of final-state r.The first three lines categorize the three different types of momentum recoils that appear in our scheme.The fourth and third are identical with respect to the recoil, but we decided to categorize them separately already for later convenience.
We subtract a momentum q, with subscripts depending on the subtraction term, from the arguments of the L-function.This is allowed as long as this momentum vanishes in the limit for which the respective subtraction term is supposed to cure a singularity.In our situation of k T -factorization, this momentum can have both non-vanishing longitudinal and transverse components.For the initial-state soft case, this could for example be p r itself.For the initialstate collinear cases, this momentum cannot be equal to p r , since it does not vanish in those collinear limits.Then, it can be "at most" q r,χ = xr P + k ⊥ and q r,χ = x r P + k ⊥ .For the other cases, the q-momentum can be identical to the recoil subtracted from the initial-state momenta in the matrix element.
Using the result of Appendix D, we find that the subtraction terms integrate to , where with and From these formulas it appears that the best choice for the q-momenta is it to be equal to the recoil, so the L-function does not depend on the integration variables of p r anymore.As mentioned before, this is not possible for the initial-state collinear cases, and we cannot avoid arguments n .This of course also happens in the well-known subtraction formulas for collinear factorization, leading to what in [3] is called the "P-operator".
While the ϵ −2 -and cancelling ϵ −1 -parts can be obtained analytically by applying a subtraction approach on these integrals again, the finite remainder of the integrals cannot be calculated completely analytically.This is relevant also when we choose the q-momentum to be equal to the recoil, seemingly removing the L-function from the integral, because the extraction of the ϵ poles will typically not follow the integration borders dictated by Θ(q) in order to be performed analytically.
Another issue is that, while formally correct, subtracting the recoil from Q inside the L function may not work very well in practice.The reason is that, although the recoil may be small, it is not guaranteed that its transverse components are small compared to the initial-state k ⊥ .The longitudinal variables x, x are upper limits on the size of the longitudinal components of the recoil, so the latter can indeed be considered small, but for the transverse components this is not the case.
Whether this becomes a problem depends on the actual behavior of the k ⊥ -dependent PDF.If it does, then the solution is to not subtract the transverse components of the recoil from Q inside the L function, with the consequence for all integrated subtraction terms the finite part requires numerical integration with the L function as part of the integrand.Realize, however, that these numerical integrals can be performed in a Monte Carlo approach, by including them into the phase space integral.
In our actual numerical studies, we used the parton-branching k T -dependent PDFs [30,31], and the aforementioned problem did not occur.Still, in the following, we will work out the integrals for the safe, "conservative", approach.We will express them in terms of divergent integrals that are calculated analytically, and are labelled 'div', and finite integrals labelled 'fin' that must be integrated numerically.In the "progressive" approach, the non-trivial integrands reduce to theta-functions only enforcing integration limits.More precisely, the difference between the two approaches is expressed by the function that will appear in the integrals of the subtraction terms.

Definition of the subtraction terms
Now we present the actual subtraction terms more precisely.These follow the same structure as in [1,4], but with different phase space mappings, and organized into different groups.It turns out to be beneficial to apply the well-known technique of splitting the soft factor into two terms each with only one collinear singularity For each of these with final state a = i, we can then apply the final-state collinear type of mapping, adding the energy of the radiation to the radiator.We also mention that for the finalstate collinear terms, we avoid double counting with the help of a selector function θ(E r < E i ) demanding that a radiator has larger energy than the radiation.

Final-state terms
Firstly, there are the terms in Eq. ( 37) with arguments Q − p r + z ri p i ; {p} r /;i n .While they do involve initial-state "spectators" in the soft terms, we refer to them collectively as "final-state terms".They are given by with All these matrix elements carry the aforementioned arguments.The b-sum is over all partons, initial-state and final-state (and b ̸ = i, but then again n i • n i = 0).The third line contains the soft-collinear counter term correcting the double subtraction of soft-collinear singularities.We use the notation and later we will also use if exactly one of i, j refers to a gluon, T R if none of i, j refers to a gluon. (51) Color conservation implies The soft and soft-collinear terms vanish if r does not refer to a gluon.The restriction θ(E r < E 0 ) appears both for the soft and the soft-collinear term, so these clearly match in the collinear limit n r → n i , undoing the soft-collinear double subtraction.In the soft limit E r → 0, the restriction θ(E r < E i ) = 1, so then the collinear and the soft-collinear term match.The final result for σ R should not depend on the the parameters E 0 , ζ 0 which in practice is a powerful check of correctness.Furthermore, it opens the possibility to adjust the behavior of the integrand under numerical integration, and we find that it is improved by the restrictions.Also, for smaller values of the restricting parameters, fewer subtraction terms contribute to a given phase space point, making the evaluation per phase space computationally point cheaper.In particular, by choosing E 0 = E i (or min(E 0 , E i )), at most 1 of the 2 assignments "radiator" versus "radiation" for i, r in the soft terms contributes per phase space point.

Initial-state terms
The terms in Eq. ( 37) with arguments Q − p r ; {p} r / n have initial-state radiators (χ, χ) instead of the final-state i in the previous terms.We only present the χ-terms.The χ-terms are completely analogous with exchange of the bar-and non-bar variables.Due to universality they largely fit the same formulas as the final-state ones, but we find it more convenient to write the collinear terms using the variables xr = 2P•p r /S and x r = 2 P•p r /S: There is only a collinear singularity for the space-like χ case if the radiation is a gluon, so no subtraction term with a quark or antiquark as radiation is required.For the on-shell χ-terms, they are included.

Calculation of the integrated subtraction terms
In the following, we present the quantities L X,Y Z in the decomposition of the integrated subtraction terms of Eq. (38) as (57) Each of the L X,Y Z (ϵ) will then be decomposed into a divergent part that can be expanded into a Laurent series in ϵ analytically, and a finite part that must be integrated numerically Just to be clear, we mention that L Q; {p} n contains the PDFs for the original radiative process, which may indeed be the "wrong" ones for processes of M ar 2 , in which say an initial-state gluon may have been replaced by a quark.This is usually the case for any subtraction scheme, but not always so explicitly visible if the formulas are written only at the parton-level.

Final-state terms
We need to combine Eq. (40) and Eq. ( 46).We remind the reader of the functions Θ(q) in Eq. ( 43) and ℓ ⊥ (q) in Eq. ( 44), which are defined applying the decomposition of the momentum q in Eq. (1).

Collinear terms
Realize that the explicit p i and z ri in Eq. ( 47) are the ones before changeing integration variables from p i to pi = (1 + z ri )p i , so for the integrated subtraction terms we need to insert z ri = zri /(1 − zri ) and p i = (1 − zri ) pi , and then just pi rename with The function S F,col is singular both in the collinear limit ⃗ n r → ⃗ n i , and in the soft limit E r → 0 because of the 1/z ri hidden in P ir (1 − z ri ).In either limit, however, p r⊥ − z ri p i⊥ vanishes.Thus we see that the integral is finite, and can be calculated analytically.

Soft terms
Again because p r⊥ − z ri p i⊥ vanishes both in the soft and the collinear limit, the integral is finite, while can, in principle, be calculated analytically.We find it however too complicated.Instead, we introduce which vanishes in the soft limit, and becomes equal to E r in the collinear limit, so we can define which is easier to calculate analytically.

Soft-collinear terms
Analogously, we have with

Initial-state soft and soft-collinear terms
We need to combine Eq. ( 42) with Eq. ( 54) and Eq. ( 55).We will only handle the χ-terms.The χterms are completely analogous with exchange of the bar-and non-bar variables.The arguments x q = x r , xq = xr of Θ(q) in Eq. ( 43) now stay positive.Furthermore, we will pick this function apart into its two factors and write those explicitly instead of the symbol Θ.

Soft
Let us introduce the abbreviation We need to calculate with ℓ ⊥ (q) as in Eq. (44).Since p r⊥ vanishes both in the soft and the initial-state collinear limit, we see again that a single subtraction suffices and that we can define We gave the quantities an extra label since they are not the final ones.We find the integral for L I,soft,div,1 χb (ϵ) too complicated because it both depends on n b and includes the restriction θ(x r < 1 − x).We split it further into with The θ functions in L I,soft,fin,2 χb prevent p r from becoming soft, while they become identical in the collinear limit n r → n χ ⇔ xr → 0. Notice that L I,soft,div,2 χb (ϵ) does not depend on the spectator b, and will prove to combine nicely with a similar soft-collinear contribution.The final finite soft contribution we denote L I,soft,fin χb = L I,soft,fin,1 χb + L I,soft,fin,2 χb . (79)

Soft-collinear
For the soft-collinear case, we abbreviate We need to calculate and split it as L I,soco,div,1 with By explicit calculation we see that, up to the color factor, the 1/ϵ pole in L I,soco,div,2 χ (ϵ) is identical to the one in L I,soft,div,2 χ (ϵ), so with the help of color conservation it will cancel in the sum over spectators for the latter.We collect the remaining finite pieces in (88)

Initial-state collinear terms
We need to combine Eq. ( 41) and Eq. ( 53).Realize that the variables x, x in Eq. ( 53) are the ones before shifting the initial-state variables.We consider here the χ-case, so Q → Q + x r P + p r⊥ and we need to transform x → x + x r to get with and now The χ-case is completely analogous, and we keep Q χ general without any assumptions.It is useful to already write the integral over p r explicitly in terms of the Sudakov variables In principle we could define L I,col,div χr such that L I,col,fin χr vanishes in the progressive case.However, we prefer to keep L I,col,div χr simpler at the cost of non-vanishing L I,col,fin χr , by removing the condition xr < 1 − x for the former.Using The following variable substitutions now make sense leading to where we used the fact that −zQ χr (z − 1) = P χr (z) , see Appendix A .The possible flip gq ↔ qq for the χ case is understood as explained there.
A factor 1/z 2 appears instead of the 1/z one might have expected because L includes the flux factor.
The splitting function has a term 2C χr /(1 − z), causing a soft singularity.(For a χ-case without such a term we just set C χr = 0 in the following.)We can isolate the singularity via and where the plus-distribution only acts towards the right:

Final results for the divergent parts
We calculate the divergent contributions L X,Y,div Z (ϵ) in Appendix F (there in terms of ε = −2ϵ), and bring the expressions for L X,Y,fin Z to a suitable form for numerical integration in Appendix G. Here, we present the final results for the former.The final-state collinear ones are where The soft and soft-collinear contributions are L I,soco,div The initial-state collinear contributions are where P reg χr (z) is defined in Eq. (101).The splitting functions are evaluated at ϵ = 0, and P χr is the O(ϵ) coefficient.The results for the χ case are completely analogous with x ↔ x, P ↔ P, and E ↔ Ē.If the radiation is a quark or antiquark, we simply have C χr = 0 in the formula above.The function ℓ χ was defined in Eq. (92).
The expressions above belong to the integrals of the subtraction term with label r in Eq. ( 36).If r does not refer to a gluon, then the terms labelled 'soft' and 'soco' vanish.Let us consider the terms for which r refers to a gluon (r → g), and let us abbreviate Using color conservation, it is straighforeward to check that where Clearly, the divergences are independent of E 0 , ζ 0 , ξ 0 .The sum in Soft a (ϵ) is only over b.
Summing also over a and applying color conservation again, we see that a Soft a (ϵ) gives the well-known universal expression for the soft and soft-collinear divergences (cf.[1]): using we find The collinear left-over in Eq. ( 115) is also exactly as expected, and is the one that needs to be cancelled against a collinear counter term as prescribed by factorization.This is well-known for the on-shell initial state, the one with χ instead of χ, and in literature often identified as part of the "P-operator", but shown in [29] to be also present for the space-like initial state.There is a 1/z 2 instead of the usual 1/z because ℓ χ includes the flux factor.Except the initial-state collinear left-over, the divergences should cancel against those in the virtual contribution at NLO coming from the one-loop amplitudes.The virtual contribution involves the same phase space integral as the Born contributions, but has 2Re M † M one−loop (ϵ) instead of the Born-level matrix element |M| 2 .It was argued in [29] that the cancellation should happen exactly the same way for the case with a space-like initial-state as for the completely on-shell case.More precisely, the divergences in the one-loop amplitudes that do not fit in the universal formula (as given for example in [32]) where identified in [29], and confirmed in [28], and found to be universal themselves, so they can be treated within a renormalization prescription.In the following, we show that the divergences from the real radiation we found here indeed match the one-loop formula as in [32].
The formulas for the divergent contributions we found above are the result of concentrating on a given radiative process.To cancel all divergences, one must include all sub processes contributing to a given jet process.The formula for the total divergent soft contribution Eq. ( 119) however already matches the universal formula for the soft contribution in one-loop amplitudes.This comes out nicely, because removing a soft gluon from the radiative process exactly results in the Born process for which both the one-loop amplitude and the radiative process form the NLO correction.Including all collinear divergences for which at least one gluon is involved is easily achieved by adding the contribution when a radiative quark becomes collinear to a gluon from Eq. ( 106) to the γ (g) q of Eq. ( 117), leading to We see that also γ q now matches the expression for the one-loop amplitude.
In order to get the fermion-loop contribution for γ g included, one must go beyond concentrating on a single Born process.Then one sees that a radiative process for which final-state quarkantiquark pair becomes collinear, there is a collinear divergence of −T R /(3ϵ) (from Eq. ( 107)) to the Born process for which that pair is replaced by a gluon.The same for a final-state quark and initial-state quark.Then there is a factor 2 for the contribution for which the quark and anti-quark reverse roles, and a factor counting the number of quark families included, leading to γ g = 11  6 C g − 2 3 T R n f from the one-loop formula.

Numerical checks
We implemented the presented scheme using the framework of KATIE.In particular, we augmented the phase space generator KALEU ( [33]) with the channels matching the subtraction terms within the adaptive multi-channel method [34], in order to achieve convergent Monte Carlo estimates of the phase space integrals.In this section, we present some numerical results in support of the viability of the proposed scheme.We show that it indeed leads to convergent subtracted-real integrals, and that their combination with the O ϵ 0 part of the integrated subtraction terms is independent of the parameters E 0 and ζ 0 (Eq.( 47) and Eq. ( 48)).We link ξ 0 of Eq. ( 53) to ζ 0 via Eq.( 236).The phase space and renormalization/factorization scale are defined by where the sum is over all final-state partons.
In Fig. 1 we illustrate convergence and independence for a selection of partonic radiative processes relevant to dijet production.The parton branching PDFs PB-NLO-HERAI+II-2018-  The phase space and renormalization/factorization scale are defined in Eq. ( 121).set2 [30], taken from TMDLIB [31], were employd for the space-like side, and CT18nlo [35], taken from LHAPDF [36], for the on-shell side.While Fig. 1 present results for the full cross section for separate processes, Fig. 2 presents results for differential cross sections for a sum of processes, again relevant to dijet production.The results clearly indicate independence of the parameters E 0 , ζ 0 despite the rather substantial cancellation between the individual contributions.The computation were performed on a single modern laptop with 12 cores.

On-shell limit
The on-shell limit |k ⊥ | → 0 for tree-level matrix elements with a space-like gluon with momentum k µ χ (Eq.( 3)) is given by Here, we only make the dependence on k ⊥ explicit as argument of M, and denote by M ν (0) the on-shell amplitude with the polarization vector for the, now on-shell, gluon stripped off.The on-shell limit is "smooth" only including integration over the remaining azimuthal angle dependence, which is turned into a sum over helicities For tree-level cross sections, this limit is also smooth, and can be imagined to be established with the help of a special k ⊥ -dependent PDF that has a free parameter β such that Beyond tree-level, the smooth limits do not hold anymore, as explicated in [29].What we call real-radiation contribution in this write-up is not the whole NLO real contribution, but only the part that cannot be calculated analytically and that in [29] is referred to as the "familiar" real contribution.Still, one could imagine to construct an artificial real integral that would lead to the on-shell case, or one could try to construct space-like subtraction terms to be applied to the one-shell real integral.
Unfortunately, this cannot work, and what is worse, even in the k ⊥ -dependent case the subtraction method presented here strictly speaking fails at |k ⊥ | = 0.The reason is that at |k ⊥ | = 0 there is no point-wise cancellation between the radiative matrix element and the subtraction terms.Let us denote the radiative matrix element by M(k ⊥ , r ⊥ ) 2 , so with 2 arguments, where r ⊥ refers to the transverse part of the recoil of the radiation, and by M(k ⊥ − r ⊥ ) 2 , with only 1 argument, the matrix element appearing in a subtraction term.The latter has one final-state parton fewer, and the radiative recoil subtracted from the initial state.For the situation in which |k ⊥ | → 0 before |r ⊥ | → 0, we schematically have and we see that the two cases do not match.Only after integration over the remnant azimuthal angle dependence of both k ⊥ and r ⊥ , the cases match.So while the integrated subtraction terms correctly represent the divergences, at |k ⊥ | = 0 the terms fail at the task of point-wise cancellation of the singularities in the subtracted real integral.The reason why the subtraction method presented in this write-up still works is that the measure of phase space where this problem appears vanishes.Both space-like matrix elements and k ⊥ -dependent PDFs are in practice defined such that they are finite for |k ⊥ | → 0, while the 2-dimensional integration over k ⊥ provides a factor |k ⊥ |, so the differential cross section vanishes.To illustrate this, we present these differential cross sections for gg ⋆ → gg and gg ⋆ → ggg at tree level in the upper row of Fig. 3.We present results for both k T -dependent PDFs PB-NLO-HERAI+II-2018set2 and KS-2013-linear [37].While LO distribution of the latter does exhibit the, in this discussion dangerous, growth with k T → 0 associated with linear BFKL evolution, it eventually still vanishes.The lower row shows the O ϵ 0 -part of the real-radiation contribution to gg ⋆ → gg from the gg ⋆ → ggg channel, for both k T -dependent PDFs.We remind the reader that while the full real-radiation contribution must be positive, the O ϵ 0 -part on its own does not have to be positive.

Conclusion
We presented a subtraction scheme for the calculation of the real-radiation contribution in hybrid k T -factorization.It is different from existing schemes for collinear factorization in that the so-called momentum recoil is subtracted from the initial-state variables, leaving the final-state momenta untouched as much as possible.We calculated the divergent part of the integrated subtraction terms, and show that they have the expected universal form.The remaining finite part of these terms must be integrated numerically, and we presented the integrals in a form directly applicable for such a procedure.
We implemented the scheme and performed calculations for all processes relevant for 2jet production as NLO, and found that the subtracted real-radiation integrals indeed converge.Furthermore, we found that the sum of these integrals and the integrated subtraction terms is independent of the parameters restricting the phase space on which the subtraction terms are set not to vanish.Despite the fact that the phase space integration region for the integrated subtraction terms has 3 extra dimensions because of the finite remainders, their computation time to reach a given precision is still only a fraction of the computation time for the subtracted real radiation integral.

A Splitting functions
Consider a collinear limit of momenta with arbitrary real momentum fractions x, y where p µ and q µ are light-like with p • q ̸ = 0, while p • k T = q • k T = 0.The corresponding collinear splitting functions are defined as and are given by with where P ab are the usual collinear splitting functions, but labelled following the splitting products rather than the parent, and such that there is a soft pole at z = 0 if a = g, and at z = Again for z ∈ [0, 1] this leads to

B Correlated matrix elements
We present the correlated matrix elements that appear in the soft and collinear limits.We present the color-correlated ones for completeness, and the spin-correlated ones in a form that is closer to what one would like to use in the context of helicity amplitudes.Matrix elements involve a sum over color.We can write the dependence of an amplitude on the color indices explicitly as where ω l is the adjoint index a l if l refers to a gluon, ω l is the fundamental index i l if l refers to a quark, and ω l is the anti-fundamental index j l if l refers to an anti-quark.Implying the usual summation over repeated indices, the squared matrix element summed over color becomes The color-correlated squared sum where if l refers to a quark, −T c j l ȷl if l refers to an anti-quark. (144) With this definition, the equal-index correlator M 2 color(k,k) is obtained, by only replacing δ ωk ω k with T c ωk ω T c ωω k .For point-wise cancellation between collinear singularities in the real-radiation integral, spin correlations need to be included in the subtraction terms.That is, we need the limit of Eq. ( 129) before integration over the angular part of k T , but only at ϵ = 0.The spin-correlated matrix elements necessary for the 4-dimensional subtraction terms can be written in terms of helicities as follows.We can write the dependence of an amplitude on the helicity indices explicitly as The squared amplitude summed over helicities becomes The 1-helicity-fixed matrix element is obtained by replacing δ λl λ l → δ λl + δ λ l − .
Let the helicity amplitude be constructed with polarization vectors ε (l) λ l for gluon l.We denote

C Integration limits
For a function f(x, y) that has support 0 < x < 1 for the first argument, and vanishes outside this region we have To see this, firstly we observe that with z ri = E r /E i and with the implied integration limits on p r .First, we write where we introduced the notation in order to single out the momentum p i , and where {p} r /,i / n−1 is {p} n+1 with the momenta p r , p i removed.We can write the integration over the momenta p r and p i explicitly via

Now we change integration momentum
Notice that z ri p i = zri pi .
We get The integration restriction E r < E r +E i = Ẽi ⇔ zri < 1 is from now on understood to be implied by the Jacobian (1 − zri ).Now we rename pi Finally, we substitute Q ← Q + p r − z ri p i and use the result from Appendix C to find Eq.(154).This was the "collinear" case.The "soft" case is simpler, and we can see immediately that

E Radiative integral representations
We will encounter integrals of the type If f depends on a an external momentum p i , it is useful to use a frame with p i along the z-axis: with and where Rot i is the inverse of the rotation that rotates ⃗ n i to the z-axis: We also encounter the (4 + ε)-dimensional case, but with f only depending on p i via p i •p r : When there are two inner products p i •p r , p j •p r involved, it can be useful to apply the Sudakov decomposition p µ r = y i p µ i + y j p µ j + y µ T , with so If the integrand does not depend on the direction of y T , we can evaluate further to

F.1.3 Soft
Also here it is understood that E 0 is in fact the minimum of E 0 and E i .In the following, i refers to a final-state momentum, while j can both be initial-state and final-state.

F.2 Initial-state terms
In the following, i refers to an initial-state momentum, while j can both be initial-state and finalstate.

G Massaging the finite integrals
We bring the integrals to a form where f, g are finite and ρ 1 = 0, ρ 2 = 0 represent the soft and collinear singularities.One frequently applied ingredient will be the variable substitution when the function g found at some stage turns out to behave as G.1 Final-state terms

G.1.1 Collinear
We need to calculate and where Pir (z) = z P ir (1 − z) is finite at z = 0. Simply using a frame in which p i is along the z-axis, we can write and see explicitly that there is only one singular variable, and one subtraction is enough to regularize the integral.

G.1.2 Soft
We need to calculate with

G.2.3 Soft-collinear
Using Eq. ( 4), the restriction xr < ξ 0 x r can be cast in the form of the restriction on the variable ζ < ζ 0 in the lab frame via Also, we have We need to calculate (243)

Figure 1 :
Figure 1: Statistical distributions of the Monte Carlo estimate of the O ϵ 0 part of the realradiation contribution to the NLO cross section for dijet production, for different values of the parameters E 0 , ζ 0 .The essence of the plots is that the distributions have finite width, proving convergence, and that they overlap, proving independence of the parameters E 0 , ζ 0 .The values of the integrated subtraction terms (which are precise below the per mille level) are written separately, to indicate that the same final results indeed come from different individual terms.The phase space and renormalization/factorization scale are defined in Eq. (121).

2 0 2 Figure 2 :
Figure 2: Cross sections differential in the final-state momentum imbalance (or initial-state) k T = |k ⊥ | (upper plots), and the p T of the hardest jet (lower plots) for dijet production.On the left are the individual plots for the subtracted-real integral and the O ϵ 0 part of the integrated subtraction terms, for two values of the parameters E 0 , ζ 0 .Notice that the individual contributions flip sign going from one parameter set to the other.On the right are the sum of subtractedreal and integrated subtraction, clearly independent of the parameters within statistical accuracy.The processes included form a representative selection, and are given by gg ⋆ → 3g, gg ⋆ → u ūg, ug ⋆ → ugg, ug ⋆ → ud d, ug ⋆ → uu ū, and those with u ↔ d.The phase space and renormalization/factorization scale are defined in Eq. (121).

- 1 Figure 3 :
Figure 3: Cross sections differential in the final-state momentum imbalance (or initial-state) k T = |k ⊥ |.The essence of the plots is that the distributions vanish for k T → 0. The upper plots show LO distributions, for gg ⋆ → gg and gg ⋆ → ggg, and for two different k T -dependent PDFs.The lower row shows the O ϵ 0 -part of the real-radiation contribution to gg ⋆ → gg from the gg ⋆ → ggg channel, for both k T -dependent PDFs.The phase space and renormalization/factorization scale are defined in Eq. (121).The k T -dependent PDF KS-2013-linear is set to zero beyond its range of validity x < 0.01.