NNLO final-state quark-pair corrections in four dimensions

We describe how NNLO final state quark-pair corrections are computed in FDR by directly enforcing gauge invariance and unitarity in the definition of the regularized divergent integrals. We give details of our approach and show how virtual and real contributions can be merged together without relying, explicitly or implicitly, on dimensional regularization. As an example, we recompute the H→bb¯+jets\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${ H} \rightarrow b {\bar{b}} + {jets}$$\end{document} and γ∗→jets\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma ^*\rightarrow {jets}$$\end{document} inclusive rates at the NNLO accuracy in the large NF\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_F$$\end{document} limit of QCD. This demonstrates, for the first time, that physical results with intermediate infrared singularities can be obtained at NNLO using a fully four-dimensional procedure.


Introduction
It is well known that calculations of observables in quantum field theory (QFT) are complicated by divergences in intermediate steps. On general grounds these divergences fall into two categories -ultraviolet (UV) divergences, associated with short wavelengths, and infrared (IR) divergences, associated with large wavelengths and/or collinear configurations. They show up in the form of integrals which do not exist in the four-dimensional physical space.
The customary way to deal with UV infinities is a two step procedure. Firstly, one regularizes the divergent integrals. Secondly, one reabsorbs the UV pieces into the free parameters of the Lagrangian. As for the IR infinities, they cancel when considering sufficiently inclusive observables, or can be reabsorbed in the parton densities. Depending on the approach, this cancellation can be achieved before or after integration. In the latter case, IR divergent integrals also need to be regulated.
When performing calculations of observables in QFT, the exact method of regulating the UV/IR divergences is arbia e-mail: bpage@ipht.fr b e-mail: pittau@ugr.es trary. Nevertheless, this freedom is not absolute as the chosen method should not interfere with two core tenets of QFT, that are • gauge invariance; (1a) • unitarity. (1b) These general principles have concrete consequences in perturbative calculations. Gauge invariance implies a set of graphical identities (see e.g. [1]) that need to be fulfilled to all perturbative orders by the Feynman diagrams of the QFT. Unitarity, meanwhile, demands the validity of the cutting equations [2,3] corresponding to the relation for the T matrix [4]. Both (1) are essential to preserve the Ward/Slavnov-Taylor identities (WI) at the regularized level, known as the regularized quantum action principle [5].
The most commonly used technique is dimensional regularization [6,7] (DReg). DReg exploits the fact that gauge invariance and unitarity are naturally preserved as they hold for the theory in all values of the dimensionality d of the space-time. Hence, the divergent integrals are analytically evaluated in d dimensions, and the asymptotic d → 4 limit is eventually taken. By doing so, UV/IR divergences are parametrized in terms of negative powers of a Laurent expansion in (d −4). In this framework, one still has some freedom to define objects in intermediate steps. Hence, several variants of DReg exist such as conventional dimensional regularization [8], 'tHooft-Veltman [6], four-dimensional helicity [9] and dimensional reduction [10]. We refer to [11] for an exact definition of all of them.
In recent years, a considerable effort has been pursued by several groups to introduce more four-dimensional ingredients in the definition of regularization. The main motivation being an attempt to simplify analytical and numerical methods, as well as to try to consider different theoretical perspectives. This four-dimensional program has resulted in a number of methods such as the four-dimensional formulation of FDH [12], implicit regularization [13,14], fourdimensional unsubtraction [15][16][17] and FDR [18][19][20][21][22]. They are described and compared in [23]. Whilst all approaches are different, FDR is the only method that does not rely on the customary UV renormalization procedure. Instead, the result of an FDR calculation is directly a renormalized quantity.
FDR treats UV divergences by performing a subtraction, extracting from the loop integrands the divergent parts which do not contain physical information -the so-called vacua. In the case of IR finite amplitudes, the validity of the FDR strategy has been explicitly demonstrated in [22]. In the presence of IR divergent configurations, the IR regulator should not interfere with principles (1). If this is achieved, the correct physical result is obtained for IR safe quantities. At NLO, it is known how to match real and one-loop contributions in the presence of final state IR singularities [21,23], while, so far, no FDR NNLO calculation has been performed involving infrared divergent configurations. In this paper we bridge this gap and give the first example of such a computation. We describe how NNLO final state quark-pair corrections can be computed in FDR in a way that automatically respects the principles (1). In particular, we reproduce the MS results for the N F part of H → bb + jets and γ * → jets. Since DReg is never used, explicitly or implicitly, this represents, to our knowledge, the first example of a realistic fully fourdimensional NNLO calculation.
The structure of the paper is as follows. In Sect. 2 we present our definitions of the virtual and real integrals appearing when computing NNLO quark-pair corrections. In Sect. 3 we discuss the relevant renormalization issues. Sections 4 and 5 give a detailed description of the H → bb + jets and γ * → jets calculations. Finally, in Sect. 6 we summarize our findings and discuss perspectives and directions opened by the procedures introduced in this work.

