Master Integrals for the mixed EW-QCD corrections to the Drell-Yan production of a massive lepton pair

We showcase the calculation of the master integrals needed for the two loop mixed QCD-Electroweak virtual corrections to the neutral current Drell-Yan process $(q\bar{q}\rightarrow l^+ l^-)$. After establishing a basis of 51 master integrals, we cast the latter into canonical form by using the Magnus algorithm. The dependence on the lepton mass is then expanded such that potentially large logarithmic contributions are kept. After determining all boundary constants, we give the coefficients of the Taylor series around four space-time dimensions in terms of generalized polylogarithms up to weight four.


Introduction
The Drell-Yan (DY) processes, depicted at leading order in figure 1, have big cross sections and clean experimental signatures, and thus are one of the best-studied processes at the LHC. In particular they can be used to determine important parameters in the electroweak (EW) sector like the weak mixing angle and W boson mass [2,3]. The DY processes are also playing an important role as standard candles for the LHC in the form of luminosity measurements and detector calibrations. Furthermore the abundance of clean data make the DY processes a perfect place to determine parton distribution functions and search for Beyond Standard Model (BSM) physics. All these applications rely on a precise theoretical description of the DY processes, making it crucial for the physics program at the LHC.
The perturbative corrections to the Drell-Yan processes can be divided into two classes. The pure QCD corrections only occur in the initial state of the DY processes, due to the colorless nature of the leptonic final state. These corrections are known differentially up to next-to-next-to leading order (NNLO) [4,5] and inclusively at next-to-next-to-next-toleading order (NNNLO) [6]. In contrast the EW corrections can involve both the quarkonic initial and leptonic final state. These corrections have been computed at next-to-leading order (NLO) [7][8][9][10][11][12][13][14][15][16][17][18] and there is an ongoing effort to extend the computation to NNLO [19][20][21][22].
Starting at NNLO the EW and QCD corrections start to mix and are currently assumed to be the largest unknown correction in the high energy region [23]. These mixed corrections can be further divided according to the number of vector boson exchanges. While the factorizable contributions are characterized by a single vector boson exchange between the initial and final state, the non-factorizable contributions involve the exchange of two or more bosons. Among the factorizable Feynman diagrams, the mixed double virtual corrections  were computed for the W and Z boson decay [24,25] and the Z boson production [26]. The double real contribution to the total cross section for the on-shell single gauge boson production has been presented in [27] and the O(αα s ) corrections to the total partonic cross section of the process qq → Z + X is calculated in [28].
For the non-factorizable contributions, significant work has been done in the QCD-QED sector [29][30][31] and by adopting the pole approximation [32][33][34][35]. Nevertheless BSM physics might also show up in regions outside the resonance and therefore it is important to have control over the non-factorizable corrections beyond the pole approximation. Currently all ingredients, including the master integrals [36][37][38], for the construction of the amplitudes in the zero lepton mass limit are known. However there are observable effects due to the lepton masses, when considering collinear radiation of photons. While the KLN theorem ensures that for a fully inclusive observable, all these effects cancel out, the use of lepton identification cuts breaks the fully inclusive nature of the measured cross sections [15]. These analysis cuts give rise to large logarithmic contributions of the form log(s/m 2 l ) and capturing these logarithms requires the calculation of master integrals with a non zero lepton mass, which are the main subject of this publication.
Fortunately the huge number of multi scale integrals appearing in loop amplitudes are linearly dependent through integration-by-parts identities (IBPs) [39][40][41]. Therefore this huge number of integrals can be expressed in terms of a much smaller set of basis integrals called master integrals. Interestingly the IBPs can be also employed for the solutions of these master integrals, since their derivatives in respect to the kinematic invariants lay within the same space that is spanned by the IBPs. The resulting first order differential equations [42][43][44] can then be integrated, in order to obtain solutions for the sought-after master integrals.
Recently it has been noticed that the freedom of choosing different sets of master integrals can be exploited to find particular simple differential equations, the so called canonical forms [45]. Such forms are characterized by a total differential in d log form and by the factorization of the dimensional regularization parameter from the kinematics. The arguments of the d log's are called letters and together form the alphabet of our problem. While canonical forms greatly simplify the solution of differential equations and expose many of the underlying mathematical structures, finding these forms can be challenging and has l l Z q q g γ Figure 2: An example Feynman diagram for the non-factorizable mixed QCD-EW correction to the neutral current DY process.
been the subject of several publications [46][47][48][49][50][51][52]. In this work we follow the algorithm outlined in [47,53] and first find a form that is linear in the dimensional regularization parameter, which then is brought to the -factorized form through the Magnus algorithm [47]. This approach has been shown to work in a variety of applications [36,[53][54][55][56][57][58], including cases with many kinematic invariants, non-planar integrals and non-rational alphabets. The resulting -factorized form is then expanded in the ratio of lepton over Z boson mass, leading to a rational alphabet. This new alphabet agrees with the one presented in the massless calculation [36,37] up to the desired logarithms in the lepton mass. Due to the rational nature of the alphabet our results can be conveniently written in terms of generalized polylogarithms up to weight four. The results and the corresponding rotation matrix found by the Magnus algorithm are given in the ancillary files of the arXiv version of this publication.
Throughout this computation we made use of publicly available codes Kira [59] and Reduze [60] for the generation of the IBPs and the differential equations, LiteRed [61] for the dimensional reduction identities, SecDec [62] for the numerical validation of our results and GiNaC [63] for the numerical evaluation of the generalized polylogs.

Notation
This article considers the mixed QCD-EW two-loop corrections to the Drell-Yan process (2.1) specified by the following kinematics where the lepton mass m l was kept non-zero throughout the calculation. In particular our work is concerned with the quantum corrections involving the exchange of a photon and a Z boson between the initial and final state, depicted in figure 2. The appearing integrals are of the form with the normalization factor (2.4) and the propagator definitions (2.5) In the above definitions k i denote the loop momenta and the normalization factor C( ) has been chosen such that the canonical integral I 9 is set to one.

Differential equations
Integration-by-parts identities allow us to express the derivative of a complete set of Feynman integral in respect to some kinematic invariant as a linear combination of the initially chosen Feynman integrals. This leads to a coupled first order differential equation, which can be solved in order to determine those Feynman integrals. For the process under consideration it is convenient to combine the invariants into three dimensionless ratios resulting in a set of three partial differential equations The solution of these differential equations can be written as series of iterated integrals W over the matricesÃ and a vector of boundary constants F 0 By rescaling the master integrals F with the appropriate powers in the dimensional regularization parameter = 4−D 2 we can ensure that the matrix W exhibits a Taylor series in The matrices W can be obtained by integrating one of the partial differential equations and then fixing the resulting integration constant by matching the derivative of the obtained solution successively to the other partial differential equations. In our case we choose to integrate first in x, then y and finally z, resulting in the following recursive formulas where at the first step W (−1) is replaced by the identity matrix. An important step in the integration of a differential equation is the identification of a functional basis that includes all integrals encountered during those integrations. For the presented integrals this was achieved by choosing a particular basis of master integrals, where the dimensional regularization parameter factorizes from the kinematics, which are encoded in a d logform. After an expansion for small z all arguments of the d log's are simple rational functions, which enable us to express the integrals in terms of the well known generalized polylogarithms [63][64][65][66]

9)
G(a n , . . . , a 1 ; u) = u 0 dt t − a n G(a n−1 , . . . , a 1 ; u) , (3.10) with weights a i and argument u. In our case the weights are drawn from the sets for arguments x, y and z respectively. Due to the expansion for small z, some integrals can not be directly expressed as generalized polylogartihms. Nevertheless after repeatedly applying integration by parts identities of the form we either obtain an integral corresponding to a generalized polylogarithm or a purely rational function.

