The real-virtual antenna functions for $S \to Q\bar{Q} X$ at NNLO QCD

We determine, in the antenna subtraction framework for handling infrared divergences in higher order QCD calculations, the real-virtual antenna functions for processes involving the production of a pair of massive quarks by an uncolored initial state at NNLO QCD. The integrated leading and subleading color real-virtual antenna functions are computed analytically in terms of (cyclotomic) harmonic polylogarithms. As a by-product and check we compute $R_Q=\sigma(e^+e^-\to \gamma^*\to Q\bar{Q}X)/\sigma(e^+e^-\to \gamma^*\to\mu^+\mu^-)$ and compare with existing results. Our result for $R_Q$ is exact to order $\alpha_s^2$.


Introduction
In this paper we report, within the antenna subtraction framework [1][2][3][4], on the calculation of the real-virtual antenna functions for processes involving the production of a pair of massive quarks by an uncolored initial state S at next-to-next-to leading order (NNLO) QCD, where S denotes, for example, an e + e − pair or an uncolored boson. Antenna subtraction is a method for handling infrared (IR), i.e., soft and collinear divergences in higher order QCD calculations. The general features of the method at NNLO QCD were presented in [3].
Applications at NNLO QCD involving massless partons include e + e − → 2 jets, 3 jets [5][6][7][8] and pp → di-jets [9]. For QCD processes involving massive quarks the method was worked out to NLO in [10,11]. Partial results exist for NNLO QCD processes with colored initial states and massive quarks in the final state [12][13][14]. For processes of the type (1) the un-integrated and integrated NNLO real radiation subtraction terms for the QQqq and QQgg final states were determined in [15] and in [16], respectively. The NNLO real-virtual antenna functions for (1), which were missing so far, are the subject of this paper. At this point it seems appropriate to recall that other NNLO subtraction techniques exist and have been successfully applied. A method was presented in [17][18][19] which can be used for massless and massive partons and was applied in the computation of the total hadronic tt cross section to order α 4 s [20,21]. Other techniques for handling the IR divergences of the individual contributions to partonic processes at NNLO QCD include the sector decomposition algorithm [22][23][24][25][26] and the subtraction methods [27][28][29][30][31][32]. Very recently, a NNLO QCD generalization of the phase-space slicing method has been presented by [33,34] for e + e − → γ * → QQX. Coming back to the issue of this paper, we outline in the next section the construction of separately infrared (IR) finite contributions to the NNLO differential cross section of (1) from the four-parton, three-parton, and two-parton final states within the antenna subtraction method. In Sec. 3 we describe our computation of the integrated massive real-virtual antenna functions.
As a first application and check of our antenna subtraction terms we compute in Sec. 4 the contribution of order α 2 s to the ratio R Q for inclusive massive quark-pair production by e + e − annihilation via a virtual photon and compare with existing results in the literature. We conclude in Sec. 5.

The real-virtual subtraction term for reactions (1)
The order α 2 s term dσ NNLO in the strong-coupling expansion of the differential cross section of (1), dσ = dσ LO + dσ NLO + dσ NNLO , receives the following contributions: i) the double virtual correction dσ V V NNLO associated with the second-order matrix element of S → QQ (i.e., 2-loop times Born and 1-loop squared), ii) the real-virtual cross section dσ RV NNLO associated with the second-order matrix element of S → QQg (1-loop times Born), iii) the double real contribution dσ RR NNLO associated with the squared Born amplitudes S → QQgg, S → QQqq (where q denotes a massless quark), and above the 4Q theshold, S → QQQQ. The latter contribution is IR finite and is of no concern for the purpose of this section. We will come back to it in Sec. 4, where we choose the initial state in (1) to be a virtual photon, S = γ * . All the terms given below denote renormalized quantities. The ultraviolet divergences in the loop amplitudes are removed by on-shell renormalization of the external quarks and gluonsin the following, m Q denotes the on-shell mass of Q -and by MS renormalization of the strong coupling.
The terms i), ii), iii) are, apart from the QQQQ contribution, separately IR divergent. Within the subtraction method, dσ NNLO is given schematically by The subscripts Φ n denote n-particle phase-space integrals. The integrands dσ S NNLO and dσ T NNLO denote the double-real subtraction terms (for QQqq and QQgg) and the real-virtual subtraction term, respectively. The former are constructed such that the phase-space integrals where ǫ = (4 − d)/2, are finite in d = 4 dimensions in all single and double unresolved limits and can be evaluated numerically. The real-virtual subtraction term dσ T NNLO which is to be constructed such that it reproduces both the implicit singularities in single unresolved phase space regions and the explicit poles of the real-virtual cross section dσ RV NNLO , is the topic of this paper. The sum of the last three terms in (2) is also IR-finite, see below. In order to make the cancellation of IR singularities explicit in (2), the integrals of these subtraction terms, denoted by Φ 3 dσ T NNLO and Φ 4 dσ S NNLO in (2), must be computed over the phase-space regions where IR singularities arise. In the antenna subtraction formalism, which is applied here, subtraction terms are constructed from antenna functions and reduced matrix elements with remapped momenta. The antenna functions are universal building blocks and can be derived from the respective physical color-ordered squared matrix elements.

Double real-radiation corrections
The antenna subtraction terms dσ S,QQqq NNLO and dσ S,QQgg NNLO for the final states QQqq, QQgg and their integrals over appropriate phase space regions were computed in [15] and [16], respectively.
While some of these integrals, collectively denoted by Φ 4 dσ S NNLO in (2), contribute to the realvirtual subtraction term and are therefore integrated over those regions of phase space which correspond to single-unresolved emission, others are combined with the double-virtual cross section and, hence, are integrated over the unresolved phase space of two partons. In order to keep track of these connections the subtraction terms are decomposed into a sum of three contributions [4]: where the integrands of the first two terms on the r.h.s. cover, respectively, the singularities due to single and double unresolved parton configurations in dσ RR,QQgg

NNLO
, whereas the third term removes spurious single unresolved limits from dσ S,b,2,QQgg NNLO . A splitting analogous to (4) is made for the subtraction term dσ S,QQqq NNLO . Throughout this section the notation n indicates the analytic integration over the phase space of n unresolved partons in d = 4 dimensions.

Real-virtual corrections
Now, we move to the order α 2 s contribution of to eq. (2) and construct the appropriate antenna subtraction terms following the lines of [3][4][5][6].
The real-virtual contribution to the differential cross section associated with the process (5) reads where N c denotes the number of colors, n f the number of massless quarks, andC(ǫ) = 8π 2 C(ǫ) = (4π) ǫ e −ǫγ E . The normalization factor N 0 comprises all non-QCD couplings as well as the spin averaging factor for the initial state and the flux factor. Furthermore, we have introduced a shorthand notation for the interference of the tree-level and renormalized 1-loop amplitudes of (5) (with all couplings and color matrices factored out): Here the superscript X denotes the leading color (lc), subleading color (sc), and the massless (f ) and massive (F ) fermion contributions. For the sake of brevity, we have dropped the dependence on the initial state momenta. In (6) and in the formulae below summation over all helicity states is implicit.
The measurement function J (n) 2 in (6) ensures that only configurations are taken into account where n outgoing partons form two heavy quark jets. We emphasize that the following discussion applies to any infrared safe observable associated with the reactions (1). The renormalized one-loop amplitudes M X (3,1) contain explicit IR poles due to the exchange of virtual massless partons. On the other hand, integrating (6) over regions of the 3-particle phase space which correspond to soft gluon emission leads to additional infrared singularities. Both types of singularities have to be subtracted with appropriate counterterms in order to perform the integration over the three-parton phase space numerically. It is a well-known fact from NLO QCD calculations that the explicit IR poles in (6) can be removed by adding the subtraction terms dσ S,a,QQgg NNLO and dσ S,a,QQqq NNLO , integrated over the oneunresolved parton antenna phase space, which govern the singularities due to single unresolved emission in the squared tree-level matrix elements of S → QQgg and S → QQqq: The terms on the r.h.s. are given in [15,16].

One-loop single-unresolved subtraction term
The singular behavior of (6) in the limit where the external gluon becomes soft requires another subtraction term which has the following structure [4,5]: where and M (2,i) , i = 0, 1, denote the tree-level and 1-loop matrix elements of S → QQ (with all couplings and color factors stripped off). These reduced amplitudes are evaluated at redefined on-shell momenta p 13 , p 23 (i.e. p 13 2 = p 23 2 = m 2 Q ) which are defined by Lorentz-invariant mappings p 1 , p 2 , p 3 → p 13 , p 23 . These mappings are given in [11,12].
The massive tree-level antenna function A 0 3 and the massive one-loop antennae A 1 3 ,Ã 1 3 ,Â 1 3,f , andÂ 1 3,F in (9) depend only on the original momenta p 1 , p 2 , p 3 . They can be derived once and for all from the tree-level and 1-loop matrix element of γ * → QQg. The corresponding 1-loop amplitudes are shown in Fig. 1. We define and likewise forÃ 1 3 with δM γ * ,lc (3,1) replaced by δM γ * ,sc (3,1) . The superscript γ * indicates that the quantities defined in (7) and (10) have to be determined specifically for S = γ * . We note that M γ * (2,0) and M γ * (2,1) only depend on q 2 (and m Q ). The antenna functionsÂ 1 3 2 arise from ultraviolet counterterm contributions to the renormalized one-loop term (6) and are therefore proportional to the tree-level antenna function A 0 3 , which was derived in [10]. The computation of δM γ * ,X (3,1) , δM γ * (2,1) , and M γ * (2,0) 2 is standard. The resulting expressions for 1 3,F are quite long and will be presented elsewhere [55]. Because (9) is not induced by an integrated double-real subtraction term, it has to be added back to the double virtual cross section in its integrated form. The integrated subtraction term can be cast into the form where we have introduced the integrated antenna functions and likewise forÃ 1 3 ,Â 1 3,f ,Â 1 3,F , and A 0 3 . In (13) we have made use of the fact that the antenna phase-space measure dΦ X QgQ is defined as the ratio of the usual three-particle phase-space dΦ 3 (associated with a massive quark-antiquark pair and a gluon) and the phase-space volume P QQ ≡ P QQ (q 2 , m 2 Q ) of a heavy QQ pair with total four-momentum q. Consequently, the integrated massive one-loop antenna functions In (12) the squared tree-level matrix element |M (2,0) | 2 and the interference term δM (2,1) depend on the heavy (anti)quark momenta p 1 and p 2 .
Because the antenna functionsÂ 1 3,f andÂ 1 3,F are proportional to the massive tree-level antenna A 0 3 , the integrated antenna functionsÂ 1 3,f andÂ 1 3,F are proportional to the integrated massive antenna function A 0 3 , which was already determined in the context of NLO antenna subtraction [10,11]. The calculation of A 1 3 andÃ 1 3 will be described in the next section.

Compensation terms for oversubtracted poles
The 1-loop antenna functions contain explicit IR poles which arise from the loop integrations in δM γ * ,X (3,1) and δM γ * (2,1) . The same holds for the reduced 1-loop matrix element M (2,1) in (9). These poles do not coincide with respective poles of the real-virtual contribution (6) to the differential cross section, because those have already been subtracted with the help of dσ T,a,QQg NNLO . In order to remove these spurious singularities one has to introduce an additional subtraction term which we denote by dσ T,c,QQg NNLO . This additional subtraction term is given by (cf. eq. (4) and [4,5]) To conclude, combining Eqs. (6), (8), (9), and (14) yields an expression which is free of (explicit and implicit) singularities in the entire three-parton phase space in d = 4 dimensions: As discussed above, the introduction of dσ T,a,QQg NNLO and dσ T,c,QQg NNLO is counterbalanced by integrated double-real subtraction terms. Hence, only dσ T,b,QQg NNLO has to be added back in its integrated form (12) to the double virtual contribution.

Double-virtual corrections
In order to cancel the explicit IR poles of the double virtual contribution dσ V V,QQ NNLO to (2), we have to add the integrated subtraction terms resulting from dσ S,b,2,QQgg NNLO and dσ S,b,2,QQqq NNLO (cf. eq. (4)) and (12): The integrations over the respective phase spaces of the unresolved partons have to be done in d = 4 − 2ǫ. Only after summing up all the terms in (16) and analytic cancellation of the IR poles, one can take the limit ǫ → 0 and perform the remaining integration over the two-parton phase space. Using (12) the subtraction term for the double-virtual correction has the following structure: We recall that the integrated antenna functions A 0 4 ,Ã 0 4 , B 0 4 , A 1 3 ,Ã 1 3 ,Â 1 3,f ,Â 1 3,F , and A 0 3 depend on µ 2 /q 2 and y = (1 − β)/(1 + β). The massive tree-level quark-antiquark antenna A 0 3 was computed in [10,11], and the integrated massive four-parton antenna functions A 0 4 ,Ã 0 4 , and B 0 4 which result from the second and third integrals in the first line of (17), were determined in [16] and [15], respectively. Finally, explicit results for the integrated massive 1-loop antenna A 1 3 , A 1 3 ,Â 1 3,f , andÂ 1 3,F are presented in the next section. In conclusion, the expressions (3), (15), and (16) constitute, within the antenna subtraction method, the IR-finite four-parton, three-parton, and two-parton contributions to (2). Let us stress again that the antenna subtraction terms that we have introduced in this section are applicable to any reaction of the type (1) for an arbitrary colorless initial state S.

Computation of
The integrated antenna functions A 1 3 andÃ 1 3 are computed as follows. We represent the d-dimensional three-particle phase-space measure in terms of cut-propagators [35]. We perform an integration-by-parts reduction [37] using the computer implementation FIRE [38] of the Laporta algorithm [39] in order to express these two integrated antenna functions by 22 master integrals which correspond to three-particle cuts through three-loop scalar self-energy type Feynman integrals, involving massive and massless scalar propagators, associated with 11 different topologies, see Fig. 2. The terms in the curly brackets listed below each topology (s ij ≡ 2p i · p j ) denote factors which multiply the integrand of the corresponding integral. For instance, the two integrals which are represented by the diagram Fig. 2k are, with A ∈ {1, s 13 }: where D(k, m, P ) = (k − P ) 2 − m 2 + iη and dΦ d 3 denotes the 3-particle phase-space measure in d dimensions.
The master integrals are computed by the method of differential equations [43][44][45][46], i.e., we derive inhomogeneous first order differential equations in the variables q 2 and y = (1 − β)/(1 + β) for each master integral I A can be expressed in terms of cyclotomic harmonic polylogarithms [49][50][51][52][53]. The solutions of the first-order differential equations involve integration constants which we fix by determining the boundary condition that each master integral I A (α) must satisfy at the QQ production threshold β = 0. For this purpose we expand the d-dimensional 1-loop integral and the phase-space measure which make up each I A (α) im powers of ǫ and β. Phase-space integration then yields the coefficients of the Laurent series in ǫ of the I A (α) at y = 1. The integration constants which involve HPL at y = 1 can be straightforwardly expressed in terms of transcendental numbers [41,56]. Those integration constants that involve cyclotomic HPL at y = 1 can be computed either by numerical evaluation of the respective integral representation of the cyclotomic HPL or by (computer) algebraic reduction to transcendental numbers [49][50][51][52][53][54]. In this way, the 22 master integrals I A (α) and, with these integrals, the integrated antenna functions are calculated analytically. The complete expressions of the individual I A (α) will be given elsewhere [55]. Below we list the IR-divergent pieces of the integrated one-loop antenna functions defined in eq. (13). For the integrated leading and subleading color one-loop antennae, A 1 3 andÃ 1 3 , we find: − 25 12(1 + y) + 2y + 1 2 (y 2 + 4y + 1) H(0; y) − 4H(1; y) Figure 2: The 11 topologies which correspond to the 22 master integrals which determine the antenna functions A 1 3 andÃ 1 3 . Bold (thin) lines refer to massive (massless) scalar propagators. The dashed line represents the three-particle cut. The external double line represents the external off-shell momentum q. The terms in the curly brackets below each topology denote additional factors in the corresponding integrand of the combined phase-space and loop integral.

Computation of R Q
As a first application and, especially, as a check of the IR-singularity structure as well as of the finite parts of our results for the integrated antenna functions, we compute the order α 2 s contribution to the inclusive production of a massive quark-antiquark pair by e + e − annihilation via a virtual photon for arbitrary squared center-of-mass energy above the pair production threshold, s > 4m 2 Q . The ratio R Q is defined by To order α s , this ratio has been known for a long time [63,64]. The second order term R (2) may be decomposed into gauge-invariant pieces associated with the different color structures as follows: Here, e Q is the electric charge in units of the positron charge, and R LC , R SC , R f , and R (2) F denote the leading color, subleading color, and massless fermion contributions, and the heavyquark contribution from the QQ, QQg, and QQQQ (above s > 16m 2 Q ) final states, respectively. A complete calculation of the contribution R (2) f was made first in [65]. Within the antenna framework it was calculated in [15] in terms of harmonic polylogarithms and agreement was found between the analytical results of [15] and [65]. Therefore, we do not further consider this term here.
In [66] the leading and subleading color contributions to R (2) were computed from the imaginary part of the order α 2 s photon vacuum polarization function Π (2) (q 2 ), using analytically known expansions of Π (2) at q 2 = 0, near the QQ threshold and for −q 2 → ∞, and Padé approximations in the range above threshold and below the high-energy region. With the same techniques, the contribution of order α 3 s to (23) was calculated in [68]. Very recently, the term R (2) was computed also in [34]. Here we compute R (2) LC and R (2) SC by applying the framework devised in Sec. 2 to the case S = γ * . Up to a trivial normalization R (2) can be obtained by adding up the corresponding subtracted cross sections (3), (15), and (16) associated with the four-parton, three-parton, and two-parton final states, respectively. More precisely, these contributions have to be evaluated with the specific choice J (n) 2 = 1, which corresponds to inclusive QQ production. In this particular case, the subtracted three-and four-parton cross sections (3) and (15) yield a vanishing contribution because, by construction, the original matrix elements squared and the subtraction terms coincide. Therefore, all information on R (2) is encoded in the subtracted two-parton contribution (16). The matrix element which determines the second-oder QCD contribution dσ V V,QQ NNLO to (16) was computed in [57] (and confirmed in [58]) and is used here. The sum of the integrated subtraction terms in (16) is given in (17). Adding up the ǫ poles of dσ V V,QQ NNLO [57], of A 0 4 ,Ã 0 4 [16], of B 0 4 [15], of A 0 3 [10], of A 1 3 ,Ã 1 3 ,Â 1 3,f ,Â 1 3,F (cf. eqs. (19), (20), (21), and (22), respectively), and of δM γ * (2,1) we find that all IR poles cancel analytically. Therefore, our result for the NNLO QCD correction (24) to R Q is indeed finite. In addition to (16), the contribution from the QQQQ final state has to be taken into account, which has not been discussed in Sec. 2. The four-particle phase-space integral of the squared tree-level matrix element of γ * → QQQQ is completely finite and can be evaluated numerically.
(The respective phase-space integrals cannot be expressed solely by polylogarithms.) Denoting this contribution to (23) by R (2) 4Q , its color structure is where R E is an interference term which results from the fact that there are two identical (anti)quarks in the final state. This term makes a contribution to the subleading color term R (2) SC in (24). It becomes logarithmically divergent in the high-energy limit m 2 Q /s → 0, (cf. [62]): where ζ(x) denotes the Riemann zeta function. In [62] the constant c was extracted from a numerical computation of R E in the region of small m 2 Q /s and found to be c = −8.7190±0.0013. In the limit m 2 Q /s → 0, the contribution from (16) to R SC exhibits the same logarithm as in (26) but with opposite sign. As a result, R SC remains finite for m Q → 0 as it should be. For the constant c in (26) we obtain which is in good agreement with the numerical value of Ref. [62] cited above.
In order to compare our results for the leading and subleading color contributions with the results and the asymptotic expansions given in [66,[69][70][71][72] and with the threshold expansion of [59][60][61], we switch to the color decomposition used in these papers and define SC , R SC .
The solid lines in Fig. 3 show our results for R NA and R A in the range 0 < β ≤ 1. These results are exact -the usual harmonic polylogarithms which enter the expressions for R (2) NA and R (2) A are evaluated numerically with the Mathematica package HPL [56], whereas the remaining cyclotomic polylogarithms are evaluated by utilizing their defining integral representations [52]. For β → 0, the behavior of R (2) A is dominated by the Coulomb singularity ∼ 1/β, while R (2) NA diverges only logarithmically for β → 0. In the high-energy limit x = m 2 Q /s → 0 the terms R NA/A are known as asymptotic expansions in x. In [72] these expansions are given up to and including terms of order x 6 . These asymptotic expressions are shown by the dotted curves in Fig. 3. We have expanded our exact results for R (2) NA/A around x = 0 with the help of the Mathematica packages HPL [56] and HarmonicSums [52,53] and compared with the expansions given in [72] and find agreement. In case of the leading color contribution R (2) LC this comparison is performed in fully analytic fashion. Because we have not computed the higher-order terms in the expansion (26) in analytical fashion, we checked the small-x behavior of R (2) SC only numerically. In the literature, R NA and R (2) A are known analytically near the QQ production threshold as an expansion in β up to and including terms of order β [59] (cf. also [60,61]). These expansions are shown as dashed lines in Fig. 3. To this order in β both the non-abelian and abelian term R (2) NA/A is determined by the second order QCD matrix element of γ * → QQ. As stated above we use this matrix element from [57]. The threshold expansion of this matrix element agrees with [59]. We checked by small-β expansion of the polylogarithms which enter our expressions for (16) that the three-and four-parton final states contribute to R (2) NA/A only with higher powers in β. Explicit expansions of R (2) NA/A to higher orders in β will be presented elsewhere [55]. Ref. [66] computed R (2) NA/A in the whole physical region 0 < β ≤ 1 by using the asymptotic and the small q 2 expansions of the photon vacuum polarization function Π (2) (q 2 ) as input for Padé approximation of Π (2) (q 2 ) (whose imaginary part yields R (2) ) away from the threshold and the asymptotic region. Figs. 5 of [66] show the result of this calculation, in particular R  A plotted against β (solid line). The renormalization scale is chosen to be µ = m Q . For comparison, the expansions at the kinematic threshold (dashed curves) [59,61] and in the asymptotic region (dotted curves) [72] are included as well.
Thus our antenna functions pass this non-trivial test.

Summary and Outlook
We have computed the integrated real-virtual antenna functions which are required, within the antenna subtraction framework, for the calculation of the differential cross section for massive quark-pair production by an uncolored initial state at NNLO QCD. This completes, together with our previous results [15,16] on the subtraction terms for the matrix elements of the final states QQqq and QQgg, the set of antenna subtraction terms which are required for this type of processes. As a check of our result, we have computed the inclusive cross section, respectively the ratio R Q (s) for e + e − → γ * → QQX at order α 2 s by adding up the antenna-subtracted cross sections associated with the two-parton, three-parton, and four-parton final states. Our result is infrared-finite as it should, and our second order correction to R Q (s), which we obtained without any approximation, agrees with the known threshold [59] and asymptotic [72] expansions in the respective kinematical limits.
Future applications of these subtraction terms include the computation of differential distributions for S → QQX at NNLO QCD. Our (integrated) antenna functions are, in addition, also important ingredients for the antenna-subtraction approach to NNLO QCD analyses of hadro-production of massive quark pairs.