Inclusive quark-pair corrections and FDR
Our aim is to compute the large N F limit of cross sections including NNLO quark-pair corrections. The relevant contributions are the Born, Virtual and Real reactions given by Fig. 1 The lowest order amplitude (a), the IR divergent final-state virtual quark-pair correction (b) and the IR divergent real component (c). The blob stands for the emission of n−1 particles. Additional IR finite corrections are created if the gluons with momenta q 1 and k 34 are emitted by an off-shell particle in the blob n represents the amplitude for the emission of n partons computed at the j th perturbative order, while dΦ m is the m-particle phase-space n is drawn in Fig. 1a, where the line with momentum p is an on-shell QCD parton and the blob denotes n-1 additional final-state particles. The NNLO amplitudes can be split into IR divergent and finite parts A (2) n := A (2) n,IR + A (2) n,F , A The infrared singular contributions are depicted in Fig. 1b, c, while the finite pieces are created when the splitting gluons with momenta q 1 and k 34 are emitted by an off-shell particle contained in the blob. Due to the presence of IR parts in (5), both σ V and σ R are IR divergent, although, as is well known, their combination in infrared safe quantities is IR finite. In addition, UV infinities are present in σ V that are renormalized away when bare parameters are determined at the perturbative order appropriate to match the NNLO accuracy. In conclusion, the fully inclusive sum is a physical quantity, although its parts are separately plagued by IR and UV divergences. Indeed, it is the sim- Fig. 2 The indicesρ andσ are external to the divergent sub-diagram disconnected form the restσρ plest case of an infrared safe observable, which must give the same result when computed in any consistent scheme used to deal with the divergences. In this section we use the four-dimensional FDR framework, and describe the procedures which allow one to compute σ NNLO . We put particular emphasis on the steps needed to cope with the simultaneous presence of IR and UV infinities at two loops, and on how to merge virtual and real components. Section 2.1 presents the steps needed to define σ V , while 2.2 deals with σ R . Section 2.3 contains an explicit example which guides the reader across both real and virtual procedures.

