Resonance-improved parton-shower matching for the Drell-Yan process including electroweak corrections

We use the POWHEG method to perform parton-shower matching for Drell-Yan production of W and Z bosons at the LHC at NLO QCD and NLO electroweak accuracy. In particular, we investigate an improved treatment of the vector-boson resonances within the POWHEG method. We employ an independent implementation of the POWHEG method and compare to earlier results within the POWHEGBOX. On the technical side, we provide the FKS formalism for photon-radiation off fermions within mass regularization.


Introduction
The Drell-Yan process is a standard-candle process at the CERN Large Hadron Collider (LHC). Due to the high statistics for the production of W and Z bosons and the clean signature of their leptonic decays, it allows for precision measurements and a thorough test of our theoretical understanding of hadron-collider physics (see e.g. Refs. [1,2] for recent reviews). For example, the Drell-Yan process at the LHC is used to measure the W-boson mass [3] with unprecedented accuracy and Z-boson production is sensitive to the effective weak mixing angle. Moreover, the Drell-Yan process is an excellent tool to extract information on the partonic content of the proton and is employed in PDF fits. Concerning new physics, Drell-Yan production is an important background for Z' or W' searches. For many other searches, the modeling of weak-boson production at large transverse momentum is essential.
To match the experimental precision, the theory predictions have been improved and refined continuously. Due to its simplicity, the Drell-Yan process has often served as a first application for new theoretical techniques. Concerning higher-order predictions in fixed-order perturbation theory, QCD corrections are known at NLO for a long time [4] and the NNLO corrections [5][6][7] are available in fully differential form and cast into flexible Monte-Carlo tools [8][9][10][11][12]. The NLO electroweak (EW) corrections [13][14][15][16] are also available in public tools [14][15][16][17][18][19][20][21]. NNLO QCD corrections have been combined additively with NLO EW corrections in Ref. [22]. A complete calculation of mixed QCD and EW corrections of O(αα s ) is not available, yet. Approximate schemes to include these contributions have been studied for example in Refs. [23,24]. The dominant O(αα s ) corrections have been recently computed in the pole approximation [25,26]. At the same order of couplings, the EW corrections to W/Z-boson production at high transverse momenta recoiling against a hard jet are available [27][28][29][30][31], as well as the corresponding NNLO QCD corrections [32,33]. To predict the transverse-momentum distribution of the vector-bosons also at low transverse momenta, the inclusive NNLO fixed-order calculations have been combined with resummation results [34].
The topic of this work is the parton-shower matching of the NLO predictions for the Drell-Yan process including EW corrections. While EW corrections are generically expected at the percent level due to the small fine-structure constant, logarithmic enhancements due to photon radiation in peaked distributions and due to weak loop diagrams in the high-energy tails of distributions (so-called Sudakov logarithms) [27,28] are well-known. The Powheg method [36] allows to make these EW corrections available along with QCD corrections at the level of unweighted events. Within the framework of the Powheg-Box [38], several implementations of the electroweak corrections for the neutral-current [42] and the charged-current [43,44] Drell-Yan process already exist. We use an independent implementation of the Powheg idea to scrutinize the previous implementations which keep track of photon radiation in the unweighted events [42,44]. In particular, we present new resonance-improved results following the ideas first presented in Ref. [45]. While the Powheg method is NLO exact, one is sensitive to relatively large, unphysical O(αα s ) effects if the W/Z-boson resonance is not properly treated concerning the generation of real radiation. In particular in the context of the W-mass measurement, a previously observed discrepancy between the Powheg result and other tools including EW corrections [24], which might have led to an overestimate of the theoretical uncertainty, disappears. During the completion of this work, an alternative solution to the observed discrepancy has been presented in Ref. [46]. The two approaches are compared at the end of Section 5.
As an additional feature in our approach, the differential K-factor due to EW corrections, encoded in the PowhegB-function, turns out to differ from 1 only at the level of 1% or below around the W/Z-boson invariant-mass peak. Hence, O(αα s ) corrections tō B can be expected to be negligible.
The resonance-improved parton-shower matching is discussed in Section 2. It is convenient to use FKS subtraction [47] to implement the Powheg method [36]. As a technical addition, in Section 3, we present the results for the photon-emission off massless fermions using FKS subtraction in mass regularization instead of dimensional regularization. Some details of the calculation are described in Section 4, in particular those which differ from earlier implementations [42,44]. The phenomenological results are discussed in Section 5 before we conclude in Section 6. Technical issues are detailed in the Appendices.

