Double-real corrections at O(alpha alpha_s) to single gauge boson production

We consider the O(alpha alpha_s) corrections to single on-shell gauge boson production at hadron colliders. We concentrate on the contribution of all the subprocesses where the gauge boson is accompanied by the emission of two additional real partons and we evaluate the corresponding total cross sections. The latter are divergent quantities, because of soft and collinear emissions, and are expressed as Laurent series in the dimensional regularization parameter. The total cross sections are evaluated by means of reverse unitarity, i.e. expressing the phase-space integrals in terms of two-loop forward box integrals with cuts on the final state particles. The results are reduced to a combination of Master Integrals, which eventually are evaluated in terms of Generalized Polylogarithms. The presence of internal massive lines in the Feynman diagrams, due to the exchange of electroweak gauge bosons, causes the appearance of 14 Master Integrals which were not previously known in the literature and have been evaluated via Differential Equations.


Introduction
The electroweak (EW) production of a pair of leptons, each with large transverse momentum, in hadron-hadron collisions, known as Drell-Yan (DY) process [1], is one of the historical testgrounds of perturbative quantum chromodynamics (QCD). The charged-current (CC) and the neutral-current (NC) processes are relevant not only to put stringent constraints on the proton parton density functions (PDFs), but also to perform high-precision measurements of fundamental EW parameters such as the masses and decay widths of the W and Z bosons or the EW mixing angle. Furthermore, they represent an important background to many new physics searches (for a recent review see Ref. [2]). All these studies require precise calculations of higher-order radiative effects and a corresponding implementation in simulation tools that can be used to analyze the experimental data (for a discussion on the status of simulation codes for DY processes see Ref. [3]).
In specific cases like the weak mixing angle or the W boson mass measurements, with a final precision goal in the 1 − 2 · 10 −4 range, all the elements entering the theoretical predictions need to be scrutinized. For instance, in the W mass case it is necessary to assess the uncertainty due to a still inaccurate representation of non-perturbative QCD effects parameterized in the proton PDFs [4][5][6] or in the models present in the QCD Parton Shower, or stemming from the incomplete knowledge of higher-order perturbative QCD, EW, or mixed QCD×EW contributions [7]. These measurements require an excellent control not only on the absolute normalization of the observables, but also on their shape. In this respect a major role is played by final state QED radiation as well as by the interplay of the latter with QCD corrections. A detailed study of this interplay requires the exact evaluation of the next order of perturbative corrections, namely those of O(αα s ) , which is not available yet.
A kinematical limit where the EW corrections play an important role is the so-called Sudakov regime, when the observables are characterized by values of the kinematical invariants (large invariant/transverse masses or large transverse momenta) much larger than the gauge boson masses, yielding large logarithmic factors. The EW O(α) corrections are responsible for the first large correction of this kind [8][9][10][11][12], but it has been shown [13,14] that also O(α 2 ) terms may still be sizeable. The O(αα s ) corrections represent the first QCD correction to these large EW factors and their explicit evaluation is thus needed to get the predictions in the Sudakov regime under control.
The DY cross sections can be expressed as a double perturbative expansion in the strong and electromagnetic couplings, respectively α s and α, which can be formally written as follows, with all the phase-space factors understood in the definition of the coefficients dσ: dσ = dσ 0 + α dσ α + α 2 dσ α 2 + · · · + α s dσ αs + α 2 s dσ α 2 s + · · · + αα s dσ ααs + . . . (1.1) In Eq. (1.1) we recognize terms purely due to the strong or the EW corrections, and also terms where the mixed combined effect of the two interactions is present. QCD corrections to the total cross section have been computed at next-to-leading-order (NLO) in Ref. [15] and at next-to-next-to-leading-order (NNLO) in Refs. [16,17]. Recently the next-to-nextto-next-to-leading-order (N3LO) corrections to the Higgs production gluon fusion process became available [18,19], allowing in turn the estimate of the N3LO corrections in the soft approximation also for EW gauge boson production [20,21]. The NLO-EW corrections have been computed separately for the CC-DY in Refs. [8,9] and for the NC-DY in Ref. [10]. Preliminary steps towards the evaluation of the full NNLO-EW corrections have been accomplished with the discussion of the renormalization of the full two-loop amplitudes [22][23][24][25].
The evaluation of the differential distributions of the final-state products is available in the codes described in Refs. [26][27][28][29] and in those of Refs. [8][9][10][30][31][32][33] respectively with NNLO-QCD and NLO-EW accuracy for the cross section. The inclusion of subsets of dominant higher-order corrections, going beyond the fixed-order description of Eq. (1.1), has been implemented in many codes that match exact matrix elements with a Parton Shower (PS). Focusing on the strong interactions, Refs. [34,35] provide the matching with (NLO+PS)-QCD accuracy, Refs. [36,37] with (NNLO+PS)-QCD accuracy, and Ref. [38] performs the matching in the framework of effective theories. On the EW side, the consistent matching of fixed-and all-orders effects is performed for instance in Refs. [30,31,39]. The resummation to all orders of terms enhanced by logarithms of the lepton-pair transverse momentum is available with next-to-next-to-leading-logarithm (NNLL) accuracy in the codes of Refs. [40][41][42].
The full set of exact O(αα s ) corrections to the total cross section is not available yet due to difficulties in the evaluation of the relevant virtual and phase-space integrals and only subsets of corrections are available. In Ref. [43] the authors considered the QCD×QED contributions to the production of a lepton pair in the qq channel. The O(αα s ) corrections to the decays of Z and W bosons have been computed respectively in Refs. [44] and [45]. In Ref. [46] the mixed two-loop corrections to the form factors for the production of a Z boson have been presented. Very recently, in Ref. [47] the authors evaluated all the two-loop virtual master integrals contributing to the O(αα s ) partonic processes of production of a l − l + or l − ν pair. Moreover, the Altarelli-Parisi splitting functions have been computed with O(αα s ) accuracy in Ref. [48] thus allowing for a consistent subtraction of all the initial-state collinearly divergent terms. NLO-EW corrections to V +jet and NLO-QCD corrections to V + γ final states have been computed in Refs. [49][50][51][52], including the leptonic decay of the vector boson. These results are based on the matrix elements describing the production of a gauge boson (and its subsequent decay) accompanied by one additional hard parton; they therefore include terms of O(αα s ) , but are divergent in the limit of vanishing vector boson transverse momentum.
The absence of an exact calculation of the O(αα s ) corrections to the DY processes has been partially compensated, in the past, by the use of different approximations: the restriction, for the EW corrections, to the subset of final-state QED corrections allowed the factorized combination of QCD and QED corrections [53][54][55]; an additive recipe for the NNLO-QCD and NLO-EW results has been proposed in Ref. [56]; the combination of NLO-QCD and NLO-EW matrix elements, consistently matched with (QCD+QED)-PS, has been described in Refs. [57][58][59].
A calculation of the O(αα s ) corrections to the DY processes near the resonance region has been performed in Refs. [60][61][62]. The calculation was done in the pole approximation, namely retaining all the leading terms contributing to the W (Z) boson resonance. Among the various contributions that the authors analyze, the non-factorizable terms due to softphoton exchange between the production and decay processes result to be negligible for current phenomenological purposes. The conclusion is, therefore, that the treatment of the process in the resonance region, which effectively decouples the production from the decay processes, is sufficient for the level of accuracy needed by current experiments. In particular, the factorizable contributions due to initial-state QCD with final-state QED corrections (emission of photons from the final state) turn out to be the most phenomeno-logically relevant. A comparison is in progress between these analytical results and the approximation of the mixed QCD×EW effects implemented in the Shower Monte Carlo of Refs. [58,59]. However, in the analysis of Refs. [60][61][62] the double corrections to the initial state are not calculated; they are estimated to be negligible.
In this paper we face the problem of the exact evaluation of the O(αα s ) corrections to the total cross section for the production of an on-shell weak boson (W or Z). The importance of this calculation is two-fold. From one side, an exact calculation can give a solid ground and a quantitative check to the estimation of Refs. [60][61][62]. From the other side, individual pieces of our calculation can be important for guiding and checking other ingredients necessary for the treatment of more exclusive observables, such as the gauge boson rapidity distribution, or for the calculation of the mixed QCD×EW infrared subtraction terms.
The evaluation of the O(αα s ) corrections to the production of an on-shell vector boson from qq initial-state annihilation requires the study of four different subprocesses, with 0, 1, or 2 additional partons (gluon, quark, photon) in the final state. The respective contributions to the total cross section for on-shell gauge boson production are obtained by computing the two-loop virtual corrections to the lowest-order amplitude or by integrating the relevant squared matrix elements over the full phase space of the additional partons. In the latter cases we adopt a technique called reverse unitarity, developed for the evaluation of the total cross section for Higgs production [63][64][65]. The standard phase-space integration is turned into the evaluation of "cut" two-loop integrals, namely with the additional condition that the final state particles fulfill the on-shell relation. Integrals with up to three internal massive lines appear in the calculation; some of them were not previously available in the literature and required a dedicated study. The calculation of the total cross section is done by reducing the dimensionally regularized scalar integrals coming from the squared amplitude to a set of Master Integrals (MIs) via integration-by-parts (IBP) identities [66][67][68][69][70][71][72][73]. The MIs are then computed using the differential equations method [74][75][76][77][78][79][80][81][82][83][84]. Their expressions in terms of Harmonic Polylogarithms (HPLs) [85] and their generalizations [86][87][88][89][90] can be found in Refs. [47,91,92].
In this paper we focus on the evaluation of the double-real contribution to the O(αα s ) corrections to the total cross section for on-shell single gauge boson production. We consider all possible channels involved at this order in perturbation theory. This includes qq-initiated process as well as qg-, qγ-, and γg-initiated processes. Since the W boson is charged, it can emit a photon. As a consequence, we need to consider diagrams in which a massive propagator is present along with the massive cut external particle. While the diagrams relevant for Z production give rise to MIs that were already computed in the literature, those for W production introduce additional MIs that are presented here, to our knowledge, for the first time. The cross sections corresponding to the channels under consideration are expressed as Laurent series of ε = (4 − d)/2, where d is the space-time dimension. The coefficients of the series are given in terms of generalized polylogarithms up to weight 3.
The paper is organized as follows. In section 2 we present the partonic processes under consideration in more detail and we define their cross sections as linear combinations of a limited number of MIs. Moreover, we briefly discuss the prescription of the γ 5 matrix employed in this computation. In section 3 we describe how the MIs are computed. In particular, we focus on the evaluation of the soft limits of the MIs, which are used to fix the boundary conditions of the differential equations. In section 4 we present the analytic expressions of the partonic cross sections of all the relevant processes. In section 5 we draw our conclusions. In appendix A we provide the reader with the analytic expressions of the MIs and with the expressions of the soft limits with exact dependence on the regulator ε. The complete set of cross sections and MIs is also given in an ancillary file that we include in the arXiv submission. According to the collinear factorization theorem, the inclusive total cross section for the production of a single gauge boson in hadron-hadron collisions can be written as where V = W ± , Z, the sum over i and j runs over all partons present in the proton (quark, gluons, photons), f i,h are the proton PDFs for a parton i inside hadron h, and each partonic cross sectionσ tot (ij → V + X) admits a double perturbative expansion as depicted in Eq. (1.1). The lowest-order non-vanishing contribution to inclusive single gauge boson production is due to quark-antiquark annihilation, with a cross section of O(G µ ) (G µ is the Fermi constant). At higher orders, for a subprocess initiated by a given pair of partons, one has to consider the virtual corrections to the lower-order amplitudes as well as the contribution of the radiative processes with additional emitted partons in the final state. The cancellation of the soft infrared divergences occurs after the combination of these different partonic cross sections with the same initial state. For instance, in the case of O(αα s ) corrections to quark-antiquark annihilation we have four, separately divergent contributions: with the superscripts a and b in σ ab ααs representing the correction due to a virtual (V) or real (R) exchange in the EW or in the strong interactions respectively. In Eq. (2.2) the sum is free of soft IR divergences and the inclusion of initial-state collinear subtraction terms makes eventually the result IR finite. Moreover, at a given higher perturbative order, more initial states with different combinations of partons ij have to be considered. Focusing on the O(αα s ) contributions, we need to include processes initiated by qg: initiated by qγ:σ and by gγ:σ In this work we study the partonic subprocesses that contribute at O(αα s ) to the inclusive hadronic cross section for the production of a gauge boson with two additional partons in the final state (double-real corrections), i.e. all the processes labeled byσ RR in Eqs. (2.2)-(2.5): We note that the squared matrix elements of processes (2.7)-(2.9) and (2.11)-(2.13) can be obtained by crossing those of processes (2.6) and (2.10) respectively. However, in the evaluation of their total cross sections new MIs, absent in the first two cases, appear, making a dedicated calculation necessary.

