NLO QCD+EW predictions for HV and HV+jet production including parton-shower effects

We present the first NLO QCD+EW predictions for Higgs boson production in association with a l nu or l+ l- pair plus zero or one jets at the LHC. Fixed-order NLO QCD+EW calculations are combined with a QCD+QED parton shower using the recently developed resonance-aware method in the POWHEG framework. Moreover, applying the improved MINLO technique to H l nu and H l+ l- production at NLO QCD+EW, we obtain predictions that are NLO accurate for observables with both zero or one resolved jet. This approach permits also to capture higher-order effects associated with the interplay of EW corrections and QCD radiation. The behavior of EW corrections is studied for various kinematic distributions, relevant for experimental analyses of Higgsstrahlung processes at the 13 TeV LHC. Exact NLO EW corrections are complemented with approximate analytic formulae that account for the leading and next-to-leading Sudakov logarithms in the high-energy regime. In the tails of transverse-momentum distributions, relevant for analyses in the boosted Higgs regime, the Sudakov approximation works well, and NLO EW effects can largely exceed the ten percent level. Our predictions are based on the POWHEX BOX RES + OpenLoops framework in combination with the Pythia 8.1 parton shower.


Introduction
The discovery of the Higgs boson [1,2] has opened the door to the direct experimental investigation of the Higgs and Yukawa sectors of the Standard Model. While present measurements of Higgs boson properties and interactions are consistent with the Standard Model [3], the full set of data collected during Run II and in subsequent runs of the LHC will provide more and more stringent tests of the mechanism of electroweak symmetry breaking. In this context, the associated production of a Higgs and a vector boson, pp → HV with V = W and Z, plays a prominent role. In spite of the fact that the total cross sections for these so-called Higgsstrahlung processes are subleading as compared to Higgs boson production via gluon fusion and vector-boson fusion, the possibility to reconstruct the full HV final state and the clean signatures that result from leptonically decaying vector bosons offer unique opportunities of testing Higgs boson interactions with vector bosons and heavy quarks (see Refs. [4][5][6] and references therein). The associated HV production makes it possible to disentangle Higgs boson couplings to W and Z bosons from one another and to measure them in a broad kinematic range. In addition, the presence of the associated vector boson allows for an efficient suppression of QCD backgrounds. In particular, pp → HV is the most favorable channel for measurements of the H → bb branching ratio, and thus for determinations of the bottom Yukawa coupling. In HV production with H → bb decay, the boosted region, with Higgs boson transverse momentum above 200 GeV, plays a particularly important role, both in order to achieve an improved control of the QCD backgrounds [4] and for the sensitivity to possible anomalies in the HV V couplings. Higgsstrahlung processes permit also to probe invisible Higgs boson decays, both through direct measurements of pp → HZ with invisible Higgs decays and through indirect bounds based on measurements of the H → bb branching ratio.
The accuracy of present and future measurements of HV production, at the level of both fiducial cross sections and differential distributions, calls for increasingly accurate theoretical predictions. The inclusion of higher-order QCD corrections is crucial, both for total rates and for a precise description of the QCD radiation that accompanies the production of the HV system. The role of QCD corrections can be particularly important in the boosted regime or in the presence of cuts and for observables that are sensitive to QCD radiation.
In general, in order to account for experimental cuts and observables, higher-order QCD and EW predictions should be available for arbitrary differential distributions, and experimental analyses require particle-level Monte Carlo generators where state-of-the-art theoretical calculations are matched to parton showers. Finally, when QCD and EW higherorder effects are both sizable, also their combination needs to be addressed.
Theoretical calculations for the associated-production processes are widely available in the literature. Among the numerous studies on HV production at next-to-leading order (NLO) QCD we quote here Refs. [7][8][9]. Predictions for inclusive HZ and HW production at next-to-next-to-leading order (NNLO) QCD have first been obtained in Refs. [10,11] and are implemented in the VH@NNLO program [12]. Besides contributions of Drell-Yan (DY) type, where the Higgs boson results from an s-channel V * → HV subtopology, Hig-implement EW corrections.
In this paper, for the first time, we present NLO QCD and NLO EW calculations for the production of a Higgs boson in conjunction with a ν or + − leptonic pair, plus zero or one jet, at the LHC. While, for convenience, the above-mentioned processes will often be denoted as HV /HVj production (with V = W ± and Z) in the rest of the paper, all the results we are going to present always correspond to the complete decayed final-state processes, with spin effects, off-shell and non-resonant contributions taken into account.
At NLO QCD we include the full set of O(α S ) contributions to pp → HV and O(α 2 S ) contributions to pp → HVj. Although terms of non-DY type are implemented in our codes, we have not included them in our simulations. In addition, we do not include NNLO-like loop-induced contributions to HZ plus 0 and 1 jet production.
Besides showing fixed-order NLO QCD+EW predictions at parton level for typical observables, we also present full NLO+PS simulations for HV and HVj production. To this end, we have implemented our NLO calculations for HV and HVj production into four separate codes (HW ± , HW ± j, HZ and HZj) in the POWHEG BOX framework. In this way, we have consistently combined the radiation emitted at NLO QCD+EW level with a QCD+QED parton shower. In this context, photon radiation from the charged leptons can lead to severe unphysical distortions of the Z-and W -boson line shapes, if not properly treated. This problem was first pointed out in the context of NLO QCD+PS simulations of off-shell top-quark production and decay, and was solved in the context of the POWHEG BOX framework by means of the so-called resonance-aware method [58]. The first application of this method and its variants, in the context of electroweak corrections, has appeared in Refs. [59,60]. In this paper, we exploit the flexibility of the resonance-aware method to perform a fully consistent NLO QCD+EW matching in the presence of non-trivial EW resonances. To this end, our NLO calculations and generators are implemented in the new version of the POWHEG BOX framework, known as POWHEG BOX RES. In this recent version, the hardest radiation generated by POWHEG preserves the resonance virtualities present at the underlying-Born level. At the same time, the resonance information can be passed on to the parton shower, which in turn preserves the virtualities of intermediate resonances of the hard process in subsequent emissions.
Similarly to what was done in Ref. [33] for HVj production at NLO QCD, we have applied the improved MiNLO [34,35] approach to HVj production in order to get a sample of events that has simultaneously NLO QCD accuracy for HV plus 0 and 1 jet. In the MiNLO framework, also the NLO EW corrections to HV and HVj production have been consistently combined in the same inclusive sample. This can be regarded as an approximate treatment of O(α S α EM ) corrections in observables that are very sensitive to QCD radiation and receive, at the same time, large EW corrections. Moreover, although we do not present a rigorous proof, based on considerations related to unitarity and factorization of soft and collinear QCD radiation, we will argue that our MiNLO predictions should preserve full NLO QCD+EW accuracy in the phase space with zero or one resolved jets. As we will see, this conclusion is supported by our numerical results.
While our NLO EW results are exact (apart from the treatment of photon-initiated contributions), we also present approximate NLO EW predictions in the so-called Sudakov limit, where all kinematic invariants are well above the electroweak scale. Specifically, based on the general results of Refs. [49,61], we provide explicit analytic expressions for all logarithmic EW corrections to pp → HV + 0 and 1 jet in next-to-leading-logarithmic (NLL) approximation. Based on the observed accuracy of the NLL Sudakov formulas, this approximation can be exploited both in order to speed up the evaluation of EW corrections at NLO and in order to predict the dominant EW effects beyond NLO.
All needed matrix elements for pp → HV + 0 and 1 jet at NLO EW have been generated using the OpenLoops program [62,63], which supports the automated generation of NLO QCD+EW scattering amplitudes for Standard Model processes [64][65][66]. The implementation in the POWHEG BOX RES framework was achieved exploiting the generic interface developed in Ref. [67]. For what concerns NLO QCD corrections, on the one hand we implemented in-house analytic expressions for the virtual corrections. On the other hand, following the approach of Ref. [33], for real-emission contributions we used MadGraph4 [68] matrix elements, via the interface described in Ref. [69].
The paper is organized as follows. In Sec. 2 we introduce the various ingredients of HV and HVj production at NLO QCD+EW. In particular, in Sec. 2.2 we present a schematic proof of the NLO QCD+EW accuracy of MiNLO predictions for inclusive observables. Further technical aspects of the calculation as well as input parameters and cuts are specified in Sec. 3. Fixed-order NLO QCD+EW predictions are discussed in Sec. 4, while in Secs. 5 and 6 we present NLO+PS QCD+EW results for HV production and MiNLO QCD+EW results for HVj production, respectively. The predictions of the NLO+PS HV and MiNLO+PS HVj generators are compared in Sec. 7. Our main findings are summarized in Sec. 8. In the appendices we document the validation of EW corrections in HV production against HAWK (App. A), detailed NLO EW formulas in the Sudakov approximation (App. B), a reweighting approach that we employ in order to speed up the evaluation of EW corrections (App. C), and technical aspects of the interface between the POWHEG BOX RES and Pythia 8.1 (App. D).

