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. The whole procedure is strictly four-dimensional and renormalization is performed by simply fixing bare parameters in terms of physical measurements. We give details of our approach and demonstrate how virtual and real contributions can be merged together without relying on dimensional regularization. As an example, we recompute the ${ H} \to b \bar b + {jets}$ and $\gamma^\ast \to {jets}$ inclusive rates at the NNLO accuracy in the large $N_F$ limit of QCD.


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 arbitrary. Nevertheless, this freedom is not absolute as the chosen method should not interfere with two core tenets of QFT, that are • gauge invariance; (1.1a) • unitarity. (1.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.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], four-dimensional 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.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.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 section 2 we present our definitions of the virtual and real integrals appearing when computing NNLO quark-pair corrections. In section 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 section 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 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 in which P is the initial state momentum. The amplitude A (0) n is drawn in figure 1-(a), 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 3) The infrared singular contributions are depicted in figure 1-(b,c), while the finite pieces are created when the splitting gluon is emitted by the blob. Due to the presence of IR parts in (2.3), 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 simplest 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 1 -dependent propagators where k is the number of propagators in the blob of figure 1-(b), and N is the numerator of the integrand. The hats on Lorentz indices means that they are external to the divergent sub-diagram, as in figure 2. The difference between hatted and un-hatted indices can be ignored until equation (2.20).
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 (2.5) given in (2.32). 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 gluon is attached to a second external massless parton.
In FDR, an 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 (1.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 that 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 (2.5). Subsequently, one adds an unphysical scale µ 2 to all such self-contractions We denote this procedure by the symbol → GP . Thus, the replacement A → GPĀ , applied to 1 See appendix B. 2 This, combined with the shift invariance of FDR integrals is sufficient to prove all Ward identities graphically [20]. any function A(q 2 1 , q 2 2 , q 2 12 ), produces a new functionĀ with arguments replaced as in (2.7), A(q 2 1 ,q 2 2 ,q 2 12 ) := A(q 2 1 ,q 2 2 ,q 2 12 ). (2.8) In the case at hand one has (2.11) 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. NeverthelessJ(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 12) 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 followsFρσ [J(q 1 , q 2 )] GV := Fρσ (q 1 ,q 2 1 ) D VḠ ρσ . (2.14) [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 (2.13) gives rise to a subtracted integrand which is globally UV convergent but still divergent when q 2 → ∞: (2.15) This is fixed by subtracting the Sub-Vacuum (SV) fromḠ ρσ by means of the expansion (2.16) The final result has the formḠ ρσ = Ḡ ρσ SV + Ḡ ρσ F , (2.17) so that the fully UV subtracted integrand 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 (2.9) is defined as follows 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 (2.20) 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 ρ 2 q σ 1 + q ρ 1 q σ 2 in (2.11) 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 figure 3. Nevertheless, we keep this piece in (2.11), as we will explicitly show in our calculation that it never contributes.
Let us now consider the unitarity properties of this prescription. Equation (1.2) relates different orders of perturbation theory, so a given divergent subgraph D can appear embedded in loop diagrams of different orders ℓ. If the result of D depends on ℓ then the relation (1.2) will in general not survive. For a more concrete example relevant to the calculations of this paper, consider the consequence of (1.2) shown in figure 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 figure 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 figure 2, q ρ 2 q σ 12 + q ρ 12 q σ 2 − g ρσ (q 2 · q 12 ), which appears in (2.5). 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ρσ 3 (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 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 figure 1-(b), so that GP is needed forq 2 2 as well which is the prescription used in (2.20). It is the mismatch among (2.24), (2.25) and (2.26) which violates SIC. Our solution to this problem is modifying the integrand in (2.20) as follows: • we do not apply GP toq 2 2 terms whose origin is a contraction with indices external to the UV divergent sub-diagram; • 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, (2.24) 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 (2.20), and how one applies our solution. We enforce SIC by rewriting (2.20) 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 (2.24). In practice,Z(q 2 1 ,q 2 2 ) is constructed from N in (2.5) as follows: • globally prescribe N , N → GPN , leavingq 2 2 unbarred; • perform GV subtraction and determine, for each term T inN , the appropriate func-tionH to be used in (2.29). The result of this will always have a factorized form. For instance, if (2.12) 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; We denote the last two operations with the symbol → SIC . Thus, the changeĀ → SICÃ , applied to any globally prescribed and GV subtracted functionĀ(q 2 1 ,q 2 2 ), produces a new function defined asÃ Thus, the SIC compatible version of (2.29) reads To continue to ensure gauge cancellations, we unbar also the propagators, that leads tõ (2.32) Equation (2.32) defines the SIC preserving four-dimensional integration over the integrand in (2.5). 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 wave-function 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 figure 1-(c), the integrands contributing to σ R in (2.1) have the following form where 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 four-dimensional integration that achieves this. 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, . Figure 5. Virtual and real cuts contributing to the IR divergent parts of σ V (a,b) and σ R (c,d).
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 figure 1-(b,c) and are depicted in figure 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 (1.2), the cutting equations [2,3] ensure that the last term in (2.34) 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 Σ c = σ R + σ V . The reason for this second requirement is the different origin of the potential numerators multiplying the two sides of (2.35a). 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 / k : must hold to guarantee the validity of (1.2). Note that (2.35b) 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 a few propagators; • SV and GV subtraction from the integrand.
We first study how FDR preserves (2.35a). We start dealing with the effect of the q 2 2 → GPq 2 2 and q 2 12 → GPq 2 12 replacements in cuts (a,b). Equation (2.35a) is preserved if Thus, theq 2 2 andq 2 12 propagators in σ V must correspond to external particles in σ R obeying Hence, we replace in (2.1) Φ n+2 →Φ n+2 , where the phase-spaceΦ n+2 is such that k 2 3 = k 2 4 = µ 2 and k 2 i = 0 when i = 3, 4. However, this is not enough. One also needs to show that (2.36) survives the SV subtraction of (2.18). We prove this explicitly in the case of the last term in (2.11). The proof is unchanged for the other contributions. The relevant expansion is where SV is the term to be subtracted. We consider a piece of the finite part of (2.38) as a numerator factor f and observe that f → 1 whenq 2 12 → 0. We therefore first put theq 2 12 propagator on-shell, giving (2.40) When alsok 2 3 goes on-shell, one obtains the same result as applying (2.36) before subtracting the vacuum. Hence, the SV subtraction is "invisible" from the point of view of (2.36), and (2.35a) is fulfilled if k 3,4 obey (2.37).
As for (2.35b), in order to preserve it, one must treat / k 3 and / k 4 in the numerator N R of (2.33) using the same prescriptions imposed on / q 2 and / q 12 in N . This means replacing in N R in accordance to the SIC preserving requirement we have used to constructZ(q 2 1 , q 2 2 ). We denote all of this by introducing a globally prescribed and SIC preserving version of N R Only a dependence on µ 2 is left in the r.h.s. of (2.43) because of the deltas inΦ n+2 . The only remaining propagator modified by our definition of the virtual component is in the top-left line of cut (a), that must correspond to the cut k 1 propagator in (d) through the relation 4 Everything is massless in this case, so that (2.35b) is fulfilled and it is sufficient to check that GV subtraction does not alter (2.45). 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 inZ(q 2 1 ,q 2 2 ), this term gets multiplied by a factor f m inZ(q 2 1 , q 2 2 ), where (2.47) But f = 1 when the propagator goes on-shell. So that (2.45) survives GV subtraction.
In summary, we define the four-dimensional integration over the integrand in (2.33) as followsR (2.48) 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 (2.32) and (2.48) 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 Figure 6. A cut contributing to H → bb at NNLO. Only the term proportional to g ρσ is considered in (2.49).
both phase-spaces, the IR cancellation must occur between integrals constructed one from the other via the replacement in (2.35a). 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 sections 2.1 and 2.2 work in practice.
The two-loop contribution to σ V is depicted in figure 6. We focus on the term generated by the g ρσ piece of the internal trace. It reads , (2.49) 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 examplē q 2 2 andq 2 12 generate vacua andD 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 (2.49) is symmetric. Figure 7. A cut contributing to H → bbqq. Only the g ρσ term is considered in (2.57).
The GV in (2.49) is fully removed by the subtraction of the scale-less integral. Thus, so that the physically relevant two-loop contribution reads . (2.55) Γ V develops IR divergences in the form of powers of L = ln(µ 2 /s). 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 (2.56) These logarithms are cancelled by a term contributing to the four-particle cut-diagram in figure 7 where the IR behavior is now regulated by the two external massive lines k 2 3 = k 2 4 = µ 2 . N R is obtained by applying the GS operation defined in (2.44) to the denominator of the diagram N R = 64(k 3 · k 4 )(k 1 · k 234 )(k 2 · k 134 ). We now can check that the integrand in (2.57) is the correct object to cancel the logarithms in (2.56). In fact, it can be obtained from (2.55) by means of the Cutkosky replacement in (2.35a), together with the relabellings inferred by comparing figures 6 and 7, and the substitution which is necessary because of the gluon propagator appearing on the r.h.s. of figure 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 (2.59) give zero when evaluated at the two-particle cut: they cannot be "seen" byΓ V . This leads us to the conclusion 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 (2.19) 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 [24].
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 section 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 The Yukawa coupling is renormalized by using its proportionality to the bottom mass. The corrected bottom propagator at the pole is proportional to (3.7)

H → bb + jets
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 finalstate 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 (3.7).
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 sections 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. Squaring V (0) 2 gives the LO H → bb decay width in which N C is the number of colors. The one-loop correction is depicted in figure 8. Following reference [21] one obtains Figure 8. The one-loop QCD correction to the Hbb vertex.
The globally prescribed integral needed to compute the two-loop correction is given by left part of figure 6. It reads , (4.6) withD i written in (2.49). 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 (2.5) 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 section 2.3, with N V given in (2.50). 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 which generates scale-less integrals. That explicitly proves the WI in (2.22). 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 5 By using this method one finds N ′ B does not induce the appearance of global UV divergences in (4.6), hence the numerator function directly readsZ(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 (4.11) In terms of the master integrals listed in appendix D it reads (4.12)

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 figure 9 and integrating over a phase-space in which all final-state particles acquire a small mass µ, as described in reference [21]. The result is where L := ln(µ 2 /s). (4.14) A for H → bbqq, two diagrams contribute to the amplitude V 4 . They can be read from figure 9 by allowing the gluon to split into a qq pair. As described in section 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 (2.43). As a result of this, the function to be integrated is a rational combination of the invariants s 34 , s 134 , s 234 and µ 2 : S(s 34 , s 134 , s 234 , µ 2 ), (4.15) where the µ 2 dependence is induced by (2.42). It is interesting to note this µ 2 dependence factorizes. In particular, one finds S(s 34 , s 134 , s 234 , µ 2 ) = S ′ (s 34 , s 134 , s 234 )w(µ 2 ), with w(µ 2 ) = 1 + 2 µ 2 s 34 .
(4. 16) 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 (4.1). The correction factor δΓ N F receives contributions from processes with up to four partons that are obtained by inserting the renormalization equations (3.1) and (3.7) in the amplitudes given in the previous sections. One finds (4.20) Equations (4.19) are written in a form that highlights the contributions generated by renormalization. Collecting all the pieces gives the IR finite result Equation (4.21) is written in terms of the pole mass m. It is possible to reabsorbe 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 [25]. Using the N F part of it gives

γ * → jets
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 cross-section 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 (3.2) and (4.5), respectively. The globally prescribed two-loop integral we need to compute V The unbarred N β is obtained from (2.5) with Fρσ = γρ(/ q 1 + / p 1 )γ β (/ q 1 + / p 2 )γσ by using the fact that, according to the WI in (2.22), the term proportional to q ρ 2 q σ 1 + q ρ 1 q σ 2 in the fermion trace does not contribute 6 . Hence Using tensor decomposition and the p 1 ↔ p 2 symmetry gives N β → GPN β ≃M β with M β (q 2 1 ,q 2 2 ) = +2γ β 4(q 2 · p 1 )(q 2 · p 2 ) + 4 s (q 2 · P ) (q 1 · p 1 )(q 2 · p 2 ) − (q 1 · p 2 )(q 2 · p 1 ) where the factor −2(q 1 · p 1 )/q 2 1 multiplying the last term subtracts its GV. Thus, the SIC preserving numerator function isZ β (q 2 1 , q 2 2 ) =M β (q 2 1 , q 2 2 ), giving 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 (4.14). As for σ (0) 4 , it is obtained by computing the amplitude squared, modifying it according to the prescription in (2.43) 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 (3.1) in the results of the previous sections. One has Gathering all the pieces gives which reproduces the MS result [27].

Conclusion and outlook
In this paper we have demonstrated that a fully four-dimensional framework to compute NNLO quark-pair corrections can be constructed based on the requirement of preserving the two principles given in (1.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. (4.4) and (5.4)); • renormalization is equivalent to the process of expressing (finite) bare parameters in terms of measurable observables (e.g. (3.1) and (3.7)).

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 (4.19) and (5.13)); • one-to-one integrand correspondences can be written down between virtual and real contributions (see section 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 [28], or exploiting the virtual/real integrand correspondence to construct four-dimensional local counter-terms directly from the virtuals. 7 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.
A Sub-integration consistency with and w/o IR divergences When no IR infinities are present, the mismatch between equations (2.25) and (2.26) is cured by adding the so called extra-extra integrals (EEI) introduced in [22]. Their exact definition is not needed here. It suffices to say that terms proportional to the difference are included. They multiply UV 1/µ 2 poles and generate logarithms of µ 2 that restore the correct renormalization properties of the two-loop amplitude. Such contributions are missed by (2.26).
In the presence of IR divergences an additional complication is generated by the GP q 2 1 →q 2 1 in (2.25) and (2.26). After GV subtraction, the difference also hits 1/µ 2 poles of IR origin. This gives rise to different renormalization constants for processes with or without IR divergences, which is unacceptable. This leads to the choice of letting q 2 1 unbarred, as discussed in section 2.1. For the sake of consistency, also the EEIs part needs to be modified accordingly. The problem is that the EEIs become unregulated when unbarringq 2 1 at the integrand level. The solution to this is replacing EEIs with the difference of two ordinary FDR integrals, generated by the combination which is sometimes referred as an extra-integral (EI). One shows that EEIs and EIs share the same logarithmic content, 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 section 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.

B Massless wave-function corrections
Wave function corrections are generated when the lower gluon in figure 1-(b) reconnects to the emitting massless parton. In this appendix, we use the results of section 2.1 to demonstrate that they vanish. The relevant integrand is obtained by taking Fρσ = γρ(/ q 1 + / p)γσ in (2.5), that gives N = N A + N B with N A = / q 1 (/ q 1 + / p)/ q 2 + / q 2 (/ q 1 + / p)/ q 1 and N B = 2/ q 2 (/ q 1 + / p)/ q 2 + 2(q 2 · q 12 )(/ q 1 + / p). (B.1) Furthermore, D = q 4 1 D p , so that the integrals we have to consider arē 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 .

(B.4)
I A only depends on p 2 = 0. In addition, it is both UV divergent and logarithmically IR divergent, so that it is a scale-less integral. Such integrals vanish in FDR as a consequence of an exact cancellation between UV and IR singularities, thereforeĪ A = 0, as required by (2.22). In the same way,Ī B is fully scale-less so that self-energy correctionsĪ A +Ī B vanish. 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 (2.32). We concentrate on first term of (B.5)Ī 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 (2.12) Theq 2 2 / p piece is less UV divergent, so that a single subtraction is sufficient 1 The vacua are subtracted by the integral operator. That defines the numerator function associated withN CZ which produces with D p = q 2 1 + 2(q 1 · p). The unbarred N (/ p) is given by (2.5) with Fρσ = γρ(/ q 1 + / p + m)γσ. The result reads N (/ p) = N A (/ p) + N B (/ p), where N A (/ p) = / q 1 (/ q 1 + / p + m)/ q 2 + / q 2 (/ q 1 + / p + m)/ q 1 = 2/ q 2 D p − (/ p − m)/ q 1 / q 2 − / q 2 / q 1 (/ p − m), N B (/ p) = 2(q 2 · q 12 )(/ q 1 + / p − 2m) + 2/ q 2 (/ q 1 + / p + m)/ q 2 .

D The virtual master integrals
In this appendix we sketch out the computation of the two-loop integrals appearing in our calculation. The q 2 integration is performed first. To determine B α we use tensor decomposition . (D.4) The first two terms cancel each other due to the q 2 2 ↔ q 2 12 symmetry of the integral. Thus When inserting these results in equations (4.11) and (5.8) one finds that the two-loop vertex corrections can be expressed in terms of three master integrals L ′ 3 + 5L ′ 2 + L ′ 56 3 + π 2 − 12ζ 3 + 5 3 π 2 + 328 9 , I 2 = − 1 4 L ′ 2 + 16 3 L ′ + π 2 3 + 104 9 , with L ′ given in (4.5). Finally, the two-loop integrals in (C.9) arẽ (D.10) Their asymptotic expansions read with L written in (4.14).