One-loop electroweak Sudakov logarithms: a revisitation and automation

In this work we revisit the algorithm of Denner and Pozzorini for the calculation of one-loop electroweak Sudakov logarithms and we automate it in the MadGraph5_aMC@NLO framework. We adapt the formulas for modern calculations, keeping light-quarks and photons strictly massless and dealing with infrared divergences via dimensional regularisation. We improve the approximation by taking into account additional logarithms that are angular dependent. We prove that an imaginary term has been previously omitted and we show that it cannot be in general neglected for $2\rightarrow n$ processes with $n>2$. We extend the algorithm to NLO EW corrections to squared matrix-elements that involve also QCD corrections on top of subleading LO terms. Furthermore, we discuss the usage of this algorithm for approximating physical observables and cross sections. We propose a new approach in which the QED component is consistently removed and we show how it can be superior to the commonly used approaches. The relevance of all the novelties introduced in this work is corroborated by numerical results obtained for several processes in a completely automated way. We thoroughly compare exact NLO EW corrections and their Sudakov approximations both at the amplitude level and for physical observables in high-energy hadronic collisions.


Introduction
After more than ten years of operation of the Large Hadron Collider (LHC), we have tremendously improved our knowledge of the fundamental interactions of elementary particles. Above all, the long-sought Higgs boson has been observed [1,2] and especially its properties have been studied in detail: they have been found to be compatible with those predicted by the Standard Model (SM) [3]. In general, no clear and unambiguous sign of beyond-the-SM (BSM) physics has been found at colliders. In parallel, at the LHC the SM itself has been stress-tested in all its sectors: e.g., electroweak (EW) interactions, QCD dynamics and flavour physics. However, the BSM search programme at colliders is only at its initial phase. At the LHC, 20 times more data will be collected in the next years, large part of it during the High-Luminosity (HL) runs [4][5][6][7][8][9]. Moreover, several options are possible for future colliders (see e.g. Ref. [10] for a recent review), involving collisions at higher energies between a pair of protons or leptons (both electrons/positrons and muons).
The success of this ambitious research programme is interconnected with the availability of precise and reliable SM predictions. A plethora of new calculations and techniques have already appeared in the literature for improving both SM and BSM predictions. QCD radiative corrections at fixed order, going from Next-to-Leading-Order (NLO) to Next-to-NLO (NNLO) or even Next-to-NNLO (N 3 LO), have become available and techniques for the resummation of large logarithms appearing at fixed order have also been improved. On the other hand, an enormous effort has been done for the calculations of NLO QCD and also NLO EW corrections for processes with high-multiplicity final states. In particular, such corrections have been implemented in Monte Carlo generators and they have been even automated [11][12][13][14][15][16][17], at different levels in the different frameworks, using various one-loop matrix-element providers [18][19][20][21][22][23][24].
A particular feature of EW corrections are the so-called "Sudakov enhancements" or "Sudakov logarithms" [25], which enhance O(α n ) fixed-order corrections, the so-called N n LO EW corrections, with terms of order ∼ −α n log k (s/M 2 W ) with k in the range 1 ≤ k ≤ 2n. In high-energy collisions, these logarithms involve two separate ranges of energy scales: the W -boson mass M W and the centre-of-mass energy √ s. Most importantly, at variance with QCD, EW Sudakov logarithms do not cancel in IR-safe physical observables [26][27][28][29], which typically are not inclusive on the additional emission of neither W nor Z bosons. Therefore, they induce large and negative corrections, especially in the tails of the distributions. With the future runs of the LHC, and especially with future colliders, higher energies will be probed at higher precision and therefore a reliable evaluation of such corrections and its automation in modern Monte Carlo generators is necessary.
It is an established fact (see e.g. section 17 of Ref. [90]) that EW Sudakov logarithms from both virtual corrections and real emissions are sizeable at high energies (a few TeV's) for several processes and, especially for the former class, resummation is necessary in order to achieve precise predictions or even just sensible results. Indeed, fixed-order corrections, precisely due to the EW Sudakov logarithms, can approach the size of −100% of the LO. In other words, in order to reach the percent accuracy, not only the exact O(α) (NLO EW) corrections have to be calculated, but the Sudakov-enhanced component must be identified and resummed. Although the computation of exact NLO EW corrections is technically more involved than the case of its virtual Sudakov-logarithm subclass, only the former has been (fully) automated and implemented in Monte-Carlo's by different collaborations [11][12][13][14][15][16][17]. Recently, the pioneering algorithm of Denner and Pozzorini [39,43], which allows to calculate both single and double one-loop virtual EW Sudakov logarithms at O(α), has been automated for the first time [91] in the Sherpa framework [92]. Previously, only specific (classes of) processes had been considered, as done e.g. in Ref. [79] within AlpGen [93].
In this work we automate the algorithm of Denner and Pozzorini [39,43] in the Mad-Graph5 aMC@NLO framework [16,94], which already allows for the fully-automated calculation of NLO QCD and EW corrections and more in general Complete-NLO predictions [16,17]. Thus, with MadGraph5 aMC@NLO, it is now possible not only to calculate in a completely automated approach NLO EW corrections, but also their subcomponent that is typically dominant: the double and single virtual Sudakov logarithms. On the one hand, this work opens up the possibility of further building an automated framework for resumming Sudakov EW logarithms and match the result to NLO EW calculations in the MadGraph5 aMC@NLO framework. On the other hand, since virtual Sudakov logarithms are the dominant component of EW corrections and their evaluation is much faster and stable of the exact O(α) result, this implementation leads also to a very good and fast approximation of NLO EW corrections at high energy. Before implementing the algorithm of Denner and Pozzorini [39,43], however, we have revisited it. This work therefore does not only describe the technical steps underlying the implementation in MadGraph5 aMC@NLO, but it also provides an algorithm that is based on the one of Denner and Pozzorini and introduces relevant novelties w.r.t. it. First, we have reframed the algorithm by setting the mass of the photon and light-fermion masses exactly to zero, regularising IR divergences by mean of Dimensional Regularisation (DR), as in modern NLO EW calculations and in general Monte Carlo implementations. Second, we have identified an imaginary term that was omitted in Ref. [39], which however cannot be in general neglected for 2 → n processes with n > 2. Third, we have modified part of the expressions in order to take into account additional angular dependences, without assuming that all the invariants are of the same size of s. Fourth, since (virtual) NLO EW corrections originate from both "genuine" EW corrections on top of the dominant LO contributions and QCD corrections on top of the subdominant ones, we provide additional terms for taking into account the logarithmic dependence not only of the former, as done in the original work of Denner and Pozzorini, but also of the latter. All these modifications concern the approximation of the matrix elements or, more precisely, the interference of tree-level and renormalised one-loop amplitudes leading to the ultra-violet (UV) finite virtual corrections. The fourth point mentioned before was not considered in the pioneering work of Ref. [39] precisely because therein the focus was on the one-loop amplitudes and not their interferences with tree-level amplitudes, which is what yields the virtual contribution to NLO EW corrections. The effects of all these modifications and their validation are then showcased presenting numerical results obtained via the implementation in Mad-Graph5 aMC@NLO. In a fully automated way, for several different processes, we compare exact results for virtual contributions obtained via MadLoop [18], one of the modules of Besides purely virtual contributions, which are unphysical and IR divergent, we also show comparisons between the Sudakov approximation and the exact NLO EW corrections for production processes in proton-proton collisions, taking into account also the necessary additional terms to achieve IR finiteness: real emission of massless particles (photons, quarks and gluons) and PDF counter-terms. We show how, for a large class of processes and IR-finite observables, at variance with the case of virtual amplitudes, the exclusion of the contribution of photons from the Denner and Pozzorini algorithm [39,43] leads to a better approximation of NLO EW corrections. We describe in detail how the algorithm has to be altered in order to exclude the QED component (photons) and keep the purely weak one (W and Z bosons).
The article is organised as follows. In Sec. 2 we revisit the work of Denner and Pozzorini [39,43] at the pure amplitude level. We set the notation, using as much as possible the same one of Ref. [39], and we introduce three of the main novelties mentioned before: the formulation with strictly massless light-fermions and photons, the missing imaginary term, and additional terms that better take into account angular dependences and differences among the invariants. In Sec. 3 we move to the virtual NLO EW level, considering the interference of tree-level and one-loop amplitudes. Therein, we provide the additional term for taking into account logarithms of QCD origin. In Sec. 4 we present a modification of the algorithm such that the contribution of photons and gluons is excluded. We discuss how this can be a superior approach for approximating physical (IR-safe) cross sections at NLO EW accuracy. In Sec. 5 we describe the technical steps for the implementation of the algorithm in the MadGraph5 aMC@NLO framework. We explain in detail the procedure for the generation of all the additional amplitudes that are necessary for the evaluation of the Sudakov logarithms and especially for the evaluation of the amplitudes themselves. This procedure requires new features of the code, such as the possibility to evaluate the interference between amplitudes with different external legs or numerical derivatives of the matrix elements. In Sec. 6 we provide numerical results comparing NLO EW virtual contributions obtained via MadLoop and via the new implementation of the Sudakov approximation. We show the relevance of the novelties introduced w.r.t. Ref. [39] and at the same time we validate the new implementation. In Sec. 7 we repeat a similar comparison for physical observables from selected hadronic production processes in highenergy collisions. We discuss the numerical results and show how the exclusion of the contribution of photons leads to better approximations of the NLO EW results. Both in Sec. 6 and Sec. 7 results are obtained in a completely automated way via the new implementation in the MadGraph5 aMC@NLO framework. We give our conclusions and outlook in Sec. 8.

The Denner and Pozzorini algorithm revisited
We start by revisiting the pioneering work of Denner and Pozzorini [39,43], which provides an algorithm for calculating one-loop EW double-logarithmic (DL) and single-logarithmic (SL) corrections of the form for any individual helicity configuration of a generic SM partonic processes. We introduce three novel features w.r.t. the algorithm of Ref. [39], which we will denoted in the following as DP algorithm: 1. In all formulas light-fermions, photons and gluons are strictly massless, as in modern higher-order calculations. In other words, IR divergences are regularised via DR, introducing a IR-regularisation scale Q.
2. We correct the expressions for the Subleading Soft-Collinear terms (see Secs. 2.2 and 2.4), taking into account imaginary contributions that are relevant for 2 → n processes with n > 2.
3. We keep track of terms that are proportional to the (squared) logarithm of the ratio between two invariants that can be built via two different pairs of momenta among those of the external particles. Namely, using the notation that will be introduced later in this section, the terms proportional to log r kl r k l and log 2 r kl In this section, we will try to use as much as possible the notation of Ref. [39], where the DP algorithm has been formulated. In this way the reader can easily detect the differences introduced in our work. Moreover, we will try to avoid unnecessary repetitions of the content of Ref. [39], but we will also define all the quantities that are entering the actual implementation in MadGraph5 aMC@NLO, which is then discussed in Sec. 5.

Range of validity and conventions
The DP algorithm strictly relies on the assumption that processes with on-shell external legs are considered and, especially, that all invariants are much larger than the gauge-boson masses. In other words, if k and l are two generic external particles with momenta p k and p l respectively, then It is interesting to notice that the condition (2.2) still allows for kinematic configurations with r kl r k l M 2 W , where the quantities r kl and r k l represent a generic pair of the many possible invariants that one can build with two external momenta. However, since the required formal accuracy consists of the DL and SL in (2.1), although logarithms of the form α 4π log 2 r kl r k l and α 4π log r kl r k l , (2.3) are present at O(α) and can be non-negligible for configurations with r kl r k l M 2 W , they are not taken into account. In other words, the algorithm assumes the regime (2.2), but large logarithms may be anyway not captured unless the condition is satisfied for any possible pair of r kl and r k l invariants. In fact, condition (2.4) is quite unrealistic for actual calculations in collider physics, since cross sections are dominated precisely by regions where one or more r kl invariants tend to be much smaller than s ≡ r 12 M 2 W . Indeed, the r kl invariants are related to those entering the propagators. Even if cuts are devised in order to maximise any possible value of r kl for a given s, the fulfilment of condition (2.4) is strictly impossible. For instance, if (2.2) is valid, one has that min(r kl /s) < 0.5 for a 2 → 2 process. This bound is even tighter and tighter for a generic 2 → n process with n growing. 1 It is worth to remind the reader an important limitation of the DP algorithm. For a given process, at least one helicity configuration of the matrix element must not be mass suppressed, i.e., it must not vanish in the limit M 2 W /s → 0. 2 Indeed, such an assumption is one of the hypotheses under which the algorithm has been derived. On the other hand, most of the processes do satisfy this hypothesis, having at least one helicity configuration that is not mass suppressed 3 . Moreover, thanks to the condition (2.2), helicity configurations that 1 The determination of the configuration where the smallest invariant is maximal in a 2 → n process is related to the determination of the largest possible value for the minimum angle between any two of the n final-state momenta. This is the typical example of a mathematical problem that it is easy to define and with a solution that is far from trivial. See for example http://neilsloane.com/packings/index.html#I.
2 An equivalent formulation of this condition is that the scaling of the matrix element M with the centre- are not mass suppressed are by definition also dominant in the kinematic regime considered. The condition (2.2) also implies that processes including unstable particles and their decays cannot be treated in this approximation if physical observables are dominated by resonant configurations. Rather, the process without decays should be first considered and the decays should be then taken into account only after applying the DP algorithm. Being aware of all the possible limitations given by the conditions (2.2) and (2.4), we describe the DP algorithm and some modifications we have introduced in order to achieve the formal leading and subleading logarithmic accuracy (LA), i.e., taking into account only enhanced DL and SL terms of the form (2.1), for one-loop EW virtual corrections to any SM amplitudes, in DR and therefore with possibly massless particles. The problems related to the validity of condition (2.4) will be also addressed, giving a pragmatic solution.
The starting point of the DP algorithm is that since all the terms considered are logarithmic, they can be expressed via the quantities where r kl can be any of the invariants 4 and M any of the masses among M W , M H , m t and M Z , depending on the associated Feynman diagrams. Moreover, in the case of massless particles, the regularisation of the divergences will lead to logarithms of the form (2.5) where M → Q and Q is the IR-regularisation scale. The most important point, in order to understand the novelties introduced in this section, is that the DP algorithm splits twice the logarithms of the form in (2.5); both splittings are connected to the modifications of the DP algorithm that we present in this work. First, logarithms of the form in (2.5) are split into two classes: a symmetric and solely energy-dependent class, which is associated to the scales M W and √ s and parametrised by the quantities L(s) ≡ L(s, M 2 W ) and l(s) ≡ l(s, M 2 W ) , (2.6) and a remaining class of logarithms involving mass ratios and ratios of invariants. This splitting involves the imaginary component that we are going to introduce in the formulas and that is not present in Ref. [39]. It also involves the modifications that take care of the violation of condition (2.4). Second, while above the scale M W all one-loop EW contributions are treated in an unified approach, without separating purely QED from purely weak effects, below the M W scale only the QED component is present, involving logarithms between M W and the IR scale. In other words, for the contribution from QED loops M W works as a technical separator. Above M W we have for example (see eq. (2.19)) quantities parametrised via the electroweak Casimir operator C ew , which involves the entire SU(2) × U(1) group, while below M W we have only quantities that involve the charges Q k of the external particles.
The latter class of contributions is denoted by the apex "em", standing for electromagnetic, and in Ref. [39] it arises from the energy hierarchy M H , m t , M W , M Z m f =t λ, where λ is the mass of the photon. In this separation the logarithms l(M 2 W , M 2 Z ), l(m 2 t , M 2 W ), and l(M 2 H , M 2 W ), similarly to the quantities log (|r kl |/s) and log 2 (|r kl |/s), are neglected when they do not multiply the term l(s).
In Ref. [39] all the quantities denoted as electromagnetic ("em") depend on m f =t and λ, which here are considered exactly equal to zero, The consequence of (2.7) is that DR becomes necessary. In D = 4 − 2 dimensions, electromagnetic logarithms transform into 1/ poles plus finite terms and logarithms involving the IR-regularisation scale Q. In this context we are not interested in the structure of the IR poles, which for NLO EW corrections is discussed, e.g., in Refs. [16,17,24,95]; we are interested only in the logarithmic dependence of the finite part. This can be simply derived via the substitutions in the expressions of Ref. [39]. We want to comment on the dependence on Q, the IR-regularisation scale, which is introduced here and it is not present in the original DP algorithm. The derivation of the formulas in Ref. [39] depends on the assumption that µ 2 = s, but therein µ is the UVregularisation scale since all the IR divergences are regularised via m f =t and λ. However, similarly to this work, therein formulas have been derived assuming an on-shell renormalisation scheme, such as the α(M Z ) or G µ ones. With such a renormalisation scheme, no renormalisation-scale dependence is present for one-loop renormalised amplitudes, both if exactly calculated or using the LA. Therefore, the DP algorithm, although derived assuming a specific value for µ (µ 2 = s), returns results that do not depend on µ. The substitution in (2.8), that we perform due to the condition (2.7), does not depend on the condition (2.2). Moreover, it affects only the regularisation of IR divergences and does not concern the UV ones. Therefore, this substitution introduces the correct dependence on Q even if a common regulator for UV and IR divergences (Q = µ) is used. Exceptions are discussed in Sec. 3.
Before providing the expressions necessary for automating one-loop EW Sudakov logarithms, we introduce further conventions and notations according to Ref. [39]. Amplitudes are assumed with n arbitrary external particles and all momenta p k incoming. Needless to say, any 2 → n − 2 amplitude can be rewritten into a n → 0 amplitude via crossing symmetry. Processes are denoted as where the (anti)particles ϕ i k are the components of the various multiplets ϕ of the SM: • f κ σ andf κ σ : chiral fermions and antifermions, with chiralities κ = R, L and the isospin indices σ = ±, • V a = A, Z, W ± : gauge bosons transversely (T) or longitudinally (L) polarised. Neutral gauge bosons are also denoted as N = A, Z, • Φ i : the scalar doublet containing the Higgs particle H and the neutral and charged Goldstone bosons χ, φ ± .
An important technical point of the DP algorithm is that, since the high-energy limit is assumed, the Goldstone-boson equivalence theorem can be used. In fact, with this algorithm, contributions from longitudinal gauge-bosons are always evaluated via the Goldstone-boson equivalence theorem. We will return to this point in Sec. 5.1.
Following the same notation of Ref. [39], the couplings of each external field ϕ i k to the gauge bosons V a is denoted by ieI Va (ϕ), namely, ieI Va ϕ i ϕ i (ϕ) is the coupling corresponding to the V aφi ϕ i vertex, with all fields that are incoming. For simplicity, in the formulas the components ϕ i k are replaced by their indices i k , namely, I a i k i k (k). All the values and formulas for the quantities I a i k i k (k), as many other terms appearing in the next sections are reported in detail in the appendices of Ref. [39]. We do not repeat them here, but we want to warn the reader that the same exact conventions for Feynman rules have to be used in order obtain consistent results.
For any process denoted as in (2.9), the Born matrix element reads The O(α) corrections to M 0 in LA, δM, has the form Equation (2.11) means that the result can be written in a factorised form, but that involves Born amplitudes for different processes. The contributions to δM have different origins: The quantities δ LSC and δ SSC are respectively the leading and subleading soft-collinear logarithms. They both emerge from the DL, which in turn originate from the eikonal approximation of one-loop diagrams where gauge bosons are exchanged between external legs and are soft-collinear. The former represents the symmetric and solely energy-dependent class of logarithms, while the latter involves mass ratios and ratios of invariants. The quantity δ C consists of the collinear logarithms, originating from virtual collinear gauge bosons from external lines and field renormalisation constants. The logarithms resulting from parameter renormalisation, which can be determined by the running of the couplings, correspond to the term δ PR . In the case of longitudinally polarised bosons the equivalences  In the following subsections we provide the formulas entering the implementation in MadGraph5 aMC@NLO, which is described in Sec. 5. We will discuss in detail only the aspects concerning the differences w.r.t. Ref. [39].

Logarithm splittings
As already mentioned, the DL corrections come from loop diagrams with virtual gauge bosons V a = A, Z, W ± connecting two external legs. In particular, they originate from regions where the gauge boson is soft and collinear to one of the external legs. Their expressions can be derived by evaluating them in the eikonal approximation. In Ref. [39], DL have been in general identified as where the invariant r kl depends on the angle between the momenta p k and p l . Equation (2.14) precisely represents the first of the logarithm splittings that has been mentioned before. In the first line of (2.14) the quantity L(|r kl |, M 2 ) is split into L(s, M 2 ), which is symmetric and energy dependent, and other two terms, of which the second can be neglected in the approximation (2.4). Moving to the second line, the remaining terms are further rearranged such that if M = M W , the mass-ratio logarithm log M 2 W M 2 is kept only when multiplying l(s). The dots at the end stand for the terms that are dropped in the splitting of the logarithms. In Ref. [39], the first two terms in the second line of eq. (2.14) are identified as the leading soft-collinear (LSC) contribution, which as already mentioned is angular-independent and involves only the s/M 2 W ratio in the logarithms. The remaining term leads to the angular-dependent subleading soft-collinear (SSC) contribution.
When loop diagrams with virtual gauge bosons V a = A, Z, W ± connecting two external legs are evaluated in the eikonal approximation, the logarithmic dependence can be derived by the expansion of the C 0 function in the high-energy limit, namely condition (2.2). The expression can be found in Ref. [96]. If the gauge boson V with mass M is exchanged between the external particles φ k and φ l , the relevant quantity is, following the conventions of Ref. [96], where Θ is the Heaviside step function. It is then clear that rather than starting from L(|r kl |, M 2 ) as in eq. (2.14) the correct quantity to be taken into account is The difference is an imaginary component that involves a term proportional to l(s). For 2 → 2 processes, this is completely irrelevant and therefore all the results presented for specific processes in Ref. [39] are not affected by this additional term. Indeed, since 2 → 2 tree-level amplitudes are always real (as a consequence of the optical theorem), the imaginary part of the one-loop (or Sudakov-approximated) amplitude drops out when the real part of the loop-tree interference is considered. However, this is no longer the case starting from 2 → 3 processes, and indeed we do find that this imaginary part is not irrelevant. We therefore repeat the procedure of eq. (2.14) in order to identify how the impact of the term 2iπΘ(r kl ) translates into the DP algorithm. Moreover, we keep track of the terms that would be otherwise discarded assuming condition (2.4). Starting from (2.16) we obtain where we have dropped in the splitting of the logarithms only terms involving neither s nor r kl . 5 In the third line of eq. (2.17) there are terms that are relevant for the formal expansion in LA, i.e., the correct expression to be used instead of (2.14). The first two terms in the sum give the LSC logarithms, while the third one contributes to the SSC ones. On the contrary in the fourth line there are further terms that become relevant when s r kl M , i.e., departing from condition (2.4). Formally, they do not enter the LA so they cannot be identified neither as LSC nor as SSC. On the other hand, since they depend on r kl , we will take into account their contribution in the expression of the SSC logarithms (Sec. 2.4). For this reason we have denoted them in eq. (2.17) as SSC s→r kl .
As we will discuss in more detail in Sec. 6.2, even taking into account the SSC s→r kl contribution, the full control of logarithms involving the ratios of |r kl | invariants and s cannot be achieved via the DP algorithm. We will discuss the case of a specific process for which this limitation is manifest. On the other hand, several numerical results in Sec. 6.2 and Sec. 7 clearly show how the inclusion of the SSC s→r kl terms substantially improves the approximation of such class of logarithms and in turn of EW virtual one-loop corrections at high energy.

LSC: Leading soft-collinear contributions
The LSC logarithms can be rearranged as a single sum over the external legs,  where δ LSC i k i k (k) reads (2.19) In this case, besides the term L em (s, Q 2 , m 2 k ), the expression is the same as Ref. [39]. The expressions for the electroweak Casimir operator C ew i k i k (k), the squared Z-boson coupling (I Z (k)) 2 i k i k , and the charge Q 2 k for a generic particle k and a specific polarisation can be found in Ref. [39]. It is important to note that the first two quantities have indexes and can be non-diagonal. We will return to this point discussing the implementation in MadGraph5 aMC@NLO. Using DR the electromagnetic DL reads (2.21)

SSC: Subleading soft-collinear contributions
Unlike the LSC terms, the SSC ones remain a sum over pairs of external legs of the form This part is the one with the largest differences w.r.t. Ref. [39]. The exchange of soft neutral gauge bosons contributes with and charged gauge bosons yields The quantity ∆ s→r kl (r kl , M 2 ) is set equal to zero when the condition (2.4) is assumed and the LA is applied in a strict sense, as done in Ref. [39]. Taking instead into account the fact that s r kl M 2 , this quantity reads 25) and precisely corresponds to the SSC s→r kl logarithms of eq. (2.17). The quantities I A , I Z and I ± are the couplings with respectively the photon, the Z boson and the W ± boson, where we have omitted the indices i j i j . While I A is always diagonal in these indices, I Z can be non-diagonal and I ± (k) is always off diagonal. The impact of the new imaginary terms proportional to iπΘ(r kl ) on results obtained with the DP algorithm is directly connected to the aforementioned off-diagonal structures. Indeed the virtual contribution to NLO EW corrections involves terms of the form 2 (M 0 δM * ), where are real, like in 2 → 2 processes, the contributions of imaginary terms proportional to iπΘ(r kl ) vanish, otherwise they formally contribute. It is also interesting to note that with DR and massless photons, setting Q 2 = s the entire δ A,SSC contribution vanishes if we also set ∆ s→r kl (r kl , M 2 W ) = 0. This can be seen from the definition of δ A,SSC in eq. (2.23). This argument will also be recalled in Sec. 3.1, where the QCD contribution to NLO EW corrections to squared matrix-element is discussed.