Resonance-improved Powheg method
The improved treatment of resonances within the Powheg method has recently been discussed in Ref. [45]. Within the Powheg method, a real-emission phase-space point is attributed to the various possible particle emissions or particle splittings using the so-called S functions. To reflect the soft/collinear enhancement in the real-emission matrix-elements, the S functions are usually based on angles and particle energies in a given reference frame. In the collinear limits, the S functions are unambiguously defined. Their specific form outside the collinear limits is, however, a matter of choice. At NLO, the result is independent of the choice for the S functions. However, within the Powheg framework, the S functions should reflect the underlying physics to minimize unphysical artifacts within the higher-order corrections. In particular, in the presence of resonances, the S functions should include information on the resonance structure of the process. For the Drell-Yan process including EW corrections, it is the objective to distinguish between final-state photon radiation (FSR) and initial-state photon radiation (ISR) using the S functions.
For example, the diagrams for the emission of a hard isolated photon (with momentum k) from the leptonic decay products (with momenta l ± ) of an on-shell Z boson (Q 2 = (l + + l − + k) 2 ∼ M 2 Z ) in Drell-Yan are enhanced compared to the diagrams for an emission from the initial-state quarks for which the Z-boson propagator is off-shell (Q 2 = (l + + l − ) 2 < M 2 Z ). Analogously, a radiative return event, where a photon emission from the initial state leads to an on-shell Z boson (Q 2 = (l + + l − ) 2 ∼ M 2 Z ) has only little contributions from FSR diagrams (Q 2 = (l + + l − + k) 2 > M 2 Z ). While FSR and ISR emission beyond the logarithmically enhanced terms cannot be unambiguously defined in general, reflecting the enhancement of different diagrams due to the resonant propagator within the S functions improves the physical description as we will explicitly show in this work for the Drell-Yan process.
The Drell-Yan process with EW corrections is the simplest resonance process which can be discussed in this context. Here, the final state at Born level is unambiguously connected to the Z-boson resonance. Hence, the problem to separate the cross section into contributions with different intermediate resonances (the so-called resonance histories in Ref. [45]) does not apply. Moreover, the center-of-mass frame of the Born process coincides with the center-of-mass frame of the resonance itself. Therefore, the conservation of the resonant momentum when emitting a bremsstrahlung particle, as discussed in Ref. [45], is trivially fulfilled by the standard Powheg momentum assignment. Nevertheless, the assignment of an emitted photon to ISR or FSR via the S functions, i.e. the assignment of the real resonance history according to Ref. [45], has interesting consequences for the Drell-Yan process. A wrongly assigned real resonance history leads to large corrections to the differential NLO K factorB(Φ n ), which is associated to each Born phase-space point Φ n in the Powheg method. Generating real emission at a given phase-space point using the same S functions, these large corrections are compensated at the NLO level, since the Powheg method has NLO accuracy independent of the choice of the S functions. Hence, as mentioned above, the discussion aims at avoiding unphysical artifacts within the higher-order corrections. In particular, this question is relevant for the mixed O(αα s ) corrections which are important for precision predictions for the Drell-Yan process, see [24] and references therein.
In the following, we discuss our choice for the resonance-aware S functions for the Drell-Yan process. The corresponding phenomenological results are discussed in Section 5. Concerning photon radiation, the neutral-current Drell-Yan process has three singular regions. One region corresponds to photons radiated collinearly to the beam with an associated S function denoted by S 0 . Moreover, there are two regions for photons radiated collinearly to the charged leptons, which we approximate to be massless (see Section 4), denoted by S ± . The default S functions [36] can be chosen as where E γ and E ± denote the energies of the photon and the charged leptons in the partonic center-of-mass frame, θ 0 the angle of the photon and the beam, and θ ± the angle between the photon and one of the charged leptons. In the charged-current case, there is only one final-state region with the obvious changes applied. For the Drell-Yan process, the general formalism introduced in Ref. [45] boils down to introducing the Breit-Wigner factors where Q 2 i is the invariant mass of the W/Z-decay products, i.e. Q 2 FSR = (l + + l − + k) 2 for FSR and Q 2 ISR = (l + + l − ) 2 for ISR. The resonance-aware S functions are obtained by using in the definitions of the resonance-improved S functions S res , i.e.
For example, on-shell FSR kinematics (Q 2 FSR = M 2 Z ) leads to i.e. the ISR S function is suppressed by the off-shellness of the ISR diagram.
To investigate the resonance improvement further, we introduce two alternative resonance-aware S functions for the neutral-current Drell-Yan process. The first one is based on the fermion charges which enter the photon-fermion coupling, i.e. we define where q i denotes the charge of the final-state leptons for P q FSR and the charge of the initial-state quarks for P q ISR . The corresponding S functions using P q i instead of P i are denoted by S q i . As already discussed in Ref. [45], the S functions can also be based on explicit matrix-element information. In particular, for the neutral-current Drell-Yan process, the squared amplitude for photon emission can be separated in a gauge-invariant way into an ISR, an FSR, and an interference contribution according to where M ISR includes only the diagrams with photon emission off the quarks and M FSR includes only the diagrams with photon emission off the charged leptons. Due to this physical decomposition, the FKS method can be applied separately to each piece, i.e. the three terms in (8) are treated as three separate real-emission processes. For the ISR/FSR piece, there is only the ISR/FSR singular regions and the correct resonance assignment is guaranteed. For the interference piece, we use the same S functions as for the full Drell-Yan process before. The corresponding phenomenological results are discussed in Section 5, where we refer to the matrix-element splitting as S M .
For the charged-current Drell-Yan process, the S functions based on (7) and (8) cannot be defined. Since the W boson is charged itself, one cannot assign a charge to the initial or final state and a gauge-invariant decomposition of the matrix element is not available either. Hence, we only use the S functions S def and S res in our phenomenological analysis. Since there is only one charged lepton in the final state, there is only one final-state region, e.g. all terms involving a d − are absent and there is no S def − or S res − for the production of W + bosons. As demonstrated explicitly in the neutral-current case in Section 5, the S functions S res are a working solution to achieve resonance-improved results. For the charged-current Drell-Yan process, also the W boson can radiate photons. For hard photons, either the W-boson propagator before or after the photon emission may be on-shell. In this case, the corresponding phase-space point is associated to ISR or FSR by S res in the sense that ISR/FSR means radiation in the production/decay of an on-shell W boson. Hence, the kinematic imprint of the underlying physics is captured as in Z-boson production. Concerning soft-photon emission, diagrams, where the photon is emitted from the W boson, cannot be associated to FSR or ISR in any way. The resonance-improved S functions S res reflect this fact since they do not affect the distinction of FSR and ISR for soft photons where P ISR and P FSR are of similar size. Since S res gives a decent description of the underlying physics for the neutral-current case, we expect it to also improve the results for the charged-current Drell-Yan process. However, further improvements should be searched for in the future and might be particularly needed for distributions, where the results of S res and S M differ in a non-negligible way for the neutral-current Drell-Yan process (see also Section 5).

