Mixed EW-QCD two-loop amplitudes for $q\bar{q} \to \ell^+\ell^-$ and $\gamma_5$ scheme independence of multi-loop corrections

We perform a dedicated study of the $q \bar{q}$-initiated two-loop electroweak-QCD Drell-Yan scattering amplitude in dimensional regularization schemes for vanishing light quark and lepton masses. For the relative order $\alpha$ and $\alpha_s$ one-loop Standard Model corrections, details of our comparison to the original literature are given. The infrared pole terms of the mixed two-loop amplitude are governed by a known generalization of the dipole formula and we show explicitly that exactly the same two-loop polarized hard scattering functions are obtained in both the standard 't Hooft-Veltman-Breitenlohner-Maison $\gamma_5$ scheme and Kreimer's anticommuting $\gamma_5$ scheme.


Introduction
A precise theoretical understanding of the neutral-and charged-current Drell-Yan processes [1] in Standard Model perturbation theory is essential to the success of the Large Hadron Collider physics program. It is therefore unsurprising that significant resources have been allocated, over the course of a multi-generational effort, to bring the higherorder corrections to the Drell-Yan scattering processes under control. Very recently, the total cross section in pure Quantum Chromodynamics (QCD) was successfully calculated through to next-to-next-to-next-to-leading order in the strong coupling constant for both the neutral-and charged-current processes [2][3][4][5][6][7]. On the electroweak side, neutral gauge boson production and decay has been calculated through to next-to-next-to-leading order in pure Quantum Electrodynamics (QED) [8,9], but full electroweak (EW) corrections are known only at next-to-leading order [10][11][12]. Due to the size of the next-to-leading order electroweak corrections, it has long been of interest to calculate the mixed EW-QCD corrections of relative order αα s as well. Indeed, many partial results in this direction have been put forth over the last decade [13][14][15][16][17][18][19][20], but a complete calculation in the regime of large lepton pair invariant mass is not available yet.
Despite the significant attention these two-loop EW-QCD calculations have received, it has long remained a challenge to include the non-factorizable two-loop virtual corrections. For many years, the primary technical problem was the evaluation of the massive two-loop box-type master integrals which emerge from the most complicated Feynman diagrams (see [21][22][23] for various results). Recently, some of us successfully evaluated the complete set of two-loop box-type master integrals relevant to the non-factorizable corrections as linear combinations of standard multiple polylogarithms in the physical regions of phase space [24]. In this paper, we take an important step towards complete relative order αα s corrections by elucidating the treatment of γ 5 for the box-type diagrams and deriving the order α 2 α s two-loop helicity amplitudes for qq → + − .
As has been known since the earliest higher-order studies of the Drell-Yan process, it is important to have a solid understanding of the γ 5 problem of dimensional regularization [25][26][27][28] if one wishes to calculate all Feynman diagrams at a given order, including those where a W or Z connect the initial-and final-state fermion lines. As will be discussed in detail in what follows, care must be taken in practical calculations to consistently include scheme-dependent higher-order-in-terms, see also [29] for a recent review. In order to establish confidence in our results, we carried out two independent calculations, one where 't Hooft-Veltman-Breitenlohner-Maison's (HVBM's) γ 5 scheme [25][26][27]30] was employed and one where Kreimer's γ 5 scheme [31][32][33][34] was employed. For our calculations, Kreimer's γ 5 scheme was practically superior because it does not require the introduction of finite counterterms to restore the chiral symmetry of the Standard Model.
This article is organized as follows. In Sections 2.1 -2.4, a general discussion of the subtleties involved with the application of dimensional regularization in the presence of chiral couplings is given, both for HVBM's γ 5 scheme and Kreimer's γ 5 scheme. In Sections 2.5 -2.7 our renormalization scheme for the coupling constants, the wavefunctions, and the particle masses is described. In Section 2.8, we discuss the finite renormalizations we perform in HVBM's γ 5 scheme. In Sections 3.1 -3.4 we define the scattering process we consider more precisely and relate further technical details of the numerator algebra scripts we used to carry out our Feynman diagram calculations. In Section 3.5, we elucidate the infrared structure of the one-and two-loop scattering amplitudes we calculate and define hard scattering functions which, at leading order in the parameter of dimensional regularization, encode physical information about the kernel of the virtual corrections which remains after infrared subtraction. In Section 6.1, we provide a recipe to pass from our form factors to scattering amplitudes for states of definite helicity. In Sections 4.1 -4.3, supplemented by Appendices A.1 and A.2, explicit results for the one-loop Standard Model corrections to the Drell-Yan process are given, together with a detailed comparison to the original literature. In Sections 5.1 -5.3, we provide details specifically relevant to the twoloop mixed EW-QCD calculation and, in Section 6.2, visualizations of our final results for the order αα s , order α 2 , and order α 2 α s polarized hard scattering functions. Finally, in Section 7, we summarize our findings.

Universal features of dimensional regularization
In this section, our aim is to review some background and motivation for the more technical, scheme-dependent discussions of chiral couplings in dimensional regularization which follow in Sections 2.2 and 2.3. The simultaneous dimensional regularization of ultraviolet and infrared divergences [30,[35][36][37][38] has proven to be an essential tool for the calculation of higher-order perturbative scattering amplitudes in gauge theories. For a gauge theory Lagrange density defined in four spacetime dimensions, the method purports to provide a prescription for the continuation of all off-shell Feynman rules from four to d dimensions. All higher-order correlation functions and scattering amplitudes in the model under consideration are then taken to be meromorphic functions of d in a complex domain containing the point d = 4. Accordingly, the Fourier integrals which appear at higher orders in perturbation theory are continued from four dimensions to d = 4 − 2 dimensions as where is the parameter of dimensional regularization and µ is the 't Hooft scale, a unit of mass introduced to maintain a dimensionless action. The ultraviolet and infrared singularities of higher-order correlation functions and scattering amplitudes manifest themselves as poles in which ultimately cancel only in higher-order corrections to physical observables. While some of the prescribed continuations of the Lorentz and Dirac algebra such as and γ µ γ µ = 1 2 g µν {γ µ , γ ν } = g µν g µν = 4 − 2 (2.4) feel very natural, other aspects of the program are far more subtle and require careful explanation. Before proceeding to a discussion of these subtleties, we state a universal identity for the d-dimensional contraction of two Dirac matrices with n d-dimensional Dirac matrices interspersed which will be useful later on: (−1) i+n γ ν i γ ν 1 · · · γ ν i−1 γ ν i+1 · · · γ νn + (4 − 2 )(−1) n γ ν 1 · · · γ νn . (2.5) This identity is simply obtained by repeatedly applying Eq. (2.3) to move the γ µ factor to the left until it is finally possible to apply Eq. (2.4). The reader can easily check for the first few values of n that the well-known textbook results (see e.g. [39]) are all encoded in Eq. (2.5).
In both schemes we consider, Eqs. (2.12)-(2.16) force us to apply Eq. (2.10) to split indices,ḡ µν k µ i k ν j = k i · k j −k i ·k j . Split loop momenta lead to well-known technical complications: in sufficiently-complicated Feynman diagrams, Dirac traces produce d-dimensional Feynman integrals with numerator insertions ofk i ·k j . Fortunately, general approaches which involve the introduction of dimensionally-shifted Feynman integrals [41][42][43] have been worked out at one loop [44] and at two loops [45]. For the convenience of the reader, we provide details of a modern formulation of these developments in Section 2.4 below.
As reviewed in [46], many authors have in the past adopted the so-called naive γ 5 scheme of [28] (see also [47]) for one-loop electroweak calculations. The naive scheme amounts to a set of ad hoc rules for γ 5 manipulations designed to lead to correct results while avoiding the technical complications which might arise in a more mathematicallyrigorous handling of chiral couplings. Historically, this was justified at the one-loop level by checking a posteriori that the relevant Ward identities check out. Naturally, in order to calculate with confidence at higher orders in Standard Model perturbation theory, we feel that it is of considerable importance to adopt a mathematically consistent framework. It should be stressed that the naive γ 5 scheme is equivalent to Kreimer's γ 5 scheme (introduced below) for what concerns the calculation of the fermion self-energies and vertex form factors considered in this paper. As we shall see in what follows, the Feynman diagrams of box type are really the only ones in 2 → 2 scattering processes with sensitivity to the γ 5 problem.
1 It actually turns out that the algebraic form of the standard family of Dirac trace identities which generalizes Eq. (2.8) is the same in both schemes, but this requires some explanation.
As advertised above, this modification of the anticommutation relation has an immediate impact on the form of the Standard Model Feynman rules in dimensional regularization. As an example of particular relevance to this work, let us consider the axial part of the coupling of the Z boson to the light quarks and leptons in four dimensions:  [40]. HVBM's γ 5 scheme allows for a very familiar treatment of the usual families of Dirac traces. As pointed out already in [30], the trace of an even number of Dirac matrices with d dimensional indices may be conveniently evaluated using the recursive formula 3) to move the γ ν 1 factor all the way to the right. Finally, the cyclic property of the trace is required to identify the final term, −tr {γ ν 2 · · · γ νn γ ν 1 }, with the original expression on the left-hand side. The trace tr {γ ν 1 · · · γ νn γ 5 } (2.27) for even n may be evaluated in much the same way by replacing γ 5 by the covariant version of its definition, (2.31) However, we must simplify the trace of the summand of the first term on the right-hand side of Eq. (2.5) using the induction hypothesis before we can finish the proof: it is easy to see that is a valid identity under the assumption tr γ ν 1 · · · γ ν n−2 = 0. Eq. (2.32) provides the desired reduction of the right-hand side of the trace of Eq. (2.5), tr {γ ν 1 · · · γ νn γ 5 } = 0 (odd n). (2.35) HVBM's prescription for γ 5 has the very unfortunate drawback that the Standard Model's chiral symmetry is violated in all bare perturbative calculations involving chiral couplings. As was originally pointed out by Breitenlohner and Maison in [25][26][27] and reviewed in [40], it is possible to restore the chiral symmetry of the Standard Model as part of the renormalization program by introducing various finite counterterms on top of the usual divergent ones. Let us emphasize here that, in general, one expects all interactions in the Standard Model to require finite counterterms at some order in perturbation theory, not just those which explicitly involve γ 5 . Furthermore, these finite counterterms must generically be calculated beyond O 0 to consistently restore chiral symmetry at still higher orders. We discuss some elementary examples relevant to this article in Section 2.8.
For what concerns the particular one-and two-loop calculations considered later on in this work, Larin's principle [48] allows us to bypass many finite counterterm computations for fermion self-energies and vertex form factors involving chiral couplings. Larin pointed out that, for these building blocks, chiral symmetry implies the finite renormalizations do nothing but turn results obtained in HVBM's γ 5 scheme into analogous results obtained in an anticommuting γ 5 scheme. However, as we will discuss in greater detail in Section 2.8, the HVBM scheme requires the insertion of explicit finite counterterms and a careful treatment also of the higher-order-in-terms if box-type Feynman diagrams with chiral couplings contribute to the two-loop or higher-order corrections of interest (see also [49] for a related recent study). Overall, one comes away with the feeling that there must exist a pure dimensional regularization scheme better suited to the study of higher-order perturbative corrections in the Standard Model. In fact, the arduous and subtle character of perturbative calculations carried out in HVBM's γ 5 scheme is what led to the development of Kreimer's prescription for γ 5 .

Kreimer's prescription for γ 5 in dimensional regularization
In Kreimer's γ 5 scheme, chiral couplings are treated in a consistent way in d dimensions by retaining the four-dimensional definition of γ 5 , Eq. (2.6), at the expense of Dirac trace cyclicity. Kreimer argued in [31] that Dirac trace cyclicity is not a natural property to retain in d dimensions, due to the fact that the Dirac algebra necessarily becomes infinite dimensional. The standard trace of classical linear algebra, tr, is replaced in Kreimer's γ 5 scheme with one of several possible exotic trace operations which differ only in their reading point prescriptions. The reading point prescription effectively specifies with which Dirac matrix to start the non-cyclic trace. In our case, the Lorentz projectors we define in Section 3.4 provide a natural anchor to start all Dirac traces involving the external fermion lines. In order to clearly distinguish the trace operations in what follows, we will write Tr rather than tr when using a reading point prescription. While not relevant to the specific higherorder calculations discussed in this paper, 2 we would average over all possible choices of first entry for internal closed fermion loops. Other reading point prescriptions, such as one which mandates an average over all possible choices of first entry in all cases, could have been chosen as well, but we find our choice particularly convenient.
Kreimer's γ 5 scheme has practical advantages over HVBM's γ 5 scheme. Crucially, it does not force a change of algebraic form for any Standard Model Feynman rule when considering interactions in d spacetime dimensions. Employing a consistent reading point prescription then automatically preserves the chiral symmetry of the Standard Model in bare perturbative calculations and the complete absence of finite counterterms makes renormalization straightforward. In addition, although box-type diagrams do ultimately force the introduction of split loop momenta also in Kreimer's γ 5 scheme, the bulk of the numerator algebra can be carried out before the split enters. This stands in stark contrast to HVBM's γ 5 scheme, where split loop momenta should be introduced immediately in order to progress with the necessary algebraic simplifications efficiently (this assertion will be clarified in Section 3.3).
As alluded to above, the standard families of odd Dirac traces still vanish identically when tr is replaced by Tr, Tr {γ ν 1 · · · γ νn } = 0 (2.36) and Tr {γ ν 1 · · · γ νn γ 5 } = 0 . (2.37) For the standard families of even Dirac traces, [32] defines where σ is the set of permutations of the 2m indices which appear in the sums on the right-hand sides, subject to the restrictions (2.28), one immediately arrives at the same conclusion for the Kreimer trace of a string of Dirac matrices with a γ 5 factor appended (see also [34] for a supplementary discussion). Despite the fact that Kreimer's γ 5 scheme was introduced a long time ago [31][32][33][34], it has not been universally accepted and embraced as a superior alternative to HVBM's γ 5 scheme; researchers continue to express doubts about the applicability of Kreimer's γ 5 scheme beyond one loop (see e.g. [46] for a current review of the status of the γ 5 problem). In this article, we present an explicit application of Kreimer's γ 5 scheme to a non-trivial calculation of two-loop box-type diagrams involving chiral couplings and show that it does yield physically sensible results. Also worth mentioning is concurrent work by one of us and others, where the two-loop box contributions to gg → ZZ featuring a closed top-quark loop have been calculated in Kreimer's γ 5 scheme [50]. The results were compared against calculations in a more commonly used γ 5 scheme and found to agree after ultraviolet and infrared subtraction. While the exact top mass corrections to gg → ZZ in next-to-leading order QCD couple the ZZ pair to the massive top quark line and involve anomalous double triangle contributions, the γ 5 dependence of the mixed two-loop EW-QCD corrections to the Drell-Yan process at the level of Feynman diagrams is significantly more complicated.

d-dimensional Feynman integrals with four-dimensional scalar numerators
In this section, we describe how we deal with particularly awkward terms which often occur in the numerators of d-dimensional Feynman integrals in Standard Model perturbation theory: four-dimensional contractions of d-dimensional loop momenta,ḡ µν k µ i k ν j . As mentioned above, we decompose our loop momenta into four-dimensional and (−2 )-dimensional components, k i =k i +k i . This decomposition maps the awkward four-dimensional numerators to ordinary inverse propagators and "µ-term" insertions of the formk i ·k j , as introduced long ago by Bern and Morgan [44]. Below, we review algorithms which allow for a direct treatment of µ-term insertions at one and two loops with dimensionally-shifted Feynman integrals. While the two-loop algorithm we review was, to our knowledge, the first to appear in the literature [45], other approaches such as that of [51] are certainly possible. It seems that some computational effort is simply unavoidable, as both approaches require the calculation of auxiliary integration by parts reductions [52][53][54]. 3 We found it very appealing to work with dimensionally-shifted Feynman integrals, as all of the necessary auxiliary integration by parts identities were already calculated by us some time ago for numerical cross-checks on our master integrals [22,24], facilitated by passing to a finite basis of master integrals [56,57] through Tarasov dimension shifts [43].
As has been clear at least since the appearance of [44], non-zero one-loop Feynman integrals with µ-term insertions are of the form where the D i are the independent propagators of the integral family. Crucially, Eq. (2.41) can be rewritten in terms of standard dimensionally-shifted one-loop Feynman integrals with the same propagator structure. By absorbing the (k 1 ·k 1 ) r numerator insertion into the integration measure, it is shown in [44] that (2.42) Although Eq. (2.42) involves dimensionally-shifted integrals in an essential way, its form does not depend on the number of legs, the complexity of the kinematics, or the details of the ν i , observations which suggest a connection to the dimension shift of Tarasov [43]. Since, to our knowledge, a manipulation of the Schwinger representation furnishes the simplest possible proof of Tarasov's relation, it is natural to expect the Schwinger representation to also allow for an efficient alternative proof of Eq. (2.42). To proceed, it is useful to recall that the Tarasov dimension shift is derived in the Schwinger parametrization by simply multiplying and dividing by the first Symanzik polynomial. Furthermore, for an integral topology with t propagator denominators, the first Symanzik polynomial at one loop has the especially simple, universal form due to the uniqueness and triviality of the one-loop vacuum (tadpole) integral topology. The general Schwinger representation is most simply expressed if one first orders the Feynman propagators such that the → 0 limits of their exponents form a monotonically decreasing sequence. We assume that, in this limit, the first n + propagators appear raised 3 Note that one can circumvent the discussion of this section altogether by performing a classical Passarino-Veltman tensor reduction [55] on all one-and two-loop integrals. In fact, we found it very productive to carry out such an analysis as a comprehensive cross-check of our µ-term implementation. to positive powers and the remaining n − propagators appear raised to negative powers (we may assume without loss of generality that the final n − propagator exponents are in fact non-positive integers). In terms of the Symanzik polynomials of the integral family, U and F, we have where d 0 is an even positive integer (d 0 = 4 is of course the main case of interest to us in this article). In the absence of inverse propagators, the usual result given in terms of the Symanzik polynomials of the integral topology under consideration, U and F, is recovered. Now, as emphasized in [44], generic d-dimensional Feynman propagators at one loop have the formk 1 ·k 1 + 2P ext ·k 1 + P ext · P ext − m 2 +k 1 ·k 1 , where P ext denotes some sum of four-dimensional external momenta and m some particle mass. This form motivates a split of the loop momentum integral in the derivation of the ordinary Schwinger parametrization into four-dimensional and (−2 )-dimensional pieces. Due to the fact that, by construction, µ-term insertions are expected to lead to Feynman integrals explicitly suppressed by powers of , the authors of [45] deduced that the (−2 )-dimensional part of the momentum integral must be proportional to the -dependent part of Eq. (2.44), and therefore serve as a generating function for all µ-term insertions via differentiation with respect to U n . Indeed, by differentiating the right-hand side of the above relation r times with respect to U (1) n , we find and, by recognizing that the shift → − r corresponds precisely to a shift d 0 → d 0 + 2r in Eq. (2.44), immediately recover the form of Eq. (2.42) as expected. Fortunately, the approach described above may be fruitfully applied to the study of generic µ integrals at the two-loop level as well. This is primarily because there is a unique non-factorizable two-loop tadpole topology (see Figure 1). The main new feature at the two-loop level is that shifts in the propagator exponents appear alongside the various dimension shifts. At two loops, generic non-zero µ-term insertions can be constructed by taking products ofk 1 ·k 1 ,k 2 ·k 2 , and (k 1 −k 2 ) 2 . Due to the uniqueness of the non-factorizable two-loop tadpole topology, the first Symanzik polynomial of an arbitrary non-factorizable two-loop integral family is of the form (2.47) where T k 1 ·k 2 is the sum of the Schwinger parameters associated to propagators which depend on the dot product k 1 · k 2 , T k 2 1 is the sum of the Schwinger parameters associated to propagators which depend on k 2 1 but not on k 1 · k 2 , and T k 2 2 is the sum of the Schwinger parameters associated to propagators which depend on k 2 2 but not on k 1 · k 2 . It follows that the two-loop analog of Eq. (2.45) ought to be , and T k 1 ·k 2 . As at one-loop, the above relation suggests that arbitrary two-loop µ integrals may be rewritten as a linear combination of ordinary Feynman integrals by simply comparing the result of each term of the prescribed differentiation of the right-hand side of (2.48) with the form of Eq. (2.44).
Starting at two loops, a subtle complication can arise in our formulation due to the fact that we prefer to write all numerator polynomials as linear combinations of power products of inverse propagators and particle masses. If a propagator depending on e.g. k 1 alone appears in the numerator of a Feynman integral alongside an explicit insertion ofk 1 ·k 1 , it will not work to use the above line of reasoning to eliminate the explicitk 1 ·k 1 factor. This is simply because the k 2 1 =k 1 ·k 1 +k 1 ·k 1 term of the propagator itself depends onk 1 ·k 1 . In practice, we found it most convenient to temporarily delay rewriting all scalar products of the form k i · p j in terms of inverse propagators. Due to the fact that the k i · p j scalar products live in four dimensions, (2.48) may be applied to eliminate all explicit µ terms before rewriting the remaining scalar products as linear combinations of propagators.
To clarify the procedure described above, we consider a non-trivial two-loop µ integral from integral family A of [22,24] which may be directly rewritten using the technology of [45]: where we have introduced the notation D 7 ≡ (k 2 − p 1 − p 2 ) 2 from [22,24] on the left-hand side for later convenience. Note that, in this example, there is no clash between the explicit factor of (k 1 ·k 1 ) 2 in Eq. (2.49) and the implicit factor ofk 2 ·k 2 present in The Schwinger parameter subsets of interest are and differentiating the right-hand side of Eq. (2.48) twice with respect to T k 2 1 gives ( − 1) times an insertion of into Eq. (2.44) specialized to our two-loop example. The factor in the denominator induces the dimension shift 4 − 2 → 8 − 2 . Powers of Schwinger parameters associated to the propagator denominators induce positive shifts in the corresponding propagator exponents, and, for h powers of parameter i, an overall factor of (−1) h Γ(ν i + h)/Γ(ν i ). Powers of the Schwinger parameter associated to the numerator also induce positive shifts in the exponent, but only an overall factor of h! from the action of the derivatives in Eq. (2.44). Due to the fact that α 7 is merely an auxiliary parameter set to zero at the end of the calculation, the insertion of α 2 7 gives a vanishing contribution for the specific example considered. We find: (2.53) The integrals on the right-hand side of Eq. (2.53) can be straightforwardly evaluated in Feynman parameters to all orders in due to the fact that they are iteratively one loop. After some simplifications, we obtain To check Eq. (2.54) and our treatment of the two-loop µ terms, we invite the reader to first integrate out k 2 in integral (2.49) above and then treat the µ-term insertion in the context of the remaining integration over k 1 (i.e. by applying Eq. (2.42)). Finally, let us stress again that a numerator insertion of D 7 (k 2 ·k 2 ) 2 would not be directly covered by our method.
To be clear, factorizable two-loop µ integrals also appear in our calculations. They could be evaluated as products of one-loop µ integrals, but we prefer to treat them along the lines described above. In the factorizable case, (k 1 −k 2 ) 2 →k 1 ·k 1 +k 2 ·k 2 and we have for our generating function. This is nothing but a degenerate case of Eq. (2.48) above with T k 1 ·k 2 set to zero. Note that, while we believe that it should be possible to make a similar construction to handle µ-term insertions in d-dimensional Feynman integrals at higher loops, more effort would be required to identify the underlying vacuum topologies.

MS renormalization of α and α s
In this work, we renormalize both α and α s in the MS scheme. Due to our neglect of all fermion masses, this choice is particularly convenient for our calculation of the two-loop mixed EW-QCD corrections to the neutral-current Drell-Yan process. As explained in greater detail in Section 3.1, we neglect the gauge-invariant subset of Feynman diagrams with internal closed fermion loops throughout this work for simplicity. This includes in particular the two-loop Standard Model vacuum polarization diagrams featuring a single quark loop and a single gluon exchange. As a corollary, the form of the contributing two-loop Feynman diagrams (see Section 5.1) immediately implies that the order αα s corrections to the electric charge renormalization constant must vanish identically.
Throughout this work, we follow the convention of [13] for fixed-order results. Therefore, quantities in Standard Model perturbation theory are written in the form to all orders in the coupling constants. In this notation, neglecting contributions with closed fermion loops which are proportional to the number of light leptons (n l ) or light quarks (n q ) or involve the mass of the top quark (m t ), we have for the order αα s corrections to the electric charge renormalization constant.
Since the tree-level Drell-Yan amplitude is of order α, the trivial renormalization of the strong coupling constant is all that is required for our order α 2 α s two-loop calculation. On the other hand, the oneloop electric charge renormalization constant does make an appearance. It has been known since the work of [58] in the standard on-shell scheme where the charge renormalization constant is fixed order-by-order in the Thomsen limit (see e.g. [59] for a review). In the MS scheme it is given by, the pole part of Eq. (5.40) of [58] once one discards the first, fermionic term inside the bracket. For notational convenience, we shall suppress the qualifiers provided on the lefthand sides of Eqs. (2.57) and (2.59) throughout the rest of this work. When necessary, appropriate caveats will be provided in the text.

On-shell electroweak gauge boson wavefunction and mass renormalization
In this section, we review the on-shell renormalization of the electroweak gauge bosons in a generalized 't Hooft-Feynman gauge. Our conventions for the electroweak sector of the Standard Model Lagrange density and Feynman rules are identical to those of [58], upon which the useful resources [46,59,60] were modeled. Our electroweak gauge boson renormalization constants are also defined with exactly the same renormalization conditions, but we have chosen to factor the squared electroweak gauge boson mass, m 2 v , out of the mass renormalization constants, Z m 2 v , to render them dimensionless. As usual, a particularly convenient choice for the gauge parameter renormalization constants is because it simplifies the all-orders form of the electroweak gauge boson kinetic terms in the renormalized Lagrange density. Subsequently, our generalized 't Hooft-Feynman gauge is defined by setting all of the renormalized gauge parameters to one at the outset, This is justified because, due to Eqs. (2.60), the gauge parameters in this scheme receive no radiative corrections. Apart from the exactly-transverse photon-photon self-energy, the electroweak gauge boson self-energies naturally decompose into a sum of transverse and longitudinal components: and The longitudinal components of the self-energies, Σ L V V , are unphysical and need not be calculated explicitly.
Working in d dimensions, we see that the relevant Lorentz projectors for the transverse gauge boson self-energies can be cleanly summarized as For all of the Lorentz projection operators considered in this work, a sum over the available open Lorentz and spin indices is always implied for their action on Feynman diagrams. We find the following results for the transverse projections of the Fourier-transformed electroweak gauge boson kinetic terms,L µν V V , of order α: Analyzing Eqs. (2.67) -(2.70) and taking into account the renormalization conditions of [58], we see immediately that Apart from light fermion and top quark contributions which we consistently neglect throughout this work, explicit expressions for the one-loop self-energies and counterterms relevant to the higher-order corrections of interest to us will be provided to all orders in in Section A.2. Note that we do not provide an explicit expression for δZ (1,0) W ± because it plays no role in the renormalization of any of the Feynman diagrams we calculate.

On-shell wavefunction renormalization for massless fermion fields
In this section, we review the systematics of on-shell renormalization for massless fermion fields in the Standard Model. The on-shell scheme for massless fermions is widely used in QCD because, due to the scalelessness of the contributing gluon-exchange diagram, the one-loop quark wavefunction counterterm vanishes identically in dimensional regularization. Although we will ultimately take the q 2 → 0 limit, it is nevertheless instructive to begin with the off-shell setup of [58], where the momentum transfer q 2 in the Lorentz decomposition of the fermion self energy is taken to be different from zero. 4 In writing Eq. (2.72), we have implicitly assumed Kreimer's γ 5 scheme; analogous calculations in HVBM's γ 5 scheme will be discussed in Section 2.8 and, accordingly, we have writtenΣ f (q) instead of Σ f (q) in order to allow for a side-by-side comparison in the next section.
We employ the notation of [58] for the interactions of the electroweak gauge bosons with matter. The Z interaction is parametrized by flavor-dependent axial vector and vector coupling coefficients, The W interaction, on the other hand, is parametrized by a universal coupling coefficient We will make extensive use of these aliases throughout the rest of this paper. At first sight, it might seem counterproductive not to employ Eq. (2.74) to eliminate a w and v w , as they are equal and have no dependence on the fermion flavor. However, as will become clear later on in this section, it is useful to retain the dependence on a w and v w at intermediate stages of fermion self-energy and vertex calculations because it generically allows for a determination of axial vector components from vector components through chiral symmetry. Actually, the fermion self-energies considered in this section are particularly simple and have no independent W -exchange Feynman diagrams; all Wexchange diagrams can be obtained from the Z-exchange diagrams for free by making the replacements a f → a w and v f → v w (see Figure 2).
Thus, for the relatively low perturbative orders of interest to us, it suffices to consider the calculation of the Z-exchange contributions to the vector component of the fermion selfenergy. At this stage, one is faced with a conceptual hurdle. Taken at face value, the selfenergy in the on-shell scheme for massless fermions would seem to be most appropriately defined as the q 2 → 0 limit of the fermion self-energy calculated off-shell. 5 Of course, it also seems that one ought to be able to simply set q 2 = 0 in all Feynman diagrams at the very beginning, significantly simplifying all subsequent calculations. In fact, as we now demonstrate at the one-loop level, both approaches are possible and yield the same answer.
If we begin with q 2 generic, one can immediately verify that the Lorentz projector of interest is and our Kreimer reading point prescription instructs us to begin reading our Dirac traces from the projector insertion. Therefore, using the Feynman rules of [58] and the anticommutation relation satisfied by Kreimer's γ 5 , Eq. (2.7), the contribution of the one-loop Z-exchange diagram to the vector component of the fermion self-energy is As Eq. (2.76) does not depend on γ 5 , it can be trivially simplified using the contraction identity and even trace formulae of Section 2.2. After carrying out the numerator algebra, we obtain: , (2.78) which can immediately be rewritten as a linear combination of standard scalar integrals using Passarino-Veltman reduction [55]; the result is where, using Feynman/Schwinger parameters, In writing Eq. (2.81), q 2 was assumed to be less than m 2 z . Modulo minor typos, 7 our results agree as expected with those of [58] through to O 0 in the limit of vanishing fermion mass.
From Eqs. (2.80) and (2.81) it is a simple matter to work out the q 2 → 0 limit of Eq. (2.79). From l'Hôpital's rule, we find: The above derivation is certainly valid, but let us stress that attempting to calculate the q 2 → 0 limit ofΣ V, f (q 2 ) in this way is not recommended at higher orders in perturbation theory; while we were able to confirm numerically that the correct q 2 → 0 limit exists at order αα s for fixed, finite , it is not amenable to analytic evaluation. In general, our opinion is that it is far simpler to implement the on-shell condition q 2 = 0 from the get-go. Of course, Eq. (2.75) is ill-defined at q 2 = 0; a valid projector when q 2 = 0 is given by where η is an arbitrary four-vector such that η · q = 0. Setting q 2 = 0 from the beginning of the calculation, we havē . (2.84) Proceeding directly from Feynman/Schwinger parameters for the integral evaluations, we find 7 While comparing to [58], we noticed a typo in the expression for Σ iσ A on the 3rd line of Eq. (5.27): the first term inside the bracket should actually read −2 viσaiσ 2 B1 k 2 ; miσ, Mz + 1 . Additional typos were discovered in Eq. (B.2); the tensor reduction formula needed for B1 k 2 ; miσ, Mz should read Note that, in the degenerate q 2 = 0 kinematics under consideration, all contributing oneloop Feynman integrals must ultimately turn out to be proportional to the one-loop tadpole. Finally, we obtain and see immediately that as claimed.
We now consider the ramifications of chiral symmetry. In the Standard Model, the right-and left-handed components of the Z-exchange contribution to the one-loop fermion self-energy are exchanged under γ 5 → −γ 5 . This implies that the axial vector component of the Z-exchange contribution to the fermion self-energy ought to be obtained by a re- In order to fix the overall sign, recall that the coupling of the W to matter is exactly left-handed, that the W -exchange contribution is obtained by making the replacements a f → a w and v f → v w , and that a w = v w . Taken together, these considerations allow us to conclude that with no additional calculation.
In the on-shell scheme, massless fermion wavefunction counterterms are defined to exactly cancel the perturbative corrections to the corresponding fermion self-energies orderby-order. Therefore, we find for the vector and axial vector components of the one-loop fermion wavefunction counterterm.
We now turn to the vector and axial vector components of the two-loop quark selfenergy of order αα s . To our knowledge, these mixed two-loop corrections to the self-energy were first considered in [61], though not in the on-shell scheme with massless fermions. It is straightforward to calculate the desired quantities by acting with the appropriate Lorentz projector in Kreimer's γ 5 scheme, Eq. (2.83), on the eight non-zero two-loop Feynman diagrams summarized by Figure 2. Remarkable simplifications are achieved by setting q 2 = 0 at the beginning of the calculation. In the end, we find that the vector component of the quark wavefunction counterterm is given by and the axial vector component of the quark wavefunction counterterm is given by As at one loop, the simple form of Eqs. (2.92) and (2.93) results from the fact that only the one-mass two-loop tadpole integral can survive in the q 2 → 0 limit.

Chiral symmetry restoration in HVBM's γ 5 scheme and Larin's principle
As mentioned in Section 2.2, HVBM's γ 5 scheme is especially challenging to work with in practice due to the fact that, in general, a plethora of finite counterterms must be introduced to restore the chiral symmetry of the Standard Model. We shall see that Larin's principle [48] provides a practical recipe which helps tremendously to minimize the number of finite counterterms which must be computed explicitly for a given application. While we will be most interested in the standard one-loop example of the order α s gluonexchange correction to the Zqq vertex shown in Figure 4 [40], let us first revisit the one-loop fermion self-energy calculation of Section 2.7 in HVBM's γ 5 scheme to illustrate the utility of Larin's principle with minimal computational effort. In HVBM's γ 5 scheme, the Lorentz decomposition for the massless fermion self-energy reads where we have assumed the q 2 = 0 Lorentz projector, Eq. (2.83), will be employed to compute Σ V, f (0) and, as in Section 2.7, the color indices are suppressed. It is once again sufficient to consider the Z-exchange contribution: (2.95) In order to simplify the term proportional to a 2 f , it is useful to resolve the d-dimensional momenta and Dirac matrices into four-dimensional and (−2 )-dimensional components using Eqs. (2.10) and (2.11). Crucially, the linear dependence on k 1 in the numerator of Eq. (2.95) implies that the Dirac traces cannot produce non-zero µ integrals. This observation allows for the immediate neglect of all terms involvingk 1 .
Using the HVBM γ 5 anticommutation relations of Section 2.2, Eq. (2.95) simplifies to Finally, the trace on the second line of Eq. (2.96) is nothing but the → 0 limit of the trace on the first line of Eq. (2.96); from the analysis of Section 2.7, we see that (2.97) Starting at O 0 , Eq. (2.97) clearly exhibits a chiral mismatch, indicating that the vector components of the fermion kinetic terms of the Standard Model Lagrange density receive finite counterterm corrections already at order α. In other words, a correction must be made because the coefficient of a 2 f is not equal to the coefficient of v 2 f in Eq. (2.97) at O 0 . As Larin emphasized in [48], the form of the finite counterterm may be deduced by simply demanding that the a 2 f and v 2 f coefficients be equal. For the present calculation, Larin's principle dictates that the finite counterterm is given to all orders in by the differencē In fact, we could have simply skipped the calculation of the bare self-energy in HVBM's γ 5 scheme altogether; requiring that our result respect the chiral symmetry of the Standard Model to all orders in uniquely leads to the result obtained in Kreimer's γ 5 scheme in Section 2.7 after finite renormalization. 8 Of course, Larin's principle can be successfully applied to the fermion self-energies and vertex form factors because they have a direct correspondence to terms in the renormalized 8 In [48], it is suggested that one ought to perform additional renormalizations for the purpose of chiral symmetry restoration only after the removal of all ultraviolet divergences with ordinary counterterms. In our view, this is not necessarily the most transparent approach. However, at higher orders in perturbation theory, it may be that restoring the chiral symmetry first in the way we propose leads to divergent "finite" counterterms. Therefore, in our approach, "finite" may not be the best moniker to generally describe the counterterms responsible for the restoration of the Standard Model chiral symmetry in HVBM's γ5 scheme. Lagrange density. For more complicated contributions to higher-multiplicity scattering amplitudes at two loops and beyond, it is expected that non-trivial corrections will arise from the finite counterterms, beginning with insertions into lower-loop box-type diagrams. Indeed, as we shall see later on in our discussion of the two-loop mixed EW-QCD corrections to the neutral-current Drell-Yan process, certain two-loop box-type diagrams which involve axial vector couplings to quarks receive corrections from finite counterterm insertions into corresponding one-loop box diagrams (see Figure 3). The bare calculation in HVBM's γ 5 scheme needed to derive the finite counterterm featured in Figure 3 is a familiar one [40,48] and, fortunately, the only one relevant to this work which we cannot bypass.
In HVBM's γ 5 scheme, the Zqq vertex has the Lorentz decomposition where s = (p 1 + p 2 ) 2 denotes the virtuality of the Z and color indices are suppressed for brevity. Here, we have assumed an incoming massless quark of momentum p 1 and an incoming massless antiquark of momentum p 2 for later convenience. The Lorentz projectors  whereV (0,1) Zqq (s) is nothing but v q times the bare time-like one-loop quark vertex form factor of massless QCD: The derivation of the analytic expression forF q 1 ( ) given above is reviewed in many places in the literature, see e.g. [62]. Thus,Ā Clearly, the remaining Dirac traces in Eq. (2.107) may be evaluated in a straightforward manner by applying Eqs. (2.5), (2.25), and (2.26) with set to zero and the d-dimensional metric tensor replaced by the four-dimensional metric tensor throughout. Discarding all contributions which vanish due to scalelessness and settingk 1 ·k 1 = k 2 1 −k 1 ·k 1 , we find (2.109) It is worth pointing out that, by construction, such explicit µ terms could not occur in a direct calculation ofĀ (0,1) Zqq (s) in Kreimer's γ 5 scheme. Finally, taking into account the algorithms of [56] and their implementation in Reduze 2 [63][64][65] to treat the dimensionallyshifted integral on the right-hand side of Eq. (2.109), it follows that after a straightforward integration by parts reduction. As is well-known (see e.g. [66]), and therefore Now, as explained above, the finite counterterm of interest is obtained as the difference between the Kreimer's γ 5 scheme expression deduced above and Eq. (2.112): Eq. (2.113) is manifestly finite at = 0; as expected, we see that δZ (0,1) Zqq = 4a q C F + O( ), in complete agreement with Eq. (9) of [48]. 9 9 In making this comparison, note that Eq. (3) of [48] was not written down for the Zqq interaction verbatim and therefore does not have the overall factor of −aq present in interaction (2.24). Note also that the cited equation numbers refer to the journal version of [48].
Before leaving this section, let us stress that we do not employ in this paper what is commonly referred to as "Larin's scheme" in the QCD literature. It is clear that the goal of [48] was to provide a recipe in FORM for the consistent treatment of multi-loop QCD corrections to the axial vector components of Standard Model vertex form factors in HVBM's γ 5 scheme; Larin did not discuss a prescription for the treatment of multiloop box-type diagrams, for example. Indeed, in contrast to the prescription of [48] for the treatment of multi-loop QCD corrections to Standard Model vertex form factors, we found it essential for the restoration of chiral symmetry in the two-loop box-type diagrams treated below to consistently keep also higher-order terms of the expansions of our finite counterterms in order to avoid any premature truncation.

Process definition
In this work, we study lepton pair production in quark-antiquark annihilation, at one and two loops. We obtain complete results for scattering amplitudes at order αα s , order α 2 , and order α 2 α s in the approximation of vanishing light quark and lepton masses, i.e. p 2 1 = p 2 2 = p 2 3 = p 2 4 = 0. As the present paper represents only a first decisive step towards the assembly of phenomenological results for Drell-Yan lepton production including complete two-loop EW-QCD effects in all channels, we defer the calculation of the gauge-invariant contributions proportional to the number of light quark or lepton flavors and all top quark mass corrections to future work. Our expectation is that the most important two-loop corrections we omit, coming from the light fermion contributions to the electroweak gauge boson self-energies, are already known (see e.g. [67]).
As usual, the 2 → 2 kinematics of (3.1) is conveniently described by the Mandelstam variables In total, five independent kinematic scales, s, t, m w , m z , and m h , appear in our results.
Here, m w is the mass of the W boson, m z is the mass of the Z boson, and m h is the mass of the Higgs boson. For simplicity, we present our final results in a physical phase-space region where t is allowed to take arbitrary physical values while s is constrained to be above all two-particle thresholds present in the contributing Feynman diagrams. This amounts to the condition imposed by the Z-Z self-energy. Working in the kinematic region defined by (3.3) allows us to do without the complex mass scheme, thereby avoiding further technical complications for our proof-of-concept study. As discussed in [24], this restriction also simplifies the explicit expressions for the two-loop master integrals.

Form factor decomposition
We decompose the one-and two-loop scattering amplitudes for our process into basic Lorentz structures ("standard matrix elements") multiplied by form factors, which depend only on the kinematic invariants. Since we are considering massless external fermions in the Standard Model, the amplitude in conventional dimensional regularization can be written up to two loops in terms of the Lorentz structures [68] where we suppressed color and spin indices and ignored insertions of γ 5 . For our EW-QCD corrections, structures T 5 and T 6 can be omitted even in a general R ξ gauge because the gluon couples only to the initial-state quark line.
Allowing also for insertions of γ 5 at the end of each Dirac chain, we obtain sixteen different Lorentz structures. For the sake of our argument, let us consider an anticommuting γ 5 to illustrate the general principle of our projector method. The described form of the amplitude can be obtained either using Passarino-Veltman reduction of the tensor integrals or using projectors. In d = 4 dimensions, only four out of the sixteen Lorentz structures are independent, as can be seen using Schouten and Fierz identities. In order to facilitate a smooth d → 4 limit, we can define a d-dimensional basis of the four structures derived fromT 1 by inserting γ 5 matrices: (3.10) and twelve further structures (derived fromT 2 ,T 3 , andT 4 in a similar way), which are strictly "orthogonal" in d dimensions and smoothly vanish for d → 4. The d dimensional projection onto the structures (3.10)-(3.13) obtained in this way are identical to those calculated by ignoring the evanescent directions at the outset. This procedure may be justified by considering the fact that, even in d dimensions, helicity is conserved and, as is clear from Eq. (2.26), the spin space retains its usual (integer) number of degrees of freedom, see also [69][70][71]. We note that the four independent structures in d = 4 are of course directly related to the four independent helicity amplitudes; alternatively, for a given γ 5 scheme, we could have written polarized Lorentz structures by replacing e.g. u(p i ) with For finite, renormalized and infrared-subtracted scattering amplitudes, O( ) deviations from the four-dimensional projectors will be irrelevant for the O 0 result. For infrared divergent amplitudes, the O 0 term depends on the exact choice of the projector. However, as long as we use the exact same prescription for the projector for all divergent contributions, those ambiguities will exactly cancel for a finite sum. This argument assumes that we use an analytic continuation which gives a well-defined meaning to the higher-order-interms of the divergent quantities and that we do not truncate their Laurent expansions in prematurely. In the schemes we consider here, for contributions involving Levi-Civita tensors, we expect that consistency is ensured because Levi-Civita tensors effectively truncate the higher-order-in-terms of finite quantities. We will discuss the details of how we construct our projectors for the γ 5 schemes we consider in the following.

Calculations in HVBM's γ 5 scheme
The Lorentz decomposition of the Drell-Yan scattering amplitude in HVBM's γ 5 scheme is with the Lorentz structures In the manner described in Section 3.2, we employ Lorentz projectors to calculate the form factors C VV , C AA , C VA , and C AV in d dimensions. We request C = spin,color and obtain In practice, we use these projectors only for box-type diagrams; for diagrams with selfenergy or vertex insertions, we choose to sew effective propagators and vertices into appropriate tree-level diagrams. Note that, as we plan to perform a calculation of all Feynman diagrams in Kreimer's γ 5 scheme anyway, Larin's principle allows us to completely bypass the direct calculation of almost all vertex Feynman diagrams in HVBM's γ 5 scheme; as explained in Section 2.8, the vertex form factor projector defined in Eq. (2.101) above was employed for the calculation of the finite counterterm δZ (0,1) Zqq of Eq. (2.113), but that is all. 10 After carrying out the renormalization program outlined in Sections 2.5-2.8, the looplevel amplitude form factor coefficients still diverge in the infrared. These singularities in the virtual amplitudes manifest as poles in which are canceled by taking real radiative corrections into account and, in the end, considering perturbative corrections to physical observables. Fortunately, we have a priori knowledge of the infrared singularity structure of the scattering amplitudes we calculate in this work due to the well-known dipole formula [72][73][74] and a simple generalization thereof. In Section 3.5, we review the two-loop mixed EW-QCD neutral-current Drell-Yan dipole singularity structure derived in [13] and show how subtraction functions may be defined which allow for the comparison of finite hard scattering functions in different γ 5 schemes.
We now turn to a technical discussion of our approach to Dirac trace manipulation and numerator algebra in Mathematica. First of all, note that all terms which involve an odd number of γ 5 matrices in total must vanish due to the fact that we consider a 2 → 2 scattering process with purely massless external momenta which cannot support a contraction with a Levi-Civita tensor. One possible approach to calculate the non-vanishing terms in HVBM's γ 5 scheme, popularized by [48], would be to make use of the identity While the associated performance penalty could be ameliorated to some extent by exploiting the total antisymmetry of the Levi-Civita tensor, we prefer instead to follow [44] and split indices from the very beginning. We use the anticommutation relations, Eqs. (2.20) and (2.21), to write 11 Once all explicit γ 5 matrices have been anticommuted away and traces with an odd number of hatted variables have been discarded, we apply Eq. (2.5) and its companion identity, (3.26) to shorten all Dirac traces at an early stage to the maximal extent possible. 10 Strictly speaking, there is another finite counterterm which is needed as well, δZ (0,1) W ±q q , but it may be trivially obtained from Eq. (2.113) via the replacement aq → aw. 11 To be clear, we distilled Eq. (3.25) from [40]; axial vector couplings were not considered in [44].
In our view, this strategy is most sensible because, due to the presence of box-type diagrams, it is anyway inevitable that we must split indices at some stage. As a bonus, Eq. (3.26) also allows us to exploit the on-shell kinematics to discard vanishing terms at an early stage. While it may not be obvious at first sight, Eqs. (2.5), (2.25), and (3.26) can be applied verbatim to traces which contain both barred and hatted variables, provided that mixed contractions of indices are discarded as soon as they are produced by the formulae. Although it was not really needed for the two-loop calculations considered in this work, it is worth mentioning that an aggressive use of memoization for the Dirac traces appearing in our calculations could have significantly decreased the runtime required by our code to process all Feynman diagrams.
Once all Dirac traces have been evaluated and all compatible Lorentz indices contracted, contractions and products of Levi-Civita tensors from the box-type diagrams must be dealt with using Eqs. (2.12) -(2.16). Subsequently, all scalar products of the formḡ µν k µ i k ν j must be processed using the methods reviewed in Section 2.4. Finally, taking into account the algorithms of [56] and their implementation in Reduze 2, we obtain the amplitude form factor coefficients via integration by parts reduction as a linear combination of master integrals. At both one and two loops, we find it particularly convenient to rotate to normal form bases of Feynman integrals, characterized by the simple, -decoupled differential equations they satisfy [75,76]. Our one-and two-loop normal form bases are defined in Sections 4.2 and 5.2.

Calculations in Kreimer's γ 5 scheme
The Lorentz decomposition of the Drell-Yan scattering amplitude in Kreimer's γ 5 scheme is essentially Eq. (B.10) of [58]. Suppressing the color indices, we have PC iĀ DY , forC =C VV ,C AA ,C VA ,C AV , (3.28) and are given by  (3.35) are useful. 12 Of course, a decomposition completely analogous to (3.33) exists for the lepton vertex form factors, and P μ  We calculate our one-and two-loop scattering amplitudes in Mathematica, using a code pipeline similar to that described in Section 3.3; as mentioned in Section 2.3, the algebraic form of the non-vanishing families of even and odd Dirac traces is the same in Kreimer's γ 5 scheme as in HVBM's γ 5 scheme. A key advantage of Kreimer's γ 5 scheme, however, is that, due to the standard anticommutation relation, Eq. (2.7), Dirac traces need not be calculated with split indices. Rather, one can first evaluate all traces and introduce scalar products of the formḡ µν k µ i k ν j only at a later stage, when contractions and products of Levi-Civita tensors are eliminated using Eqs. (2.12) -(2.16). As a cross-check, we also calculated the interference of the one-and two-loop scattering amplitudes with the tree-level Drell-Yan amplitude in Kreimer's γ 5 scheme using an independent setup in FORM 4 [77] which, in particular, employed Passarino-Veltman reduction [55] instead of the technology of Section 2.4 to simplify scalar products of the formḡ µν k µ i k ν j . Finding complete agreement between our independent implementations served as a crucial cross-check on our calculations.

Infrared dipole singularity structure and subtraction functions
The dipole formula [72][73][74] provides a particularly concise and straightforward recipe for the generation of infrared subtraction terms for virtual gauge theory scattering amplitudes calculated in pure dimensional regularization. While, historically, the singularity structures of scattering amplitudes in particularly simple models such as massless QED, massless QCD, and N = 4 super Yang-Mills theory were discussed most frequently in the literature, broader applications to theories with massive particles are certainly possible and have long been known (see e.g. [78]). Indeed, it was shown in [13] that with straightforward modifications, the dipole formula can be extended to describe the singularity structure of the two-loop mixed EW-QCD neutral-current Drell-Yan scattering amplitude of interest.
We begin by introducing the building blocks required to describe the infrared singularities of the order αα s , order α 2 , and order α 2 α s neutral-current Drell-Yan scattering amplitudes. As has long been clear, the leading infrared singularities of gauge theory scattering amplitudes are governed by cusp anomalous dimensions [79]. The one-loop quark cusp anomalous dimension of massless QCD may be extracted, for example, from the expression for the one-loop quark form factor of massless QCD given in Eq. (2.103). In QCD, we have (3.39) whereas in QED we obtain the result from Eq. (3.39) by replacing the quadratic Casimir invariant C F with the squared charge of the fermion flavor, However, as pointed out in [13], the mixed quark cusp anomalous dimension vanishes: The next-to-leading infrared singularities are more complicated and receive contributions from both soft anomalous dimensions and resummation functions derived from the −1 poles of massless vertex form factors with an appropriate number of gluon and/or photon exchanges. The one-loop Drell-Yan soft anomalous dimensions are well-known and have been calculated many places in the literature, see e.g. [80,81] for a thorough discussion. 13 Rewriting the results to make all imaginary parts explicit in the physical kinematic region of interest, we find and S (1,0) (3.43) 13 Note that, for general QCD scattering processes, the soft anomalous dimension becomes a mixing matrix in color space. For processes with just two partons in the initial state such as the Drell-Yan process, the color space structure trivializes and all matrices in the dipole formula may be replaced with functions.
Once again, as shown in [13], the two-loop mixed EW-QCD Drell-Yan soft anomalous dimension vanishes identically: We define our massless QCD resummation functions in the framework of [82]. In practice, following [13], all results may be derived from the one-and two-loop quark resummation functions of massless QCD by making simple replacements. In what follows, f [k] denotes the coefficient of k in the Taylor series expansion of f ( ), (3.45) To the relevant orders, we have Finally, the two-loop mixed EW-QCD quark resummation functions are obtained from Eqs. (3.50) and (3.51) by setting C A to zero and replacing C 2 F with Q 2 q C F : and Given the ingredients discussed above, we are now in a position to present the predictions of the generalized dipole formula for the ultraviolet-renormalized, infrared-divergent scattering amplitudes A DY andĀ DY in, respectively, HVBM's γ 5 scheme and Kreimer's γ 5 scheme. At tree level we have (3.58) where the dots stand for terms of higher order in the coupling constants. As the notation suggests, the tree-level hard functions, H in the notation of [58] for the couplings of the Z boson to matter, see Eqs. (2.73).
As is clear from the form of Eqs. (3.59), the tree-level amplitudes in the two γ 5 schemes we consider trivially coincide in the four-dimensional limit because all of the differences between the schemes disappear in four spacetime dimensions: For the remainder of this work, it will sometimes be convenient to view the hard scattering functions in the two schemes we consider as vectors of coefficients with respect to where the one-and two-loop hard scattering functions which appear above are determined by demanding equality between Eqs. (3.62) and the corresponding explicit calculations at one and two loops. All other things being equal, it is generally believed that HVBM's γ 5 scheme and Kreimer's γ 5 scheme produce identical results through to one loop and O 0 [46]. This proposition was implicitly affirmed at order αα s in Section 2.8. At order α 2 , our explicit one-loop calculations, discussed in detail in the next section, do indeed show that However, we also find that for k > 0 . In the Standard Model, it is expected that all terms of virtual and real radiative corrections sensitive to the γ 5 problem eventually cancel out of combinations which furnish complete fixed-order perturbative corrections to physical observables [83]. Therefore, it is expected that (3.65) actually is an equality. We will test this expectation with our explicit calculations below.

Diagrammatic structure
In this section, we present the diagrammatic structure of the one-loop perturbative corrections to the neutral-current Drell-Yan process. Due to their familiarity, the renormalized self-energy corrections to the photon and Z propagators, , will not be visualized in explicit detail. In what follows, we will also use a generic wavy line for both the photon and the Z, keeping in mind that some contributions involving e.g. a coupling to neutrinos can occur only for the Z. As is well-known, just a single class of diagram contributes to the one-loop corrections of relative order α s , .
The remaining one-loop corrections of relative order α fall into three categories: initial state vertex diagrams, final state vertex diagrams, and box diagrams. Both the initial state vertex diagrams and the final state vertex diagrams require renormalization, , due to the presence of ultraviolet divergences. The box diagrams have infrared divergences only.

One-loop integral definitions
All of the one-loop integrals defined below are pure functions taken in the standard MS normalization. Therefore, all of our one-loop Feynman integrals are understood to have a multiplicative factor of in the integration measure. In what follows, thin solid lines denote massless propagators or massless external momenta, thick, dotted lines denote propagators of mass m w , thick, dashed lines denote propagators of mass m z , and thick solid lines denote propagators of mass m h . In total, 31 linearly-independent one-loop integrals appear in our calculations: where λ(x, y, z) is the Källén function, For clarity and later convenience, we have made all dependence on the kinematic variables s, t, u = −s − t, m 2 w , m 2 z , and m 2 h completely explicit on the right-hand side. Moreover, we have included integrals for all relevant permutations of the kinematic variables separately. In order to actually calculate the master integrals, it is enough to find explicit expressions for I 1 , I 2 , I 3 , I 6 , I 7 , I 8 , I 13 , I 14 , I 16 , I 17 , I 19 , I 22 , I 24 , and I 28 , as all other integrals may be obtained by making simple replacements. For example, I 31 is obtained from I 28 by replacing m 2 z with m 2 w and t with u. Due to the fact that we consider the physical region s > (m z + m h ) 2 (see (3.3)), the same form of the analytic result can be used in a straightforward way.
Our calculation requires knowledge of the most complicated one-loop box-type master integrals (i.e. I 24 and I 28 ) expanded through to fourth order in . We could not find suitable analytic solutions for all integrals in the literature. Several of the explicit results we could locate in the Loopedia database [84] were either not expanded to a sufficiently high order in for our purposes or not provided in a form convenient for numerical evaluations in the physical region of interest to us.
We therefore computed all of the integrals from scratch, either by direct integration for generic followed by an expansion 14 in or by using the method of differential equations [41,42,[86][87][88][89][90]. The integral definitions given above lead to a -decoupled form for the differential equations [75,76], which we integrate in terms of multiple polylogarithms. To proceed, we construct real-valued multiple polylogarithms in the region of phase space of interest to us, employing the functional basis presented in [22,24], together with a few additional logarithms and polylogarithms required to integrate e.g. I 6 to a sufficiently high order in . We provide our results in the ancillary file oneloopmasters.m attached to the arXiv submission of this paper.

Assembly of one-loop results
Assembling the components provided in Appendix A.1 and A.2, we find (1,0) V, 14 For some of the more complicated expansions, we used the Mathematica package HypExp [85].
for the renormalized one-loop scattering amplitude coefficients in HVBM's γ 5 scheme. In writing the above expressions, we made use of the tabulated one-loop vertex counterterms collected in Appendix A of [59]. Results in Kreimer's γ 5 scheme may be simply obtained from Eqs.

Diagrammatic structure
In this section, we present the diagrammatic structure of the two-loop perturbative corrections to the neutral-current Drell-Yan process of relative order αα s . As at one loop, we do not explicitly identify the photon and Z. Due to their familiarity, diagrams with one-loop renormalized self-energy insertions of the form will not be visualized in explicit detail. Several other two-loop diagrams are also trivial, in that they essentially consist of simple products of one-loop contributions: . Note that, in the above, the sum of the final four diagrams vanishes identically. For the sake of brevity, we draw only one half of the two-loop vertex diagrams with a single massive vector boson stretched across the quark line: .
The remaining diagrams of this type may be trivially recovered from the above by exchanging the gluon and the electroweak gauge boson attached to the quark line.
Of course, there are also diagrams with γW + W − and ZW + W − interactions: .
Finally, the order αα s vertex receives a correction from a two-loop vertex counterterm insertion, .
The most complicated two-loop diagrams are those of box type: .
Our evaluation of these box-type diagrams is one of the most important new results of this work. In HVBM's γ 5 scheme, further diagrams, , need to be included to ensure that the order α 2 α s hard scattering function in HVBM's γ 5 scheme respects the chiral symmetry of the Standard Model (see Section 2.8 for more details).

Two-loop integral definitions
All of the two-loop integrals defined below are pure functions, converted from the idiosyncratic notation of [22,24], where they were originally evaluated in the physical region above all two-particle thresholds, to standard MS normalization. To be explicit, the m i integrals from [22] and the m i integrals from [24] must be multiplied by (5.1) That is to say, the integral measure of our two-loop integrals is exactly: integrals not treated explicitly in references [22] or [24] may be directly obtained through to weight four from the integral evaluations carried out in [22,24] via the replacement t → u. For the sake of completeness, we present a simple expression for J 29 in the physical region above its one-mass threshold. Indeed, the integral admits a compact evaluation in Feynman/Schwinger parameters to all orders in : Expanding Eq. (5.4) in , we find results equivalent to those obtained from J 62 by performing a crossing from t to s and then carefully rewriting the result in terms of the basis of multiple polylogarithms suitable for the expansion of J 28 in the physical region above the one-particle threshold.
Of course, it is crucial to bear in mind that the renormalization program in HVBM's γ 5 scheme is complicated by the fact that insertions of the finite counterterms 17 δZ (0,1) Zqq 16 We remind the reader that explicit all-orders-in-results for the components of the quark wavefunction counterterm at order ααs are provided in Eqs. (2.92) and (2.93). 17 Expressions for the finite counterterms relevant to this work are provided in unintegrated form in Eqs. is a highly non-trivial relation.
As the equality sign in (5.10) suggests, after applying the computational methods and master integral evaluations described above to the Feynman diagrams of Section 5.1, our explicit two-loop calculations confirm that (3.65) does in fact hold: the hard scattering functions for the two-loop mixed EW-QCD corrections to the neutral-current Drell-Yan process we consider in HVBM's γ 5 scheme and in Kreimer's γ 5 scheme coincide. In our view, the rigorous establishment of Eq. (5.10) constitutes a very important cross-check on our entire formulation and a significant theoretical result. In the future, our detailed comparison of HVBM and Kreimer at the two-loop level should benefit further calculations of still-unknown multi-loop perturbative corrections to more complicated Standard Model processes in pure dimensional regularization.
In future studies of Standard Model scattering processes, it would certainly be possible to perform two parallel calculations for each process in the two γ 5 schemes considered in this work, but there may be a more practical way to rigorously cross-check perturbative calculations going forward. During the course of our investigation of qq → + − in Kreimer's γ 5 scheme, we discovered that we could effectively check the correctness of the two-loop hard scattering functions within Kreimer's γ 5 scheme itself by setting up our code pipeline to allow the user to make different choices of reading point prescription for the Dirac traces. For example, we found that the fully-symmetrized reading point prescription mentioned in Section 2.3 led to different higher-order-in-one-loop hard scattering functions, but the same order α 2 α s hard scattering function at O 0 .
The two-loop hard scattering functions calculated in this work are too long to be included in the article. Once we have implemented support for the remaining kinematic regions of the physical phase space, we plan to provide public access on HEPForge to an efficient C++ program for the evaluation of the polarized hard scattering functions discussed in the next section.

Helicity amplitudes
Having arrived at finite quantities in four dimensions, we can directly cast our results into expressions for the scattering of fermions with definite helicities using explicit fourdimensional representations.
Due to the form of the massless fermion propagator and the manner in which the electroweak gauge bosons couple to matter, the helicity quantum number is conserved along the initial and final fermion lines (see e.g. [91] for a review). This implies that there are just four non-vanishing helicity amplitudes, those where both the quark and the antiquark have opposite helicity and the lepton and antilepton have opposite helicity. We parametrize the center-of-momentum frame according to sin θ, 0, cos θ), and use an explicit basis of spinor states in the chiral representation [39] u + (p 1 ) = s 1/4 0, 0, 1, 0 , Plugging these representations 18 into our Lorentz structures directly gives us results for the helicity amplitudes. Decomposing the hard functions (or finite remainders) according to for coupling orders m and n. In Eqs. (6.5) -(6.8), due to the observed γ 5 schemeindependence of the hard scattering functions at zeroth order in the expansion (see Section 5.3), we drop the dual notation employed in Eq. (6.4).

Final results
In this section we present visualizations of the polarized hard scattering functions H (0,0)

and H
(1,1) λ 1 λ 2 λ 3 λ 4 defined in Section 6.1 for all four non-trivial helicity 18 The phase conventions of (6.3) above are in line with the spinor helicity formalism, as in [92] where u±(pi) = v∓(pi) for all i. Let us stress, however, that we do not adopt the all-outgoing convention for the external four-momenta. We made use of Eqs. (6.1) and (6.2) to derive Eqs. (6.3). configurations. In the numerical analysis, we set the gauge boson and Higgs masses to their on-shell values as listed in [93], m w = 80.379 GeV, m z = 91.1876 GeV, m h = 125.10 GeV, (6.9) and the renormalization scale µ R = √ s. Since we keep powers of α and α s factored out of the expressions we plot, we do not need to specify their values here. For illustrative purposes, we found it sufficient to consider up-type quarks in the initial state. Our figures do look rather different for down-type quarks, but our impression is that, on the whole, they do not introduce completely new features which would be of paramount importance to discuss here. We elected to plot only the real parts of our final results, as the real parts of H (1,1) λ 1 λ 2 λ 3 λ 4 contain the most complicated weight four multiple polylogarithms appearing in our calculations.
In our numerical analysis, we focus on larger values of √ s where the previously unknown non-factorizable two-loop box-type Feynman diagrams of Section 5.1 are expected to be important. In order to compare the different orders in α and α s , we find it convenient to include in our plots the 4π suppression factors taken out of each higher-order term in Eq. (3.62); that is, we consider H (m,n) λ 1 λ 2 λ 3 λ 4 /(4π) m+n , where the relative orders in the electroweak and strong couplings are given by (m, n) (see also Eq. (2.56)).
In all plots, we include the tree-level and relative order α s results for reference, they show a rather simple behavior. The one-loop QCD corrections lie almost on top of the tree level results. Indeed, they fully factorize from the tree results; we have for all helicity configurations and both up-and down-type quarks in the initial state. For the real part we see that π/3 ≈ 1.05, thus explaining the observed similarity. Figure 5 shows the dependence of the hard functions on the center-of-mass energy √ s for fixed central scattering angle. We observe that the absolute values of all plotted real parts of relative order αα s increase as a function of √ s, justifying the expectation that the calculation of the mixed two-loop EW-QCD corrections for the leptonic final state is well-motivated in the kinematic regime depicted in the figure. We note that the real parts of the one-loop EW and the two-loop EW-QCD corrections are not aligned for all helicity configurations. Figure 6 shows the dependence of the hard functions on the cosine of the scattering angle for fixed center-of-mass energy √ s. The electroweak corrections show a rather complex angular dependence. While the one-loop EW and the two-loop EW-QCD corrections show a similarity in their angular dependence, they differ in the details.
We also compared the curves of Figure 6 to analogous ones for the pure QED-QCD model studied in [13]. 19 While maybe not instructive from the phenomenological point of view, it was interesting to see a somewhat similar qualitative angular dependence emerge when comparing the two-loop polarized hard scattering functions for QED-QCD and EW-QCD normalized by their respective tree level contributions.

Summary and outlook
In this article, we calculated the relative order αα s mixed EW-QCD corrections to Drell-Yan lepton pair production, qq → + − . We performed the calculation in both the 't Hooft-Veltman-Breitenlohner-Maison γ 5 scheme and in the Kreimer γ 5 scheme using the projector method for chiral fermions. While our two-loop scattering amplitudes were found to be scheme dependent starting at O −1 , unique polarized hard scattering functions in d = 4 were obtained after infrared subtraction. In the 't Hooft-Veltman-Breitenlohner-Maison scheme, we restored chiral symmetry using local counterterms, where we found it essential to consistently calculate our renormalization constants to higher orders in . To the best 19 To our knowledge, the authors of [13] calculated unpolarized hard scattering functions only. of our knowledge, our application of Kreimer's γ 5 scheme to the calculation of genuine two-loop 2 → 2 scattering amplitudes with non-trivial electroweak effects is the first of its kind.
Our calculation provides a major building block for the calculation of the relative order αα s corrections to off-shell Drell-Yan production in the high energy region, which is particularly relevant to new physics searches and the establishment of constraints on possible extensions of the Standard Model.

A.2 Relative order α results
Corrections to the neutral-current Drell-Yan process of relative order α are generated by electroweak gauge boson self-energy diagrams, vertex diagrams, and box diagrams. For the sake of completeness, we also present expressions for the counterterms of order α discussed in Sections 2.6 and 2.7.
The fermion wavefunction counterterms, calculated to all orders in in Section 2.7, are and δZ Due to their excessive length, we abbreviate the integral coefficients of the electroweak gauge boson self-energies and associated counterterms as r i and present explicit expressions for them in the ancillary file twopointcoeffs.m attached to the arXiv submission of this paper. In the notation of Section 2.6, we havē ZZ (s) = r 5 I 1 + r 6 I 2 + r 7 I 4 + r 8 I 5 + r 9 I 6 , (A. 13) andΣ (1,0) W + W − (s) = r 10 I 1 + r 11 I 4 + r 12 I 5 + r 13 I 8 + r 14 I 9 + r 15 I 10 (A.14) for the transverse parts of the bare electroweak gauge boson self-energies of relative order α. Inserting explicit expressions for our master integrals, we find complete agreement through to O 0 with the results of [58,59]  for the relative order α electroweak gauge boson mass and wavefunction counterterms appearing in our calculations. As it was not possible for us to find Eqs. (A.11)-(A.20) in the literature verbatim, it may be that we present explicit all-orders-in-results for the above quantities for the first time.
terms of Eqs. (5.28) and (5.30) of [58] modulo minor typos. 20 Finally, for HVBM's γ 5 scheme, we let B (1,0) VV = Q 2 q Q 2 (a 1 I 13 + a 2 I 17 + a 3 I 18 + a 4 I 19 + a 5 I 20 ) + Q q Q a q a (a 6 I 4 + a 7 I 17 + a 8 I 18 + a 9 I 21 + a 10 I 22 + a 11 I 23 + a 12 I 24 + a 13 I 25 )   where, due to their excessive length, the a i , b i , c i , d i ,ā i ,b i ,c i , andd i are provided in the ancillary files HVBMboxcoeffs.m and Kreimerboxcoeffs.m attached to the arXiv submission of this paper. In fact, it is possible to compare our results for the one-loop box diagrams to the original calculation of [58], where a photon mass λ was employed to regulate infrared singularities, by making the identification ln λ 2 /s → 1/ . Inserting explicit expressions for our master integrals into Eqs. (A.29)-(A.36), we find, modulo typos, complete agreement through to O 0 with Eqs. (B.10) and (B.12) -(B.14) of [58]. To find agreement with [58], we first set the factor of α/(2πs) in Eq. (B.10) to one, the factors of α/(2π) in Eq. (B.12) to 2/s, the factors of α/(2π) in Eq. (B.13) to 2/ s − M 2 Z , and the factors of α/(2π) in Eq. (B.14) to 2/s. In addition, we find that the term inside the brackets on the second line of Eq. (B.12) should read that the term inside the parentheses on the right-hand side of the second equality of Eq.
(B.13) should read 21 and that the factor of x 2 + x 1 in the denominator on the third line of Eq. (B.14) should be replaced by x 2 − x 1 . From the ingredients presented above, we distill complete results for C