Mixed EW-QCD two-loop amplitudes for qq¯→ℓ+ℓ−\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q\overline{q}\to {\mathrm{\ell}}^{+}{\mathrm{\ell}}^{-} $$\end{document} and γ5 scheme independence of multi-loop corrections

We perform a dedicated study of the qq¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q\overline{q} $$\end{document}-initiated two-loop electroweak-QCD Drell-Yan scattering amplitude in dimensional regularization schemes for vanishing light quark and lepton masses. For the relative order α and α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 γ5 scheme and Kreimer’s anticommuting γ5 scheme.


JHEP05(2021)213 1 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 higher-order corrections to the Drell-Yan scattering processes under control. The total cross section in pure Quantum Chromodynamics (QCD) was successfully calculated through to next-tonext-to-next-to-leading order in the strong coupling constant for both the neutral-and charged-current processes [2][3][4][5][6][7][8]. 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) [9,10], but full electroweak (EW) corrections are known only at next-to-leading order [11][12][13][14][15][16][17][18][19]. 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 [20][21][22][23][24][25][26][27], 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 [28][29][30] 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 [31]. 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 most complicated unknown parts of the order α 2 α s two-loop helicity amplitudes for qq → + − . Here, we do not include contributions involving the top-quark, light fermionic loops, or bb initial states, and focus on the region of phase space with center-of-mass energy above all two-particle thresholds.
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 [32][33][34][35][36] 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 [37] 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 [32][33][34]38] was employed and one where Kreimer's γ 5 scheme [39][40][41][42] 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 JHEP05(2021)213 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 section 3.1, we define the scattering process we consider more precisely. In section 3.2, supplemented by appendix A, we describe our general approach to the form factor decomposition of the amplitude, and, in sections 3.3 and 3.4, we give further technical details of the numerator algebra scripts we used to carry out our Feynman diagram calculations in the different γ 5 schemes. 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 sections 4.1-4.3, supplemented by appendices B.1 and B.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 give details specifically relevant to the two-loop mixed EW-QCD calculation. In section 6.1, we provide a recipe to pass from our form factors to scattering amplitudes for states of definite helicity 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 [38,[43][44][45][46] 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.
As has been clear from its inception [38,44], it is technically challenging to handle chiral couplings in dimensional regularization. This is primarily because there is no canonical way to continue the usual four-dimensional definition of γ 5 , into d dimensions. It was emphasized in [38] by 't Hooft and Veltman that one cannot expect to maintain the usual anticommutation relation {γ µ , γ 5 } = 0 (2. 7) and, at the same time, work with traces of Dirac matrices which enjoy all of the usual mathematical properties. In particular, a cyclic Dirac trace was explicitly assumed alongside eq. (2.7) in the algebraic proof that the trace of γ µ γ ν γ ρ γ σ γ 5 must either vanish in dimensional regularization or be subject to a modified set of algebraic rules (see e.g. [48] for a review). 't Hooft and Veltman influenced most later authors to modify the γ 5 anticommutation relation, but, as was pointed out by Kreimer [39], retaining eq. (2.7) while abandoning the cyclicity of the Dirac trace is another path towards mathematical consistency. It would at this stage be reasonable to expect a concrete prescription for the d-dimensional continuation of classical trace formulae involving γ 5 such as tr {γ µ γ ν γ ρ γ σ γ 5 } = −4i ε µνρσ (2.8) in some refined dimensional regularization scheme. In fact, there is no d-dimensional replacement for γ 5 or the Levi-Civita tensor, ε µνρσ , in the dimensional regularization schemes we consider; rather, one simply accepts that eq. (2.6), eq. (2.8), and the familiar properties γ 2 5 = 1 and γ † 5 = γ 5 (2.9)