Treatment of γ 5
The squared matrix element of each subprocess, averaged over initial spin polarizations and color and summed over final spin polarizations and color, must be computed in an arbitrary number of dimensions d = 4 − 2ε, in order to include all the finite contributions due to the interplay of the squared amplitude with the divergent phase-space integration treated in dimensional regularization. In this respect, to perform our calculation we need to adopt a prescription for the manipulation of the Dirac matrix γ 5 , as it is not defined in a non-integer number of dimensions. Therefore, in the present work we consider the proposal of Ref. [93], and take γ 5 anticommuting with all the other γ µ matrices in arbitrary d dimensions. Before the evaluation of the traces, the product of Dirac γ matrices is rearranged by shifting all the γ 5 to the utmost right position, using the anticommuting property. We do not rely on the possibility of a cyclic permutation of the matrices inside the trace, because under the assumption of anticommuting γ 5 in all d dimensions the cyclicity property of the trace does not hold. Moreover, since there are four independent momenta in the process, it is possible to saturate all the indices of a Levi-Civita tensor resulting from the computation of the traces and thus yielding non-vanishing factors. These terms containing Levi-Civita tensors are responsible for the two following problems, after evaluation of the traces: the presence of gauge-dependent terms, when the polarization sum is done with an arbitrary gauge vector, and the presence of purely imaginary terms out of a squared matrix element, which should obviously be real-valued. The solution is found, consistently with the prescription of Ref. [93], by promoting also the Schouten identity to be valid in an arbitrary number d of dimensions; all the problematic terms thus exactly cancel.