C: Collinear and soft single logarithms
In this section we provide the results obtained in Ref. [39], adapting them for the case with massless light-fermions and photons. The formula for the collinear and soft single logarithms can be written as a sum over the external particles and polarisations, that depends on the external particle and polarisation ϕ i k . We provide the results in the following. The expressions for all the new terms introduced in the formulas can be found in Ref. [39].

Transverse charged gauge bosons W
The result is

Transverse neutral gauge bosons A,Z
The results for symmetric and antisymmetric parts are expressed in terms of the coefficients b ew N N of the β-function. The result is The quantity Z em AA in DR reads

Longitudinally polarised gauge bosons
By means of amplitudes involving Goldstone bosons, the complete collinear corrections (2.27) for longitudinal gauge bosons is

Higgs boson
The complete correction is

PR: Logarithms connected to the parameter renormalisation
The last ingredient is the logarithms related to the UV renormalisation. In Ref. [39] they have been identified via the formula where the quantities are related to the top-quark Yukawa coupling and to the scalar self coupling, respectively. All the δ's are the logarithmic part of the renormalisation counter-terms of the corresponding dimensionless quantities. In the δ's, regardless of the value of Q chosen for the regularising the IR divergences in the other contributions (LSC, SSC, C), the UV regularisation-scale µ must be set as µ 2 = s. Indeed, although renormalised amplitudes in an on-shell scheme do not depend on the value of the unphysical UV-regularisation scale µ, the DP algorithm has been derived assuming µ 2 = s. Therefore, in order to preserve the cancellation of the µ dependence related to the UV poles, in the logarithmic part of the UV counter-terms it is necessary that µ 2 = s.
Here, we rearrange the formula in eq. (2.38) for practical purposes related to the implementation in MadGraph5 aMC@NLO, discussed in Sec. 5, but the results are fully equivalent with those of Ref. [39]. In practice we rearrange it into It is worth to recall that the renormalisation of masses in propagators or in couplings with mass dimension is not relevant, because those contribute only to mass-suppressed amplitudes. The parameter n tad is a technical parameter that has the only purpose of keeping track of the appearances of the tadpole counter-term. In practice what we do is to modify Feynman rules for three-scalar and four-scalar vertices by rescaling their value by the parameter n tad , which is then set equal to one in the numerical evaluation. We use the following formulas: where the purely electromagnetic part reads In this work, all the results are presented by adopting the G µ scheme, where in the place of α the input parameter is G µ , which is related to α via the tree-level relation . This translates into the substitution in eq. (2.40).
The remaining terms are where on-shell renormalisation for the mass is assumed, and finally with the contribution from the tadpole renormalisation reading

