Two-Loop QCD-EW Master Integrals for Z Plus Jet Production at Large Transverse Momentum

The production of electroweak $Z$ bosons that decay to neutrinos and recoil against jets with large transverse momentum $p_\perp$ is an important background process to searches for dark matter at the Large Hadron Collider (LHC). To fully benefit from opportunities offered by the future high-luminosity LHC, the theoretical description of the $pp \to Z+j$ process should be extended to include mixed QCD-electroweak corrections. The goal of this paper is to initiate the computation of such corrections starting with the calculation of the Feynman integrals needed to describe two-loop QCD-electroweak contributions to $q \bar q \to Z+g$ scattering amplitudes. Making use of the hierarchy between the large transverse momenta of the recoiling jet, relevant for heavy dark matter searches, and the $Z$ boson mass $m_{Z}$, we present the relevant master integrals as a series expansion in $m_{Z}/p_\perp$.


Studies of electroweak vector bosons play an important role in experiments at the Large
Hadron Collider (LHC). Large cross sections of processes with one and even two electroweak gauge bosons and their clean leptonic decay signatures allow many precision tests of the Standard Model (SM) [1,2]. However, these same large cross sections imply that processes with vector bosons are important backgrounds to searches for physics beyond the Standard Model (BSM). For dark matter searches specifically, the most important among them is pp → Z(→ νν) + j (with j denoting a hadronic jet), where the produced neutrinos escape detection and lead to the missing energy signature [3,4]. It was pointed out in Ref. [5] that understanding this process to the degree required for dark matter searches at the LHC can not rely on purely data-driven techniques and that theoretical input is required.
Interestingly, in spite of the fact that theoretical studies of Z + j production have a long and successful history, it was argued recently [6] that one important ingredient is still missing.
Indeed, to search for heavy O(1 TeV) dark matter particles, we look for events with high-p ⊥ jets and large missing energy; for the main background process pp → Z + j, this O(1 TeV) missing energy is comparable to the transverse momentum of the vector boson p ⊥ [7]. It is well-known that for processes with large p ⊥ , electroweak corrections get enhanced by electroweak Sudakov logarithms log 2 p ⊥ /m Z . Hence, although the QCD corrections play a much more important role at small p ⊥ than the electroweak ones, this hierarchy becomes less obvious at high transverse momentum. As a consequence, a recent state-of-the-art theoretical analysis of vector boson plus jet production [6] came to the conclusion that unknown mixed QCD-electroweak corrections are among the largest sources of uncertainty in a theoretical description of Z + j production at high p ⊥ at the LHC. Developing a better understanding of these corrections is the long-term physics goal of this paper.
The theoretical understanding of Z + j production in hadron collisions is very advanced.
Indeed, the QCD corrections to Z +j production are known to an impressive next-to-next-toleading order (NNLO) [8][9][10]; at high p ⊥ , these corrections increase the next-to-leading-order (NLO) predictions by about 50%. The electroweak corrections have been computed up to NLO [11][12][13] and electroweak Sudakov logarithms, relevant at high p ⊥ are known through next-to-leading logarithmic terms at NNLO in the weak coupling constant [14][15][16].
The relevance of mixed QCD-electroweak corrections for dark matter searches was quan- tified in Ref. [6] by comparing a multiplicative and additive prescription for combining NLO QCD and NLO electroweak corrections. It was found that the difference between the two prescriptions leads to a 5 − 10% ambiguity in the total cross section at high p ⊥ [6]. Since the NLO QCD K-factors are large and the EW corrections can amount to up to a few tens of percents due to the EW Sudakov logarithms [14,17], it is crucial to accurately compute the NNLO mixed QCD-EW corrections to achieve few-percent accuracy for future highluminosity LHC (HL-LHC) measurements. These corrections involve unknown two-loop virtual Feynman integrals with massive propagators and an off-shell leg.
In Figure 1 we show some representative 2 → 2 Feynman diagrams that contribute to Z + j production at the LHC at various orders in QCD and EW. The dominant partonic contributions at LO are the qq → Zg and qg → Zq processes. The majority of the missing two-loop mixed QCD-EW diagrams are those where the LO graphs are supplemented with an extra gluon, and either a massive V = W ± , Z or massless photon exchange. There are also contributions from closed quark loops (cf. bottom middle graph in Figure 1), with either a massive (top) or massless quark-loop.
In this paper we will include contributions of all massless quarks but systematically neglect the contributions that involve the top quark, since these contributions contain another scale and it is very difficult to compute them. However, at (very) large p ⊥ , the top-quark contribution may, to some extent, be approximated by that of a massless quark and we can check to what degree the top quark is relevant by simply including it in the massless ap-proximation. While not ideal, this will provide us an estimate of how relevant the top-quark contribution is for mixed QCD-EW corrections to Z + j production at high p ⊥ .
In what follows, we compute the master integrals (MI) contributing to mixed QCD-EW corrections to Z + j, treating all quarks as massless, in the limit of large p ⊥ . Note that in the massless quark approximation, exchanges of a Higgs boson (cf. bottom-right graph in Figure 1) do not contribute. Furthermore, we take the Z and W boson masses to be equal, m W ± = m Z . A rough estimate of the error introduced by such an approximation is on the QCD-EW corrections, which corresponds to an error of order 1% or less on the full Z + j prediction. The only place where this approximation may be questioned is in the arguments of the Sudakov logarithms since their mass difference is not p ⊥ -suppressed. However, even in this case, it leads to a tiny relative modification of QCD- Since this expansion improves at higher p ⊥ , our results are certainly valid within 1% in the relevant p ⊥ ∼ 600 − 1500 GeV range where the statistical error is expected to be about 1 − 10% at the HL-LHC and the theory errors associated to the missing QCD-EW corrections are about 2 − 5% [6].
The method of differential equations (DE) has been very fruitful for computing MI [18][19][20][21][22][23][24]. We will follow Ref. [25], where an algorithmic method for expanding in m 2 top /(p H T ) 2 1 was used to compute the two-loop virtual amplitudes to the H+jet process at large Higgs transverse momentum p H T including finite top-mass m top effects (see also [26][27][28][29][30] for similar expansions). In this method, the DE satisfied by the MI are expanded directly in the small parameter m 2 Z /p 2 ⊥ . This method has consistently shown to reduce the complexities of computing MI with massive propagators [31]. The MI presented in this paper are expanded through next-to-leading-power in the small parameter m 2 Z /p 2 ⊥ ∼ 10 −2 for p ⊥ ∼ 1 TeV. The remainder of the paper is organized as follows. In Section II, we introduce the notation used in this paper and describe the integral families required for computing the mixed QCD-EW virtual Z + j amplitude (with all quarks massless). Section III summarizes the method of using differential equations to compute the MI as an expansion in m 2 Z /ŝ ij , while the calculation of the boundary constants to fix the MI is explained in Section IV.
In Section V we give some brief details on how to analytically continue and cross our MI solutions to other production channels and in Section VI we report on various numerical checks that were performed for the MI solutions. Finally, we conclude in Section VII.
Alongside this paper, we include ancillary files that contain the solutions for our MI.