JHEP05(2021)213
remain unchanged in d dimensions. That γ 5 and ε µνρσ are intrinsically four-dimensional objects inevitably forces one to distinguish a four-dimensional subspace and sacrifice somewhat for what concerns the construction of a manifestly Lorentz-covariant formalism. In this work, it is frequently convenient to explicitly distinguish the four-dimensional and (−2 )-dimensional subspaces of d-dimensional Lorentz indices in perturbative calculations. Following [48], we write k µ =k µ +k µ (2.10) and γ µ =γ µ +γ µ , (2.11) where bars denote four-dimensional objects and hats denote (−2 )-dimensional objects. Similarly, we write the usual four-dimensional identities,
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 [49][50][51] have been worked out at one loop [52] and at two loops [53]. For the convenience of the reader, we provide details of a modern formulation of these developments in section 2.4 below.
As reviewed in [54], many authors have in the past adopted the so-called naive γ 5 scheme of [35] (see also [55]) 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

JHEP05(2021)213
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.
HVBM's γ 5 scheme [32][33][34]38] and Kreimer's γ 5 scheme [39][40][41][42] are the two variants of dimensional regularization we work with in this paper. The precise treatment of γ 5 in these schemes has a profound impact on both the Feynman rules and the handling of Dirac traces which we discuss at greater length in sections 2.2 and 2.3. 1 Indeed, fourdimensional Dirac trace identities which rely on γ 5 manipulations must be, depending on the γ 5 scheme, either carefully reformulated or abandoned. Consider, for example, the algebraic proof given in [47] that the trace of an odd number of Dirac matrices must be zero. For n odd, we have which implies tr {γ ν 1 · · · γ νn } = 0 as desired. As should by now be clear, the above proof is valid only in four spacetime dimensions; different lines of reasoning must be adopted in d = 4 − 2 dimensions. As we shall see, in HVBM's γ 5 scheme, eq. (2.5) facilitates a simple proof of the proposition by induction, whereas, in Kreimer's γ 5 scheme, the trace of an odd number of Dirac matrices is argued to be zero in d dimensions by appealing to a mathematical property of the infinite-dimensional Dirac algebra [39].
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.

JHEP05(2021)213
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: HVBM's γ 5 scheme allows for a very familiar treatment of the usual families of Dirac traces. As pointed out already in [38], the trace of an even number of Dirac matrices with d dimensional indices may be conveniently evaluated using the recursive formula Eq. (2.25) is derived in a manner very similar to that of eq. (2.5), but by repeatedly applying eq. (2.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 for even n may be evaluated in much the same way by replacing γ 5 by the covariant version of its definition, and then evaluating the resulting trace of n + 4 Dirac matrices using eqs. (2.25) and (2.26). As expected, the trace of an odd number of Dirac matrices is still zero in d dimensions. As mentioned above, this can be proven by induction on the trace of eq. (2.5). Let us consider the argument for n = 1. On the one hand, by trace cyclicity and eq. (2.4), we have (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), Finally, the difference of eqs. (2.31) and (2.33) implies that tr {γ ν 1 · · · γ νn } = 0 (odd n) (2.34) by induction. Similarly, eq. (2.28) and the line of reasoning used to establish eq. (2.34) together reveal that 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 [32][33][34] and reviewed in [48], 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 [56] allows us to bypass many finite counterterm computations for fermion self-energies and vertex form factors involving chiral couplings. Larin pointed JHEP05(2021)213 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 [57] 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 [39] 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 higher-order 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).

JHEP05(2021)213
As alluded to above, the standard families of odd Dirac traces still vanish identically when tr is replaced by Tr, and Tr {γ ν 1 · · · γ νn γ 5 } = 0 . (2.37) For the standard families of even Dirac traces, [40] defines where σ is the set of permutations of the 2m indices which appear in the sums on the right-hand sides (m = n and m = n + 2 in (2.38) and (2.39), respectively), subject to the restrictions  5 . As a corollary, due to eq. (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 [42] for a supplementary discussion). Despite the fact that Kreimer's γ 5 scheme was introduced a long time ago [39][40][41][42], 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. [54] 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 [58]. 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 JHEP05(2021)213 to ordinary inverse propagators and "µ-term" insertions of the formk i ·k j , as introduced long ago by Bern and Morgan [52]. 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 [53], other approaches such as that of [59] are certainly possible. It seems that some computational effort is simply unavoidable, as both approaches require the calculation of auxiliary integration by parts reductions [60][61][62]. 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 [29,31], facilitated by passing to a finite basis of master integrals [64,65] through Tarasov dimension shifts [51].
As has been clear at least since the appearance of [52], 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 [52] 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 [51]. 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 [63] 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