Sudakov logarithms and NLO EW corrections
In the previous section we have revisited the DP algorithm, which allows the calculation of electroweak DL and SL in LA for virtual scattering amplitudes. On the other hand, for collider results and in general for the calculation of physical observables, the relevant quantities are amplitudes that are either squared or interfered among them. In particular in this work our final goal is the NLO EW corrections to LO cross sections. For any differential or inclusive cross section Σ, adopting the notation already used in Refs. [12,16,17,82,[97][98][99][100][101][102][103][104], the different contributions from the expansion in powers of α S and α can be denoted as: with k being process dependent and k ≥ 1. Each Σ LO i denotes a specific α n S α m perturbative order that can be present at LO, i.e., arising from tree-level diagrams only. On the contrary, each Σ NLO i denotes a specific NLO perturbative order to which the interferences between different classes of tree-level and one-loop diagrams can contribute. For a given process, the values of n and m vary for each 3) implies something that is very well known and, e.g., has been discussed in detail in Ref. [82]. If Σ NLO i involves EW corrections (i > 1) and it is not the term with the possibly highest α power at NLO (i < k + 1), then both QCD and EW loops on top of tree-level amplitudes can enter into the game. Even worse, this separation into "QCD loops" and "EW loops" is artificial and especially cannot be rigorously defined. Since one of the main features of our implementation of the DP algorithm is the possibility of comparing the DL and SL terms in LA against the exact result for NLO EW corrections, the contribution of such "QCD loops" cannot be ignored.
In LA, the contribution from one-loop corrections to the quantity Σ NLO i , denoted as Σ virt NLO i can be written in the form In the following section we provide the necessary ingredients for taking into account single and double-logarithmic enhanced contribution of QCD origin in the computation of Σ NLO 2 , what is typically dubbed in the literature as "NLO EW corrections". The case of the Complete-NLO, i.e. the complete set of Σ NLO i contributions, is left for future work.
In practice, what is discussed in this work is sufficient for both the case of Σ NLO 2 and Σ NLO k+1 , since the latter never receives contributions from "QCD loops", as can be seen from eq. (3.3).
The quantity δ EW LA is what is calculated via the DP algorithm revisited in Sec. 2 and summarised in eqs. (2.10)-(2.12). For the case of Σ NLO 2 , or equivalently the case Σ NLO k+1 , if M 0 is the amplitude that once squared leads to Σ LO 1 , or equivalently Σ LO k , then As we said, we leave the case of the Complete-NLO for future work. In that case also eq. (3.5) would receive modifications since a generic Σ LO i−1 term with 1 < i < k can itself arise from the interference of amplitudes factorising different powers of α S and α.

Contributions from QCD loops
or an interferenceM 0,1M * 0,2 of two amplitudesM 0,1 andM 0,2 with with j being an integer in the range 0 < j ≤ min(n − 1, m + 1). Similarly to eqs. (2.10)-(2.12), where starting from the amplitude M 0 the logarithmic-enhanced EW corrections are denoted as δM, we can denote the logarithmic-enhanced QCD corrections toM 0 as δM. This implies that and respectively In principle, following the same steps of Sec. 2, one could derive a general algorithm for obtaining δM starting form a genericM 0 . Indeed, besides the case of non-abelian gluon vertices, the DL and SL logarithms can be identified by looking at the purely QED part of the expressions of Sec. 2. In practice, this would lead to non-trivial terms involving colour-linked amplitudes, which are not per se problematic, but still avoidable via two simple assumptions on the value of Q 2 and ∆ s→r kl (r kl , M 2 ).
As already mentioned at the end of Sec. 2.4, if one sets Q 2 = s and ∆ s→r kl (r kl , M 2 ) = 0, as in the formal derivation of Ref. [39], then the SSC contribution from purely QED origin vanishes. Similarly, simplifications in the rest of the expressions of Sec. 2 happen. Setting Q 2 = s and ∆ s→r kl (r kl , M 2 ) = 0, these simplifications are present also for the case of QCD. Especially, the SSC contribution vanishes also in the case of QCD corrections. As we will see later in Sec. 4.1, these two assumptions are innocuous for what concerns δ QCD LA in LA when physical observables are considered.
With these two assumptions, where n t and n g are the number of top quarks and gluons in the external legs, respectively. The quantities L t (s), l α S (µ 2 ) and (δm t ) QCD are defined as (3.14) with C F = 4/3, and they have a very different origin, as explained in the following. The terms proportional to L t (s) can be obtained by performing the substitution  .28)). The reason why the top quark is special is that we are understanding the use of the five-flavour scheme. If other fermions f are treated as massive (m f = 0), the corresponding logarithms with t → f should be also taken into account. This is true also for the remaining contributions discussed in this section. Conversely, for all the other massless quarks, if one sets Q 2 = s and ∆ s→r kl (r kl , M 2 ) = 0, not only the SSC but also the LSC and C contributions to δ QCD LA vanish. The term proportional to l α S (s) can be derived from the diagonal C contribution for the photon by applying the substitution (3.15). These logarithms are the virtual counterpart of the quasi-collinear logarithms emerging from the g → tt splittings. The term proportional to l α S (µ 2 R ) has instead a different origin; it is connected to the MS renormalisation of α S . While the renormalisation of the EW sector can be performed without introducing a renormalisation-scale dependence, this is unavoidable in QCD. With five active flavours, the logarithmic-enhanced part of the α S counter-term reads where µ R is the renormalisation scale and the quantity β 0 = 11 − 2 3 n f is the leading term of the QCD β function in the SM (n f = 6). We are assuming Q 2 = s and it is reasonable to assume also µ 2 R ∼ s, which let us to ignore the term proportional to log µ 2 R Q 2 in eq. (3.11). While the contribution from α S -renormalisation to the LA can be expressed via algebraic formulas, this is not possible in the case of m t -renormalisation (or equivalently for any other quark that would be considered as massive), where the derivative δM 0 /δm t is entering the expression. 6 The formula in eq. (3.14) has been obtained in the on-shell scheme, consistently with the EW case.
Via eqs. (3.11) and (3.10) we can finally write a compact formula for δ QCD