FKS subtraction using mass regularization
In contrast to NLO QCD calculations, EW corrections can also be performed using mass regularization instead of dimensional regularization, i.e. fermion masses m f regularize collinear singularities and a photon mass m γ is introduced to regularize soft singularities. There is no advantage or disadvantage using one regularization scheme or the other. It is merely a matter of convenience to consider FKS in mass regularization if the virtual corrections are already available in this regularization scheme.
For Catani-Seymour subtraction [48], the results for mass regularization have been presented in full generality in Ref. [49,50]. For collinear-safe observables, the FKS subtraction using mass regularization can be formulated such that only the soft-virtual contribution has to be modified. Here, we present results for photon emission off massless fermions, i.e. we keep the fermion masses only as a regulator in mass-singular logarithms and neglect the fermion masses otherwise.
In complete analogy to the soft-virtual contribution in Section 2.4.2. of Ref. [36], the result for EW corrections in mass regularization reads where α is the fine-structure constant and i, j run over all initial-state and final-state particles. The full virtual one-loop contribution is denoted by V 1−loop . It is assumed to be calculated in mass regularization and, thus, includes mass singular logarithms. Moreover, with the Born matrix element M and following Ref. [49] where q i is the electric charge of particle i and σ i = +1 for incoming fermions and outgoing anti-fermions and σ i = −1 for incoming anti-fermions and outgoing fermions. The functions I ij and Q read and In turn, we define where the index k refers to final-state particles, l to initial-state particles and i again to both. Moreover, s denotes the center-of-mass energy squared, E i the energy of particle i in the partonic center-of-mass frame, p i its momentum, and µ F the factorization scale. The functions I ij are derived from the eikonal integrals calculated in Ref. [51]. The result for the usual plus-distributions regularizing the real-emission contributions in FKS are obtained for the choice ξ c = 1 and δ O = 2. Other values of ξ c and δ O refer to the modified plus-distributions defined in Ref. [36] (see Appendix A). Since ξ c and δ O only influence the treatment of non-singular regions, the dependence of these parameters is the same as in dimensional regularization. The integration of the real-emission contribution over the real phase space and the initial-state collinear remnants which are regularized via plus-distributions are unchanged with respect to the results in Ref. [36] up to the trivial replacement of coupling constants and color factors by electric charges [44] (see also Appendix A). Hence, the results for V using dimensional regularization or mass regularization have to be identical for each phase-space point. We have checked that this is indeed the case for the Drell-Yan process. More details concerning the calculation are given in Appendix A.
To start from an implementation in dimensional regularization, it is convenient to express the above results in terms of the difference of the two regularization schemes. This difference can be written in compact form and reads where Q 2 in the notation of Ref. [36] plays the role of the parameter in dimensional regularization usually called µ 2 . Since only the soft-virtual contribution changes, this result can also be interpreted as a translation rule for the translation of the corresponding virtual amplitude V 1−loop in mass regularization into the finite part of the virtual amplitude V fin in dimensional regularization. In particular, the difference of the virtual amplitudes is of course independent of the subtraction scheme. The same result has already been derived in a different context in Appendix A of Ref. [52].

