NNLO QCD corrections to $pp \to \gamma^* \gamma^*$ in the large $N_F$ limit

We compute the NNLO QCD corrections for the hadroproduction of a pair of off-shell photons in the limit of a large number of quark flavors. We perform a reduction of the two-loop amplitude to master integrals and calculate the latter analytically as a Laurent series in the dimensional regulator using modern integration methods. Real radiation corrections are evaluated numerically with a direct subtraction of infrared limits which we cast in a simple factorized form. The results presented here constitute a gauge invariant part of the full NNLO corrections but are not necessarily dominant. We view this calculation as a step towards a complete computation. Our partial corrections to the total cross-section are about $1\%-3\%$ and vary with the virtuality of the two off-shell photons.


Introduction
The Tevatron and the LHC have performed studies on a wide spectrum of processes which probe the electroweak sector of the Standard Model. In particular, the production processes of a pair of electroweak gauge bosons [1][2][3][4][5][6][7][8][9][10][11][12] are of great interest as they allow to test the electroweak theory, constrain physics beyond the Standard Model and are background to signals of the Higgs boson decaying into H → W W, H → ZZ. While the bulk of the cross-sections is due to on-shell production of the W or Z bosons, off-shell production is interesting especially for the background estimation in Higgs searches.
In this publication we make a first step towards the computation of NNLO corrections for diboson production in the case of two off-shell electroweak gauge bosons. We restrict ourselves to computing the NNLO cross-section for an idealized process pp → γ * γ * in the limit of a large number of massless quark flavors N F . While the large-N F limit is not necessarily dominant it provides the opportunity of obtaining a gauge invariant part of the cross section and serves as an excellent means to treat and develop analytic and numeric methods. We generate and reduce the required amplitudes to master integrals using established methods [47][48][49][50]. We evaluate the latter by directly performing the integrations over the Feynman parameter following methods similar to the ones introduced in refs. [51][52][53][54][55][56][57][58][59]. As a by-product, we construct a set of basis functions up to transcendental weight four with the correct branch cut structures which are sufficient to write down the answer for the class of integrals studied in this paper. Moreover, the master integrals presented here have been computed independently and agree numerically with the results of refs. [60][61][62]. For the calculation of real radiation corrections we apply a subtraction scheme based on a hierarchical parameterization of the phase-space and the universal collinear and infrared limits of the squared matrix-elements. All singularities cancel after adding the partonic cross-sections together and performing UV renormalization.
This article is organized as follows. In section 2 we present our notation and setup of the calculation. In section 3 we present the calculation of the two-loop amplitude in the large N F limit and we outline the computations of the relevant master integrals in section 4. The computation of corrections due to real radiation and our subtraction scheme are presented in sections 5 and 6. We demonstrate the numerical impact of the contributions that we have computed here in section 7. We conclude in section 8.

Setup and notation
In this article, we compute the fully differential cross-section at the LHC for the process of producing two idealized off-shell photons, where P denotes a proton and X is a shorthand notation for the associated QCD finalstate radiation. In parentheses we indicate the momenta of the external particles.
We compute cross sections which are fully differential in the momenta p 3 and p 4 of the photons, as well as in the momenta of the associated QCD jet radiation. The hadronic cross section for a generic observable J is given by with σ ij→γ * γ * X [J ] denoting the differential cross section for the process where i and j run over the parton flavors g, u,ū, d,d, . . . relevant to this process, p 1 = x 1 P 1 and p 2 = x 2 P 2 are the momenta of the initial-state partons and f b i (x) the bare parton distribution functions (PDFs). The function J depends on the finalstate momenta and restricts the phase-space to the desired infrared-safe observable.
The partonic cross sections are computed as a perturbative expansion in the bare strong coupling constant α b s , where k and l run over the final-state parton flavors. The partonic cross sections with definite final state γ * γ * , γ * γ * q, γ * γ * q q , etc, are given by: where s = 2p 1 · p 2 is the partonic center-of-mass energy squared and |M ij→γ * γ * ... | 2 (m) is the m−loop contribution to the ij → γ * γ * . . . amplitude squared, summed over spin and colour and averaged over initial state quantum numbers. We compute the matrix elements using conventional dimensional regularization in d = 4 − 2 spacetime dimensions. We assume that the photons do not decay and use the polarization sum: where p denotes the photon-momentum. We consider N F = 5 light quark flavours and we ignore the effects of the top-quark both in the loops and the evolution of the strong coupling.
In the present article, we compute the complete O(α s ) corrections, while at O(α 2 s ) we retain only the gauge-invariant terms which contribute in the N F → ∞ limit. Some tree and two-loop diagrams that contribute to the NNLO large-N F Figure 1: Sample tree and two-loop diagrams contributing to the NNLO corrections for qq → γ * γ * in the large-N F limit. correction are shown in figure 1. The two-loop diagrams contributing to the large N F limit are in one-to-one correspondence with the one-loop diagrams appearing at NLO, by replacing the gluon propagator by its one-loop self energy graph. At NNLO, the partonic processes which contribute to the correction are qq → γ * γ * q q and qq → γ * γ * qq. In the latter process, we retain only the interference terms with two spin lines.
The Lorentz invariant phase space is given by and we define the following Mandelstam variables and their ratios: Ultraviolet renormalization is performed in the MS scheme. The bare strong coupling constant α b s is given in terms of the renormalized coupling α s (µ) as where β 0 and β 1 are the first and second coefficients of the QCD beta function 2Nc , N C = 3, T R = 1 2 ; and S = e (log 4π−γ E ) . Since the Born cross section is independent of α b s , only the α 2 s (µ) term of eq. (2.8) is required for renormalization. We absorb the initial-state collinear singularities into the parton densities in the MS-factorization scheme. The bare PDFs f b i (x) are written in terms of the renormalized PDFs f j (x, µ) as where implicit summation over j is understood, and the convolution integral is defined as (2.10) The kernels ∆ (1,2) ij can be written in terms of the Altarelli-Parisi splitting kernels as The splitting kernels relevant for this computation are 14) (2.15) The D n (1 − z) plus-distributions are defined as qq we need only the terms proportional to N F . We remark, however, that in the numerical evaluation of the PDFs and the strong coupling from their values at their initial scales we use the complete β−function and Altarelli-Parisi kernels and not just their N F parts.
In the rest of this article we will set the renormalization and factorization scales to be equal, µ f = µ r ≡ µ. The generic dependence on both scales can be easily restored by first setting µ = µ f and writing:

