NNLO QCD corrections to jet production in deep inelastic scattering

Hadronic jets in deeply inelastic electron-proton collisions are produced by the scattering of a parton from the proton with the virtual gauge boson mediating the interaction. The HERA experiments have performed precision measurements of inclusive single jet production and di-jet production in the Breit frame, which provide important constraints on the strong coupling constant and on parton distributions in the proton. We describe the calculation of the next-to-next-to-leading order (NNLO) QCD corrections to these processes, and assess their size and impact. A detailed comparison with data from the H1 and ZEUS experiments highlights that inclusive single jet production displays a better perturbative convergence than di-jet production. We also observe that the event selection cuts in some of the di-jet measurements of both H1 and ZEUS induce an infrared sensitivity that destabilises the perturbative stability of the predictions. Our results open up new opportunities for QCD precision studies with jet production observables in deep inelastic scattering.


Introduction
The HERA electron-proton collider produced a very precise data set on deeply inelastic scattering (DIS) processes, both fully inclusively as well as for specific hadronic final states [1,2]. Especially jet final states in deep inelastic scattering contain important information on the distribution of partons in the proton, which are a crucial ingredient for the prediction of any cross section at the LHC and at other future hadron colliders. Single jet inclusive and di-jet cross sections in deep inelastic scattering [3] are among the few precisely measured processes [4][5][6][7][8][9][10][11][12][13][14][15] that provide sensitivity to the strong coupling constant and the gluon distribution already at tree level. The importance of these cross sections for precision QCD studies has long been known. However, having the relevant hard sub-process coefficient functions only at next-to-leading order (NLO) in perturbative QCD [16][17][18][19], the HERA deep inelastic jet data could not be confronted with a sufficiently precise theory description to fully exploit their physics potential.
Calculations of collider observables to NNLO require to combine three types of partonlevel contributions, which are individually infrared divergent: double-real radiation, singlereal radiation at one loop and two-loop virtual contributions. To implement these contributions into a parton-level event generator, a method to extract and recombine their infrared (IR) singular parts is required. Our group has developed the antenna subtraction method [25,26] for this purpose. This method forms the basis of the NNLOjet code, which provides the necessary infrastructure and building blocks to implement NNLO corrections to different collider processes. Up to now, pp → Z + j [44], pp → H + j [42] and pp → 2j [51] have been implemented in NNLOjet.
The same framework is used in the calculation of NNLO corrections to jet production in DIS, and first results on di-jet final states were reported in a short letter [53]. This paper extends this calculation to single-jet inclusive production in DIS, and provides a detailed documentation of its implementation in the NNLOjet framework. Section 2 establishes the notation and describes the kinematical situation of jet production in DIS. The calculation of the NNLO corrections is described in Sec. 3, which also addresses in detail the handling of initial-state partons in the antenna subtraction method. Our results for the NNLO corrections to single jet inclusive production and di-jet production in DIS are discussed in detail in Secs. 4 and 5 where we also compare to the available data from the HERA experiments. We conclude with Sec. 6. P P ′ q = (0, 0, 0, −Q) k k ′ Figure 1. An illustration of the basic hard scattering process in the Breit frame with incoming proton momentum P , incoming lepton momentum k, virtual boson exchange momentum q, outgoing proton momentum P and outgoing lepton momentum k .

Kinematics of jet production in deep inelastic scattering
The basic interaction in deep inelastic lepton-proton scattering is mediated by a virtual gauge boson. The kinematics of the fully inclusive process can be inferred from the momenta of the incoming particles and of the outgoing lepton: such that a four-momentum q = k − k is transferred to the proton. Measurements are carried out in terms of the following variables (neglecting the proton mass): s lp = (k + P ) 2 , Q 2 = −q 2 , x = Q 2 2q · P , y = q · P k · P = Q 2 xs lp (2.1) The underlying parton-level process is the scattering of a quark from the proton with the virtual boson. At leading order, this quark carries a momentum fraction x of the proton momentum and scatters by an angle γ h , defined as cos γ h = (1 − y)xE p − yE l (1 − y)xE p + yE l , (2.2) with E denoting the energy. More detailed information on the underlying parton-level dynamics can be gained by examining the hadronic final state X. This is often analysed in the Breit frame of reference, which is defined by requiring proton and gauge boson momenta to take the form P B = (Q/(2x), 0, 0, Q/(2x)) , q B = (0, 0, 0, −Q) , (2.3) with Q = Q 2 . Momenta in the Breit frame are indicated by a subscript B. The Lorentz transformation from the laboratory frame to the Breit frame can be determined from the measured lepton kinematics. In this frame, the leading order DIS process results in an outgoing quark with vanishing transverse momentum, such that higher order contributions to DIS can be resolved by looking at final state objects with non-vanishing transverse momentum p T,B in the Breit frame.
Of particular interest is jet production in the Breit frame, which has a large cross section and allows the measurement of a variety of different distributions that can provide constraints on the parton content of the proton and on the strong coupling constant. Inclusive DIS at leading order is at O(α 2 α 0 s ) and only contains quarks in the initial-state; sensitivity to the strong coupling, α s and initial-state gluons only arises at NLO. In contrast, jet production is sensitive to α s and both initial-state quarks and gluons already at LO.
Jets are reconstructed in the Breit frame based on the transverse momenta p T,B (or equivalently transverse energies E T,B ), pseudorapidity η B = − ln tan(Θ B ) (with η B > 0 in the incoming proton direction) and azimuthal angle φ B of individual particles. Jets are formed from the particles by using an iterated clustering algorithm based on a distance measure [61]: where R 0 is a resolution parameter and δ = 1 for the k T algorithm and δ = −1 for the antik T algorithm. If the minimal distance measure is d ij , then objects i and j are recombined, if it is d iB , then object i is called a jet and removed from the list. The recombination of particle clusters is done in the E T scheme, as This scheme results in massless clusters and consequently massless jets, and p T,B = E T,B . The jets are accepted if they are within rapidity cuts and above a minimum transverse momentum. They are ordered in decreasing transverse momentum, and denoted as j1, j2 and so forth. The HERA electron-proton collider operated in its run II with beam energies E e = 27.5 GeV and E p = 920 GeV, resulting in √ s ep = 318 GeV. The HERA experiments H1 and ZEUS measured jet production cross sections [9,10,15] using the k T algorithm with R 0 = 1 as function of the Breit frame jet variables M jj = M 12 = (p j1 + p j2 ) 2 , E T,B = p T 2 = 1 2 (p T,B,j1 + p T,B,j2 ) , For di-jet production at leading order, ξ 2 can be identified with the momentum fraction of the incoming parton relative to the proton momentum.

Calculation of NNLO corrections
To calculate jet production at NNLO, all matrix elements which contribute to the process at the desired order in α s must be considered. Beyond leading order, radiative corrections can proceed via virtual loops or additional on-shell final-state particles, both of which can modify the cross section for an observable in non-trivial ways. We use dimensional regularisation with d = 4 − 2 space-time dimensions as regulator of both ultraviolet and infrared divergences in these contributions. Inclusive and semi-inclusive jet cross sections involve final-states containing up to one additional jet at NLO and up to two additional jets at NNLO, relative to the LO final-state. The theoretical prediction for a given n-jet final state at NNLO receives doublereal (RR) corrections from all relevant tree-level (n + 2)-parton matrix elements, realvirtual (RV) corrections which involve the interference between one-loop and tree-level (n + 1)-parton amplitudes, and double-virtual (VV) corrections involving the interference between two-loop and tree-level n-parton amplitudes as well as the one-loop n-parton contribution interfered with itself. For di-jet production in deep inelastic scattering, the relevant matrix elements have been known for a long time [54][55][56]. The RR and RV contributions are also part of the NLO corrections to tri-jet production in DIS [19], and the relevant matrix elements can now be generated automatically using standard tools [57][58][59][60].