Calculational setup and input parameters
We employ a process-specific C ++ implementation 1 of the Powheg method [36] which is independent of the Powheg-Box [38]. The calculational setup and the generalization to include electroweak corrections are similar to the calculation presented in Refs. [42,44]. QCD and QED radiation are treated on the same footing concerning their generation.
The QED processes simply enter as additional radiation processes. For the calculation of theB functions, we additively combine QCD and EW corrections (see also Section 5). In the following, we present some details of our implementation, in particular the differences with respect to Refs. [42,44]. As discussed in Section 3, we have employed both mass regularization as well as dimensional regularization to treat the infrared (soft and collinear) divergences of the EW loop integrals. Adding the corresponding soft endpoint of the real-emission contribution, the soft-virtual contribution V is independent of the regularization scheme at each Born phase-space point.
In contrast to the Powheg-Box implementation of the EW corrections to the Drell-Yan process [42,44], we consider all fermions to be massless (except for the regulator masses in mass regularization). In particular, the fermions are treated as massless in the real-emission matrix elements and for the generation of phase space. The massless fermion approximation is adequate as long as the transverse momentum of emitted photons with respect to their emitters is large compared to the physical lepton mass. In the Powheg method, the generation of the hardest emission is usually performed only above a given minimal transverse momentum k min T . If massive fermions are used, this cut-off can in principle be removed for the generation of photon radiation. For a massless calculation, it has to be much larger than the lepton mass. Since a simultaneous simulation of QCD and QED radiation requires a cut-off above the QCD scale anyway, we do not consider this limitation as a serious drawback. We use (k min T ) 2 = 0.8 GeV 2 for the results shown in Section 5. Photon radiation beyond the hardest emission or below k min T is then generated by the QED shower employing the physical lepton masses. To provide the QED shower with massive leptons we use a simple on-shell projection for the lepton momenta. The on-shell projection is performed in the rest frame of the lepton pair, where the combined three momentum of the two leptons vanishes. Keeping the energy of the lepton-pair fixed, we rescale the three momenta of the two leptons by the same factor so that on-shell massive four-momenta are obtained. By construction, three-momentum conservations is respected since the total three momentum of the two leptons still vanishes. Since the combined energy is not changed, also four-momentum conservation is respected. By varying k min T and simulating photon emission only (i.e. QCD emission is turned off), we have verified that the massless fermion approximation is valid within the statistical fluctuations at the per mill level. Hence, in contrast to Refs. [42,44], the FKS construction of the real-emission phase space for massive particles and its usage throughout the calculation is not needed. Moreover, when calculating theB function, the mass-logarithms in the virtual corrections are analytically canceled according to Section 3, while the cancellation has to be achieved numerically in the framework employing massive fermions. For simplicity, all results are given for so-called bare muons, i.e. no photonlepton recombination is performed (see for example [24]) on resolved photon-emissions generated by the Powheg method or the parton-shower.
There is a singular region according to the FKS prescription for initial-state radiation and for the emission off each final-state lepton. Before resonance improvement, we have used the S functions in Section 2 which are also used in the Powheg-Box version 2 with the setting "olddij=1". We do not observe qualitative differences if the default S functions of version 2 are used. We have also adapted exactly the same running for α s like the Powheg-Box, extracted from the Powheg-Box source. To scrutinize the resonance description, we also investigate the various S functions introduced in Section 2.
The virtual corrections for the charged-current Drell-Yan process are taken from Ref. [53]. The virtual corrections for the neutral-current Drell-Yan process have been extensively checked against Ref. [16]. The real-emission matrix elements have been calculated with MadGraph5_aMC@NLO [54]. For simplicity, we do not include photoninduced processes in our analysis. Their inclusion is straight-forward but does not add anything to the discussion of resonance improvements, since the associated photon splitting into a quark pair is a pure initial-state effect. For complete predictions at the percent-level, the photon-induced processes are, however, a non-negligible ingredient.
We generate unweighted events in the LHE format [55] and pass them to Pythia 8.215 [56,57] for showering. All relevant Pythia and Powheg-Box settings are given in Appendix B, where we also discuss how the actual matching can be performed in Pythia. The plots shown in Section 5 are produced with Rivet 2.4.0 [58]. Since we are interested in the effects of the parton-shower matching, we switch off hadronization and multi-parton interactions. In particular, we have checked that the differences in the results due to the choice of S functions are unaffected from hadronization and multi-parton interactions.
We generate inclusive event samples for the LHC running at a center-of-mass energy √ s = 13 TeV. According to the FKS prescription, the radiation phase space for theB function is built starting from tree-level phase-space points. To obtain a finite crosssection in the neutral-current Drell-Yan process, we use a cut on the invariant dilepton mass m ll > 50 GeV (19) applied to the Born phase space, i.e. radiated photons are treated completely inclusively and infrared safety is guaranteed. From the inclusive samples, we select events on the basis of typical identification cuts for the charged leptons. To be specific, we use cuts on the charged lepton transverse momentum p l T and rapidity η l p l T > 20 GeV and η l < 2.5 .
For the charged-current Drell-Yan process, a missing transverse momentum cut is applied which is equivalent to a cut on the transverse momentum of the neutrino at parton level. For the neutral-current Drell-Yan process, we also employ the m ll > 50 GeV cut in the event selection. As mentioned above, we present results for bare leptons in the event selection. Since the parton shower is based on massive fermions, the fermion-mass logarithms which are present in the results in this case are appropriately taken into account.
Concerning the electroweak input-parameter scheme, we employ the so-called G µ scheme. The numerical value for the fine structure constant α Gµ is based on the Fermi constant and the gauge-boson masses in this scheme according to We use α Gµ in the calculation ofB to absorb various higher-order corrections into the leading-order result (see for example the discussion in Ref. [53]). Moreover, the corresponding charge-renormalization constant in this scheme includes the radiative corrections to muon decay and it is independent of the light fermion masses. Concerning the radiation generation, we use the fine-structure constant α(0) defined in the Thomson limit. This setup is also supported by the Powheg-Box implementation Refs. [42,44]. To consistently describe the gauge-boson resonance, we use the complex-mass scheme [59,60], where the gauge-boson masses and all related quantities like the weak mixing angle are treated as complex quantities, where Γ V is the gauge-boson width. As numerical input for the results in Section 5, we use Within the Powheg framework, the quark masses m c and m b only enter as thresholds in the α s -running. The lepton masses only enter in the on-shell projection described above. The CKM-matrix elements V ij are applied as global factors in the Born matrix elements for the charged-current Drell-Yan process. The CKM matrix is taken as block diagonal, i.e. mixing with the third generation is neglected. We use the NNPDF23_nlo_as_0118_qed PDF set (ID: 244600) [61][62][63] with the LHAPDF(6.1.6) [64] interface. Concerning the factorization with respect to EW corrections, we assume the DIS scheme [65] for the given PDF set. The factorization scale µ F and the renormalization scale µ R are set to the W/Z-boson mass M W/Z , i.e. µ F = µ R = M W/Z .

Phenomenological results
For the specified setup and input parameters, we obtain 1890 ± 1 pb for the inclusive neutral-current Drell-Yan cross section (no-cuts applied other than m ll > 50 GeV at Born level) and 10941 ± 2 pb for the inclusive charged-current Drell-Yan cross section (only a technical cut m ll > 1 GeV applied as in Ref. [24]). The error of these benchmark numbers only includes the statistical error due to Monte-Carlo integration of the cross sections. The Powheg-Box results [42,44] agree at the per mill level. All differential results discussed in the following are normalized to our total cross section results.
In Figure 1, we compare our implementation of the Powheg method for the neutralcurrent Drell-Yan process with the Powheg-Box implementation [42]. For the dilepton invariant-mass spectrum m ll without lepton-photon recombination (bare muons), we find good overall agreement with respect to QCD as well as EW corrections, where also the Pythia shower contains either only QCD or only QED radiation. We also show the fixed-order NLO QCD and EW correction to the m ll distribution. For all the absolute predictions shown in the following, there are QCD uncertainties due to missing higher-order corrections (typically assessed via scale variations), PDF uncertainties or non-perturbative uncertainties due to intrinsic transverse momentum of the partons within the proton, hadronization or the underlying event. In Figure 1, we show, for example, the uncertainty due to the scale variation between µ F = µ R = 2M Z and µ F = µ R = M Z /2. As can be seen, QCD uncertainties are rather flat and affect the overall normalization. In contrast, the shape of the leptonic invariant-mass or transverse mass distributions close to their peak is affected much more by EW corrections. Concerning EW corrections, QCD uncertainties for the relative corrections are negligible and not included in Figure 1. While QCD uncertainties have been extensively analyzed for a thorough understanding of the Drell-Yan process and in particular for the measurement of the W-boson mass (see for example Ref. [3,46] and references therein), they are not relevant for the following comparison of shape differences in the distributions due to different resonance improvements and are not discussed further.
In Figure 2, we show the impact of the EW corrections onB as a function of the dilepton invariant mass m ll . For the standard choice of the S function without resonance improvement, the difference between B andB has a pronounced percent-level structure around the resonance peak and is larger than 10% in the low invariant-mass tail. The results agree well with the Powheg-Box prediction. SinceB depends on the choice of the S functions, it has no physical meaning on the level of a fixed-order NLO calculation, where the differences inB are exactly canceled by including real radiation. However,B impacts the results within the Powheg method at the level of higher-order corrections beyond NLO which are included in the prediction. The relatively large EW corrections inB can be attributed to the resonance issue discussed in Section 2. Using S def , the wrongly assigned resonance history is responsible for moving events from the resonance peak to the tails. FSR events which are erroneously assigned as ISR events are moved below the resonance peak. ISR events which are erroneously assigned as FSR events are moved above the resonance peak. Since FSR is dominant due to the larger lepton charges, more events from the peak are moved to lower values of m ll .
Indeed, definingB by distinguishing between ISR and FSR on the basis of the matrix element (labeled by S M in Figure 2), the EW corrections inB are small and resemble the weak virtual corrections (see Figure 12 in Ref. [16]) which form a gauge-invariant subset of corrections. They grow at low dilepton invariant masses since the photon-exchange at LO starts to impact the distribution for which the G µ scheme is not the optimal choice. Subtracting weak corrections, the remaining photonic corrections turn out to be at the per mill level only.
Using S res (see Section 2) to define resonance-improved S functions, the results agree with the matrix-element method at the one-percent level. Including the relevant fermion charges in the definition of S q , the matrix-element based result S M is reproduced within statistical fluctuations.
At the level ofB, NLO QCD and EW corrections are simply added in our calculation. A factorized approach, where the EW corrections are multiplicatively applied to the best QCD prediction, would lead to essentially the same result around the resonance peak, since the EW corrections are so small. The resonance-improved (physical) definition of the underlying S functions seems to imply very small O(α s α) corrections since already the O(α) corrections toB are tiny around the resonance peak. The bulk of the mixed O(α s α) is then contained in the interplay of photon radiation off the final-state leptons with the QCD corrections as found in Ref. [26]. Higher-order EW corrections at the level ofB are certainly negligible around the Z resonance and will not exceed a few per mill even in the low invariant-mass tail (see also Ref. [16]). Figures 3 and 4 finally quantify the impact of the choice of S function on the final Powheg result including QCD and EW effects inB as well as for the generation of radiation. Pythia8 is used as a parton shower. The agreement with the Powheg-Box is at the per mill level at the resonance and differences up to 1% are observed in the tail of the distribution. As expected, the differences due to the choice of S functions are reduced in the full result compared to the results forB. Part of the difference is compensated by generating hardest emissions according to the Powheg method, where real photons are emitted based on the same S function. However, there is a percent-level distortion of the resonance peak as shown in Figure 4, where the relative differences with respect to the default Powheg result are displayed for the different S functions discussed in Section 2. Since the hardest emission generated by Powheg is mostly of QCD type, not all photonic emissions, needed to cancel theB distortion up to O(α 2 ), are generated by the Powheg method. They are added by the QED parton shower instead. The higher-order O(α s α) differences observed in Figure 4 are beyond the formal accuracy of the calculation and could serve as an error estimate of the Powheg method. However, since the physical reason for the deviations can be associated to the resonance treatment, the error would be clearly overestimated if a proper resonance-improved S function is used. The dominant error of the full prediction is not intrinsic to the Powheg part of the computation but depends on the accuracy of the parton-shower tool which is used. For Pythia, the difference between the matched calculation and a LO plus parton-shower calculation for the pure EW corrections only exceeds 5% of the EW corrections at invariant masses below 70 GeV. Close to the Z peak, it is approximately 1% of the LO result. Around 80 GeV, where the EW corrections are largest, it reaches 5%. Since, due to hard QCD radiation, not all first photon emissions are guaranteed to be matrix-element improved in our approach, a part of this uncertainty is also present in the full prediction. It can thus be estimated to be several per mill around the Z peak and can reach the percent level in the tails of the distribution. Our resonance improvement can be combined with a first photon-emission which always has matrix-element accuracy [46], as discussed below. In this case, the uncertainty should be further reduced to approximately 1-2% of the EW correction for a given observable.
For the charged-current Drell-Yan process, the transverse-mass distribution m T is one of the most crucial observables for the W-boson mass measurement. In contrast to the neutral-current process, the matrix element for W-boson production cannot be separated in a gauge-invariant way into ISR and FSR parts. Also, due to the W-boson charge, there is no unambiguous definition of the fermion charges in S q for W-boson production. Hence, for the charged-current case, we only use S res to improve the resonance description. In the left plot of Figure 5, for comparison, we show the transverse-mass distribution for the neutral-current process. The distribution based on S res reproduces the best result, based on S M , within the statistical uncertainty up to m T = M Z . Like for the invariant-mass distribution, there is a percent-level discrepancy for m T > M Z . In analogy, the theoretical uncertainty due to the resonance modeling in W-boson production, which is shown in the right plot of Figure 5, can be assumed to be quite small for m T < M W . Above M W , the relative difference between the results using S def and S res is probably a good estimate for the theoretical uncertainty of the prediction according to the observations in the POWHEG QCDxEW + Pythia8 S M S res S q BOX Figure 4: Same as the lower plot in Figure 3, but the relative differences are shown with respect to the result obtained with S def . We show the results for the different resonanceimproved S functions S res and S q along with S M and the result of the Powheg-Box.
neutral-current case. Hence, a refined resonance improvement could further improve the prediction for the charged-current Drell-Yan process for m T > M W , e.g. by combining our method with the findings of Ref. [46] as discussed below. Since there is only one radiating charged lepton in the charged-current case, the size of the resonance improvement is slightly reduced compared to the neutral-current case. At the Jacobian peak of the transverse-mass distribution, it is about 3 per mill. This difference can be compared to the results in Figure 34 of Ref. [24], where the results of the (not resonance-improved) Powheg-Box matched to Photos [66] as a QED parton shower are compared to the results of the standard QCD only Powheg implementation with Photos. The difference is attributed to higher-order O(α s α) corrections missing in the latter approach. However, the observed difference is at least partly due to the missing resonance improvement as shown in Figure 5. Hence, the true higher-order O(α s α) corrections can be estimated to be much smaller.
These finding have also been addressed in Ref. [46], where a different resonance improvement has been employed which is briefly discussed in the following. If FSR is associated wrongly to ISR due to the resonance insensitive S functions S def , the corresponding energy loss of the leptons is encoded already in theB function, e.g. events are moved from the Z-boson peak to lower invariant masses as shown in Figure 2. This part of the radiation is also not added as FSR in the Powheg generation of radiation, i.e. the leptons do not loose energy again due to photon emission, and NLO accuracy is reached. If, however, the first radiation is dominated by QCD emissions, the parton shower takes over and may add this FSR radiation again such that there is double counting of the energy loss of the leptons. In Ref. [46], as proposed already in Ref. [45], an extension of the Powheg method is employed which allows to generate the first emission in a resonance decay independently of the emissions in the rest of the process. Hence, the first FSR emission in Drell-Yan is always generated by Powheg and the double counting for the leptonic observables is avoided. On the other hand, the corresponding photon radiation which should be generated by Powheg ISR will still be missing since ISR is dominated by QCD radiation. However, precision measurements are based on the leptonic observables, so this drawback is phenomenologically not severe.
As an advantage, the first FSR emission is always described with Powheg precision in contrast to our approach. On the other hand, theB function reflects the physics much better if a resonance-improved S functions is used such that unwanted artifacts in the higher-order corrections can be avoided from the start. In particular, our approach points towards negligible O(α s α) corrections at the level of theB function. Finally, we want to stress that the two solutions to achieve resonance improvement can be easily combined to achieve all the advantages of each method and to avoid all the disadvantages of the other. In particular, the resonance-improved S function discussed in detail in this work can be easily incorporated in the public Powheg version presented in Ref. [46].