Virtual corrections
Ingredients of the NLO and NNLO corrections are the one-loop and two-loop amplitudes for the partonic process qq → γ * γ * . We generate the required Feynman diagrams using QGRAF [47] and then compute the interference of the one-loop amplitude and the tree-amplitude as well as the interference of the two-loop amplitude and the tree amplitude, summing over external-state colours and polarizations. We perform the Dirac and colour algebra with programs implemented in the FORM [48] programming language.
From the interference of the tree and two-loop amplitudes, we keep only the terms which contribute to the large N F limit. These are expressed in terms of two-loop integrals of the form: where we have used the shorthand notation q 1···n ≡ q 1 + · · · + q n and the external momenta q i take the values: (q 1 , q 2 , q 3 ) ∈ {(p 1 , p 2 , p 3 ), (p 1 , p 2 , p 4 )}. The powers n i take integer values in the range n i ∈ [−4, 2]. These integrals are not independent and they can be reduced to a basis of six master integrals. We use the program AIR [50] based on the Laporta algorithm [49], and obtain the following two-loop master integrals: T 2 (0, 1, 0, 0, 1, 0, 1, 1, 1, p 1 , p 2 , p 4 ) ≡ (3.6) T 2 (0, 1, 0, 0, 1, 0, 1, 1, 2, p 1 , p 2 , p 4 ) ≡ (3.7) The same integrals with p 3 and p 4 exchanged also appear in the two-loop amplitude. Similarly, the interference of the tree and one-loop amplitudes can be expressed in terms of integrals of the form: where the integer powers n i range in [−4, 1]. The one-loop integrals are reduced to the following master integrals : The master integrals T 1 (1, 0, 1, 1, p 1 , p 2 , p 3 ), T 1 (1, 1, 1, 1, p 1 , p 2 , p 3 ) also appear in the one-loop amplitude.
In the following section, we present a computation of the required master integrals, as well as of some master integrals which are needed for the full calculation beyond the large N F limit. The complete set of master integrals contributing to diboson production at two-loop order was recently computed in ref. [62]. We have performed an independent computation and confirm these results.

Master Integrals
In this section we present the analytic results for all master integrals that enter the N F -part of the amplitude for qq → γ * γ * up to two-loop order.