Sudakov logarithms and physical cross sections
What has been discussed up to this point concerns the LA of one-loop "EW corrections" (Sec. 2) and one-loop "QCD corrections" (Sec. 3.1) to amplitudes and their combination for the LA of the virtual contribution to the perturbative orders Σ NLO i with i = 2 or i = k + 1 (eq. (3.4)), in the perturbative expansion of the cross section Σ. Both cases, amplitudes or virtual contributions, are unphysical and cannot be directly used for theoretical predictions of physical quantities. Since electromagnetic contributions are included (one-loop QED corrections), the DP algorithm leads to the LA of a quantity that is IR divergent and must be combined at least with the LA of the IR-divergent real-emission contributions. Alternatively, the DP algorithm can be slightly modified by excluding the QED contribution, which is the only one leading to unphysical quantities. While virtual QED SL and DL involve the unphysical quantity Q, the remaining contributions of purely weak origin involve the physical masses. Clearly, also these logarithms can be partially canceled by their real-emission counter part, the heavy-boson-radiation (HBR) of an extra W , H or Z boson, but these cancellations strongly depend on the specific set-up and the degree of inclusiveness of the observable considered, see e.g. Ref. [73]. In other words, while photon and gluon emissions and real radiation of light quarks are unavoidable contributions for obtaining IR-finite predictions, the HBR is not necessary for the sake of IR finiteness. The contribution of HBR may be also relevant for the LA of the entire Σ NLO i prediction, but this critically depends on the process and the set-up considered.
In the literature, e.g. in the recent work of Ref. [91] or in Ref. [79], this problem has been circumvented by dropping the contributions tagged as "em" in the DP algorithm, namely those involving the ratios M 2 W /λ 2 or M 2 W /m 2 f , in the original formulation, or equivalently the ratio M 2 W /Q 2 in this work. We believe this approach is artificial and based on a wrong interpretation of the role of M W in the DP algorithm. While the DL and SL induced by W and Z boson loops (the L(|r kl |, M 2 ) and l(|r kl |, M 2 ) terms with M = M W , M Z ) are physical, in the case of QED M W is only a technical separator used in order to split the logarithms according to the logic discussed in Sec. 2.2. Bypassing the problem of IR finiteness by simply removing the logarithms involving M W and the IR scale is therefore an approach mostly driven by simplicity.
We propose a more rigorous approach for avoiding IR sensitivity in the implementation of the DP algorithm for physical cross sections. We will denote this approach as SDK weak , where SDK is an abbreviation for Sudakov. The SDK weak approach is based on the idea of selecting only the DL and SL of purely weak origin, excluding the contributions of QED corrections. Actually, for a large class of processes this approach leads to predictions that are much closer to the exact NLO EW corrections than in the case of approaches based on the removal of all "em" terms, what we will denote from now on also as SDK 0 approach. Indeed, in sufficiently inclusive observables, most of the logarithms of QED origin cancel against their real-emission counterparts.
The DP algorithm has been formulated in Ref. [39] for one-loop amplitudes and generalised in Sec. 3 for their interference with tree-level amplitudes. From now on, we will also denote the latter case as simply the SDK approach. One should keep in mind that, at variance with the SDK 0 and SDK weak approaches, the SDK one leads to IR-divergent quantities, which approximate correctly the virtual contribution to NLO EW cross sections, but that cannot be used for physical observables. We describe in the following how expressions of Secs. 2 and 3 should be modified in order to adapt the DP algorithm, which has been formulated so far for the SDK approach, to the SDK weak approach.

SDK weak : Purely weak LA for cross sections
In general, when W bosons are not involved in a process, virtual EW corrections can be divided into a QED and a purely weak component in a gauge-invariant way. QED corrections consist of all loops involving QED interactions between fermions and photons, excluding the vacuum-polarisation diagrams. 7 The purely weak part consists of all the rest of contributions, including also the vacuum-polarisation diagrams and the renormalisation of the photon wave-function and of the fine-structure constant α.
In order to isolate purely weak effect in DL and SL, we exclude all contributions induced by fermions and photons interactions in all formulas of Sec. 2, with the exception of those from parameter renormalisation (PR). Moreover, we exclude also DL and SL related to photons interacting with W bosons as external legs. While the classification of the W -γ interaction as either a purely weak or QED effect is ambiguous, the identification of the terms in the expressions of Sec. 2 that originate from such interaction is unambiguous.
The purely weak version of the DP algorithm, SDK weak , can be obtained following these steps: 1. Calculate the δ PR in eq. (2.12) as in the standard SDK approach.
2. For each external particle ϕ i k in (2.9), set This step alone has the effect of eliminating all the terms tagged as "em", with the exception of δZ em AA . It also eliminates all the SSC terms and C terms that lead to SL originating from photons, with the exception of those related to transverse W bosons.
3. For each external particle ϕ i k in (2.9), perform the replacement with the value of Q 2 k before enforcing eq. (4.1). This, in combination with eq. (4.1), has the effect of eliminating the DL due to photons.
This has the effect of eliminating for the transverse W bosons the C terms that lead to SL originating from photons.

Set
and perform the replacement This has the effect of eliminating, for the photons, the C terms that lead to SL originating from light fermions.
6. Calculate the remaining terms in eq. (2.12) with the new redefinitions of steps 2-5.
We want to stress that, thank to the step 1, the redefinitions of steps 2-5 do not apply to all the PR contributions discussed in Sec. 2.6; for them any QED-like contribution is retained. We remind the reader that also in this context we assume the use of either the α(M Z ) or G µ -scheme, which both have an IR structure that is MS-like, namely, IR poles are not present in the α counter-term, δα. 8 This difference of treatment for the PR terms, besides the definition of purely weak and QED introduced before, can also be understood in a different way. Logarithms from PR are related to UV renormalisation and do not involve IR sensitivity, as can be seen from all the equations in Sec. 2.6, which do not depend on the infrared regulator Q. Therefore, in order to achieve physical predictions and eliminate the Q dependence, there is no need to exclude their components that are related to QED interactions. This is in contrast with the LSC, SSC and C contributions, where the QED contributions are Q-dependent (see eq. (2.19) and eq. (2.23) for respectively LSC and SSC contributions, and eq. (2.29), which defines the term l em entering most of the C contributions) and therefore IR-sensitive. Moreover, while LSC, SSC and C contributions have real-emission counterparts that, together with PDF counter-terms, lead to IR finiteness and the (partial) cancellation of Q dependence, this is not the case for PR contributions.
As already mentioned, in sufficiently inclusive observables, most of the logarithms of QED origin cancel against their real-emission counterparts. Thus, as we will show in Sec. 7, the purely weak version of the DP algorithm, the SDK weak approach, reproduces very well the logarithmic dependence of NLO EW predictions for differential distributions in several cases. For example, if we consider leptons or massive particles in the final state, it is easy to understand how the purely weak version of the DP can be superior when predictions for physical observables are considered.
If the real emission of photons is considered in the eikonal approximation and integrated up to the energy E γ = √ s/2, where E γ is the energy of the photon, large part of the QED logarithms cancel. For example, all the QED contributions from LSC and SSC terms are canceled exactly. There are two classes of collinear SL that are left uncanceled: those associated to final-state radiation and those associated to initial-state radiation. 9 The former class of SL can be eliminated by clustering photons with charged particles. In the case of massless charged particles, the clustering is anyway unavoidable for IR safety (e.g. the case of dressed leptons) and eliminates the Q dependence. In the case of massive charged particles, namely top quarks and W bosons, clustering is also very reasonable since in a realistic experimental set-up the separation of collinear real radiation from very boosted objects is not feasible and, from a theoretical point of view, it leads to larger uncertainties. The clustering eliminates the physical collinear SL of QED origin, which have the form l(s, M ), where M is the mass of the massive charged particle. The latter class of SL (initial-state) are those related to the PDFs, which therefore in an exact NLO EW calculation are subtracted by the corresponding PDF counter-terms.
Both for the case of initial-and final-state collinear SL, a special case is given by the photon. Since the top quark and the W boson are massive, the corresponding collinear SL associated to their contribution to δZ em AA are not canceled. Indeed, in the case of the initial state, no corresponding PDF counter-terms are present, since massive particles do not enter the PDF evolution. In the case of the final state, no γ → tt or γ → W + W − radiation is generated for a final state giving the same signature. 10 In order to be on the same ground, we also remove from eq. (3.11) the term proportional 9 We analytically verified this statement for the process e + e − → W + W − combining the results of the DP algorithm together with the eikonal approximation of the real emission of photons [105]. 10 We remind the reader that we are assuming in this work that the α(MZ ) or Gµ scheme is used. The use of α(0) and therefore the case of isolated photons in the final state is not considered in this context.
to L t (s). Indeed this is canceled by the clustering of real emissions of gluons and top quarks into recombined top-quarks. At this point it is also easy to understand why in Sec. 3.1 we have said that setting Q 2 = s and ∆ s→r kl (r kl , M 2 ) = 0 is innocuous for what concerns δ QCD LA in LA when physical observables are considered. If we lifted these two assumptions, we would obtain many more contributions to eq. (3.11), but they would be all canceled by the real-radiation counterparts, together with the PDF counter-terms. The presence of the term n g l α S (s) can also be seen now as the gluon QCD-counterpart of the photonic uncanceled SL that have been mentioned in the previous paragraph: being massive, top quarks do not enter in the PDF counter-terms and no g → tt is generated. We leave to future work the exploration of these effects from QCD corrections in numerical results for physical cross sections.
All in all, the modifications to eq. (3.11) that have to be implemented in the SDK weak approach in order approximate physical cross sections is simply: in eq. (3.11).
Finally, we want to stress that the real emission of photons from electrically charged particles and of gluons from coloured particles cannot be neglected for IR-safe observables, but at high energy even the HBR may be considered and taken into account in inclusive calculations. In that case, the SDK weak approach should be further modified. For instance, taking into account the Z radiation, also the contribution of the Z boson should be removed from the DP algorithm. We leave the exploration of this approach for future work.

Implementation in MadGraphaMC@NLO
The theoretical framework described in the previous sections has been implemented in MadGraph5 aMC@NLO [94], specifically in the part of the code that is deputed to the calculation of NLO QCD and EW corrections and more in general Complete-NLO predictions [16,17]. This allows a direct comparison of results in LA and at exact fixed-order, both at amplitude level and for physical observables.
We remind the reader that in MadGraph5 aMC@NLO the IR singularities are dealt with via the FKS method [106,107], automated for the first time in MadFKS [108,109]. In MadGraph5 aMC@NLO one-loop amplitudes can be evaluated via different types of integral-reduction techniques (the OPP method [110] or the Laurent-series expansion [111]) and techniques for tensor-integral reduction [112][113][114], all automated within the module MadLoop [18]. Moreover, the codes CutTools [115], Ninja [116,117] and Collier [118] are employed within MadLoop, which has been optimised by taking inspiration from Open-Loops [20] for the integrand evaluation.
As already possible in the code, NLO QCD and EW corrections can be invoked via the syntax [QCD] [QED], see Refs. [16,17] for more details. However, now the code allows also for the evaluation of virtual one-loop Sudakov logarithms by adding after the command generate or add process the flag --ewsudakov. As we have said, the code works for the moment for O(α) corrections to the Σ LO i contribution with i = 1 and i = k, according to eqs. (3.1) and (3.2). In order to implement the DP algorithm in MadGraph5 aMC@NLO, three main technical features had to be implemented: 1. The generation of all the amplitudes that are necessary for the computation of the DL and SL.
2. The evaluation of the amplitudes, especially the interferences of amplitudes involving different external legs.
3. The evaluation of the derivatives of the amplitudes, which enter the formulas concerning the PR terms.
In the following subsections we address each of the previous points.

