Antenna subtraction at NNLO with hadronic initial states: double real radiation for initial-initial configurations with two quark flavours

The antenna subtraction formalism allows to calculate QCD corrections to jet observables. Within this formalism, the subtraction terms are constructed using antenna functions describing all unresolved radiation between a pair of hard radiatior partons. In this paper, we focus on the subtraction terms for double real radiation contributions to jet observables in hadron-hadron collisions evaluated at NNLO. An essential ingredient to these subtraction terms are the four-parton antenna functions with both radiators in the initial state. We outline the construction of the double real subtraction terms, classify all relevant antenna functions and describe their integration over the relevant antenna phase space. For the initial-initial antenna functions with two quark flavours, we derive the phase space master integrals and obtain the integrated antennae.


Introduction
Jet production processes constitute an important tool for precision studies due to their large cross sections at high energy colliders. Reliable theoretical predictions for these observables require the calculation of at least the next-to-leading order QCD corrections. For these observables, the inclusive cross section with two incoming hadrons H 1 , H 2 can be written as dσ = a,b where ξ 1 and ξ 2 are the momentum fractions of the partons of species a and b in the incoming hadrons, f being the corresponding parton distribution functions and dσ ab (ξ 1 H 1 , ξ 2 H 2 ) is the parton-level scattering cross section for incoming partons a and b.
The partonic cross section dσ ab has a perturbative expansion in the strong coupling α s such that theoretical predictions for a hadronic process at a given order in α s are obtained when all partonic channels contributing to that order of the partonic cross section are summed and convoluted with the appropriate parton distribution functions as in eq. (1.1).
In general, beyond the leading order, each partonic channel contains both ultraviolet and infrared (soft and collinear) divergences. The ultraviolet poles are removed by renormalisation in each channel. Collinear poles originating from the radiation of initial state partons are cancelled by mass factorisation counterterms and absorbed in the parton distribution functions. The remaining soft and collinear poles cancel among each other when all partonic contributions are summed over [1]. As jet observables depend in a nontrivial manner on the experimental criteria used to define them, they can only be calculated numerically. The computation of hadronic observables including higher order corrections therefore requires a systematic procedure to cancel infrared singularities among different partonic channels before any numerical computation of the observable can be performed.
For the task of next-to-leading order (NLO) calculations, the infrared divergences present in real radiation contributions can be systematically extracted by process-independent procedures, called subtraction methods. The purpose of any subtraction method at NLO is to provide a subtraction term which has the same singular behaviour as the real radiation squared matrix element and is sufficiently simple to be integrated analytically over the radiation phase space which has been factorised from the (m + 1)-particle phase space. The actual form of this subtraction term depends on the subtraction formalism used. Several successful subtraction formalisms have been proposed in the literature [2][3][4][5][6]. Most notably, the FKS [3] subtraction by Frixione, Kunszt and Signer and the dipole formalism of Catani and Seymour [2] have been implemented in an automated way, the former in [7], the latter in [8][9][10][11][12][13]. The major challenge for NLO calculations is the computation of oneloop amplitudes for multiparticle processes. The evaluations of 2 → 4 processes at the next-to-leading order represent the current frontier [14][15][16][17].
Nevertheless, for some hadronic processes, in particular 2 → 1 or 2 → 2 scattering processes such as Drell-Yan, Higgs production, dijet production, vector-boson plus jet, vector-boson pair production or heavy quark pair production, the accuracy of the next-toleading order predictions is not sufficient to match the anticipated experimental accuracy, expected to reach the order of a few percent or better. Accurate precision studies enabling the extraction of fundamental parameters of the theory will require that the theoretical predictions have the same precision. Those need therefore to be evaluated up to the nextto-next-to-leading order (NNLO) in perturbative QCD .
The calculation of observables with m jets in addition to other objects (like for example vector bosons) at the NNLO requires three distinct contributions: the double real radiation dσ R N N LO with (m + 2) final state partons, the mixed real-virtual radiation dσ V,1 N N LO with (m + 1) final state partons and the 2-loop virtual radiation dσ V,2 N N LO , with m final state partons. Those build the NNLO cross section which is given by 2) The individual contributions in the m-, (m + 1)-and (m + 2)-parton final states are all separately infrared divergent. After renormalisation and factorisation, their sum is finite, though. For most massless jet observables of phenomenological interest, the two-loop matrix elements have been computed some time ago, while the one-loop matrix elements are usually known from calculations of NLO corrections to (m + 1)-jet production [18] . The one-loop and two-loop matrix-elements contain explicit infrared divergences from the loop integration. Those cancel with divergences which are implicit in the real radiation for (m + 1)-and (m + 2)-parton processes. These real radiation divergences become explicit only once the phase space integration is carried out. The main issue of these calculations is therefore to find a method to extract and cancel the infrared divergences among these three contributions in order to finally evaluate numerically the finite remainders to obtain the NNLO contribution to the cross section. Like at NLO, a subtraction formalism is needed in order to extract the infrared divergences from these contributions. At parton level, the general form of the cross section for an m-particle final state at NNLO including subtraction terms is given by [19]: Here, dσ S N N LO denotes the subtraction term for the (m + 2)-parton final state which behaves like the double real radiation contribution dσ R N N LO in all singular limits. Likewise, dσ V S,1 N N LO is the one-loop virtual subtraction term coinciding with the one-loop (m + 1)particle contribution dσ V,1 N N LO in all singular limits. The two-loop correction to the mparton final state is denoted by dσ V,2 N N LO . In addition, when there are partons in the initial state, there are two mass factorisation contributions, dσ M F,1 N N LO and dσ M F,2 N N LO , for the (m + 1)-and m-particle final states respectively. Like at the next-to-leading order level, the subtraction terms are needed in their unintegrated as well as in their integrated forms.
For processes with initial-state partons, the antenna subtraction formalism has been so far fully worked out only to NLO in [59]. It has been extended to NNLO for processes involving one initial state parton relevant for electron-proton scattering in [60] while an extension of the formalism to include two initial state hadrons at NNLO is under construction [61,62]. An essential step towards this aim is performed in [62] where an explicit derivation of the subtraction terms needed for the double real contributions to the sixgluon process is presented. The general structure of the unintegrated subtraction terms relevant for the double real contributions to any hadronic observables evaluated at NNLO is presented there as well.
In this paper, we will focus on the integrated form of the subtraction term relevant for double real radiation for processes involving two partons in the initial state. In the antenna subtraction formalism, the subtraction terms are built with so-called antenna functions. The latter describe all unresolved partonic radiation off a hard pair of colour-ordered partons, the radiators. Depending on where the two hard radiators are located, three cases need to be distinguished: both radiators are in the final state (final-final), only one radiator parton is in the initial state (initial-final) or both radiator partons are in the initial state (initial-initial). The subtraction terms and the antenna functions building them are separated corresponding to these three cases. In the most general hadronic process (two partons in the initial state, two or more partons in the final state), all three configurations have to be taken into account.
As discussed in [19,60,62] the subtraction terms are separated according to the colour connection of the unresolved partons. In this paper, we will specialise on the subtraction term for the case where two unresolved (soft or collinear) partons are colour connected to the two incoming partons. This is indeed the only case where new ingredients, namely the four-parton initial-initial antennae are needed in unintegrated as well as in integrated forms. On a longer term we are aiming to evaluate the whole set of integrated four-parton initial-initial antennae, those are universal building blocks for the subtraction terms for any hadronic process evaluated at NNLO.
In a first step towards this aim, in this paper, we have focused on the crossings of two partons from a set of three 4-parton final-final antennae involving two quark flavours. More precisely, the paper will be organised as follows: In Section 2, we present the general formulae for the subtraction terms related to double real radiation for initial-initial configurations while the colour-connected case is treated explicitly in Section 3. Section 4 establishes a list of all non-identical initial-initial four-parton antenna functions relevant to construct the subtraction term for the double real radiation off two initial-state partons. In Section 5, the phase space mapping appropriate for initial-initial configurations is presented.
Finally, Section 6 contains our results for the integrated initial-initial antenna functions with two quark flavours and Section 7 our conclusions. Since the results are lengthy, we show only the leading pole terms in the manuscript, and attach the complete results as a Mathematica file.