Conclusions
We have used a resonance-improved implementation of the Powheg method to investigate the NLO QCD and EW corrections to the Drell-Yan process matched to a parton shower. The importance of a proper treatment of the resonance with respect to associating the emission of photons to ISR or FSR by means of FKS S functions is highlighted. Without resonance improvement, higher-order O(α s α) corrections, which are beyond the formal accuracy of the Powheg method but which are artifacts of the method, enter the prediction. For the Z-boson line shape, they are at the percent level in the resonance region. For the W-boson transverse-mass distribution, they only reach several per mill. However, they are particularly important for the W-boson mass measurement at the level of a few MeV.
A proper resonance improvement of the Powheg method also leads to small NLO EW corrections in the PowhegB function. Hence, the mixed O(α s α) corrections in B can also be expected to be rather small. The dominant O(α s α) effects are then contained in the kinematic effects due to photon emission interleaved with QCD ISR which are particularly well described by the matrix-element improved parton shower within the Powheg method. This picture is consistent with the findings of Ref. [26], where the dominant O(α s α) corrections in the pole approximation are found to arise from EW corrections in the final state combined with QCD corrections in the initial state. The resonance description can be further improved by using the resonance-improved S functions within the approach discussed in Ref. [46].
On the technical side, we provide the results to use the Powheg method for calculations including massless fermions using mass regularization instead of dimensional regularization.
In the future, resonance-improved predictions for the production of weak gauge bosons at high transverse momenta recoiling against a jet are a natural next step to further improve the theoretical description of weak boson production at the LHC.