Structure of the NNLO cross section
The colourless ingredients of the matrix element as well as the QCD coupling and colour factors can be pulled into an overall factor for each sub-process; this leaves the colourstripped matrix element as the basic object on which to perform the subtraction, which is appropriate as this is the function which contains all the IR singularities. The overall factor at leading-order for initial-state parton flavour i, is given by where V = N 2 − 1, N = 3 the number of colours, α s is the strong coupling and α is the fine structure constant. The initial-state colour averaging factor for the incoming parton is given by The initial-state spin averaging factors are given by where the number of spin states for the various initial-state particles are given by s e = s q = 2, s g = 2(1 − ), and the factors C( ) = (4π) 8π 2 e − γ ,C( ) = 8π 2 C( ) (3.4) which are arising from dimensional regularisation are introduced. The nomenclature for squared matrix elements we adopt in this paper is to denote the various types of colour stripped squared matrix elements with the shorthand:  Table 1. Summary of the scattering channels and matrix elements contributing at leading order. The argument with a caret denotes the initial-state parton, and the lepton momenta are omitted from the argument list.
where M can stand for A (all gluons), B (one quark-pair), C (two quark pairs of distinct flavour), D (interference terms for two quark pairs of identical flavour), n denotes the number of gluons in the process and the the number of loops. The case of two quark pairs of the same flavour can be decomposed into distinct-flavour matrix elements (C-type) and interference-only contributions (D-type). In the list of momentum arguments, the two outermost momenta denote the primary quark pair. A secondary quark pair is enclosed by semi-colons in the momentum list. In the case of matrix elements containing two quark pairs it is also necessary to distinguish which quark line the exchanged electroweak boson couples to. Sub-leading orders in colour are indicated byM , and closed fermion loop contributions byM . The different contributions to the cross section at leading order are given in Tab. 1. The antiquark-initiated processes are implemented explicitly in the NNLOjet program but considered only implicitly in this paper as they are closely related to the quark-initiated processes by relabelling.
To simplify the equations in this paper we absorb all additional powers of the coupling relative to leading-order into overall factors suitable for the correction being considered. We purposefully do not absorb additional colour factors into these factors so it is always clear which colour factor we are considering. To this end we define the following NLO factors for the real (R) and virtual (V) contributions, which are used to define the contributions to the cross section at NLO summarised in Tabs. 2 and 3. At NNLO we define the set of overall factors, level sub-process factor notation R qγ → qgg level sub-process factor notation Table 3. Summary of the gluon-channel partonic sub-processes and colour factors for the NLO calculation. The matrix element B γ,0 2 (3, i, j, 4) includes a summation {i, j} ∈ P (1, 2), where1 is the initial-state gluon. Each • denotes a closed quark loop.
which are used to consistently define the various contributions to the NNLO cross section in Tabs. 4 and 5. Occasionally, in the subtraction terms listed in App. A and App. B, we retain some information about the coupling of the vector boson to the quark line. In the case of a matrix element containing two distinct flavour quark lines (of flavours q and Q) in an unresolved limit where two of the quarks become collinear to form a gluon, it is often necessary to retain the information about which quark line remains. The reduced two-quark matrix element in the subtraction term is then denoted by either of the two terms, B γ, n,q , B γ, n,Q (3.11) to distinguish the different flavoured quark lines. Similarly it is sometimes necessary to use a matrix element which has been symmetrised over the final-state quark pair's momenta; level sub-process factor notation RR qγ → qggg Table 4. Summary of the quark-channel partonic sub-processes and colour factors for the NNLO calculation and the notation used for the colour-stripped squared matrix elements. Non-numeric arguments denote summation over permutations of partons, e.g. B γ,0 3 (1, i, j, k, 2) is summed over the six permutations of the labels {i, j, k} ∈ P (3,4,5). The lepton momentum labels are suppressed. Each • denotes a closed quark loop. level sub-process factor notation RR gγ → qqgg FB γ,2 1 (2,1, 3) Table 5. Summary of the gluon-channel partonic sub-processes and colour factors for the NNLO calculation. Non-numeric arguments denote summation over the set of gluons, which may include the initial-state gluon. e.g. the five-parton matrix element B γ,0 3 (4, i, j, k, 5) represents a sum over the permutations {i, j, k} ∈ P (1, 2, 3). The four-parton matrix element B γ,1 2 (3, i, j, 4) represents a sum over {i, j} ∈ P (1, 2). Each • denotes a closed quark loop.
we denote these matrix elements like so, Aside from the colour factors given in Tabs. 1-5 there are also contributions carrying the charge-weighted sum over flavours (3.13) These occur in multi-quark interferences between amplitudes where the vector boson couples to lines of different quark flavour or 1-loop matrix elements where the vector boson couples to the closed quark loop. The contributions coming from this overall factor are finite and only amount to a negligible correction to the cross section and are not included in our calculation.

Application of the antenna subtraction method
With the notation defined in Sec. 3.1 we can write the NNLO correction to the cross section, with incoming parton species a, in the form, where dσ RR a , dσ RV a dσ V V a are the contributions from tree-level, one-and two-loop squared matrix elements. They are decomposed into individual quark-initiated and gluon-initiated sub-processes in Tabs. 4 and 5. The dσ M F 1,2 a are the counter-terms arising from the massfactorisation of the physical PDFs.
It is well known that each of these contributions contains singularities which render them individually ill-defined: dσ RR a contains single-and double-unresolved divergences when integrated over the (n + 2)-parton phase space. dσ RV a + dσ M F 1 a contains explicit poles as deep as O( −2 ) in the dimensional regularisation parameter as well as singleunresolved divergences when integrated over the (n+1)-parton phase space. dσ V V a +dσ M F 2 a contains explicit poles as deep as O( −4 ).
To ensure the cancellation of all explicit poles in and render each phase space integral finite, one constructs three local subtraction terms and recasts the cross section in the form, Various methods for the construction of these subtraction terms at NNLO have been put forward in the literature [23][24][25][26][27][28][29]. We use the antenna subtraction method at NNLO [25,26]. In this method, antenna functions are introduced to account for all unresolved partonic radiation off a pair of two radiator partons [62]. The antenna subtraction term for each colour-ordered squared matrix element is then constructed as a sum (over all pairs of radiators) of products of antenna functions with reduced squared matrix elements of lower multiplicity. Each antenna function is integrated analytically over the sub-space of unresolved parton momenta in d = 4 − 2 dimensions and added back into a lower multiplicity final-state contribution as a Laurent series in , thus ensuring the analytic cancellation of all explicit poles. The following subsections provide a detailed description of the construction of the antenna subtraction terms in Eq. (3.15).

Phase-space mappings
To ensure that each antenna subtraction term converges to the unresolved limits of the full squared matrix element that it is intended for, the momenta in its reduced squared matrix element must be appropriate to this limit. However, by its construction, the antenna subtraction term is defined over the full phase space. To ensure the correct factorisation behaviour in all unresolved limits consequently requires a factorisation of the full phase space into a reduced phase space and a so-called antenna phase space corresponding to the unresolved radiation in each antenna subtraction term. The momenta of the reduced phase space are constructed from the original momenta by a phase space mapping, which is typically a non-linear transformation. This ensures the correct factorisation behaviour of the subtraction term in the unresolved limits, but also the factorisation of the phase space itself into disjoint spaces so that the analytic integration can be achieved over each antenna sub-space. The form of the mapping depends on the kinematic configuration of the two hard radiators involved in the antenna function, final-final (FF) [25], initial-final (IF) [63] or initial-initial (II) [63]. Since DIS processes have only one parton in the initial state, only the FF and IF mappings can occur. In each configuration there exists a mapping appropriate to single-unresolved (n + 1 → n particles) and double-unresolved (n + 2 → n particles) configurations. The subtraction terms are then constructed such that the antenna function depends on potentially unresolved unmapped momenta, whereas the reduced-multiplicity squared matrix element (and any event selection criteria) depends only on the mapped momenta.