Antenna subtraction for double-real radiation at NNLO in the initialinitial configuration
Antenna subtraction has been derived explicitly for final-final and initial-final configurations at NLO and NNLO in [19] and in [59,60] respectively. For the initial-initial case, it has been derived at NLO in [59] and is under construction at NNLO. At NLO, the subtraction term is introduced to extract and cancel the infrared divergences present in the real contributions. The general forms of the subtraction terms required at NLO in any of the three configurations (final-final, initial-final and initialinitial) have been given in [19,59] and summarised in [62]. At this order, only tree-level three-parton antenna functions involving one unresolved parton are needed to build the subtraction terms. Those functions are usually denoted by X 0 ijk , X 0 i,jk , X 0 ik,j in the three configurations. Their definitions will be recalled in Section 4.
At NNLO, two contributions, double real and mixed real-virtual, require the introduction of a subtraction term. For final-final and initial-final configurations both types of subtraction terms have been constructed; at this order tree-level four-particle antennae involving two unresolved partons and one-loop three-parton antenna functions are needed respectively. Those functions have been derived and integrated over the corresponding factorised phase space in [19] and [60]. For initial-initial configurations, the general structure of the unintegrated subtraction terms relevant for the double real contributions to pp → m jets is presented in [62], we will recall it in this section.
The double real radiation contribution to the m-jet cross section in pp collisions reads (2.1) In this equation, |M m+2 (p 1 , p 2 ; k 1 , . . . , k m+2 )| 2 stands for the colour-ordered 2 → m + 2 matrix-element squared. The symmetry factor S m+2 accounts for identical partons in the final state. The normalisation factor N includes all QCD-independent factors as well as the dependence on the renormalised QCD coupling constant α s . m+2 denotes the sum over all configurations with m + 2 partons. The initial state momenta are labelled as usual as p 1 and p 2 whereas the m + 2 momenta in the final state are labeled k 1 , . . . , k m+2 . dΦ m+2 is the 2 → m + 2 particle phase space where we have introduced the abbreviation [dk 1 ] = (d d−1 k 1 )/(2E 1 (2π) d−1 ). The jet function J (m+2) m (k 1 , . . . , k m+2 ; p 1 , p 2 ) ensures that out of (m + 2) final state partons, an observable with m jets is built. The incoming parton momenta p 1 , p 2 serve as reference directions to define transverse momenta and rapidities of the jets.
The double real radiation contribution given in eq. (2.1) can become singular if either one or two final state partons are unresolved (soft or collinear). Consequently, when constructing the corresponding subtraction term dσ S N N LO in eq. (1.3), which shall correctly subtract all those single and double unresolved singularities we must distinguish the following configurations according to the colour connection of the unresolved partons: (a) One unresolved parton but the experimental observable selects only m jets.
(c) Two unresolved partons that are not colour connected but share a common radiator (almost colour-unconnected).
(d) Two unresolved partons that are well separated from each other in the colour chain (colour-unconnected).
(e) Large angle soft gluon radiation.
This separation among subtraction contributions according to colour connection is valid in final-final, initial-final or initial-initial configurations and in any of those the subtraction formulae have a characteristic structure in terms of required antenna functions. This antenna structure has been derived for the final-final and initial-final cases in [19,60] and in [62] for the initial-initial case.
In here, we focus only on the initial-initial case with the kinematical situation where two unresolved partons are colour-connected to the two incoming partons. This is the only case where new ingredients are needed, namely the initial-initial four-parton antenna functions denoted by X il,jk . The four-particle initial-initial antenna functions will be defined explicitly in Section 4.