NLO QCD and EW corrections to HV and HVj production
In this section we describe the QCD and EW NLO corrections to the production of a Higgs boson in association with a ν or + − leptonic pair plus zero or one additional jets. For convenience, these Higgsstrahlung processes will be denoted as associated HV and HVj production, with V = W ± or Z. However, all results presented in this paper correspond to the complete processes including all spin-correlation and off-shell effects. The combination of HW + /HW + j and HW − /HW − j Higgsstrahlung will be denoted as HW /HWj production. In our calculations, we have considered only one leptonic generation, and all leptons are treated as massless.

NLO QCD+EW matrix elements
In this section we describe the various tree and one-loop amplitudes that have been assembled to form a NLO QCD+EW Monte Carlo program based on the POWHEG BOX framework [30]. Figure 1. A sample of QCD (a, b), and EW (c), real radiation diagrams contributing to HW − j production. While only one photon or gluon at a time is present at fixed order, for illustration purpose, in (a) and (c), we have shown various possible gluon and photon emissions.
Associated HV production proceeds through quark-antiquark annihilation at leading order, which corresponds to O(α 3 EM ). In HVj production, where the leading order corresponds to O(α S α 3 EM ), additional (anti)quark-gluon initiated processes contribute. All O(α S α 3 EM ) NLO QCD corrections to HV production have been computed analytically, since they simply affect the V qq vertex, and the calculation of the real and virtual corrections is trivial. In HVj production, the virtual O(α 2 S α 3 EM ) NLO QCD corrections have been computed analytically [70]. The color-and spin-correlated Born amplitudes and the real contributions at O(α 2 S α 3 EM ) have been computed using the automated interface [69] between the POWHEG BOX and MadGraph4 [68]. The real contributions involve tree diagrams with either an additional gluon or an external gluon replaced with a qq-pair. Example diagrams are shown in Fig. 1 (a, b). Figure 2. A sample of virtual EW diagrams contributing to HW − j production.
The virtual EW corrections to HV and HVj production comprise loop amplitudes up to pentagon and hexagon configurations, respectively. Example diagrams for HVj production are shown in Fig. 2. All the internal resonances have been treated in the complex-mass scheme [71,72] throughout.
As pointed out in the Introduction and illustrated in Fig. 2 (c), the virtual NLO EW amplitudes induce a dependence on the Higgs trilinear coupling λ HHH . This dependence arises both from the bare virtual amplitudes and from the Higgs boson self-energy entering the Higgs boson wave-function renormalization. In view of the possibility of exploiting precision measurements of Higgsstrahlung processes for an indirect determination of λ HHH , we allow λ HHH to be set independently of the Higgs boson mass. 1 The real NLO EW corrections to HV and HVj production comprise QED radiation off all charged particles, i.e. they have an additional photon in the final state, as illustrated in Fig. 1 (c). Photon-induced real radiation contributions, where the photon is crossed to the initial state, are, on the other hand, not considered here, as they are suppressed by the small photon density in the proton. These corrections for HV production have been computed for the first time in Ref. [44] and are included in the HAWK [45] Monte Carlo generator. Interestingly, they reach several percent for inclusive HW production, but remain at the 2% level when leptonic selection cuts are applied, and are negligible for HZ production [73]. For HVj production, photon-induced contributions enter already at Born level: however, they are of O(α 4 EM ) and thus formally subleading with respect to the O(α S α 3 EM ) leading order. Still, the NLO QCD corrections to these photon-induced processes are of O(α S α 4 EM ) and thus formally of the same order as the NLO EW corrections to the quark-antiquark and (anti)quark-gluon initiated channels in HVj production. Also not considered here are mixed QCD-EW bremsstrahlung contributions to HVj production at O(α S α 4 EM ). These tree-level contributions are finite and can easily be investigated separately. Similar contributions in the NLO EW corrections to V +jet production are known to yield relevant contributions only in jet observables at very large transverse momentum [65,74]. Finally, also virtual QCD corrections to HV γ production contribute formally at O(α S α 4 EM ) and are thus of the same perturbative order as the NLO EW corrections to HVj production. However, if a photon isolation is applied, as is done in this paper (see Sec. 3.4), HV γ production can be considered as a separate process and thus excluded from the definition of HVj production.
All the electroweak real and virtual corrections have been computed using a recent interface of the POWHEG BOX RES to OpenLoops [67].
In this study we combine NLO QCD and EW corrections in an additive way, i.e. corresponding perturbative contributions are simply added. At fixed order, an improved description can easily be obtained via a factorized ansatz, where differential NLO QCD cross sections are multiplied with relative EW correction factors. Such a multiplicative combination can be motivated from the factorization of soft QCD radiation and EW Sudakov logarithms, which can be tested comparing relative NLO EW corrections for HV and HVj production.

MiNLO approach for HVj production at NLO QCD+EW
In order to obtain an optimal description of QCD radiation, both in the hard and soft regime, all NLO QCD+EW calculations for pp → HVj have been performed using the "Multiscale improved NLO" (MiNLO) [34] method. This approach effectively resums logarithmic singularities of soft and collinear type to NLL accuracy, thereby ensuring a finite HVj cross section in all regions of phase space, even when the extra jet becomes unresolved. In the MiNLO approach, NLL resummation is achieved by means of a CKKW scale setting [75,76] for the strong coupling factors associated with each QCD vertex, together with an appropriate factorization-scale choice and NLL QCD Sudakov form factors. These are applied to all internal and external lines corresponding to the underlying-Born skeleton of each event. In addition, improving the MiNLO resummation as described in Ref. [35], we have obtained a fully inclusive description of HV production with NLO QCD accuracy in all phase space regions. In other words, besides providing HVj kinematic distributions that are NLO accurate and also finite when the hardest jet goes unresolved, the improved MiNLO predictions for pp → HVj are NLO accurate also for distributions in inclusive variables such as the rapidity or the transverse momentum of the HV pair.
All NLO QCD+EW predictions for pp → HVj presented in this paper, both at fixed order and including matching to the parton shower, are based on the MiNLO approach, which is applied to all contributions of NLO QCD and NLO EW type. Technically, the MiNLO Sudakov form factors and scale choices are implemented at the level of pp → HVj underlying-Born events that correspond to so-calledB terms in the POWHEG jargon. 2 Note that the MiNLO procedure resums only logarithms associated with soft and collinear QCD singularities that result form the presence of QCD radiation at Born level, while QED radiation is not present at Born level. Thus there is no need to introduce NLO EW effects in the MiNLO Sudakov form factors. This implies that, in contrast to the case of NLO QCD, the NLO EW corrections to pp → HVj do not need to be matched to the MiNLO form factors. In practice, for what concerns the EW corrections, the MiNLO procedure is applied in a way that is equivalent to Born level.
For observables where QCD radiation is integrated out, the MiNLO improved NLO EW contributions assume the form where Φ HV and Φ j denote the factorized phase spaces of the HV system and the jet, respectively. The termB EW HVj (Φ HV , Φ j ) includes O(α EM ) corrections 3 of virtual and real type, and the latter are integrated over the corresponding emission phase space. The MiNLO approach is implemented through an implicitly understood CKKW scale choice for the α S term inB EW HVj , and through the NLL Sudakov from factor ∆ k T (Φ j ) in Eq. (2.2). For later convenience, together with the Sudakov form factor, we introduce a corresponding emission 2 Real-emission events of NLO QCD and NLO EW type are related to underlying-Born events of type pp → HVj via FKS mappings [77]. 3 Since Born contributions are part of the usual QCDB term, in theB EW term we only include O(αEM) corrections.
kernel K(Φ j ) that is formally related to ∆ via In the following, based on the factorization properties of soft and collinear QCD radiation, encoded in the kernel K(Φ j ), and using the unitarity relation where MiLO denotes the Born (or LO) version of the MiNLO approach. The above identity can be written as Here the singularities associated with QCD radiation in the soft and collinear limits are factorized 4 into the pp → HV Born term times the NLL kernel K(Φ j ), while the B fin HVj remainder is free from singularities. Thus, upon integration over the jet phase space, the B fin HVj remainder yields only O(α S ) suppressed contributions with respect to B HV , while using the unitarity relation (2.4) it is easy to show that the singular term in Eq. (2.9) leads to Eq. (2.8).
Thanks to the fact that applying the MiNLO approach to NLO EW contributions is largely equivalent to applying MiNLO at Born level, the NLO EW accuracy property (2.5) can be proven along the same lines as for the LO accuracy property (2.7). As sole additional ingredient, the NLO EW proof requires certain factorization properties of soft and collinear QCD radiation. More precisely, the factorization properties of Eq. (2.9) must hold also in the presence of EW corrections, i.e.
Here the remainderB EW, fin HVj should be free from QCD singularities, so that it yields only O(α S )-suppressed contributions relative toB EW HV , when the extra jet is integrated out. Based on this natural assumption, in full analogy with the LO case, we easily arrive at 11) which is equivalent to the hypothesis (2.5).
In summary, based on unitarity and factorization properties of QCD radiation, we expect that the improved MiNLO procedure applied to NLO QCD+EW matrix elements for pp → HVj should preserve its full QCD+EW accuracy when the jet is integrated out. As we will see, this conclusion is well supported by our numerical findings in Secs. 4-7. Nevertheless, due to the schematic nature of the presented derivations and related assumptions, the above conclusions should be regarded as an educated guess that deserves further investigation.