Definition of the total cross section and Reverse Unitarity
We define the total partonic cross sections of the processes under consideration as: (2.14) where z = M 2 V /ŝ is the ratio between the gauge boson mass squared and the partonic centerof-mass energy squared and we conventionally assign the momentum p 3 to the massive gauge boson. The Reverse Unitarity (RU) technique relies on the remark that, in terms of distributions, the following replacement (Cutkosky rule) holds The phase-space measure of each final state particle can thus be rewritten as the difference of two propagators with opposite prescriptions for their imaginary part (with η an infinitesimal positive real number). The integral over the full phase space of the two additional partons, necessary to compute the total cross section, is transformed into the evaluation of the imaginary part of two-loop integrals with the additional constraint that lines corresponding to the final-state particles are cut, i.e. are on-shell (optical theorem). The calculation of the total cross section of processes (2.6)-(2.13) can therefore be accomplished by means of the techniques developed for the study of virtual corrections. After computing the squared amplitude and applying the Cutkosky rule, the phasespace integral of Eq. (2.14) consists of a very large number of terms. Most of these terms, however, are not independent. By means of algebraic relations, Lorentz (LI) and IBP identities (in our case implemented in the codes Reduze [68,69] and FIRE [70][71][72]), it is possible to simplify the sum of these phase-space integrals and express it as a combination of a limited number of irreducible MIs. For the processes under consideration, the number of the independent MIs that eventually have to be explicitly computed is of O(10). The expression of the total cross section of a given process X can therefore be cast as: where the coefficients c X i are rational functions and are process dependent. The cross section is a combination of MIs I i , whose precise number and expressions depend on the process and on some choices applied in the reduction procedure. In our case, the total partonic cross sections of the processes (2.6)-(2.13) have been expressed in terms of (11,13,13,11,7,9,9,7) MIs respectively, with a total of 30 distinct integrals, of which 16 already known in the literature and 14 new. In Section 3 we discuss the techniques developed to compute the new MIs and in Appendix A we list the explicit expressions of all the new integrals (as well as the expressions of the others for completeness), written in terms of HPLs.
The total cross sections of the processes (2.6)-(2.13) are IR divergent. In dimensional regularization the highest-order singularity can be at most an ε −4 pole due to the simultaneous soft and collinear divergences of both additional partons (e.g. photon and gluon in the qq-initiated process). The rational coefficients and the MIs in Eq. (2.17) depend in a non-trivial way on the regularization parameter ε. The explicit expressions of the cross sections are obtained by expanding both in powers of ε, keeping all the terms of the product that are non-vanishing in the limit ε → 0. The total cross sections can therefore be written as Laurent series in ε: We remark that in the qq-initiated processes, in order to extract the soft singularity z → 1 (thus obtaining the ε −4 pole), the following identity is used: with the so-called "plus" distribution defined as Lastly, we note that the assumption that the final-state W boson is on-shell yields additional IR soft divergences with respect to the off-shell case. The production of a W boson differs with respect to the case of a Z boson because of its electric charge: since a photon can be radiated off each charged leg, in the case of W production the amplitude receives a contribution from additional Feynman diagrams. In the case of quark-antiquark annihilation, the additional Feynman diagrams are those in the last row of Fig. 1. From the point of view of strong interactions, the amplitude for W production can be thus divided into two gauge-invariant subsets: the first two rows of Fig. 1, common to W and Z production, and the last one. The invariance under electromagnetic gauge transformations requires instead the sum of all the diagrams, and it can be checked by writing explicitly the charges of up-type quarks, down-type quarks, and of the W boson, respectively Q u , Q d , and Q W .