Subtraction terms for two colour-connected unresolved partons in the initial-initial configuration
When two unresolved partons j and k are adjacent and colour-connected to two initialstate partons, the subtraction term related to the double real contribution dσ R N N LO given in eq. (2.1) reads : where the sum runs over all colour-adjacent pairs j, k and implies that the hard momenta i, l are chosen accordingly. X 0 il,jk denotes a four-particle tree-level initial-initial antenna function defined explicitly in Section 4. By construction those contain all colour connected double unresolved limits of the 2 → m + 2 parton matrix element associated with partons j and k being unresolved between radiators i and l. However this antenna can also become singular in single unresolved limits associated with j or k only where it does not coincide with limits of the matrix element. To ensure that this subtraction term is only active in the double unresolved limits of the real matrix elements squared we remove these single unresolved limits of the four-particle antennae. Those limits are products of two tree-level three-particle antennae, namely products of an initial-final antenna X a,bc and an initial-initial antenna X Ad,C . In these antennae we have replaced the original hard radiators with new particles,Î andL. When both radiators are in the initial state as it is the case here, pÎ = x i p i , pL = x l p l . The product of these three-particle antenna functions in dσ S,b,(ii) N N LO then subtracts the single unresolved limits of the associated four-particle antenna.
The 2 → m matrix element |M m (pÎ ,pL;k 1 , . . . ,k i ,k l , . . . ,k m+2 )| 2 is evaluated with new on-shell momenta which are Lorentz boosted as a result of the mapping required to ensure factorisation of matrix element and phase space. The NNLO initial-initial mappings have been discussed in [59] and will be recalled in this paper in Section 5.
Using a factorised form for both matrix element and phase space as shown in figure 1, one is able to obtain an integrated form of the subtraction term with m partons in the final state which can be combined with the virtual two-loop contributions having also m final state partons as defined in eq. (1.3).
More explicitly, the factorisation of the phase space with m + 2 final state particles denoted by k 1 , . . . , k m+2 , reads, Using this factorised form of the phase space we can rewrite the integrated (colourconnected) subtraction term involving only the four-parton antennae in the form, Moreover, the integrated antennae are defined as the antenna functions integrated over the antenna phase space as defined in eq. (3.2) including a normalisation factor to account for powers of the QCD coupling constant by, where C(ǫ) is given by, These integrations are performed analytically in d dimensions to make the infrared singularities explicit. The integrated initial-initial antennae are presently unknown and the calculation of a sub-set of those involving two quark flavours is the subject of Section 6 and constitutes our main result in this paper. A list of all non-identical initial-initial four-particle antennae will be presented in Section 4.