-factorized form
As a first step towards the calculation of the relevant master integrals we identify a special set of master integrals, for which the dimensional regularization parameter factorizes from the kinematics. This can be achieved in a two step process through the Magnus algorithm, presented in [47,53]. First we identify a special set of master integrals which are depicted in figure 3. These integrals satisfy a precanonical differential equation, which has a linear dependence on the dimensional regularization parameter . Secondly we follow the Magnus algorithm to build a rotation matrix, which identifies the corresponding factorized master integrals 12 m 2 z t F 31 + t − m 2 l 2 2 F 15 + F 16 + 6 F 31 + 4m 2 z F 32 + 2 t − m 2 l F 9 + 4 m 4 z t F 32 , where we introduced the abbreviations r 1 = √ −s 4m 2 l − s and r 2 = m 2 z − 4m 2 l . The transformations given in eq. (4.2) and (4.1) are provided in the ancillary files accompanying the arXiv version of this publication. The presented integrals satisfy a set of three partial differential equations where the dimensional regularization parameter has been factorized from the kinematic dependence encoded in A.

Solution
The vastly different size of the lepton mass compared to the Z boson mass allows us to expand the differential equations for small values of z = m 2 l m 2 z . This expansion captures the important logarithms in the lepton mass, while keeping the analytic complexity at a manageable level. In particular for the leading term in the z-expansion the matrices A x and A y are independent of z and are completely comprised of rational entries. Combining the expanded differential equations into a total differential exposes the d log form of our problem with an alphabet that up to the logarithm in the lepton mass is identical to the massless case [36,37] A precise approximation of the integrals appearing in the amplitude is only guaranteed if the integrals T i are well described within the expansion. For this reason the kinematic factors defined in eq. (4.2) need to be considered, when determining the appropriate point at which we can drop all higher order terms. Expanding the kinematic factors for small lepton masses (small z) we find divergences in those factors For this reason we have expanded our canonical integrals I j up to the first order in z such that all finite pieces are captured for the integrals T i . These higher order terms can be integrated as described in section 3 and lead to rational terms in our canonical master integrals. Since the rational functions do not contribute to the alphabet and are irrelevant for the analytic continuation, the differential equations were integrated under the same constraints as the massless case which correspond to the unphysical region Furthermore the analytic continuation to the physical region follows the same recipe that was already laid out in [36].

