Improving NLO QCD event generators with high-energy EW corrections

In this work we present a new approach for the combination of electroweak (EW) corrections at high energies, the so-called EW Sudakov logarithms (EWSL), and next-to-leading-order QCD predictions matched to parton-shower simulations (NLO+PS). Our approach is based on a reweighting procedure of NLO+PS events. In particular, both events with and without an extra hard emission from matrix elements are consistently reweighted via the inclusion of the corresponding EWSL contribution. We describe the technical details and the implementation in the MadGraph5_aMC@NLO framework. Via a completely automated procedure, events at this new level of accuracy can be obtained for a vast class of hadroproduction processes. As a byproduct we provide results for phenomenologically relevant physical distributions from top-quark pair and Higgs boson associated production ($t\overline{t}H$) and from the associated production of three $Z$ gauge bosons ($ZZZ$).


Introduction
After the first two runs of the Large Hadron Collider (LHC), our knowledge of the fundamental interactions of elementary particles has tremendously improved.Above all, the Higgs boson has been observed [1, 2] and its properties have been found compatible [3] with those predicted by the Standard Model (SM), further corroborating the current theory of fundamental interactions.On the other hand, no clear sign of beyond-the-SM (BSM) physics has been found at the LHC so far, as it has been the case for previous colliders.In 2022, ten years after the discovery of the Higgs boson, the Run-3 has started and during this period and the subsequent High-Luminosity (HL) runs [4][5][6][7][8][9], the total amount of recorded data by the LHC will increase by a factor of 20 w.r.t. the Run-1 and Run-2 data sets combined.Moreover, several options have been proposed for future colliders, involving collisions at higher energies between protons and/or leptons (including also muons).It is therefore clear that the quest for new physics at colliders is only at its initial stage.
The success of this quest relies on the availability of precise and accurate SM predictions.In other words, the possibility of calculating QCD and electroweak (EW) higher-order corrections.Regarding fixed-order perturbative expansion, QCD radiative corrections at Next-to-Leading-Order (NLO), Next-to-NLO (NNLO) or even Next-to-NNLO (N 3 LO) accuracy are nowadays available for several processes.In fact, both NLO QCD and NLO EW corrections can be calculated for processes with high-multiplicity final states, with limitations given only by computing power.This is possible since such corrections have been implemented in Monte Carlo generators and their calculation has been even automated [10][11][12][13][14][15][16][17], at different levels in the different frameworks, using different one-loop matrix-element providers [18][19][20][21][22][23][24].
Another unavoidable ingredient for the simulation of events at (hadron) colliders is the modelling of the multiple emission of (QCD) partons and their hadronisation, i.e., Parton Shower (PS) simulations.However, while the matching of NLO QCD corrections and parton shower effects has already been achieved [25][26][27] (also for NNLO [28][29][30][31] and recently even N 3 LO accuracy [32] for specific processes) and automated since a long time, in the case of NLO EW corrections a process-independent approach still needs to be formulated and either only approximations or case-by-case exact solutions have appeared in the literature so far [33][34][35][36][37][38][39][40][41].It is therefore desirable that the lack of a general algorithmic procedure for the exact matching of NLO QCD+EW predictions and PS simulations is solved as soon as possible, especially since it is well known that NLO EW corrections can strongly depend on the kinematics, and the naive estimate of their relative impact (NLO EW∼ O(α) ∼ O(1%) in absolute value) can be easily violated by one or more order of magnitudes.On the one hand, the origin of this violation can be due to a specific mechanism for a specific process and/or observable [42][43][44][45][46][47][48].On the other hand, the origin of this violation (NLO EW ≫ 1% in absolute value) is typically related to two different kinds of effects, which are universal.First, the final-state-radiation (FSR) of photons from light fermions, which is of QED origin and, e.g., distorts the Breit-Wigner distributions of the Zboson decay products.The modelling of FSR is available within modern PS simulators, such as Pythia8 [49][50][51], Herwig7 [52], and Sherpa [53][54][55].Second, the EW Sudakov logarithms (EWSL) [56] of the form α n log k (s/M 2 W ) with k ≤ 2n, which are mainly of weak origin and become relevant at high energies (s ≫ M 2 W ).An algorithmic procedure for their evaluation at one- [57,58] and two-loop [59][60][61][62] accuracy is available since long time, the so-called Denner and Pozzorini (DP) algorithm for the Sudakov approximation.It has been automated for the first time [63] in the Sherpa framework, and, after the revisitation and improvement of particular features [64], in the MadGraph5_aMC@NLO framework [10,16].
It is clear therefore that having automated tools for generating events at NLO accuracy matched to PS simulation (NLO+PS) where not only NLO QCD corrections (NLO QCD +PS) but also the dominant NLO EW ones, FSR and EWSL, are taken into account is very useful for current and future experimental analyses, especially if the addition of the EW contributions does not slow down the generation of the events.An example of a tool of this kind based on Ref. [63] and the Meps@Nlo [65] method has already appeared in the literature [66] and it has been applied to a specific process: ZZ and ZZj merged production.This work has shown the relevance of such studies and the advantages of a general-purpose automation for (at least) SM processes in general.
The current work precisely presents the automation, in the MadGraph5_aMC@NLO framework, of combined EWSL and NLO QCD + PS accuracy, including QED FSR, for the event generation of SM processes.We have implemented it in this framework, since it offers all the capabilities for achieving this goal.Indeed, in MadGraph5_aMC@NLO the following three features are already available: • The automated calculation of matched NLO QCD + PS simulation, via the interface with external PS simulators.
• The automated evaluation of EWSL at one-loop accuracy.
• The possibility of reweighting events (e.g.changing model parameters) from LO and NLO simulations [67].
Our strategy therefore is based on the reweighting of NLO QCD + PS events taking into account the EWSL contribution.FSR can then eventually be simulated directly via the PS.In doing so, we do not reweight with the EWSL only the LO contribution from the hard process, but also the QCD one-loop virtual contribution as well as the first QCD real emission.For the latter, we take into account that both the kinematics and the external states are different.In this way, especially for high-energies (s ≫ M 2 W ), a good approximation of NLO EW corrections is correctly taken into account both for the Born-like process and the one with an extra hard jet.Automatically, the NLO QCD prediction for the inclusive production is correctly reweighted via the EWSL, and QCD shower effects are taken into account.Moreover, by adopting the so-called SDK weak scheme [64] for the EWSL, which consists of a complete removal of QED effects of infrared (IR) origin, FSR or in general QED effects can be included in the PS simulation avoiding their double counting.
Since the evaluation of EWSL involves only tree-level matrix elements and compact analytical formulas (see Ref. [64]), one of the advantages of reweighting via the Sudakov approximation is the speed of this procedure and especially the numerical stability of the results.However, especially for the real emission contributions, it is crucial that the EWSL are damped in phase-space regions where any of the kinematical invariants involving two external states is smaller than M 2 W .This is necessary not only because in such phase-space regions the Sudakov approximation is not valid, but also because the soft and collinear limits relating n + 1 and n final states in QCD must be preserved for a correct matching of NLO QCD predictions and PS simulations also after the reweighting.The approach we have adopted for solving this problem is completely general and, although we will discuss in the paper its implementation in MadGraph5_aMC@NLO, could be in principle extended for other matching schemes.Since the approach is based on the reweighting, i.e. a step happening after the event generation, it does not rely on the strategy for the event generation itself.
In this paper we focus on the technical implementation and its validation, leaving phenomenological studies and comparison with exact NLO EW accuracy predictions to dedicated works.Nevertheless, we show results for two representative SM processes: the top-quark pair and Higgs associated hadroproduction (pp → t tH) and the hadroproduction of three Z gauge bosons (pp → ZZZ).We consider for both processes the final-state particles as stable, while for the latter we also consider the case of Z bosons decaying into e + e − pairs, enlightening the relevance of considering both EWSL and QED FSR contributions.In doing so, we use Mad-Spin [68] for performing the decays and therefore, while the reconstruction of the tree-level spin correlations are taken into account, we do not preserve the information of the correlation of the helicity-dependent EWSL with the helicities and angular distributions of the decay products. 1he paper is organised as follow.In Sec. 2 we discuss the motivations for this work and the general structure of our implementation in MadGraph5_aMC@NLO of an automated framework for performing NLO QCD +PS simulations including also the effects of EWSL.In Sec. 3 we give the technical details of the implementation of the reweighting approach in MadGraph5_aMC@NLO.
In Sec. 4 we present results for t tH and ZZZ hadroproduction at 13 TeV collisions.In Sec. 5 we give our conclusion and outlook.Finally, in Appendix A we briefly summarise the DP algorithm as revisited in Ref. [64], focusing on the concepts and the formulas that are relevant for the discussion of Sec. 3.