Sudakov approximation at NLO EW
In the Sudakov high-energy regime, where all kinematic invariants are of the same order and much larger than the electroweak scale, the NLO EW corrections are dominated by soft and collinear logarithms of Sudakov type. Based on the general results of Refs. [49,61] we have derived analytic expressions for the NLO EW corrections to HV and HVj production in NLL approximation. Details and scope of this approximation are discussed in App. B.
The Sudakov approximation at NLO provides us with qualitative and quantitative insights into the origin of the dominant NLO EW effects. Moreover, it can be easily extended to the two-loop level [51,52], thereby opening the door to approximate NNLO EW predictions based on the combination of exact NLO EW corrections with Sudakov logarithms at two loops. From the practical point of view, the Sudakov approximation at NLO permits to obtain the bulk of the EW virtual corrections at much higher computational speed as compared to an exact NLO EW calculation.
In Sec. 5, we will assess the quality of the Sudakov approximation 5 through a detailed comparison against exact NLO EW corrections. Finally, in App. C, we show how the NLL EW approximation can be used in order to speed up the Monte Carlo integration, while keeping full NLO EW accuracy in the final predictions.

The POWHEG BOX RES framework at NLO QCD+EW
The QCD+EW NLO calculations for HV and HVj production have been matched to parton showers using the POWHEG method. To this end, we used the recently-released version of the POWHEG BOX framework, called POWHEG BOX RES. The major novelty of this new version is the resonance-aware approach [58], which guarantees a consistent treatment of intermediate resonances at NLO+PS level. This is achieved by generating the hardest radiation in a way that preserves the virtuality of resonances present at the underlying-Born level. At the same time, the resonance information can be passed on to the parton shower, which in turn preserves the virtuality of intermediate resonances of the hard process in subsequent emissions. This method was introduced in order to address the combination of NLO QCD corrections with parton showers in the presence of top-quark resonances. However, since it is based only on general properties of resonances and infrared singularities, the resonance-aware approach is applicable also to the combination of EW corrections with QED parton showers. In fact, this method has already been applied in the context of electroweak corrections in Refs. [59,60].
In the POWHEG BOX RES jargon [58], a radiated parton (or photon) can be associated to one or more "resonances" present in the process, or to the "production" part, if it cannot be associated to a particular resonance. The POWHEG BOX RES framework automatically finds all the possible so-called "resonance histories" for a given partonic process. For the processes at hand, considering QED radiation, only two resonance histories are detected: a production history, where the photon can be emitted by any quark (both in the initial and in the final state), and a vector-boson decay history, where the photon is radiated off a final-state charged lepton, and the virtuality of the intermediate vector boson needs to be preserved. Soft photons that are radiated from a W resonance are attributed either to the production subprocess or to the W decay, consistently with the virtualities of the quasi-resonant W propagators "before" and "after" the photon emission.
The treatment of QED radiation was first introduced in the POWHEG BOX for the calculation of the EW corrections to Drell-Yan processes [59,60,78,79]. In this context, leptons were considered as massive particles, and QED subtraction in the POWHEG BOX was implemented accordingly. In the study at hand, leptons are treated as massless, and we have implemented the treatment of photon radiation off massless charged particles (both leptons and quarks). To this end, we have adapted the QCD soft and virtual counterterms already present in the POWHEG BOX to the QED case. Moreover, we have computed a new upper-bounding function for the generation of photon radiation with the highest-bid method, as described in Ref. [29].
By default, in the POWHEG BOX RES framework, only the hardest radiation out of all singular regions is kept, before passing the event to shower Monte Carlo programs like Pythia or Herwig. In this way, for each event, at most one of the decaying resonances (or the production part of the process) includes an NLO-accurate radiation. Moreover, in case of combined QCD and EW corrections, QED emission occurs in competition with the QCD one. The POWHEG BOX RES uses the highest-bid method to decide what kind of radiation (QED or QCD, initial-or final-state) is generated. Due to the larger centerof-mass energy available in the production stage, initial-state radiation is enhanced with respect to final-state radiation, and since the QCD coupling is larger than the QED one, initial-state quarks tend to radiate gluons rather than photons. Thus, QED emission from the decay of a resonance would hardly be kept at the Les Houches event (LHE) level, and the QED radiation would mainly be generated by the shower Monte Carlo program.
The resonance-aware formalism implemented in the POWHEG BOX RES framework offers the opportunity to further improve the POWHEG radiation formula. With this improvement, first introduced in Ref. [80], radiation from each singular region is generated and, instead of keeping only the hardest overall one, the hardest from each resonance is stored. As a result, the LHE file contains a radiated particle for each decaying resonance, plus possibly one emission from the production stage. In this way NLO+LL accuracy is ensured for radiation off each resonance. The subsequent shower from each resonance generated by the Monte Carlo shower program has to be softer than each corresponding POWHEG radiation. 6 All NLO+PS results presented in this paper are based on this multiple-radiation scheme.
As a final remark, we note that in the POWHEG BOX RES framework both the HV and HVj processes can be computed at NLO or NLO+PS level with only QCD corrections, with only EW corrections, or with combined NLO QCD+EW corrections. 7

OpenLoops tree and one-loop amplitudes
All needed amplitudes at NLO EW have been generated with OpenLoops [62,63] and implemented in the POWHEG BOX RES framework through the general interface introduced in Ref. [67]. Thanks to the recursive numerical approach of Ref. [62] combined with the COLLIER tensor reduction library [81], or with CutTools [82], the OpenLoops program permits to achieve high CPU performance and a high degree of numerical stability. The amplitudes employed for the EW corrections in this paper are based on the recently achieved automation of EW corrections in OpenLoops [64][65][66].
Within OpenLoops, ultraviolet and infrared divergences are dimensionally regularized in D dimensions. However, all ingredients of the numerical recursion are handled in four spacetime dimensions. The missing (4 − D)-dimensional contributions, called R 2 rational terms, are universal and can be restored from process-independent effective counterterms [83][84][85]. The implementation of the corresponding Feynman rules for the complete EW Standard Model in OpenLoops is largely based on Refs. [86][87][88][89]. Relevant contributions for HV and HVj production have been validated against independent algebraic results in D = 4 − 2 dimensions. UV divergences at NLO EW are renormalized in the on-shell scheme [90] extended to complex masses [71].

