Single-jet inclusive rates with exact color at O ( α 4 s )

Next-to-next-to-leading order QCD predictions for single-, doubleand even triple-differential distributions of jet events in proton-proton collisions have recently been obtained using the NNLOjet framework based on antenna subtraction. These results are an important input for Parton Distribution Function fits to hadron-collider data. While these calculations include all of the partonic channels occurring at this order of the perturbative expansion, they are based on the leading-color approximation in the case of channels involving quarks and are only exact in color in the pure-gluon channel. In the present publication, we verify that the sub-leading color effects in the single-jet inclusive doubledifferential cross sections are indeed negligible as far as phenomenological applications are concerned. This is the first independent and complete calculation for this observable. We also take the opportunity to discuss the necessary modifications of the sector-improved residue subtraction scheme that made this work possible. ar X iv :1 90 7. 12 91 1v 1 [ he pph ] 3 0 Ju l 2 01 9


Introduction
Pure jet observables are not only interesting within the portfolio of Standard Model (SM) measurements but also as tools for New Physics searches. When paired with high-precision predictions within Quantum Chromodynamics (QCD), they are an important input for Parton Distribution Function (PDF) fits. Consistency of the latter application requires theoretical predictions at the same order of perturbation theory as the specified accuracy of the given PDF set. In consequence, it is nowadays indispensable to determine differential distributions of jet events at the next-to-next-to-leading order (NNLO) in QCD. Although not the subject of the present publication, non-perturbative effects from underlying event and hadronisation [1] as well as electroweak corrections [2] should also be included for a realistic comparison with measurement data.
A calculation of jet rates at NNLO in QCD remains a very challenging project. Until now, results could only have been obtained with the help of the antenna subtraction scheme [3][4][5][6][7][8][9][10][11][12][13][14][15][16] as implemented in the NNLOjet framework [17]. While early studies concentrated on the pure-gluon case [18,19], several publications [20][21][22][23] appeared recently that present results including the complete set of partonic channels. It turns out that quark-gluon and quark-quark scattering dominates jet rates at high transverse momenta and/or rapidities. These contributions are, therefore, necessary for a complete description. On the other hand, the current NNLOjet implementation is based on the leading-color approximation for all channels but the pure-gluon channel. It has been argued in the past that inclusion of the sub-leading color effects at O(α 4 S ) while keeping exact color dependence for the lower order contributions will have a negligible impact on the predictions. Nevertheless, it is necessary to verify this statement by explicit calculation. This is the main purpose of the present publication.
The published results on jet rates correspond to several different setups. In particular, [20,21] correspond to 7 TeV ATLAS data, while [22] to 8 TeV and [23] to 13 TeV CMS data. Single-jet inclusive cross sections (every identified jet in an event is accounted for in the histogram) that are double-differential in the jet transverse momentum, p T , and rapidity, |y|, have been published for 7 and 13 TeV Large Hadron Collider (LHC) center-of-mass energies and various jet radii, R, defined with the anti-k T jet algorithm [24]. Furthermore, for the same center-of-mass energies, there are available results for di-jet cross sections that are double differential in the di-jet invariant mass, m jj , and the jet rapidity difference, |y * | ≡ |y 1 − y 2 |/2. Finally, there is a published result for di-jet cross sections corresponding to 8 TeV CMS data that is triple-differential in the average transverse momentum, p T,avg ≡ (p T,1 + p T,2 )/2, rapidity difference, |y * |, and the boost, y b ≡ |y 1 + y 2 |/2. For convenience of the reader, we summarise the available results in Appendix A.
Di-jet cross sections have consistently been evaluated with central renormalisation and factorisation scales µ R = µ F = µ = m jj . However, it turned out that the scale choice for single-jet inclusive cross sections is a non-trivial issue, since seemingly well-justified choices lead to large differences in the predictions as illustrated by the results for 7 TeV center-ofmass energy using the customary scales µ = p T (each jet input into the histogram with the cross section evaluated with the scale equal to that particular jet transverse momentum) -2 -and µ = p T,1 (cross section scale set to the hardest-jet transverse momentum). The question of scale setting has been very thoroughly studied in [23]. The final recommendation of that publication is to use either the jet-based scale µ = 2p T or the event-based scale µ =Ĥ T (the scalar sum of the transverse momenta of the partons in the event).
The high computational cost of the calculation of jet rates with our software (see Section 3) enabled us to only perform one complete Monte Carlo simulation for the present publication. Since the setups of the different theoretical predictions described above differ in energy, jet transverse momentum cuts and jet radii, we had to choose a single specific setup to study the sub-leading color effects. Our goal was to compare results for a classic observable used in PDF fits that has been previously evaluated with one of the recommended scales. In view of these considerations, we have chosen to evaluate the single-jet inclusive cross section for 13 TeV center-of-mass energy available from [23] (jet radius R = 0.7) using the jet-based scale µ = 2p T . We should also point out that there are no available numerical values for cross sections or cross sections ratios in either of the publications [20][21][22][23] -only histograms have been provided. Appendix D contains our results for the K-factors, i.e. ratios of cross sections evaluated at NNLO and NLO in QCD with the same PDF set. This is a first step towards an easy inclusion of jet data in NNLO PDF fits. In the future, it would certainly be desirable to provide fastNLO [25][26][27] or APPLGRID [28] tables for all setups. We intend to undertake this task once more computational resources become available.
Besides the study of sub-leading color effects in jet rates, there is another aspect to the present publication. The calculation of a cross section at NNLO in QCD requires a method to handle infrared (IR) singularities occurring in contributions of different final state multiplicities. Apart from the already mentioned antenna subtraction scheme, there are several other methods currently being developed for this purpose: the CoLoR-fulNNLO scheme [29][30][31][32][33], q T -slicing [34,35], N -jettiness slicing [36][37][38][39][40], sector-improved residue subtraction [41][42][43] and its spin-off called nested soft-collinear subtraction scheme [44][45][46], the projection-to-Born method [47,48], local analytic sector subtraction [49] and geometric IR subtraction [50]. The results that we report in Section 3 have been obtained with our implementation of the sector-improved residue subtraction in the C++ library Stripper (SecToR Improved Phase sPacE for real Radiation). However, we have introduced several improvements w.r.t. to Ref. [43]. These improvements imply for instance a minimal number of subtraction terms per phase space point. They also require a modified reduction of the construction in Conventional Dimensional Regularisation (CDR) to four dimensions ('t Hooft-Veltman scheme). In the present publication, we discuss these issues in detail.
The paper is organised as follows. In the next section, we discuss modifications of the sector-improved residue subtraction. These consist of an improved phase space parameterisation and a new approach to the dimensional reduction of the formulation of the scheme. The section is closed with details on the implementation and tests. Subsequently, we present our results for the single-jet inclusive cross sections at 13 TeV. The main text is closed with an outlook. Appendices provide and overview of published results on jet cross sections, define the notation for cross sections contributions, provide a list of expressions -3 -necessary for the implementation of the subtraction scheme in four dimensions, and, finally, provide numerical values of the NNLO K-factors.

Minimal sector-improved residue subtraction
We define a subtraction scheme to be minimal, if it has the minimal number of subtraction terms (defined by final state kinematics) for a given phase space point. Consider a phase space configuration with n+2 final state particles contributing to a cross section at next-tonext-to-leading order of perturbation theory. A configuration corresponding to a single soft (one gluon with vanishing energy) or collinear (two partons collinear) limit will have n + 1 resolved final state particles (single-unresolved configuration), while a configuration corresponding to a double-soft (two gluons or a quark-anti-quark pair with vanishing energy) or triple-collinear (three partons collinear) limit will have n resolved final state particles (double-unresolved configuration). Let us divide the phase space into sectors according to collinear limits (see Section 3 of Ref. [43] for more details). A sector that only allows for one singular configuration (divergent cross section contribution) in a specific single collinear limit is called single-collinear. A sector that allows for one singular configuration in a specific limit of two pairs of collinear partons is called double-collinear. Finally, a sector that allows for one singular configuration in a specific limit of three collinear partons is called triple-collinear. It is easy to convince one-self that the minimal number of subtraction terms in a single-collinear sector is one. Similarly, the minimal number of subtraction terms in a double-collinear sector is three (two for single-unresolved configurations and one for the double-unresolved configuration). Finally, the minimal number of subtraction terms in a triple-collinear sector is four (three for single-unresolved configurations and one for the double-unresolved configuration). The phase space construction of the sector-improved residue subtraction scheme as defined in Ref. [43] generates more subtraction terms. Here, we present an alternative phase space that, due to the additional sector decomposition in the triple-collinear sector, never requires more than three subtraction terms for a given phase space point. It turns out, however, that the four-dimensional formulation of the scheme requires modifications with this phase space. On the other hand, the methods presented here allow for numerical checks of pole cancellation in CDR for a fixed Born phase space point [51,52]. In contrast, an implementation of Ref. [43] yields finite results at the level of distributions only.
We note, finally, that a single subtraction configuration in single-collinear parameterisations has also been obtained in the FKS subtraction scheme [53] implementations of Refs. [54,55]. The approach presented here is conceptually different, and allows to cover next-to-leading and next-to-next-leading order cases on the same footing. It might be extensible to even higher orders. Let us introduce the following notation for the phase space measure corresponding to a single particle with mass m ≥ 0 and momentum k in d dimensions The complete phase space for a process involving n q = 0 1 (not explicitly parameterised) resolved momenta q i , n u = 1, 2 unresolved momenta u i , and 0 < n r ≤ n u reference momenta r i , n f r of which in the final state, can be decomposed as follows where n = n q + n f r + n u and P is the total initial state momentum. In the second line, we have inserted an intermediate momentum q with invariant mass Q ≥ nq i=1 m i . In the third line, however, we have performed the integration over q, leaving an integration measure in the square brackets, which only depends on the fixed invariant mass of q.
The reference and unresolved momenta are the momenta that are allowed to correspond to a singular configuration in a given sector. Reference momenta and unresolved momenta have different behaviour w.r.t. the soft limit: a soft reference momentum does not generate a singularity.
We now introduce a mapping from the full phase space to the Born phase space which is invertible for fixed unresolved momenta and which conserves the invariant mass of the intermediate state q The mapping only involves a rescaling of the reference momenta. Specifically, for a final state reference momentum we require where x is given by a function f x of the full kinematics, in particular of r. The phase space measure is modified . (2.7) For an initial state reference, we require similarly where z is given by a function f z of the full kinematics, in particular of r. We consider the phase space measure together with the integration over the parton momentum fraction with the parton distribution function φ where p h is the initial state hadron momentum. The Born configuration is only properly defined, i.e. physical, if the rescaling parameters are always non-negative andq 0 ≥ 0. The phase space measure may be rewritten as where θ {u l } ∈ U represents constraints on the unresolved momenta, and J is the Jacobian of the transformation. Furthermore, we have used Eq. (2.5) and the implied existence of a Lorentz transformation Λ, q = Λq, together with the Lorentz invariance of dµ m i (q i ). The Born phase space measure is clearly singled out, which leads to the following algorithm for the construction of the full phase space 1. generate a Born configuration; 2. generate unresolved momenta subject to constraints; 3. determine the rescaling parameters and, by the same, the full reference momenta; 4. determine the Lorentz transformation yielding q fromq, and apply it to the final state momenta of the Born configuration; 5. multiply the weight by the Jacobian.
One advantage of this procedure is that it allows for the use of a multi-channel phase space generator for the Born configuration, which is particularly useful in case of intermediate resonances. Furthermore, electroweak decays should not affect efficiency, since intermediate invariant masses are not modified in the absence of QCD radiation from the decay products. The rescaling parameters are fixed by Eq. (2.5) and an additional constraint in the case of two references (see section 2.1.4). The relations always involve a sum (final state reference) or a difference (initial state reference) of a reference momentum and an unresolved momentum. Thus, if u = α r (collinear or soft limit), for some unresolved, u, and reference, r, momenta, then the constraints fix r ± u (the observable or initial state momentum) independently of u. This is the reason for the minimal number of subtraction kinematics in all cases but the single-unresolved configurations of the triple-collinear parameterisation (see Section 2.1.3).
In the case of initial state references, the energy and rapidity of the initial state is modified. In order to have a minimal number of configurations, it is necessary to choose a frame with constant boost w.r.t. the laboratory, e.g. the laboratory frame or the center-ofmass frame of the underlying Born configuration. On the other hand, the center-of-mass frame of the (n + n u )-configuration cannot be used as it does not satisfy this constraint.
The same (n + n u )-configuration corresponds to different Born configurations for different parameterisations (single-collinear, triple-collinear and double-collinear). Thus, if the angles and energies of the unresolved momenta are required to have the same meaning across parameterisations, it is not possible to use the center-of-mass frame of the Born configuration either. This is important, since using the same angle and energy definition for all configurations yields simpler 't Hooft-Veltman corrections. For this reason we use the laboratory system in the construction of the phase space.

One reference momentum
We explicitly treat the triple-collinear parameterisation. The single-collinear case is recovered by setting u 2 = 0. The following constraints determine x or z final state reference: initial state reference: They have the following properties 1. the Born configuration is well defined for any full configuration, i.e. x, z,q 0 ≥ 0, and x, z ≤ 1, In consequence, an iterative energy parameterisation with u 0 1 determined in the range 0, u 0 1 max , followed by u 0 2 in the range 0, u 0 2 max (u 1 ) covers the full phase space. The maxima of the energies, u 0 i max , are obtained at x = 0 or z =x. Due to the further requirement u 0 1 ≥ u 0 2 necessary to factorise double-soft limits, we introduce the parameterisation

Energy parameterisation for single-unresolved configurations
The parameterisation of the angles of the unresolved momenta in the case of a single reference momentum follows Section 4.2 of Ref. [43]. In particular, there iŝ The phase space is further decomposed in the variablesη 1 andη 2 as shown in Fig. 1, where we have merged sector 2 and 3 defined in Ref. [43] as suggested in Ref. [44]. We consider the four sectors separately and verify the number of subtraction terms in the single-unresolved configurations Sector 1 involves two independently vanishing variables (η 2 or ξ 2 ) in the single-unresolved kinematics, which corresponds to u 2 = α r. The relevant resolved momenta are u 1 and r + u 2 (r/z − u 2 ) for a final-state (initial-state) reference. u 1 is specified completely without any reference to u 2 thanks to the iterative energy parameterisation. The single-unresolved configuration is thus unique by the arguments of Section 2.1.1. Sector 2,3 involves one vanishing variable for each of the two partons in the singleunresolved kinematics (η 1 for u 1 , ξ 2 for u 2 ) implying that there are two singleunresolved configurations.
Sectors 4 and 5 involve two independently vanishing variables (η 2 for sector 4, η 1 for sector 5, or ξ 2 in both sectors) in the single-unresolved kinematics, which corresponds to u 2 = α u 1 . The relevant resolved momenta are r (if in the final state) and u 1 + u 2 .
In the iterative energy parameterisation, the energy of u 1 + u 2 depends on the energy of u 2 . For this reason the single-unresolved configuration is not unique.
To make the single-unresolved configuration unique, we introduce an alternative energy parameterisation in terms of the sum of the energies and their relative proportion . (2.13) While ξ 2 ≤ 2, restricting its variation range to ξ 2 ∈ [0, 1] implies u 0 1 ≥ u 0 2 . For any value of ξ 2 , x and z are monotonically decreasing functions of u 0 12 . Let -9 -For a final state reference there is while for an initial state reference there is (2.16) With the energies parameterised as follows the integration measure is If u 2 = α u 1 , u 0 12 max does not depend on u 0 2 . Thus ξ 1 uniquely determines the resolved momentum u 1 +u 2 . At the same time the rescaling of the reference momentum is also unique. In consequence, we have only one single-unresolved configuration in sector 4. There is also only one single-unresolved configuration at η 1 = 0 independently of ξ 2 in sector 5. However, there is a second single-unresolved configuration at ξ 2 = 0, η 1 = 0 as discussed below.
The remapping of the energy variables can be introduced into the parameterisation of the phase space from the start, since the sectors 4 and 5 are defined independently of the energy of the partons. In principle, both energy parameterisations may be used in the sector 2,3. However, sector 1 requires the original parameterisation in order not to introduce additional subtraction kinematics.
Sector 5. At this point, the configurations corresponding to η 1 = 0 and ξ 2 = 0, η 1 = 0 are different. This is due to the fact that the direction of the soft momentum u 2 in the latter case influences the direction of u 1 which is resolved. Thus, there is a second single-unresolved configuration.
By the above considerations, we have one single-unresolved configuration in sectors 1 and 4, and two single-unresolved configurations in sectors 2,3 and 5.

Two reference momenta
The following constraints allow to determine x 1,2 , z 1,2 (the classification corresponds to the position of the references r 1 and r 2 in that order, p is the second initial state momentum) final-final: -10 - final-initial: initial-final: initial-initial: They have the following properties 1. the Born configuration is well defined for any full configuration, i.e. x 1,2 , z 1,2 ,q 0 ≥ 0, and x 1 ≤ 1, z 1,2 ≤ 1, 2.
x 1 , z 1 are monotonically decreasing functions of u 0 1 independent of u 2 , 3. x 2 , z 2 are monotonically decreasing functions of u 0 2 at fixed u 1 .
In consequence, an iterative energy parameterisation with u 0 1 determined in the range 0, u 0 1 max , followed by u 0 2 in the range 0, u 0 2 max (u 1 ) covers the full phase space. The maxima of the energies, u 0 i max , are obtained at x i = 0 or z i =x i . Imposing an ordering of the energies, u 0 1 ≥ u 0 2 , for a given phase space parameterisation is compensated by adding the contribution of the phase space for the parameterisation corresponding to swapped references. This covers the full phase space if and only if the condition u 0 1 ≥ u 0 2 is applied in the same frame in both cases. Unfortunately, if at least one of the references is in the initial state, the chosen parameterisations lead to different Born frames upon swapping the references. Therefore, the ordering of the energies cannot be imposed in the Born frame. Instead, we apply the parameterisations in a fixed frame with respect to the lab. The full phase space is then correctly covered with the energy parameterisation The parameterisation of the angles of the unresolved momenta follows Section 4.3 of Ref. [43].

Reduction to four dimensions
The construction of local subtraction terms following the strategy of sector decomposition (see Ref. [43] for details) yields integrable expressions in CDR. The different cross section contributions are Laurent-series expansions with poles in (for the CDR parameter d = 4 − 2 ) whose sum is finite due to the finiteness of the (next-to-)next-to-leading order cross section. The construction in d dimensions is straightforward but computationally cumbersome due to: 1) the necessity of including higher order -expansion terms of matrix elements; 2) the growth with multiplicity of the number of effective dimensions of the phase space parameterisation. In Ref. [43], it has been shown that a four-dimensional formulation of the subtraction scheme can be given by introducing additional corrections such that the single-unresolved (SU ) and double-unresolved (DU ) contributions to the next-to-next-toleading order cross section are separately finite. In order to determine these corrections, -12 -separately finite sets of contributions must be identified. For the sake of completeness we first review the necessary notation. The hadronic cross section is given by the collinear factorisation expression where P 1,2 are the momenta of the hadrons h 1,2 , while p 1,2 = x 1,2 P 1,2 are the momenta of the partons.
Here, the superscript "R" denotes emission of an additional parton w.r.t. to the Born final state, "V " denotes a virtual-loop integration, while "C" a convolution with Altarelli-Parisi splitting kernels. Precise definitions are given in Appendix B. After introduction of sectors followed by the derivation of the subtraction and integrated subtraction terms, the next-to-leading order real-emission contribution is decomposed as followsσ whereσ R F contains the (n + 1)-particle tree-level matrix elements together with appropriate subtraction terms, whileσ R U contains the respective integrated subtraction terms and nparticle tree-level matrix elements.
In general, infrared divergences can be factorised from virtual amplitudes as follows This decomposition translates into a decomposition of the virtual contribution at next-to-leading order whereσ V F contains the n-particle one-loop finite remainders, whileσ V U contains the nparticle tree-level matrix elements of the Z (1) operator.
In consequence, there are three separately finite contributions in the next-to-leading order caseσ In each separately finite contribution, it is allowed to take the → 0 limit by removing higher-order terms in the -expansion of the matrix elements and reducing the dimension of the resolved momenta to four. This is the essence of the four-dimensional formulation of the subtraction scheme. This construction can be extended to next-to-next-to-leading order, which yields the following decompositionsσ The different contributions are identified as follows.σ RR,RV,V V F contain the same multiplicity and number-of-loops finite-remainder matrix elements asσ RR,RV,V V together with appropriate subtraction terms to make them integrable (unnecessary forσ V V ). The subscript "F R" (Finite Remainder) denotes the remaining contributions containing at most n-particle one-loop finite-remainder matrix elements. The subscript "SU " (Single Unresolved) denotes the remaining contributions containing at most (n + 1)-particle tree-level matrix elements together with appropriate subtraction terms to make them integrable. Finally, the subscript "DU " (Double Unresolved) denotes the remaining contributions containing only n-particle tree-level matrix elements. The F R, SU and DU contributions contain explicit poles in due to Z (1,2) -operator insertions and integration over subtraction terms of the F -contributions (and SU -contributions in the case of DU -contributions).
By construction, three contributions are separately finitê The finiteness of the next-to-next-to-leading order cross section implies that (2.38) Following the argument of Ref. [43],σ F R is separately finite due to the finiteness of the next-to-leading order cross section. Indeed,σ F R is given by the same expressions asσ U once tree-level amplitudes are replaced by one-loop finite remainders in the latter. This leavesσ DU +σ SU to be finite. Unfortunately, it turns out thatσ DU andσ SU are both separately divergent despite having different multiplicity resolved final states. A four-dimensional formulation of the subtraction scheme is only obtained under the assumption that aσ HV , the 't Hooft-Veltman corrections contribution, linear in F n exists such that are separately finite. An appropriateσ HV has been constructed in Ref. [43]. Here, we present a different construction which exploits the fact thatσ SU is finite for a next-toleading order measurement function, i.e. for F n = 0. The approach relies on the idea that, since the cut on the additional radiation in F n+1 is arbitrary, the divergences inσ SU with a next-to-next-to-leading order measurement function, i.e. for F n = 0, can, in fact, be described with n-particle kinematics and matrix elements. It should thus be possible to systematically identify them and subsequently shift them toσ DU .

't Hooft-Veltman corrections
The replacement of the next-to-next-to-leading order measurement function with a nextto-leading order measurement function, i.e. one with F n+1 = 0 and F n = 0, turnsσ SU into theσ U of an (n + 1)-particle next-to-leading order cross section witĥ At next-to-leading order, the contributionsσ R U ,σ V U andσ C may be written in the following formσ Notice that I R n+1 is given by a sum of contributions of individual sectors, which makes the factorisation of the d-dimensional phase space measure dΦ (d) n+1 non-trivial. With the phase space parameterisations of Section 2.1, this factorisation can, nevertheless, be achieved explicitly in each sector.
Let us now consider the structure of the -expansions of the integrands I c n+1 , c ∈ {R, V, C}, while keeping the exact -dependence of the occurring matrix elements. We have where each I c(i) n+1 is a function of through the -dependence of the tree-level matrix elements only. The simple structure of I V n+1 is due to the fact that it is given by matrix elements of the Z (1) -operator which we chose to be defined in the MS scheme, where it consists of pure poles. Even though collinear factorisation is also performed in the MS scheme, I C n+1 does have a non-trivial -expansion for µ R = µ F because of the expansion of the pre-factor (µ 2 R /µ 2 F ) in (B.3). The finiteness of the next-to-leading order cross section implies This analysis translates directly toσ SU with an appropriate change of superscripts.
Parameterised measurement function. Let us introduce a family of measurement functions F α m , m ∈ {n, n + 1, n + 2} with the following properties • F α m is infrared safe; • F α =0 n = 0 and F α=0 n = 0.
Hence, α = 0 corresponds to a next-to-leading order calculation within a next-to-nextto-leading order setup, while α = 0 corresponds to the general next-to-next-to-leading order calculation. Since we only consider single-and double-unresolved contributions, it is not necessary to define F α n+2 . In order to identify the 't Hooft-Veltman corrections, a particularly useful realisation is given by with F α n obtained from F α n+1 by taking soft and/or collinear limits. α η and α ξ are a set of global infrared sensitive variables -16 - where θ ij is the angle between two parton momenta, and p 0 i is a parton energy. The variable α η measures the minimal angle between any two partons i and j, while α ξ measures the minimal energy of the partons with respect to some arbitrary energy scale E norm . For α = 0, F α m is a well-defined next-to-leading order measurement function. For α = 0, it corresponds to the original next-to-next-to-leading order measurement function.
Identification of 't Hooft-Veltman corrections. Using the parameterised measurement function, the next-to-next-to-leading order version of Eq. (2.48) takes the form with c ∈ {RR, RV, C1}. Considering the full next-to-next-to-leading order case, we can schematically write the contributionsσ c SU in the following form represent the subtraction terms that regulate the n-particle limit of I c n+1 . Here, following the discussion of the next-to-leading order case, I c(i) n+1 contain the unexpanded (n + 1)-particle matrix elements. Hence, I c(i) n consist of the appropriate unexpanded factorisation formulae for (n + 1)-particle matrix elements in the single-soft and collinear limits.
Consider now the differenceσ c SU − I c . By reshuffling terms and neglecting O( ) contributions, it can be written aŝ The integral C c neither depends on the parameter α nor contains poles in . The integrand of Z c (α) is integrable while the phase space volume is restricted by 1 − θ α . The phase space volume thus vanishes in the α → 0 limit and so does Z c (α). Finally, the phase space integral in contains integrations over the angle and energy variables which might give rise to singularities regulated by α. In particular, after sector decomposition, the only singularities in a given phase space sector are due to soft and collinear limits of the unresolved parton momenta. In consequence, we can safely take the limit α → 0, if neither α η nor α ξ correspond to unresolved partons. Hence, the general contribution to N c (α) contains terms regular at α → 0 and integrals of the form where x is one of η 1 , η 2 , ξ 1 , ξ 2 . Expansion in yields the power-log series 58) where N c k (α) are regular at α → 0. The modified SU contributions are used as followŝ The left-hand side is independent of α and, therefore, the right-hand side has to be independent as well. Since Z c (α) are regular functions of α which vanish in the limit α → 0, the logarithms appearing in (2.58) have to cancel across the different contributions c. In the limit α → 0, we thus find that the poles that do not cancel withinσ SU are given by Thus, subtractingσ HV fromσ SU yields a finite quantity where all poles cancel. Since all terms inσ HV are proportional to F n ,σ HV can be added toσ DU . By the finiteness of the next-to-next-to-leading order cross section it follows that are separately finite. The formal manipulations of the different contributions must be performed in d dimensions. However, after this procedure, the 't Hooft-Veltman regularisation discussed in Section 8 of Ref. [43] can be applied yielding the desired four-dimensional formulation of the subtraction scheme. A collection of the required 't Hooft-Veltman corrections can be found in Appendix C.

Implementation and tests
We have implemented the complete subtraction scheme in the C++ library Stripper. In principle, the library provides sufficient functionality to evaluate NNLO QCD corrections to any process in the Standard Model. In practice, it requires appropriate matrix elements at tree-level (including up to double correlations in color and/or spin), one-loop level (including single correlations in color or spin), and two-loop level. Tree-level matrix elements for arbitrary Standard Model processes are provided by default by the Fortran library [56] introduced in Ref. [57]. The code generates amplitudes on-the-fly for arbitrary polarisations (helicities) and color configurations of external states and evaluates them numerically in double precision. We note that, while completely general, this is slower than dedicated analytic expressions for parton-scattering processes at low multiplicity. The fivepoint one-loop matrix element values that were required for our computation of jet rates have been obtained using the NJet C++ library presented in Refs. [58,59]. However, the general numerical-unitarity algorithm implemented in the library turned out to be too slow and unstable for our purposes. Instead, we have used the analytic formulae [60][61][62] for five-parton matrix elements implemented in NJet. The two-loop amplitudes have been taken from Ref. [63], which is based on Refs. [64][65][66][67][68].
The correctness of the results obtained with Stripper depends on the correctness of the matrix elements, splitting and soft functions, d-dimensional and four-dimensional phase spaces, and, finally, the 't Hooft-Veltman corrections. Apart from matrix elements, the majority of these contributions is involved in the evaluation of the NNLO QCD corrections to hadronic top-quark pair production. With the new implementation, we were able to reproduce the results of Refs. [69,70], which have recently been confirmed in Ref. [71]. On the other hand, our most recent results [72] for this process including Narrow-Width-Approximation top-quark decays in the di-lepton channel have only been obtained with the new version of Stripper.
A final test of the phase space implementation and the 't Hooft-Veltman corrections has been performed by calculating the single-jet inclusive double-differential (p T , |y|) cross section in the pure-gluon channel for proton-proton collisions at 7 TeV center-of-mass energy, jets defined with the anti-k T algorithm with R = 0.7 and with the scale µ R = µ F = µ = p T,1 /2, p T,1 , 2p T,1 for the MMHT2014nnlo68cl PDF set. In this case, we have found perfect agreement within the statistical errors estimated at below 1% with results obtained with NNLOjet (private communication). This test covers all aspects of our software necessary for the evaluation of the cross sections in the remaining channels involving quarks.

Single-jet inclusive rates for LHC @ 1TeV
In this section, we present our results for the single-jet inclusive differential distributions in the jet transverse momentum (p T ) in several jet rapidity (|y|) slices. We assume an initial state corresponding to proton-proton collisions at 13 TeV center-of-mass energy and use the central PDF4LHC15 nnlo PDF set to obtain the parton densities. The strong coupling constant running corresponds to this PDF set as well. Jets are identified with the anti-k T jet algorithm with R = 0.7, p T > 114 GeV and |y| < 4.7. Every jet identified in a given final-state parton configuration (event) is input into the appropriate rapidity-slice histogram with a weight that corresponds to a cross section contribution evaluated with the scale µ R = µ F = µ ∈ {p T , 2p T , 4p T }. The weight corresponding to µ = 2p T is taken as the central value of the prediction, while the remaining values are used to estimate the scale uncertainty. In Fig. 2  The scale uncertainty of the latter is also shown. We do not include non-perturbative or electroweak corrections. It is expected that the non-perturbative effects are smaller for R = 0.4 than for R = 0.7. However, we are not interested in having the most complete prediction, and thus choose somewhat arbitrarily one of the R values. Furthermore, with the current experimental precision, neither non-perturbative nor electroweak effects are necessary to obtain a good description of the data. In Fig. 3, we compare our results with those obtained with NNLOjet as presented in Ref. [23]. Since no numerical values are given in the latter publication, we superimpose our values, including the estimated Monte Carlo integration error as listed in Appendix D, on the respective plot from Ref. [23]. The results appear to be compatible within their respective errors 2 . The largest significant differences are observed in the first rapidity slice at low jet transverse momenta. However, this is the phase space region, where pure-gluon contributions dominate. The latter have been compared separately (see Section 2.3) and agree within one percent. We also note that even though the bulk of the events are in the low-p T /central-rapidity region, our calculation is not optimised to yield very small integration errors there. More interesting is the comparison of the results for higher p T and in the first four rapidity slices (|y| < 2.0). There, our calculation has an estimated integration error at the level of about one percent, and is still compatible with the NNLOjet result. This implies that sub-leading color effects missing in the contributions of channels involving quarks in NNLOjet are indeed at most at this level. The integration errors of our calculation in the fifth rapidity slice (2.0 < |y| < 2.5) are still less than about five percent and the two calculations remain compatible, although NNLOjet results are clearly more precise. While the integration errors in the sixth rapidity slice (2.5 < |y| < 3.0) remain below ten percent, the results can hardly be used as a precise indicator of sub-leading color effects. We also provide the outcome of our calculation in the last rapidity slice (3.2 < |y| < 4.7) to illustrate the limits of reasonable convergence within our setup. Let us finally comment on the convergence of the Monte Carlo integration in Stripper for this process. The results presented in Figs. 2 and 3 required about 350.000 CPU hours. In particular,σ RR F was evaluated with 200.000 CPU hours,σ SU with 100.000 CPU hours andσ DU with 50.000 CPU hours. A further improvement of the integration errors would require doubling the evaluation time forσ SU . It is important to note that the computation cost of the integrated subtraction terms present inσ DU andσ SU is still less than that of the pure real radiation. Hence, even if one could reduce it substantially by performing analytic integrations of the subtraction terms, the calculation would be at most twice faster for the same quality of the results. To put the performance into perspective, we point out that results for typical top-quark distributions as published recently require less than a twentieth of the quoted running times. -22 - In the present publication, we have performed a first independent and also the first complete calculation of single-jet inclusive rates for LHC @ 13 TeV with NNLO QCD accuracy. After comparing with results obtained with NNLOjet, we concluded that the sub-leading color effects not included in NNLOjet are negligible for phenomenological applications of the studied observable. One obvious extension of the present work is to evaluate fastNLO and/or APPLGRID tables for all measured jet observables in order to allow for inclusion of the experimental data in PDF fits. In view of the current computational costs, this requires either a substantial improvement of the efficiency of Stripper or the acquisition of substantial computational resources. Apart from the calculation of jet rates, we have also discussed a further evolution of the sector-improved residue subtraction scheme which allows for a minimal number of subtraction terms per phase space point. This approach does improve the convergence of cross sections, but to quantify the improvement requires further studies. There are still several avenues to explore in order to optimise the subtraction scheme. They range from pure Monte Carlo techniques such as better sampling of the initial state, through speedups of matrix element evaluation by using analytic formulae, and finishing with further modifications of the phase space treatment. We hope to be able to substantially reduce the cost of calculations in the future.

Double-real contributions
The 't Hooft-Veltman corrections to double-pole contributions contained in S 1 , S 4 , S 5 and S 6 (apart from the special case discussed below) are identical and given by (C.1) together with the corresponding subtraction terms obtained by expanding in η 1 . Here, h(η) = E norm /u 0 1,max (η) and where dµ (u i ) = i=0 dµ (i) (u i ) i is the integration measure for the unresolved parton momentum u i , S denotes the selector function, and a η 2 = 1 for S 1,6 and a η 2 = 2 for S 4,5 . The corrections to single-pole contributions depend on the sector and can be written as The single-pole contribution correction has the same structure as (C.4) with

Final remark
Through the parameterised measurement function, a dependence on the arbitrary energy scale E norm has been introduced. The final result for the next-to-next-to-leading order cross section, however, does not depend on this scale. The independence from E norm can be used either as a check on the implementation or to steer the cancellation of the arising logarithms.

D Single-jet inclusive NNLO QCD K-factors
The following tables correspond to K-factors depicted in Figs. 2 and 3.