Evaluation of the Master Integrals
The MIs necessary to compute the total cross sections of processes (2.6)-(2.13) involve at least one massive line (the EW gauge boson in the final state) and possibly an additional one from those diagrams where a photon is emitted off a W leg. For the processes under consideration, we found a total of 30 MIs, of which 16 with one massive line, and 14 with two massive lines. All the integrals with one massive line were already available in the Figure 1. The Feynman diagrams that contribute to the probability amplitude of the processes qq → Zγg and qq ′ → W γg. Diagrams in the first two rows are common to both processes, whereas the two at the bottom are typical of qq ′ → W γg because the photon only couples to charged particles. The Feynman diagrams contained in this article were drawn with Jaxodraw [94].
literature after the evaluation of the NNLO-QCD corrections to the inclusive Higgs boson production in gluon fusion [63,64]. In order to validate the routines developed for the present calculation, we recomputed them and found complete agreement. The computation of all the necessary MIs has been performed using the differential equations method. The system of equations has been written with the help of the package Reduze Specifically, the MIs contributing to processes (2.6)-(2.8) and (2.10)-(2.13) can be written, as a function of the variable z, in terms of HPLs with weights {0, ±1}. The process (2.9) requires instead an enlargement of the basis of functions and the use of non-linear weights, while keeping the variable z. In this case, following Ref. [89], we use the set of weights {0, ±1, − 1 4 , − r 0 4 }. In some cases, HPLs with non-linear weights can be transformed into combinations of HPLs with linear weights at the price of introducing new weights ("letters") in the set ("alphabet") (see e.g. Ref. [90]). In our case, by performing the change of variable ξ = ξ(z) defined through the equations and by introducing two new linear weights a 1 , a 2 defined as HPLs of variable z that contain the non-linear weight − r 0 4 can be expressed in terms of HPLs of variable ξ and linear weights {0, ±1, a 1 , a 2 }. In particular, the additional weights a 1 , a 2 have to be introduced only for those HPLs of variable z that simultaneously contain the weights 1 and − r 0 4 . We further observe that HPLs with the latter combination of weights exactly cancel in the final results for the partonic cross sections of the process (2.9). Two explicit examples of the aforementioned transformations are: In the ancillary Mathematica file we list all the transformations needed for the MIs that contribute to process (2.9). The advantage of this type of transformations is that the HPLs that appear in the final expressions can be easily converted into ordinary logarithms and polylogarithms and evaluated numerically.