Overview of the problem and proposed solution
As already mentioned in the introduction, in this section we discuss the motivations for this work and the general structure of the implementation in MadGraph5_aMC@NLO of an automated framework for performing NLO QCD +PS simulations including also the effects of EWSL, i.e., NLO EW corrections in the Sudakov approximation.For the very interested readers, the technical details can be found in Sec. 3. In the following, we start by introducing the notation, which will be also used in Sec. 4 where we present numerical results for selected processes.

Notation
At fixed order, adopting the notation already used in Refs.[70, 12, 71-73, 46, 16, 74-77, 17], the different contributions from the expansion in powers of α S and α to the differential or inclusive cross section Σ of a generic (SM) process up to NLO can be denoted as: where k ≥ 1 and is process dependent.At LO, meaning tree-level diagrams only, there can be more than one perturbative order α n S α m and each one of the orders is associated to a different . Similarly, including one-loop corrections, there is more than one perturbative order and each one is associated to a different Σ NLO i , and if . The full set of LO and NLO orders is the so-called Complete-NLO order.
In this paper we are interested, from the fixed-order side, in the quantities ) ) ) ) ) The rest of the Complete-NLO prediction, Σ LO i with i > 1 and Σ NLO i with i > 2 is not considered and is not relevant for the discussion in this paper.Still it is worth to mention that, as it is already very well known (see e.g.Ref. [70]), the NLO EW corrections, i.

State of the art and general considerations
The reason why at high energies, and when all the invariants are large, EWSL are a very good approximation of NLO EW corrections (the NLO 2 term) is that which is a non-negligible contribution only if we are considering a process for which a precision of very few percent is relevant.EWSL are precisely the ingredient that is missing in many present and future experimental analyses at the LHC that are targeting particles with transverse momenta larger than or equal to a few hundreds of GeV.Due to statistics and other uncertainties of both experimental and theoretical origin, these analyses are sensitive to effects of O(10%).Thus, in those cases, EWSL cannot be ignored.The same analyses cannot ignore also the effects due to NLO QCD corrections and clearly PS simulations.
While the matching of NLO QCD+EW with PS (NLO QCD+EW + PS) is in general far from trivial, the case of NLO QCD+EWSL (NLO QCD+EWSL + PS) is indeed trivial once a NLO QCD + PS matching is already available.In the NLO QCD + PS matching, PS can simulate only effects of O(α n S ) with n > 0 on top of NLO QCD .The matching consists in avoiding the double counting of the case n = 1 on top of the LO QCD component.If now we start with the NLO QCD+EWSL at fixed order, PS effects on top of the EWSL components will induce effects of O(αα n S ) with n > 0 on top of the LO QCD , which cannot be double-counted by construction.
In the previous argument we have implicitly assumed that PS involves only QCD effects, but in modern PS simulators also QED effects such as FSR, possibly involving also photons splitting into fermions, are taken into account.In those cases, the PS simulations will induce effects of O(α m α n S ) with m + n > 0 on top of the LO QCD component, which can instead lead to a possible double-counting of the EWSL in the case m = 1, n = 0.The solution in this case is using the scheme denoted as SDK weak in Ref. [64].This scheme completely neglects effects of pure-QED origin in the evaluation of EWSL and therefore takes into account only the purely weak ones, avoiding the double-counting (more details are given in Appendix A).The QED component of the EWSL would anyway vanish in observables inclusive in the photon radiation, but otherwise large QED effects are simulated directly by the shower. 2These effects include not only the QED component of the EWSL but also FSR effects such as the already mentioned distortion of the Breit-Wigner distribution from the Z boson decay.This means that it is possible to generate NLO QCD+EWSL events for the production of heavy objects (W, Z, H bosons or top quarks) and let them decay directly via the shower or programs like MadSpin [68] including QED FSR effects from the shower in NLO QCD+EWSL + PS simulations. 3We stress the fact that with the notation "PS" we understand the presence of both QCD and QED effects in the shower. 4When we will refer to the purely QCD effects in the PS we will use the notation PS QED .The reader may wonder what would be the problem if in the previous argument instead of considering EWSL, the exact NLO EW corrections were considered, namely the NLO QCD+EW + PS case.First of all, it is important to note that claiming NLO EW accuracy means having under control the exact O(α) effects, together with the advantages of shower simulations, meaning e.g. the possibilities of setting hard jet-vetoes or studying the transverse momentum of the total final-state system without obtaining the typical pathological results of fixed-order simulations, where large logarithms are not resummed.For this reason not only the QED shower on top of LO QCD but also the standard QCD part of PS on top of the LO 2 must be taken into account and matched to the NLO EW corrections.Especially the latter contribution poses non-trivial 2 In Ref. [64] the SDK weak has been devised in order to take into account the cancellation between the virtual and real QED components.This cancellation has not been formally proven but it is supported by several examples that are reported and discussed in Ref. [64].One should notice that here we do not want to take into account this cancellation, instead to completely remove the QED component and leave its simulation to the PS.Therefore, the choice of the SDK weak approach is not motivated by the cancellation between virtual and real QED contributions. 3The approach described here would fail in case of the inclusion of purely weak effects directly in the shower [78][79][80][81].However, the multiple emission of heavy bosons (denoted later in the text as HBR) are not relevant for 10-100% precision even at a 100 TeV proton-proton collider [82].Similar consideration applies for the evolution of proton PDFs involving weak splittings [83,84].We note that the case of a high-energy lepton collider is a completely different scenario [85][86][87][88].
challenges, since the colour flow is not defined as the LO 2 contribution is typically an interference and not a squared amplitude.Intermediate solutions to these problems have been proposed (see Refs. [35,37,36] for applications to phenomenology, possibly extended to multi-jet merged processes, and Ref. [89] for more formal aspects related e.g. to the definition of colour flows for interferences), and results with an exact matching at NLO QCD+EW + PS have been presented only for cases where the LO i , with i > 1 are not present [33,34,[38][39][40], recently even at NNLO QCD+EW + PS accuracy [41].However, a general method has not appeared so far in the literature.
Even if a general NLO QCD+EW + PS method were available, the possibility of obtaining NLO QCD+EWSL + PS is still desirable for practical reasons: speed and stability.EWSL can be calculated via compact analytical formulas and tree-level matrix elements only, as explained in Appendix A. Therefore, unlike exact NLO EW corrections, they can be evaluated in a much faster and numerically stable way, without numerical cancellations between real and virtual contributions.
That said, we want to stress that the EWSL are an approximation of the exact NLO EW corrections and therefore there could be non-negligible effects at high energies that cannot be captured, such as photon-initiated contributions (see, e.g., Ref. [90]).In general, given the possibility of performing at fixed order both the calculations of EWSL and of exact NLO EW corrections, for any phenomenological study, we strongly suggest to compare the two approaches beforehand.

Proposed solution: NLO QCD ⊗ EWSL + PS
In this work we present the automation not only of the NLO QCD+EWSL + PS, and also of its simpler version at LO denoted as LO QCD+EWSL + PS, but also of what we will denote as which will be our best prediction and will be described briefly in the following, leaving the technical details to Sec. 3. First of all, predictions at NLO QCD ⊗ EWSL + PS accuracy are obtained by showering events at NLO QCD ⊗ EWSL accuracy: at variance with the NLO QCD+EWSL case, not only the LO QCD contribution but also the NLO QCD corrections (Σ NLO 1 ) receive corrections from EW Sudakov logarithms of O α log k (s/M 2 W ) with k = 1, 2. It is important to note that the NLO QCD corrections originate from both virtual and real contributions.While the former have the same external states and kinematics as the Born contribution, the latter are different.The EWSL have to be therefore evaluated separately for the virtual and real contributions.
Using a notation that will be exploited in Sec. 3, we denote with δ EWSL (S) the relative impact of the EWSL on top of the LO QCD , Similarly, considering the same process plus the radiation of a Hard QCD parton, the corresponding quantity is instead denoted with δ EWSL

(H)
. As guiding principle, in the NLO QCD ⊗ EWSL predictions we want that, similarly to the LO QCD contribution, also NLO QCD virtual and real contributions in the Soft/collinear regions receive δ EWSL (S) corrections.Conversely, the contribution from hard and non-collinear real emissions should be corrected by δ EWSL

(H)
. At the same time, we need to ensure that in the soft/collinear limits .
(2.12) Condition (2.12) is unavoidable for two different reasons: 1. Virtual and real IR poles need to receive the same corrections from EWSL, so that the cancellation of the divergences is preserved.
2. In the soft and/or collinear limits at least one of the kinematical invariants involving two external states is by definition smaller than M 2 W , invalidating the applicability of the Sudakov approximation and the sensibility of δ EWSL (H) .Condition (2.12), leaves freedom on how to implement the mapping between δ EWSL (H) and δ EWSL (S) and we will discuss the practical implementation in Sec. 3.
The strategy adopted for correcting the different contributions by either δ EWSL (S) is very similar to the one used in, e.g., Ref. [91] and relies on the general framework introduced in Ref. [67]: reweighting NLO events before showering them.We will give more details in Sec. 3, but the idea is the following.In the MC@NLO formalism two kinds of events are generated, namely the S and H events.The latter class corresponds to the contribution from hard real emission to NLO QCD corrections.It takes into account the contribution of the Monte Carlo (MC) counter term, which is precisely added in order to avoid the double counting from PS effects on top of the LO QCD .On the contrary, the rest of the contributions entering the NLO QCD predictions corresponds to S events.Given a process pp → X, with X having multiplicity n, S events are of the kind 2 → n, while H events are of the kind 2 → n + 1. Denoting the weights of the former as w S and the weight of the latter as w H5 , the events generated at NLO QCD accuracy with the MC@NLO matching scheme can be promoted to NLO QCD ⊗ EWSL accuracy performing the following reweighting before the parton shower: After the reweighting, events can be showered obtaining predictions at NLO QCD ⊗ EWSL + PS accuracy.Again, we will give many more details on the procedure in Sec. 3.

Technical details of the EW Sudakov reweighting strategy
In Sec. 2 we described the general features of the NLO QCD ⊗ EWSL + PS approximation and the motivations behind it.In this section we provide the technical details of the reweighting procedure, which we have implemented by extending the general-purpose reweighting module of

MC@NLO matching and reweighting
The structure of a fixed-order NLO calculation of a cross section dσ, as performed within Mad-Graph5_aMC@NLO, for a 2 → n production process can be summarised by the following equation dσ The terms B, V, R are respectively the Born, virtual and real emission contributions.The term C is the local counterterm that renders the integral over the dϕ n+1 phase-space finite, where dϕ n+1 ≡ n+1 k=1 d φk and d φk is the differential of the phase-space integration associated to the particle k.The term C int is the integrated form of C over dϕ n+1 /dϕ n , such that C int − C d φn+1 = 0.The specific form of the counterterms depends on the subtraction scheme that is used, e.g.FKS [92] or CS [93], where the former is the one on which the Mad-Graph5_aMC@NLO implementation is based.
In the case of matching of NLO QCD computations with PS in the MC@NLO formalism, on top of the local counterterm C one has to also include the so-called Monte-Carlo counterterm C MC [25], in order to avoid the double counting of PS effects on top of the B contribution.The counterterm C MC accounts for the cross section one obtains from PS simulations by truncating the perturbative expansion at O(α m+1 S ), where LO 1 is of O(α m S ).The MC counterterm depends on the specific PS simulator one interfaces the calculation to 6 , but since the leading IR behaviour of any PS simulator is the same as the one of R (or equivalently −V after integrating over d φn+1 ), the analogue of Eq. (3.1) for NLO QCD + PS simulation is where dσ (S) and dσ (H) are the cross sections associated to the S and H events, respectively. 7nlike fixed-order calculations (see Eq. (3.1)), MC counterterms are such that the dσ (S) and dσ (H) subtracted cross sections are separately finite and therefore Born-like (S) and real-emission (H) events can be unweighted.
The reweighting prescription of Eqs.(2.13) and (2.14) corresponds to )dσ (S) , (3.4) )dσ (H) . (3.5) As can been easily seen in Eq. (3.3), the exact cancellation between the term C and C int is preserved.The cancellation of the C MC dependence between Eqs. (3.2) and (3.3) is instead more subtle and relies on the condition (2.12).Before giving details on the implementation of δ EWSL , and especially the functional form of the mapping between them for ensuring this condition, we discuss the implications of Eq. (3.4) and (3.5).
First of all, since O(δ ) = α, all features of pure-QCD origin are exactly preserved.Considering EW interactions, the term EWSL in Eq. (2.6) is given by Bδ EWSL .Similarly, the LO QCD+EWSL + PS approximation is obtained by keeping only the term Bδ EWSL