Generation of the amplitudes
We start discussing the case of a generic partonic process and at the end we return to the case of proton-proton collisions. The formulas of Sec. 2, which are given for n → 0 processes, can be easily reframed in terms of more common 2 → n − 2 amplitudes via crossing symmetry, 3) with 0 ≤ k W ≤ n W and 0 ≤ k Z ≤ n Z is used for the following steps in the generation of the amplitudes. As discussed in Sec. 2, the formulas for the different contributions leading to DL and SL involve amplitudes with external particles that are different from the original ones in M 0 . In particular, starting from the process in (2.9) it is necessary to generate the amplitudes for all the processes with 1 ≤ k ≤ n that can be obtained applying the substitution ϕ i k → ϕ i k of the form: H ←→ χ .
With the symbol ←→ we understand that the substitution works in the two directions. Substitution (5.7) is necessary for the off-diagonal components of C ew entering the LSC terms and of b ew N N entering the C terms. Substitution (5.8) is necessary for the off-diagonal components of (I Z ) 2 entering the neutral SSC terms. Moreover it is necessary to generate also the amplitudes for the processes that can be obtained either applying two substitutions ϕ i k → ϕ i k and ϕ i l → ϕ i l of the form (5.8), again for the off-diagonal components of (I Z ) 2 in the neutral SSC terms, or two different ϕ i k → ϕ i k and ϕ i l → ϕ i l substitutions that together do not violate charge conservation, each one of them of the form:
For hadronic calculations (protons in the initial state, jets, etc.) different partonic subprocesses can contribute at the Born level. Therefore, the procedure described so far has to be separately repeated for each of them.

Evaluation of the amplitudes
The evaluation of the amplitudes follows the standard procedure of the MadGraph5 aMC-@NLO framework, which relies on the helicity routines supplied by Aloha [119]. Here, the additional complication consists in the evaluation of interferences of amplitudes that can have different particles in the respective initial and/or final state, as shown in eq. (2.26). As discussed in the previous section, there can be one or even two different external particles between the two interfering amplitudes. Consequently, without altering the external momenta, external particles cannot be in general on-shell in both amplitudes. In order to preserve the on-shell conditions of external legs, external momenta have to be modified for one of the two amplitudes that are interfered. We follow this approach, modifying the external momenta of the amplitude with different external states w.r.t. the Born one.
From a technical point of view, this approach is very similar to the momentumreshuffling techniques discussed in Ref. [120], in the context of the so-called "Simplified Treatment of Resonances" (STR) that are needed to perform, e.g., computations in supersymmetric theories. 11 In both cases, on-shell conditions are enforced by modifying part of the external momenta and the remaining ones (possibly a subset) are reshuffled in order to preserve momentum conservation. On the other hand, in this context, not only this procedure has to be applied at the amplitude level and not at the squared-amplitude level, but it is also intrinsically more articulated. While in the case of Ref. [120] only one on-shell condition is enforced and involves the invariant mass of two final-state particles, in the present case one or more on-shell conditions have to be enforced and involve the kinematic mass, p 2 0 − | p| 2 , of one or two individual external momenta, from the initial and/or final state. Our implementation is based on the one described in detail in Sec. 5.2 of Ref. [120], on which we base the following discussion. We use the same notation for describing the technical details.
The case of only one on-shell condition for a particle in the final state can be directly derived from Sec. 5.2 of Ref. [120]. Using the same notation, it can be summarised as: given a set of momenta k i , generate a new setk i where a given (final-state) particle, denoted as β, changes its mass from m to m β = M . The two particles labeled as δ and γ in Ref. [120] are irrelevant for our case. Among the infinite number of solutions, two options, dubbed as A and B in Ref. [120], have been considered. In the former the energy-momentum conservation is imposed by modifying all the other final-state particle momenta, while in the latter, which is the default option in the code, by changing the momenta of the initialstate particles. The case of two on-shell conditions for particles in the final state can be achieved by applying the procedure iteratively.
The case of β being a particle in the initial state was not relevant for Ref. [120] and we will briefly present it here. In this case, we have chosen to change the momentum of only the other initial-state particle, leaving the final-state ones untouched (k i = k i for i ≥ 3). This implies that the centre-of-mass energy s is conserved and therefore the procedure is very simple. We start with the original initial-state momenta k i , with i = 1, 2, where the one with i = β is going to have a new mass. Since s must be conserved and we want the new momentak i collinear to the beam pipe, in the partonic centre-of-mass frame one has to simply derive the new quantity |k β,z | = |k 3−β,z | enforcing momentum conservation and on-shell conditions. We conclude by commenting on the fact that, when masses are modified for both an initial-state and a final-state particle, the procedure can again be performed iteratively. However, when the default option B is used for the final-state case, since it assumes massless initial-state momenta, the case of a new mass in the final state should be considered first, and only afterwards one should consider the initial-state one.
Before moving to the next section we want to clarify that all this procedure would be unnecessary if an analytical calculation were performed and all the mass-suppressed term were discarded. This is on the other hand not possible in an automated approach. The procedure outlined here leads to a correct evaluation of all the terms that are not mass suppressed. Indeed, all the modifications of the momenta and subsequent reshuffling operations involve only scales connected to the mass of the SM particles. The differences in the kinematics before and after the procedure outlined in this section are mass suppressed themselves, leading to smaller and smaller effects when the energy inrceases. Ambiguities related to the choice of a specific reshuffling technique and to the order in which the reshuffling is performed are also mass suppressed.

Derivative of the amplitudes
As can be seen in eqs. (2.40) and (3.11), in order to compute the logarithmic contributions induced by the parameter renormalisation, the derivatives of the amplitudes w.r.t. part of the input parameters have to be calculated. Although one may in principle use dedicated Feynman rules, such as those used for generating UV counter-terms in an NLO computation, we have opted to calculate the derivatives via numerical methods. In other words, for each phase-space point, we evaluate the quantity where x is any of the variables for which the derivative has to be performed (M W , M Z , etc.),x is its numerical value when the amplitudes are evaluated and δx is a small value, which has been set to δx = 10 −5 for the results presented in this work. The same procedure is done forM in eq. (3.11).
We have checked that this procedure has a mild impact on the speed of the code and the choice δx = 10 −5 is excellent in terms of both stability and precision. The use of the numerical derivatives allows also to easily adapt the calculation of PR terms for possible BSM scenarios, where additional particles and couplings would be present. Moreover, at variance with what has been done in the recent automation [91] in the Sherpa framework, the SL from PR terms are calculated exactly at O(α) as all the other type of DL and SL logarithms, without including spurious terms from higher orders in the α expansion. This fact is crucial for the systematic comparisons we are going to carry out in Secs. 6 and 7 between NLO EW exact results and their Sudakov approximations.

Numerical Results: matrix-element level
In this section we present numerical results obtained via the revisitation and implementation of the DP algorithm in MadGraph5 aMC@NLO, which has been described in the previous sections. We focus here on the Sudakov approximation of one-loop amplitudes and in particular on their interferences with the corresponding tree-level ones; results for cross sections of processes that are relevant at colliders are discussed in Sec. 7. We compare exact results for the finite part of the virtual contribution to the NLO EW corrections with their Sudakov approximation, what is denoted as the "SDK approach" following the notation introduced in Sec. 4. After having specified the input parameters in Sec. 6.1, we start in Sec. 6.2 by discussing the effect of the SSC s→r kl terms in eq. (2.17), which are relevant when the condition (2.4) is not satisfied. Next, in Sec. 6.3 we show the numerical relevance of the terms proportional to 2iπΘ(r kl ), discussed in Sec. 2.2. Then, in Sec. 6.4 we show the relevance of the QCD corrections to the subleading LO (see Sec. 3.1) in order to compare NLO EW corrections with their Sudakov approximation.
Throughout this section, we will consider only the finite part of the virtual contribution to the NLO EW corrections. It is worth recalling that, in general, the virtual contribution is IR divergent and non-physical by itself. Since we regularise IR divergences in DR, the finite part depends on the IR-regularisation scale Q, which we set here always as Q 2 = s.

Input parameters
The results presented in this section are obtained using the following input parameters: All the other SM particles are treated as massless and all the decay widths are set equal to zero. Consistently with the input parameters, EW interactions are renormalised in the G µ -scheme and masses and wave-functions in the on-shell scheme. QCD interactions are renormalised in the MS-scheme, with the renormalisation-group running at two loops and the renormalisation scale µ R set to µ 2 R = s.