Input parameters, scales choices and other aspects of the setup
In our pp → HV (+jet) simulations at NLO QCD+EW, we have set the gauge-boson masses and widths to the following values [91] M Z = 91.1876 GeV, The latter are obtained from state-of-the-art theoretical calculations. Assigning a finite width to the Higgs boson in the final state would invalidate EW Ward identities: we then consider the Higgs boson as on shell with Γ H = 0 and set its mass to M H = 125 GeV. The top-quark mass and width are set respectively to m t = 172.5 GeV and Γ t = 1.5083 GeV. All other quarks and leptons are treated as massless. In the EW corrections, the top-quark contribution enters only at loop level, the dependence of our results on Γ t is thus completely negligible.
For the treatment of unstable particles we employ the complex-mass scheme [71,72], where finite-width effects are absorbed into complex-valued renormalized masses The electroweak couplings are derived from the gauge-boson masses and the Fermi constant, G µ = 1.16637 × 10 −5 GeV −2 , and the electromagnetic coupling is set accordingly to where µ 2 W and the squared sine of the weak mixing angle The absolute values of the CKM matrix elements are set to  Our default set of parton-distribution functions (PDF) is the NNPDF2.3_as_0119_qed set [92], that includes QED contributions to the parton evolution and a photon density. 9 The value of the strong coupling constant corresponding to this PDF set is α S (M Z ) = 0.119. 8 By default we use the Gµ scheme throughout. However, in the POWHEG BOX RES framework, there is the option to evaluate the virtual EW corrections using αEM computed in the Gµ scheme, and use the Thomson value αEM(0) = 1/137.035999 in the evaluation of the contribution due to photon radiation. 9 It corresponds to the PDF set 244800, in the LHAPDF6 [93] numbering scheme.
Finally, in HV production, the renormalization and factorization scales are set equal to the invariant mass of the HV pair at the underlying-Born level, where 1 and 2 are the final-state leptons, while in pp → HVj the improved MiNLO [34,35] procedure is applied, and the scales are set accordingly.
Predictions at NLO+PS generated with the POWHEG method are combined with the Pythia 8.1 QCD+QED parton shower using the "Monash 2013" tune [94]. Effects due to hadronization, multi-particle interactions and underlying events are not considered in this paper.

Physics objects and cuts in NLO+PS simulations
In the following we specify the definition of physics objects and cuts that are applied in the phenomenological NLO+PS studies presented in Secs. 5-7.
All leptonic observables are computed in terms of dressed leptons, which are constructed by recombining the collinear photon radiation emitted within a cone (in the (y, φ) plane) of radius R γ = 0.1 from charged leptons, and the recombined photons are treated as unresolved particles. Observables that depend on the reconstructed vector bosons are defined by combining the momenta of the dressed charged leptons and the neutrino associated with their decay. The latter is taken at Monte Carlo truth level.
Jets are constructed with FastJet using the anti-k T algorithm [95,96] with R = 0.5. The jet algorithm is applied in a democratic way to QCD partons and non-recombined photons, with the exception of photons that fulfill the isolation criterion of Ref. [97] with a cone of radius R 0 = 0.4 and a maximal hadronic energy fraction h = 0.5. The hardest of such isolated photons is excluded from the jet algorithm and is treated as resolved photon.
The following standard Higgsstrahlung cuts are applied. For every dressed charged lepton we require In HW /HWj production, we also impose where / E T is the transverse momentum of the neutrino that results from the W -boson decay at Monte Carlo truth level. In HZ/HZj production, the invariant mass of the dressedlepton pair is required to satisfy Besides these inclusive selection cuts, we also present more exclusive results in the boosted regime. In this case, we impose the following additional cuts on the transverse momentum of the Higgs and vector bosons Such a selection of events with a boosted Higgs boson improves the signal-over-background ratio in the H → bb decay channel.

Results for HV and HVj production at fixed NLO QCD+EW
In this section we present fixed-order NLO QCD+EW predictions for pp → HV and pp → HVj at 13 TeV. For HVj production the improved MiNLO approach [34,35] is applied. Higgs boson production in association with W and Z bosons is discussed in Secs. 4.1 and 4.2, respectively. Predictions based on exact NLO EW calculations (apart from photoninitiated contributions that have been neglected) are compared against the Sudakov NLL approximation (see App. B), which includes virtual EW logarithms supplemented by an exact treatment of QED radiation. The fixed-order results presented in this section are not subject to the cuts and definitions of Sec. 3.4. No acceptance cut is applied, and differential observables are defined in terms of the momenta of the Higgs and vector bosons. The latter are defined in terms of the momenta of their leptonic decay products at the level of underlying-Born events, i.e. before the emission of NLO radiation. Photons and QCD partons are clustered in a fully democratic way using the anti-k T algorithm with R = 0.5. Effectively this procedure corresponds to an inclusive treatment of QED radiation.
Besides total cross sections, we consider various differential distributions, focusing on regions of high invariant masses and transverse momenta, where EW corrections are enhanced by Sudakov logarithms. Such phase-space regions play an important role for experimental analyses of HV production in the boosted regime.
An in-depth validation of our fixed-order NLO EW results for HV production against the ones implemented in the public Monte Carlo program HAWK [44,45], is presented in App. A.

HW and HWj production
In this section we focus on NLO results for pp → HW and pp → HWj. In Tab. 1 we report inclusive NLO cross sections. In the case of HWj production, the improved MiNLO approach yields finite cross sections without imposing any minimum transverse momentum on the hardest jet. For comparison, we report also HWj MiNLO cross sections for the case where a minimum p T of 20 GeV is required for the hardest jet. In this case, the MiNLO Sudakov form factor plays hardly any role, since it damps the cross section only at p T of the order of a few GeV, i.e. far below the imposed cut. Thus, at fixed order, the MiNLO procedure only affects the choice of scales, as described in Sec. 3.3. The EW corrections lower the inclusive NLO QCD cross section by roughly −7% for HW ± production and −5% for HW ± j production, while they amount to only −2% when a resolved jet with p j 1 T > 20 GeV is required in the HW ± j calculation. Inclusive cross sections in the NLO QCD+NLL EW approximation differ by several percent from the exact NLO QCD+EW results. This is expected, since the NLL approximation is only valid in the high-energy regime.
In the following we investigate the impact of EW corrections and the validity of the NLL approximation in differential distributions for HW − and HW − j production. Results for HW + (j) production (not shown) are very similar.
In Fig. 3 we plot the invariant mass of the reconstructed HW − pair, both for HW − and HW − j production. The three curves represent predictions at NLO QCD, NLO QCD+EW    and in NLO QCD+NLL EW approximation. While EW corrections have a moderate impact on the total cross sections, they affect the tail of the M HW distribution in a substantial way. At large M HW we observe the typical Sudakov behavior, with increasingly large negative EW corrections that reach the level of −25% (−30%) for HV (HVj) production at 2 TeV. The Sudakov NLL approximation captures the bulk of these large EW corrections as expected.
In the tail it agrees at the percent level with the exact result for both processes, while for moderate invariant masses it overestimates EW correction effects by up to 5%.  In Fig. 4 we investigate the transverse momentum of the Higgs boson. Also in this case EW corrections become negative and large in the tail, exceeding −20% in the TeV region. For both processes the Sudakov approximation agrees at the percent level with exact NLO EW results for p H T > 300 GeV. The EW corrections have a sizable impact also on the missing transverse momentum distribution, shown in Fig. 5. Size and shape of these corrections are very similar to the ones observed for the Higgs boson p T distribution.
In Fig. 6 we present HW − j predictions for the distribution in the p T of the leading jet. At low p j 1 T (left plot) the MiNLO Sudakov form factor damps soft and collinear singularities at zero transverse momentum yielding finite cross sections below the Sudakov peak, which is located around 3 GeV. Concerning EW effects, the NLL approximation converges to the exact NLO results already for values of p j 1 T around 200 GeV. In the region of moderate transverse momentum, NLO EW corrections are nearly constant, and in the limit of vanishing jet-p T they converge towards an EW K-factor that is very close to the one of the NLO QCD+EW calculation for the inclusive pp → HW − cross section (see Tab. 1). This observation is consistent with the theoretical considerations presented in Sec. 2.2, namely with the fact that EW corrections are insensitive to soft and collinear QCD radiation, and that MiNLO predictions for HVj production preserve NLO QCD+EW accuracy when the extra jet is integrated out. In fact, in the inclusive distributions of Figs. 3-5, we observe TeV dσ/dp j1 TeV dσ/dp j1 TeV that the EW corrections obtained from HW − and HW − j calculations are very similar, with small differences that can be attributed to NNLO effects.
Finally, in Fig. 7 we see that the EW corrections affect the rapidity of the leading jet in a rather uniform way over the whole phase space. We have observed a similar behavior of EW corrections in several angular distributions.

HZ and HZj production
In line with the discussion of HW and HWj production, we present in this section fixedorder results for HZ and HZj production. In Tab. 2 we collect the inclusive cross sections 13 TeV 13 TeV  In Figs. 8 and 9 we show distributions of the invariant mass of the reconstructed HZ pair and of the transverse momentum of the Higgs boson. Similarly as for HW (j) production, the inclusion of EW corrections is essential in the tails of these distributions, where the NLL Sudakov approximation agrees well with the exact NLO EW predictions.