Subtraction for double-real (RR) contributions
The double-real matrix elements contain both single-and double-unresolved divergences and the subtraction term therefore has to be able to regulate both types of limits. The single-unresolved limits of the matrix element, M γ,0 n+2 , when considering at least n jets, are dealt with by subtraction terms of the general form, where X 0 3 denotes an antenna function, which may in general be in the FF, IF or II configuration and is a function of three of the momenta from the (n+2)-parton momentum set, {p n+2 }. J (n+1) n represents the jet function which builds at least n jets from n + 1 partons. The reduced matrix element and the jet algorithm depend only on the (n + 1)parton mapped momentum set, {p n+1 }.
Colour-connected double-unresolved limits in the matrix elements are regulated using a four-parton antenna function and a double-unresolved momentum mapping. The subtraction terms for these limits take the form, Colour-disconnected and almost-colour-connected double-unresolved divergences (where unresolved partons share at most one hard neighbour in the colour ordering) are removed using subtraction terms with iterated single-unresolved mappings, These three types of subtraction term are sufficient to remove all potential divergences in the (n + 2)-parton phase-space integral and allow it to be computed in d = 4 using numerical techniques. The explicit forms of the RR subtraction terms for quark-and gluon-initiated channels are given in App. B.1 and App. B.2, respectively.

Subtraction for real-virtual (RV) contributions
The real-virtual matrix elements contain only single-unresolved divergences but unlike the double-real matrix elements, they also contain explicit poles in which have to be cancelled. The explicit poles of the (n + 1)-parton real-virtual matrix elements, M γ,1 n+1 , are cancelled analytically by a subtraction term of the form, where the function J 1,kl 2,ij is constructed [26] from the integrated form of the antenna function introduced in Eq. (3.16) and possibly also the kernels coming from the mass factorisation counterterms. The labels i and j denote the flavours of the hard partons and belong to the set {QQ, QG, GQ, GG} for quark-antiquark, quark-gluon, gluon-quark and gluon-gluon dipole functions. The labels k and l denote the kinematical configuration and belong to the set {F F, IF, II} for final-final, initial-final and initial-initial configurations respectively. In DIS kinematics only FF and IF configurations are present.
The single-unresolved divergences of the real-virtual matrix elements are regulated using two distinct subtraction terms, built to reflect the factorisation behaviour of oneloop amplitudes. The first term reflects the divergence associated with a tree-like singular function and a one-loop reduced matrix element and has the form, The second type of subtraction term reflects the divergence associated with a one-loop singular function and tree-like reduced matrix element and so involves a one-loop antenna function, The subtraction terms in Eq. (3.19) remove all explicit poles from the matrix elements yet contain their own single unresolved divergences, whereas the subtraction terms in Eqs. (3.20), (3.21) remove all single-unresolved divergences from the matrix elements yet contain their own explicit poles in . To cancel these remaining poles and divergences we introduce subtraction terms of the form, which render the RV contribution finite and integrable using numerical techniques. The explicit form for the RV subtraction terms in quark-and gluon-initiated channels are given in App. B.1.7 and App. B.2.7, respectively.

Subtraction for double-virtual (VV) contributions
The double virtual matrix elements, M γ,2 n , are integrated over the n-parton phase space and exhibit no IR divergences; however they do contain explicit poles in which are cancelled analytically using subtraction terms of the form, and The function J 2,kl 2,ij contains the integrated form of the four-parton antenna from Eq. (3.17), the integrated one-loop antenna from Eq. (3.21) and several other elements of lower complexity (essentially products of three-parton antenna functions) which conspire to cancel the poles of the two-loop matrix elements. The integrated antenna functions up to NNLO were derived in Ref. [25] for FF kinematics, in Ref. [64] for IF kinematics and in Ref. [65] for II kinematics. In DIS processes, only FF and IF kinematics contribute to the subtraction terms.
Once the RR and RV subtraction terms have been defined, the form the VV subtraction terms is completely fixed and so for brevity we do not present these explicitly in an appendix.

Initial-state identity changing collinear limits
In scattering processes with initial-state partons, there is an additional complication that arises when attempting to capture the singularity structure of the matrix element where the identity of the initial-state parton changes due to a collinear limit with a final-state parton: for example an initial-state gluon becoming collinear with a final-state quark, combining to form an initial-state quark. In this situation the matrix element will factorise, as usual, into a product of a splitting function and a reduced multiplicity matrix element, e.g. consider the gluon-initiated two-quark-two-gluon DIS sub-process (of which some indicative Feynman diagrams are shown in Fig. 2) in a collinear limit between the initial-state gluon and the final-state quark, where the momentum of the initial-state quark in the reduced matrix element is given by (3.28) We would like to regulate this kind of divergent limit using an appropriate subtraction term constructed from antennae and reduced matrix elements. A candidate subtraction term would be, . (3.29) This subtraction term does indeed mimic the matrix element faithfully in the 1 g ||3 q limit, which can be seen by considering: An attractive feature of the antenna subtraction method is that one antenna function regulates several unresolved limits. This is possible because of two facts: 1. the antenna function tends to the appropriate universal function in each unresolved limit.
2. the composite momenta for the reduced matrix element tend to the appropriate resolved momenta in each unresolved limit.
The former is ensured by a proper definition of the antenna functions, the latter achieved by an appropriate phase-space mapping which interpolates between several unresolved limits. For soft or collinear limits between final-state partons this procedure is unproblematic, as it is for collinear limits between initial-state partons and final-state partons where the species of the initial-state parton is unchanged (for example, a final-state gluon becoming collinear with an initial-state quark); we refer to these limits as identity-preserving (IP) limits. We refer the initial-final collinear limits that change the species of the initialstate parton as identity-changing (IC) limits. When a subtraction term, such as that in Eq. (3.29), attempts to regulate both IP and IC limits then an inconsistency arises due to the fact that we have to make a choice about the momentum assignment in the reduced matrix element. This determines in turn the momentum crossing in which the reduced matrix element is evaluated. For example, consider the matrix element in Eq. (3.27) and the subtraction term in Eq. (3.29). The matrix element contains, among others, divergences in the1 g ||3 q and1 g ||2 g initial-final collinear limits. As we have seen, the subtraction term in Eq. (3.29) successfully regulates the1 g ||3 q limit. In the1 g ||2 g the matrix element factorises according to, where now, The subtraction term in Eq. (3.29) also contains a divergence in the1 g ||2 g limit, This clearly does not regulate the matrix element in this limit because the reduced matrix element in the subtraction term has the quark as its initial-state parton. This choice is appropriate for the1 g ||3 q limit, but inappropriate for the1 g ||2 g which requires an initialstate gluon in the reduced matrix element. The phase-space mapping can interpolate between resolved momenta in different unresolved limits but cannot interpolate between different crossings of the reduced matrix element. There are several solutions to this problem of mixed IP and IC limits. A simple solution is to partial fraction the antenna such that it only contains either the IP or IC singularities, thus fixing the choice of crossing for the reduced matrix element. Following this approach we would divide the subtraction term into two terms, one handling each unresolved limit, e.g.
where the first line only contains the1 g ||3 q limit and the second line only contains the1 g ||2 g limit and they carry different reduced matrix elements. The disadvantage of this approach is that it relies on our ability to partial fraction the antenna function and successfully integrate the sub-antenna analytically. The three-parton antennae relevant for NLO calculations have been partial fractioned in this way and successfully integrated [63] so all ingredients are available to apply this method at NLO. However, the problem of mixed IP and IC limits is also present in three-parton one-loop and four-parton antennae used for NNLO subtraction terms. Partial fractioning these more complicated antennae is a non-trivial task and integrating the resulting sub-antennae analytically poses a substantial technical challenge. An alternative to this method, which requires no additional analytic integrations, is to construct subtraction terms which remove all IP limits from the matrix element. The antennae in these subtraction terms will generally contain IP and IC collinear limits, but the reduced matrix element is chosen to be appropriate for IP limits only. The strategy is then to use quark-antiquark antennae to remove the spurious IC limits from these subtraction terms. Once this is achieved, the genuine IC limits of the matrix element can be removed using the same set of quark-antiquark antennae but now carrying the crossing of the reduced matrix element appropriate to the genuine IC limit. For example, consider the combination of matrix elements, These contain the initial-final IP collinear limit1 g ||2 g , initial-final IC limits1 g ||3 q ,1 g ||4q, as well as several final-state soft and collinear limits which are irrelevant to this discussion. We can regulate these matrix elements with the subtraction term, The first three lines carry crossings of the reduced matrix element appropriate for IP limits (gluon initiated). The fourth line carries the crossing of the reduced matrix element relevant to the IC limits (quark initiated). The first two lines regulate the1 g ||2 g IP limits but contain spurious IC1 g ||3 q and1 g ||4q limits. These spurious singularities are then removed by the quark-antiquark antenna in the third line, which by necessity also carries a gluon-initiated matrix element. The fourth line is similar to the third insofar as it captures IC limits, but now uses a quark-initiated reduced matrix element to regulate the genuine limit. Taken together, the block of terms in Eq. (3.38) successfully disentangles and removes all IP and IC divergences for the block of matrix elements in Eq. (3.37). The advantage of this method is that although it is more cumbersome at NLO, as in the example considered here, the same method can be applied to combinations of threeparton one-loop antennae and four-parton antennae to disentangle the NNLO IP and IC limits. This method requires no new integrals and all ingredients are readily available such that the problem is reduced to constructing an appropriate subtraction term. For example, the subset of RR matrix elements for gluon-initiated DIS, contain the IC triple-collinear limits (1 g ||3 q ||i g ), (1 g ||3 q ||j g ), (1 g ||4q||i g ), (1 g ||4q||j g ), as well as the IP triple-collinear limits (1 g ||i g ||j g ), (3 q ||i g ||j g ), (4q||i g ||j g ). Following the strategy outlined above at NLO, we construct a subtraction term that regulates the IP limits using full antennae, Of course, these subtraction terms contain spurious IC limits which can be removed using quark-antiquark four parton antennae, The genuine IC triple-collinear limits of the matrix elements in Eq. (3.39) can then be regulated by the same set of quark-antiquark antenna, but now carrying a quark-or antiquarkinitiated reduced matrix element, 1 The alternative would be to partial fraction the D 0 4 antenna into sub-antennae which is significantly more difficult than the partial fractioning of the D 0 3 antenna due to the presence of overlapping single-and double-unresolved limits. The integration of such a sub-antenna would also present significant technical challenge and additional master integrals.
In practice, we employ a mixture of approaches to the subtraction terms constructed for this calculation. We use the known NLO partial fractioned antennae where possible to simplify the subtraction term and use a carefully constructed combination of antennae for the genuine NNLO IC limits such as the triple collinear or one-loop single collinear IC limits. In principle all partial fractioned antennae could be eliminated from the subtraction terms, but this would only serve to make the construction more baroque so we take an optimal approach and use partial fractioned NLO antennae where possible.

Implementation into NNLOjet and validation
The NNLOjet code is a parton-level event generator that provides the framework for the implementation of jet production processes to NNLO accuracy, using the antenna subtraction method. Besides containing the event generator infrastructure (Monte Carlo phase-space integration, event handling and analysis routines), it supplies the unintegrated and integrated antenna functions and the phase-space mappings relevant to all kinematical situations. The implementation of processes in the NNLOjet framework requires the availability of the matrix elements for all RR, RV and VV processes, and the construction of the antenna subtraction terms. NNLOjet provides testing routines to verify the point-wise convergence of the subtraction, as documented for example in Ref. [66]. Processes included in NNLOjet up to now are Z and Z + j production [44], H and H + j production [42] as well as di-jet production in hadron-hadron collisions [51].
Our NNLOjet implementation [53] of jet production in DIS uses the same matrix elements [54][55][56] as were used in Z + j production [44], however in different kinematical crossings. While the phase space for Z + j production corresponds to a single crossing region of the matrix elements, three different crossing regions are required [16,67] to describe each DIS process, depending on the relative size of Q 2 compared to the partonparton invariants. The matrix element contributions are summarised in Tabs. 4 and 5. The antenna subtraction terms relevant to each of these contributions are collected in the Appendix.
To validate our implementation of the tree-level and one-loop matrix elements, we compared the NLO predictions for di-jet and tri-jet production against SHERPA [68] (in DIS kinematics [69]), which uses OpenLoops [57] to automatically generate the one-loop contributions at NLO. The antenna subtraction is then verified by testing the convergence of subtraction terms and matrix elements in all unresolved limits and by the infrared pole cancellation between the integrated subtraction terms and the two-loop matrix elements.

Scale dependence of the NNLO cross section
The coupling constant renormalisation and mass factorisation procedures are performed at a renormalisation scale µ r and factorisation scale µ f , where the strong coupling constant α s (µ r ) and parton distribution functions f i (x, µ f ) are evaluated, respectively.
Beyond leading order, the fixed-order contributions to the parton-level cross sections contain an explicit dependence on µ r and µ f , which compensates the dominant scale dependence from the coupling constant and parton distributions at the previous orders. These scale-dependent terms can be predicted from the renormalisation group equations (our normalisation conventions are summarised in the appendix of Ref. [42]). Starting from the evaluation of the DIS di-jet cross section at fixed default values µ f = µ r = µ 0 (which can be chosen dynamically event-by-event): the full scale dependence of the cross section can be predicted in terms of It takes the form: In the above expressions, summation over the parton indices is implicit. Equation (3.45) can be used to compute the cross section at multiples of an initially chosen scale, and was also employed to perform detailed validations of our implementation of the different contributions to the NNLO corrections to di-jet production.

Inclusive jet production
Inclusive jet production in deep inelastic scattering has been widely studied by the H1 [4, 5, 7-10] and ZEUS [12][13][14] experiments at DESY HERA. The jet measurements are preformed in the Breit frame, where the transverse momentum requirement on the jet ensures the existence of a partonic recoil, even if only a single jet is reconstructed in the kinematical acceptance.
In this section, we use the kinematic criteria (see Tab. 6 below) used in the final H1 measurements [9,10] to discuss several generic features of the NNLO corrections to inclusive jet production, followed by an in-depth comparison of the newly derived predictions to the H1 data [9,10].