A Derivation of the soft-virtual contribution for FKS subtraction using mass regularization
In this appendix we derive the soft-virtual contribution to the cross section as given in Equation (9) in Section 3. The partonic part of the real cross section is given bŷ with where M n+1 is the real-emission matrix element with massive fermions and a photon mass is introduced in potentially singular propagators. The phase space is given by where k is the momentum of an emitted photon with a regulator mass m γ , E = m 2 γ + k 2 , k 1,2 are the initial state momenta, and p i the final-state momenta of the particles present at Born level with corresponding energies E i .
We are concerned with the massless limit, i.e. we extract the mass singular logarithms in the real cross section and neglect the masses otherwise. Concerning soft singularities, we subtract and add the soft or eikonal approximation for the matrix element where i, j runs over all external particles, M n is the Born matrix element, and q i , σ i are defined in Section 3. In contrast to dimensional regularization, in mass regularization also the terms with i = j contribute (k i · k i = m 2 i ). Neglecting the soft-photon momentum in the momentum-conservation δ-function, the soft approximation can be integrated analytically [51] and yieldŝ with the definitions given in Section 3. To obtain this result for the massless limit, the results for arbitrary masses in Ref. [51] are expanded, keeping only the mass-singular logarithms and mass-independent terms. The real-emission contribution can now be written aŝ where ξ = 2E/ √ s denotes the relative photon energy and we have defined the plus- for 0 < ξ c ≤ 1. Note, that the second term in (30) does not contain a soft singularity due to the plus-distribution such that the photon mass is no longer needed as a regulator. Concerning the collinear singularities, according to the FKS prescription, we introduce S functions to decompose the real-emission phase space such that eachσ i contains either only the initial-state collinear singularities (i = 0) or one final-state collinear singularity for photon emission collinear to particle i.
To extract the initial-state singularity, one subtracts and adds the collinear approximations for the matrix element |M (initial,i) n+1 where the splitting function is given by P f f (z) = 1+z 2 1−z and z = 1 − E E i = 1 − ξ denotes the momentum fraction of the initial-state fermion i after the emission of the collinear photon. The argument of M n (zk i ) is given in order to state that the reduced momentum zk i enters the Born matrix element. The subtracted cross section is free of collinear singularities, so that the regulator masses can be taken to zero. Hence, the subtracted cross section can be expressed in terms of plus-distributions in complete analogy to dimensional regularization. In the subtraction term, the collinear limit has to be used in the momentum-conservation δ-function. Hence, the angular integrals over the collinear approximation can be again analytically solved to obtain where s = (zk 1 + k 2 ) 2 is the center-of-mass energy squared of the Born matrix element and the integral over the angle θ between the emitter and the photon has been restricted by Θ(δ I − 1 + cos θ) with 0 < δ I ≤ 2. The choice of δ I corresponds to using the plus-distribution [36] 1 −1 with y = cos θ to define the subtracted cross section. For the corresponding hadronic cross section one finds where x 1 = zx 1 , the parton luminosity L(x 1 , x 2 , µ F ) = f (x 1 , µ F )f (x 2 , µ F ) is given by the product of the corresponding PDFs, and s is again given by the center-of-mass energy squared of the Born matrix element M n , i.e. s = (k 1 + k 2 ) 2 . The logarithm of the fermion mass is absorbed by the renormalization of the partondistribution functions (PDFs). In mass regularization, the PDF renormalization for a quark q or the corresponding anti-quark amounts to [16,65] where the usual plus-distribution (ξ c = 1) enters. Moreover, µ F is the factorization scale and K f f = 0 in the MS scheme and K f f = K DIS f f in the DIS scheme with After PDF renormalization, we find σ (initial,1) where is equivalent to the result in dimensional regularization and the second line of (40) contributes to Q initial . For the second initial-state fermion, the equivalent result is found. Concerning the electroweak corrections, the DIS scheme is usually assumed to be the best scheme choice for the current PDF sets.
Concerning the extraction of final-state collinear singularities, we follow the same strategy. In each singular region due to the emission of a photon from a final-state fermion, we subtract and add the collinear limit of the matrix element in mass regularization where z = E i E i +E = ξ i ξ i +ξ . Again the subtracted cross section is free of singularities such that the regulator masses can be taken to zero. The subtracted cross section can be written in terms of the same plus-distributions as the result in dimensional regularization. Since the collinear limit has to be used in the momentum-conservation δ-function of the subtraction term, one can analytically integrate the angular integrals to obtain where δ O plays the same role as δ I in restricting the range of the angular integration. After rescaling the momentum p i to absorb the explicit factor of z in the delta-function and the matrix element, solving the ξ-integral yieldŝ corresponding to Q final defined in Section 3.