Initial-initial antenna functions
In this section we shall recall how the antenna functions are defined in any of the configurations and how they enter in the construction of the subtraction terms. We also present a list of all non-identical four-parton initial-initial antenna functions.

Definitions
In the antenna formalism, in any of the three configurations (final-final, initial-final or initial-initial) the subtraction terms are constructed from products of antenna functions with reduced matrix elements (with fewer final state partons than the original matrix element). The integrated subtraction terms are obtained after an integration over a phase space which is factorised into an antenna phase space (involving all unresolved partons and the two radiators of the antenna) multiplied with a reduced phase space (where the momenta of radiators and unresolved radiation are replaced by two redefined momenta). These redefined momenta can be in the initial state, if the corresponding radiator momenta were in the initial state as we saw in Section 3 in eq. (3.1) for the subtraction term σ An antenna function is determined by the external states it contains and the pair of hard partons it collapses to in the unresolved limits. In general we denote the antenna function as X. For antennae that collapse onto a hard quark-antiquark pair, X = A for qgq. Similarly, for a quark-gluon antenna, we have X = D for qgg and X = E for qq ′q′ final states. Finally, we characterise the gluon-gluon antennae as X = F for ggg, X = G for gqq final states.
At NLO, we only need to consider tree level three-particle antennae involving only one unresolved parton. At NNLO we will need four-particle antennae involving two unresolved partons and one-loop three-particle antennae.
In all cases the antenna functions are derived from physical matrix elements associated to the decay of a colourless particle into partons: the quark-antiquark antenna functions are derived from γ * → qq + (partons) [63], the quark-gluon antenna functions fromχ → g + (partons) [64] and the gluon-gluon antenna functions from H → (partons) [65]. The tree-level antenna functions are obtained by normalising the three-and four-parton treelevel colour sub-amplitudes squared to that of the basic two-parton process: The final-final three-and four-particle antennae are respectively defined by: where S denotes the symmetry factor associated with the antenna, which accounts both for potential identical particle symmetries and for the presence of more than one antenna in the basic two-parton process. It is chosen such that the antenna function reproduces the unresolved limits of a matrix element with identified particles. The initial-final tree-level three-and four-parton antennae denoted by X 0 i,jk and X 0 i,jkl are in principle obtained by crossing one-parton to the initial state starting from the corresponding final-final antennae. However this crossing might be ambiguous as was first noticed in [59] for the quark-gluon type antenna. Depending on the unresolved limit considered, the pair of hard partons it collapses to may be different.
The initial-initial tree-level three-and four-parton antennae denoted by X 0 ik,l and X 0 il,jk are obtained by crossing two partons to the initial state, starting from the final-final antennae. This crossing procedure, unlike in the initial-final case, is free of ambiguity as the pair of hard partons the initial-initial antenna collapses to is always uniquely defined. More explicitly, the three-parton initial-initial antenna function X 0 ik,j is defined as where i denotes that particle i is crossed to the initial state. Therefore we can see that the initial-initial antenna is connected to the final-final antenna (where all coloured particles are outgoing) by with ∆ F the difference in the number of fermion crossings between the three-particle and the two-particle subamplitude. The tree-level four-particle initial-initial antenna X 0 il,jk is defined as and the relation to the final-final antenna is (4.5) The four-particle tree-level antenna functions are not determined by the species of the particles alone but also by the colour-connection. We distinguish leading-colour antennae, denoted by letters without tilde, where the particles are colour-connected in the order they are listed and subleading colour antennae, denoted by letters with tilde, where the gluons are photon-like. This notation has been used in [19,60] and we will further use it in this section in the tables below.
The unresolved limits of the initial-initial antennae can be obtained from those of the final-final antennae by crossing. The crossing of the triple-collinear splitting functions is explained in [66]. Those limits will be reported elsewhere.

Catalogue of tree-level four-particle initial-initial antenna functions
Any two particles of a four-particle final-final antenna can be crossed to the initial state to obtain an initial-initial antenna; therefore one final-final four-particle antenna gives rise to six initial-initial antennae. Due to symmetries, at most four of these initial-initial antennae are different. The independent crossings are listed below. To make the colour connection clear, in this list we write out the arguments of the antennae explicitly, i.e we write X 0 4 i, j, k, l (where i denotes an incoming particle) instead of X 0 ik,jl . Some initialinitial antennae are free of singular limits. These finite antennae are not needed for the construction of subtraction terms, but their integrated form could be needed for crosschecks. They are marked with * * below.

quark-antiquark antennae
is symmetric under the interchange of the two photon-like gluons. The nonidentical-flavour antenna B 0 4 is separately symmetric under interchange of q ′ with q ′ and of q with q.

quark-gluon antennae
Due to the cyclic colour connection, D 0 4 is symmetric under interchange of the second and fourth gluon.
is symmetric under cyclic interchange of its arguments. G 0 4 is symmetric under the interchange of the two gluons as well as under the interchange of q with q. H 0 4 has three symmetries, q ↔ q, q ′ ↔ q ′ and the flavour renaming q ↔ q ′ , q ↔ q ′ .

Phase space factorisation and mappings
The construction of subtraction terms requires a mapping from the original set of momenta onto a reduced set. The mapping interpolates between the different soft and collinear limits which the subtraction term regulates. An appropriate mapping for the initial-initial case, both for single and double unresolved configurations, has been discussed in [59]. By requiring momentum conservation and phase space factorisation, the phase space mapping is strongly constrained. The remapping of initial state momenta can only be a rescaling, since any transversal component would spoil the phase space factorisation. For two unresolved partons j and k, a complete factorisation of the phase space into a convolution of an m-particle phase space depending on redefined momenta only and the phase space of the unresolved partons j and k can be achieved with a Lorentz boost. This boost maps the momentum q = p 1 + p 2 − k j − k k , with q 2 > 0 and p 1 , p 2 being the momenta of the hard emitters, into the momentumq = x 1 p 1 + x 2 p 2 , where x 1 and x 2 are fixed in terms of the invariants as follows: These two definitions guarantee the overall momentum conservation in the mapped momenta and the correct soft and collinear behaviours . The two momentum fractions x 1 and x 2 satisfy the following limits in double unresolved configurations: 1. j and k soft: 2. j soft and k k = z 1 p 1 : 3. k j = z 1 p 1 and k k = z 2 p 2 : and all the limits obtained from the ones above by the exchange of p 1 with p 2 and of k j with k k . The construction of NNLO antenna subtraction terms also requires that all single unresolved limits of the four-parton initial-initial antenna functions X il,jk , with radiators i and l, have to be subtracted, such that the resulting subtraction term is active only in its double unresolved limits. A systematic subtraction of these single unresolved limits by products of two three-parton antenna functions can be performed only if the NNLO phase space mapping turns into an NLO phase space mapping in its single unresolved limits. A detailed discussion of the corresponding translation between these two momentum mappings can be found in [59].
The factorisation of the (m + 2)-parton phase space into an m-parton phase space and an antenna phase space involving the unresolved partons j and k given in eq. (3.2) can equivalently be written as where J is the Jacobian factor defined by This phase space parametrization can also be used to give an equivalent definition of the integrated initial-initial antennae first given in eq. (3.4). Those are given as, where We shall use this definition of the integrated antennae X 0 il,jk (x 1 , x 2 , ε) to compute them in Section 6.

Integration of the four-parton initial-initial antennae
In the first part of this section we describe how the initial-initial four-parton antennae are integrated over the antenna phase space while in the second part of this section we restrict ourselves to the evaluation of integrated antennae involving two quark flavours.