Impact of SSC s→r kl terms
In this section we show numerical results for 2 → 2 partonic processes both varying the value of s and the angle θ between the first and third particle, which in turn parametrises the value of t. We select representative processes for which the relevant plots are displayed in Figs. 1 and 2. In both figures, the plots of each column refer to the same partonic process and the upper plots show the dependence of several quantities on the center-of-mass energy √ s, while the lower plots show their θ dependence. In the following, we describe the layout of the plots and how they should be interpreted.
In the first panel we show the value of the LO squared matrix-element, separately for each leading-helicity configuration and possibly their sum if there is more than one. In order to improve the readability of the legends in the plots, therein we display not only the helicity of any external particle, but also a conventional number associated to the ordering of the helicity configurations within MadGraph5 aMC@NLO. Conventionally, leading-helicity configurations have been identified as those with a value, for their squared amplitudes, that is at least 10 −3 times the one of the dominant helicity configuration. The main purpose of this conventional choice is to probe and select via a numerical method all the helicity configurations that are not mass suppressed. 12 In the first inset we show the 12 If at least one helicity configuration is not mass suppressed and it is the dominant, the ratio of its squared amplitude and the one of another helicity configuration that is not mass suppressed asymptotically converges to a positive constant at high energies. Therefore, leading helicities can be present over the entire  ratio between the O(α) virtual corrections and the LO in different approximations. We display as separate dots the exact results obtained via MadLoop (Virt) for selected values of s, while as lines 13 the LA approximation of Sudakov logarithms that are obtained via the new implementation of the modified DP algorithm described in this work. Dashed lines refer to the pure LA (SSC s→r kl terms not included), denoted in the plots as "SDK, s → r kl OFF", while the solid lines to the case in which SSC s→r kl terms are taken into account, denoted in the plots as "SDK, s → r kl ON". As expected, the values of the ratio over LO for s range if and only if they are not mass suppressed and this ratio is larger than 10 3 . In other words, if an helicity configuration is mass suppressed, it is for sure not tagged as leading, while if it is not mass suppressed can be not tagged as leading, but it means its contribution is at less than per-mill level of the dominant-helicity squared amplitude. 13  both dots and lines are negative and grow in absolute value for large values of s. A correct implementation and evaluation of the LA of Sudakov logarithms implies that the differences between each line and the dots converge to a constant value for s → ∞. Indeed, since all the mass-suppressed terms of O(α) corrections go to zero for large s, the terms that survive are either logarithmic enhanced, those that have to be exactly captured by the LA (lines), or constant for t/s fixed. We therefore separately display the interpolation of the difference between the dots and the solid line (second inset) and between the dots and the dashed line (third inset). These quantities are denoted as (Virt-SDK)/LO in the plots. The layout of the lower plots of Figs. 1 and 2 is very similar to the one of the upper plots, however, in this case the x-axis refers to the angle θ between the first and third particle, which in turn parametrises the value of t, in the range 10 −2 θ π/2. We have fixed the value of s to √ s = 10 TeV for all lower plots. In order to produce the upper plots, the scan in √ s with t/s fixed, we have performed the following procedure. We start by generating the momenta for a phase-space point with √ s = 10 3 GeV and t/s = −1/20 for the specific process considered. Then, we iteratively repeat the following steps for increasing the value of √ s by keeping fixed the t/s ratio within an error of the order of permille. First, we rescale the trimomenta of the outgoing particles by a common factor. Second, we impose on-shell conditions for the outgoing particles in order to obtain their energies. Finally, we impose momentum conservation for determining the momenta of the initial state. In this way, we can generate several phasespace points by scanning the √ s range and keeping the ratio t/s very stable. Each one of the phase-space points obtained is then used as input for evaluating the exact virtual NLO EW corrections of O(α) as well the LA with and without the inclusion of the SSC s→r kl terms. The SDK, s → r kl ON and the SDK, s → r kl OFF lines are the interpolation of these LA results.
As can be seen in both Figs. 1 and 2, all the second and third insets of upper plots show perfectly horizontal lines for large values of s, for each individual helicity configuration. We have shown here only representative processes, but we did not see any exception in all cases that we have checked. This is a clear sign of a correct implementation of the LA of Sudakov logarithms.
In order to rigorously check the last statement, we have fitted the quantities (Virt-SDK)/LO via a function of the form with the method of least squares. While the coefficient B has been found in general of the order of few percents for the plots shown here, the quantity A is in general of the order of 10 −4 and compatible with 0 due to the associated statistical error, 14 therefore supporting our previous statement about the correct implementation of the LA of Sudakov logarithms. If we consider Fig. 1, comparing the second and third inset of the upper plots we can also appreciate how SSC s→r kl terms further reduce the gap between the LA and the exact 14 We remind the reader that statistical errors also include effects induced by the numerical method that is used for performing the derivatives, which is discussed in Sec. 5.3, as well as by possible instabilities of the evaluation of exact virtual amplitudes.
value of the virtual. Being |t|/s = 1/20, these terms are constant when varying s, but still non-negligible. In other words, their inclusion preserves the LA and further improves the approximation of the exact result at high energy. In the lower plots, the impact of the SSC s→r kl terms can be better understood. By varying θ, the difference between the strict LA and the exact result is not expected to be a constant times the LO, because the condition (2.4) is assumed. This can be observed in the second and third insets. However, although in both insets lines are not perfectly horizontal, one can notice how the SSC s→r kl terms substantially reduce the slope. In other words, as can be also seen in the first inset, they lead to a better approximation of the exact results also at small angles. Clearly, at very small angles power-suppressed terms may not be negligible, since the M 2 W /|t| ratio increases and this in general leads to non-horizontal lines even with the inclusion of the SSC s→r kl terms. These effects are enhanced for the left and central plots of Fig. 2 and even more for the right one, the LO 3 ∝ α 2 of the dd → dd process, which deserves some more comments.
By looking at the right-upper plot, it is clear that the LA works well and both the second and third insets show perfectly horizontal lines, however, by looking the lower plots it is manifest that both with or without the inclusion of SSC s→r kl terms the θ dependence cannot be correctly captured, although better approximated in the former case. First, it is important to note, as can also be seen in the main panel of the bottom-right plot, that the LO cross section diverges for θ → 0, indeed it is proportional to 1/t 2 . Actually, if condition (2.2) is not assumed, not only terms proportional to 1/t 2 (the photon t-channel propagator) are present, but also terms proportional to 1/(t 2 −M 2 Z ) (the Z-boson t-channel propagator) appear. Approaching the regime with |t| = (1−cos(θ))s/2 M 2 Z , which implies θ →∼ 0.02 for √ s = 10 TeV, condition (2.2) is not valid and therefore even including SSC s→r kl terms a good approximation is not expected. 15 Second, the SSC s→r kl correctly takes into account the value of the different r kl invariants entering eq. (2.16), but this equation derives from eq. (2.15), which by itself does not depend on ratios of different invariants. This means, for instance for 2 → 2 processes, that additional corrections that involve the (s/t) ratio and have a functional form different from the terms included in the DP algorithm cannot be recovered, even with the inclusion of SSC s→r kl terms. In other words, while the correctness of LA is guaranteed and SSC s→r kl terms further improve the approximation keeping track of the correct dependence on the different r kl invariants in eq. (2.16), this does not mean that the full dependence is retained. In order to do that, the information on the internal structure of the Feynman diagrams would be necessary. Indeed, the starting point in the derivation of the SSC s→r kl terms is the C 0 function in (2.15), which is associated to simply the masses and the invariant mass of two external particles involved in the process. However, already with 2 → 2 processes, D 0 functions can appear in virtual corrections, involving also at high energies more than one invariant and leading to additional terms 15 We have noticed an additional interesting behaviour in the θ dependence, which is not shown here in the plot. When even smaller angles are considered, |t| M 2 Z , while the SDK, s → r kl OFF predictions depart even more from the exact result, the SDK, s → r kl ON approximation considerably improves. This is due to the fact that, for the LO3 ∝ α 2 of the dd → dd process, the QED contribution is dominant in this regime and furthermore it does not involve neither MZ nor MW . when the condition (2.4) is not satisfied.

Impact of the imaginary component
As explained in Sec. 2.2, in the original work of Ref. [39] an imaginary component has been omitted in the formulas. On the other hand, this component affects results only for 2 → n processes with n > 2. In this section we show numerical results about this aspect, for 2 → n partonic processes with n = 3, 4. 16 Again, we select representative processes for which the relevant plots are displayed in Fig. 3. Each plot shows the dependence on s of several quantities, and the layout is very similar to the one of the upper plots of Figs. 1 and 2. Here, the LA always includes the SSC s→r kl terms, 17 but we distinguish the case in which the terms proportional to iπΘ(r kl ) in eqs. (2.23)-(2.25) are excluded, as in the original DP algorithm in Ref. [39], or retained. The former are displayed as dashed lines (iπΘ(r kl ) OFF) and the latter as solid lines (iπΘ(r kl ) ON). For each leading helicity configuration, we also show in the second and third inset the difference between the LA and the exact result both normalised to the LO, respectively with and without taking into account the imaginary component.
In order to produce the plots, scanning in √ s, we have performed a procedure similar to the one explained in the previous section for the upper plots in Figs. 1 and 2. The only difference here is the starting point. For 2 → n partonic processes with n > 2, besides s, there is more than only one independent kinematic invariant that can be built with the external momenta. In order to avoid pathological configurations with an |r kl | M 2 W , we randomly generate the first set of external momenta setting √ s = 10 4 GeV and requiring We remind the reader, as already explained in footnote 1, that eq. (6.3) is a condition that can be satisfied for 2 → 3 or 2 → 4 processes, but not in general for 2 → n, for which this lower bound has to be lowered more and more increasing the value of n, further departing from the condition of eq. (2.4). Looking at Fig. 3, it is manifest how the case including terms proportional to iπΘ(r kl ) correctly catches the LA, while the other one does not; perfectly horizontal lines are present in the second inset, while in the third inset a dependence on s is clearly visible. For some of the processes considered, such as dd → Zdd, this dependence seems to cancel out for the sum over the different helicity configurations. In large part this is correct, but a small dependence is still present and it is simply not visible from the plot. We in general see this feature also for individual helicity configurations, namely the iπΘ(r kl ) is often formally relevant but sometimes the numerical effect is very small. For other processes, such as e + e − → e + e − µ + µ − or ud → Zud, even for the helicity-summed result the lack of the terms proportional to iπΘ(r kl ) leads to sizeable numerical effects.
In order to provide a more quantitative statement, we list in Tab. 1 the results of the fit of (Virt-SDK)/LO for each leading-helicity configuration (and their sum) of the process dd → Zdd. We have used again the method of least squares and the functional form of eq. (6.2). As can be seen in the third column of Tab. 1, all helicities exhibit a non-vanishing slope when the terms proportional to iπΘ(r kl ) are turned off. Notably, as anticipated before, this happens also for the sum over the helicities, which for this particular process and kinematic configuration (condition (6.3)) leads to a cumulative error of 2.6% in the LA for every factor of 10 in increase of the energy. The error is process dependent and can also be larger, as can be seen in the bottom-left plot of Fig. 3 for the e + e − → e + e − µ + µ − process. Finally, given Fig. 3, we would like to stress how well the LA of Sudakov terms can work in the high-energy regime. All these processes receive corrections of the order of −200% and the difference between the LA and the exact result is always (well) below the 10%.

Impact of the corrections of QCD origin
As explained in Sec. 3, NLO EW corrections can receive contribution of both EW and QCD origin (see discussion therein for details). In this section we show numerical results on this aspect, for 2 → n partonic processes with n = 2, 3, 4. We again select representative processes for which the relevant plots are displayed in Figs. 4. Each plot shows the dependence on s of several quantities, and the layout is very similar to the one of the plots of Fig. 3 in the previous section. Phase-space points have also been generated following the same procedure described in the previous section and according to the condition (6.3). As in Sec. 6.3 the LA always includes the SSC s→r kl terms in δ EW LA , however here we distinguish the cases in which the term Σ LO i δ QCD LA in eq. 3.4 is excluded or retained. The former are displayed as dashed lines (SDK, δ QCD LA OFF) and the latter as solid lines (SDK, δ QCD LA ON). For each leading helicity configuration, we also show in the second and third inset the difference between the LA and the exact result both normalised to the LO, respectively with and without taking into account the δ QCD LA component. As can be seen in Fig. 4, perfectly horizontal lines are displayed only in the second inset. 18 For each plot, lines in the third inset have a slope, which largely depends on the process considered. The slope is directly connected to the term Σ LO i δ QCD LA , therefore its size depends on both the size of δ QCD LA (see eq. (3.17)) and, since in the plots shown in Figs. 4 we have i = 2, on the LO 2 /LO 1 ratio. 18 Some lines do not span the entire √ s-range. These are related to helicity configurations that are actually mass suppressed but still contribute more than 10 −3 of the dominant helicity configuration for sufficiently "low" energies. See also the definition of leading helicity configurations at the beginning of Sec. 6.2 and in footnote 12. Considering simple 2 → 2 processes, one can see how different is the impact of Σ LO 2 δ QCD LA in the top-left and top-center plots; in the case of bb → tt both δ QCD LA and the LO 2 /LO 1 ratio are larger. The top-right plot refers to the process uū → ttgh, the simplest process for which all the terms of eq. (3.17) are non-vanishing. The lower plots refer to different partonic processes entering the process pp → tttt. As already discussed in Ref. [100] a large part of the NLO EW corrections are of QCD origin, and this can be observed also in the lower plots of Fig. 4.
Given the large number of leading helicity configurations for processes with tttt in the final state, we did not list them in the legend of the corresponding plots of Fig. 4. On the other hand, in Tab. 2 we provide, as an example, the results of the fit of (Virt-SDK)/LO summed (9.9 ± 8.4) · 10 −4 (−2.7 ± 0.4) · 10 −2 (−8.9 ± 1.2) · 10 −1 (3.1 ± 0.6) · 10 0 Table 2: Result of the fit of the quantity (Virt-SDK)/LO using the method of least squares and the function (6.2) for the representative process gg → tttt. The case including(excluding) the δ QCD LA contribution corresponds to the quantities shown in the second(third) inset of the lower-central plot of Fig. 4.
for all leading-helicity configurations (and their sum) for the process gg → tttt. We have used again the method of least squares and the function (6.2) as done for Tab. 1, but in this case we have considered two different scenarios: the inclusion or the exclusion of the δ QCD LA contribution. By comparing the numbers in the first and third column of Tab. 2, one can further see how the δ QCD LA contribution is essential in flattening the (Virt-SDK)/LO curve.

