Analytic results for deep-inelastic scattering at NNLO QCD with the nested soft-collinear subtraction scheme

We present analytic results that describe fully-differential NNLO QCD corrections to deep-inelastic scattering processes within the nested soft-collinear subtraction scheme. This is the last building block required for the application of this scheme to computations of NNLO QCD corrections to arbitrary processes at hadron colliders.


Introduction
In this paper, we apply the nested soft-collinear subtraction scheme for NNLO QCD computations, introduced by some of us in Ref. [1], to the deep-inelastic scattering (DIS) proa e-mail: raoul.rontsch@cern.ch cess Pe → e + X . We note right away that our goal is not DIS phenomenology; rather, we would like to extend this subtraction scheme to processes that involve QCD partons in both the initial and the final state.
Compared to the cases of color singlet production [2] and decay [3] that we studied earlier, the DIS process requires us to deal with the situation where QCD partons in the (leading order) hard processes are not back-to-back. This makes the computation of NNLO QCD corrections to deep-inelastic scattering an important next step in the development of the nested subtraction scheme (NSS).
In spite of the fact that the computation that we report in this paper is new, we would like to emphasize that we can reuse significant parts of the analytic computations described in Refs. [2,3]. This is so because collinear singularities in QCD factorize on external lines so that their treatment, including analytic integration of respective subtraction terms [4], is process-independent. Hence, everything that needs to be known about collinear singularities in DIS and their regularization can be inferred from the treatment of the collinear singularities in color-singlet production and color-singlet decays, see Refs. [2,3].
At variance with collinear limits, important differences arise in the treatment of the (double-) soft radiation which is sensitive to the relative orientation of three-momenta of hard emittors. The integrated double-soft subtraction term for the case when the momenta of hard emittors are at an angle to each other was analytically computed in Ref. [5]. The computation of NNLO QCD corrections to the DIS process that we report in this paper is the first application of that result.
The main result of this paper is the set of analytic formulas that, in conjunction with the fully-resolved regulated contribution, provides a fully-differential description of DIS at NNLO QCD. It is our long-term goal to employ these for-mulas as ingredients to describe "initial-final dipole" contributions when computing NNLO QCD corrections to generic processes. Because of that, it is important to ensure that the analytic results for initial-final dipoles reported in this paper are correct. Studies of DIS are advantageous from this perspective since analytic results for DIS coefficient functions at NNLO are available [6][7][8] and we can use them to check our computations to a very high precision.
We note in passing that, in the past decade, a large number of subtraction schemes and slicing methods appeared ; they enabled a large number of NNLO QCD computations for important LHC processes . Nevertheless, in spite of all successes, the construction of a fully local, analytic, physically transparent and scalable subtraction scheme remains an interesting challenge. We believe that further development of the NSS, that we describe in this paper, will contribute to finding an answer to this challenge.
The remainder of the paper is organized as follows. In Sect. 2 we describe how leading order (LO) and next-toleading order (NLO) DIS cross sections are computed. In Sect. 3 we discuss the NNLO computation. In Sect. 4 we validate our results against analytic ones. We conclude in Sect. 5. Useful formulas are collected in several Appendices. Analytic results for NNLO QCD DIS subtraction terms in computer-readable format are provided in an ancillary file attached to this submission.