The NNLO definition of the virtual contribution
A generic two-loop integrand in A (2) n,IR has the form with internal loop momenta q 1 , q 2 and q 12 := q 1 + q 2 . The denominator D collects all q 2 -independent propagators where k is the number of propagators in the blob of Fig. 1b, and N is the numerator of the integrand. The hats on Lorentz indices means that they are external to the divergent subdiagram, as in Fig. 2. The difference between hatted and un-hatted indices can be ignored until equation (23). Here we just point out that hatted indices can be identified in more general cases as described in [22]. We now analyze all possible divergences generated upon loop integration. We do this to determine the form of the FDR UV subtractions and the structure of the IR regulator needed to define the FDR integration over (7) given in (35). J (q 1 , q 2 ) is quadratically divergent when q 2 → ∞ at fixed q 1 . Given the presence of many propagators in D, this is the only possible UV sub-divergence. Depending on where the lower gluon reconnects to the blob, it may also generate global UV infinities. In addition, due to the on-shell condition p 2 = 0, a double collinear IR singularity arises when both q 1 and q 2 are proportional to p. A second potential IR divergent configuration is q 2 1 → 0 but q 2 2 = 0. Nevertheless, this divergence is cancelled by the UV behavior of the q 2 integration. This can be understood by noticing that a scaleless q 2 -type sub-integral is generated in this case, that vanishes in FDR. 1 An additional double collinear IR singularity is created if the lower gluon is attached to a second external massless parton.
In FDR, a common unphysical scale μ 2 is added to all propagators in order to regulate divergences. As is well known, performing this operation only in the denominator is contrary to the core principle (1a) of gauge invariance. As such, in FDR one performs a Global Prescription (GP) [18], also making the replacement in the numerator such that all integrand cancellations between numerator and denominator take place at the regulated level. 2 These cancellations are called "gauge cancellations". In practice, one first determines the dependence of J (q 1 , q 2 ) on q 2 1 , q 2 2 and q 2 12 generated by Feynman rules. This is the reason for the explicit q 2 1 as an argument of F in (7). Subsequently, one adds an unphysical scale μ 2 to all such self-contractions We denote this procedure by the symbol → GP . Therefore, the replacement f → GPf , applied to any function f (q 2 1 , q 2 2 , q 2 12 ), produces a new functionf with arguments replaced as in (9), Note thatf has the same functional form as f . Nevertheless, we keep the bar to make explicit that GP has been implemented on f . In the case at hand one has with D :=D pq and In our calculation, we simplify theq 2 i s inN with the denominators of (11), After this is done, μ 2 only appears in the propagators. The IR singularities of J (q 1 , q 2 ) are now regulated by the addition of μ 2 in the propagators. As we shall see, the asymptotic limit μ 2 → 0 will be eventually taken after integration. That generates logarithms of μ 2 of IR origin. Nevertheless J (q 1 , q 2 ) is still UV divergent. The global UV infinities are subtracted by separating physical and non-physical scales in D p . That means using the identity and noticing that the second term is more UV convergent than the original propagator. The same expansion has to be applied to the other propagators in 1/D, untilF/D is written as follows where [F/D] V do not depend on physical scales. Since alsō G ρσ does not contain physical scales, [F/D] V defines the global UV divergent behavior ofJ (q 1 , q 2 ): [J (q 1 , q 2 )] GV is called a Global Vacuum (GV) and is written between square brackets, that is the standard FDR notation to indicate the vacuum part of an object. Note that (F/D) F in (16) gives rise to a subtracted integrand which is globally UV convergent but still divergent when q 2 → ∞: This is fixed by subtracting the Sub-Vacuum (SV) fromḠ ρσ by means of the expansion The final result has the form so that the fully UV subtracted integrand is integrable in four dimensions. 3 Upon integration, the vacua subtracted in (21) induce the appearance of logarithms of μ 2 of UV origin, so that both IR and UV singularities are regulated by the same regulator.
The procedure leading to (21) is conveniently encoded in a linear integral operator whose action on a two-loop integrand is defined by three subsequent operations: -subtract the vacua; -integrate over q 1 and q 2 ; -take the asymptotic limit μ 2 → 0.
The last operation means retaining only the logarithmic pieces in the asymptotic expansion, neglecting O(μ 2 ) terms. Thus, the FDR two-loop integration overJ (q 1 , q 2 ) in (11) is defined as follows 4 In the following, we will often omit terms that integrate to zero in globally prescribed numerators. Hence, it is convenient to introduce a notation for that, that isN N if both numerators give the same result upon FDR integration. Equation (23) defines a gauge-invariant object, in which the necessary gauge cancellations are preserved by the GP operation. By "gauge-invariant object" we mean that a calculation ofĪ in a different gauge will give the same result. It is instructive to check that no change inĪ is produced if one shifts the numerator of the gluon propagator as Thus,Ī gives the same result when computed in any gauge. Another consequence of the WIs is that the term proportional q (13) should not contribute toĪ when contracted withFρσ . That is, After taking into account the vanishing of scaleless integrals, this corresponds to the WI depicted in Fig. 3. Nevertheless, we keep this piece in (13), as we will explicitly show in our calculation that it never contributes.  (2) for H → bb + jets in the large N F limit Let us now consider the unitarity properties of this prescription. Equation (2) relates different orders of perturbation theory, so a given divergent subgraph S G can appear embedded in loop diagrams of different orders . If the result of S G depends on then the relation (2) will in general not survive. For a more concrete example relevant to the calculations of this paper, consider the consequence of (2) shown in Fig. 4. This specifies that the discontinuity of a two-loop graph is given in terms of a one-loop graph. As such, the higher loop regularization must be consistent with the lower loop one. In FDR this places strong constraints on the GP. In our example of Fig. 4 we see that on the left hand side, the momenta associated to the gluon lines take part in the GP, but this is not true for the right hand side. In FDR, this tension is resolved by enforcing "sub-integration consistency" (SIC) [22].
We now consider this procedure in detail for the integrand J (q 1 , q 2 ). We analyze a subset of the full integrand, specifically the numerator of the gluonic self-energy of Fig. 2 12 ), which appears in (7). This is the piece that needs to be treated carefully to maintain SIC. We consider, in particular, the result of the contraction with a gρσ 5 (potentially) contained in Fρσ . It reads where we have explicitly separated the contributionq 2 2 := gρσ q ρ 2 q σ 2 . From the point of view of the sub-diagram disconnected from the rest,q 2 2 should not be globally prescribed, hence the GP replacement to be performed is 5 The analysis for terms like γρ γσ is also required, but equivalent.
On the other hand, embedding this in the full diagram requires integrating over q 1 , therefore q 2 1 needs to be barred Nevertheless, ρ and σ are internal indices of the two-loop diagrams of Fig. 1b, so that GP is needed forq 2 2 as well which is the prescription used in (23). It is the mismatch among (27)-(29) which violates SIC. Our solution to this problem is modifying the integrand in (23) as follows: -we do not apply GP toq 2 2 terms whose origin is a contraction with indices external to the UV divergent subdiagram; 6 -we replace back everywhereq 2 1 → q 2 1 after GV subtraction.
Note that the last operation is possible because barring q 2 2 and q 2 12 is sufficient to regulate the IR divergences. This is a consequence of the fact that the only IR divergent configuration is the double collinear limit. Furthermore, this solution does not affect any of the WIs and so still defines a gauge invariant object. In summary, after GV subtraction, (27) is maintained as it is also when embedded in a two-loop calculation. A comparison between this solution and the IR-free case in given in Appendix A.
Let us now return to the matter of enforcing SIC in the entirety of (23), and how one applies our solution. We enforce SIC by rewriting (23) as where we have used the fact that [F/D] V is subtracted by the integral operator. The structure of the expansions needed to extract the GV is such thatD can be always pulled out from the rest. Thus, it is possible to rewrite Next, we introduce the numerator function where the explicit dependence onq 2 2 is generated by the l.h.s of (27). In practice,Z(q 2 1 ,q 2 2 ) is constructed from N in (7) as follows: -globally prescribe N , N → GPN , leavingq 2 2 unbarred; -perform GV subtraction and determine, for each term T inN , the appropriate functionH to be used in (32). The result of this will always have a factorized form. For instance, if (15) has to be used once to subtract the global vacuum, the contribution of T toZ(q 2 1 ,q 2 2 ) is −2(q 1 · p)/q 2 1 × T ; -eliminates the bars from theq 2 1 s; -identifyq 2 2 with q 2 2 .
We denote the last two operations with the symbol → SIC . Thus, the changef → SICf , applied to any globally prescribed and GV subtracted functionf (q 2 1 ,q 2 2 ), produces a new function defined as In the equation above,f has the same functional form asf , but we use the tilded notation to denote that the SIC operation has been imposed onf . Thus, the SIC compatible version of (32) reads To continue to ensure gauge cancellations, we unbar also the propagators, that leads tõ Equation (35) defines the SIC preserving four-dimensional integration over the integrand in (7). Note that, since the GV has been subtracted in it, [d 4 q 1 ] is replaced by a customary integration d 4 q 1 . The asymptotic limit μ 2 → 0 is understood after taking the two integrations.
A first consequence of this definition is that external wavefunction corrections vanish for massless particles, so that they can always be neglected in actual calculations. The proof is given in Appendix B.

The NNLO definition of the real component
Given the propagator structure of Fig. 1c, the integrands contributing to σ R in (3) have the following form where 0 ≤ α, β ≤ 2. The symbol S collects all the remaining propagators and N R is the numerator of the amplitude squared. Depending on the value of the exponents α and β, J R becomes IR divergent under integration over Φ n+2 . These IR singularities must be regulated coherently with our treatment of the virtual component, without violating unitarity and gauge invariance. In this section, we determine a fourdimensional integration that achieves this.ρ .

Fig. 5
Virtual and real cuts contributing to the IR divergent parts of σ V (a, b) and σ R (c, d). The special indices needed to enforce (27) and (38b) are denoted byρ andσ Our starting point is the representation of σ V and σ R in terms of cut diagrams, in which we put the complex conjugate amplitudes on the right side. With this convention, normal Feynman rules are assumed on the left and complex conjugate ones on the right. The cuts generating IR divergent configurations are obtained by squaring the amplitudes in Fig. 1b, c and are depicted in Fig. 5. IR divergences manifest themselves as pinch singularities of the loop integrals in (a,b) and endpoint phase-space singularities in (c,d). These two kinds of singularities are related to each other by the identity in which the poles of the propagator on the l.h.s. may create a pinch in the complex plane of the loop integration and the δ + (k 2 ) on the r.h.s. may induce a non-integrable end-point configuration. Dubbing Σ c the sum over all cuts that appears in the r.h.s. of (2), the cutting equations [2,3] ensure that the last term in (37) does not contribute to the singular part of each cut in Σ c , and that Σ c is IR finite. This theorem implies a unitarity-preserving cancellation of the IR singularities if the Cutkosky relation giving the possibility of a one-to-one integrand level identification of the infrared divergent parts contributing to different cuts in Σ c , is preserved. However, one should also prove that The reason for this second requirement is the different origin of the potential numerators multiplying the two sides of (38a). In the case of a fermion line, that is the only cut relevant for this paper, the l.h.s. gets multiplied by the numerator of the propagator Hence, in addition to (38a), the identity must hold to guarantee the validity of (2). Equation (38b) states that the numerators of the loop propagators in Fig. 5a, b should be treated consistently with the sum over the spin polarizations of the cut lines in Fig. 5c, d. Note that (38b) also guarantees consistent gauge cancellations in all terms contributing to Σ c . Let us consider how to make these relations consistent with the procedure from the previous section where we make two modifications to the integrand: -adding μ 2 to the propagators with momenta q 2 and q 12 ; -SV and GV subtraction from the integrand.
As for (38b), in order to preserve it, one must treat k 2 3 , k 2 4 and (k 3 + k 4 ) 2 = s 34 in the numerator N R of (36) using the same prescriptions imposed on q 2 2 , q 2 12 and q 2 1 in N of (7). This means replacing in N R k 2 3,4 →k 2 3,4 = 0, where the last equalities are induced by the delta functions in (39). These changes should be performed everywhere in N R except in contractions induced by the external indicesρ andσ in cuts (c,d). In this casê in accordance to the SIC preserving requirement we have used to constructZ(q 2 1 , q 2 2 ) from N . 7 We denote all of this by introducing a globally prescribed and SIC preserving version of N R where the action of → GS on a function f (k 2 in which the tilde onf expresses the fact that GS has been imposed. Only a dependence on μ 2 is left in the r.h.s. of (46) because of the deltas inΦ n+2 . The only remaining propagator involved in our definition of the virtual component is in the top-left line of cut (a), that 7 Although it does not occur in our calculation, it is worth mentioning what happens when loop propagators and cut lines refer to spin-1 particles. In this case (38b) is replaced by g μν prop = g μν spin , which states that the metric tensor in the numerator of a loop propagator should be handled consistently with the one generated by the sum over polarizations. This is again achieved by treating the self contractions k 2 3 , k 2 4 , (k 3 + k 4 ) 2 in σ R exactly as their counterparts q 2 2 , q 2 12 , q 2 1 in σ V , as done in (44) and (45). must correspond to the cut k 1 propagator in (d) through the relation 8 1 Everything is massless in this case, so that (38b) is fulfilled and it is sufficient to check that GV subtraction does not alter (48). The proof is similar to the one used for the SV. In fact, if m expansions are needed to subtract the vacuum in front of a term in Z(q 2 1 ,q 2 2 ), this term gets multiplied by a factor λ m iñ But λ = 1 when the propagator goes on-shell. So that (48) survives the GV subtraction.
In summary, we define the four-dimensional integration over the integrand in (36) as follows By doing that, unitarity preserving IR cancellations occur by construction between σ R and σ V , without violating gauge invariance.

An example of cancellation
The integrals in (35) and (51) can be computed independently, and this is the strategy we adopt in this paper. In fact, term by term cancellations in Σ c are difficult to find. The reason is that one can add to the numerator of the virtual piece arbitrary vanishing terms that nevertheless contribute to the real part, and vice-versa. This is due to the different structure of the deltas contained in dΦ n and dΦ n+2 . However, if the numerator of a term does not change -modulo a relabelling of the momenta -when multiplied by both phase-spaces, the IR cancellation must occur between integrals constructed one from the other via the replacement in (38a). In this section we illustrate this phenomenon by considering a piece of the full H → bb + jets calculation presented in the following section. This also serves us as a concrete example on how the procedures of Sects. 2.1 and 2.2 work in practice.
The two-loop contribution to σ V is depicted in Fig. 6. We focus on the term generated by the g ρσ piece of the internal trace. It reads The complex conjugate of (48) links cuts (b) and (c).
with D i := (q 1 + p i ) 2 . The numeratorN V is obtained via GP from its unbarred version where we have neglected couplings and color factors, but not phases needed to compute the overall sign. We rewrite with s = P 2 , that gives This expression can be simplified by dropping terms which integrate to zero, for exampleq 2 2 andq 2 12 generate vacua and D 1 gives a scale-less integral. Furthermore, (P · q 1 ) does not contribute because it is antisymmetric when p 1 ↔ p 2 , while the denominator in (52) is symmetric.
The GV in (52) is fully removed by the subtraction of the scale-less integral. Thus,Z V (q 2 1 ) =N V (q 2 1 ) and so that the physically relevant two-loop contribution reads Γ V develops IR divergences in the form of powers of When splitting the result of the integration in a part which collects all terms containing powers of L, dubbed logarithmic part (L .P.), plus a remainder, one finds These logarithms are cancelled by a term contributing to the four-particle cut-diagram in Fig. 7 Γ where the IR behavior is now regulated by the two external massive lines k 2 3 = k 2 4 = μ 2 .Ñ R is obtained by applying the GS operation defined in (47) to the numerator of the diagram We now can check that the integrand in (61) is the correct object to cancel the logarithms in (60). In fact, it can be obtained from (58) by means of the Cutkosky replacement in (38a), together with the relabellings inferred by comparing Figs. 6 and 7, and the substitution which is necessary because of the gluon propagator appearing on the r.h.s. of Fig. 7. Therefore,Γ R must contain a contribution with the same singular behavior of −Γ V . We dubΓ R such a contribution, andÑ R its numerator function. The terms proportional to s 234 or s 134 in (63) give zero when evaluated at the two-particle cut: they cannot be "seen" byΓ V . This leads us to the conclusion that which corresponds toZ V (q 2 1 ) in (57), modulo the first replacement in (64). Thus,Γ R is obtained by replacing N R →Ñ R in (61). An explicit calculation confirms that

Renormalization
In this section, we discuss and implement the FDR renormalization program in the context of our calculation. To do this we need to distinguish, at least conceptually, between UV regulator, IR regulator and renormalization scale. We denote them by μ UV , μ IR and μ R , respectively. In the case of IR free observables, the FDR integral operator in (22) subtracts the UV infinities before integration. For this reason, after taking the asymptotic limit μ UV → 0, μ UV can be directly interpreted as the finite renormalization scale μ R [18]. In this sense, FDR directly produces a finite, renormalized result for the loop part: nothing needs to be subtracted from it. However, this result is arbitrary until bare parameters are fixed by experimental measurements. After doing so, if the theory is renormalizable, the scale μ R gets replaced by physical scales, leading to an unambiguous prediction. This can be understood as a finite renormalization necessary to make the theory predictive [25].
In the presence of IR divergences, no distinction is made in the virtual component between μ UV and μ IR . As a matter of fact, the procedure in Sect. 2.1 assumes μ = μ UV = μ IR , preventing one from setting μ = μ UV = μ R , as is possible in the IR free case. Our solution is fixing the bare parameters in terms of physical quantities before combining virtual and real components. After this is done, the μ R scales get automatically replaced by physical scales, hence the left over μs are the μ IR s which cancel the IR behavior of the real counterpart.
The bare parameters in our calculation are α 0 S and the Yukawa coupling y 0 b . In order to implement our renormalization program we need relations linking them to measured quantities at the appropriate perturbative order, which is one loop for α 0 S and two loops for y 0 b . We are interested in corrections proportional to N F . Hence, α 0 S can be linked to the customary α MS S (s) by using the fact that the N F contribution to the running coincides in FDR and MS [21]. As a consequence, we choose our renormalized strong coupling constant to be α S = α MS S (s), that gives the relation with where L is given in (59). The Yukawa coupling is renormalized by using its proportionality to the bottom mass. The corrected bottom propagator at the pole is proportional to where m 0 is the bare mass and the Σ ( j) s are computed in Appendix C. This gives a relation between m 0 and the pole mass m which translates into with where Equation (72) contains the bare QCD coupling. Inserting (68) gives the desired two-loop relation between bare and renormalized Yukawa coupling

H → bb + j ets
In this section, we use FDR to reproduce the physical prediction for the inclusive decay width of the Higgs into two b jets up to the NNLO accuracy in the large N F limit of QCD. That means computing the observable where Γ (0) 2 (y b ) is the tree-level decay width and δΓ N F collects all the NNLO terms proportional to α 2 S N F . The correction factor δΓ N F receives contributions from processes with up to four final-state particles, namely: -H → bb up to two loops; -H → bbg at the tree level; -H → bbqq at the tree level. The tree-and one-loop two-and three-body decays in the above list contribute to δΓ N F through renormalization. As a matter of fact, due to the scalar nature of the Yukawa coupling, Γ NNLO (y b ) is a simple process in terms of the contributing tensor structures. Nevertheless, it requires the two-loop renormalization of (75).
In the following, we compute all components in the massless limit of QCD, namely with m = 0 only in y b . Our notation is as follows. We dub V In this paper we focus on the new aspects of FDR at NNLO, namely the procedures presented in Sects. 2.1 and 2.2. For this reason we do not go into detail of the calculation of the NLO part. The corresponding expressions can be computed as described in references [21,23]. However, we emphasize that FDR NLO formulae stay the same also when they contribute to a NNLO calculation. The same holds true for LO expressions multiplying higher order corrections. This is in contrast to d-dimensional regularization methods, in which higher powers in the (d − 4) expansion must be added.

H → bb up to two loops
In our conventions, the lowest order H → bb vertex is where k and l are the color indices of the bottom quarks.
in which N C is the number of colors. The one-loop correction is depicted in Fig. 8. Following reference [21] one obtains with The globally prescribed integral needed to compute the two-loop correction is given by left part of Fig. 6. It reads withD i given right after (52). Here we have replaced bare quantities with renormalized ones, because the difference is O(α 3 S ). In the following, we computeN A,B starting form their unbarred counterparts N A,B , which can be obtained from (7) by taking Fρσ = γρ(/ q 1 + / p 1 )(/ q 1 + / p 2 )γσ . By denoting / q 2 := γρq ρ 2 = γσ q σ 2 one finds where the last term originates from the piece we have studied in detail in Sect. 2.3, with N V given in (53). N A does not contribute toV (2) 2 . In fact, aū( p 1 )(u( p 2 )) is understood on the l.h.s.(r.h.s.), so that, by virtue of the Dirac equation, one can replace N A → N A = D 1 (/ q 1 + / p 2 )/ q 2 + / q 2 (/ q 1 + / p 1 )D 2 , giving which generates scale-less integrals. That explicitly proves the WI in (25). As for the N B piece, there are several ways [19,20] to deal with strings of γ -matrices to extract the dependence on (q i · q j ) and q 2 i needed to implement GP. They are based on replacements of the type / q i → / q i − μ i , where the "masses" μ i serve as a bookkeeping tool. In this paper we find it more convenient to use Clifford algebra until we reach the configurations 9 9 Note the invariance of the last term under / Fig. 9 The LO Hbbg vertex By using this method one finds N B does not induce the appearance of global UV divergences in (81), hence the numerator function directly reads Z(q 2 1 ,q 2 2 ) =N B (q 2 1 ,q 2 2 ). Thus, the SIC preserving numerator function isZ(q 2 1 , q 2 2 ) =N B (q 2 1 , q 2 2 ), and the two-loop correction is In terms of the master integrals listed in Appendix D it reads

H → bbg and H → bbqq at the tree level
The LO H → bbg decay width is obtained by squaring the V (0) 3 vertex drawn in Fig. 9 and integrating over a phasespace in which all final-state particles acquire a small mass μ, as described in reference [21]. The result is with L defined in (59).
As for H → bbqq, two diagrams contribute to the amplitude V (0) 4 . They can be read from Fig. 9 by allowing the gluon to split into a qq pair. As described in Sect. 2.2, Γ (0) 4 is obtained by squaring V (0) 4 and integrating over a massive 4-particle phase-spaceΦ 4 such that k 2 1 = k 2 2 = 0 and k 2 3 = k 2 4 = μ 2 . Prior to integration, the integrand should be modified according to the GP and SIC replacements given in (46). As a result of this, the function to be integrated in (51) is a rational combination of the invariants s 34 , s 134 , s 234 and μ 2 : where the μ 2 dependence is induced by (45). It is interesting to note this μ 2 dependence factorizes. In particular, one finds with Therefore, (51) has the form of a massless-like matrix ele-mentJ R (s 34 , s 134 , s 234 ), obeying all relevant WIs, times a global factor w(μ 2 ) integrated over a massive phase-space. We conjecture that a similar gauge-invariant structure should naturally appear also when dealing with more complicated cases involving, for instance, IR divergences induced by vector particles like gluons.
In terms of the integrals reported in Appendix E the result reads

The large N F limit of the inclusive width
Here we gather all the calculated components and compute Γ NNLO (y b ) in (76). The correction factor δΓ N F receives contributions from processes with up to four partons that are obtained by inserting the renormalization equations (68) and (75) in the amplitudes given in the previous sections. One finds where Equations (94) are written in a form that highlights the contributions generated by renormalization. Collecting all the pieces gives the IR finite result Equation (96) is written in terms of the pole mass m. It is possible to reabsorb the large logarithms of the ratio m 2 /s in a new Yukawa coupling y MS b defined through the known two-loop relation between m and the MS mass [26]. Using the N F part of it gives hence Equation (98) coincides with the known MS result [27,28].

γ * → j ets
In this section we compute the large N F limit of the inclusive e + e − → γ * → jets production rate up to the NNLO accuracy. That is the observable where σ (0) 2 is the tree-level e + e − → γ * → qq crosssection and δσ N F contains the QCD corrections proportional to α 2 S N F . QCD renormalization only involves α S , in this case. Nevertheless, higher rank tensors contribute, so that preserving gauge cancellations and unitarity in such an environment provides a more stringent test for our procedures. In this respect, γ * → jets is complementary to H → bb + jets.
The processes which contribute to δσ N F are e + e − → qq up to two loops; e + e − → qqg at the tree level; e + e − → qqq q at the tree level, where we understand a photon mediating the reactions. We dub V ( j)β i the final-state current producing i partons computed at the j th QCD order, where β is the Lorentz index of the virtual photon. The Feynman diagrams representing the vertices are obtained from those in the previous section by replacing the Higgs with a photon. Hence, we do not draw them. Contracting V ( j)β i with the initial-state current, squaring and integrating over the phase-space gives the corresponding cross section, denoted by σ ( j) i . In the following, we describe the FDR computation of the various components.

e + e − → qq up to two loops
The lowest order vertex is where Q q is the electric charge of the quark. The corresponding cross section reads in which α is the fine-structure constant. The computation of V (1)β 2 is described in [23]. The result is with a 0 and L defined in (69) and (80), respectively. The globally prescribed two-loop integral we need to compute V The unbarred N β is obtained from (7) with Fρσ = γρ(/ q 1 + / p 1 )γ β (/ q 1 + / p 2 )γσ by using the fact that, according to the WI in (25), the term proportional to q ρ 2 q σ 1 + q ρ 1 q σ 2 in the fermion trace does not contribute. 10 Hence Using tensor decomposition and the p 1 ↔ p 2 symmetry where the factor −2(q 1 · p 1 )/q 2 1 multiplying the last term subtracts its GV. Thus, the SIC preserving numerator function 10 The proof is analogous to the one given in Sect. 4.1.
In terms of the two-loop integrals in Appendix D one finds 5.2 e + e − → qqg and e + e − → qqq q at the tree level A NLO computation produces with L is given in (59). As for σ (0) 4 , it is obtained by computing the amplitude squared, modifying it according to the prescription in (46) and integrating over theΦ 4 phase-space. In terms of the integrals in Appendix E the result reads

The large N F limit of the inclusive jet production rate
Here we collect all components needed to compute σ NNLO . The correction can be split as follows where the various contributions are obtained by inserting (68) in the results of the previous sections. One has with δV (1) Gathering all the pieces gives which reproduces the MS result [29,30].

Conclusion and outlook
In this paper we have demonstrated that a fully fourdimensional framework to compute NNLO quark-pair corrections can be constructed based on the requirement of preserving the two principles given in (1). The FDR idea of enforcing gauge invariance and unitarity at the level of the UV subtracted integrands is at the base of the procedures we have used to define UV and IR divergent integrals. A few advantages of such an approach that have appeared in our calculation are, for the UV part -no (explicit or implicit) UV counterterms have to be included in the Lagrangian; -lower-order substructures are used in higher-order calculations without any modification (see e.g. (79) and (102)).
For the IR sector -infrared divergences in the real component directly show up in terms of logarithms of a small cut-off parameter μ IR , with no need for a prior subtraction of 1/(d − 4) poles (see, for instance, the four-parton rates in (94) and (111)); -one-to-one integrand correspondences can be written down between virtual and real contributions (see Sect. 2.3).
In this paper we have focused our attention on a special class of NNLO corrections. However, we believe that the basic principles that have guided us towards a consistent treatment of all the pieces contributing to the final NNLO answer will remain valid also when considering more complicated environments, with the final aim of constructing a completely general procedure including also initial state IR singularities. This is certainly the main subject of our future investigations. Other possible directions are: using μ IR as a separation parameter in slicing-based subtraction methods at NNLO [31], or exploiting the virtual/real integrand correspondence to construct four-dimensional local counter-terms directly from the virtuals, 11 in the same spirit of the procedure presented at NLO in the FDR section of [23].
On a more general ground, we envisage that the intrinsic four-dimensionality of FDR can pave the way to new numerical methods and that there is room for fully exploiting its potential in NNLO calculations.
which fixes the correct UV behavior. In addition, EIs admit theq 2 1 → q 2 1 limit that matches the rest of the calculation. In summary, the solution presented in Sect. 2.1 is equivalent to the following procedure: -apply GP; -subtract GV; -downgradeq 2 1 → q 2 1 in the result; -identify the EEIs to be added (using the same algorithm as in the IR-free case); -replace each EEI with the corresponding EI.
It would be interesting to establish whether this strategy works also for IR finite two-loop calculations. That would make unnecessary the use of the EEIs. We leave this to further investigations.

Appendix B: Massless wave-function corrections
Wave function corrections are generated when the lower gluon in Fig. 1b reconnects to the emitting massless parton. In this appendix, we use the results of Sect. 2.1 to demonstrate that they vanish.
The relevant integrand is obtained by taking Fρσ = γρ(/ q 1 + / p)γσ in (7), that gives N = N A + N B with Furthermore, D = q 4 1 D p , so that the integrals we have to consider arē (B.5) One finds Only the third term contributes toĪ A . All the others generate vacua or result from contractions of antisymmetric combinations of γ -matrices with symmetric integrals. Thus I A only depends on p 2 = 0. In addition, it is both UV divergent and logarithmically IR divergent, so that it is a scaleless integral. Such integrals vanish in FDR as a consequence of an exact cancellation between UV and IR singularities, thereforeĪ A = 0, as required by (25). In the same way,Ī B is fully scale-less The proof that scale-less integrals do not contribute can be found in [20]. Here we prove that they vanish also when q 2 1 is unbarred, as in (35). We concentrate on first term of (B.8) The proof is unchanged for all the other contributions. A GV subtraction is needed in front ofq 2 2 / q 1 . This is achieved by using twice the identity in (15) 1 (B.10) Theq 2 2 / p piece is less UV divergent, so that a single subtraction is sufficient The vacua are subtracted by the integral operator. That defines the numerator function associated withN C . (B.14) I C diverges logarithmically in the double collinear configuration in the absence of regulator. The barred q 2 -type denominators are sufficient to regulate this. That is a consequence of the fact that (B.10) and (B.11) do not alter the IR power counting. By tensor decompositionĨ C ∼ / p p 2 /μ 2 + O( p 4 /μ 4 ) , so that it vanishes on-shell. In summary, the GV subtraction does not leave finite pieces in scale-less integrals.

Appendix D: The virtual master integrals
In this appendix we sketch out the computation of the twoloop integrals appearing in our calculation. The q 2 integration is performed first. 12 This means computing (D. 28) The coefficients are obtained by subtracting the SV by means of (19), and integrating over the finite part. The result reads When inserting these results in Eqs. (86) and (106) one finds that the two-loop vertex corrections can be expressed in terms of three master integrals whereĨ 2 is UV finite because p 2 1 = 0. Integrating over q 1 and x and neglecting O(μ 2 ) terms gives −12ζ 3 + 5 3 π 2 + 328 9 ,  with L written in (59).