JHEP05(2021)213
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 [52], 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 [53] 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 JHEP05(2021)213 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 of k 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 [29,31] which may be directly rewritten using the technology of [53]: where we have introduced the notation D 7 ≡ (k 2 − p 1 − p 2 ) 2 from [29,31] 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

JHEP05(2021)213
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 where we have set (p 1 +p 2 ) 2 = −1. 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.

JHEP05(2021)213
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 section 3.1, we neglect the gauge-invariant subset of Feynman diagrams with closed fermion loops and any contributions involving top quarks throughout this work. 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 twoloop 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 [20] 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 [66] in the standard on-shell scheme where the charge renormalization

JHEP05(2021)213
constant is fixed order-by-order in the Thomsen limit (see e.g. [67] for a review). In the MS scheme it is given by, the pole part of eq. (5.40) of [66] 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 [66], upon which the useful resources [54,67,68] 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:

JHEP05(2021)213
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 [66], 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 B.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

JHEP05(2021)213
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 [66], 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 [66] for the interactions of the electroweak gauge bosons with matter. The Z interaction is parametrized by flavor-dependent axial vector and vector coupling coefficients, w . (2.74) 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 4 Throughout this section, we suppress the color indices of the quarks to streamline our discussion and allow for a unified description of propagator corrections at one loop. The two-loop mixed EW-QCD scattering amplitudes of primary interest to us in this work always come with a color factor of CF , due to the fact that all contributing Feynman diagrams feature exactly one gluon exchange between quark lines. 5 Here, the q 2 → 0 limit should be taken before expanding in . 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 (2.75) and our Kreimer reading point prescription instructs us to begin reading our Dirac traces from the projector insertion. Therefore, using the Feynman rules of [66] 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

JHEP05(2021)213
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 [63]; the result is

JHEP05(2021)213
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 [66] 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 [66], 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

JHEP05(2021)213
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 replacement of the form 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 [69], 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 JHEP05(2021)213 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 [56] 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 gluon-exchange correction to the Zqq vertex shown in figure 4 [48], 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)

JHEP05(2021)213
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 [56], 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 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 [48,56] and, fortunately, the only one relevant to this work which we cannot bypass.

JHEP05(2021)213
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 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. [70]. Thus,Ā

JHEP05(2021)213
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.108) From eq. (2.42), we have (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 [64] and their implementation in Reduze 2 [71][72][73] 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. [74]), 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) in complete agreement with eq. (9) of [56]. 9 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 [56] 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 9 In making this comparison, note that eq. (3) of [56] 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 [56].

JHEP05(2021)213
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 [56] 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 inclusion of the gaugeinvariant contributions proportional to the number of light quark or lepton flavors and all top quark mass corrections to future work. The most important two-loop corrections we omit are those which involve a closed fermion loop; these contributions arise due to electroweak gauge boson self-energy insertions and are already known, see e.g. [26,75]. We also do not consider bb initial states, which lead to double W exchange diagrams involving the top quark; in other words, we restrict ourselves to four light flavors in the initial state.
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 [31], 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 [76] 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 additional directions at the outset. This construction is described in detail in appendix A, see also [77][78][79][80][81] for related work. 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 1 2 (1 ± γ 5 ) u(p i ). 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,

JHEP05(2021)213
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 where the dots denote additional structures which decouple in the limit d → 4 and can effectively be ignored in our projector based setup. The relevant Lorentz structures for our calculation are 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 P C iA DY , for C = C VV , C AA , C VA , C AV , (3.19) to hold in d dimensions 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 JHEP05(2021)213 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 [82][83][84][85][86][87][88] 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 [20] 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 [56], 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 [52] 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, to shorten all Dirac traces at an early stage to the maximal extent possible. 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, 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 [48]; axial vector couplings were not considered in [52].

JHEP05(2021)213
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 [64] 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 [89,90]. 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 [66]. Suppressing the color indices, we have iĀ DY = i C VVTVV +C AATAA +C VATVA +C AVTAV + . . . , (3.27) where the dots denote additional Lorentz structures which can be ignored in our projector based setup, and the Lorentz structures of relevance are defined in (3.10)-(3.13). The corresponding Lorentz projectors in Kreimer's γ 5 scheme fulfill in d dimensions and are given by with D ≡ (2 − )s 2 + 2t(2s + t).