II. DEFINITIONS AND TOPOLOGIES
As discussed in the introduction, we consider Feynman integrals that are needed to describe the mixed QCD-EW two-loop corrections to Z + j production, in the approximation when all quarks are massless. We also take the vector boson masses to be equal, m V := m Z = m W and consider the external vector boson to be on the mass shell.

A. Kinematics
We choose the kinematics such that the three partons have (incoming) momenta p 1 , p 2 , p 3 , and then the kinematic quantities are defined as Thus there are three scales in our problem: s, t, and m V . From these we can define an overall dimensionful scale, and two dimensionless quantities In the following we will mostly consider u-channel kinematics (corresponding to the particles with momenta p 2 and p 3 being incoming). This implies or correspondingly σ, χ, µ > 0.
In the u-channel we have and high transverse momentum corresponds therefore to χ ∼ O(1) and |µ| {1, |χ|}, justifying the expansion we use in our computations.

B. Integral families
In order to find the classes of integrals that need to be evaluated, we first look at the positioning of the external Z boson. We find that there are six options (up to crossings) for genuine seven-propagator sectors, shown in Fig. 2. These cases can be fitted into two 9propagator families, a planar and a non-planar. We choose the momenta of the propagators in these families as for the planar sectors and for the non-planar sectors.
We note that this discussion has not regarded internal masses. For a mixed QCD-EW contribution, there has to be at least one internal electroweak boson in any diagram. Either there can be one internal electroweak boson connecting two quark-lines directly, or there can be two vector bosons that couple to the external vector boson and to quark-lines through Since the internal vector boson can be a photon, it is possible for a diagram to have no internal masses. For cases with one internal massive vector boson, it can be anywhere in the diagram except on a line connecting to the external Z, and since the momenta can always be chosen in such a way that the external Z connects to the corner formed by momenta q 4 and q 5 (both in the planar and non-planar case) the single mass cases can have the mass With this in mind, our convention for the Feynman integrals is the following dimensionless where γ E is the Euler Mascheroni constant, = (4 − d)/2 is the dimensional regularization parameter, a ≡ 9 i=1 a i , and P f i is the ith propagator in a family f . Our convention for the propagators is such that where m 2 may be either 0 or m 2 V depending on the mass of the virtual particle.