B Pythia and Powheg-Box setup
We use the default settings for Pythia 8 if not specified otherwise in the following. Concerning the matching of the Powheg results with the Pythia shower, Powheg and Pythia are both based on transverse-momentum as a shower-evolution scale. However, the precise definition of the evolution scale differs. Ignoring this difference, i.e. assuming that the Pythia and Powheg scales are the same, the simplest matching solution is to use the Powheg scale of the hardest radiation (given by SCALUP in the LHE file) as starting scale for Pythia. One can use the Pythia settings to ensure that Pythia uses the starting scale provided by Powheg also when adding radiation to resonance decays. Since there is only one scale for the hardest radiation in our setup, using this flag is equivalent to implementing the Pythia function UserHooks::scaleResonance as recommended in Ref. [67].
To take the mismatch of the precise definition of the transverse-momentum based evolution scales into account, the Pythia 8 plug-in PowHegHook is provided by the Pythia authors [68]. An adapted version of this plug-in is also used in the Powheg-Box POWHEG QCDxEW + Pythia8 Figure 6: The plot shows the relative difference in the prediction of the invariant-mass distribution using the simple matching, where the Powheg scale is the starting scale of the shower, with respect to the vetoed shower as described in the text. Here, we have employed S M , however, the result does not depend on the S functions in a significant way.
implementations Ref. [42,44] and we follow the usage employed there. If this plug-in is used, Pythia starts its evolution at the kinematic limit and vetoes all radiation with a Powheg evolution scale above the scale of the hardest Powheg emission. This procedure is possible because the Pythia and Powheg scales are similar. A truncated shower [69] is not needed in this case according to Ref. [36,68]. The default PowHegHook plug-in does not apply any veto in resonance decays. Therefore, one has to modify the code to enable vetoing in resonance decays. Moreover, the correct Powheg scale has to be calculated by the plug-in. In our case, the ISR scale is the transverse momentum and the FSR scale is based on the soft-collinear limit of the transverse momentum, i.e.
where the index i denotes the radiated particle and j denotes the emitter after radiation. The necessary changes in the code can be also found in the public implementation of the neutral-current Drell-Yan process in the Powheg-Box. Our implementation is equivalent.
The results presented in Section 5 employ the matching based on PowHegHook. The simpler matching procedure described before, shows differences at the level of 1% in the tails of the invariant-mass and transverse-mass distributions (see Figure 6). Changes in the ratios shown in Figure 4 and Figure 5 are negligible, since the parton shower matching affects the Powheg events in the same way. Moreover, we use the fine-structure constant in the Thomson limit: TimeShower:alphaEMorder 0 SpaceShower:alphaEMorder 0 Finally, we use an internal NNPDF set of Pythia:

PDF:pSet 15
For the left panel of Figure 1, we use only the Pythia QCD shower. We switch off the QED shower by the following settings: To compare EW effect in the right panel of Figure 1, we disable the QCD shower in Pythia by the following settings: Checks:event off PartonLevel:Remnants off SpaceShower:QCDshower off

TimeShower:QCDshower off
For performance reasons we restrict the analysis to parton level by switching off hadronization and multi-parton interactions: HadronLevel:all off PartonLevel:MPI off We checked that hadronization and multi-parton interaction do not alter the results. In particular, we checked that there are no differences for the ratios of the results employing different (resonant-improved) S functions. Concerning the Powheg-Box, we use the default settings. In addition, we choose the G µ scheme, µ F = µ R = M Z and photon-induced processes are switched off: For the plots shown in Section 5, we use: However, the default choice for the S functions or the treatment of lepton masses does not lead to noticeable differences in the differential distributions. Moreover, we choose the complex-mass scheme for the charged-current Drell-Yan process [42] by setting the corresponding flag complexmasses = .true. in init_couplings.f. If svn version r3337 or earlier is used for the neutral-current Drell-Yan process [44], one should make sure to change nubound to a value of 10 5 or higher.