Results for HV production at NLO+PS QCD+EW
In this section we present NLO QCD+EW predictions for HV production completed by the Pythia 8.1 QCD+QED parton shower using the "Monash 2013" tune [94]. All predictions are subject to the cuts and physics object definitions specified in Sec. 3.4, and NLO EW corrections are treated exactly throughout, except for photon-initiated processes, that have  been neglected. The NLL Sudakov approximation is only used in order to speed up the Monte Carlo integration, as detailed in App. C. In Sec. 5.1 we compare predictions at fixed-order NLO QCD+EW against corresponding predictions at the level of Les Houches events, which include only the hardest emission generated in the POWHEG BOX RES framework, and at NLO+PS level, where the full QCD+QED parton shower is applied. The effect of EW corrections is studied in Sec. 5.2 in the case of fully showered NLO+PS simulations. By default, at NLO+PS level, the full QCD+QED parton shower is applied, both for NLO QCD+EW and for pure NLO QCD simulations. Occasionally, we also present NLO QCD simulations with a pure QCD shower, where QED radiation is switched off. Such predictions are labeled "QCD (no QED shower)".
The consistent combination of the NLO radiation to the parton shower requires the vetoing of shower emissions that are harder than the radiation generated in the POWHEG BOX RES framework. Since no standard interface is available in a multi-radiation scheme, we have implemented a dedicated veto procedure on the Pythia 8.1 showered events, as described in App. D. This veto procedure is applied in case of NLO QCD+EW simulations. Instead, in case of NLO QCD simulations combined with the Pythia 8.1 QCD+QED shower, only QCD radiation is restricted by the POWHEG BOX RES hardest scale, while arbitrarily hard QED radiation can be generated by the shower.
We have verified that inclusive cross sections at NLO+PS QCD and NLO+PS QCD+EW agree within statistical uncertainties with the corresponding fixed-order results reported in Tabs. 1 and 2. Thus, in the following we will focus on differential distributions.

From fixed NLO QCD+EW to NLO+PS QCD+EW
In this section we compare NLO QCD+EW predictions at fixed order with NLO+PS ones at LHE level and completed with the Pythia 8.1 shower. Since the various Higgsstrahlung processes behave in a very similar way, we will focus on HW − production.
In Fig. 10 we plot the rapidity of the reconstructed HW − pair, which is NLO accurate, and its transverse momentum, which is only LO accurate. Due to the inclusiveness of the rapidity of the HW − pair, we find, as expected, very good agreement, within the integration errors, among the three predictions. The fixed-order curve for the transverse momentum displays the typical divergent behavior at low p T . At LHE level, instead, the divergence is tamed by the Sudakov form factor. The effect of the parton shower is modest in the tail of this distribution, while at low p T it slightly shifts the position of the Sudakov peak.    Figure 11. Pseudorapidity (left) and transverse-momentum distribution (right) of the Higgs boson in HW − production. Same curves and labels as in Fig. 10.
In Fig. 11 we plot the pseudorapidity and the transverse momentum of the Higgs boson: thanks to the inclusiveness of this variable, we find again very good agreement among the three predictions.

Impact of the EW corrections in NLO+PS events
In this section we investigate EW correction effects at the level of fully showered NLO+PS predictions. In Fig. 12 we show the rapidity (left) and the transverse momentum (right) of the charged dressed lepton in HW − production. In the rapidity distribution, the impact of NLO EW effects is constant and amounts to about −7%. The shape of the p T distribution, instead, changes drastically due to EW Sudakov logarithms in the high-p T region, where differences with respect to the pure QCD predictions reach −30% around 1 TeV.  In Fig. 13 we plot the transverse mass of the reconstructed W − boson where ∆φ is the azimuthal angle between the charged lepton and the missing transverse momentum. Similarly, as for the lepton rapidity, the EW corrections do not change the shape, but lower the differential cross section by about 7% with respect to the pure QCD corrections. If no QED shower is activated when Pythia 8.1 showers QCD-corrected events, the curve that is obtained is very similar to the QCD one, i.e. the impact of the QED shower is small for this distribution and no radiative tail can be observed.
In Fig. 14 we show the rapidity and the transverse momentum of the Higgs boson in the boosted regime, as defined by the cuts of Eq. (3.10). The EW corrections have a constant negative impact around 10% on the rapidity distribution, and reach −25% around 1 TeV. Similar conclusions can be drawn for the rapidity and transverse momentum of the W − boson.
We conclude this section by presenting kinematic distributions for HZ production in Figs. 15-17. In Fig. 15 we show the distribution in the rapidity and the transverse momentum of the dressed electron. The EW corrections give a constant contribution of about −5% in the plotted rapidity range, while in the high-energy tail of the p T distribution the EW corrections decrease the differential cross section by roughly 30% due to Sudakov logarithms.  In Fig. 16 we plot the invariant mass of the reconstructed leptonic pair in the region around the Z resonance. In spite of the fact that the shape of the Z resonance it known to receive very large O(α EM ) radiative corrections (see e.g. Refs. [79,98]), NLO EW effects turn out to be almost constant and as small as −5% when we compare showered NLO+PS predictions at NLO QCD+EW versus NLO QCD. This is due to the fact that the bulk of the O(α EM ) radiation is correctly described by the QED shower in Pythia 8.1. The importance of O(α EM ) radiation becomes evident when we switch off the QED shower ("no QED shower") in the NLO QCD simulation. This results in a radiative tail with distortions of up to 40% in the region below the Z peak. In HZ production, the momentum of the vector boson can be fully reconstructed. Thus, in Fig. 17 we display the rapidity and transverse-momentum distributions of the Z boson in the boosted regime, as defined by the cuts of Eq. (3.10). These results are very similar to the ones obtained for the Higgs boson in HW − production in Fig. 14. While EW corrections have a constant impact of about −8% on the rapidity distribution, the tail of the p T distribution is dominated by large negative EW Sudakov logarithms, and we observe differences with respect to the pure QCD result of the order of −25% for p T ∼ 1 TeV.

Results for HVj production at NLO QCD+EW with MiNLO+PS
In this section we study pp → HVj at NLO+PS accuracy in the MiNLO approach, denoted in the following as MiNLO+PS. Similarly as in the previous section, in Sec. 6.1 we first compare NLO QCD+EW predictions obtained with MiNLO at fixed order against corresponding results at the LHE level or including also the full QCD+QED parton shower. The effect of EW corrections is studied in Sec. 6.2 in the case of fully showered MiNLO+PS simulations. The cuts and physics object definitions of Sec. 3.4 are applied throughout, and we do not impose any cut that requires the presence of jets.

From fixed-order MiNLO to MiNLO+PS at NLO QCD+EW
In Fig. 18 we analyze the rapidity and the transverse momentum of the reconstructed HW − system. As a result of the MiNLO prescription, the rapidity distribution, as well as any other inclusive observable, is finite. For the rapidity distribution we observe that the three predictions are very close to each other. At variance with Fig. 10, here fixed-order predictions for the p T distribution are finite at small p T , since soft and collinear divergences are suppressed by the MiNLO Sudakov form factors. Moreover, the NLO accuracy in the spectrum of the HW − system leads to an improved agreement between fixed-order and NLO+PS results. In Fig. 19 we show the pseudorapidity and the transverse-momentum distributions of the Higgs boson, finding again very good agreement among the three predictions.
We refrain from presenting results for HW + j and HZj production since they behave qualitatively very similar as the results shown here for HW − j production.