C. Integral reductions and master integrals
It is well known that families of Feynman integrals can be expressed through a minimal basis of independent integrals, referred to as "master integrals". We perform the reductions to master integrals for each family, using IBP technology [32,33], with the program Kira [34].
Since Kira is able to identify identical integrals in different families, we treat our problem as one big coupled system of linear equations rather than treating the families individually.
In total we use a basis of 468 master integrals, constructed in such a way that the integrals in the 7-propagator sectors are independent under any permutation of s, t, and u.

III. MASTER INTEGRALS FROM THE DIFFERENTIAL EQUATIONS
Integration by parts identities can be used to construct differential equations for the master integrals. To solve these equations in the kinematic regime |µ| {1, |χ|}, we use a series ansatz [35]. The ansatz for an integral f n is of the form We call each individual expansion term that appears in the ansatz in Eq. (11) a branch.
The accuracy of this approximate expression is formally determined by the maximal value i max included in the i-sum. Once i max is chosen, all other limits in the sum are fixed from the differential equations discussed below. For most of the integrals, the i-sum starts from i = 0 but for some, divergences are present which make the i-sum start from a value as low For the j sum only negative values appear, with all integrals having values in the interval from −4 to 0. For the log-terms only k = 0 and k = 1 appear. Finally for a few integrals i and j also appear with certain half-integer values [25].
In summary for the most general case we have to consider the following values of i, j, k For the results provided with this paper, we choose i max = 1 for the top-sector integrals, which requires 1 us to choose at least i max ≥ 1 for the lower sectors in order to achieve this level of accuracy.
We solve for the the c-coefficients of Eq. (11) using the method of differential equations [18]. We combine the ν = 468 MIs into a vector f and write Note that we can choose the order of the master integrals in such a way that the matrix A takes a (lower) block triangular form, with the blocks consisting of integrals with the same 1 The reason for this is that mass-suppressed coefficients of the lower sectors enter the DE of the top-sector at i max = 1. We therefore compute the µ-series expansion of the lower sectors up to an i max -value required for reaching i max = 1 for the top-sector integrals.
set of propagators. Since all the entries of A are rational functions of µ, it follows from Eq. (11) that terms with different j or k decouple; hence expanding the RHS of Eq. (13) in µ allows us to extract relations between the coefficients c n,i,j,k in Eq. (11) for each block, j, and k at a time. At the end we are left with only ν = 468 unfixed coefficients.
At this point we have a valid solution for the µ differential equation Eq. (13). The next step is to insert that result into the corresponding χ differential equation for f , which leads to a set of differential equations for the undetermined c-coefficients. Introducing a new vector g for these unfixed coefficients, we get the differential equation As this is a one-scale system, it can easily be 2 brought into a canonical form [24] using automatic public tools such as Fuchsia [36]. We find along with rules for mapping between the old (g) and new (g) vectors. As the entries ofB are in a d log-form, this equation system can easily be integrated, leading to results of the The solutions g (w) (χ) are given in terms of generalized polylogarithms (GPLs) [37] of χ with weights w and entriesā, G(ā, χ). Finally, to obtain a complete solution, boundary values for the integrals are needed. Their computation is discussed in Section IV.
As mentioned in the previous Section, integrals can be separated into families with 0 or 1 internal masses, and integrals with 2 internal masses adjacent to the external massive vector-boson. For integrals in the first category, all polylogarithms can be expressed as GPLs with entries taken from 3 the list {0, 1, −1}, making them equivalent to harmonic polylogarithms [39]. We note that the entry 1 is spurious, in the sense that our integrals 2 Only for the integrals in the seven-propagator sectors, it was necessary to do this. For the lower sectors the equations were simple enough that it was possible to integrate up Eq. (14) with more traditional methods. 3 In our results as they appear in the ancillary files, we also have GPLs with the entries −2 and − 1 2 appearing. However, integral transformations of the form discussed in ref. [38] allow for the mapping of these GPLs to the HPL set, at the cost of introducing an argument different from χ. have no divergence in the point χ = 1. In the two-mass category, the entries are taken from are third roots of unity, and the roots of the cyclotomic polynomial x 2 +x+1. Polylogarithms with such entries are discussed in Refs. [40,41]. In addition 4 four special, transcendental constants appear in the result. They are The results have been checked numerically as discussed in Section VI.