(S)
. Both approximations can be achieved in a much easier way without reweighting, but rather accounting directly for the effects of EWSL in the event generation.This is precisely how we perform such simulations.
The NLO QCD ⊗ EWSL + PS accuracy consists instead of showering events generated at NLO QCD ⊗ EWSL accuracy, which in turn consists of the reweighting procedure of Eqs.(3.4) and (3.5) applied on events generated via the MC@NLO approach (Eqs.(3.2) and (3.3)).From the arguments of the previous paragraph it is clear that the NLO QCD + PS, LO QCD+EWSL + PS and NLO QCD+EWSL + PS accuracies are still valid within the NLO QCD ⊗ EWSL + PS one, which therefore is superior.More specifically, ) where in the previous three equations we have specified only the leading term in the combined α S and α expansion.Thus, what is left for discussion is the consistency of our approach for higher orders and in particular the combination of NLO QCD corrections, EWSL and PS effects: the terms of O(α S α) (and higher) indicated in Eq. (3.6).Implicitly in the C MC terms there is a dependence on the shower scale µ S .Roughly speaking, emissions of partons at an energy smaller than µ S are dealt by the PS simulator, while the first emission at energy larger than µ S is given by the matrix element R. 8 An NLO QCD + PS simulation preserves NLO QCD accuracy matching it to Leading-Log (LL) accuracy in QCD for soft and collinear emissions.Still, a µ S dependence, beyond the aforementioned accuracy, is left and it can be exploited for estimating higher-order effects.
Since the NLO QCD ⊗ EWSL + PS accuracy is an ad hoc approximation for accounting for the dominant EW corrections together with PS effects, the question "At what order is the µ Sdependence emerging?" is rather academic.We want to elaborate anyway on that in the following, since it will help to understand the consistency of our approach and the relevance of the 1 st motivation for the prescription (2.12).Concerning the pure-QCD contributions, it is exactly the same situation of NLO QCD + PS: it appears beyond NLO QCD accuracy.Taking into account EW corrections, the µ S -dependence can emerge only at one order of α S beyond the EWSL in Eq. (2.6).
In addition to this, if µ S ∼ M W ≪ √ s, in the relevant region of the matching where the transition between S and H events take place, condition (2.12) actually takes the form δ EWSL . As it will be explained in Sec 3.2, the condition (2.12) is ensured if any of the invariants is smaller or equal to M 2 W . Therefore no dependence on µ S related to EWSL is present at all.If instead µ S ∼ √ s > M W , a dependence on µ S can be present at one order α S beyond the EWSL in Eq. (2.6).It is important to note that this dependence is often due to the unbalance between the δ EWSL (S) for a given process of the form (A.2) and δ EWSL (H) of the same form with an additional gluon (either in the initial or final state), which does not interact electroweakly.Thus, these effects are in fact expected to be even smaller than their naive estimate: O(α S )×O(EWSL).We show a concrete example of what we have discussed in this section.In Fig. 1 we consider the case of ZZ hadroproduction where we have set the cuts p T (Z) > 600 GeV , m(ZZ) > 1200 GeV , (3.9) in order to probe the region µ S ∼ √ s > M W .We show results for two different values of the shower scale µ S , in particular µ S = k × H T /2 with k = 1 and k = 0.2.The case k = 1 is the standard value 9 , while k = 0.2 is an ad hoc value which has been chosen just for our purpose.
In the main panel of Fig. 1 we show the NLO QCD + PS predictions for the two different values of k, denoted in the plot as S + H (blue) and separately the results for the S (orange) and H (green) events alone.The case k = 1 corresponds to the solid lines, the case k = 0.2 to the dashed ones.It is important to notice that different values of k return very different contributions from the S events and the H ones and a non-negligible dependence on k, as expected, is left also in the total prediction S + H.In the first inset we show the ratio for the three different sets of events and two different shower scales.It is manifest how the impact of the EWSL is very different for the S, or H, events alone when k = 1 or k = 0.2, while in the case of the full set of events S + H there is almost no dependence on the shower-scale choice.This supports our previous argument regarding the fact that although a dependence on µ S of O(α S ) × O(EWSL) can be present, it is actually expected to be even smaller.

EWSL EWSL EWSL
Figure 1: Technical test that the relative impact of EWSL on the S + H samples depends very mildly on the value of µ S .The representative case of ZZ hadroproduction with the cuts in (3.9) is shown.
In the last inset we show the ratio between the predictions for k = 1 and k = 0.2 for the three different sets of events.The case of the NLO QCD ⊗ EWSL + PS predictions correspond to the dotted lines while the NLO QCD + PS case to the solid ones.First, we can see that for each set of events (also for the set S + H) there is a visible difference between the case k = 1 and k = 0.2.Second, one can notice how the ratio is unaffected by the presence of the EWSL contribution. .In this work, we employ the following approach:

Implementation of δ EWSL
where δ EW LA is the quantity defined in Eq. (A.11) and we have specified that it is evaluated for an S event, denoted as e S , which is associated to a process of the form e S : φ i 1 (p 1 ) . . .φ in (p n) → 0 . (3.12) In Eq. (3.12) we have used the same notation of Eq. (A.2) and understood that for a 2 → n process n = n + 2. Equation (3.11) is actually refining the definition that was given in Eq. (2.11).The EWSL are calculated in the SDK weak scheme, as described in Appendix A. This scheme was conceived in Ref. [64] in order to reproduce as close as possible NLO EW corrections.Here, the final goal is the same, but the SDK weak is actually employed in order to not double-count QED effects from PS simulations.Equation (3.11) implies also that we assume δ QCD LA = 0.This assumption has clearly no effect for all the processes for which the LO 2 is zero or anyway smaller than α/α S ∼ 0.1, i.e., the naive expectation for O(LO 2 /LO 1 ).However it is also a reasonable assumption for a much larger class of processes.Indeed, even if O(LO 2 /LO 1 ) ∼ α/α S , according to Eq. (A.12) and the related discussion, at least one of the following conditions must be satisfied for δ QCD LA to be in practice relevant: • Σ LO 2 has a sizeable dependence on m t , e.g., due to the Yukawa interaction of the top quark, and therefore there is a dependence on the parameter renormalisation of m t in QCD, (δm t ) QCD .
• The LO 2 involves matrix elements for partonic processes with external gluons.
• s ≫ µ 2 R and n − 1 ̸ = n g .As examples, any purely EW process such as multi-boson production is free of these issues since LO 2 is not present in those cases.The processes involving top quarks in the final state are also typically exhibiting small contributions from the perturbative order LO 2 . 10On the other hand, we reckon that this approximation may miss non-negligible contributions for (multi-)boson production in association with more than one jet, for instance Z + 3j studied in Ref. [63], since LO 2 contribution is not negligible in the tails of the distributions.
We discuss now the case of the H events, denoted in the following as e H .As we mentioned multiple times we wish to ensure that condition (2.12) is valid if at least one of the r kl invariants, defined as r kl ≡ (p k + p l ) 2 (see also Eq. (A.1)), is such that |r kl | < M 2 W . Actually, since this is a prescription for matching δ EWSL preserving the EWSL accuracy in both the n and n + 1 final states, a more general condition is preferable, where c H→S should be chosen of O(1) and can be varied around the default value c H→S = 1 in order to test the dependence on it.Analogously to Eq. (3.12), an e H event is associated to a process of the form where we understood again that for a 2 → n process n = n + 2. In the FKS language, e H is one of the r ∈ R n+1 partonic processes with n + 1 particles in the final state. 11If one considers all possible j → kl branchings (g → gg, g → q q, and q → qg, but also Q → Qg, with Q being a quark with non-zero mass) a list of new processes with n particles in the final state is obtained by removing any possible (k, l) pair and substituting it with j in its place.In doing so, one also obtains for a given e H the (k, l) pairs that are associated with a soft and/or a collinear singularity, the set of FKS pairs P FKS (e H ) [102], and at the same time the S events e (k,l) S associated to processes of the form: where we remove the φ i k and φ i l particles and we add the particle φ i j in the position n.In the case of hadronic collisions, for a given production process the partonic process (3.14) that can be associated to an event e H is not unique.Thus, for each event e H , the FKS pairs P FKS (e H ) and the processes associated to the e (k,l) S events can be different.It is important to note that the mapping from the n+1 momenta {p} of e H to the n momenta {p ′ } of e (k,l) S is not uniquely defined.We have used for this purpose techniques for the momentum reshuffling analogous to those of Ref. [103]. 12In particular, given an FKS pair (k, l), one particle (which we associate to the index l) always belongs to the final state, while the other one (k) can be either initial or final.
• If k is a final-state particle, the pair k, l is first replaced by its mother particle j with p j = p k + p l .Then, the energy component of p j is changed in order to fulfil the mass-shell condition.Finally, the initial-state momenta are changed so that the new set of momenta satisfies momentum conservation.
• If k belongs to the initial state, then l is removed from the event.The remaining finalstate particles are boosted to their total-momentum centre-of-mass frame and, again, the initial-state momenta are changed in order to satisfy momentum conservation.
We now specify the quantity δ EWSL