LO and NLO calculation
We consider deep-inelastic scattering of an electron on a proton P(P 1 ) + e( p 2 ) → e( p 3 ) + X, (2.1) mediated by a neutral current. The cross section of this process is computed as a convolution of parton distribution functions with the partonic cross section that describes partonelectron scattering. Schematically, we write In Eq. (2.2) we denote a parton of type i ∈ [−5, . . . , 5] as f i , with f 0 = g and f ±1,±2,±3,±4,±5 = {d/d, u/ū, s/s, c/c, b/b}. With a slight abuse of notation, we also use f i (x) to denote the parton distribution function of parton f i . The partonic cross section dσ f i +e→X can be computed in QCD perturbation theory as an expansion in the strong coupling constant α s . We write At leading order, electron-quark and electron-anti-quark scattering processes q( p 1 ) + e − ( p 2 ) → e − ( p 3 ) + q( p 4 ), q( p 1 ) + e − ( p 2 ) → e − ( p 3 ) +q( p 4 ), (2.4) contribute. For the purpose of computing QCD corrections, there is no difference between these two processes and we focus on the electron-quark scattering.
To compute the partonic cross section of this process, we employ the notation that has been used in earlier papers on the NSS [1][2][3], and define (2.5) where N includes normalization and symmetry factors, d = 4 − 2 is the space-time dimensionality, is the phase-space volume of a parton f i , A tree is the tree-level matrix element and O is a generic observable that depends on momenta p 1,.., 4 . E max is a sufficiently large but otherwise arbitrary 1 parameter that provides an upper bound on energies of individual partons; its role will become clear later.
We will also use the notation F LM (i, j, . . .) δ to indicate that the corresponding cross section is fully-differential with respect to momenta that are shown as arguments of the function F LM . In this notation, the fully differential LO cross section for quark-electron scattering reads 2s · dσ LO q = F LM (1 q , 2 e , 3 e , 4 q ) δ , (2.7) where s = 2 p 1 · p 2 is the partonic center-of-mass energy squared. We now discuss NLO QCD corrections. As we already emphasized in Refs. [1][2][3], at this order in the α s expansion, our subtraction scheme is equivalent to the FKS one [121,122]. In spite of that, it is useful to discuss the NLO QCD computation of DIS here, if only to develop a familiarity with our notation. At NLO QCD, both the quark q/q +e → q/q + e + g and the gluon g + e → q +q + e channels contribute to the DIS cross section. We consider the quark channel first, and start by discussing the real emission contribution. For the sake of definiteness, we focus on the following process (2.8) In analogy with Eq. (2.5), we define The scattering amplitude is singular when the gluon g 5 is soft or when it is collinear to the incoming or outgoing quark. Following our earlier work on the NSS [1][2][3], we introduce operators S 5 and C 51 , C 54 to extract the leading soft and collinear behavior of scattering amplitudes squared. 2 We use these operators to isolate non-integrable singularities in differential cross sections by systematically rewriting the identity operator I as (2.10) In both of the above equations, the first term describes a singular contribution and the second is free from soft and collinear singularities.
In the spirit of FKS subtraction [121,122], we partition the phase space using 14) 2 The definition of the S 5 and C 5i operators was described in detail in Refs. [1][2][3]. Here, we remind the reader that if f ( , while if f is Taylor-expandable around E 5 = 0 then S 5 f (E 5 ) = 0. Similar remarks apply to the collinear operators. 3 For the choice of Eq. (2.12), this follows immediately from C 5i ρ 5i = 0. that leads to simplifications in the collinear limit; in particular, the only collinear singularity present in the partition defined by w 5i is C 5i . This allows us to write 4 [1,4] C 5i (I − S 5 ) F LM (1 q , 4 q , 5 g ) + i∈ [1,4] (2.16) In Eq. (2.15), the first term on the right hand side describes the soft limit of the process Eq. (2.8), the second term describes two soft-subtracted collinear limits and the last term describes fully-regulated contributions that can be calculated in four dimensions. We continue with the discussion of the different terms in Eq. (2.15), starting with the soft contribution. We have where g s,b is the bare QCD coupling, C F = 4/3, and the definition of F LM (1 q , 4 q ) follows from Eq. (2.5). Since F LM (1 q , 4 q ) does not depend on p 5 anymore, we can integrate over the energy and angles of the unresolved gluon. We obtain where and η i j = ρ i j /2 = (1 − cos θ i j )/2. This allows us to write where we have introduced . (2.21) Note that at variance with cases of color-singlet production and decay that were studied in Refs. [2,3], the soft contribution depends non-trivially on the angle between the two hard emittors. Next, we discuss the soft-subtracted collinear terms in Eq. (2.15). We begin with the term proportional to C 51 that describes the situation when the collinear gluon is emitted by an incoming quark. We parametrize the gluon energy as and find where z min = 1 − E max /E 1 . The soft and collinear limits of the matrix element squared are well known. They read where z · 1 q denotes a quark with momentum z · p 1 and is the splitting function for this limit.
Since the matrix elements in Eq. (2.24) are independent of the emission angles, we can integrate over them. The relevant integral reads .
Putting everything together, we obtain (2.27) We note that by construction E max ≥ E 1 (see Footnote 1), so that z min ≤ 0. This implies that F LM (z · 1 q , 4 q ) = 0 for z ∈ [z min , 0] but the integration of the second term in angle brackets on the right hand side of Eq. (2.27) extends all the way to z min . We isolate the term in P qq that is singular in the z → 1 limit, and write (2.29) Note that in Refs. [2,3], we have chosen E max = E 1 but we prefer to keep E max generic in the current computation. Indeed, since final results are supposed to be E maxindependent, the possibility to vary this parameter provides a useful check on the implementation of the subtraction formulas. The plus distribution in Eq. (2.29) is defined as usual The discussion of the final-state collinear singularity, extracted by applying the operator C 54 to the matrix element squared, is very similar. In this case we define and repeat steps similar to the ones that led to Eq. (2.29). We obtain (2.32) We note that further details about final-state collinear splittings can be found in the discussion of QCD corrections to color-singlet decays, see Ref. [3].
To facilitate the -expansion of Eqs. (2.29) and (2.32), we write (2.33) The various splitting functions and anomalous dimensions are reported in Appendix A. We also define and write the real contribution to the NLO cross section Eq. (2.15) as follows [1,4] As the next step, we consider virtual corrections. Using notation analogous to Eq. (2.5), we define (2.36) We employ the Catani's representation [123] for the renormalized amplitudes to write the NLO contribution as 5 where (2.38) and We note that γ q can be found in Appendix A [see Eq. (A.4)] and that F fin LV is a finite remainder, free of any singularities. To obtain the final result for the NLO QCD cross section, we combine the real-emission contribution Eq. (2.35) with virtual corrections Eq. (2.37) and the contribution that originates from the collinear renormalization of parton distribution functions qq is the Altarelli-Parisi splitting function, see Appendix A. We find 2s · dσ NLO q = i∈ [1,4] (2.41) where the various (generalized) splitting functions and anomalous dimensions are defined in Appendix A. We have also defined where γ i = γ q (γ g ) and C i = C F (C A ) if particle i is a quark(gluon), and λ i j = 1 if both particles i and j are in the initial or in the final state, and zero otherwise. We also remind the reader that in our notation η i j = (1 − cos θ i j )/2, where θ i j is the angle between the directions of particle i and j. Comparing Eq. (2.41) to similar results for the production and decay of a color singlet, considered in Refs. [2,3], we note two main differences. First, Eq. (2.41) depends nontrivially on the relative angle between the incoming and out-going hard quarks. Second, subtraction terms in Eq. (2.41) explicitly depend on E max . This explicit dependence is supposed to be canceled by an implicit dependence contained in theÔ NLO F LM (1 q , 4 q , 5 g ) terms. Checking the E maxindependence provides a useful cross-check of the correctness of the implementation of Eq. (2.41) in a numerical program. Furthermore, we note that E max controls the amount of (soft) subtractions; by varying E max , we move contributions from the regulated hard emission termÔ NLO F LM (1 q , 4 q , 5 g ) to integrated subtractions. In this sense, E max is closely related to the so-called ξ cut parameters of the FKS formalism, c.f. Refs. [121,122].
In addition to the quark-electron scattering, at NLO we have to consider the gluon-electron scattering (2.43) The matrix element that describes this process is singular when the quark or anti-quark becomes collinear to the incoming gluon. These singularities are physically equivalent, so we find it convenient to treat both of them at once. To this end, we introduce the following partitioning and define

(2.45)
This effectively remaps both the gq and the gq singularities onto the p 1 || p 5 collinear configuration. Since a final state quark does not induce soft singularities, a subtraction formula for the gluon channel is simpler than the formula for the quark channel. Repeating the same steps that led to Eq. (2.41), it is straightforward to obtain (2.46)

NNLO calculation
In this section, we discuss the calculation of NNLO QCD corrections. Many details of the calculation are very similar to the color-singlet production and decay cases discussed in Refs. [2,3] and we do not repeat them here. Rather, we skim through the derivation of the subtraction formalism and concentrate on new features that arise in the DIS case.
As we already remarked in the previous section, the most important new feature is the fact that hard partons are not back-to-back. As the result, the integrated subtraction terms become functions of the opening angle between these partons. The second new feature is that we work with an arbitrary E max . As we explained in the previous section, the E maxindependence of the final result arises through a non-trivial interplay of subtractions and fully-regulated contributions. As such, E max provides both a powerful tool to check the correctness of the implementation of the subtraction framework in a numerical program and also allows us to shuffle contributions from numerical subtractions to analytically-integrated subtraction terms.
We find it convenient to deal with quark-and gluoninitiated processes separately. The primary reason for that is that only the former ones contain double-soft singularities, while in the latter case the only genuine NNLO singularities are of the collinear type. We start by discussing the quark channel.

Quark channel
We consider collision of an electron and a quark and write the differential NNLO partonic cross section as where • dσ VV q is the double-virtual contribution, which requires the one-loop squared and the two-loop amplitudes for the q + e → e + q process; • dσ RV q is the real-virtual contribution, which requires the one-loop amplitude for the q + e → e + q + g process; • dσ RR q is the double-real contribution, which requires the tree level amplitudes for the q + e → e + q + g + g, q + e → e + q + Q +Q, with Q = q, and q + e → e + q + q +q processes; • dσ PDF q originates from the collinear renormalization of parton distributions.
To efficiently manage the flavor structure and to follow what is commonly being done, we arrange the different partonic contributions into non-singlet and singlet pieces. We now briefly describe how this is done.
We consider the double-real contribution. Schematically, we write 6 We note that scattering amplitudes that involve an additional quark-anti-quark pair can be constructed from the "master" amplitude A qq;QQ shown in Fig. 1. For example, the amplitude for the process q( and the amplitude for the process q( By systematically re-labeling partons, it is straightforward to re-write dσ RR q Eq. (3.2) in terms of the amplitude A qq;QQ as (3.5) 7 Initial-and final-state crossings in A qq;QQ are understood. Crossings are reflected in a q vs.q mismatch in the flavor index between the subscripts and arguments of A.
We stress that, contrary to Eq. (3.2), the sums in Eq. (3.5) run over all quark flavors. It is conventional to split the real-emission processes into non-singlet, singlet and the interference contributions. To that end, we write The individual contributions in Eq. (3.6) are written as integrals of the corresponding amplitudes squared These amplitudes read The different contributions have a distinct structure of infra-red and collinear singularities. The strongest singularities are present in the non-singlet contribution |A ns | 2 which exhibits non-vanishing soft and collinear singular limits. More specifically, |A(1 q , 4 q , 5 g , 6 g )| 2 is singular if either one or both gluons are soft, or when they are collinear to quarks q 1,4 or to each other. Other contributions to the nonsinglet amplitude squared are less singular. For example, |A qq;QQ (4 q , 1 q ; 5 q i , 6q i )| 2 is singular when quarks q 5 and q 6 are both soft, or when they are collinear to each other, or when they are simultaneously collinear to q 1 or q 4 .
The singlet contribution |A s | 2 only contains collinear singularities. Double-collinear singularities arise when q 5 is collinear to q 1 or whenq 6 is collinear to q 4 . Triple-collinear singularities arise when q 1 , q 5 andq 6 are collinear or when q 1 , q 4 and q 5 are collinear. The pure interference contribution |A int | 2 is finite. Furthermore, for photon-mediated DIS, that we consider in this paper, the integral of this contribution vanishes due to Furry's theorem.
Since double-virtual and real-virtual corrections only contribute to the non-singlet cross section, we write where the three terms on the right hand side, defined as are separately finite. The collinear renormalization counterterms in Eq. (3.10) are explicitly given by 8 As usual, the sign "⊗" stands for the convolution product andP (0/1) are leading and next-to-leading order Altarelli-Parisi splitting functions, see e.g. [124]. We report them in Appendix A for convenience. The leading-order QCD β function, which appears in Eq. (3.11), reads where C A = 3, T R = 1/2 and n f is a number of massless quark flavors.

NNLO corrections in the non-singlet channel: derivation
The goal of this section is to describe the calculation of NNLO QCD corrections to neutral current DIS in the nonsinglet channel. Our goal in this discussion is two-fold. On the one hand, we aim to show that many ingredients of the current computation are similar to cases of color-singlet production and decay discussed in Refs. [1][2][3] and, for this reason, can be directly borrowed from these references. On the other hand, we want to emphasize new elements required for the construction of the nested subtraction scheme when a process involves color-charged initial-and final-state partons. We start by discussing the double-real contribution dσ RR q,ns . It contains double-soft singularities which arise when partons f 5,6 become soft, E 5 ∼ E 6 → 0. Because of this, we find it convenient to order energies of partons f 5,6 , see Ref. [1]. Using the notation of Sect. 2, we then define (3.14) To extract soft and collinear singularities from F LM,ns (1, 4, 5, 6), we closely follow the procedure described in Refs. [1][2][3]. First, we extract the double soft E 5 ∼ E 6 → 0 singularity. Similar to Refs. [1][2][3], we introduce an operator S that extracts the leading soft behavior of the matrix element, and sets E 5,6 to zero in both the momentum-conserving δfunction and the observable O in Eq. (3.13). We write  The second term on the r.h.s. of Eq. (3.15) is free of doublesoft singularities. In the first term, partons f 5,6 completely decouple from the hard matrix element and any infra-red safe observable. Explicitly, we have (3.16) where the function Eik (1,4;5,6) is the square of the double eikonal currents computed e.g. in Ref. [125]. Compared to the color-singlet production and decay cases described in Refs. [1][2][3], the integral of Eik(1, 4; 5, 6) depends on the relative angle between directions of hard partons f 1,4 . This integral was computed in Ref. [5] for a generic opening angle between f 1 and f 4 , so we can directly take the result from there.
The second term on the l.h.s. of Eq. (3.15) is free from the double-soft singularity, but still contains single-soft E 6 → 0 singularity. To extract it, we use I = S 6 + (I − S 6 ) and write (I − S)F LM,ns (1, 4, 5, 6) δ = S 6 (I − S)F LM,ns (1, 4, 5, 6) We deal with the first term on the right hand side of Eq. (3.17) following the discussion of the NLO computation, c.f. Eq. (2.20). The only differences with respect to that case are a) a more involved color structure and b) the maximal allowed energy of f 6 is now E 5 , because of the energy ordering. Taking into account that the S 6 singularity is only present if parton f 6 is a gluon, we obtain ( 3.18) The right-hand side in Eq. (3.18) still contains unregulated singularities that arise when g 5 is collinear to q 1 or q 4 . To extract them, we proceed as in the NLO case. To this end, we again use the partition of unity shown in Eq. (2.11) and write I = i∈ [1,4] ( , the first term on the right hand side of Eq. (3.19) leads to a fully-regulated contribution, while the second term extracts the collinear divergences. Its integration over the unresolved phase space proceeds similar to the NLO case except for two differences.
• When compared to NLO calculations, Eq. factor. This leads to modified powers of (1 − z) −2 in collinear limits of differential cross sections. To efficiently deal with this case, we find it convenient to define a generalised version of Eq. (2.33). It reads (3.20) where P qq is the splitting function defined in Eq. (2.25).
• The pre-factor in Eq. (3.18) leads to terms of the form This slightly changes the angular integral that needs to be computed. This point was already discussed in Refs. [1][2][3], and we refer the reader to these references for details.
Taking into account these two differences, and repeating steps explained in the context of the NLO calculation, we obtain for initial-state singularities, and for final-state ones. Combining the various terms discussed above, we arrive at the following formula where ellipses stand for various contributions from which all soft and collinear singularities have been extracted, as described above. However, the first term on the right hand side of Eq. (3.23) still contains unregulated collinear singularities. To proceed with their extraction, we follow the FKS approach [121,122] and its NNLO extension [18,19], and partition the phase space in the following way where w 5i,6 j are designed to dampen all but a few collinear singularities. More specifically, we ask that the damping factors behave in the following way: (3.25) We also find it convenient to construct the functions w 5i,6 j in such a way that w 51,64 and w 54,61 vanish in the limit when f 5 and f 6 become collinear to each other and that they are symmetric under the 5 ↔ 6 exchange. Apart from these requirements, the explicit form of the damping factors is mostly immaterial. An explicit construction of these factors, valid also for the DIS case, was discussed in Refs. [1][2][3]. We report it in Appendix B for convenience.
In the "double-collinear" partitions w 51,64 , w 54,61 , collinear singularities are effectively NLO-like. In the "triplecollinear" partitions w 51,61 and w 54,64 genuine triplecollinear singularities occur when partons f 5,6 become collinear to f 1 or f 4 , respectively. However, since these triplecollinear limits can be reached in a variety of ways, it is useful to introduce additional angular ordering to isolate physicallyinequivalent configurations [18,19]. For the partition w 5i,6i , we write We note in passing that the angular ordering Eq. (3.27) can also be enforced by constructing appropriate damping factors [36,37]. Nevertheless, we find it practical to employ Eq. (3.27) to isolate and extract collinear singularities and to integrate the subtraction terms analytically. In particular, a phase space parametrization that naturally implements the sector decomposition as in Eq. (3.27) and that we employ in this paper can be found in Refs. [18,19].

30)
• and, finally, the fully-regulated term reads F s r c r LM,ns (1, 4, 5, 6) δ = i∈ [1,4] θ (a) [14,41] ( (3.31) We note that in Eqs. (3.29,3.30,3.31) we used the notation introduced in Refs. [1][2][3]. In particular, C i denote triplecollinear limits, and collinear operators act on everything that appears to the right of them. For example, by writing C i j [d f j ] we indicate that the phase space for parton j has to be taken in the corresponding collinear limit, see Refs. [1][2][3] for details. A detailed discussion of double-and triple-collinear sectors, both for initial-and final-state collinear singularities, can be found in Refs. [2,3]. The fact that the discussion of these limits does not change is the consequence of the fact that collinear singularities only depend on color charges and types of external partons; as the result, once the production and decay of color singlets are understood, the description of similar limits in DIS naturally follows. For this reason we do not repeat the discussion of collinear limits per se but, instead, illustrate new features that arise in the DIS case by focusing on a few representative examples.
We start by discussing the contribution shown in Eq. (3.30). Compared to the cases studied in Refs. [2,3], the collinear limits now have an explicit E max dependence. For the triplecollinear limits, relevant results were obtained in Ref. [4], and we refer the reader to this reference for details. For the double-collinear case, we need to evaluate 9 (1,4,5,6) . (3.32) In the non-singlet channel, DC is non-vanishing only if both partons 5 and 6 are gluons. Schematically, Eq. (3.32) reads We note that apart from the operator S 6 and the energyordering condition, this expression is symmetric under the exchange f 5 ↔ f 6 . Accounting for that and renaming partons appropriately, it is easy to show that Eq. (3.33) can be written as where we used the short-hand notations C s r 5(6)i = C 5(6)i (I − S 5(6) ) and We note that each of the three terms in the curly bracket in Eq. (3.34) does not contain unregulated soft divergences. The first term is just the product of two NLO-like contributions; for this reason, it can be immediately read off from Eqs. (2.29,2.32). In the second term, the C 51 S 5 soft-collinear limit leads to 9 We note that to go from Eqs. (3.30) to (3.32), we used C 5i C 6 j w 5i,6 j = 1, which follows from the definition of the w 5i,6 j damping factors. (3.35) where the function is defined in Eq. (2.34). According to Eq. (3.34), we now have to take the soft-regulated collinear limit C s r 64 of Eq. (3.35). For the term proportional to E max , everything proceeds as in the NLO case. The term proportional to E 6 , however, induces an extra power of (1 − z) −2 after performing the change of variables E 6 = (1 − z)E 64 , see Eq. (2.31). As we already explained when discussing single-soft singularities, this leads to a term proportional to γ 24 qq , c.f. Eq. (3.20). More precisely, we obtain An analogous result can be found for the last term in the curly bracket of Eq. (3.34) that describes the regulated initial-state radiation and the soft final-state radiation. The last double-real contribution that we need to discuss is the soft-regulated single-collinear term Eq. (3.29). Once again, an in-depth discussion of this term can be found in Refs. [2,3], and we do not repeat it here. Rather, we illustrate the new features arising when considering the DIS process by focusing on the initial triple-collinear sectors w 51,61 θ (a,c) . Once these cases are understood, generalization to other sectors does not pose conceptual challenges and can be obtained following the discussion in Refs. [2,3]. Schematically, we write θ (a) where we used the fact that in the non-singlet channel the C 51 and C 61 limits are only singular if both partons 5 and 6 are gluons. To proceed further, we follow the discussion of the double-collinear contribution. We use thesymmetry of F LM,ns with respect to an interchange of g 5 and g 6 and re-write Eq. (3.37) as We now analyze the two terms in the curly bracket. We note that both of them only contain collinear divergences. Indeed, the soft singularity is always regulated, either by the explicit subtraction, as in the first term, or by the energy-ordering condition, as in the second one.
It is easy to see that the structure of the first term in Eq. (3.38) is nearly identical to the NLO case except that in the soft-collinear limit the upper bound on the energy of the parton f 6 is E 5 (contrary to E max in the NLO case). This observation allows us to immediately write the result for this contribution following the discussion in Sect. 2. We find (3.39) We note that functions P qq,R2 and , that appear in Eq. Finally, we note that the factor (ρ 51 /2) − , which is not present in the NLO case, arises from the ordering condition θ(ρ 61 < ρ 51 /2) in Eq. (3.38). Equation (3.39) still contains an unregulated collinear singularity that occurs when g 5 and q 1 become collinear. 10 We extract it in the usual way by inserting I = C 51 + (I − C 51 ). The regularization of the term proportional to P qq was 10 We note that the collinear p 5 || p 4 singularity is removed by the w 51,61 6||1 damping factor. explained in detail in Ref. [2]. The regularization of the term proportional to is analogous to what we just discussed for the double-collinear contribution.
We then move to the second term in curly brackets of Eq. (3.38), which corresponds to the nested soft-collinear limit. Since in the limit when g 6 is collinear to q 1 , the emission of the soft gluon g 5 can not resolve g 6 and q 1 , we are allowed to write (3.41) Similar to the NLO case, the dependence on the momentum of gluon g 5 has disappeared from the hard matrix element. However, a residual dependence on p 5 in w 51,61 6||1 and in the pre-factor (ρ 51 /2) − appeared. The dependence on the kinematics of gluon g 5 is described by the following integral where we defined (3.43) The contribution from the second term in the curly bracket in Eq. (3.38) is then 11 (3.44) At this stage, we can treat the collinear singularity as before. We note that the pre-factor E −2 6 will lead to additional overall factors of (1 − z) − , as discussed above. Taking this into account and repeating steps similar to what is done in the NLO case, we obtain This completes our discussion of the regularization of the triple-collinear sectors w 51 ,61 θ (a,c) . Finally, we note that this procedure is generic and that one can regulate all the remaining sectors following it. Before discussing the real-virtual and double-virtual contributions, we comment on the explicit dependence of Eq. (3.45) on the damping factor w. In general, one expects that 1/ poles of the double-real contribution are independent of the damping factors, since they do not appear in any other part (double-virtual, real-virtual etc.) of the calculation. Eq. (3.45) seems to contradict this assertion. We will now show that this is not the case. To this end, we note that the sum over double-collinear partitions in Eq. (3.29) leads to a contribution which is almost identical to Eq. (3.45 (3.45). The term 61 S 5 enters the differential cross section in the combination 61 S 5 / 2 , c.f. Eq. (3.45); this implies that we should prove that the dependence of 61 S 5 on the partitioning w only starts at O( 2 ). Since the integrand in Eq. (3.46) is singular in two collinear limits ρ 51 → 0 and ρ 54 → 0, we need to regularize and extract these singularities. Following the (by now) standard practice, we write (3.47) In the first term on the r.h.s. of Eq. (3.47) the dependence on the partitioning is absent because The second term on the r.h.s. of Eq. (3.47) is fully regulated and we can expand it in . We find Finally, we note that since the damping factors are process-dependent, it is not possible to analytically compute partitioning-dependent finite contributions once and for all. We do not view this as a problem. Indeed, it is simple to see that the 61 can be re-absorbed in a slightly different definition of the collinear operator in Eq. (3.31). We do not pursue this avenue further since, in the DIS case, we were able to compute the required integrals 61 S 5 analytically, using the damping factors given in Appendix B. This computation is outlined in Appendix C.
We now briefly discuss the real-virtual and double-virtual contributions to the partonic cross section Eq. (3.10). To discuss dσ RV q , we define where A 1−loop is the UV-renormalized one-loop amplitude. We then write (3.52) The function F LV (1 q , 4 q , 5 g ) contains a soft singularity that arises when the energy of g 5 vanishes, and collinear singularities that arise when the momenta of g 5 and g 1 or the momenta of g 5 and g 4 are parallel.
We now sketch the procedure to extract these singularities, and refer the reader to Refs. [1][2][3] for additional details. To extract the soft singularity, we write I = S 5 + (I − S 5 ). The soft limit of a generic masslessone-loop scattering amplitude was studied in Ref. [125]. Adapting the general result to our case, we write (3.53) where the term proportional to β 0 appears because we deal with UV-renormalized amplitudes. (3.55) For a particular choice of E max , the regularization of collinear singularities was discussed in detail in Refs. [1][2][3]. Generalization to arbitrary E max can easily be done following steps similar to the ones discussed for the double-real contribution; for this reason we won't discuss it further. At the end, the RV contribution is written as  (3.57) where (3.58) and F fin LV is a finite μ-dependent remainder. In Eq. (3.58), we use the notation The last term we need to discuss is the double-virtual contribution. We define To extract all 1/ poles, we use results of Ref. [123] and write where F fin LVV and F fin LV 2 are finite remainders, see Refs. [1-3] for more details. 12 In Eq. (3.61), we used the I 14 defined in Eq. (2.38), K = (67/18 − ζ 2 ) C A − 10/9 T R n f , and (3.62) We combine the double-real, real-virtual and double-virtual contributions described in this section with the PDF collinear renormalization Eq. (3.11) and obtain a fully-regulated finite final result. We present it in Sect. 3.1.2.

NNLO corrections in the non-singlet channel: results
The procedure outlined in the previous section allows us to rewrite the non-singlet NNLO differential cross section as dσ NNLO q,ns = dσ NNLO q,ns,3j + dσ NNLO q,ns,2j + dσ NNLO q,ns,1j , (3.63) where the three terms on the r.h.s. are individually finite and integrable in four dimensions. The term dσ NNLO q,ns,3j requires tree-level amplitudes with up to two additional partons relative to the Born configuration. It only receives contributions from double-real emission processes, and it is given by The second term, dσ NNLO q,ns,2j , requires tree and loop amplitudes with at most one additional parton relative to the Born configuration; it can be written as [1,4] [1,4] where the one-loop finite remainder F fin LV (1 q , 4 q , 5 5 ) is defined in Eq.  [1,4] w 5i,6i 6||5 ln , (3.66) with w 5i,6 j 6||k = lim η 6k →0 w 5i,6 j . Finally, the transverse metric is defined as where is the polarization vector of the gluon, and F μν LM is defined indirectly through the following equation (3.68) The vectors r (i) μ in Eq. (3.65) are remnants of spin-correlations, and can be thought of as a particular choice for the gluon polarization vector. Indeed, they satisfy r (i) · p 5 = 0 and r (i) · r (i) = −1. Their role in the subtraction framework and their explicit construction is discussed in details in Refs. [1][2][3], and we refer the reader to those references for more details. Here, we mention that if the momentum of gluon 5 is parametrized as p 5 = E 5 (1, sin θ 5i cos ϕ 5 , sin θ 5i sin ϕ 5 , cos θ 5i ), (3.69) the r vector reads r (i) = (0, − cos θ 5i cos ϕ 5 , − cos θ 5i sin ϕ 5 , sin θ 5i ). (3.70) The last term in Eq. (3.63) that needs to be consider is dσ NNLO q,ns,1j . It describes the exclusive process q + e → q + e, and only requires tree-level and loop amplitudes with Bornlike kinematics. It reads 13 (3.71) We note that finite parts of loop amplitudes have been defined in Eqs. (2.37, 3.61). The function S E i j can be found in Eq. (2.42). The terms i j and r μ r ν are the only contributions where the explicit dependence on the choice of partition functions appear. They are discussed in Appendix C, see Eqs. (C.14, C.16). They are multiplied by the generalized splitting function and anomalous dimension The two functions NS and T NS that appear in Eq. (3.71) contain the bulk of the NNLO (integrated) subtractions. They contain both standard and harmonic polylogarithms, and can be found in an ancillary file provided with this submission. All the other splitting functions and anomalous dimensions used in Eq. (3.71) are reported in Appendix A.

NNLO corrections in the singlet channel
We turn to the discussion of the singlet contribution to the cross section dσ NNLO q,s defined in Eq. (3.10). The singlet channel is much simpler than the non-singlet one discussed previously. This is so because (a) it only receives contributions from the double-real emission and the collinear renormalization of PDFs and (b) the singularity structure of the double-real contribution is very simple. In particular, it does not contain soft singularities and no genuine final-state collinear singularities. Indeed, the matrix element squared |A s (1, 4, 5, 6)| 2 defined in Eq. (3.8) is singular when p 1 || p 5 and/or when p 1 || p 5 || p 6 or p 1 || p 5 || p 4 . Since the two triple-collinear configurations are physically equivalent, we find it convenient to treat both of them at once. To this end, we first introduce a partition which is analogous to the one we used for computing NLO corrections in the gluon channel, c.f. Eq. (2.44), and write Then, we use the symmetry of the amplitude and the phase space to write Since |A 156 s (1, 4, 5, 6)| 2 does not contain any soft singularity, it is not necessary to order partons f 5 and f 6 in energy. Nevertheless, we find it practical to treat all contributions to the quark channel in the same way. We then define F LM,s (1,4,5,6) = i∈ [1,4] θ (a) [14,41] ( (3.80) where θ (a,b,c,d) i are defined in Eq. (3.27). We note that the θ (b,d) sectors don't require regularization since there is no single p 5 || p 6 collinear singularity. Because of this, one could easily do away with the sectors. Nevertheless, as we have already mentioned, we find it convenient to use the same parametrization for both the singlet and non-singlet quark channel, so we keep the sectors for the singlet contributions.
The second term on the r.h.s. of Eq. (3.79) can be written as where F LM,g and˜ 61 have been defined in Eq. (2.45) and Eq. (3.66) respectively, while all the splitting functions can be found in Appendix A. Finally, the fully-unresolved contribution reads where the (universal) transition function T S is reported in an accompanying ancillary file.

Gluon channel
In this section, we discuss NNLO QCD corrections in the gluon channel e + g → e + X . Similar to what we did in the quark case, c.f. Eq. (3.1), we write We note that, at variance with the quark channel, there are no double-virtual corrections in this case. We now briefly discuss the three terms on the right hand side of Eq. (3.83). Employing notation familiar from the discussion of the quark channel, we write the real-virtual contribution as The function F LV (1 g , 4 q , 5q ) contains unregulated collinear singularities when a quark or an anti-quark becomes collinear to the incoming gluon. To consider both singularities at once, we proceed as we did in the NLO case and define where the damping factor w 41 g is given in Eq. (2.44) and A 1−loop is UV-renormalized one-loop amplitude. We note that F LV,g , defined as in Eq. (3.85), is only singular when p 5 || p 1 . Since there are no soft singularities, E max does not play any role in the regularization procedure, which therefore follows the discussion in Ref. [2]. We refer the reader to that reference for details. The final result can be written as where (3.88) and L i j is defined in Eq. (3.59). We now discuss the double-real contribution dσ RR g . Similar to the quark channel, we write it as 2s · dσ RR g = F LM (1 g , 4 q , 5q , 6 g ) δ . (3.89) The matrix element for the process g( p 1 )+e( p 2 ) → e( p 3 )+ q( p 4 ) +q( p 5 ) + g( p 6 ) is singular in the following kinematic configurations: • q orq are collinear to the incoming gluon; • the outgoing gluon is collinear to the incoming one, or to the outgoing (anti)quark; • the outgoing gluon is soft.
Similar to the case of real-virtual corrections, the p 1 || p 4 and p 1 || p 5 singularities are equivalent. We then define 90) which is regular in the p 4 || p 1 limit. The regularization of the remaining singularities in the function F LM,g proceeds similarly to what we discussed in the case of the quark channels. There is only one main difference: since in this case there are no double-soft singularities, we do not order partons f 5 and f 6 in energy. This slightly changes the construction of the subtraction terms, as described in Refs. [2,3]. Tak = i∈ [1,4] θ (a) 14 We note that sometimes the action of the C i j operators is zero. For example C 54 F LM,g (1, 4, 5, 6) = 0, since the configuration where the two outgoing quarks are collinear to each other is not singular in this channel. Nevertheless, we retain the symmetric notation of Eq. (3.93) for convenience. +θ (d) [14,41] ( where F fin LV has been defined in Eq. (2.37), the various splitting functions can be found in Appendix A and T g is reported in an ancillary file that accompanies this paper.

Validation of results
In this section, we validate the results for the NNLO corrections to the deep-inelastic scattering process P + e → e + X obtained with our subtraction scheme by comparing them to analytic results for the NNLO DIS coefficient functions [6][7][8], as implemented in Hoppet [22,126,127]. Since the goal of this paper is to validate fully differential formulas for NNLO corrections to DIS, rather than to perform phenomenological studies of this process, we consider the simplest possible setup that allows for a thorough cross-check.
To this end, we only consider photon-induced neutral-current DIS. Furthermore, we only consider contributions proportional to either the gluon or the up-quark PDF. In other words, we define the non-singlet and singlet quark distributions as q s = q ns = u, which is sufficient for validating the results presented in this paper.
We now describe the setup of our computation. We consider photon-induced DIS collisions with hadronic centerof-mass energy equal to √ s = 100 GeV, and consider the total DIS cross section where the momentum transfer from an electron to a proton q 2 = −Q 2 is restricted to the interval 10 GeV < Q < 100 GeV. We include contributions of 5 massless flavors (2 up, 3 down) in the final state. We always use the NNPDF3.0 NNLO set [128] as implemented in Lhapdf [129] for both the parton distribution functions and the strong coupling.
In order to check the scale dependence of our result, we set the renormalization and factorization scales to μ R = μ F = Q max = 100 GeV instead of a more natural choice μ = Q. In order to study the robustness of our framework, we did not devise a specific parametrization for the phase space of the underlying DIS process. Specifically, we did not use a phase space that naturally accommodates the t−channel vector boson exchange. Hence, our phase-space parametrization is clearly not optimal. We believe that by not optimizing it we stress-test the numerical performance of our subtraction scheme. In general, we find that we can get per mill precision on the NNLO total cross section, corresponding to a few percent precision on the NNLO coefficient, in a few hours on an 8-core machine.
We where the subscript indicates whether the result has been obtained from our fully exclusive calculation ("NSS") or from the direct integration of the analytic coefficient functions ("an") over Q 2 and z. The Monte Carlo integration error for the former is shown in parentheses; for the analytic case, this error is always negligible, so we don't show it here. It follows from the above results that we can compute the NNLO DIS coefficients with a few per mill precision, and that the agreement between numerical results and analytical predictions is excellent. We have checked that this also holds true for other values of the factorization and renormalization scales. As we explained in Sects. 2 and 3, our framework contains a parameter E max which allows us to control the amount of (soft) subtraction. As such, one can view this as a prototype for a ξ cut parameter in the FKS formalism [121,122]. We have explicitly checked that our results are E maxindependent.
Finally, we note that we performed other checks by splitting numerical and analytic results into contributions of individual color factors. This allows us to cross-check subtle interference effects, which are color-suppressed and, hence, largely invisible in the full result for NNLO coefficients. We have found good agreement between numerical and analytic results for all such cases as well.

Conclusion
In this paper, we presented analytic results for NNLO QCD corrections to deep-inelastic scattering within the nested softcollinear subtraction scheme introduced by some of us in Ref. [1]. These results allow us to extend the nested subtraction scheme to processes involving partons both in the initial and in the final state. We have validated our calculation by computing NNLO QCD corrections to inclusive neutral-current DIS and comparing them against predictions obtained from a direct integration of analytic DIS coefficient functions. We found that despite a sub-optimal parametrization of the DIS phase space in the numerical routines, our formalism performed well and allowed us to check individual NNLO coefficients to a few per mill precision.
Apart from their relevance for processes like DIS or vector boson fusion in the factorized approximation, the results presented here constitute the last building block for applying the nested subtraction scheme to generic collider processes. Indeed, the nested subtraction scheme has been previously formulated for processes involving two hard partons both in the initial [2] and in the final [3] state. Since at NNLO the structure of infrared singularities is basically dipole-like, those results combined with the ones presented in this paper provide all the necessary building blocks to deal with arbitrary collider processes.
In practice, there are still two small issues that must be confronted when dealing with higher multiplicity reactions. First, the framework, as currently formulated, involves some partitioning-dependent contributions that must be dealt with in an efficient way, see the discussion around Eq. (3.49). We are confident that this issue can be dealt with by using a small modification of the subtraction operators. Second, for processes involving 4 or more partons, non-trivial color correlations appear. Although we have not studied such effects in detail yet, we do not anticipate that they would prevent us from extending the nested subtraction scheme to generic processes. We leave the investigation of these issues to the future.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical article and no experimental data has been listed.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .