The Two-Loop Master Integrals for $q \bar{q} \to V V$

We compute the full set of two-loop Feynman integrals appearing in massless two-loop four-point functions with two off-shell legs with the same invariant mass. These integrals allow to determine the two-loop corrections to the amplitudes for vector boson pair production at hadron colliders, $q \bar{q} \to V V$, and thus to compute this process to next-to-next-to-leading order accuracy in QCD. The master integrals are derived using the method of differential equations, employing a canonical basis for the integrals. We obtain analytical results for all integrals, expressed in terms of multiple polylogarithms. We optimize our results for numerical evaluation by employing functions which are real valued for physical scattering kinematics and allow for an immediate power series expansion.


Introduction
Precision studies of the electroweak interaction at the LHC are based on a wealth of observables derived from vector boson pair production, γγ, Zγ, W γ, ZZ, W W , W Z, which allow to test the SU (2) L × U (1) Y gauge structure and the field content of the Standard Model. Anomalous contributions to these interactions can probe physics beyond the Standard Model at energy scales well beyond direct searches. To fully exploit these observables, precise theoretical predictions are of crucial importance, including especially higher order perturbative corrections. To reach an accuracy in the per-cent range, thus matching the experimental precision at the LHC, corrections to next-to-leading order (NLO) in the electroweak theory and to next-to-next-to-leading order (NNLO) in QCD are to be included.
The full set of NLO QCD corrections [1][2][3][4] and large parts of the NLO electroweak corrections [5][6][7][8][9][10] have been derived for vector boson pair production. NNLO QCD corrections were calculated up to now for γγ [11] and Zγ [12] production. Key ingredient to the NNLO calculations are the two-loop matrix elements for qq → V 1 V 2 , which are known for γγ [13,14] and V γ [15,16] production, and for qq → W W in the high energy approximation [17]. The full calculation of two-loop matrix elements for the production of two massive vector bosons is still an outstanding task, and requires the derivation of a new class of two-loop Feynman integrals: two-loop four-point functions with internal massless propagators and two massive external legs. First results on these were obtained already, with the derivation of the full set of planar two-loop integrals for vector boson pairs with equal mass [18] and two different masses [19]. In the present paper, we extend our earlier calculation [18] to compute the full set of two-loop four point functions relevant to vector boson pair production with two equal masses. A subset of these, the three-point functions with three massive external legs and all massless internal propagators, has already been known for some time in the literature [20,21]. The problem kinematics and notation are described in section 2. Working in dimensional regularization with d = 4 − 2 space-time dimensions, we identify the relevant master integrals (MI) and derive differential equations for them employing integrationby-parts (IBP) [22,23] and Lorentz-invariance (LI) [24] reductions through the Laporta algorithm [25] implemented in the Reduze code [26,27]. The master integrals are then determined by solving these differential equations [24,[28][29][30] and matching generic solutions to appropriate boundary values obtained in special kinematical limits. Improving upon our earlier results [18], we are now transforming the differential equations to a canonical form [31] which renders their integration trivial after an expansion in . The algorithm applied for this transformation is described in detail in section 3, similar procedures have been put forward most recently in [32,33]. With this, the remaining non-trivial step in the calculation of the master integrals is the determination of the boundary terms, which we describe in section 4. We use a similar setup for our parametrization and treatment of functions as in the first calculation [34] of non-planar double boxes with this type of external kinematics. In particular, our solutions are described in terms of multiple polylgarithms. We fix the boundary terms of all complicated integrals by imposing a simple set of regularity conditions. The implementation of these conditions and further processing of the multiple polylogarithms relies on computer-algebra implementations of the coproduct augmented symbol formalism [35][36][37][38][39] and other techniques, which we described in [18,40]. Section 5 contains a discussion of our solutions and the checks we performed on them. In section 6 we describe the final form of our analytical results in terms of a particular set of real valued Li 2,2 , Li n (n = 2, 3, 4) and ln functions, optimized for fast and stable numerical evaluations. We complement our exact results by expanding them at the production threshold and in the high-energy region. We conclude in section 7 and specify the exact definition of our canonical basis in appendix A. For all algebraic manipulations we made extensive use of FORM [41] and Mathematica [42].

Notation and reduction to master integrals
We consider the production of two vector bosons V of mass m in the scattering kinematics: where p 2 1 = p 2 2 = 0 and q 2 1 = q 2 2 = m 2 . The Mandelstam invariants are so that in the physical region relevant for vector-boson pair production we have: We choose to work with dimensionless variables x, y and z defined by The Landau variable x absorbs a square root s(s − 4m 2 ) in the differential equations which is associated with the two massive particle threshold.
Following the work started in [18] we organize all Feynman integrals required for the computation of qq → V V into three different integral families named Topo A, Topo B and Topo C, where the first two topologies are needed to represent respectively the doubleboxes with adjacent and non-adjacent massive legs, while the third contains all non-planar integrals. We choose the propagators of the three topologies as listed in Table 1. Table 1. Propagators in the three different integral families used to represent all two-loop 4-point integrals with two massless and two massive legs with the same mass.
As it is well known, using IBPs, LIs and symmetry relations all Feynman integrals described by these three integral families can be reduced to a small subset, the master integrals. We performed this reduction for all integrals relevant for our process using the automated codes Reduze 1 and Reduze 2 [26,27,43,44]. After the reduction we find that all integrals can be expressed in terms of 75 MIs, some of which are actually not genuinely independent, but can instead be related to each other through a permutation of the external legs p 1 ↔ p 2 .
In [18], we described the computation of the MIs embedded in Topo A and Topo B. In the present work we conclude the computation of all non-planar MIs in Topo C.

Building up a canonical basis
It has been recently suggested [31] that a suitable choice of basis for the MIs renders their computation in the differential equation method more transparent. In particular it has been conjectured that, if all MIs for a given topology can be integrated in terms of Chen iterated integrals only [45,46] (the commonly used generalized harmonic polylogarithms, GHPLs, are a special case of these [47][48][49][50][51][52]), then there must exist a basis choice, with m = { m j }, such that the differential equations with respect to the external invariants can be cast in the canonical form: where x j are the external invariants, the differential d acts on all external invariants, and the dependence on the dimensional regularization parameter is completely factorized from the kinematics. Moreover, if the matrix A(x j ) can be written as: where A k are constant matrices and r k are simple rational functions of the external invariants x j then, by their very definitions, upon integration of (3.1), the result will only contain GHPLs of alphabet { r k } n k=1 . While casting the differential equations in the canonical form (3.1) is not strictly necessary for their integration, it is still very desirable for different reasons. In particular the MIs computed in the canonical basis m, once expanded as Laurent series in , end up having a particularly compact representation in terms of pure functions of uniform transcendentality, i.e. order by order in each MI is given only by a combination of transcendental functions of uniform weight [31]. Having a result in this form, in particular in the case of multi-scale and/or multi-loop problems, helps to handle the largeness of the intermediate expressions and also the complexity of the final result.
The issue of the existence of such a basis for any multi-loop problem remains, in particular for those cases which cannot be expressed in terms of GHPLs or general Chen Iterated integrals [53][54][55][56]. Moreover, even in those cases where it is known that the final result will contain only GHPLs, no algorithm for finding such basis is known, while only some general criteria have been pointed out recently [32,33,57,58].
In what follows we will describe in detail the procedure that we exploited in order to build up a canonical basis in the explicit case under study. While we claim no generality in this approach and no proof can be given that such approach would work in more involved cases, we found it particularly elementary and algorithmically straightforward to implement, so that its extension to more difficult cases should not present particular conceptual difficulties.

Building up the basis bottom-up in t
Before describing in detail the method we used to find our canonical basis, let us recall some notation and definitions which will be useful in the following. We start off by considering a topology (or sector) given by a set of t different propagators (matching the loop integrals of some Feynman diagram). Its sub-topologies (or sub-sectors) are defined as the set of all possible arrangements of propagators obtained from the original topology by removing one or more propagators in all possible ways. In the case of two-loop corrections to vector boson pair production in massless QCD, where all tadpoles are identically zero, the first non-zero sectors will be those with t = 3 (corresponding to the sunrise topologies), while the highest sectors will contain at most t = 7 different propagators.
As it is well known, by generating and solving all IBPs, LIs and symmetry relations for a given integral family, some sectors will be reduced to one or more MIs while some other, the so-called reducible sectors, will be completely reduced to their sub-topologies. In what follows we can completely neglect these reducible sectors.
Let us consider now a sector with a given value of t and which is reduced to n ≥ 1 MIs. As it is well known, the differential equations for the latter will in general contain all their sub-topologies as inhomogeneous terms. Therefore, it is natural to try and follow a bottomup approach in t, such that, when studying the differential equations for the MIs of a given sector with a given value of t, we can assume that all its sub-topologies fulfil differential equations already in canonical form. For the MIs of the sector under consideration we use the notation: f a (x k ; ) with a = 1, ..., n. The differential equations for the n MIs read in total generality: al (x k ; ) are at most rational functions of x k and , while the m l (x k ; ) are the sub-topologies, whose differential equations are, by construction, already in the canonical form: Let us consider now the n × n matrix of the coefficients of the homogeneous equation ab (x k ; ) } ab . Our method relies on the assumption that we can find a starting basis of MIs f a (x k ; ) such that: Obviously in the case where n = 1 the matrix C(x k ; ) reduces to a scalar and the condition 2. is always trivially satisfied. On the other hand there is no real restriction on the dependence on of the functions D al (x k ; ). As an exemplification we can assume a typical situation where they contain terms of the following form: where the functions α al depend only on the external invariants x k . Note that if the factor 1/(1 − 2 ) were substituted by any other linear factor 1/(u + v ), with u, v ∈ Z, the argument would proceed in the exact same way. Moreover, as it will become clear in what follows, a more complicated dependence on in the inhomogeneous terms (for example polynomial in ) can be, at least in principle, treated with a suitable extension of the method described below.
For every sector at a given t we proceeded as follows: 1. Starting from (3.3), and using the assumption (3.5), we first attempt to solve the homogeneous system for = 0 in terms of rational functions only. While there is obviously a priori no guarantee that this can be done in general (without introducing, for example, logarithms of the external invariants x k ), we found that, in all cases we worked with, this was always the case. If this is possible, then it is equivalent to finding a rotation f a (x k ; ) → g a (x k ; ) such that the system (3.3) becomes: where the functions C 2. Once the differential equations are in form (3.8), only the sub-topologies need to be fixed in order to achieve a complete canonical form. Assuming an -dependence as in eq.(3.9), we start removing first all subtopologies proportional to the coefficients γ al . This can be attempted performing a shift in the MIs basis as follows: where the Γ al (x k ) are rational functions of the external invariants and whose explicit form will be determined in the following. Note that, since the differential equations for the sub-topologies are already in canonical form (3.4), this ensures that upon performing this shift and partial-fractioning in we will only produce terms proportional 0 , , or 1/(1 − 2 ). Upon performing the shifts in (3.10), in fact, we are left with: The explicit form of the functions Γ al (x k ) can be, at least in principle, determined by imposing that all terms proportional to 1/(1 − 2 ) are cancelled, in other words that: Eq. (3.12) is a linear system of first-order coupled differential equations for the unknown Γ al (x k ) whose solution can be, at least in principle, as difficult as the solution of the original system (3.3). Nevertheless, in all cases that we encountered, the system could be easily solved with an Ansatz. In particular, assuming that a basis which realizes the canonical form (3.1) exists and assuming that such basis can be reached through a rotation which only involves rational functions 2 , we can write the most general Ansatz for the functions Γ al (x k ) as linear combination of all possible linearly independent rational functions 3 which appear in the original differential equations (3.8).
Collecting for the independent rational functions, and requiring their coefficients to be zero, we are left with a large system of linear equations with numerical coefficients whose solution is now, at least in principle, completely straightforward.
Note that, in a typical case, there will be more equations than unknowns and the system will be over-constrained with many equations being linearly dependent from each other. For the same reason, it is in no way guaranteed that a solution to such a system exists. Nevertheless, once more, for all cases that we worked with, a solution could always be found.
2 Note that this requirement is perfectly sensible as long as we assume that such basis can be reached from any other basis through IBPs, LIs and symmetry relations only and potential roots coming from the solution of the homogeneous equations can be rationalized, similar like in our case where the Landau variable x absorbs the root s(s − 4m 2 ). 3 Where here "linearly independent" has to be intended in the sense of a complete partial fractioning in all external invariants x k .
3. Once this step has been performed, we are left with a new system of equations for the new MIs h a (x k ; ) which reads: al (x k ) are again simple rational functions of the external invariants x k . We can now proceed removing the remaining terms which are not proportional to . This can be achieved in the same way as before by performing the shift: where again the G al (x k ) are rational functions of the external invariants.
Note that, since the MIs depend in general on many external invariants x j , for every a, l fixed, there has to exist a single function G al (x k ), such that the terms not proportional to in (3.13) cancel under the shift (3.14) for all differential equations in all external invariants x j . This condition can be rephrased as: With the same assumptions as before we can solve these equations with an Ansatz imposing that the solution must be a linear combination of rational functions in the external invariants only. Again, in all cases where we applied this method, a solution could always be found.
4. After the final shift (3.14) is performed the canonical form is reached: where the H We would like to emphasize that there is no guarantee that, given any sets of MIs to start with, all steps described above can be always successfully carried out. These require in each instance to find a shift which only involves rational functions and which eliminates at every step the un-wanted terms in the differential equations. It should also be noted that the two steps 2. and 3. must be performed in this order. It is clear in fact from (3.11) that the first step will produce in general terms in the equations which are proportional to 0 , and which will be removed only during the following step.

Extension to polynomial dependence on
In the very same way as discussed above we can treat the more general case where the differential equations (3.8) admit a polynomial dependence on in the sub-topologies and/or higher powers of factors 1/(u + v ). Let us consider for simplicity a sector with only one master integral f (x; ), which depends on one single external invariant x and which has only one sub-topology m(x; ). Let us assume again that the differential equation for the sub-topology is in canonical form Moreover, let us suppose that the differential equation for the MI f (x; ) is almost in canonical form except for a term proportional to n , where the C j (x) are all rational functions of the external invariant x. Using again (3.17), we can perform the shift: Inserting this expression in (3.18) and using the fact that the differential equations for the sub-topology m(x; ) are in canonical form, it is clear that we will produce terms proportional to n and n−1 only. We can then fix the function G(x) imposing that the shift (3.19) removes all terms proportional to n , being in this way left only with terms proportional to n−1 . Proceeding in this way, starting from the highest power of , we can tentatively remove all undesirable terms from the differential equations and bring them to the canonical form (3.1). The very same idea applies to higher powers of factors 1/(u + v ) which multiply any subtopology whose equations are already in canonical form. Starting from the highest powers we can tentatively remove all terms one after the other until we are reduced to the case treated in the section above.

The basis
Applying the algorithm described above we could find a canonical basis for all MIs contributing to both planar and non-planar corrections to qq → V V . In total there are 75 MIs, some of which are not truly independent but instead are related by a permutation of the external legs p 1 ↔ p 2 (or equivalently q 1 ↔ q 2 ).
Following the ideas described above, we started building up our basis with the following initial choice of MIs: Thick lines denote the massive external particles, dots additional powers of the corresponding propagator. We introduced the short-hand notations p ij = p i + p j and p ijk = p i + p j + p k , where p 3 = −q 1 . In some cases we denote additional numerators of the integrand, where the definition of the loop momenta k, l is implicitly given by the corresponding integral family in table 1. The familiesĀ andB emerge from A and B by swapping the two incoming momenta p 1 and p 2 .
Given any of the integral families T = A, B, C (and the crossed variants), every integral is defined with the integration measure: where the indices T j run over the propagators in the different integral families, d = 4 − 2 and The transition between (m 2 ) > 0 and (m 2 ) < 0 is understood to be taken with (m 2 ) > 0.
Let us note here that, as already discussed above, the main point of our bottom-up construction of the canonical basis is that we do not need to look at the global properties of the 75 × 75 matrix, but we can move step by step for increasing values of t, treating separately the differential equations for different sectors which are topologically disentangled from each other. In this sense we need to verify separately for every topology that the two requirements of (1) triangularity as → 0, and (2) linearity in , are satisfied only for the homogeneous part of the sub-system 4 .
Deriving the differential equations for the basis above one immediately sees that: 1. All sub-systems of differential equations, topology by topology, have the property that their homogeneous part is triangular (or even decouples) in the limit → 0.
2. On the other hand almost all sub-systems fulfil the property of linearity in , except six, which are: In particular the homogeneous systems contain terms proportional to 2 and 1/(1−2 ).
In order to put these sub-systems in the right form (note that four of them are equal two-by-two under a permutation of the external momenta) one can proceed in different ways. Since the systems are very simple (in all cases 2 × 2 systems) one can study directly the homogeneous parts of the latter and, one-by-one, look for appropriate linear combinations of the two MIs which remove the terms in 2 and 1/(1 − 2 ). In particular in the case of sectors 181 and 182 (and their crossing) it is easy to see that these unwanted terms are removed by a simple rescaling of the two masters.
For the first two sectors (A166 and A198), this is not enough and a proper linear combination of the two MIs in the sector is needed. Nevertheless, having the exact solution for all planar MIs at hand [18] it is clear that, for both sectors, the second MI (the one with a dotted propagator f A166 19 , f A198 21 ) is already in the right form (i.e. it is a function of uniform transcendentality multiplied by a single rational factor), and therefore must not be changed. This implies that only the first MI must be substituted by a linear combination of the two, which can be easily found imposing that also the latter becomes a function of uniform transcendentality multiplied by a single rational factor (or, equivalently, that the homogeneous part of the system contains only a linear dependence on ).
Once the differential equations for these four sectors have been put in the right form we can easily apply our algorithm to all remaining MIs. As a result we get a new basis m = {m j } with j = 1, ..., 75 whose differential equations with respect to both external invariants (x, z) are in canonical form. We enclose the explicit definition of the basis expressed in terms of the initial choice of MIs depicted above in Appendix A.

Comments on the basis change
Following the construction of the canonical basis as described in the previous section, it appears clear how, at least in the case under consideration, the issue of finding a canonical basis for a set of master integrals can be identified to that of being able to integrate out the homongeneous part of the system in = 0 in terms of rational functions only 5 . In this sense, having a basis of MIs whose differential equations are in canonical form is, for practical purposes, almost equivalent to having a basis whose differential equations are triangular for = 0, and whose homogeneous parts can be integrated in terms of rational functions only. From this point of view, casting the system of differential equations into canonical form consists in separating into two different steps the two conceptually different issues of: (a) integrating the homogeneous parts of the equations (which provide the somehow "trivial" rational prefactors of the MIs), and (b) integrating the "non-trivial" dependence of the master integrals on transcendental functions (and in particular on GHPLs). It looks reasonable to think that such a "factorization" can be achievable as long as the master integrals are expressed in terms of GHPLs only (and Chen Iterated integrals in general). 5 The same seems anyway to apply also to more involved cases where a preliminary analysis has already been carried out. Moreover note that similar conclusions on the behaviour of the differential equations for = 0 are drawn in [33].
On the other hand, it is not clear how and if this structure could be preserved or generalized in the case of more complicated functional behaviours.
It is also important to stress that, both if the equations are in canonical form, and if the equations are triangular for = 0 (with the homogeneous part solvable in terms of rational functions only), the integration of the differential equations becomes straightforward. A task which remains is the determination of the boundary terms, which is nontrivial in particular in the case of non-planar integrals. In this regard, having the equations in canonical form does not solve any conceptual difficulties by itself. In either case, the canonical basis gives not only a clearer view on the structure of the problem, it also has great practical advantages due to the simpler and more compact expressions which need to be handled in this approach.

Differential equations
Given the basis in Appendix A we can derive differential equations in both independent variables (x, z). As already anticipated the equations take the canonical form (3.1): where the differential d acts on the two variables x, z. The matrix A(x, z) does not depend on and it can be decomposed as where the A k are constant matrices, whose entries are in particular just rational numbers, and the r k are polynomial functions of (x, z) and constitute the so-called alphabet of the GHPLs which will be needed in order to integrate the equations: Expanding in , the canonical form ensures full decoupling of the differential equations (3.22) order by order in . The integration up to a boundary term becomes trivial and can be carried out entirely algebraically.

Integration and boundary conditions
We consider the full system of differential equations for all 75 master integrals in a uniform manner. Our normalization is such that the solutions for our master integrals have a Taylor expansion, where the weight 0 contributions start at 0 . We solve the full vector of coefficient functions m (i) order by order in up to and including weight 4.
For master integrals depending on z we choose to integrate the partial differential equation in z for fixed x implied by (3.22). This gives us the solution up to a function of x, which needs to be fixed by additional constraints discussed below. For master integrals independent on z we integrate the partial differential equation in x, which determines the solution up to a constant. It is obvious that this procedure naturally leads to iterated integrals. The d ln form of the differential equations ensures that the iterated integrals can be expressed in terms of Goncharov's multiple polylogarithms Here, the w i are complex rational functions of the indeterminants. To handle non-linear letters we also employ generalized weights [59].
where f (o) is an irreducible rational polynomial and o is a dummy variable.
In order to fix the boundary terms we use two ingredients. For some of the simplest integrals, namely a small number of tadpole, bubble and triangle integrals, we use their known analytic solutions from the literature [18,60]. For all other integrals we require the absence of logarithmic divergencies for the solutions in certain kinematical limits. This requires the linear combinations of master integrals multiplying the corresponding d ln terms in the differential equations to vanish in the respective limit. This completely fixes the remaining boundary terms, i.e. the unknown functions of x respectively the constants.
Since we consider also non-planar integrals it is unavoidable to deal with cuts in s, t and u at the same time, i.e. we have to handle uncrossed and crossed kinematics. From the Feynman parameter representation it is clear that there is a Euclidean region with s, t, u, m 2 all less than zero, such that the integrals are real. From the on-shell relation (2.4) one sees, however, that it is not possible to parametrize this region employing real valued parameters x, z and m 2 . Note that in the scheme [34] employed here, the solutions develop explicit and implicit imaginary parts already during the iterative integration procedure.
We require regularity of each integral in some of the following collinear and, depending on its cut structure, threshold limits: anymore, we anticipate its former presence with (m 2 ) < 0. This translates to (small) imaginary parts (x) > 0, (z) < 0 and (y) = ((1 + x 2 − xz)/x) > 0. The limits were computed using in-house Mathematica packages for multiple polylogarithms [61,62], where we employed both, coproduct based and non-coproduct based limit algorithms. In practice we employ small but finite imaginary parts, such that the complex parameters fulfil the on-shell relation (2.4). We match the constants appearing in the limit computations by a numerical fitting procedure. This step utilizes the numerical evaluation of multiple polylogarithms [63] in GiNaC [43].

Solutions and checks
We obtain the solutions in terms of GHPLs of argument z and weights {0, −1, x, 1/x, (1 + We performed several checks on the results. First of all, we integrated the whole 75×75 system of differential equations at once, fixing consistently all boundary conditions using the limits described above. We explicitly verified that our solutions fulfil the partial differential equations both in x and in z. This is a non-trivial check for integrals depending on both variables, for which we fixed x-dependent boundary terms by regularity conditions.
As a subset of the integrals considered here, we re-calculated all non-trivial planar master integrals presented in [18]. Taking z → z + i0 and swapping masters 3 and 4 of sector B213 in the result of that reference, we translate these expressions to our new functional basis and find perfect agreement at the analytical level.
For the previously unknown non-planar master integrals we compared our results against numerical samples obtained with the sector decomposition program SecDec2 [64,65]. We found the program particularly useful since it allowed us to perform checks of our results both in the Euclidean and in the physical region. In the Euclidean region we set x to a truly complex number. Note that to obtain a real number at all consists already in a very non-trivial check of our solution. In particular, we could verify all our master integrals in the Euclidean region with a typical precision of at least 4 digits for the weight 4 coefficients. On the other hand, the numerical evaluation by sector decomposition in the Minkowski region was much more cumbersome. Using SecDec2 we could evaluate all corner integrals (integrals with no dots nor scalar products) up to weight 4, finding good agreement with our result. For sectors involving more than one master integral we additionally considered integrals with dots and/or scalar products. For them we could check at least the weight 3 contributions, in some cases also the weight 4 parts. The combination of these checks in the Euclidean and in the Minkowski region provides stringent evidence for the correctness of our results.
-16 -For the purpose of numerical evaluation in the physical region the primary form of our solutions is not optimal yet, e.g. because the multiple polylogarithms are not single valued and their numerical evaluation is not straightforward. We follow the procedure described in [37] and project onto a new functional basis which consists of Li 2,2 , classical polylogarithms Li n (n = 2, 3, 4) and logarithms. The Li-functions are related to the G-functions via Li n (x 1 ) = −G(0, · · · , 0, 1 In our new functional basis we allow for rather complicated rational functions of x and z.
We choose them such that the functions are real valued and the imaginary parts of the solutions are explicit over the entire physical domain.
In [40] it was demonstrated that this method works also in the presence of generalized weights, which could in fact be eliminated at the level of the amplitude. In the present case we work at the level of the master integrals. Also here, we successfully apply this projection onto real valued functions and eliminate all generalized weights {[1 + o 2 ], [1 + o + o 2 ]} using a coproduct based algorithm.
We can actually go one step further and restrict the target function space even more. For the functions Li n (x 1 ), Li 2,2 (x 1 , x 2 ) we select real arguments with In this way, the multiple polylogarithms are not only real valued but correspond directly to a convergent power series expansion see e.g. eq. (20) of [63]. While it is not a priori obvious that such a restricted set of functions is sufficient to represent our master integrals, we find that this is indeed the case. Our choice of functions drastically improves the numerical evaluation time, since it avoids additional transformations which would be required otherwise to map to an appropriate expansion. Evaluating all master integrals discussed in this paper takes only fractions of a second in a generic phase space point on a single core. For completeness, we also expand our solutions both at the production threshold and in the small mass region. The threshold region is characterized by β → 0 for fixed cos θ, where β = 1 − 4m 2 /s is the velocity of each vector boson and θ the scattering angle in the center-of-mass frame, such that z = 1 + 2β(β + cos θ)/(1 − β 2 ). We find it convenient to directly expand our full solutions in the real-valued function representation, rather than the individual G-functions. The expansion contains β, ln β, GHPLs of argument cos θ and weights {−1, 1} as well as the constants ln(2) and Li 4 (1/2). Similarly, we consider the small mass limit m 2 /s → 0 for fixed φ = −(t − m 2 )/s. The expansion contains m 2 /s, ln(m 2 /s) and GHPLs of argument φ and weights {0, 1}. The first couple of orders for both expansions as well as our results in terms of real-valued functions are provided via ancillary files on the arXiv.

Conclusions
In this paper, we computed the full set of master integrals relevant to the two-loop QCD corrections to the production of two vector bosons of equal mass in the collision of massless partons. These two-loop four-point functions are computed using the differential equation method [24,[28][29][30]. We describe in detail how we find a canonical basis [31] for the master integrals. In this basis, the differential equations for the master integrals can be solved in an elegant and compact manner in terms of iterated integrals. These general solutions are then matched onto appropriate boundary values, requiring non-trivial transformations of the iterated integrals. Our analytical results for all master integrals are expressed in terms of multiple polylogarithms, they are provided with the arXiv submission of this article. We find that it is possible to employ a restricted set of multiple polylogarithms, which allows for a particularly fast and precise numerical evaluation. We validated our solutions against numerical samples obtained using sector decomposition.
With the full set of master integrals derived in this paper, it is now possible to derive the two-loop corrections to the amplitudes for qq → W + W − and qq → ZZ, and to compute the NNLO corrections to the pair production of massive vector bosons. Combined with precision measurements of these observables at the LHC, these results will allow for a multitude of tests of the electroweak theory at unprecedented precision. m52 = 3 m 4 (1 + x + x 2 − xz) We remark here that even if the formulas look in some cases rather cumbersome, they are always at most linear combinations of the starting basis f j with rational coefficients. Obviously, choosing differently this starting basis can simplify or even complicate substantially these relations. On the other hand the main point of the derivation given in Section 3 is to show how, starting from a basis whose differential equations fulfil some initial requirements, a canonical basis (if existing!) can be built in an almost algorithmic way.