Structure of the inclusive jet production cross section at NNLO
Inclusive jet production in the Breit frame receives leading order contributions both from quark-initiated and gluon-initiated processes. The relative magnitude of these depends on the final-state kinematics. Figure 3 shows the relative contributions from quark and gluon initial states to inclusive jet production, evaluated at NLO and NNLO for a representative range in Q 2 , using the cuts from the H1 low-Q 2 analysis (see below). We observe that at low transverse momentum, p T ≤ 10 GeV, inclusive jet production is mainly gluon-initiated (to almost 80%). With increasing p T , the quark-initiated contribution becomes more and more important, reaching 50% at around p T ≈ 30 GeV, and further increasing towards higher p T . NNLO corrections affect the relative importance of the different initial states only in a moderate manner, and only at high-p T . Compared to NLO, the gluon-initiated fraction decreases more slowly for larger values of p T . The overall behaviour of the parton fractions and their modification between NLO and NNLO remains largely unchanged at higher Q 2 . The inclusive jet production cross section receives contributions from different jet multiplicities. At NNLO, final states with up to four identified jets contribute. The jets are ordered in transverse momentum. Figure 4 displays the contribution of the first, second, third and fourth jet to the inclusive jet distribution at NNLO for a given bin in Q 2 . It can be seen that the leading jet and subleading jet contribute about 85% and 12% to the distribution over the full p T range. This behaviour can be understood from the fact that the jet production is measured in the Breit frame, where the leading order process will always yield a p T -symmetric two-jet final state. The jets are not necessarily both found inside the rapidity cut, such that in some fraction of the events, only the leading jet is observed. Furthermore, real radiation from higher order corrections can lower the transverse momentum of the second jet compared to the first one, such that the same event will enter the p T distribution at a larger value with the first jet than with the second jet. Due to the sharp decrease of the distribution with increasing p T , the relative importance of the   Table 6. Kinematical cuts used to define the inclusive jet phase space in the measurements of H1 at high-Q 2 [9] and low-Q 2 [10], and ZEUS [13] .
second jet is lower than of the first jet. The contributions from the third and fourth jet are at the level of a few per-cent and a few per-mille respectively at low p T . Their importance decreases to higher p T , which can be understood from the limited final-state phase space that is available for multi-jet production at large transverse momenta.

Comparison to HERA data
Inclusive jet production in the Breit frame (using the inclusive k T algorithm with a massless E T recombination scheme) was measured double differentially in Q 2 and p B T by the H1 experiment, which distinguishes a low-Q 2 [10] and a high-Q 2 sample [9] and by the ZEUS experiment [14]. The event selection criteria for all three studies are summarised in Tab. 6. We compute theoretical predictions at LO, NLO and NNLO, always using the same set of parton distribution functions (NNPDF3.0 NNLO) with α s (M Z ) = 0.118. Our predictions use a dynamical central scale choice µ 2 r = µ 2 f = (Q 2 +p 2 T,B )/2, and the theory uncertainty is determined from a seven-point scale variation with rescaling factors [1/2, 2]. The theoretical predictions are multiplicatively corrected for hadronization effects, using the correction tables from the experimental papers [9,10]. Figure 5 compares our NNLO predictions to the H1 data. We observe that the NNLO corrections are very substantial at low-Q 2 and low-p B T , with an up to 60% enhancement with    respect to NLO. These large corrections are within the NLO uncertainty band (close to the upper edge), and result in a residual theory uncertainty of 20% even at NNLO. Especially at low Q 2 , the shape and normalisation of the theory prediction changes significantly going from NLO to NNLO, and results in a considerably improved theoretical description of the data, as already statistically quantified in the experimental H1 study [10]. A very similar pattern is also observed for the ZEUS measurement shown in Fig. 6. With increasing Q 2 , the size of the NNLO corrections decreases, accompanied with very small residual theoretical uncertainties (decreasing from 10% at Q 2 = 150 GeV 2 to 2% at 5000 GeV 2 ). In this region, the combination of precision data with the newly derived NNLO corrections has clearly the potential to provide important new constraints for precision QCD studies.

Di-jet production
In the Breit frame, di-jet production and single jet inclusive production in deep inelastic scattering are mediated by the same Born level processes and are closely related. In contrast to single jet inclusive production, where only the rapidity and transverse momentum of the jet can be studied, di-jet final states allow for more kinematical observables to be reconstructed (see section 2 above). Typically, di-jet cross sections are measured inclusively based on the two leading jets in an event, i.e. including events with more than two reconstructed jets. Inclusive di-jet production was measured by the H1 [4, 6-10] and ZEUS [11,13,15] experiments at DESY HERA.
In this section, we adapt the event selection (see Tab. 7 below) used in the final H1 measurements [9,10] to discuss several generic features of the NNLO corrections to inclusive jet production, followed by an in-depth comparison of our NNLO predictions to the ZEUS [15] and H1 [9, 10] data.

Scale setting in the di-jet production cross section at NNLO
The dependence of the NNLO cross section on the renormalisation and factorisation scales has been derived in Sec. 3.5, where it can be seen that each order in the perturbative series compensates the dominant scale dependent terms from the previous order. Consequently, the residual scale dependence of a fixed-order prediction is commonly used to estimate the error on the prediction resulting from missing higher order corrections. The scale dependence of the theoretical prediction is quantified by choosing central values for µ f and µ r , and then evaluating the cross section for variations around these central scales, typically by a factor two. These central scale values should reflect the dynamics of the process. In processes involving only a single mass scale (such as for instance inclusive deep inelastic scattering, depending only on the photon virtuality Q 2 ), the central scale choice is unambiguous (at most up to a constant factor). For processes involving multiple physical scales, several choices are possible (and a priori equally well motivated). The only guiding principle for choosing the central scales in this case is the occurrence of large logarithmic terms in each order in perturbation theory, which spoil the convergence of the perturbative expansion and indicate an inappropriate choice of the central scale.
Di-jet production in DIS depends on two scales: the photon virtuality Q 2 and the average transverse momentum of the two jets p B T 2 . Using the kinematical cuts of some of the bins in the H1 low-Q 2 di-jet measurement (which are discussed in detail in the following subsection) as an example, we study different choices for the central scales in Fig. 7. Table 7. Kinematical cuts used to define the inclusive di-jet phase space in the measurements of H1 (high-Q 2 [9] and low-Q 2 [10]) and ZEUS [15].
We compare the following three options: All three options were considered previously in comparisons of NLO predictions to the H1 and ZEUS jet data [4][5][6][7][8][9][10][11][12][13][14][15]. Option (a) represents the average of both hard scales in the jet production process, and is used as default throughout this paper; option (b) is used frequently in the experimental studies, with the argument that the partonic structure of the proton (µ f ) is resolved by the virtual photon, while it is hard QCD emission (µ r ) that yields the two-jet final state; finally option (c) is entirely based on the photon virtuality to describe the hardness of the interaction.
We observe that the scale choices (a) and (b) yield similar predictions, except in the region of large transverse momentum. Especially at large Q 2 , the choice (b) results in slightly higher cross section predictions, accompanied with larger scale uncertainties. In contrast, scale choice (c) results in unphysical predictions (negative cross sections) if applied at low Q 2 and large transverse momentum. These results confirm earlier observations made at NLO [16][17][18][19]. By examining the analytical form of some of the NLO contributions [17,18] for scale choice (c), the emergence of large logarithmic corrections in the jet transverse momenta could be established. These corrections are largely compensated in the hard coefficient functions for scale choices (a) and (b), which are clearly preferable in terms of reliability and perturbative stability. It should be pointed out that the large logarithmic terms alone (which can be inferred from threshold resummation [70]) do not properly account for the bulk of the NNLO corrections, as observed in Ref. [10].