Impact of EW corrections at MiNLO+PS level
The impact of EW corrections at the level of fully showered MiNLO+PS predictions for HW − j production is illustrated in Figs. 20-22. For the distributions in the rapidity and transverse momentum of the Higgs boson in the boosted regime (Fig. 20) we find that the EW corrections induced by the boosted cut, p H T ≥ 200 GeV, are nearly independent of y H and around −10%, while they grow up to −20% and beyond when p H T enters the TeV regime. These results closely agree with the corresponding ones shown in Fig. 14 for the NLO+PS simulation of inclusive HW − production. Consistently with the fixed-order findings discussed in Sec. 4, also this observation supports the theoretical considerations of Sec. 2.2, where we have argued that MiNLO improved predictions for HVj production should preserve NLO QCD+EW accuracy when the extra jet is integrated out. Also other inclusive observables, such as the distribution in the missing transverse momentum shown in Fig. 21, confirm this observation.
The EW corrections to the leading-jet p T distribution, shown in Fig. 22, do not feature the standard Sudakov behavior. In this distribution, EW effects remain rather small, at the level of −5%, in the entire plotted range, i.e. from very low jet-p T up to 400 GeV. This is not surprising, since a similar "non-Sudakov" behavior for inclusive jet spectra was already observed in Ref. [65,74] for the case of V + jets production. Another important feature of Fig. 22 is that EW corrections are nearly constant in the region where the jet p T approaches zero. Again, this confirms the considerations made in Sec. 2.2 regarding the factorization of EW corrections in the presence of soft or collinear QCD radiation, and the NLO QCD+EW accuracy of inclusive MiNLO simulations. To be more precise, in the left panel of Fig. 22 we see that EW corrections effects are nearly constant at small p T with the exception of the first bin. This effect can be attributed to photonic contributions to the jet transverse momentum, and to the fact that the Sudakov peak associated with the damping of QCD radiation is located well above the one associated with the damping of QED radiation. This mismatch tends to enhance the relative importance of QED radiation in the region between the QCD and QED Sudakov peaks. In any case, this effect cancels upon integration over the soft region of the jet spectrum. Thus, it should not spoil the expected NLO QCD+EW accuracy of inclusive MiNLO predictions.  We conclude this section by discussing the impact of NLO EW effects in HZj production, illustrated in Figs. 23 and 24. The distribution in the Z-boson p T in the boosted regime (Fig. 23) features the typical Sudakov EW behavior, with negative EW corrections that exceed −20% in the tail. In the leading-jet p T distribution (Fig. 24)    atively small and rather constant EW corrections. Both distributions behave similarly as the corresponding distributions for HW − j production.
We refrain from showing further plots for HZj or HW + j production, since EW correction effects are quite similar to the ones already discussed.

Comparison between the HV and HVj generators
In this section, we discuss and compare NLO+PS predictions for pp → HV against MiNLO+PS predictions for pp → HVj, both at NLO QCD+EW accuracy. A similar comparison at NLO QCD accuracy was presented in Ref. [33]. Since the various Higgsstrahlung processes behave in a very similar way, we will focus on the associated HW − production.
The comparison between the HV and HVj generators is motivated by the fact that the improved MiNLO prescription [35] applied to HVj production provides NLO accuracy also for inclusive HV quantities, i.e. for observables where the associated jet is not resolved. While this is well known at NLO QCD level, in Sec. 2.2 we have argued that also NLO EW accuracy should be preserved when the jet is integrated out. We also study the dependence of our results on scale variations. To this end we apply standard seven-point variations obtained by multiplying the central value of the renormalization and factorization scales µ R and µ F , defined in Eq. (3.6) for HV production, by the factors K R and K F , respectively, chosen among the seven pairs Scale variations are in general larger in HVj with respect to HV production. This is due to the fact that, in standard POWHEG BOX simulations, the scale associated to the emission of the hardest jet is kept fixed at the corresponding transverse momentum, while scale variations are applied only at the level of the so-calledB term, where QCD and QED radiation are integrated out. For this reason, scale variations in MiNLO+PS simulations provide a more realistic estimate of scale uncertainties associated with QCD radiation.
Figures 25-27 display differential distributions subject to the cuts of Sec. 3.4. Red bands correspond to scale variations for HV and HVj production. We do not show the statistical uncertainties associated to the integration procedure on these bands, since they are much smaller than the width of the bands. When plotting instead the blue curves for the distributions computed at the central scales, we display the statistical uncertainties of the integration procedure as an error bar. The plots on the left-hand side show the uncertainty band for the HV process, while the ones on the right-hand-side show the uncertainty band for HVj production. In Fig. 25 we display the rapidity distribution of the HW − system. Since this inclusive quantity is predicted at NLO QCD accuracy by both simulations, we find very good numerical agreement between the two curves at NLO QCD+EW level. The uncertainty band is larger in the HW − j case. This is due to the fact that for HW − production there is no renormalization-scale dependence at LO, while in HW − j such dependence is already present at leading order.
In Figs. 26 and 27 we compare the transverse momentum of the HW − pair in two  different p T ranges. Here we observe significant differences due to the fact that this distribution is only computed at leading order in the HW − simulation, while it is NLO accurate in the HW − j case. Since we included also EW corrections, in our plots these differences are slightly more pronounced than in the pure QCD implementation of Ref. [33]. The fact that such differences emerge in the region below the QCD Sudakov peak (Fig. 26) is consistent with the observation of enhanced EW effects in that region (Fig. 22) as discussed in Sec. 6.2. We also note that the uncertainty band for the HW − generator is smaller than the HW − j one. This is due again to the fact that, at Born level, HW − production does not depend upon α S , while HW − j production does, and this dependence amplifies the scale-variation band.

Summary and conclusions
In this paper we have presented the first NLO QCD+EW calculations for HV and HVj production, with V = W ± , Z, at the LHC. Specifically, we have considered complete Higgsstrahlung processes corresponding to Higgs boson production in association with off-shell ν or + − leptonic pairs plus zero or one jet. In addition to fixed-order predictions we have presented realistic simulations obtained by combining NLO QCD+EW calculations with a QCD+QED parton shower. This was achieved by means of the POWHEG BOX RES generator, a recent extension of the POWHEG BOX V2 framework, that allows for consistent NLO+PS simulations in the presence of resonances. In the case of HVj production, using the improved MiNLO approach, we have extended the applicability of NLO QCD+EW predictions to the full phase space, including regions where the hardest jet is unresolved. This is the first application of the MiNLO and POWHEG BOX RES approaches in combination with NLO EW corrections. We have studied several kinematic distributions for HV and HVj production in protonproton collisions at 13 TeV, and we have discussed predictions at fixed-order NLO, at the level of POWHEG BOX RES Les Houches events, and at NLO+PS level using Pythia 8.1. Particular care has been taken in combining the QCD+QED shower of Pythia 8.1 with the POWHEG BOX-generated events, since no standard interface is available, at present, to deal with multiple NLO emissions that can arise at production and decay level in resonant processes.
Electroweak corrections typically lower NLO+PS QCD predictions by 5 to 10% at the level of integrated cross sections and in angular distributions. We have observed quantitatively similar and rather constant EW corrections also in the jet-p T spectrum, as well as in the reconstructed Z-mass and transverse W -mass in the vicinity of the corresponding resonances. In contrast, due to Sudakov logarithms, EW corrections can be much more sizable in the tails of transverse-momentum and invariant-mass distributions. For example, in the Higgs and vector-boson p T distributions, EW corrections reach up to −25% around 1 TeV. In this respect, the HV and HVj Higgsstrahlung processes behave similarly, i.e. the emission of a jet does not have a sizable impact on EW corrections.
We have studied theoretical uncertainties associated with standard factor-two variations of the renormalization and factorization scales. In the context of the POWHEG formalism, scale variations are performed only at the level of the underlying-Born cross section, while the scale of the strong coupling constant associated with NLO radiation is kept fixed at the corresponding transverse momentum. Thus the resulting scale-variation bands are typically smaller as compared to the ones obtained in fixed-order NLO calculations. In the total cross sections for HV and HVj production we have found scale uncertainties around 1-2% and 2-4%, respectively, while scale variations in kinematic distributions are typically at the 10% level.
Thanks to the improved MiNLO prescription, simulations based on NLO QCD+EW matrix elements for HVj production can be applied to inclusive observables and compared against more conventional simulations based on NLO QCD+EW matrix elements for HV production. At NLO QCD, the observed agreement between HV and HVj predictions confirms that, as is well known, the improved MiNLO approach guarantees NLO QCD accuracy also when the extra jet is integrated out. A similarly good level of agreement was found also at NLO QCD+EW level in a variety of observables. In this regard, based on unitarity and factorization properties of soft and collinear QCD radiation, we have sketched a proof of the fact that the improved MiNLO approach, applied to QCD jet radiation computed with NLO QCD+EW matrix elements, should provide NLO QCD+EW accuracy in the full phase space.
All relevant matrix elements at NLO EW have been computed using a recent interface of the POWHEG BOX RES framework with the OpenLoops matrix-element generator. The other QCD amplitudes have been computed in part analytically and in part using the standard interface to MadGraph4. We have also presented simple analytic expression that approximate the virtual EW amplitudes in the Sudakov regime at next-to-leading-logarithmic accuracy. This approximation captures the bulk of EW corrections and reproduces exact NLO EW results with reasonable accuracy. Moreover it can be exploited in the combination of the reweighting approach that permits to speed up NLO QCD+EW simulations while providing full NLO EW accuracy in the final results.
The POWHEG BOX RES code together with the generators that we have implemented for HV and HVj production can be downloaded following the instructions at the POWHEG BOX web page: http://powhegbox.mib.infn.it