Calculational method
The initial-initial antenna functions have the scattering kinematics where q is the momentum of the outgoing colourless particle. The momenta satisfy: The four-parton initial-initial antennae defined in Section 4 need to be integrated over the phase space of the unresolved partons j and k. This integration yields a result which depends only on s 12 , x 1 and x 2 . From dimensional counting, one can immediately conclude that the dependence on s 12 is only multiplicative, according to the mass dimension of the integral.
The initial-initial antenna phase space integrals are derived from squared matrix elements and can be represented by forward scattering diagrams as in the following figure: The two delta functions in eq. (5.2) can be represented as mass-shell conditions of fake particles and are shown in the previous picture as a thick solid line (representing a massive particle with mass M 2 = x 1 x 2 s 12 ) and a dashed line (representing the other constraint) . This allows us to use the optical theorem to transform the initial-initial antenna phase space integrals into cut two-loop integrals and therefore use the methods developed for multiloop calculations [67][68][69] . Up to eight-propagator integrals with four cut propagators are generated in this way . Using the reduction techniques, the calculation of the integrated antennae can be related to the evaluation of a reduced set of master integrals . For the complete set of non-identical initial-initial four-parton antennae tabulated in Section 4, we find 32 such integrals, obtained using integration-by-part (IBP, [70]) and Lorentz invariance (LI, [71]) identities, following the Laporta algorithm [72]. A private implementation as well as a public one [73] have been used. For those four-parton antenna functions there are 13 different propagators, including the four that are cut in the phase space integration (D j , D k , D jk12 , D jk123 ): where p 3 = x 2 p 2 − x 1 p 1 . To perform the reduction to master integrals, we drop any integral where D j , D k , D jk12 , D jk123 are not in the denominator and impose momentum conservation. The integrands of the 32 master integrals found can all be written as rational polynomials of the denominators above. These master integrals are calculated using either the method of differential equations or by a direct evaluation of the phase space integrals in terms of hypergeometric functions. The simplest master integral is a two loop bubble with all the internal lines cut, it is obtained from eq. (5.3) by replacing the jacobian J and the antenna X 0 il,jk with unity: This integral can actually be expressed as a hypergeometric function: it has been checked against the master integral I[0] appearing in the calculation of the gauge-boson rapidity distribution at NNLO in [68], the notation translates as u = x 2 /x 1 , z = x 1 x 2 . The set of master integrals which we denote by I i (x 1 , x 2 , ε) are functions of x 1 , x 2 and ε.
We begin by factoring out the leading behavior of the master integrals I i (x 1 , x 2 , ε) in the limits x 1 → 1 and x 2 → 1, keeping the exact ε-dependence: The integers m 1 , m 2 are characteristic to each master integral. The functions F i (x 1 , x 2 , ε) are regular at x 1 = 1, at x 2 = 1, and at x 1 = x 2 = 1 and can be calculated as Laurent series with, at most, second order poles in ε. The integrated antennae given by X (x 1 , x 2 , ε) are linear combinations of these master integrals I i (x 1 , x 2 , ε), with coefficients containing poles in ε, as well as in (1 − x 1 ) and (1 − x 2 ). After the masters have been inserted into the integrated antennae, those take the form where R(x 1 , x 2 , ε) is regular at the boundaries x 1 = 1, x 2 = 1, and at x 1 = x 2 = 1. The ε-expansion of the singular factors (1 − x i ) −1−2ε is done in the form of distributions: To evaluate the integrated antennae, we decompose the phase space into four regions depending on the values of x 1 and x 2 . Those regions are given by: • x 1 = 1, x 2 = 1, which we refer to as the hard region • x 1 = 1, x 2 = 1, and x 1 = 1, x 2 = 1, referred to as collinear regions • x 1 = 1, x 2 = 1, which we denote the soft region .
In the hard region (x 1 = 1, x 2 = 1), harmonic polylogarithms of weight two appear in the O(ε 0 ) term of R. Therefore, the ε-expansion of the master integrals in the hard region is needed at least up to the order at which terms of transcendentality 1 two appear for the first time generally.
In the collinear regions (x 1 = 1 or x 2 = 1), since the expansion in distributions (6.6) generates additional 1/ε factors, the function R is required up to O(ε) where harmonic polylogarithms of weight 3 appear. The masters evaluated in the collinear region need therefore to be expanded at least up to the order at which terms of transcendentality 3 appear generally.
Finally, in the soft region (x 1 = x 2 = 1), since the expansion of the distributions (6.6) generates additional 1/ε 2 coefficients, the function R is required up to O(ε 2 ) where polylogarithms of weight 4 appear. The masters evaluated in the soft region need therefore to be expanded up to transcendentality 4 at least.

Integrated antennae with two quark flavours
In a first step towards the calculation of all integrated four-parton initial-initial antenna functions for the double real radiation case, in this paper, we focus on all the crossings of two partons from the following four-parton final-final antennae: B 0 4 (q, q ′ ,q ′ ,q),Ẽ 0 4 (q, q ′ ,q ′ , g) and H 0 4 (q,q, q ′ ,q ′ ) defined in [19] . We found that the reduction involving only those initial-initial antennae leads to 12 master integrals. Those without numerators are shown in Fig. 2 . We have performed the calculation of these integrated antennae with two choices of master integral bases differing by four master integrals. In the first basis, the definitions of the master integrals involved in the calculation are as follows: In the other choice of basis, the masters with scalar products in the numerator (I 1 , I 3 , I 4 , I 5 ) are replaced by alternative master integrals: The relation to the first basis is: Before we proceed with the details of calculating the masters, we present in Table 1 a summary of which regions contribute to the crossings of the antenna functions B 0 4 (q, q ′ ,q ′ ,q), E 0 4 (q, q ′ ,q ′ , g) and H 0 4 (q,q, q ′ ,q ′ ). We then present results of the required masters in each of these regions. As explained in Section 4.2, not all of the six crossings of a final-final antenna are different. Labelling the final-final antenna functions as B 0 4 (1q, 3q ′ , 4q ′ , 2q), E 0 4 (1q, 2q ′ , 3q ′ , 4g) and H 0 4 (1q, 2q, 3q ′ , 4q ′ ), the identical crossings are the following: 14) As mentioned in Section 4.2, B 34 is free of singular limits. It will therefore not be needed for the construction of subtraction terms and its integrated form is free of poles in ε. It might be needed for checks of the integrated antennae, however, which is why it is included here.

Master integrals
In the following, we present the results for the masters in the hard, collinear and soft regions while restricting ourselves to those which are explicitly involved in the calculation of the integrated antennae with two quark flavours.

a) The hard region
The master integrals defined above were computed in the hard region mainly with the differential equations technique [71,74,75]. The only masters that were calculated directly are I 1 , I 2 , and I 14 . Since the dependence of the integrands on x 1 and x 2 is via the constraints, C 1 and C 2 , as shown in eq. (5.4), we derive the differential equations for each master integral by employing the following operators at the integrand level: The boundary conditions required for the solution of the differential equations are either obtained from self-consistency conditions on the integrals, or by explicit evaluation in the collinear or soft limits. The solution of the system of differential equations yields twodimensional generalized harmonic polylogarithms (GHPL, [76,77]) of up to weight two, or products of weight one harmonic polylogarithms (HPL, [78]) of argument x 1 or x 2 . The definition of the HPL and GHPL functions involved in the solution of the master integrals in the hard region is recalled below: The harmonic polylogarithms of higher weight are defined recursively (we group weights into vectors, b 1 , . . . , b w = b): with weight functions This results in the derivative formula The two-dimensional generalized harmonic polylogarithms are defined in a very similar fashion: In the results for the hard region below, the GHPL's will have x 1 as their argument and 0, ±1 and −x 2 as their weights, the conversion to HPL's has been inserted where possible. An important check on our results in this region is to compare the soft and collinear limits of the hard result against the ones derived in the soft and collinear regions calculated by a direct integration of the phase space. We found agreement on all the powers of the εexpansion in the hard region using these two different methods. The results of the masters expanded to the needed order in ε are given below. The common prefactor N Γ is defined as: (2) )  (2) )   (2)) (2)) where The master I 14 was calculated using differential equations and is written in terms of logarithms and dilogorithms of complicated arguments instead of GHPL's and HPL's. The coefficients of the differential equations for I 14 involve a new denominator, cubic in x 1 and quadratic in x 2 , which does not occur in the differential equations for any of the other masters presented in this paper. Because of this denominator, I 14 has this more complicated form. The differential equations for I 14 are solved by first considering the corresponding homogeneous equations. A solution of these equations, for ε = 0, is given by Note that it is singular along the curve Q 3 = 0, part of which is inside the physical region 0 < x 1 < 1, 0 < x 2 < 1. On the other hand, the full solution for I 14 , which we write as cannot be singular inside the physical region. Therefore, the coefficient C 14 must vanish along the curve Q 3 = 0. In order to solve the differential equations for I 14 , it is convenient to consider them keeping x 12 = x 1 + x 2 fixed, so that Q 3 = x 1 (x 2 12 + 4) − 4 x 12 becomes linear in x 1 . Using the boundary condition that C 14 must vanish at the point x 1 = x 0 1 = 4x 12 /(x 2 12 + 4), where Q 3 = 0, the solution of the differential equation of C 14 in x 1 can then be written as a one-dimensional integral: The integration can be performed analytically in terms of logarithms and dilogarithms. The final result for I 14 is shown in eq. (6.43).
We have presented the result for I 14 in a form where all (di)logarithms are real in the region where Q 3 > 0 and s is real. By construction, I 14 has no singularity at Q 3 = 0. For Q 3 < 0, s is imaginary and the arguments of the (di)logarithms that depend on s become complex. Nevertheless, I 14 remains real in this region because the factor C 14 is now imaginary.
Finally, it is worth noting that I 14 is only needed up to finite order in the hard region. It is required up to order ǫ in the collinear x 1 region where it has been calculated separately using differential equations and direct integration.

b) The collinear regions
The real-real master integrals were calculated in the collinear regions by deriving a hypergeometric integral representation starting from their definition as phase space integrals. This allows us to give the results in closed form. We also have derived these results as an expansion in ε using the differential equations method, which provides an important check on our results. We list below the results of the masters contributing to the antennae B 0 4 ,Ẽ 0 4 , H 0 4 in these regions.
• collinear x 1 In this region, only I 1 , I 2 , I 7 and I 8 require to be expanded up to O(ε 3 ) since they contribute toẼ 12 which starts its power expansion at O(ε −3 ).
× (−(k j · p 1 )(k j · p 3 )) = I 3 = (s 12 ) 1−2ε N Γ 2 −3−ε (1 − x 1 ) 1−2 ε × (6.50) • collinear x 2 In this region, only the following four masters (all of them are needed up to O(ε 2 )) con-tribute: × k k · p 1 = I 1 = (s 12 ) 1−2ε N Γ 2 −3−ε (1 − x 1 ) 1−2 ε × (6.57) c) The soft region As discussed previously, the master integrals in the soft region are expected to be needed at most two orders in ε higher than in the hard region. It turns out here that the masters are only needed up to O(ε 3 ) at most in the soft region. The reason behind this is the absence of double soft gluon configurations in the antennae B, E and H. Furthermore, as can be seen from Table 1, only B 12 will have a contribution from the soft region at all. This is due to the fact that only B 12 allows for a singular double soft configuration, a soft q −q pair. The master integrals which are required in this case up to O(ε 3 ) are either I 2 or I ′

Conclusions
Within the antenna subtraction formalism, allowing the calculation of higher order QCD corrections to jet observables, subtraction terms are constructed from antenna functions. Those functions describe all unresolved radiation between a pair of hard radiator partons. At NNLO, this formalism has been fully developed and applied so far only for colourless initial states. In this paper, we have focussed on the extension of this formalism to evaluate NNLO corrections to jet observables at hadron colliders and concentrated on the construction of subtraction terms for the double real radiation contributions. More precisely, we have considered the subtraction terms needed to account for the radiation of two colour-connected unresolved partons off two initial state partons. For these subtraction terms, four-parton tree-level initial-initial antenna functions are required in unintegrated and integrated form. The integration over the phase space associated with two unresolved partons has to be performed analytically.
In this paper, we have given a catalogue of all non-identical four-parton initial-initial antenna functions. Furthermore, after applying standard reduction techniques, we found that 32 master integrals are necessary to obtain their integrated form. As a step towards the integration of the full set of integrated initial-initial antenna functions, in this paper we have focussed on the initial-initial antennae obtained from crossing two partons in the final-final antenna functions characterised by the presence of two quark flavours. After reduction, 12 masters were required to obtain those. We presented the decomposition of the calculation according to four phase space regions: hard, collinear and soft and we gave the master integrals in these regions. Finally, we presented the results for those integrated initial-initial four parton antenna functions themselves. Since the results are lengthy, we have shown only the leading pole terms in the manuscript and have attached the complete results as a Mathematica file.

Acknowledgements
This work is supported in part by the U.S. Department of Energy, Division of High Energy Physics, under contract DE-AC02-06CH11357 and by the Swiss National Science Foundation (SNF) under contract PP0022-118864.