Evaluation of the soft limits
We use the soft limit (i.e. z → 1 limit) of the MIs as boundary conditions to the differential equations. We compute the soft limit of all the MIs relevant for the present calculation with the method described in Ref. [65]. The main idea of this method is to rescale the momenta of the final-state partons in the propagators of the MIs by a factor (1 − z) and to perform an expansion of the integrals around the threshold z = 1. The coefficients of these expansions are integrals with simpler structures, e.g. eikonal propagators. By means of the IBP identities it is then possible to reduce these "soft" integrals to a combination of a very small number of "soft Master Integrals", which have to be computed explicitly. By construction, the first term in the threshold expansion of the MIs is the leading behavior as z → 1, i.e. their soft limit. For the processes under consideration, we found that the soft limits of all the necessary MIs can be expressed as combinations of three soft MIs, two of which were already known in the literature while one, to our knowledge, was not available yet. We also observe that in the case of I 21 (z; ε), according to the indexing of Appendix A, the integration constants of the differential equation can be fixed only by computing also the next-to-leading term in the threshold expansion of the soft limit of the MI.
The first soft MI is the pure phase-space integral. It can be computed using the energy-angles parameterization of Refs. [64,96] and reads where we defined the normalization factor common to all MIs . (3.5) The second soft MI appears in the soft limits of some of the MIs relevant for the qγ-and qg-initiated subprocesses. Its expression has been discussed in Refs. [64,65] and reads 1 The third soft MI is peculiar of W production. In this case, the presence of an additional internal massive line spoils the factorization of the different integrations over the energy/angles variables discussed in Refs. [64,96]. More specifically, The integration contour can be chosen such that all the poles of the Γ(a + u) are located to the left of the vertical line defined by u 0 and all the poles of the Γ(b − u) are located to the right. The integration can then be solved using the residue theorem by choosing a finite closed rectangular contour to the left of the vertical line at u = u 0 and then taking the limit of an infinitely extended contour. In this limit, the contribution of the additional lines vanishes and the result of the integral is thus given by the infinite sum of the residues of the integrand. Explicitly, we find: (3.10) In Appendix A.2 we collect the expressions, exact in ε, of the soft limits of all the MIs appearing in this calculation. 1 The expressions Y(z; ε) and X18 in Eq. (3.4) of Ref. [65] differ by a normalization factor, namely -11 - We now present the analytic expressions of the total partonic cross sections of the processes under consideration. For each subprocess we indicate which MIs contribute to the partonic cross section according to the indexing proposed in Appendix A. We present the results expressed as Laurent series in the dimensional regulator ε = (4 − d)/2 and as functions of the dimensionless variable z = M 2 V /ŝ with V = W ± , Z. We remark that all the cross sections are expressed in terms of HPLs up to weight 3, as the coefficients that contain HPLs of weight 4 in the expansion in ε of the individual MIs do not contribute to the cross sections up to O(ε 0 ). The only exceptions to this are the cross sections for the qqinitiated processes, where HPLs of weight 4 coming from integrals I 1 and I 2 (according to the indexing of Appendix A) enter the cross section at O(ε 0 ), but eventually exactly cancel among each other. We remark that some of the MIs contributing to the process gγ → W ± q iqj are represented in terms of generalized HPLs with non-linear weights, as it is explicitly shown in the results of Appendix A. Nevertheless, we observe that the generalized HPLs that are eventually part of the cross section can all be transformed into HPLs with linear weights {0, ±1} and variable ξ(z) defined in Eq. (3.1).
Finally, in order to facilitate the numerical evaluation of the results, we convert all the HPLs appearing in the cross sections into ordinary logarithms and polylogarithms of variables z and ξ(z), and Riemann zeta functions. We perform this conversion with the HPL package [97] and with in-house Mathematica routines.
Throughout this section we use the following normalization factor for the cross sections: . (4.1)