Analytic results in the Euclidean region
We start by giving the analytic results for the master integrals in the Euclidean region where all consecutive Mandelstam invariants are negative. Note that in this region the variables u, v and w defined in section 2 are all positive. The results with the two virtualities p 2 3 and p 2 4 exchanged can easily be obtained from the replacement Before presenting our results, we first discuss some general properties of the integrals. In dimensional regularization with d = 4 − 2 , every master integral is computed as a Laurent series in , whose coefficients are expressed in terms of polylogarithmic functions. The simplest possible representatives of this class of functions are the ordinary logarithm and classical polylogarithms, defined by with Li 1 (x) = − log(1 − x). However, more general functions can also appear. These are the multiple polylogarithms [63,64], defined by G( 0 n ; r) ≡ 1 n! log n r and G(a 1 , . . . , a n ; r) = r 0 dt t − a 1 G(a 2 , . . . , a n ; t) , with G(r) = 1 and the arguments a i , r ∈ C. The number of elements of the vector a = (a 1 . . . , a n ) is called the weight of the multiple polylogarithm. Note that up to weight three, all multiple polylogarithms can be expressed in terms of classical polylogarithms and ordinary logarithms. In particular, the two-loop amplitude for qq → γ * γ * in the large N F limit only involves polylogarithmic functions up to weight three (up to O( 0 )), and hence we can always express our two-loop amplitudes in terms of classical polylogarithms only. This greatly facilitates the numerical evaluation. This point will be discussed in more detail in section 4.2 when discussing the analytic continuation from the Euclidean region to the Minkowski region. The arguments of the polylogarithms are in general algebraic functions of the Mandelstam invariants, and in particular they involve the square root λ (1, u, v), where λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc denotes the Källén function. A convenient parameterization which rationalises this square root is given by or equivalently This choice of parameterization is inspired by ref. [51], where it was argued that the variables (r,r) define a natural set of variables for parameterizing the kinematics of a massless three-point function with all external legs off shell. These integrals naturally appear as master integrals in our case. Furthermore, it was shown in ref. [51] (see also refs. [65,66]) that in the region where λ(1, u, v) < 0, such that r andr are complex conjugate to each other, massless three-point functions are described by single-valued functions in the complex r plane. Indeed, it is well-known that massless loop integrals can only have branch cuts starting at points where one of the Mandelstam variables vanishes. The single-valuedness condition is equivalent to the condition that these functions have the correct physical branch cuts. The advantage of this approach is that for every weight, there is only a very limited set of single-valued functions. In ref. [51] a method was presented to construct these functions explicitly up to weight four in the case of massless three-point functions (see also refs. [67] for similar ideas). In particular, up to weight three only three functions can appear besides the ordinary logarithms, log u = log(rr) and log v = log(1 − r)(1 −r). Following ref. [51], we denote these functions by P 2 (r), P 3 (r),P 3 (1 − r) and Q 3 (r). The functions P n (r) are closely related to the so-called Bloch-Wigner function, where B k denotes the k-th Bernoulli number and R n denotes the real part for odd n and the imaginary part otherwise. Note that the function defined by eq. (4.7) is a combination of classical polylogarithms without branch cuts for r ∈ C, and it is therefor natural to call the functions (4.7) the single-valued versions of the classical polylogarithms. The function Q 3 (r) is defined by Up to weight three and two loops, all massless three-point functions can be written as linear combinations of (products of) these functions [51] (with coefficients that are Q-linear combinations of ζ values). The previous considerations, however, only apply to massless three-point functions. It is nevertheless straightforward to generalise these ideas to four-point functions with two adjacent off-shell legs. In appendix A we present a way to construct a set of basis functions up to weight four with the correct physical branch cuts contributing to the large N F limit of the the qq → γ * γ * amplitude at two loops. In the following we only concentrate on the set of basis functions up to weight three, which is relevant in the present case. Besides the functions defined in eq. (4.6 -4.8), we find six possible classical polylogarithms, and two new functions R ± 3 (r, w) ≡ R ± 3 (r,r, w), where the superscript '±' refers to the parity of the functions under the exchange r ↔r, R + 3 (r, w) = G(0, v,r(−1 + r); w) + G(0, v, (−1 +r)r; w) First, we emphasise that each of these functions has the correct branch-cut structure corresponding to a massless four-point function with two adjacent off-shell legs, i.e., they have branch cuts at most starting at points where one of the external Mandelstam invariants vanishes. Second, this set of functions is linearly independent, i.e., it is not possible to express any of these functions as a linear combination of (products of) all the others. It is therefore justified to call these functions a set of basis functions. As a consequence, all master integrals contributing to the large N F part of the qq → γ * γ * two-loop amplitude can be expressed as a unique linear combination of (products of) basis functions. The construction of these functions, as well as the proof that they form a basis, is given in appendix A.
In the rest of this section we collect our results for the master integrals contributing to the large N F part of the qq → γ * γ * two-loop amplitude. Details about the computation can be found in appendix B. All the expressions are valid in the Euclidean region, and the results are given in terms of the basis functions we have just defined. We explicitly show the results up to weight three. Analytic results up to weight four are provided as ancillary files with the arXiv submission.
We checked that our results satisfy the differential equations for the master integrals. Moreover the results were checked numerically with FIESTA [68], which is based on the method of sector decomposition [69] (the multiple polylogarithms were evaluated using GiNaC [70,71]. In addition, we have compared our results with existing results in the literature whenever available [51,[60][61][62][72][73][74][75].
One-loop integrals. We start by summarising the one-loop integrals. The relevant one-loop two, three and four-point functions are given by where γ E = −Γ (1) denotes the Euler-Mascheroni constant and we introduced the usual normalization factor Note that all results are entirely expressed in terms of the basis functions defined at the beginning of this section, as expected. The one-loop box has been previously computed up to the finite part in the -expansion in ref. [13].
Two-loop integrals. In this subsection we give the analytic expression for the two-loop integrals. Besides the loop integrals necessary for the amplitudes presented in this work, we also display all the boxes with bubble insertions with two adjacent off-shell legs. The master integrals are presented up to the order in that corresponds to coefficients of weight up to three. The full results including coefficients of weight four can be found in the file attached as ancillary files to the arXiv submission of the paper.
After we completed the computation of these integrals, a complete basis of planar master integrals for the production of off-shell vector bosons was presented in ref. [62].

Analytic continuation into the physical region
The results of the previous section are only valid in the Euclidean region, where s, t, p 2 3 , p 2 4 < 0, such that u, v, w > 0. In this section we perform the analytic continuation into the region defined by The analytic continuation can be performed via the usual replacements This implies that the ratios u, v and w are analytically continued according to the prescription In ref. [62] it was shown how to analytically continue the multiple polylogarithmic functions to the physical region using this prescription. In the following, we present an alternative way of performing the analytic continuation, which will allow us in the end to express the large N F part of the amplitude entirely in terms of classical polylogarithms of weight three at most and with arguments in the interval [0, 1] everywhere in the physical phase space. Consequently, the classical polylogarithms are real and admit a convergent power series representation. The advantage of this representation is very fast and stable numerical evaluation. It turns out, however, that in order to obtain such a representation, we need to split the phase space into three different regions, such that a representation of the desired type exists in each region. We first discuss these different regions, and present our procedure to perform the analytic continuation in each region at the end of this section.
Definition of the regions. In the following we describe how to identify the parts of the physical region in which the results can be expressed in terms of classical polylogarithms with arguments inside the range [0, 1], where all the functions are convergent. The results of the previous section were valid in the Euclidean region where λ(1, u, v) < 0, and thus r andr complex conjugate to each other. Without loss of generality we may assume that in that region we have Im r > 0 and Imr < 0. (4.29) It is easy to check that the physical phase space, however, corresponds to λ(1, u, v) > 0, i.e. r andr real. In ref. [51] it was shown that the correct prescription for the analytic continuation from λ(1, u, v) < 0 to λ(1, u, v) > 0 while keeping u and v real is It is sufficient to work out the analytic continuation for the basis functions. Moreover, since in the physical region √ s > m 3 + m 4 , we must have 0 < u, v < 1, which implies 0 <r < r < 1 [51]. In the following we show that we must furthermore havē r <w + rr < r in the physical region.
In order to show this inequality, we parameterize the external momenta as where e 3 = (0, 0, 1) T and (4.32) We thus obtain for p 3 and so for some θ ∈ [0, π], and where we used rotational invariance to remove the dependence on the azimuthal angle. At this point we can already conclude that λ(1, u, v) > 0, i.e. r andr are indeed real and moreover we see from eq. (4.32) that 0 <r < r < 1.
Using this parameterization, we find and sor <w + rr < r.
In the end we assert that the only non-trivial part in switching to the physical region is the analytic continuation of the basis functions depending on w. Besides the analytic continuation in w, some of the functions appearing in the basis functions are not defined for arbitrary values of r,r and w. Consider for example which develops an imaginary part if −w >r(1 − r). We find that we have to split some of the basis functions for physical values into three different regions where τ is defined throughw In deriving the analytic continuation of functions depending on w we have to keep in mind these different regions. Note that the physical phase space corresponds to 0 ≤ τ ≤ 1.
Analytic continuation of the functions. In this section we demonstrate how to perform the analytic continuation (4.27). Our main goal is to obtain a representation of the amplitude in the physical region in terms of classical polylogarithms up to weight three with arguments lying in the range [0, 1], such that the polylogarithms admit a convergent power series representation. Technically speaking, we are looking for a functional equation which allows us to express the amplitude in terms of functions that are real in the physical region, and where all the imaginary parts are explicit. Functional equations among multiple polylogarithms are most conveniently described in terms the Hopf algebra of multiple polylogarithms (see appendix A). In a nutshell, multiple polylogarithms admit a coproduct structure which allows to decompose a polylogarithm of weight n into a sum of pairs of polylogarithms of weight (k, n − k). It is then possible to find functional equations among multiple polylogarithms of weight n recursively by first decomposing them into functions of lower weight, for which all relations are assumed to be known. Let us illustrate this on some simple examples. First, we know that there are only three basis functions of weight one, and their analytic continuation follows immediately from eq. (4.27), Next, consider one of the basis functions of weight two, and let us consider Li 2 1 − u w as a representative example. In the physical region where w = −w < 0, the argument of the dilogarithm becomes greater than 1, and so the dilogarithm develops an imaginary part. Acting with the coproduct, we obtain, where in the second line we used eq. (4.27). Note that by construction every basis function of weight n will only contain log u, log v or log w in the first factor of its (1, n − 1) component of the coproduct (see appendix A for details). The imaginary part of this dilogarithm can immediately be read off from the second term, At this point we need to find a real function whose coproduct matches the real part of eq. (4.40). It is easy to check that Hence, we can conclude that (4.43) We can thus conclude that the arguments of ∆ 1,1 are equal, up to (constant) terms that vanish when acting with ∆ 1,1 . In order to determine this constant, we expand both side close to the branch point at w = 0, where in the first line we have used the fact that log w = logw + iπ. Equating the two expressions, we see that Analogously we obtain the analytic continuation of all the basis function depending on w.

Single-real contributions
We now turn our attention to the phase-space integrations over tree-level matrix elements for partonic processes, where the photon pair is produced in association with an additional parton in the final state: As before, we denote in brackets the momenta of the partons. These processes contribute at NLO to the hadronic process, and via renormalization and massfactorization also at NNLO.

Quark-antiquark channels
We first consider the channels with a qq-pair in the initial state. The corresponding tree-level cross section is given by where |M qq→γ * γ * g | 2 (0) is the qq → γ * γ * g tree matrix element squared, summed over spin and colour and averaged over initial-state quantum numbers. The case where the quark and the anti-quark are exchanged is identical. The phase-space measure can be decomposed into a phase space producing a gluon and an intermediate offshell particle in the final state with momentum Q of virtuality Q 2 = zs, and a phase space for the decay of the intermediate particle into two photons as where We parameterize the momentum of the gluon as such that where z, λ ∈ [0, 1] and e T is a unit vector transverse to p 1 and p 2 in d = 4 − 2 dimensions. In this section and the following ones we use the shorthand notation where we use dΩ d−2 to denote the differential solid angle generating e T . The matrix element squared in the integrand of eq. (5.1) is singular in the collinear limits p g p 1 and p g p 2 (corresponding to λ → 0 and λ → 1 respectively) and in the soft limit p g → 0 (corresponding to z → 1). The singular behaviour of matrix elements squared is universal [76], in the sense that it is independent of the process under consideration. In particular, the formulae presented here are valid for any colourless final state in place of the two off-shell photons γ * γ * . Explicitly, we have where g 2 s = 4πα b s , and the splitting kernel is given by and In the above, the squared matrix elements for the Born process qq → γ * γ * are evaluated with the momentum in the collinear direction rescaled by a factor of z.
Since we consider a colourless final state, the soft limit does not involve any colour correlations and can be simply written as Note that the sum of the collinear limits, given by equations (5.8) and (5.9), reproduces eq. (5.13) exactly in the limit where z → 1. This means that although the matrix element squared is singular in the soft limit, no explicit subtraction of this singularity will be needed. We now recast the partonic cross-section as: where σ H qq has an integrand which is finite in all singular limits (λ,λ,z → 0) as we take to zero and it is therefore allowed to perform a Taylor expansion in , while σ C 1 qq and σ C 2 qq are divergent as → 0.
The contributions read and We extract the pole in of σ C 1 qq and σ C 2 qq by integrating over the variables λ and e T . Since J , B i , and Q do not depend on these variables anymore, this step is straightforward. The result is still singular in the z → 1 limit and we use an expansion in plus-distributions to extract this last singularity. We obtain and where the Born cross section σ (0) qq→γ * γ * [J ] is evaluated with rescaled momenta in the collinear direction, and the integrated splitting kernel is with P (0) qq (z) being the Altarelli-Parisi splitting kernel (2.13). The partonic cross section can then be subtracted using eq. (5.14). Recalling that p 1 = x 1 P 1 and p 2 = x 2 P 2 , and using eqs. (5.18) and (5.19), we obtain (5.20) where we used the trivial identity with y = xz. The first term of eq. (5.20) is finite, while the second and third terms contain all the poles in .

(Anti-)quark gluon channels
The remaining channels qg → γ * γ * q, gq → γ * γ * q,qg → γ * γ * q , and gq → γ * γ * q are treated similarly, and it is only necessary to consider the channel qg → γ * γ * q. We parameterize, as before, The matrix element squared is finite in in the limit where p q p 1 but is singular when p q p 2 , with the asymptotic behaviour Note that since we consider averaged matrix elements squared (with d − 2 polarizations for the gluons), we needed to compensate for averaging factors. The matrix element squared has a simple pole at z = 1, but since the phasespace measure (5.7) vanishes linearly in this limit, the cross section is free of soft singularities.
We subtract as before, writing where with where the Altarelli-Parisi splitting kernel is given by eq. (2.14). As before, the corresponding partonic cross-section can be written as where the second term contains all the poles in .

Double-real contributions
We now consider the double-real contributions to the partonic cross sections. As explained in section 2, only the channels qq → γ * γ * q q andqq → γ * γ * q q contribute to the large N F limit. We will first consider observables that do not involve differential information about the final-state quarks separately. In the second part of this section we then present a fully differential subtraction scheme.

Semi differential subtraction
Restricting ourselves to observables that do not resolve any of the differential properties of the final state quarks allow us to write J (p 1 , p 2 , p 3 , p 4 , p q , pq ) = J (p 1 , p 2 , p 3 , p 4 , p g * ), (6.1) where p g * is the momentum of the parent off-shell gluon, p g * = p q + pq . The phase space of the final-state quarks can then be integrated out, simplifying the extraction of limits. Hence we first consider where int. indicates that we restrict ourselves to the aforementioned observables. The phase space can be factorized as with s g * = p 2 g * , and since the function J does not depend on p q and pq , we can perform the integration over the decay phase space of the off-shell gluon explicitly. We obtain where |M qq→γ * γ * g * | 2 (0) is the averaged tree-level matrix element squared for the production of two off-shell photons and an off-shell gluon 1 , and A( ) is given by such that eq. (6.2) becomes To perform the subtraction for σ (0),int.
qq→γ * γ * q q we parameterize the momentum of the off-shell gluon as where z, λ, ρ ∈ [0, 1] and e T is again a unit vector transverse to p 1 and p 2 in d = 4−2 dimensions. We obtain the invariants Using this parameterization, the phase-space measure reads where dΩ d−2 denotes the differential solid angle parameterizing e T , and eq. (6.6) becomes The singular limits of the matrix element squared are once again universal but are asymmetric as a consequence of the asymmetry of the parameterization (6.7) under the exchange p 1 ↔ p 2 . We have to consider the following singular limits: • p g * p 1 : This corresponds to λ → 0, such that p g * →zp 1 , and the matrix element squared has the asymptotic behaviour with S qq;1 (z, ρ) = C F 2z + (1 − )z 2 ρ . (6.13) • p g * p 2 : This corresponds to λ → 1, such that p g * →zp 2 , and the matrix element squared has the asymptotic behaviour • s g * = 0: This is the final state collinear singularity, when the gluon becomes on-shell, but remains in the hard region. It corresponds to ρ → 1, and the matrix element squared has the smooth limit The corresponding singularity comes from the factor s −1− g * in (6.11).
Note that in the limit where ρ → 1, both splitting kernels (6.13) and (6.15) tend smoothly to the splitting kernel we obtained in the previous section, eq. (5.10). We proceed to the subtraction by writing (6.17) where σ HH qq has a Taylor expansion around = 0 and the other terms contain all the poles in . The different contributions are: and where we have p g = lim ρ→1 p g * such that the parameterization (6.7) tends to eq. (5.4) in the limit where ρ → 1. As for the real radiation, the soft limit does not need to be subtracted explicitly. This can be checked by expanding the integrand of eq. (6.18) around z = 1. We can extract the poles in of σ CC 1 qq and σ CC 2 qq by integrating over the variables λ, ρ and e T . The integration is slightly more complicated than for the real contributions, but the result can be written in terms of hypergeometric functions 2 F 1 (a, b, c;z), where a and b depend on . We have used the HypExp [77] package to expand them in and then performed a plus-distribution expansion overz to extract the double soft singularity. We note here that the residue of the soft pole at z = 1 is the same for both counterterms σ CC 1 qq and σ CC 2 qq , such that the asymmetry due to the parameterization is limited to the regular coefficients and does not affect the delta-and plus-distribution terms. As for the real corrections, we can write the counterterms as where we have The pole in of σ R;C qq is extracted by integrating over ρ only, since the variables z and λ still parameterize the on-shell gluon in ρ → 1 limit. Using we can expand as 27) where σ H qq is our NLO expression (5.15), and we defined After summing over the final-state quark flavours q , the first term of eq. (6.27) will cancel the β 0 term coming from the renormalization of α b s in the large N F limit, given by eq. (2.8), applied to the NLO contribution σ H qq , given by eq. (5.15). Also note that σH qq has the exact the same structure as σ H qq , except for the prefactor and the term on the second line.
At the partonic level, we finally obtain (6.29) where the first term is finite and the others contain all the poles in .

Fully differential subtraction
In this section, we show how to perform a similar subtraction in the case where the phase space of the two final state quarks cannot be integrated out. Hence we consider 2s dΦ 12→γ * γ * q q J (p 1 , p 2 , p 3 , p 4 , p q , pq )|M qq→γ * γ * q q | 2 (0) , (6.30) where the function J depends now on all external momenta. We extend our parameterization in a way similar to ref. [78]. In order to construct the momenta of the q q pair, we introduce another transverse unit vector e T , with the conditions e T · p 1 = e T · p 2 = e T · e T = 0 and e T 2 = −1. 2 The full phase space measure then reads where dΩ d−2 and dΩ d−3 denote the integrals over e T and e T respectively, and Q = p 1 + p 2 − p q − pq . The new variables y 1 , y 2 ∈ [0, 1] parameterize the phase space of the decay of the off-shell gluon. The expressions for the invariants s 1q = (p 1 − p q ) 2 , s 1q = (p 1 − pq ) 2 , s 2q = (p 2 − p q ) 2 and s 2q = (p 2 − pq ) 2 can be found in the aforementioned reference.
The momenta of the two final-state quarks can now be fully reconstructed and read p q =z λ y 1 p 1 + λ y 1 ρ +ρȳ pq =z λȳ 1 p 1 + λ ȳ 1 ρ +ρ such that the momentum of the parent gluon p g * = p q + pq is still given by our previous expression (6.7). Note that pq can be obtained from p q by replacing y 1 ↔ y 1 , y 2 ↔ȳ 2 and e T → −e T . There are no new singular limits to consider in this case. The singular limits where p q + pq p 1 and p q + pq p 2 have the same structure as in eqs. (6.12) and (6.14), and can be written as

33)
|M qq→γ * γ * q q | 2 (0) = 4g 4 s µ 4 S qq;2 (z, ρ, y 1 , y 2 ) with new, fully differential, splitting kernelsS qq;1 andS qq;2 . Note that although the singularities are now quadratic at the level of the matrix element squared, the cross section diverges only logarithmically because the phase-space measure (6.31) vanishes linearly in these limits. The singular limit where s g * = 0, such that p q pq , is now non-trivial and involves spin correlations, and reads In the above we have defined where pol. denotes the sum over the polarizations of the photons and the spins of the quarks and A qq→γ * γ * g µ is the tree amplitude for the process qq → γ * γ * g where the gluon is not contracted with the corresponding polarization vector.