A. The ancillary files
Alongside the arXiv version of this paper we include two ancillary files. The first file definitions.txt contains, in Mathemathica format, the definitions of the integral families, the master integrals, and a few other definitions that are relevant to interpret the results.
The other file results.txt contains the solutions for the master integrals. They are presented in the branch structure given by Eq. (11), and the branches are expanded such that each integral includes terms up to weight 4, as described above. Each branch contains an object ORD indicating the missing -order, and branches that are not exact in µ contain an object ALARM indicating the first missing power in the µ-expansion.

IV. BOUNDARY CONDITIONS
The differential equations in µ and χ in Eqs. (13), (14) fix the coefficients c n,i,j,k (d, χ) that enter the µ-expansion ansatz (cf. Eq. (11)) up to 468 integration "constants" that only depend on the dimensional parameter = (4 − d)/2. The differential equations relate the coefficients c n,i,j,k (d, χ) in branches with different i, n-values but same j, k to each other as explained in the previous Section. For this reason, the 468 constants that still need to be computed appear in branches of several master integrals and we may compute a suitable branch, in order to deduce the corresponding constant. We computed these constants either by matching the appropriate branches to massless m V = 0 master integrals, or using the so-called regularity conditions, or lastly by direct computation of specific branches at a convenient phase space-point s, t, u. We illustrate these three methods below.