W production
In the calculation of the different subprocesses that contribute to W boson production we have retained the full dependence on the quarks and W boson electric charges, obtaining expressions that are lengthier than those for Z production. For the sake of brevity, we present here results where the explicit charge values have been inserted in the formulae, while in the ancillary files we deliver the generic expressions.

The subprocess q iqj → W ± g γ
We present here the fully inclusive partonic cross section for the tree-level processes: The cross sectionsσ RR q iqj →W ± gγ (z; ε) are obtained by summing the following combination of MIsσ RR q iqj →W ± gγ (z; ε) = 11 k=1 c q iqj →W ± gγ k (z; ε)I k (z; ε) , (4.3) where the explicit expressions of the MIs I k (k = 1, . . . , 11) can be found in Appendix A. After expanding in ε and introducing plus distributions, we recast the results aŝ which is the coefficient of δ(1 − z) in the Born cross section, with sin 2 θ W the squared sinus of the weak mixing angle and N c the number of colors. We remark that the total cross sections for the processes q iqj → W + gγ and q jqi → W − gγ are identical. In terms of ordinary logarithms and polylogarithms, the functions P (n) W ± (z) read: (4.10) We present here the fully inclusive partonic cross section for the tree-level processes: The cross sectionsσ RR q i g→W ± q j γ (z; ε) are obtained by summing the following combination of MIsσ with the explicit expressions of the MIs I k (k = 1, . . . , 23) collected in Appendix A. We observe that c q i g→W ± q j γ k = 0 for k = 3, 7−11, 14, 19−21. The cross section expressed as a Laurent series in the dimensional regulator ε has the form: where A W has been defined in Eq. (4.1) andσ 0 q iqj →W ± (z) in Eq. (4.5). We remark that the cross sections of the subprocesses initiated by a gluon and an up-or a down-type quark differ because of the different electric charge flow probed by the final state photon. For the specific process u g → W + d γ, the functions Q (n) ug,W + (z) read 2 : . (4.17)

The subprocess q
We now focus on the partonic cross section for the tree-level processes: We obtain the cross section aŝ where c q i γ→W ± q j g k = 0 for k = 3, 5 − 11. As a Laurent series, the cross section can be rewritten aŝ with A W andσ 0 q iqj →W ± (z) as earlier defined. We remark that the cross sections of the subprocesses initiated by a photon and an up-or a down-type quark differ because of the different electric charge flow probed by the photon. For the specific process u γ → W + d g, the functions G (n) uγ,W + (z) read 3 : Li 3 (z) (1 − z) 2 − 4z 2 (Li 2 (−z) + ln (z) ln (1 + z)) − 1 9 38z 2 + 142z + 193 Li 2 (z) ln (1 − z) + 147z 2 − 61z 2 + 87 8 Finally, we consider the tree-level processes: The partonic cross sections are written aŝ with c gγ→W ± q iqj k = 0 for k = 3, 5−8, 10−23. Expressed as Laurent series, the cross sections have the form: We remark that the total cross sections of the processes g γ → W + dū and g γ → W − ud are identical. The functions T (1 + 2z) 2 10ζ 2 − 10Li 2 (z) + 13 2 ln 2 (z) Li 2 (z) ln (z) where H − r 0 4 , − r 0 4 , 0; 1 explicitly reads (4.31) 4.2 Z production 4.2.1 The subprocess q iqi → Z g γ We present here the partonic cross section for the tree-level process: The cross sectionσ RR q iqi →Zgγ (z; ε) is obtained by summing the following combination of MIŝ where c q iqi →Zgγ k = 0 for k = 5−7. We can then rewrite the cross section aŝ where A Z is obtained from Eq. (4.1). We defined which is the coefficient of δ(1 − z) in the Born cross section, with C v,i , C a,i the coefficients of the vector and axial-vector couplings of the Z boson to a fermion i. The functions P (n) Z (z) read: We consider here the tree-level process: The partonic cross section is obtained by summing the following combination of MIŝ with c q i g→Zq i γ k = 0 for k = 3, 5−11, 14. The cross section expressed as a Laurent series in the dimensional regulator ε has the form: with A Z andσ 0 q iqi →Z (z) as earlier defined. The functions Q (n) Z (z) read: We present here the cross section for the tree-level process: The result is obtained by summing the following combination of MIŝ with c q i γ→Zq i g k = 0 for k = 3, 5−11, 14. The cross section can be rewritten aŝ For the functions G Since the Z boson does not couple to the photon, the two subprocesses q i γ → Z q i g and q i g → Z q i γ have the same Feynman diagrams upon exchanging the photon with the gluon. Therefore, the two cross sections are identical apart from a color factor due to the sum over final state color in one case or average over initial state color configurations in the other case.