JHEP05(2021)213
Again, we employ eqs.  (3.35) are useful. 12 Of course, a decomposition completely analogous to (3.33) exists for the lepton vertex form factors,  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 [91] which, in particular, employed Passarino-Veltman reduction [63] 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 crosscheck on our calculations.

Infrared dipole singularity structure and subtraction functions
The infrared singularities of general scattering amplitudes in two-loop QCD can be predicted in a process independent way [53,[92][93][94][95][96]. The dipole formula [82][83][84][85][86][87][88] provides a particularly concise and straightforward recipe for the generation of infrared subtraction terms at two loops. 13 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 or mixed gauge groups are certainly possible and have long been known, see e.g. [20,107]. In particular, it was shown in [20] that with straightforward modifications, the dipole formula can be extended to describe the singularity structure of two-loop mixed QED-QCD Drell-Yan scattering amplitudes. This setup covers also the case of the mixed EW-QCD corrections discussed in this paper.
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 [108]. 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 [20], 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 in many places in the literature, see e.g. [83,84] for a thorough discussion. 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. Rewriting the results to make all JHEP05(2021)213 imaginary parts explicit in the physical kinematic region of interest, we find and S (1,0) Once again, as shown in [20], the two-loop mixed EW-QCD Drell-Yan soft anomalous dimension vanishes identically: We define our massless QCD resummation functions in the framework of [109]. In practice, following [20], 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

50)
and from eqs. (3.10) and (3.11) of [109] after discarding all contributions proportional to the number of light quarks. As with the cusp anomalous dimensions, the one-loop QED results may be obtained by making the replacement C F → Q 2 f in eqs. (3.46)-(3.49). Explicitly, we have (3.54) and G 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 (0,0)

(-)
, are defined to be independent of the coupling constants. We have in the notation of [66] 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: T AA . This minor abuse of notation will allow us to rewrite e.g. eq. (3.60) as At higher perturbative orders, we have

(-)
A DY = 4πα    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.

JHEP05(2021)213
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 [54]. 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 Thus, in light of the fact that the higher-order-in-one-loop hard functions first enter eqs. (3.62) at the level of the −1 poles of the relative order αα s expressions, it is not obvious that an analog of eq. (3.67) continues to hold at higher orders in perturbation theory: While it is generally expected that all terms of virtual and real radiative corrections sensitive to ambiguities in the treatment of γ 5 eventually cancel out of combinations which furnish complete fixed-order perturbative corrections to physical observables, hard partonic scattering functions by themselves are not physical observables. We will test whether (3.69) actually is an equality 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 , .

JHEP05(2021)213
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