Boundary conditions
By the nature of determining the Feynman integrals through their derivative, a boundary constant has to be specified in order to determine the exact solution. These constants can be obtained by either matching our solutions to easier integrals through some appropriate limit or by demanding the absence of unphysical thresholds in our alphabet. For our problem at hand we used the following conditions • The integrals I 3,14 were provided as an independent input.
• In the limit m 2 l → 0, integrals I 37,43 were matched against their massless counterparts presented in [36].
• Taking the limit t → −m 2 z and demanding its regularity results in relations between boundary constants of integrals I 24,31,49 .
• Integrals I 42,48 were fixed by demanding regularity in the limit t → −m 2 z s m 2 z −s . In our normalization (2.4) the independent input integrals have the following expressions The weight structure of the boundary constants established through the input integrals persists for all master integrals. Namely the boundary constants at order i for 0 ≤ i ≤ 4 are proportional to the corresponding constants 1, 0, ζ(2), ζ(3) and ζ(4). The obtained boundary constants finalize the determination of the master integrals, whose expression are attached to the arXiv version of this publication.

Consistency checks
In order to verify the obtained solutions we performed a number of checks: • We independently computed the equivalent one-loop integrals to our process and checked that all factorizable integrals I 1,2,3,4,5 analytically agree with the actual product of one-loop integrals.
• Our expanded expressions were compared with the exact integrals obtained numerically by the package SecDec and we found agreement within the expected errors due to our approximation.
• Following the steps outlined in section 7.1 of [55], we took the lepton mass to zero for 31 1 combinations of our integrals and found agreement with the solutions presented in [36].
The zero lepton mass limit for the last check was taken by performing a Jordan decomposition of the pole matrix A z,−1 defined as The Eigenvectors v = c i I i of this matrix define linear combinations of our canonical integrals, whose behavior in the limit z going to zero is defined by the corresponding Eigenvalue where h( ) is some function that is independent of z. The 31 Eigenvectors belonging to the Eigenvalue 0 define linear combinations of our integrals that are finite in the limit z → 0. These finite combinations can be employed in two ways. Firstly we can now safely take the z → 0 limit directly on the analytic expressions of these combinations. Secondly we take the lepton mass to zero directly at the integrand level of these combinations. The resulting integrals can then be expressed in terms of the master integrals presented in the massless calculation [36]. Finally we found that the 31 expressions derived from our analytic results were in full agreement with the equivalent expressions derived by taking the lepton mass to zero at the integrand level.

Conclusions
The subject of this publication was the calculation of the previously unknown master integrals, needed for the two-loop mixed QCD-Electroweak virtual corrections to the Drell-Yan process (qq → l + l − ). The lepton mass dependence was kept up to logarithmic terms, such that the incomplete cancellation of these potentially large contributions can be studied in cross sections which are not fully inclusive due to lepton identification cuts. This study requires the knowledge of the two-loop virtual amplitudes, whose computation is now feasible through this publication. The 51 master integrals were evaluated with the method of the differential equations. In particular we found a set of MIs obeying a precanonical system of differential equations, which was brought to an -factorized form with the help of the Magnus exponential. After an expansion for small lepton masses, the boundary conditions were imposed by matching the solutions onto simpler integrals at special kinematic points, or by requiring the regularity of the solution at pseudo-thresholds. Finally the coefficients of the Taylor series around four space-time dimensions were given in terms of generalized polylogarithms up to weight four.