A. Massless solutions
A subset of the constants are fixed by the requirement that the massless branch coefficient c n,0,0,0 in the ansatz in Eq. (11), is equal to the massless master integral f n (m V = 0) that is obtained when the vector-boson mass m V that appears in the propagators (cf. Eq. (10)) and in the off-shellness of the external leg, is set to zero at the integrand level. These massless double-box (both planar and non-planar) integrals have been previously computed in [24,[42][43][44][45] to weight 4 and we take the results for the coefficients c n,0,0,0 from there.

B. Regularity conditions
For the planar master integrals, a large subset of the remaining constants is fixed by requiring that any spurious poles or branch-points that arise from solving the differential equations vanish. The physical poles and branch-points of Feynman integrals can be understood from their cuts. Unlike the non-planar integrals, the planar Feynman integrals only have cuts in two of the three s, t, u Mandelstam variables, cf. Fig. 3. As it happens, the differential equations for some of the planar master integrals allow solutions with poles and branch-points that do not correspond to their physical cuts. Requiring that these unphysical singularities vanish allows us to fix many constants in the planar sectors.

C. Direct computation of integration constants
After using the massless branch and regularity conditions described above, we are still left with a few remaining constants, all of them appearing in branches of non-planar master integrals. We fix these remaining constants by computing suitable branches of some integrals where the required constants appear, either at a regular kinematic point or taking limiting values of s, t corresponding to a physical pole of the corresponding Feynman diagrams. For [1] [9] [7] [4] [5] [8] both not corresponding to any physical pole of the master integrals. These branches were either computed by applying the method of expansion by regions [35,46] to the Feynmanparametric representation of the master integrals, or by applying the Mellin-Barnes method.
We refer to Ref. [27] for a detailed example of a computation of a specific branch of a sevenpropagator integral through the use of Feynman parameters 5 and to Ref. [25] for an example of a computation of an integral with six propagators using the Mellin-Barnes method.
Unfortunately, the application of these methods did not allow us to compute ten constants needed for non-planar integrals with seven propagators. For these cases we computed the relevant branches in either the limit −t → 0 + or u → 0 + . Since both limits correspond to physical poles of the non-planar master integrals, their relevant branches typically behave as (−t) −1+l 1 or u −1+l 2 in these two limits. The branches that we computed were chosen in such a way that the remaining unknown constants in the solutions multiplied a χ-dependent factor that contained one (or both) of the two possible t, u singularities and we chose to compute these branches in the limit where that χ-dependent factor has a pole, either at −t → 0 + or u → 0 + . We give now an example of a constant that we computed by considering the −t → 0 + limit of a relevant branch of a seven-propagator master integral, using the method of Mellin-Barnes.
Let us consider the master integral I np2 110110111 , shown in Figure 4, where the measure D d kD d l is taken as in Eq. (9). It has a branch µ − , whose coefficient c 449,0,−1,0 contains a constant that multiplies a factor ∝ t −1 . We are therefore interested in extracting that branch and computing its limit −t → 0 + . The integral in Eq. (19) may be expressed through a Feynman-parametric representation with Symanzik polynomials as follows with u = m 2 V − s − t and κ := −1 − i0. 6 The above Symanzik polynomials equal We perform the integrals over the Feynman parameters by introducing Mellin-Barnes integrals, that split up terms inside the Symanzik polynomials as follows, The contour runs parallel to the imaginary axis in the complex z-plane and is chosen such that the singularities of Γ(−z) and Γ(λ + z) are to the right (left), respectively of the integration contour. We need to introduce seven Mellin-Barnes integrals to be able to perform the integration over all seven Feynman parameters and we are left with, We put −s = 1 and first extract the µ − = (m 2 V ) − branch by closing appropriate integration contours to the left, picking up residues at −2 − z 1 − z 2 − z 3 − 3 = − . Afterwards, we close contours such that we pick up residues at z 1 = −1 + a for any real-valued a. The latter residues correspond to the leading power pole 7 at −t → 0 + . Finally, we expand the result in . All of these steps can be performed with the packages collectively known as MBTools [48]. After expanding the result in , we are left with a sum of one-fold Mellin-Barnes integrals that may be performed by closing the contours either to the left or right, picking up a ladder of residues by the virtue of Cauchy's theorem. The resulting sum over residues may be then performed with the package XSummer [49]. We note that, had we not taken the limit −t → 0 + but instead just evaluated the branch at a regular point, e.g. −t = 1, the final expression after expanding in would have contained various three-fold Mellin-Barnes integrals that would have made the final calculation much more complicated.
The final result for the µ − branch of the I np2 110110111 at s = −1, in the limit χ = −t → 0 + equals to weight four The above result is then matched to the solution for the coefficient c 449,0,−1,0 in the point s = −1 and the limit χ = −t → 0 + , which determines the required integration constant.
We used the same method explained above to compute all ten constants that appear in branches of the non-planar master integral solutions with seven-propagators. Seven of these branches are of integrals with only one massive propagator and there the above method resulted in one-fold Mellin-Barnes integrals. For four of those seven branches (including the one explained above) we could close the corresponding contours and compute the final residue sums analytically with the XSummer package. The remaining three sets of one- 7 One finds by closing contours that there is no t −2+a or higher pole.
fold Mellin-Barnes integrals were performed numerically with the package MBTools [48] and we then used PSLQ [50] to match an integration constant onto a basis of transcendental constants up to weight four. Finally, the last three of the ten constants appear in branches of integrals with two massive propagators. These three constants were expressible in terms of two-fold and three-fold Mellin-Barnes integrals. We followed the steps explained in [51] to compute them, by mapping them onto parametric Euler-type integrals, at which point, they could be calculated using standard methods.

V. ANALYTIC CONTINUATION AND CROSSINGS
The top-sector master integrals that we include in our ancillary files are independent under crossings of the external momenta p 1,2,3 and may be directly used in the region where p 2,3 are incoming, i.e. u > 0 and χ = t/s > 0. However, the physical two-loop virtual mixed QCD-EW amplitudes describing the process i 2 (p 2 )i 3 (p 3 ) → i 1 (−p 1 ) + V with coloured partons i 1,2,3 , contain, after IBP reduction, also crossings of our chosen top-sector and lower-sector master integrals. 8 As explained in Section III, all our master integrals are expressed in terms of GPLs with argument χ. In order to see any possible cancellations among different crossed master integrals in the amplitudes, it is required to express them also in terms of GPLs with the same argument χ. In this section we briefly explain how these crossings can be performed in practice.
There are three physical scattering regions, defined in terms of the Mandelstam invariants The integrals in the ancillary files included with this paper are all defined in the region (4a) + . In Ref. [44] it is explained how to perform the analytic continuation from the region (4a) + to the other two regions (2a) + , (3a) + and we refer to that paper for details.
In order to compute the physical scattering amplitudes in the region (4a) + (or the other two), one will need all possible permutations of the external momenta Consider now for example the crossing σ 12 , which maps χ → −1−χ−µ and the region (4a) + to the region (3a) + . In order to compute the crossing σ 12 of our master integrals, we first map χ → −1 − χ − µ and expand the mapped integrals in small µ, resulting in GPLs with argument −1 − χ. Then we express these GPLs in terms of GPLs with argument χ using so-called supershuffling identities, as implemented for example in the Maple-based package HyperInt [47]. This supershuffling step needs to be done assuming one is in region (3a) + and results in all GPLs having the argument χ. After this step, one follows Ref. [44] and analytically continues the final expression from region (3a) + back to (4a) + to get the final result for the σ 12 -crossed master integrals as defined in the region (4a) + . Similar steps can be performed for the other permutations in Eqs. (31), (32).
We have implemented the above crossings and, in fact, used them to make various nontrivial checks of our master integral solutions gathered in the ancillary files, that we computed by solving the differential equations as explained in section III. Namely, many of the lowersector integrals that we computed in the total list of 468 master integrals are related by crossings of the external momenta p 1,2,3 . By applying the above crossing steps, we have checked that these crossing relations are indeed satisfied by our master integral solutions.

VI. VALIDATION OF MASTER INTEGRALS
We have performed multiple checks to make sure that the computed master integrals are correct. To begin with, two "basic" checks were performed for almost all master integrals.
Firstly, the consistency of our analytic solution of a given master integral with its corresponding µand χdifferential equations was checked. Secondly, upon deriving a boundary constant analytically with the methods of Section IV, it was validated against numerical integration using the Mellin-Barnes method as implemented in MBtools [48]. In addition to these checks and the ones mentioned in the previous Section, we employ two other which are discussed below.
One straightforward and robust way of verifying our results is to compute the master integrals numerically and compare them with their analytic solutions at a kinematic point.
To this end, we have used the computer programs pySecDec [52][53][54][55] and FIESTA [56], that are based on the sector-decomposition method, to compute integrals numerically. In this way, the master integrals with up to five denominators (t = 5) were systematically compared against numerical results and the agreement at the per-mille level was obtained. However, it was already difficult to achieve this precision for some master integrals with five propagators, due to their specific kinematics. Typically, problematic integrals have threshold singularities that are hard to treat with the sector decomposition method.
Even though there is an algorithm that is implemented in PySecDec, 9 we opted for a different approach that is known in the literature [58,59]   If we choose I d+2 such that it is finite in the → 0 limit, DRRs can be used to put constraints on lower-dimensional integrals. Indeed, if I d+2 is finite, then at the right hand identities. The exact cancellation of poles on the right hand side of Eq. (33) provides a nontrivial check of the computed MI. To check finite parts of the master integrals, we computed finite I d+2 using pySecDec. We note that DRRs do not resolve the problem with threshold singularities, but it was much easier to achieve per-mille numerical precision in the case of finite Feynman integrals. An example of such a comparison is given in Fig. 5.
Various programs were used to pursue this check. To find d = 6 and 8 finite integrals, the reduction program Reduze2 [62] that implements an algorithm of Ref. [63] was used. To generate the DRRs in Eq. (33), we used the Mathematica-based package LiteRed [64]. To perform reduction to our master integrals, we employed both Kira [34,65] and Reduze2.

VII. CONCLUSIONS
In this paper we computed the Feynman integrals needed for the double-virtual mixed QCD-EW NNLO contribution to Z + j production at hadron colliders (disregarding contributions from the top and the Higgs, as discussed previously). The computation was performed as an expansion around m Z = 0, making the result valid in the high p T limit.
This makes it useful for background determination for certain types of new physics, such as weakly interacting dark matter (see ref. [6] and the references therein).
Our results are presented as an expansion in µ = −m 2 Z /s of the form given by eq. (11), keeping terms up to O(µ), and with the coefficients of the expansion being expressed in terms of generalized polylogarithms. The results are added as an ancillary file as described in Section III A. We obtain numerical agreement with SecDec (with the help of finite basis methods) as described in Section VI. If higher accuracy is needed, there is no conceptual barrier to extending the µ-expansion to higher orders.
Additionally we have computed the integrals (when normalized as in eq. (9)) up to terms of transcendental weight 4. This is expected to correspond to the terms contributing to the amplitude up to finite orders in . It could, however, happen that intricate cancellations take place in such a way as to make higher -orders for some individual integrals needed. In that case there is nothing conceptual that prevents us from continuing the expansion to the required -order.
An obvious next step, is to use the integrals computed in this paper, to perform the complete calculation of the NNLO mixed QCD-EW correction to the double-virtual scattering amplitude for Z + j production. Additionally, one has to compute the real-virtual and the double-real contributions, before being able to perform the infrared subtractions and combine it all into a final result for the scattering cross-section. Moreover, we may be able to estimate or even compute the top-mass contribution that were disregarded in this paper as discussed in the introduction. Finally, we note that the integrals computed in this paper have a large overlap with those needed for W ± + j production, with only a few extra integral families missing. All of this may be the subjects of future publications.