JHEP05(2021)213
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 [110] 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 [49,50,[112][113][114][115][116]. The integral definitions given above lead to a -decoupled form for the differential equations [89,90], 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 [29,31], 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
In this section, we give our results for the renormalized one-loop amplitude coefficients in both HVBM's γ 5 scheme (3.14) and Kreimer's γ 5 scheme (3.27). Assembling the components provided in appendices B.1 and B.2, we find 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 [67]. Results in Kreimer's γ 5 scheme may be simply obtained from eqs.  [2]. We observe analytically that the 1/ 2 and 1/ poles cancel in the hard scattering functions, leaving a finite remainder, as expected.
Let us emphasize that, upon substituting the explicit expressions for the master integrals discussed in the previous section (see also appendix B.2), we find

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: .

JHEP05(2021)213
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 -44 -

JHEP05(2021)213
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 [29,31], 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 [29] and the m i integrals from [31] must be multiplied by (5.1) That is to say, the integral measure of our two-loop integrals is exactly: due to our neglect of all contributions proportional to the number of light fermion flavors or involving the top quark. Together with eq. (2.57), eqs. (5.5) imply exceptionally-simple forms for the order αα s γqq and Zqq vertex counterterms: δZ (1,1) A,γqq = −Q q δZ (1,1) A, q , (5.7) δZ (1,1) V,Zqq = v q δZ (1,1) V, q − a q δZ (1,1) A, q , (5.8) and δZ (1,1) A,Zqq = −a q δZ (1,1) V, q + v q δZ (1,1) A, q , (5.9) where δZ (1,1) V, q and δZ (1,1) A, q are, respectively, the vector and axial vector components of the order αα s quark wavefunction counterterm. 16 Eqs. (5.6)-(5.9) are in fact the only nonzero counterterms of order αα s which directly contribute to the renormalization of the scattering amplitudes we calculate.
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 That is to say, 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. As expected, once all relations between crossed-and uncrossed-channel multiple polylogarithms were identified, we observed that the 1/ 4 through 1/ poles cancelled analytically in the two-loop hard scattering functions we derived from eqs. (3.62). 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 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.

JHEP05(2021)213
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 to a unique order α 2 α s hard scattering function at O 0 . Even when written as a linear combination of a minimal number of multiple polylogarithms such that the cancellation of infrared poles becomes manifest, 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 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. [117] 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 (1, sin θ, 0, cos θ), and use an explicit basis of spinor states in the chiral representation [47] u + (p 1 ) = s 1/4 0, 0, 1, 0 , 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 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)). 18 The phase conventions of (6.3) above are in line with the spinor helicity formalism, as in [118] 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. 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 [20]. 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 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 Decoupling of irrelevant Lorentz structures
The helicity amplitudes presented in this article have been calculated using projectors involving only a set of Lorentz structures which are independent in four dimensions. In this appendix we show how the additional Lorentz structures in d dimensions can be cleanly separated in d dimensions such that the limit d → 4 is smooth. Our construction is similar to the discussion in [80], but here we consider also explicit occurrences of γ 5 . We present the argument for the anticommuting γ 5 scheme.
The and decompose our two-loop amplitudes in d dimensions according to either by explicit Passarino-Veltman reduction of the tensor integrals, or using a projector method. In fact, we could consider a matrix with elements to construct projectors using the standard spin summation of conventional dimensional regularization, Both problems can be cured by changing to a different basis of Lorentz structures. Before we discuss this new choice we note that the matrix M above is block diagonal, and the VV, AA block (indices 1 . . . 8) is identical to the AV, VA block (indices 9 . . . 16). It is therefore sufficient to restrict the following discussion to the VV, AA block and restrict indices to 1 . . . 8. We consider a new basis of Lorentz structures T α = β R αβ T β such that and the first two structures are strictly orthogonal to the others in d dimensions, The remaining part of R is obtained according to enters, which, unlike M −1 , is free of 1/ poles. For R we find explicitly

JHEP05(2021)213
In the new basis, the matrix M takes the form where M 6×6 has an overall factor of and therefore vanishes for d → 4. Since in the new basis the first two Lorentz structures are orthogonal to the other structures in d dimensions, the projectors for the first two form factors can be constructed without explicitly considering the other structures. Moreover, from the form of M 6×6 we see that contributions due to the other structures are suppressed by a power of . This motivates to consider only the first two Lorentz structures for the calculation of physical helicity amplitudes, just ignoring the other "irrelevant" Lorentz structures (in the terminology of [80]). As pointed out before, the same construction applies to the AV, VA structures as well, leaving us with a total of four Lorentz structures, matching the number of different helicity configurations.

B.1 Relative order α s results
Corrections to the neutral-current Drell-Yan process of relative order α s are produced exclusively via gluon exchange across the initial quark line. As should be clear from the diagrams of section 4.1, the corresponding order α s vertex form factors, are nearly trivial. In HVBM's γ 5 scheme, the finite counterterms relevant to the renormalization of the order α 2 α s scattering amplitude, and δZ (0,1) are also of order α s (see section 2.8 for more details).

B.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 , (B.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 (B.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 [66,67]  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. (B.11)-(B.20) in the literature verbatim, it may be that we present explicit all-orders-in-results for the above quantities for the first time.

JHEP05(2021)213
for our master integrals, we find complete agreement through to O 0 with the relevant terms of eqs. (5.28) and (5.30) of [66] 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 [66], 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. (B.29)-(B.36), we find, modulo typos, complete agreement through to O 0 with eqs. (B.10) and (B.12)-(B.14) of [66]. To find agreement with [66], 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

JHEP05(2021)213
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 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.