A Validation of the fixed-order NLO EW corrections in HV production
In this section we compare our fixed-order NLO EW predictions for HW and HZ production with predictions obtained with the Monte Carlo program HAWK [45].

Setup for the comparison
In order to make a comparison between the results generated by HAWK and our results, we switched off photon-initiated contributions in HAWK, since these contributions are currently not included in the POWHEG BOX RES HV generators. Similarly, bb-initiated contributions have been discarded in the POWHEG BOX RES, since this sub-process is not included in HAWK. The CKM matrix elements have been set to Photons are recombined with collinear charged leptons if R γ < 0.1, where R γ is the angular separation variable in the (y, φ) plane. If more than one charged lepton is present in the final state, the eventual recombination is performed with the lepton having the smallest value of R γ . After photon recombination, we apply the following cuts on the charged dressed leptons while for HW production we also require a missing transverse momentum of

Results
In Figs. 28 and 29 we compare NLO EW predictions obtained with POWHEG BOX RES (solid line) and HAWK (dashed line) for selected observables in HW + and HW − production.  As examples for the validation of HZ production, in Fig. 30 we present the transverse momentum and the rapidity of the Higgs boson, and in Fig. 31 the rapidity of the produced electron. Again we find a perfect overlap between the POWHEG BOX RES and HAWK predictions, and the same level of agreement was found for all kinematic distributions that we have examined.   Figure 31. NLO EW predictions for the rapidity of the electron. Same curves and labels as in Fig. 30.