Numerical Results: differential cross sections
In this section we consider IR-safe collider observables and we compare exact NLO EW predictions and their LA (the approximation given by Sudakov logarithms) obtained via the DP algorithm as framed in Sec. 4, i.e., the SDK weak approach. We show differential distributions for proton-proton collisions at 100 TeV, which is one of the possible experimental set-ups that has been considered as an option for a high-energy future colliders [90,[125][126][127]. In this regime, LA is expected to be very efficient in capturing NLO EW effects. By explicit examples, we show how both the inclusion of the SSC s→r kl terms in eq. (2.17), which are relevant when the condition (2.4) is not satisfied, and the usage of the purely weak LA described in Sec. 4, the SDK weak approach, has to be in general preferred for predictions of physical observables. Indeed, these features not only improve the LA, but they are also instrumental in order to capture the correct logarithmic dependence.
As already said, we consider cross-section predictions for differential distributions of IRsafe observables. Thus, bare leptons, which are treated as massless, have to be recombined with photons into dressed leptons. A dressed lepton is here obtained by recombining a bare lepton with any photon γ that satisfies the condition where ∆R( , γ) ≡ (∆η( , γ)) 2 + (∆φ( , γ)) 2 , and ∆η( , γ) and ∆φ( , γ) are the differences of the bare-lepton and photon pseudo-rapidities and azimuthal angles, respectively. 19 In this context, however, we recombine with photons any electrically charged particle, including top quarks and W bosons, which are massive. This choice is not due to IR safety, but rather to the fact that very energetic massive particles are typically identified as tagged jets, namely a jet which contains the considered particle. The recombination of photons and heavy charged particles is inspired by precisely this procedure, and similarly the condition (7.1) leading to a large cone. In one case we will also consider the aforementioned tagged jets, where the clustering is performed not only with photons but with all the particles and afterwards tagging the considered heavy particle among the jet constituents. In that case, we will use the anti-k T algorithm [128] as implemented in FastJet [129], with R = 0.4. We remind the reader that, as already explained in Sec. 4.1, photon recombination has an effect on the LA: it cancels the virtual QED contributions associated to the collinear configuration in the final state, which is precisely what is taken into account by the SDK weak approach.

Drell-Yan
We start discussing the case of Drell-Yan production with leptons in the final state, namely the process pp → + − . 20 For all processes that will be considered in this section we have used the same input parameters already listed in Sec. 6.1. However, in the case of Drell-Yan, we have employed the complex-mass scheme [16,130,131] for the exact NLO EW corrections. 21 We use the PDF set NNPDF3.1 [132,133], in particular the   MZ and therefore passing the cuts. This configuration is not associated to any enhancement and therefore very rare, but in the on-shell scheme it leads to the evaluation of a resonant Z propagator with zero width and therefore it is inconsistent. cially a photon density following the LUXqed parameterisation [134,135]. The renormalisation (µ R ) and factorisation (µ F ) scales are both set equal to the partonic center-of-mass energy √ s. This set-up is common with all the other processes discussed in this section. In the Drell-Yan simulation the following cuts are imposed on the dressed leptons: 2) On the one hand, these cuts are imposed in order to resemble realistic experimental cuts for high-energy objects. On the other hand, they avoid additional logarithmic enhancements from collinear splittings appearing in the real radiation processes or even at the Born level. In Fig. 5 we show differential distributions for the transverse momentum of the electron, p T ( − ), for the transverse momentum of the leading (trailing) lepton, p T ( 1 ) (p T ( 2 )), and for the dilepton invariant mass m( + , − ).
The layout of each plot in Fig. 5, and in general of each plot in this section 22 , is the following. In the main panel we show the differential distribution at LO (solid blue line) and NLO EW (solid orange line) accuracy, where the exact O(α) corrections are taken into account. If the NLO EW prediction turns negative, meaning that NLO EW corrections are negative and larger than the LO in absolute value, the curve corresponds to its absolute value and is drawn as dashed. In the first inset we show the relative impact of EW corrections, δ X ≡ X/LO − 1, in different approximations. The solid orange line corresponds to the one in the main panel with the same style, i.e. the exact O(α) corrections (NLO EW), and the dotted orange line corresponds to the same case where the photon PDF has been set equal to zero (NLO EW, no γ). The other curves correspond to results in LA, with different assumptions. First, the solid curves include the SSC s→r kl contribution (SDK X , s → r kl ), while the dashed ones do not (SDK X ). Second, the green lines are obtained by simply omitting the QED and IR-sensitive terms, which are dubbed as "em" in the DP algorithm. This is analogous to the approach of e.g. Refs. [79,91] and dubbed here as SDK 0 . The red lines are instead obtained by completely removing the QED contribution, namely, following the procedure described in Sec. 4.1, the SDK weak approach. Both the SDK 0 and SDK weak predictions, similarly to the NLO EW ones in this section, include also the LO contribution. Needless to say, the closest a line is to the solid orange one, the better is the approximation of the exact NLO EW corrections. Therefore, in order to better judge this characteristic, in the second inset we zoom on the lines by simply plotting for each line in the first inset the difference with the solid orange one. Clearly, the reference prediction in LA is the solid red line, which both includes the SSC s→r kl contribution and is obtained via the SDK weak approach.
We remind the reader that neither the SDK 0 nor the SDK weak approach are equal to the approach dubbed as SDK in Sec. 6, which concerns the LA of the interference of Born amplitude and the IR-divergent virtual amplitude. The SDK prediction cannot be used for IR-safe observables. Moreover, the SDK 0 approach, even when the SSC s→r kl contributions are not taken into account, is not exactly equal to the one used so far in the literature, since we do include also in this case the terms proportional to iπΘ(r kl ). This is particularly relevant for Secs. 7.2 and 7.4, where 2 → 3 processes are considered.
Starting from the top-left plot of Fig. 5, the p T ( − ) distribution, we can see how sound is the LA in the high-energy limit. The distribution is spanning two orders of magnitude of p T ( − ), from 200 GeV to 20 TeV, and NLO EW corrections reach ∼ −50% in the tail. As can be seen in the second inset, the solid red line differs from the solid orange line by only a very few percents of the LO. The situation is opposite for the green solid line. In that case the difference with the solid orange line grows logarithmically with p T ( − ) and reaches ∼ −20% of the LO in the tail. This is a clear example of how the SDK weak approach can be superior to the SDK 0 one in the approximating the exact NLO EW corrections for IR-safe collider observables. For this observable the differences between dashed and solid lines of the same colour are negligible. This is not surprising, since large p T ( − ) implies large values of t and therefore |t|/s ∼ O(1), leading to small contributions from the SSC s→r kl terms. Finally, by looking at the difference between the solid and dotted orange line, we can also appreciate how the photon-initiated processes are relevant for this process and especially unavoidable also in the LA approximation in order to correctly approximate the exact NLO EW effects.
Moving to the top-right plot of Fig. 5, the m( + , − ) distribution, we see a very similar situation to the p T ( − ) distribution. In this case there is a visible difference between solid and dashed lines of the same colours, with the solid red line better approximating the exact NLO EW result than the dashed red line. Indeed, large m( + , − ) values do not imply large values for |t|. However, due to the cuts in (7.2) the effect of the SDK weak terms is very mild and also the dashed red line is leading to a very good approximation. The lower plots, the p T ( 1 ) (left) and p T ( 2 ) (right) distributions, display some differences w.r.t. the upper ones. For these two observable the LA is less efficient, indeed, as can be seen in the plots, the difference with the exact results can reach up to 5% of the LO also for solid red lines. However, this difference is converging to a constant value in the tail, for both distributions. Moreover the difference is of opposite sign for the two complementary observables, p T ( 1 ) and p T ( 2 ). This behaviour is due to the indirect cut that is affecting the recoiling particle when a particular value of p T ( 1 ) or p T ( 2 ) is considered, namely the trivial condition p T ( 1 ) > p T ( 2 ). This has nothing to do with EW Sudakov effects, but it is a particular feature of this kind observables: a similar pattern was observed for topquark pair production in Ref. [136]. Indeed, we also verified this explicitly for top-quark pair production, which we do not show here since results are very similar to the case of Drell-Yan that we are discussing. The bottom line is that the SDK weak approach including SSC s→r kl terms for LA can always have O(α) finite discrepancies with the exact result, especially it cannot take into account effects that are not related to a Born-like kinematic. On the other hand, also for these two observables there are not logarithmically enhanced differences with the exact NLO EW prediction, only constant terms are present. Instead, with the SDK 0 approach, the logarithmically enhanced differences are clearly visible.

ZZZ
In Fig. 6 we show plots, with the same layout of those in Fig. 5, for the process pp → ZZZ. This process has a neutral final state, so we do not expect large differences between the SDK 0 and SDK weak approaches. On the other hand, being a 2 → 3 process, the effect of the SSC s→r kl terms is supposed to be more relevant. The upper plots of Fig. 6 correspond to the transverse-momentum distributions of respectively the hardest Z-boson (p T (Z 1 )), the second-hardest Z-boson (p T (Z 2 )) and the softest one (p T (Z 3 )). The lower plots instead correspond to the invariant masses m(Z i , Z j ) of the three different Z-boson pairs. All the results have been obtained by applying the following cuts: 3) Similarly to (7.2), these cuts resemble realistic experimental cuts for high-energy objects, but they also avoid additional logarithmic enhancements from collinear splittings appearing in the real-radiation processes or even at the Born.
First of all, it is important to notice the size of the EW corrections. For most of the spectrum of all distributions, they are negative and larger than the LO in absolute value, reaching ∼ −200% of it in the tail. Since they are negative, this means that fixed-order NLO EW corrections are also negative in this regime and therefore non-physical. These distributions are a clear example of how large Sudakov logarithms, and in turn NLO EW corrections, can be at high energy. Also they clearly point to the necessity of resumming them in order to obtain sensible predictions. Here, on the other hand, we are not providing phenomenological predictions but rather showing the accuracy of the LA and testing its implementation in MadGraph5 aMC@NLO.
As expected, for all distributions, the difference between green and red lines (SDK 0 and SDK weak ) amounts to only few percents of the LO, with no clear logarithmic enhancement in the high-energy limit. Also as expected, the impact of the SSC s→r kl terms (solid versus dashed lines) is much larger for this process than for Drell-Yan production. In the upper plots of Fig. 5, the p T (Z i ) distributions, the dashed lines are differing from the solid ones by 5-10% of the LO for the full spectra, with the latter in turn differing only by a very few percents from the exact NLO EW prediction. The difference between dashed and solid lines is even larger in the lower plots, the m(Z i , Z j ) distributions, and especially a clear logarithmic trend can be observed. It is worth to stress that for all these distributions, with the exception of the far tail in the m(Z i , Z j ) ones, the inclusion of the SSC s→r kl terms leads to an accuracy of very few percents for corrections spanning from ∼-80% to ∼-200%. This is not the case for the pure LA without the SSC s→r kl terms.