(H)
. Given all the possible (k, l) pairs of external states for an event e H we define r min, abs ≡ min( with ( k, l) being the pair returning the smallest value for |r kl |.Introducing the quantity the process e H with the n + 1 kinematic.Otherwise, the Sudakov contribution is calculated via the same quantity evaluated instead for the underlying Born configuration, and the associated n-body kinematics, that is obtained via the replacement of the FKS pair giving the smallest invariant with its parent particle.In the unlikely (but possible) situation that the smallest invariant is given by a pair not corresponding to a QCD branching, Eq. (3.19) says that the EWSL are calculated directly for the process e H with the n + 1 kinematics, but within the DP algorithm the quantity r kl is replaced by M 2 W times the sign of r kl .Events of this kind with W are very unlikely, since they are not associated to any divergence.However, this replacement ensures that events are not reweighted via artificially large Sudakov contributions in a region where the approximation is not supposed to work.
The last point concerning the replacement is actually more general and we implemented it in MadGraph5_aMC@NLO as a safety feature as Indeed the EWSL approximation should be used only when invariants are large, but we want to prevent that artificially large correction may arise from simulations performed for processes with |r kl | < M 2 W already at LO QCD accuracy.The replacement is performed not only for simulations in the NLO QCD ⊗ EWSL + PS approximation, but also for the NLO QCD+EWSL one.
In conclusion, we can summarise the description of the NLO QCD ⊗ EWSL + PS predictions as follows: 1. Events are generated at NLO QCD accuracy via the MC@NLO matching scheme.
2. Events are reweighted [67] via the prescription in Eqs.(2.13) and (2.14) using the SDK weak scheme and neglecting the second term in the r.h.s. of Eq. (A.10).This leads to NLO QCD ⊗ EWSL accuracy in the MC@NLO matching scheme.
3. Events are showered via a parton shower including QED effects (possibly after heavy particles are decayed using external tools).

Numerical results
In this section we present numerical results for phenomenologically relevant physical distributions from two different production processes at hadron colliders: the top-quark pair and Higgs boson associated production (t tH), and the associated production of three Z gauge bosons (ZZZ).
Here we focus on the presentation of the EWSL-based predictions, and we do not perform any comparison with the exact NLO EW corrections.Still, we remind the reader that EWSL are an approximation for the NLO EW corrections.Hence, the quality of such an approximation should always be checked, at the differential level, before relying on it for phenomenological predictions.
For both processes considered here we show inclusive results (without any cut applied) as well as applying the following cuts: where X, Y is any of the particle in the final state at the Born level, p T is the transverse momentum and ∆R ≡ (∆ϕ) 2 + (∆η) 2 with ∆ϕ being the azimuthal angle between X and Y and ∆η is the difference between their pseudorapidities.The results have been obtained by generating events with MadGraph5_aMC@NLO and using Pythia8 as parton shower.Hadronisation is disabled in the parton shower.Input parameters are defined in the G µ scheme for what concerns EW renormalisation: and the top quark and Higgs boson masses are set to We employed the NNPDF4.0parton-distribution-functions [108], with NNLO evolution and α S (M Z ) = 0.118.The renormalisation and factorisation scales have been set equal to H T /2, where H T is the scalar sum or the transverse energies of all the particles in the final state, before showering the event.For what concerns the shower starting scale, we use the default setting in MadGraph5_aMC@NLO, which, for processes with massive or colourless particles in the final state, is a value proportional to H T [10,101].Jets are clustered via the anti-k T algorithm [109] as implemented in FastJet [110], with R = 0.4 and are required to have p T,min = 10 GeV.

t tH
We discuss results for the following differential distributions from t tH production in protonproton collisions at 13 TeV: p T (H) in Fig. 2, p T (t) in Fig. 3, the invariant mass of the top-quark pair m(t t) in Fig. 4, and p T (j 1 ) in Fig. 5, where j 1 is the hardest jet.In each of the Figs.2-5 we show inclusive results without cuts in the left plot and those with cuts (4.1) in the right one.The layout and the rationale of each plot is the following.In the main panel we show the central values for predictions with the three different accuracies: • NLO QCD + PS (grey), • NLO QCD ⊗ EWSL + PS (blue), • NLO QCD+EWSL + PS (red).
In the first inset we display the scale uncertainty band for the same three predictions normalised to the central value of NLO QCD +PS, where the uncertainty has been evaluated by independently varying the renormalisation and factorisation scale by a factor of two up and down (the usual 9point scale variation).The band for NLO QCD + PS is obviously centred around one and we show only the upper and lower bounds as dashed grey lines.In the last inset we show the ratio of the NLO QCD ⊗EWSL+PS and NLO QCD +PS predictions for different values of the parameter c H→S , introduced in Eq. (3.13) and entering Eq. (3.18).We remind the reader that c H→S parametrises how the condition (2.12) is implemented, 13 as can be seen from the aforementioned equations. 13Roughly speaking this means δ EWSL if an invariant connected to a soft/collinear limit is smaller in absolute value than c H→S M 2 W .The default case c H→S = 1 corresponds to the ratio of the blue and grey lines in the main panel.Before commenting the individual distributions we give some general considerations.We first note that the plots without cuts display both distributions for p T 's well below the M W scale (Figs.2-5) and the m(t t) distribution at the threshold (Fig. 4).It is well known that the Sudakov approximation is not expected to hold for these cases and also that additional effects such as Sommerfeld enhancements can be present [132].However, as already mentioned before, the focus of this paper is not to present phenomenological predictions for t tH (or ZZZ in Sec.4.2), but rather to document the features and the technical implementation of the NLO QCD ⊗ EWSL + PS approximation in MadGraph5_aMC@NLO.We leave discussions and comparisons with NLO EW corrections and/or data for future detailed studies.We here simply reckon that owing to the reweighting procedure, further classes of EW effects can be naturally incorporated via the redefinitions of the quantities δ EWSL two predictions are almost equal for the entire spectrum.The non-vanishing (still at the level of a very few percent) and flat effect at small p T is due to the fact that all Born invariants are at least ∼ m t in absolute value and s ≥ (2m t + m H ) 2 > M 2 W .For this reason, not only at the threshold are the EWSL non-vanishing: in this range they are positive, and are therefore clearly dominated by the single logarithms.Since (logarithms of the) invariants have very small variations for p T (t) ranging from 0 to ∼ 100 GeV, these corrections are quite flat.The fact that NLO QCD ⊗ EWSL + PS and NLO QCD+EWSL + PS are very close is due to the fact that EWSL are relatively small and the QCD K-factor, NLO QCD + PS/LO QCD + PS is very close to one and flat.Indeed, in first approximation, one expects that in other words, EW and QCD corrections combined in the multiplicative approach.An exception is the case of m(t t) in Fig. 4, where we observe the opposite trend: EWSL are not flat and NLO QCD+EWSL + PS and NLO QCD ⊗ EWSL + PS predictions are different.We verified that indeed the QCD K-factor increases at large values of m(t t).Similarly to what has been observed and discussed in Refs.[73,37,133] for the case of t t production in a similar context, one can notice that the scale uncertainty band is smaller in the case of NLO QCD ⊗ EWSL + PS predictions than in the case of NLO QCD+EWSL + PS.This is not a surprise since in the former setup the EWSL multiply NLO QCD corrections, while in the latter one they multiply only the LO QCD component, which has a LO dependence on the factorisation and renormalisation scales.This improvement is clearly not present in the p T (j 1 ) distribution (Fig. 5), since this distribution is not even present at LO QCD at fixed order.On the other hand, it is interesting to notice that for large values of p T (j 1 ), where NLO QCD + PS is dominated by hard matrix-element contributions and not PS effects, NLO QCD+EWSL + PS converges to exactly NLO QCD + PS at variance with QCD corrections are not dominated by hard emissions, such as e.g. in the case of the pT (t) spectrum, the difference between the r.h.s. and l.h.s. of Eq. (4.4) amounts to typically 1-5% of the NLOQCD + PS prediction.NLO QCD ⊗ EWSL + PS, which includes EWSL corrections also to the first real emission from hard matrix-element.If the cuts (4.1) are applied, we can clearly see that the impact of the EWSL increases and similarly the discrepancy between the NLO QCD+EWSL + PS and NLO QCD ⊗ EWSL + PS predictions increases.In the case of p T (j 1 ) distribution (Fig. 5) we see a flat contribution and a change at very large values for p T (j 1 ), where the simulation starts to be dominated by hard matrix-element contributions.One should notice, in particular for this distribution but also for all the remaining ones, that the NLO QCD ⊗ EWSL + PS is completely insensitive to the value of c H→S if varied by a factor of two up and down w.r.t. the reference value c H→S = 1.
Finally we comment on several checks that we performed and are not directly documented in the plots.The t tH cross section at LO involves contributions not only of order α 2 S α but also of order α S α 2 (and α 3 ).Therefore t tH is one of those process that potentially may involve EWSL contributions from the quantity δ QCD LA (see Eqs. (A.9)-(A.12))that we do not include (see Eq. (3.11)).We have explicitly verified at fixed order that the impact of this term is at most at the permille level and therefore can be safely ignored.We have also verified the effect of not implementing Eq. (3.20), finding no difference for the distributions, although real emission H events with an invariant smaller in absolute value than M 2 W and not associated to any QCD splitting have been identified (see Eq. (3.19)).Similarly to what has been observed in Ref. [64], the inclusion of the logarithms of the form as in Eq. (A.8) has a non-negligible impact.

ZZZ
In this section, first we discuss results for ZZZ production that are analogous to those of Sec.4.1 for ttH (Sec.4.2.1).Then, in Sec.4.2.2, we discuss the case with decays of the Z bosons, scrutinising the impact of the inclusion of QED effects in PS simulations.

Stable Z
We discuss results for the following differential distributions from ZZZ production in protonproton collisions at 13 TeV: the transverse momentum of the hardest Z boson p T (Z 1 ) in Fig. 6 and of the softest one p T (Z 3 ) in Fig. 7, the invariant mass of the two hardest Z bosons m(Z 1 Z 2 ) in Fig. 8 and of p T (j 1 ) in Fig. 9.The layout of the plots in Figs.6-9 is the same of those in Figs.2-5.
Before discussing the specific distributions we focus on the main differences with the t tH distributions in Sec.4.1.As already observed in the literature, in the case of multi-boson production EWSL are very large (see e.g.Refs.[134,16,66,135]) and much larger than in the case of t tH production.Indeed, in Figs.6-9 we observe a much larger impact of EWSL than in Figs.2-5.Moreover, at variance with t tH production, the LO cross section of ZZZ production does not depend on α S .Therefore scale uncertainties are smaller than in the t tH production.Still, NLO QCD corrections can be very large in multi-boson production (see e.g.[136][137][138]), therefore also the difference between NLO QCD ⊗EWSL+PS and NLO QCD+EW +PS is enhanced, especially when the cuts defined in (4.1) are applied.We anticipate that, similarly to the case of t tH production, we do not observe a dependence on the value of c H→S .
In the case of the p T (Z 1 ) distribution (Fig. 6) we can clearly see how EWSL are sizeable, especially in the right plot where cuts are present.The same argument applies to the discrepancy between NLO QCD ⊗ EWSL + PS and NLO QCD+EW + PS.We reckon absolute rates are smaller than what could be reasonably measured at LHC, even after the HL program, however, in the tail of the distribution, the effect of EWSL should be not only taken into account but also resummed. 15The same considerations are valid and even stronger for the case of the p T (Z 3 ) distribution in Fig. 7.
Turning to the m(Z 1 Z 2 ) distribution (Fig. 8), we can see that the previous discussion for p T distributions is valid also here, although with much weaker effects in the case without cuts 15 While the leading EWSL of the form α k log 2k (s/M 2 W ) can be in principle resummed via a simple exponentiation, the next-to-leading case α k log 2k−1 (s/M 2 W ) is not straightforward, as can be seen in Refs.[59,139] and further references therein.
(plot on the left).When we consider the case with cuts defined in (4.1) applied, we observe an additional effect for low values of m(Z 1 Z 2 ).For m(Z 1 Z 2 ) ≲ 700 GeV, the NLO QCD + PS simulation is dominated by the real emission contribution from hard matrix elements.First, this explains why the red line in the first inset of the plot on the right converges to one for small m(Z 1 Z 2 ) value, similarly to what has been discussed for Fig. 5. Second, since the dominant contribution originates from ZZZ +1 jet, it has a much larger dependence on the renormalisation scale.Indeed, the LO cross section for that process is of order α 3 α S .For this reason, the scaleuncertainty bands are larger for m(Z 1 Z 2 ) ≲ 700 GeV.
The p T (j 1 ) distribution (Fig. 9) shows the same features discussed already for the corresponding distributions in t tH production (Fig. 5).However, here the effects are magnified and clearly show the superiority of NLO QCD ⊗ EWSL + PS approximation w.r.t. the NLO QCD+EW + PS one.

Z → e + e − decays
In this section we consider the case in which Z bosons are decayed.Via MadSpin [68], the three Z bosons are decayed into e + e − pairs, after including the EWSL in the event as done in the previous section. 16We also consider the effects induced by EWSL in the NLO QCD ⊗ EWSL + PS or NLO QCD+EW + PS approximations and the impact of QED effects in PS shower simulations.As already mentioned, we denote with PS when these effects are taken into account and with PS QED when they are ignored.The process that we consider is therefore pp → e + e − e + e − e + e − , in the Breit-Wigner approximation emerging from ZZZ production.For the sake of simplicity, when the QED shower is enabled, only the photon emissions off charged particles (quarks and leptons) is allowed and the photon splitting into charged fermions is disabled.Therefore, we always have exactly six electrons/positrons (we will generally call them electrons in the following) in the final state, plus a number of photons, as well as quarks and gluon from the parton shower. 17We perform electron-photon recombination as follows: 1. Final-state electrons and photons are clustered with FastJet into jets using the Cambridge-Aachen algorithm (hence the clustering is purely geometric) [140,141], with distance parameter R = 0.1 and asking for a minimum jet transverse momentum p T > 25 GeV.
2. Out of the jets returned, we consider as leptonic jets those jets for which the sum of the charge of the constituent particles is different from zero.
3. The event is kept if exactly six leptonic jets are found, and otherwise it is discarded.
4. Negatively charged leptonic jets are sorted according to their transverse momentum and they will be dubbed as e − i , i = 1, 2, 3, where e − 1 is the hardest one. 16When performing the decay, MadSpin includes a smearing of the invariant mass of the decay products according to a Breit-Wigner distribution, for which we employ the width ΓZ = 2.49877 GeV. 17 We recall that in the simulations presented here hadronisation is turned off.

5.
To each of the e − i , the corresponding positively-charged leptonic jet e + i is assigned such that the quantity is minimised (this means that, in general, e + i will not be sorted according to their transverse momentum).
In Fig. 10 we display the following distributions: the transverse momentum of the hardest leptonjet, p T (e − 1 ) (top-left), the invariant mass of the e + 1 e − 1 pair, m(e + 1 e − 1 ), in a range close to M Z (topright), and the same distributions for the softest lepton jets, p T (e − 3 ) (bottom-left) and m(e + 3 e − 3 ) (bottom-right).
The usage of MadSpin allows for the correct reconstruction of the tree-level spin correlations, which would be lost if the decay were performed via a general PS simulator in which spin correlations are not preserved.However, the correlation of the helicity-dependent EWSL with the helicities and angular distributions of the decay products within this framework is not correctly addressed.Therefore, the impact of EWSL on observables that are sensitive to spin correlations cannot be correctly taken into account.However, we stress that this limitation is not due to the reweighting procedure for the inclusion of the EWSL presented in this work, but rather to MadSpin itself, and it affects also the case of the NLO QCD + PS predictions.
The layout of the plots in Fig. 10 is different w.r.t.those in Sec.4.2.1.In particular, the main panel is the same as in the plots of Figs.2-9, while the insets are different.In the first inset we show both the NLO QCD ⊗ EWSL + PS (blue) and NLO QCD+EW + PS (red) predictions normalised to the NLO QCD +PS one.This is similar to Figs. 2-9, but without uncertainty bands.In the second inset instead we show for both predictions, with the same colour convention of the first inset, the ratio of the PS and PS QED case.We have decided to omit from the plots the dependence on c H→S , since also in this case no visible effects have been observed.
As a general comment, the observables displayed in Fig. 10 show that both EWSL effects (first insets) and the QED effects (second insets) cannot be neglected.Also, their relative importance strongly depends on the considered observable.In the transverse-momentum distributions, the growth of the EW corrections in the Sudakov approximation is manifest, reaching −20% w.r.t.NLO QCD for p T (e − 1 ) ≃ 500 GeV and already for p T (e − 3 ) ≃ 100 GeV.Around these values, the benefits of the NLO QCD ⊗ EWSL + PS prediction over the NLO QCD+EW + PS one start to be visible, with the former displaying negative corrections about 5-10% larger than the former.For these observables, QED effects in the PS lead to ∼ −10% corrections, which are mostly related to a reduction of probability in passing the selection cuts specified in the previous bullet points.Indeed, the photon radiation reduces the leptonic-jet energy, leading also to a very mild shape distortion.We stress that, such effects strongly depends on the recombination details and that their relative impact does not depend on the employed approximation for the EWSL (NLO QCD ⊗EWSL+PS or NLO QCD+EW +PS), since the two classes of EW corrections factorise.
Turning to the invariant-mass distributions, the situation is somehow the opposite.On the one hand, as expected, the enhancement due to QED corrections in the region m(e + 1 e − 1 ) < M Z , and especially m(e + 3 e − 3 ) < M Z , is sizeable.QED corrections are important not only in the first  bins of the invariant-mass plots of Fig. 10, where they easily exceed +100% effects, 18 but also in the bins around m(e + 1 e − 1 ) = M Z and m(e + 3 e − 3 ) = M Z .These bins are the most relevant for the correct simulation of the signal region of ZZZ production.On the other hand, the impact of EWSL is flat and amounts to just −5%. 19However, should we have asked for e.g.cuts on the transverse-momenta of the e + i , e − i pairs, the EWSL effect would have been larger, as expected.For these invariant-mass distributions, no significant difference is visible among the two predictions NLO QCD ⊗ EWSL + PS or NLO QCD+EW + PS.All in all, these plots in Fig. 10 show that both effects due to QED radiation and to the EWSL have to be included, not only for achieving percent precision as shown, e.g., in Ref. [144], but also because they can in general give effects of order 10% or more.We re-stress here that using the SDK weak scheme for the evaluation of EWSL, the QED effects can be incorporated without any problem of double-counting.

Conclusion and Outlook
In this work we have presented the automation in MadGraph5_aMC@NLO of combined NLO+PS accuracy in QCD (NLO QCD + PS), Electroweak Sudakov Logarithms (EWSL) corrections, and QED final-state radiation (FSR) for event generation of SM processes at hadron colliders.Our strategy consists in the reweighting of NLO QCD + PS events for taking into account the EWSL contribution.FSR effects of QED are simulated directly via the PS.
We do not only reweight the LO contribution from the hard process, but also the QCD one-loop virtual contribution as well as the contribution from the first QCD real emission, the latter taking into account the different kinematic and external states for the evaluation of the EWSL.Moreover, since we have adopted the so-called SDK weak scheme [64] for the evaluation of EWSL, FSR or in general QED effects can be included in the PS simulation avoiding their double-counting.We have denoted the accuracy of our simulation as NLO QCD ⊗ EWSL + PS and motivated via theoretical arguments and numerical results its superiority to an approach where only the LO is reweighted with EWSL.In particular, we have shown and discussed results for physical distributions from pp → t tH production and pp → ZZZ production.Concerning the latter, we have also considered the case with Z → e + e − decays and stressed how neither EWSL nor QED FSR effects can be neglected.
In this paper we have focussed on the technical implementation of the NLO QCD ⊗EWSL+PS accuracy in MadGraph5_aMC@NLO, and also its automation and validation.The approach we have adopted is actually completely general and could be in principle extended to other tools 18 For a correct description of this process in the phase-space region |m(e + i e − i ) − MZ | ≫ ΓZ a full simulation for the complete process pp → e + e − e + e − e + e − would be necessary, or at least the contributions from γ → e + e − splittings and their interferences with Z → e + e − off-shell decays should be taken into account.This is the reason why, in the first bins of the invariant-mass plots of Fig. 10, QED effects appear even larger than what has been documented in, e.g., Refs.[142,143,77]. 19Since the invariant mass of the system is much larger than MW already at the threshold for the on-shell production, s ≥ (3MZ ) 2 , and the EW Casimir of the Z boson is particularly large, it is not surprising to observe non-vanishing EWSL also for Z bosons that are almost at rest.The value is quite flat because we calculate the EWSL for the on-shell production only.However, the variation of a few GeV for m(e + 1 e − 1 ) or m(e + 3 e − 3 ) has a negligible impact on the value of the invariants built with the momenta of the Z bosons, so we do not expect very large effects if one also takes into account the Z-boson off-shellness in the EWSL evaluation.
that use different matching schemes for NLO QCD + PS simulations.Indeed, since the approach is based on the reweighting, i.e. a step happening after the event generation, it does not rely on the strategy for the event generation itself.Moreover, since the evaluation of EWSL involves only tree-level matrix elements and compact analytical formulas, an advantage of the reweighting via the Sudakov approximation is the speed and especially the numerical stability of the results.
A natural follow up of our work is the extension of this technology to the matching and merging of NLO QCD + PS predictions with different jet multiplicities, namely, in the Mad-Graph5_aMC@NLO framework, the FxFx formalism [145], similarly to what has been done in Ref. [66].Also, an improvement of MadSpin in order to account for the information of correlation of the helicity-dependent EWSL with the helicities and angular distributions of the decay products would be beneficial for observables that are sensitive to spin correlations.
We have left phenomenological studies and comparisons with exact NLO EW accuracy predictions to dedicated works.Since EWSL are an approximation of the exact NLO EW corrections, there can be non-negligible effects at high energy that cannot be captured, such as photon-initiated contributions.These comparisons are therefore crucial before performing any phenomenological study.However, we remind the reader that a dedicated comparison of the EWSL automation in MadGraph5_aMC@NLO and exact NLO EW corrections at fixed-order has already been performed in Ref. [64].Moreover, with our approach, the reweighting factors for the NLO QCD + PS events can be augmented with further contributions on top of the EWSL.In other words, not only the percent-level mismatch with exact NLO EW corrections can be further decreased, but also higher-order EW effects or even BSM contributions can be taken into account.
for any individual helicity configuration that does not give a mass-suppressed amplitude in the high-energy limit, and for a generic SM partonic processes with on-shell external legs.First of all, the algorithm strictly relies on the assumption that all invariants are much larger than the gauge boson masses.Specifically, with k and l being two generic external particles with momenta p k and p l , The DP algorithm has been formulated for amplitudes with n arbitrary external particles, where all momenta p k are assumed as incoming.Processes are denoted as where the (anti-)particles φ i k are the components of the various multiplets φ of the SM.Moreover, contributions from longitudinal gauge-bosons are always evaluated via the Goldstone-boson equivalence theorem.If the Born matrix element for the process in (A.2) is written as ), which on the other hand involve different processes than the original one in (A.2), as can be seen in the indices.On top of that, the quantities δ i ′ where r kl denotes a generic kinematic invariant and M any of the masses of the SM heavy particles (M W , M H , m t and M Z ) or in the case of the photon the IR-regularisation scale Q, using the same notation of Ref. [64].Second, on the couplings of each external field φ i k to the gauge bosons V a and another field φ i ′ k , I a i k i ′ k (k), or associated quantities such as electroweak Casimir operators C ew , which involve the entire SU(2) × U(1) group.The exact expressions, as implemented in MadGraph5_aMC@NLO, can be found in Ref. [64] and are based on the results of Refs.[57,58].
As mentioned before, an important limitation of the DP algorithm is that, for a given process, at least one helicity configuration must not be mass suppressed, i.e., the amplitude should scale on top of LO 2 .Thus, in LA, the contribution from one-loop corrections to the quantity Σ NLO 2 , denoted as Σ virt NLO 2 can be written of the form (A.10) The quantity δ EW LA is calculated via the DP as summarised in Sec.A.1 and in particular Eqs.(A.3) and (A.4).With M 0 being the amplitude that once squared leads to Σ LO 1 , Strictly speaking, Eq. (A.10) is valid only under two simple assumptions: Q 2 = s and ∆ s→r kl (r kl , M 2 ) = 0. Otherwise, also the information on the colour-linked matrix elements would be necessary.A simple expression for δ QCD LA under the two aforementioned assumptions has been provided Ref. [64] and is reported later, in Eq. (A.12).However, as we will discuss in the following and as has already been mentioned in Sec.3.2, δ QCD LA is not so relevant as δ EW LA for the physical observables and processes considered in this work.This is ultimately the reason why the previous two assumptions are irrelevant in view of the EWSL approximation in the context of physical cross sections, which we are going to describe in the following.
Similarly to the case of δM, the quantity δ EW LA is IR-divergent and therefore non-physical.This is not a surprise, being an approximation of virtual EW corrections, which involve contributions from massless photons.Using the same notation as in Ref. [64], this is the scheme denoted as SDK, meaning the DP algorithm for the calculation of the amplitudes as in Eq. (A.4), 21and in the case of squared matrix elements contributing to the virtual component of NLO EW corrections as in Eq. (A.10).The SDK scheme is therefore a very good approximation at high energies of loop amplitude and virtual contributions, but it is not directly suitable in the case of phenomenological predictions.
In order to obtain predictions in LA that can be used for physical cross sections, the approach that has been employed often in literature is what has been denoted in Ref. [64] as SDK 0 .We stress here again that the SDK 0 is an approach mostly driven by simplicity.Indeed, it bypasses the problem of IR finiteness by simply removing some QED logarithms involving M W and the IR scale, but those logarithms arise from the conventions used in Refs.[57,58] and not from physical argument.
In Ref. [64], the SDK weak scheme has been precisely designed in order to solve this problem.The main idea behind it is that in sufficiently inclusive observables the Sudakov logarithms of QED and IR origin in the virtual contributions cancel against their real counterpart.In fact, the SDK weak scheme consists in a purely weak version of the SDK approach where almost all contributions of QED IR origin are removed 22 , while those of QED and UV origin are retained.
and (3.4).In fact, what we have denoted as NLO QCD+EWSL + PS in Eq. (2.8) corresponds to generating events by setting δ EWSL (H) = 0 and multiplying only the term B in (3.4) by δ EWSL (S) In the following we describe in detail how the δ EWSL (S) and δ EWSL (H) functions are implemented, starting with the case of δ EWSL (S)