Comparison to HERA data
Inclusive di-jet production was measured by both HERA experiments: H1 [9, 10] provides double-differential results in Q 2 and p B T 2 or Q 2 and ξ 2 , using the k T and anti-k T jet algorithms in the Breit frame. ZEUS [15] uses only the k T jet algorithm and provides single-differential results inĒ B T = p B T 2 , Q 2 , M jj , η * , ξ = ξ 2 as well as double-differential     The theoretical predictions are multiplicatively corrected for hadronization effects, using the correction tables from the respective experimental papers [9,10,15]. Figure 8 displays the inclusive di-jet cross section for the ZEUS kinematics as a function of the electron variables Q 2 (left) and of x (right). We observe the NNLO corrections to be sizeable especially at low values of Q 2 or x, where they enhance the NLO prediction by about 15%. In this region, the scale dependence of the NNLO prediction is as large as at NLO, or even larger. We also note that the shape of the data is better described by the NLO prediction than at NNLO. A similar pattern is observed in the distributions in p B T 2 and M jj shown in Fig. 9, with sizeable NNLO corrections in the lower range of the distributions. In both these distributions this range clearly correlates with the approach to the infrared limit, as can be seen even more clearly in the double differential distribution in p B T 2 and Q 2 , Fig. 10. In this limit, the QCD coupling constant increases and logarithmically enhanced contributions could deteriorate the convergence of the perturbative fixed-order expansion. This issue is further aggravated by the interplay of the M jj cut (see Tab. 7) with the transverse momentum requirements on the final state jets, as discussed previously in Ref. [53], and elaborated upon in more detail below. The relatively small scale dependence of the NLO predictions in this range is likely an artefact originating from a cross-over of the upper and lower edges of the scale-band, as investigated in detail for hadronic di-jet   Figure 11. Inclusive di-jet production cross section as a function of η * (left) and log(ξ 2 ) (right), compared to ZEUS data.
production in Ref. [51]. The scale variation at NNLO therefore provides the more realistic assessment of the theoretical uncertainty. The di-jet cross section as function of η * and of log(ξ 2 ) is shown in Fig. 11. While good perturbative convergence is observed in the plateau region η * < 0.65, NNLO corrections turn out to be very sizeable at higher rapidities. The perturbative instability in this region was already pointed out and explained by the ZEUS collaboration [15]. The log(ξ 2 ) correlates most directly with the parton distributions, indicating that the ZEUS di-jet data will likely deliver important constraints in the determination of the gluon distribution for momentum fractions in the medium range 0.01 < ξ < 0.1. The double-differential distribution in log(ξ 2 ) and Q 2 , Fig. 12, further illustrates this impact, showing a coherent pattern of discrepancy between data and NNLO theory over the full Q 2 range. It will be very interesting to further study the impact of these data in the determination of parton distributions at NNLO.
The H1 dataset is divided into a high-Q 2 and a low-Q 2 sample, see Tab. 7. For the low-Q 2 sample [10], double-differential distributions were measured in Q 2 and p B T 2 , which we compare to our NNLO calculation in Fig. 13. Compared to NLO, the NNLO corrections enhance the prediction of the di-jet cross section at lower values of p B T 2 , leading to a considerable improvement in the description of the H1 data, as already pointed out in Ref. [10]. Moreover, we observe an excellent convergence of the perturbative series and a considerable reduction of the theory uncertainty in going from NLO to NNLO. This highlights the potential of these data for future precision studies of parton distributions and the strong coupling constant.  Figure 12. Inclusive di-jet production cross section as a function of log(ξ) in bins of Q 2 , compared to ZEUS data.  Figure 13. Inclusive di-jet production cross section as a function of p B T 2 in bins of Q 2 , compared to H1 low-Q 2 data.  Figure 14. Inclusive di-jet production cross section as a function of p B T,2 in bins of Q 2 , compared to H1 high-Q 2 data. For the high-Q 2 sample, slightly different event selection criteria are applied: in particular, a minimum value of M jj is imposed. In a previous work [53], we have already studied the impact of the NNLO corrections in the double differential distributions in Q 2 and p B T 2 for the H1 high-Q 2 sample. Figure 14 collects these results. It can be seen that the improvement of the theoretical uncertainty from NLO to NNLO is less pronounced than for the low-Q 2 study, and that perturbative instabilities arise at low p B T 2 . These can be traced back to the M jj cut [53], which restricts the available LO and NLO phase space for di-jet production at low p B T 2 . Figure 15 compares the NNLO predictions for double-differential distributions in Q 2 and ξ 2 to the H1 high-Q 2 measurement (these distributions are not available for the H1 low-Q 2 study). We observe that the quantitative behaviour is very similar to the ZEUS distributions, Fig. 12. At LO, ξ 2 is directly related to the momentum fraction carried by the incoming parton, such that Fig. 15 indicates the kinematical range where the H1 data can potentially improve the determination of parton distributions. Recalling the definition of ξ 2 (2.4), we moreover observe that the H1 high-Q 2 data set typically probes lower values of ξ 2 than its low-Q 2 counterpart (which is due to the transverse momentum requirement on the final state jets, and contrasts with the kinematical correlations in inclusive DIS).
To illustrate the problematic impact of the symmetric cuts on p B T , combined with a cut on M jj , we re-evaluated the double differential distribution in Q 2 and ξ 2 for a different set of jet cuts: p B T,j1 > 5 GeV, p B T,j2 > 4 GeV. The result is shown in Fig. 16, where we observe a very substantial improvement in the perturbative convergence, compared to the cuts used in the H1 analysis [9].

Conclusions and Outlook
In this paper, we described the calculation of the second-order (NNLO) QCD corrections to jet production in deep inelastic scattering. By defining jets in the Breit frame of reference, this process requires at least two partons in the final state, thereby providing sensitivity on the gluon distribution and the strong coupling constant already at leading order. We consider inclusive production of single jets and of di-jet systems in the Breit frame, which start both at the same perturbative order.
Our calculation uses the antenna subtraction method to cancel infrared singularities among parton-level sub-processes of different multiplicity. The application of this method to processes with one hadron in the initial state requires the dedicated treatment of infraredsingular splittings that change the parton identity. Our calculation is implemented into the NNLOjet parton-level event generator framework, which was also used recently for the NNLO corrections to pp → Z + j [44], pp → H + j [42] and pp → 2j [51].
The HERA experiments H1 and ZEUS have measured inclusive single jet and di-jet production over a broad kinematical range. We observe that the NNLO corrections to inclusive single jet production modify the shape of the kinematical distributions, which are now described considerably better than at NLO. The corrections are moderate in size (ten to twenty per cent) except for very low jet transverse momenta or low photon virtuality Q 2 , and their inclusion substantially reduces the scale uncertainty on the predictions, typically well below the experimental statistical and systematical uncertainty.
The NNLO corrections to di-jet production are found to be sizeable, and often well outside the scale uncertainty band of the NLO predictions over a broad kinematical range for most of data sets from ZEUS and H1. Already in an earlier study [53], we could relate this behaviour to the interplay between the cuts on the jet transverse momentum and di-jet invariant mass, which overly restricts the final state phase space at LO and NLO. We now provide further evidence for this argument by performing dedicated comparisons of the predictions with and without invariant mass cut, with the latter displaying much-improved perturbative convergence and reduced scale uncertainty. Only one H1 di-jet data set (low-Q 2 [10]) was measured without the invariant mass cut; this data set is well-described by the NNLO theory.
Our NNLO results open up several opportunities for precision phenomenology with jet observables in deep inelastic scattering. The determination of parton distributions at NNLO from a global fit is currently dominated by processes that are only quark-initiated at leading order (inclusive DIS, Drell-Yan processes). Constraints on the gluon distribution come mainly from indirect effects (scaling violations) or from the inclusion of data from processes (like jet-production) where the full NNLO corrections are unknown. In these cases, the NNLO corrections to the hard process cross sections are either approximated using some ad-hoc assumptions or discarded altogether. A recent study [71] on the impact of LHC inclusive top quark cross section data on the determination of the gluon distribution illustrated the importance of a consistent treatment of NNLO effects (which are known for top quark production [39,40]). Our newly derived NNLO corrections to jet production in DIS enable the consistent inclusion of HERA data on this process into global parton distribution fits at NNLO. The magnitude of the corrections, as well as their kinematical dependence, makes it likely that their inclusion will lead to modifications of the gluon distribution, also leading to a substantial reduction of its uncertainty in the crucial region of medium-x.
Jet production data from HERA have been used to measure the strong coupling constant α s [7,9,12]. The error on all these measurements was dominated by the theory uncertainty inherent to the NLO predictions used in the extraction of α s . This uncertainty was found to be typically a factor two or more larger than experimental statistical or systematical uncertainties, thereby proving to be the limiting factor to further improving α s measurements from jet production in deep inelastic scattering. Given that the analysis of inclusive deep inelastic scattering structure function data typically yields values of the α s at the lower boundary of the range indicated by other determinations [72], it will be very interesting to apply the NNLO corrections derived in this paper to perform an NNLO-accurate determination of the strong coupling constant from DIS jet production data.