WZ
We now move to the case of a couple of processes where both the inclusion of the SSC s→r kl terms and the use of SDK weak is relevant. We start by showing differential distributions for the process pp → W + Z, where results have been obtained by using the following cuts |η(V i )| < 2.5 , m(W + , Z) > 1 TeV , ∆R(W + , Z) > 0.5 . (7.4) Again, these cuts resemble realistic experimental cuts for high-energy objects, but they also avoid (part of the) additional logarithmic enhancements from collinear splittings appearing in the real-radiation processes or even at the Born.
In Fig. 7 we show the transverse momentum of the hardest (p T (V 1 )) and softest (p T (V 2 )) recombined vector-bosons and their invariant mass (m(W + , Z)). Similarly to the case of leptons (7.1), the recombination is performed by recombining any charged vector boson V i with photons that satisfy the condition ∆R(V i , γ) < 0.4. We also show the transverse momentum (p T (j V,2 )) of the softest jet that is tagged as a V -jet, namely containing one of the two vector bosons. They layout of the plots is very similar to the ones of Figs. 5 and 6, but with an important difference: in the second inset we plot for each line in the first inset the difference with the dotted orange one, i.e., the exact NLO EW prediction where the photon PDF has been set equal to zero. The reason is the following. The W Z  TeV. Comparison between exact NLO EW predictions excluding photon-initiated processes (dotted orange) and their LA approximations in the SDK weak (red) and SDK 0 (green) approaches, including (solid) or excluding (dashed) the SSC s→r kl terms. More details are given in the text.
production is affected by giant K-factors not only at NLO QCD [137,138], but also at NLO EW [90,139] precisely due to the opening of the qγ → W + Zq channel. For large value of p T (V 1 ) this channel induces large positive corrections of order L(p 2 T (V 1 ), M 2 V 2 ), due to configurations in which V 1 recoils against the emitted quark q and V 2 is soft and collinear to it. This effect has nothing to do with virtual NLO EW corrections and cannot be expected to be captured by the DP algorithm, both in its SDK 0 and SDK weak adaptations for physical observables. Therefore, in order to test the LA, we use as reference the NLO EW prediction where the photon PDF has been set equal to zero. We also plot in the main panel the NLO EW prediction without photon PDF as dotted orange line and its absolute value as dashed-dotted when it becomes negative. We remind the reader that the photon PDF is entering the process starting with NLO EW real emission of quarks, while it is not present at LO.
Before commenting each plot of Fig. 7, it is again important to note how LA in general approximates very well the exact NLO EW predictions excluding photon-initiated processes, reaching almost -200% of the LO in the tail of the distributions. This can be appreciated in the first inset of each plot. There, as in the main panel, the size of the qγ → W + Zq channel can also be appreciated. Its impact is maximal in the p T (V 1 ) distribution, due to the aforementioned L(p 2 T (V 1 ), M 2 V 2 ) enhancement, and minimal in the p T (V 2 ) distribution since large values of p T (V 2 ) forbid the kinematic configurations that precisely lead to the enhancement for the p T (V 1 ) distribution.
In the top-left plot (p T (V 1 )) we can see that none of the curves is really flat. This is not surprising, since the LA works very well when the cross section is dominated by the Born-like kinematic, which is not the case for large p T (V 1 ), where V 1 can in principle recoil against a very hard photon. Still, the solid red line (SDK weak with SSC s→r kl terms included) leads to the best approximation and differs from the reference value, the dotted orange line, by less than 10% of the LO also in regions where the NLO EW corrections are ∼ −200%. The superiority of the SDK weak approach (red) over the SDK 0 one (green) and of the inclusion of the SSC s→r kl terms (solid) over their exclusion (dashed) can be better appreciated in the top-right plot (p T (V 2 )) and especially in the bottom-right one m(W + , Z). The difference between the dotted orange and solid red curves is only very few percents of the LO over the full range.
In the bottom-left plot, the p T (j V,2 ) distributions, all the curves are equal to those of the top-right one p T (V 2 ) but one: the solid orange line. In this context the only difference between a recombined V -boson and a V -tagged jet is that in the latter the V boson can be recombined also with quarks. When the qγ → W + Zq channel is included, soft+collinear configurations for (V 2 ) in the final-state splitting q ( ) → q V 2 are avoided at large p T (V 2 ), but the purely collinear ones still survive. The recombination of V 2 with quarks has therefore a sizeable impact. The bottom line is that the SDK weak approach takes into account only those emissions that are unavoidable from an IR-safety point of view, but additional realemission contributions can be present and sizeable (for instance the photon-induced quark radiation), they can lead to enhancements factorising different Born matrix elements (for instance the case discussed here where σ(qγ ) and their impact can depend on how the different particles are clustered among each other (for instance V -bosons versus V -tagged jet). These contributions cannot be taken into account via the DP algorithm. The LA and in particular the DP algorithm as implemented in MadGraph5 aMC@NLO allow for a fast and, as shown in the previous examples, very reliable approximation of   fixed-order exact NLO EW corrections. On the other hand, it cannot substitute the exact calculation. LA can be used as a starting point for improving the fixed-order NLO EW, by e.g. resumming the large Sudakov logarithms or alternatively for performing fast simulations with MadGraph5 aMC@NLO including the EW dominant effects at high-energy. The latter option, however, should be always cross-checked with an exact calculation before being used for phenomenological predictions.

WWW
As last example, we show in Fig. 8 distributions for the process pp → W + W + W − , where for the calculations the following cuts have been employed following the same rationale behind the cuts that have been used for the other processes considered in this section. In Fig. 8 we consider the same observables already considered in Fig. 6, with Z bosons replaced by W bosons. The layout of the plots is the same of those in Fig. 7 for W + Z production, where the reference line is the dotted orange one, the NLO EW prediction without photon-initiated processes. Indeed, also in this case the opening of the photon-induced channels in the NLO EW corrections (qγ → W + W + W − q ) induces additional effects that cannot be captured by the DP algorithm.
Similarly to the case of W + Z production in Fig. 7, we see that both the choices between SDK 0 and SDK weak approaches and between the inclusion or exclusion of SSC s→r kl terms are relevant. Also, the exclusion of the photon-initiated processes is unavoidable for the comparison between LA and exact NLO EW results. As any other observable considered in this chapter, once again we see very large effects from NLO EW corrections (reaching almost −200% of the LO). Moreover, the solid red line, the LA in the SDK weak approach including SSC s→r kl terms, is able to approximate within the 10% level the exact NLO EW effects excluding photon-initiated processes (dotted orange line) over the full range of the observables considered the plots of Fig. 8. The difference between the two aforementioned accuracies, the solid red line in the second inset, is in general a constant over the range of each observable considered. The only exception is m(W 2 , W 3 ) where a logarithmic trend is visible also for the red-solid line. A similar effect, although milder, could be observed also for the m(Z 2 , Z 3 ) distribution in ZZZ production in Fig. 6. To the best of our knowledge, given the cuts in (7.5), this effect is not due to any possible enhancement arising from real emission contributions. However, as explained in Sec. 6, possible additional angulardependent effects involving ratios of invariants can be present and not captured in the DP algorithm, regardless of the inclusion of the SSC s→r kl term.

Conclusions and outlook
In this work we have revisited the algorithm of Denner and Pozzorini [39] for the calculation of one-loop EW Sudakov logarithms, denoted in the previous sections as DP algorithm. We have introduced several novelties that concern different aspects of the electroweak Sudakov approximation. Moreover, we have implemented the DP algorithm, together with the novelties introduced, in the MadGraph5 aMC@NLO framework. Thanks to this new implementation we have provided in a completely automated approach several numerical results for different production processes, both for the virtual contribution in specific phase-space points and for physical (IR-safe) differential cross sections at a 100 TeV proton-proton collider. All the numerical results that we have obtained corroborate the relevance of the novelties introduced in this work and also demonstrate the correctness of the implementation of the algorithm in the code.
In particular we have introduced the following novelties: • We have reframed the DP algorithm by setting the mass of the photon and lightfermion masses exactly to zero, regularising IR divergences by mean of Dimensional Regularisation (DR), as in modern NLO EW calculations and Monte Carlo implementations. Reframing the algorithm with the language of modern calculations, it allows for further developments and compatibility with the state-of-the-art tool implementations.
• We have modified part of the expressions in order to take into account additional angular dependences, without assuming that all the invariants are of the same size of s. These modifications correspond to the terms that have been dubbed in the work as SSC s→r kl ; they have been obtained keeping track of the dependence on any invariant r kl in the derivation of the subleading soft-collinear logarithms. As already mentioned in various points in the text, even with this improvement, the full control of logarithms involving the ratios of the different |r kl | and s cannot be achieved via the DP algorithm. Information on the internal structure of the diagrams is unavoidable for this purpose. On the other hand, the SSC s→r kl terms are sensitive to the presence of a hierarchy among the various invariants that characterise a specific phase-space point and substantially improve the approximation of the aforementioned class of logarithms. In order to support this statement, we have showcased for several processes and observables the superiority of the Sudakov approximation including the SSC s→r kl terms.
• We have identified an imaginary term that was omitted in Ref. [39], which cannot be in general neglected for 2 → n processes with n > 2. This term is also present in the new SSC s→r kl terms that have been introduced. We have shown the relevance of this imaginary term for correctly capturing effects of order α log(s/M 2 W ).
• In this work we did not focus only on one-loop EW corrections to amplitudes, as done in the pioneering work Ref. [39], but we have considered also the virtual NLO EW corrections to the LO cross sections. The two cases are trivially related only when the entire LO factorises a single combination of α S and α powers. However, more in general, NLO EW corrections can also include effects due to loops of QCD origin. We provide the additional terms, denoted in the paper as δ QCD LA , that are necessary in order to take into account in Sudakov approximation also contributions from the aforementioned QCD loops. We have proved the relevance of δ QCD LA for correctly capturing effects of order α log n (s/M 2 W ) with n = 1, 2, when also QCD loops contribute to virtual NLO EW corrections to the LO cross sections.
• We describe how to modify the DP algorithm in order to exclude the contribution of photons and also the contribution of gluons from the QCD terms mentioned in the previous bullet. We dub this approach as SDK weak . In the context of differential cross sections and IR physical observables, we show how this approach is superior to the standard approach used in the literature, dubbed in this work as SDK 0 , which corresponds to simply removing the IR-divergent logarithms involving the scale M W and the IR cut-off.
We discuss also the technical steps of the implementation of the DP algorithm, together with the novelties introduced in this work, in the MadGraph5 aMC@NLO framework. The choice of the MadGraph5 aMC@NLO framework, which already automates the calculation of the exact NLO EW corrections, has been crucial for the validation of our work and for demonstrating the correctness of our implementation of the Sudakov approximation, also denoted in the paper as leading approximation (LA). Systematic comparisons between exact O(α) corrections, the NLO EW, and their LA with different assumptions (SSC s→r kl , imaginary term, δ QCD LA , SDK weak , etc.) have been performed in order to check that all the logarithmic terms are correctly captured.
Besides being essential for validation, the implementation in a framework as Mad-Graph5 aMC@NLO opens up for several new possibilities. First of all, the necessary ingredients for performing completely automated NLO EW corrections matched with resummed EW Sudakov logarithms are now available in a single tool. Then, a fast and stable method for approximating dominant effects from EW corrections is now available. Especially, the quality of the approximation can be checked via the NLO EW exact calculation, which however we remind the reader cannot in general be substituted by its LA. Sudakov effects can be also more easily integrated in Monte Carlo simulations and generations of events. These are only a few of the possible outcomes of this work, which sets the basis for future works both at the phenomenological and formal level. Another important possible development of this work is, e.g., the possibility of compute approximate EW corrections also in extensions of the SM. Indeed, only tree-level Feynman rules are needed. However, for a given BSM model, the associated Universal Feynrules Output (UFO) model [140], which is the format used in MadGraph5 aMC@NLO for encoding the Feynman rules, should provide some extra informations that are not present at the moment: how particles are arranged in SU(2) multiplets and the values of their electroweak couplings (charge, isospin component and/or hypercharge). At the moment, this information had to be written by hand in the code. We thus envisage an extension of the UFO format in this direction.
As a final remark, we want to stress once again that the Sudakov approximation cannot substitute exact NLO EW calculations. As we said, this approximation can be exploited for performing fast simulations including the EW dominant effects at high energy. However, it should be always cross-checked with an exact calculation before being used for phenomenological predictions and especially comparisons with data. Moreover, as explained many times in the text, the approximation works when all the invariants r kl satisfy the relation |r kl | M 2 W . This implies, for instance, that the Sudakov approximation can never be directly exploited for a process involving resonances. In the SM this typically means the production and decay of the heavy particles W, Z, H or t. On the other hand, given such a type of process, the Sudakov approximation can be applied for the calculation of the same process including only the production but not the decays of the aforementioned heavy particles, and the decay can be subsequently taken into account.