B The virtual EW Sudakov approximation
The calculation of EW virtual corrections is typically more complex than in the case of QCD. This is due to the nontrivial gauge-boson mass spectrum, the presence of Yukawa and scalar interactions, the fact that EW corrections enter also in leptonic vector-boson decays, as well as subtleties related to the treatment of unstable particles. For these reasons, the numerical evaluation of NLO EW virtual corrections can be time consuming. Motivated by this practical issue, in this appendix we present compact analytic formulas that provide a decent approximation of the bulk of NLO EW effects, based on the Sudakov approximation. Besides speeding up the numerics, this approximation provides also valuable insights into the origin of the bulk of the EW corrections.
The largest EW corrections originate in the Sudakov regime, where all kinematic invariants are of the same order and much larger than the electroweak scale. In this high-energy regime, the EW corrections are dominated by Sudakov logarithms [46,48,49,[99][100][101] of the form i.e. by leading (LL) and next-to-leading logarithms (NLL) involving the ratio of the partonic center-of-mass energy √ s to the electroweak-boson masses, M V = M W , M Z . Sudakov EW logarithms originate from virtual gauge bosons that couple to one or two on-shell external particles in the soft and/or collinear limits.
General factorization formulas for LLs and NLLs that apply to any Standard Model process at one loop have been derived in Refs. [49,61,102]. For a generic n-particle scattering processes with all particles ϕ i and momenta p i incoming 10 ϕ 1 (p 1 ) ϕ 2 (p 2 ) . . . ϕ n (p n ) → 0 , Here the first term is related to the running δλ i of the dimensionless coupling parameters in the Born amplitude, while the second term consists of process-independent correction factors δ ϕ 1 ...ϕn ϕ 1 ...ϕ n that contain all LL and NLL terms and multiply Born matrix elements for the process at hand. Note that the correction factors are matrices in SU(2) space. In general they act on one or two external particles, requiring the evaluation of SU(2)-transformed matrix elements M The first three terms are due to double logarithms originating from soft-collinear gauge bosons exchanged between pairs of external legs. This gives rise to angular-independent LLs proportional to L(s) (δ LSC ) and subleading angular-dependent logarithms of type l(s) log(|r kl |/s), with r kl = (p k + p l ) 2 . The latter are split into terms associated with neutral (δ SSC,n ) and charged (δ SSC,± ) soft-collinear gauge bosons. The remaining terms consist of single-logarithmic contributions from soft/collinear gauge bosons coupling to single external legs (δ C ) and from the usual renormalization-group evolution of dimensionless coupling parameters (δ PR ). In the following, the general results of Refs. [49,61,102] are applied to Higgsstrahlung processes.

B.1 NLL Sudakov approximation for HV and HVj production
In this section we discuss the application of the Sudakov approximation to resonant processes. We focus on vector-boson decays to leptons of the first generation, but our results are applicable also to µ and τ leptons. With u and d we denote generic up-and down-type quarks, with no assumptions on their generation, unless specified.
For the leading-order kinematics of HV and HVj production at particle level we use the notation where P i = q,q, g are generic partons, and k = p 4 + p 5 is the off-shell momentum of the decaying vector boson. In HV production, P 6 is not present, and the two incoming partons are always a quark-antiquark pair, while for HVj production an extra gluon can appear, both as an initial-state parton or in the final state. In order to apply the Sudakov approximation of Refs. [49,61,102], the Higgsstrahlung processes (2.1) need to be factorized into separate parts associated with the production and decay of the vector boson. This is achieved in a gauge-invariant way by using the leading-pole approximation (LPA) [103,104], which corresponds to the leading term of a systematic expansion in Γ V /M V . At leading order, the LPA for Higgsstrahlung processes reads where factorized matrix elements for vector-boson production and decay on the r.h.s. are summed over the physical polarizations of the vector boson. The propagator in Eq. (B.6) depends on the off-shell vector-boson momentum k, while, in the matrix elements on the r.h.s., an on-shell projected momentum k must be used in order ensure gauge invariance. This can be achieved with a mapping that, conserving energy and momentum, projects on shell the V and H momenta and rescales accordingly the momenta of the decay products. In our implementation, we employ such a mapping by keeping fixed the angles formed by the vector boson and by one of the leptons.
In general, in leading-pole approximation, three types of NLO EW corrections need to be considered: factorizable corrections to the production and decay parts, and nonfactorizable corrections that connect production and decay. However, the latter are typically quite small [105][106][107][108]. Moreover, vector-boson decay do not involve Sudakov EW logarithms. Thus, only the production part receives Sudakov EW corrections, i.e.
and δM P 1 P 2 →HV λ (P 6 ) as well as its decay counterpart need to be computed for both transversely-and longitudinally-polarized vector bosons. In the framework of Refs. [49,61,102], tree amplitudes with longitudinal vector bosons need to be related to corresponding amplitudes with Goldstone bosons using the Goldstone-boson equivalence theorem [109][110][111] where V a i L are the longitudinal gauge bosons, Φ a i the corresponding Goldstone bosons, ϕ i are the fermions and scalars in the process, M and E are typical scale masses and energies involved in the process, d is the mass dimension of the matrix element and Q V a k is the electric charge of the vector boson V a k .
In the following sections we present analytic results for all relevant NLL EW correction factors. These formulae contain group-theoretical quantities such as the electric charge Q of the external particles, their weak isospin T a , or the electroweak Casimir operator C ew . Their values can be found in App. B of Ref. [102]. For the sine and cosine of the Weinberg angle, we use the shorthand s w and c w , respectively.
Large logarithms of the light-fermion masses do not need to be included since we use the G µ scheme, which incorporates the running of the electromagnetic coupling up to the EW scale, and we regularize QED infrared singularities of virtual type at the EW scale by using an effective photon mass λ = M W . This approach effectively corresponds, in logarithmic approximation, to the combination of virtual EW corrections with the emission of real photons up to transverse momenta of the order of M W .
In the framework of the Sudakov NLL approximation, the Sudakov limit is applied only to virtual EW effects, while real QED radiation is treated exactly. More precisely, FKSsubtracted real-emission matrix elements are treated exactly, while only the finite part of the integrated FKS terms, defined via MS subtraction of the IR poles at the scale µ = µ R , is included. Concerning IR singularities, this MS subtraction is consistent with the cancellation of virtual QED singularities through the above mentioned λ = M W regularization approach. However, as far as QED logarithms are concerned, we note that we do not apply a fully consistent matching of the (regularized) virtual contributions to real QED radiation. In fact, the former are effectively cut off at the scale M W , while the latter are subtracted at the scale µ = µ R . This implies missing logarithmic terms of order α EM ln(µ R /M W ). Nevertheless, as demonstrated by our numerical results, such uncontrolled logarithms do not jeopardize the accuracy of the Sudakov approximation at high energies.

B.2 HW and HWj production
Here we focus on HW − production and we first consider the partonic process which involves only left-chiral quarks and leptons. Matrix elements for the charge-conjugated processd as well as crossing-related matrix elements corresponding to permutations of the initialstate quarks, can be easily obtained from the ones for the processes (B.9). For the crossed production process the Born amplitudes in the high-energy limit read where q = p 1 + p 2 . Transverse and longitudinal gauge-boson polarization states are denoted as λ = ±1 ≡ T and λ = 0 ≡ L, respectively, and For the decay of the polarized W − boson we have For the crossed HW − j production process the polarized Born amplitudes at high energy read where g s is the strong coupling, t a is the color matrix, and The gluon-initiated processes can be obtained via appropriate crossing transformations.
For longitudinally polarized W bosons, we obtain the Sudakov correction factors where the group theoretical quantities involving the charged Goldstone boson φ − arise from the Goldstone-boson equivalence theorem, and the explicit expression for the SU(2) β-function coefficient b ew W W can be found in App. B of Ref. [102]. The correction factors of Eqs. (B.22) and (B.24) are equally valid for the HW − and HW + production processes in Eqs. (B.9) and (B.10). Corresponding results for processes with the initial-state quarks interchanged are easily obtained by swapping the Mandelstam variables t and u.

Sudakov correction factors for HWj production
The Sudakov correction factors for HWj production are quite similar to the ones for HW production. In fact, the presence of an extra SU(2)×U(1) singlet gluon has only indirect effects of kinematic type on the Sudakov EW corrections. In particular, the δ LSC , δ C and δ PR factors of Eqs. (B.22) and (B.24) are directly applicable to HWj production without any modification. In contrast, the δ SSC,n and δ SSC,± factors need to be generalized by including extra angular-dependent logarithms of type log(|r 12 |/s) associated with vectorboson exchange between the initial-state quarks. For dū → HW − this kind of logarithms vanishes due to r 12 = s. However, in the case of dū → HW − g, they need to be taken into account, since they give rise to non-vanishing contributions via crossing transformations of the type r 12 ↔ r 16 , which correspond to the case of quark-gluon initial states.
In the transverse case (λ = T ) the SSC correction factors become while for longitudinal W − bosons (λ = L) they read The above correction factors are directly applicable to HW + j production as well, while processes with the initial partons exchanged require the swap r 13 ↔ r 23 and r 1k ↔ r 2k .

B.3 HZ and HZj production
One of the main differences between the Sudakov EW corrections for HZ and HW production is due to the fact that Z bosons couple to both left-and right-handed currents. As a consequence, for the HZ production process the squared Born amplitude in LPA reads Here and in the following we keep track of the quark chirality κ = R, L in the grouptheoretical quantities, while q = p 1 + p 2 , and A κ T Z = −iv κ (p 2 ) γ µ u κ (p 1 ) T µ (−k) , (B.33) A κ LZ = −iv κ (p 2 ) γ µ u κ (p 1 ) (−k + p 3 ) µ . (B.34) The interchange of the initial-state quarks modifies only the spinor part, without changing the structure of the matrix element. The Z-decay matrix element for generic λ reads be evaluated. This has the advantage that the bulk of the inclusive cross section is computed with high statistics at a reduced computational cost.
-Then, the code is run again with lower statistics by using the same importancesampling integration grids generated in the first run, and computing only the missing EW part of the virtual contributions. This is done by setting the flags virtonly and qed_qcd to 1 in the input file. In this way, if the flag select_EW_virt is set to 1, the complete virtual corrections are included. If instead the user is interested in obtaining the Sudakov approximated results, the code can be run by setting select_EW_virt to 2.
Finally, the kinematic distributions obtained in the two previous steps should be combined, by summing them together.

Les-Houches-level Monte Carlo events
If one is interested in generating Monte Carlo events at the Les Houches level, then the code can be run with some approximation of the virtual contributions (or even with the virtual corrections set to zero). At the end, the generated events need only to be reweighted, using the POWHEG BOX RES reweighting feature, with the full virtual corrections activated. In this way, when the event is generated, the use of an approximated virtual contribution (whose evaluation could be requested several times per event) considerably speeds up the code. Once the event is generated, the reweighting procedure calls the full virtual contribution only once per event. This reduces the running time in a drastic way. The different options available are the following: -The user can generate the events omitting virtual contributions of QCD and EW kind. This is obtained by setting the POWHEG BOX RES flag novirtual to 1 in the input file. Since the inclusive cross section used to generate the weight associated to the event is computed without the finite part of the virtual corrections, the weight associated with a single event can be very different with respect to the weight obtained after the reweighting procedure applied on that event. This can give rise to statistical fluctuations in the kinematic distributions, that would need a higher number of events to be smoothed.
-The user can generate the events including only the QCD part of the virtual contributions, that are quite fast to be evaluated. This is done by setting the POWHEG BOX RES flag select_EW_virt to 0 in the input file. The difference with the previous case is that an important part of the virtual corrections is included in the calculation of the inclusive cross section, and the results after reweighting tend to be smoother.
-The best option is to include the QCD part of the virtual corrections together with their EW NLL approximation. This is achieved by setting select_EW_virt to 2. This option is the one we have used to generate the events analyzed in Sec. 5. Since the EW NLL approximation of the virtual corrections captures most of the dominant Sudakov logarithms, running the code with this setting generates events whose weight is very similar to the final weight associated to each event after reweighting.
The default value for the select_EW_virt flag is 1, which corresponds to the inclusion of exact virtual EW contributions at all stages.

D Interface to Pythia 8.1 and the veto procedure
In order to generate realistic event samples at NLO+PS accuracy, including both QCD and QED corrections, the radiation of QCD partons and photons generated at the LHE level by the POWHEG BOX RES framework has to be completed by a Monte Carlo showering program. This is achieved through a dedicated interface that feeds the LH events to Pythia 8.1.

The veto procedure
In the following we discuss the veto procedure that is applied in order to guarantee a consistent combination of QCD and QED radiation generated at LHE level with subsequent parton-shower emission. Since we apply the multi-radiation mode described in Sec. 3.1, each LH event generated by the POWHEG BOX RES framework can be accompanied by both QCD and QED radiation. Radiation of QCD type arises only at the "production" level, while photon radiation can come both from "production" and from the charged leptons that arise from the decays of Z and W resonances.
For QCD radiation, the standard veto shower implemented in any Monte Carlo program is used. In practice, the highest transverse momentum of the radiation (of QCD or QED type) generated at the "production" level by the POWHEG BOX RES framework is passed to the shower Monte Carlo program through the variable scalup, in the Les Houches interface.
For what concerns QED radiation, since the Les Houches interface does not provide a standard mechanism to veto radiation from resonance decays, we have implemented a dedicated veto procedure. The POWHEG BOX RES events can have up to two photons at LHE level, one associated with the production part, and one with the decay part of the process, and the shower Monte Carlo has to be instructed to veto, separately at the level of production and decay, any photon with transverse momentum higher than the hardness of the emissions produced by the POWHEG BOX RES framework.
To this end, we first scan the Les Houches event to identify the photons that have been generated by the POWHEG BOX RES framework, determining if their mother belongs to the production or to the decay products. 12 We then shower the event with Pythia 8.1, restricting QCD radiation by means of the scalup variable as discussed above, and identifying the extra photons produced by the shower algorithm.
For photons that are generated by the shower at the production level we apply a similar veto procedure as for QCD radiation: 1. we compute the transverse momentum of each photon produced by the shower, at production level, and store its maximum value p max T for the event at hand; 2. if p max T is greater than scalup the event is vetoed. This procedure effectively amounts to requiring that, at production level, the shower does not generate any QED radiation with transverse momentum greater than the radiation of QCD or QED type produced by the POWHEG BOX RES.
3. Since, in order to ensure momentum conservation, the Monte Carlo reshuffling procedure during shower generation slightly modifies the momenta of the particles, we also check that p max T does not exceed the hardness of the LH photon after reshuffling. If this happens, the event is vetoed.
For photons associated to the resonance decay, we proceed as follows: 1. if no photon is present at the LHE level, this means that the POWHEG BOX RES has not been able to generate radiation harder than the minimum value of 10 −3 GeV, set as a minimum for the transverse-momentum of photon radiation. 13 In this case, any shower QED radiation harder than 1 MeV is vetoed in the decay.
2. If instead a photon is already present, we compute its transverse momentum with respect to the lepton emitter in the center-of-mass frame of the mother resonance, and store this value in p max T, rel . In HZ and HZj production, at Les Houches level, it is not possible to know if the photon has been emitted by the lepton or by the antilepton, and p max T, rel is set to the minimum value between the two relative transverse momenta. We then veto the event if, among the photons produced at decay level, the maximum relative transverse momentum is greater than p max T, rel .