Acknowledgments
We would like to thank Nigel Glover and Thomas Morgan for many interesting discussions during this project, Stefan Höche and Marek Schönherr for help with the NLO comparisons against SHERPA and Daniel Britzger for many discussions on the H1 jet measurements.

A NLO subtraction terms
In this appendix we include all the subtraction terms for the NLO two-jet calculation. The matrix elements and their associated factors are collected in Tabs. 2 and 3.

A.1 Quark-initiated subtraction terms
The quark channel receives contributions from the processes qγ → qgg, qγ → qqq at tree level and qγ → qg at 1-loop.

Real
Real emission corrections at NLO contain up to three final-state partons and are comprised of B γ,0 2 , C γ,0 0 and D γ,0 0 matrix elements.

A.1.1
The leading-colour contribution is composed of two-quark-two-gluon matrix elements, B γ,0 2 , which are summed over both permutations of the final-state gluons. This matrix element and the corresponding subtraction term are given by where P (i, j) denotes the permutations of the labels in the set of final-state gluons {3, 4}. The subtraction term is constructed for a single permutation and given by The two-quark-two-gluon matrix element also contributes at sub-leading colour where it is given by theB γ,0 2 matrix element. The subtracted contribution to the cross section is given byB The gluons in this function act as if they are abelian due to the particular combination of interferences used to define it. This changes the factorization behaviour compared to the leading-colour matrix element and we construct the subtraction term to reflect this: 2 ({p} 2 ) +A 0 3,q (1, i, k) B γ,0 1 (1, j, ( ik)) J With three final-state partons and an initial-state quark we must also evaluate the contribution from four-quark matrix elements of identical (qγ → qqq) and non-identical (qγ → qQQ) flavours. The full four-quark matrix element, summed over possible quark flavours, can be written as a term derived from the non-identical flavour matrix element, C γ,0 0 , and an interference term, D γ,0 0 . The C γ,0 0 function contains contributions from diagrams where the vector boson couples to each quark line and the subtraction term is given by In this subtraction term we note the appearance of the B γ,0 1,q and B γ,0 1,Q reduced matrix elements which were introduced in Sec. 3.1 and which carry over the information on which quark line the vector boson is coupling to in the reduced matrix element.
We can see that the first term of the subtraction term regulates the quark-antiquark collinear limit between partons i and j which constitute a quark line of flavour Q and so in the reduced matrix element the vector boson couples to the remaining quark line of flavour q. The remaining subtraction terms, on the other hand, regulate the collinear limit between the quarks1 and k, which are of flavour q. Consequently, the vector boson in the reduced matrix element couples to the remaining quark line of flavour Q.

A.1.4
As explained above, the four-quark matrix element also contains a term given by the interference of four-quark amplitudes with different quark assignments, arising from the identical flavour contribution, We include this contribution but note that it is finite in all unresolved limits and so requires no subtraction term.

Virtual
The virtual NLO contributions to the cross section arise from the interference of a one-loop amplitude with a tree-level amplitude for the process qγ → qg. Once all colour factors are stripped out, this leaves three separate contributions to the cross section.

A.1.5
The leading-colour contribution and its subtraction term is given by All the poles of the one-loop matrix element are cancelled analytically by this subtraction term.

A.1.6
The two-quark-one-gluon one-loop matrix element also gives a contribution at subleading colour,B γ,1 The dipole function in this subtraction term, which contains all the explicit poles, is a function of the quark and antiquark momenta but not the gluon momentum. This can be traced back to the colour connection of the virtual gluon in the loop; at subleading colour the virtual gluon is colour connected only to the quark line.

A.1.7
The one-loop diagrams which contain a closed quark loop form a distinct matrix element, B γ,1 1 , carrying a factor of N F . The poles of this matrix element are cancelled by the subtraction termB γ,1,T 1 . We must also include the integrated subtraction terms coming from Sec. A.1.3 where the reduced matrix element is gluon-initiated rather than quarkinitiated.
This type of subtraction term was discussed in Sec. 3.3 and is referred to as an identity changing (IC) term. Such a terms clearly cannot cancel the poles of the one-loop matrix element as it is proportional to the gluon-initiated reduced matrix element. Instead, the relevant integrated antennae from the real subtraction term are combined with the IC mass factorization kernels to generate a finite virtual subtraction term,B γ,1,T 1,q→g .

A.2 Gluon-initiated subtraction terms
The gluon channel is composed of the tree-level process gγ → qqg and the one-loop process gγ → qq.

Real
Due to the initial-state gluon, no four-quark matrix elements contribute to the real correction, but only two-quark-two-gluon matrix elements with one of the gluons crossed into the initial state.

A.2.1
The leading-colour contribution is, as in the quark channel, given by the B γ,0 2 function. We sum over all gluon permutations, where now one of the gluons is in the initial-state. The subtracted contribution to the cross section is given by The final two terms in the subtraction term account for IC limits where the initial-state quark becomes collinear with a final-state quark or antiquark.

A.2.2
The subleading-colour contribution to the process gγ → qqg is given by theB γ,0 2 function with one of the gluons crossed into the initial-state, The initial-state gluon cannot become soft, but can participate in a collinear limit with final-state partons. As the gluons in this matrix element behave as if they were abelian, no gluon-gluon collinear limits exist and all IF collinear limits are therefore IC limits with the final-state quark or antiquark,

Virtual
The virtual NLO contributions arise from the interference of a one-loop amplitude with a tree-level amplitude for the process gγ → qq.

A.2.3
The leading-colour contribution contains a subtraction term which cancels the poles of the one-loop matrix element and an IC subtraction term which is finite, and the finite IC subtraction term is given by

A.2.4
The subleading-colour matrix element is given by theB γ,1 1 function, whose poles are cancelled by the subtraction termB γ,1,T 1 . The IC integrated subtraction terms coming from Sec. A.1.2 are combined with the IC mass factorization kernels to form the finite subtraction term,B γ,1,T 1,g→q . The total subtracted contribution to the cross section is given bỹ The one-loop matrix element with a closed quark loop is given by theB γ,1 1 function. This function contains explicit poles, but as is clear from Tab. 3, there is no real subtraction term contributing to this colour factor. This means there are no integrated antenna functions to cancel the explicit poles of the matrix element. However, the NLO mass factorization term does contribute to this colour factor and cancels the poles of the matrix element. This is reflected in the fact that the integrated dipole for this process is formed entirely from mass factorization kernels, as is noted in Ref. [26]:

B NNLO subtraction terms
In this appendix, we include the subtraction terms that are used to regulate each colourstripped matrix element introduced in Sec. 3.1. The matrix elements associated with the subtraction terms are listed in Tabs. 4 and 5.

B.1 Quark-initiated subtraction terms
The quark channel at NNLO receives contributions from the tree-level sub-processes qγ → qggg and qγ → qqqg, the one-loop correction to the processes qγ → qgg and qγ → qqq and the two-loop corrections to the process qγ → qg.

Double real
The tree-level five-parton contributions are given by the colour decomposition of the twoquark-three-gluon and four-quark-one-gluon matrix elements with the quark crossed into the initial-state.

B.1.1
The leading-colour matrix element is summed over all six permutations of the final-state gluons. The subtraction term is summed only over the cyclic permutations of gluons. The cancellation of divergences occurs once all terms in the sum are evaluated, where P C denotes the restricted sum over cyclic permutations. The subtraction term is given by TheB γ,0 3 matrix element contributes at subleading colour, 1/N 2 relative to leading colour. The interferences of colour-ordered amplitudes which form this function result in one of the gluons behaving as if it were an abelian gauge boson, colour connected only to the quark line. The two remaining gluons are colour connected to the quark line but also to each other. The abelian-like gluon is always the first gluon argument in the functionB γ,0 3 . The subtracted contribution to the cross section is given by The abelian-like nature of the gluon i (first argument) in the subtraction term can be seen by the fact that there are no divergent collinear limits for this gluon with other gluons in the subtraction term, ×A 0 3 (1, i, ( lj)) B γ,0 1 (1, ( kj), ( i (lj))) J The four-quark-one-gluon matrix element contains contributions from identical and nonidentical quark flavours for the two quark strings. As was the case at NLO in Secs. A.1.3-A.1.4 the full four-quark contribution can be separated into a term derived from the nonidentical flavour matrix elements, given by the C γ,0 1 functions, and a remainder of interference terms given by the D γ,0 1 functions. The C γ,0 1 terms include both colour orderings and the vector boson coupling to both quark lines. The subtracted contribution to the cross section is given by − E 0 3,q →g (k, l, 1) a 0 3,g→q (j, 1, ( kl)) B γ,0 1,Q (( j (kl)), i, 1) J In this subtraction term we once again retain information on the flavour of the quark line in the reduced matrix element and employ a symmetrisation over the quark and antiquark for some reduced matrix elements indicated by aB.

B.1.5
The subleading-colour contribution to four-quark-one-gluon scattering can also be split into a term derived from non-identical flavour matrix elements,C γ,0 1 terms, and a left-over set of interferences,D γ,0 1 terms. TheC γ,0 1 terms contain all colour orderings and vector boson couplings to both quark lines. The subtracted contribution to the cross section is given by  1, ( il)) B γ,0 1,Q (k, 1, ( j (il))) J The interferences of quark orderings arising from the identical flavour contribution to fourquark-one-gluon scattering gives the subtracted contribution to the cross section, where +2 C 0 4 (1, k, j, l) B γ,0 1,q (1, i, ( jkl)) J The interferences of quark orderings also contribute to the subleading colour and require subtraction,D where the subtraction term is given bỹ 2 ({p} 2 ). (B.14)

Real-virtual
The real-virtual contributions to the quark channel receive contributions from the twoquark-two-gluon and four-quark one-loop matrix elements, where the four-quark matrix elements include identical and non-identical flavoured quarks. Each matrix element may be decomposed into its colour-stripped functions, as set out in Tab. 4.

B.1.8
The leading-colour contribution is given by the B γ,1 2 function, which is then summed over both gluon permutations. The subtracted contribution to the cross section is given by where the subtraction term for a single ordering is given by 2 ({p} 2 ) 2 ({p} 2 ) The subleading-colour matrix element contains a mixture of colour connection structures for both the virtual gluon loop and the external gluons. This leads to a set of subtraction terms which are colour connected only to the quark line, and a set which are also colour connected to gluons. The subtracted contribution to the cross section is given by The most subleading-colour contribution contains only abelian-like colour connections and therefore only contains quark-antiquark antennae and quark-antiquark integrated dipoles. The subtracted contribution to the cross section is given bỹ The two-quark-two-gluon matrix element with a closed quark loop contributes to two colour factors. The leading contribution is given by theB γ,1 2 (1, i, j, 2) function. The poles and unresolved limits of this matrix element are removed by the subtraction termB γ,1,T 2 (1, i, j, 2), and both of these functions are summed over the permutations of final-state gluons.
As was the case at NLO, there is also a contribution to the subtracted cross section from the integrated IC subtraction terms in Sec. B.1.4. These are combined with the IC mass factorization terms to generate a finite subtraction term with no unresolved limits, B γ,1,T 2,q→g . The total subtracted contribution to the cross section is given by The subleading-colour contribution to the two-quark-two-gluon matrix element with a closed quark loop,B γ,1 2 , has a subtraction terms which cancels all poles and all divergences,B γ,1,T 2 , and a finite subtraction term which is generated from the IC subtraction terms in Sec. B.1.5 and IC mass factorization counterterms,B γ,1,T 2,q→g . The total subtracted contribution to the cross section is given bŷ wherễ 2 ({p} 2 ) The contributions with four quarks and a closed quark loop can be written purely in terms of the non-identical flavour matrix elements with a closed quark loop,Ĉ γ,1 0 . The subtracted contribution to the cross section is given bŷ (1; j, k; i) = The remaining four-quark one-loop interferences between different quark orderings also contribute to two colour factors. The leading-colour term has the subtracted contribution to the cross section given by The matrix element contains explicit poles but no single unresolved divergences for the same reason that the tree-level D γ,0 0 was finite in all limits. The subtraction term is then given by The subleading-colour contribution to the identical-quark interferences similarly has no single-unresolved limits and requires only the poles to be removed. The contribution to the cross section is given byD (1, j; k, i) =

B.2 Gluon-initiated subtraction terms
The gluon channel receives contributions from the sub-processes: gγ → qqgg and gγ → qqqq at tree level, gγ → qqg at one-loop and gγ → qq at two loops. Each of these contributions can be separated into their respective colour factors which are listed in Tab. 5.

Double real
The double-real correction to the gluon channel contains two-quark-three-gluon, and fourquark-one-gluon matrix elements with one of the gluons crossed into the initial-state.

B.2.2
The subleading-colour matrix element,B γ,0 3 , is summed over all gluon permutations, including the initial-state gluon. The subtraction term is constructed in a way such that it is summed only over final-state gluon permutations.

Real-virtual
The gluon-induced real-virtual corrections have at most three partons in the final state and therefore only the two-quark-two-gluon one-loop matrix elements are necessary. These matrix elements have a colour decomposition given in Tab. 5 and we present the subtraction terms for the colour-stripped matrix elements in this section.

B.2.8
The leading-colour matrix element is given by the function B γ,1 2 , which is summed over both gluon orderings. The subtraction term B γ,1,T 2 , is constructed to cancel against the full set of orderings and is thus not summed over any permutations of gluons. The IC terms from the RR and the IC mass factorization kernels combine to form a finite subtraction term, B γ,1,T 2,g→q . The total contribution to the cross section is given by