Numerical results
We have implemented the various contributions to the differential cross section for the N F part of the process pp → γ * γ * + X up to the next-to-next-to leading order in the strong coupling expansion in two different programs. The virtual contributions are written in terms of master integrals which in turn are evaluated in terms of harmonic polylogarithms. In order to ensure the correct implementation of master integrals, various analytic and numerical checks were performed against published results in the literature as detailed in section 4. We have used the program CHAPLIN [79] for the numerical evaluation of the necessary harmonic polylogarithms in the physical region. The agreement of the poles of the one-and two-loop virtual amplitudes, as predicted by ref. [80] was checked both analytically and numerically, at the implementation level. The NLO contribution was checked against the MCFM [81] implementation 3 . The double-real contributions were implemented as described in sections 6.1 and 6.2, and double checked against another fully differential parameterization. Because the two parameterizations have different double-real counterterms, the numerical results for the double hard, the single hard and the integrated triple collinear counterterm cross sections are individually different. Only the sum of these contributions is physical, which provides a strong numerical check of our two implementations.
In the following, we present indicatively some differential distributions of interest, including their factorization and renormalization scale dependence. Since we do not include the decay of the off-shell photons to leptons, or the single-resonant diagrams in this publication, we defer a more detailed phenomenological analysis to a future publication.
In what follows, we use the central grid of the MSTW08 parton distribution functions [82], ignoring the uncertainties due to PDFs and the strong coupling constant. The strong coupling constant is run at the appropriate QCD order while the electromagnetic coupling constant is set to its value at m Z , a(m Z ) = 1/132.34.
The total cross section depends on the virtualities of the off-shell photons, m 3 = p 2 3 , m 4 = p 2 4 . First, we set the virtualities of the two photons equal and study the scale uncertainty of the NLO and NNLO K-factors as a function of the common photon virtuality, in fig. 2. For photons that are widely off-shell, i.e. with m 3,4 > 10GeV, the NNLO corrections are at the per mille level and the NNLO scale uncertainty is reduced, implying a satisfactory perturbative convergence for the process. As the limit of on-shell photons is approached the LO cross-section blows up and so does its scale uncertainty. This is expected, since we do not impose any final-state cuts on the two photons.
Next we turn to differential distributions for unequal photon virtualities. We set In fig. 3, we present the rapidity distributions of the two photons at each order in α s . The transverse momentum distributions for the two photons can be seen in fig. 4. The uncertainty due to the renormalization and factorization scales is shown as shaded regions in the figures. The scales are kept equal and varied in the interval µ r = µ f ∈ [10,40] GeV. (7.2) We note that while the NLO contribution changes the shape of the transverse momentum distributions, an effect that is more pronounced in the high transverse momentum region, the NNLO contribution does not induce any further changes. The rapidity distributions at NNLO also follow closely the NLO pattern.
Off-shell diphoton production contributes as a background, along with Z pair production, to the Higgs boson measurements in the golden channel pp → H → In that case the invariant mass of the two photons must be in a window of several GeV around the Higgs mass of 125GeV. We therefore set the virtualities of the photons to m 3 = 91.188GeV and m 4 = 27GeV, to simulate one  on-shell and one off-shell Z boson. The invariant mass distribution of the photon pair, shown at fig. 5, has its peak in the 126GeV region. The NNLO contributions are overall very small, but induce a slightly more pronounced correction at the region around the peak, and further stabilise the perturbative prediction there.
The N-jet cross section is shown in fig. 6. We have implemented the anti-k T   fig. 6. We observe that there is a migration of events away from the 1-jet bin at NNLO. On the right panel we show the dependence of the size of the 0-, 1-and 2-jet bins as a function of p min T . In general the contribution of the NNLO cross section is at the percent level or lower.
Finally the transverse momentum and rapidity distributions of the leading jet when the jet algorithm is defined with p min T = 20GeV is shown in fig. 7. This observable starts at NLO and the NNLO corrections are seen to be small and negative.

Conclusions
We have computed the NNLO corrections to off-shell diphoton production at the large N F limit, as a first step towards a complete, fully differential NNLO computation of off-shell diboson production that is necessary for improving the simulation of backgrounds for Higgs production in the four-lepton channel at the LHC.
We have provided explicit analytic expressions for the necessary two-loop master integrals in terms of classical polylogarithms, using direct integration methods along the lines of ref. [51].
We have treated the double-real radiation with a direct subtraction method where all subtraction counterterms are analytically integrated, thanks to the factorized structure of the singular limits in this process. The approach described is currently restricted to N F -type contributions, but is independent of the specific, colourless Born-level final state.
We have implemented the NNLO corrections in a fully differential partonic Monte Carlo code and provided selected differential distributions, demonstrating the numerical stability of both the double virtual and the double real contributions, in anticipation of a complete diboson computation.
where H n denotes the Q-vector space spanned by all multiple polylogarithms of weight n, and we set H 0 = Q.
Moreover, H can be equipped with a coproduct, turning it into a Hopf algebra. In the following we refrain from giving a detailed discussion of the Hopf algebra structure and only concentrate on the essentials that we will need in the following. In a nutshell (and very loosely speaking), a coproduct is a linear map ∆ : H → H⊗H that preserves the weight and the algebra structure 4 . For example, for the classical polylogarithms and the ordinary logarithms we have ∆(log r) = 1 ⊗ log r + log r ⊗ 1 , The advantage of the coproduct lies in the fact that it allows one to decompose a multiple polylogarithm of a specific weight into pairs of lower weight objects, for which properties like functional equations are already known. In addition, this decomposition can be iterated H → H ⊗ H → H ⊗ H ⊗ H → . . ., allowing one to decompose the functions into more and more combinations of functions of lower weight (which we will consider 'simpler' in the following). In the following we denote the by ∆ n 1 ,...,n k the component of the the coproduct in H n 1 ⊗ . . . ⊗ H n k . For a multiple polylogarithm of weight n this decomposition naturally stops when the function has been decomposed into an n-fold tensor product of functions of weight one, i.e., ordinary logarithms for which all identities are known. This maximal iteration of the coproduct is known as the symbol map in the literature [85][86][87][88][89].
The coproduct also encodes information on the discontinuities and the derivatives of a function. More precisely, discontinuities are encoded in the first factor of the coproduct, while derivatives only act on the second factor [84],

A.2 Construction of the basis
We now exploit the concepts reviewed in the previous section to construct a basis of functions through which all the integrals presented in this paper can be expressed. The discussion follows very closely the discussion in ref. [51], so we will be brief and online outline the main steps. Either by analysing explicit results for the integrals or by analysing the singularities of the differential equations satisfied by the master integrals, we see that the symbols of the master integrals have all their entries drawn from the set where u = rr, v = (1 − r)(1 −r) and w were defined in Section 4. Note that A 4 contains a subset A 3 = {r,r, 1 − r, 1 −r, r −r}, which corresponds to the case of the three-point functions considered in ref. [51]. Moreover, Cutkosky's rules imply that the Feynman integrals considered in this paper can only have branch cuts starting at point where u, v, w = 0. Let from now on H denote the Hopf algebra of all polylogarithmic functions whose symbols have all their entries drawn from the set A 4 , and H its subalgebra consisting of all functions having at most the branch cuts prescribed by Cutkosky's rule. Note that H is manifestly graded by the weight. Our goal is to find for every weight n (up to weight four) a basis for H n . In addition, we require this basis to be as 'simple as possible', i.e., we require that the product of two basis functions of weight m and n be an element of the basis of weight m + n.
A basis for H can now be constructed recursively in the weight. Indeed, since we know from eq. (A.4) that discontinuities are encoded in the first entry of the coproduct, we conclude that 5 Equation (A.6) is known as the first entry condition [86]. In the rest of this section we discuss how the first entry condition can be used to construct a basis for H recursively in the weight, following the procedure of ref. [51] (see also ref. [90]). Let us start with weight one. It is easy to see that a basis for H 1 is given by and a basis for the subspace H 1 is Next, we want to construct a basis for H 2 . From eq. (A.6) we know that and it is clear that a basis for H 1 ⊗ H 1 is given by However, not every element of H 1 ⊗ H 1 corresponds to a function in H 2 . Let us illustrate this with an example. Consider the element log(rr) ⊗ log r ∈ H 1 ⊗ H 1 , and suppose that there is a function f ∈ H 2 such that ∆ 1,1 (f ) = log(rr) ⊗ log r. Using eq. (A.4) and the fact that for the total differential d 2 = 0, we obtain a contradiction, because It can however be shown that This is known as the integrability condition. We can thus write down the most general linear combination of elements in B 1,1 and then solve the integrability condition. The a basis for the solution space of this problem is at the same time a basis for ∆ 1,1 (H 2 ). Every basis element of corresponds to a basis element in H 2 , and it is straightforward to find the corresponding function. Note that we also need to add all those elements ξ such that ∆ 1,1 (ξ) = 0. In our case there is just one such element, namely ζ 2 . Carrying out this procedure at weight two, we find that, besides all possible products of elements of B 1 , there are 4 new basis elements of weight two, which we choose as in agreement with the result quoted in section 4. This procedure immediately carries over to higher weight. Indeed, assume that we have constructed a basis B n−1 of H n−1 . The first entry condition and the integrability conditions imply that and a basis for H n−1 ⊗ H 1 is given by Starting from the most general linear combination of elements in B n−1,1 , we can solve the integrability condition and obtain a basis for ∆ n−1,1 (H n ). To every basis element corresponds a function in H n that can easily be constructed. We find that (up to weight four) the basis elements can be expressed in terms of multiple polylogarithms G( a; x) with functions considered here. In ref. [51] a basis up to weight four for these threepoint functions was constructed. Our basis has been chosen such that all basis functions of ref. [51] appear explicitly as basis elements in our basis.

B. Computation of the master integrals
In this appendix we illustrate how we computed the four-point master integrals defined in section 4. The method used to compute the integrals follows the algorithm introduced in ref. [52] (see also ref. [53][54][55][56][57][58][59]), which, under certain conditions which are always satisfied in the following, allows to perform the integrations one at the time.
In a nutshell, the general procedure is the following. If the integral is finite as → 0, we can expand the Feynman parameterized integral in under the integration sign. At each order in we then obtain integrands composed of (logarithms of) rational functions of the Feynman parameters and the external scales. If in addition we find an ordering of the Feynman parameters such that, after integrating over the first k Feynman parameters, all the polynomials appearing in the integrand are linear in the next Feynman parameter, then we can perform the next integration trivially using the definition of multiple polylogarithms, eq. (4.3). Several explicit algorithms to perform the integrations exist [52][53][54][55][56][57][58][59], and we refer to literature for the details. In the following we content ourselves to discuss the example of the four-point function B 2a defined in section 4.

B.1 A representative example: the integral B 2a
Let us illustrate the algorithm on the representative example of the integrals B 2a , corresponding to the integral B 2a = e 2γ d d k d d l (iπ d/2 ) 2 1 (k + l) 4 (k + p 1 ) 2 (k + p 1 + p 3 ) 2 (k + p 1 + p 3 + p 4 ) 2 l 2 (B.1) The integral over l is a massless bubble integral and can be done in closed form (B.2) After integration over the bubble, we obtain effectively a one-loop box integral where one of the propagators is raised to an -dependent power. Note that this applies to all the two-loop four-point integrals considered in section 4 and is not specific to B 2a . After Feynman parameterization of the remaining one-loop integral, we are left with the following integrals to compute In the case of B 2a , the propagator between the legs p 1 and p 2 is raised to the power 1 + , fp(1 + , 1, 1, 1, 4 − 2 ) = (−1) 2− Γ(2 + 2 ) Γ(1 + ) (B.5) Next, we would like to compute this remaining integral using the algorithm outlined at the beginning of this section. The integral is, however, divergent in d = 4 dimensions, and so are cannot naively expand in under the integration sign, but we first need to extract all the singularities. We first describe our method to extract the singularities, and then we illustrate the aforementioned algorithm on the resulting finite integrals.
Extraction of the singularities. The integral contains overlapping singularities that need to be factorized. After all singularities are factored, we can expand the singular terms in the integrand in terms of plus-distributions, obtaining a set of finite integrals that can be expanded in under the integration sign. In order to extract the singularities, we use the method of non-linear mappings introduced in ref. [78], which we review in the following. We start by considering the mapping, where the A j are constants. We then obtain for the integrand in eq. (B.5), (B.7) It is possible to remove all the kinematical dependencies from the denominator by solving the system of equations [93] A 2 A 4 = 1/s, A 1 A 3 = 1/t, . (B.10) The δ distribution can for example be solved by change of variables x 3 = y 1 , x 2 = (1 − y 1 )y 2 , x 4 = (1 − y 1 )(1 − y 2 )y 3 , x 1 = (1 − y 1 )(1 − y 2 )(1 − y 3 ), (B.11) where the Jacobian of the transformation is (1 − y 1 ) 2 (1 − y 2 ). Writingȳ i = 1 − y i we arrive at (−1) 2− Γ(2 + 2 ) Γ(1 + ) (ȳ 1 y 2ȳ2 y 3 + y 1 ) 2+2 , (B.12) where for the moment we did not apply the change of variables in the sum ( x i A i ) 3 for better readability. We obtain overlapping singularities that can be factorized completely by the following non-linear mapping y 1 → y 1 y 2 (1 − y 2 )y 3 y 1 y 2 (1 − y 2 )y 3 + (1 − y 1 ) .
The Jacobian is cancelled entirely and we end up with a integral free of overlapping singularities. Putting everything together, we obtain, − 2ȳ 3 (ȳ 1 y 2 A 2 + y 1 y 2ȳ2 y 3 A 3 +ȳ 1ȳ2 y 3 A 4 +ȳ 1ȳ2ȳ3 A 1 ) 3 , (B.14) where the A i 's are given in eq. (B.9). Substituting the functions for the A's (B.9) and trading the invariants t, m 2 3 , m 2 4 for the variables u, v, w we finally obtain the following representation for B 2a (writing the y i again as x i ), The two singularities are located in the variables x 2 and x 3 in a factorized form as intended. We then perform the expansion in with the help of the plus-distribution, i.e. by substituting and we obtain four finite integrals The sum of the four integrals represents the integral in eq. (B.15) up to order O( ). As each of these integrals is finite, they can be computed using the algorithm outlined at the beginning of this section. This will be illustrated in the rest of this section.
Doing the integrals. Let us do some of the integration explicitly to give a taste of the integration using multiple polylogarithms. The integral I 2a [1] is trivial and can be integrated directly without having to expand the integrand in . Also the integration over x 1 in I 2a [2] can be performed without any trouble, but let us not do this for the sake of illustration. The coefficient of 0 of I 2a [2] is given by where the dependence on x 1 dropped out and we are left with the integration over x 3 . The integrand can be written in terms of multiple polylogarithms and we can readily integrate over x 3 using the definition of multiple polylogarithms, eq. (4.3). We obtain All the other integrals can be done in this manner.