Conclusions
In this work we presented the analytical calculation of the total cross sections of all the partonic subprocesses that contribute at O(αα s ) to inclusive single on-shell gauge boson production, with two additional partons in the final state (double-real corrections). The results are expressed as Laurent series in the dimensional regularization parameter, contain HPLs up to weight 3, and can be cast in terms of logarithms and ordinary Euler polylogarithmic functions. This calculation required the evaluation of 14 new two-loop cut MIs with two internal massive lines, whose explicit expressions are presented in the Appendices. These results are part of the O(αα s ) corrections to the total cross section for inclusive on-shell single gauge boson production. The complete evaluation of the latter requires the calculation of the two-loop virtual corrections to the lowest-order process for gauge boson production (double-virtual corrections) and of the virtual corrections to the processes with a single emission of an additional real parton (real-virtual corrections).

Acknowledgments
We would like to thank Kirill Melnikov for interesting comments on the evaluation of the Master Integrals. We would like to thank the Galileo Galilei Institute for

A Analytical expressions of the Master Integrals
In this Appendix we present the MIs relevant for the evaluation of the total cross sections of processes (2.6)-(2.13). We recall that these are phase-space integrals with phase-space measure where the last term is the on-shellness delta function of the vector boson.
In the following, we separate the result of each MI I k (z; ε) into soft and hard part (borrowing the notation from Ref. [64]): where the soft part comes first and both terms are expanded in ε. For each MI we present the two expansions truncated at the last order in ε that is relevant for the cross sections (we note that for some MIs the last relevant order is different between soft and hard part). In addition, in the ancillary Mathematica file we present the expansion of each MI truncated at the order that contains at most HPLs of weight 4. The soft part of all the MIs is also available exact to all orders in ε in Eqs. (A.34)-(A.63) 4 . We remark that all MIs are written in terms of HPLs of argument z, linear weights {0, ±1, − 1 4 }, and the non-linear weight − r 0 4 . As discussed in Section 3, these HPLs can be converted into HPLs of argument ξ(z) and linear weights {0, ±1, a 1 , a 2 }, where a 1 and a 2 are defined in Eq. (3.2). In the ancillary file we provide the explicit transformations for all the HPLs involved in this calculation and containing the non-linear weight − r 0 4 . Lastly, in all the expressions below we extract an overall normalization factor N (ε) defined in Eq. (3.5).
A.1 Definitions and results expanded in ε