5 Figure 2 :
Figure 2: Differential distributions for p T (H) in t tH production at 13 TeV.Left: no cuts applied.Right: cuts as defined in (4.1) applied.

5 Figure 3 :
Figure 3: Differential distributions for p T (t) in t tH production at 13 TeV.Left: no cuts applied.Right: cuts as defined in (4.1) applied.

5 Figure 4 :
Figure 4: Differential distributions for m(t t) in t tH production at 13 TeV.Left: no cuts applied.Right: cuts as defined in (4.1) applied.

5 Figure 5 :
Figure 5: Differential distributions for p T (j 1 ) in t tH production at 13 TeV.Left: no cuts applied.Right: cuts as defined in (4.1) applied.

5 Figure 6 :
Figure 6: Differential distributions for p T (Z 1 ) in ZZZ production at 13 TeV.Left: no cuts applied.Right: cuts as defined in (4.1) applied.

5 Figure 7 :
Figure 7: Differential distributions for p T (Z 3 ) in ZZZ production at 13 TeV.Left: no cuts applied.Right: cuts as defined in (4.1) applied.

5 Figure 9 :
Figure 9: Differential distributions for p T (j 1 ) in ZZZ production at 13 TeV.Left: no cuts applied.Right: cuts as defined in (4.1) applied.

Figure 10 :
Figure 10: Differential distributions for the transverse momentum (left) of the hardest (top) or softest (bottom) negatively-charged electron jet, and its invariant mass with the corresponding positively-charged one (right).

1 i 1
...i ′ n in depend only on two kinds of ingredients.First, on logarithms of the form L(|r kl |, M 2 ) ≡ α 4π log 2 |r kl | M 2 and l(|r kl |, M 2 ) ≡ α 4π log |r kl | M 2 , (A.5) In a few words, Eq. (3.18) says that if the smallest invariant is larger in absolute value than the M 2 W scale the Sudakov contribution δ EWSL in 0 (p 1 , ..., p n ) ,(A.3)in LA the O(α) corrections to M 0 , δM, can be written of the form 20δM i 1 ...in (p 1 , . .., p n ) = M Equation (A.4) is the main reason why the EWSL for physical cross sections can be computed in a much faster and stable way than the exact NLO EW corrections.Indeed Eq. (A.4) means that in LA the one-loop EW corrections can be written in terms of only tree-level amplitudes (M i ′ 1 ...i ′ n 0 (p 1 , . .., p n )δ i ′ 1 i 1 ...i ′ n in .(A.4)