The automation of next-to-leading order electroweak calculations

We present the key features relevant to the automated computation of all the leading- and next-to-leading order contributions to short-distance cross sections in a mixed-coupling expansion, with special emphasis on the first subleading NLO term in the QCD+EW scenario, commonly referred to as NLO EW corrections. We discuss, in particular, the FKS subtraction in the context of a mixed-coupling expansion; the extension of the FKS subtraction to processes that include final-state tagged particles, defined by means of fragmentation functions; and some properties of the complex mass scheme. We combine the present paper with the release of a new version of MadGraph5_aMC@NLO, capable of dealing with mixed-coupling expansions. We use the code to obtain illustrative inclusive and differential results for the 13-TeV LHC.


Introduction
With more than two-thirds of the LHC Run II completed, and a collected integrated luminosity in excess of 85 fb −1 at 13 TeV, it appears that physics beyond the Standard Model is not eager to be discovered. This could be either because its characteristic scale is larger than usually assumed, and thus beyond the LHC reach; or because its signals are particularly difficult to identify by means of final-state direct searches. While we have no way at present to tell which of these two scenarios is the one realised by Nature, we know that they entail different long-term strategies. In the former instance, we must be able to explore larger scales, chiefly by increasing collider c.m. energies; in the latter instance, we must collect more statistics and fight against systematics. In both cases, indirect searches might turn out to play a vital role.
In order to cope with this situation, phenomenologists must continue in the trend which has now been established for a few years, namely to progress in the direction of increasing both the flexibility and the precision of their calculations. In terms of flexibility, the key aspects are the ability of computer codes to readily deal with the matrix elements relevant to both new-physics models and the SM, and that of embedding these matrix element results in fully realistic final-state simulations, such as those provided by parton shower Monte Carlos. As far as precision is concerned, in the vast majority of applications this is a synonym for the computation of higher orders in a fixed-order coupling-constant perturbative expansion. Which specific perturbative orders depends obviously on the theory one considers. Owing to both the large numerical values assumed by α S and the paramount importance of hadroncollision physics in the LHC era, QCD has played a particularly prominent role in recent years -fully-differential NLO results are by now standard, including the cases of processes with very involved final states, while more and more NNLO and even N 3 LO predictions are becoming available, for low-multiplicity reactions and with different degrees of inclusiveness (see e.g. ref. [1] for a recent review).
Currently, there is therefore a compelling case for considering in full generality the calculation of NLO contributions in the EW theory, and this for at least three reasons. Firstly, based on the values of α S and α, one expects NNLO QCD and NLO EW effects to be numerically comparable. Secondly, this naive scaling behaviour could actually be violated in those regions of the phase space associated with large mass scales, since there the coefficients of the EW series might grow faster than their QCD counterparts, owing to the presence of large Sudakov logarithms (see e.g. refs. [2][3][4][5][6][7] for discussions about their origin and universal nature) -large transverse momenta are a particularly important example, given the relevance of the high-p T regions in new-physics searches. Thirdly, it is likely that among future colliders there will be at least one e + e − machine, and possibly an eH one, for whose physics simulations EW computations will be crucial.
When combining the necessity of the calculation of EW NLO corrections with the requirement that they be flexible and available for arbitrary processes, one is naturally led to their automation, not least because the case for automation is strongly supported by the striking success that this strategy has enjoyed in the case of QCD, where NLO results are now mass-produced and constitute the backbone of ATLAS and CMS pp simulations. Indeed, although not as comprehensive as for QCD computations, there has been steady progress in the automation of EW NLO corrections, for both one-loop and realemission contributions, by collaborations such as RECOLA [8,9] with Sherpa [10,11], OpenLoops [12] with Sherpa, GoSam [13,14] with either MadDipole [15,16] or Sherpa, and MadGraph5 aMC@NLO [17]. Recent results obtained with these tools  clearly demonstrate how automation is allowing us to attack problems whose complexity is too great to justify their solutions through traditional approaches.
The goal of this paper is that of presenting the most important features that underpin NLO calculations performed by MadGraph5 aMC@NLO (MG5 aMC henceforth) in a mixed-coupling scenario, i.e. when two coupling constants are simultaneously treated as small parameters in a perturbative expansion. While the QCD+EW case, in the context of which MG5 aMC is able to compute the NLO QCD and the (so-called) NLO EW corrections among other things, is the most prominent in today phenomenology, the code is not necessarily restricted to that, and is structured to handle analogous situations in user-defined theories (including, in particular, EFTs). We also combine the present paper with a new release of the MG5 aMC code, which will be the first public one compatible with a mixed-coupling expansion -previous applications that included non-QCD effects [22,25,32,33,39,40,[42][43][44] had been obtained with preliminary and still-private versions. Conversely, we stress that EW-loop-induced processes could already be automatically generated by MG5 aMC, owing to the work of ref. [45]. This paper is organised as follows. In sect. 2 we lay out the scope of our work, discuss the strategy upon which MG5 aMC is based, and point out the limitations of the current version of the code. In sect. 3 we extend the short-distance FKS [46,47] formulae to the mixed-coupling case. Sect. 4 shows how one can apply the FKS subtraction formalism to those cases in which the use of fragmentation functions is required. In sect. 5 we discuss some features of the complex-mass scheme [48,49] that are particularly relevant to its automation in MG5 aMC. In sect. 6 we present some illustrative results, both inclusive and at the differential level, that have been obtained in a 13-TeV LHC configuration. Finally, in sect. 7 we draw our conclusions. Some extra technical material is reported in the appendices.

General considerations
We start by writing an observable Σ(α S , α) to all orders in α S and α by adopting the notation introduced in ref. [17]: where we have identified the LO (Σ (LO) ) and NLO (Σ (NLO) ) contributions in eq. (2.2) with the p = 0 and p = 1 terms in eq. (2.1). The integer numbers k 0 , c s (k 0 ), c(k 0 ), and ∆(k 0 ) are process-specific quantities (of which we shall give an explicit example below) that obey the constraint k 0 = c s (k 0 ) + c(k 0 ) + ∆(k 0 ). We can further decompose these contributions by defining terms that factorise a single coupling-constant combination α n S α m ; each of these terms corresponds to a single value of q in the second sum in eq. (2.1). Explicitly: In a well-behaved perturbative series, and given that α α S , eqs. (2.3) and (2.4) imply: This hierarchy suggests (see ref. [22]) to call Σ N p LO i the leading (i = 1) or i th -leading (i > 1: second-leading, third-leading, and so forth) term of the N p LO contribution to the cross section. It is customary to identify Σ NLO 1 and Σ NLO 2 with the NLO QCD and the NLO EW corrections, respectively. While this is unambiguous in the case of Σ NLO 1 , it is somehow misleading in the case of Σ NLO 2 , for two reasons. Firstly, at one loop there is no clear-cut way to define pure-EW contributions on a diagrammatic basis. Secondly, Σ NLO 2 may receive contributions from the so-called heavy-boson radiation (HBR), namely from diagrams that feature the emission of a real W , Z, or H boson (see e.g. eq. (2.3) of ref. [25]) -these are typically not included in what are conventionally denoted as NLO EW corrections. With these caveats in mind, in what follows we shall also often refer to Σ NLO 2 as the NLO EW contribution, in those cases where no ambiguity is possible. In order to give an explicit example, let us consider the case of tt production in association with a heavy EW boson. Equations (2.1)-(2.4) read: (2.6) Σ (NLO) ttV (α S , α) = α 3 S α Σ 4,0 + α 2 S α 2 Σ 4,1 + α S α 3 Σ 4,2 + α 4 Σ 4,3 , ≡ Σ NLO 1 + Σ NLO 2 + Σ NLO 3 + Σ NLO 4 . (2.7) These imply that, for this process, we have k 0 = 3, c s (k 0 ) = 0, c(k 0 ) = 1, and ∆(k 0 ) = 2.
Finally, we point out that eqs. (2.1)-(2.4), although written here for the QCD+EW case, actually apply to the perturbative expansion in any two couplings 1 α 1 and α 2 . The formal replacements α S → α 1 and α → α 2 suffice to obtain the general formulae.
As it has been anticipated in sect. 1, our final goal is the automation of the computation of eqs. (2.3) and (2.4) in MG5 aMC. We briefly recall here the basic building blocks that form the core of the code. MG5 aMC makes use of the FKS method [46,47] (automated in the module MadFKS [50,51]) for dealing with IR singularities. The computations of one-loop amplitudes are carried out by switching dynamically between two integral-reduction techniques, OPP [52] or Laurent-series expansion [53], and TIR [54][55][56]. These have been automated in the module MadLoop [17,57], which in turn exploits The implications of what has been said above are the following. The building blocks of the MG5 aMC code have been made compatible with the structure of the generic mixedcoupling expansion, eqs. (2.3) and (2.4). This implies, in particular, upgrading the handling of the inputs, and implementing a much more involved bookkeeping, both of which are due to the necessity of retaining independent control on the different Σ k,q contributions. In terms of matrix-element computations, the code is aware of the possible presence of more than one type of interactions (e.g. QCD and EW), whose details are assumed to be inherited from the model. In particular, inspection of the model allows MG5 aMC to construct all of the real-emission and one-loop diagrams that contribute to the process under study at the desired perturbative order. Crucially, it also allows the code to figure out the structure of the IR singularities, by considering the set of all possible 1 → 2 branchings with at least one outgoing massless particle, and so to set up the appropriate FKS subtractions. More details on the calculation strategies underlying MadFKS and MadLoop can be found in sect. 2.4.1 of ref. [17], and sects. 2.4.2 and 4.3 of ref. [17], respectively. Mixed-coupling expansion capabilities have been gradually added to MadFKS and MadLoop, in different stages for the two modules, during the course of the work relevant to refs. [22,25,32]. They have been completed and validated for the present paper. The bottom line is that the typical MG5 aMC usage in the case e.g. of a QCD+EW expansion may read as follows 2 : MG5 aMC> import model myNLOmodel w qcd qed MG5 aMC> generate p 1 p 2 > p 3 p 4 p 5 p 6 QCD=n max QED=m max [QCD QED] with p i (multi)particles that belong to the spectrum of the model myNLOmodel w qcd qed. The syntax above implies computing the following LO and NLO contributions: NLO : α n S α m , n ≤ n max + 1 , m ≤ m max + 1 , n + m = k 0 + 1 . (2.9) We point out that the largest power of α S in eq. (2.9) is exactly one unity larger than its LO counterpart in eq. (2.8) because of the presence of the keyword [QCD] in the processgeneration command. This instructs the code to consider all diagrams that feature two extra QCD vertices 3 w.r.t. those present at the Born level. Likewise, it is the keyword [QED] that sets the largest power of α that appears at the NLO level equal to m max + 1.
Either keyword can be omitted -we shall give explicit examples in sect. 6. More details on the syntax above, relevant to the inner workings of MG5 aMC, are given in appendix D. We now turn to underscore an issue which is specific to QED, and that stems from the fact that, in such a theory, photons and leptons can be regarded both as particles that enter the short-distance process and as observable (taggable) objects. This point has been already discussed in ref. [32], and we shall limit ourselves here to summarise the conclusions of that paper. The key point is that short-distance photons and massless leptons can be identified with the corresponding taggable objects only up to a certain Σ NLO i 0 term, beyond which (i.e. for Σ NLO i , i > i 0 ) this identification leads to IR-unsafe observables. The value of i 0 is process dependent; typically, IR unsafety manifests itself in the third-or even the second-leading NLO contribution 4 . The solution proposed in ref. [32] is as follows: 1. Short-distance photons and massless leptons are not taggable objects.

2.
A taggable photon is a photon that emerges from a fragmentation process. 3. A taggable massless lepton is either a lepton that emerges from a fragmentation process or a dressed lepton, i.e. an object whose four-momentum is equal to that of a very narrow jet that contains the short-distance lepton. 4. These rules imply that photons and massless leptons must be treated on the same footing as gluons and quarks in short-distance computations (democratic approach). 5. Short-distance computations should be performed in MS-like EW renormalisation schemes, such as the G µ or α(m Z ) ones, regardless of the initial-and final-state particle contents of the process of interest.
We point out that, since the O(α 0 ) term of both the photon and the lepton fragmentation functions is equal to δ(1 − z), the above prescriptions exactly coincide with the usual procedure followed in the calculation of NLO QCD corrections, where no distinction is made between short-distance and taggable EW objects. Furthermore, we would like to stress that, as was already mentioned in ref. [32], a direct consequence of the evolution equations for fragmentation functions at the first non-trivial order in QED is that the z → 1 dominant term in the photon-to-photon fragmentation is equal to α 0 /α(Q 2 ) δ(1 − z), with α 0 an O(α) small-scale constant, and Q 2 a scale of the order of the hardness of the process. This shows that, by retaining only such a term, a fragmentation-function approach allows one to recover naturally and in a straightforward manner the results obtained in the standard EW approach that makes use of the α(0)-scheme. On top of being fully general and backward-compatible with established procedures, working with fragmentation functions has also the appealing feature of putting QCD and QED on a similar footing, and of rendering conceptually alike the treatment of initial-and final-state photons 5 . The case of massless leptons has not been discussed in ref. [32], and we also leave it to future work; however, we point out that universal logarithmic mass effects can be included in lepton fragmentation functions through the computation of its perturbative part, for example by working in analogy to what has been done in QCD for the case of massive quarks [72]. In conclusion, there is a compelling motivation for a fully general treatment of fragmentation functions in the context of NLO computations, which is what will be done in sect. 4.
We also note that items 4 and 5 in the list above put the MC@NLO-type matching to QED showers on the same footing as its QCD counterpart, at least for those QED showers implemented in a QCD-like manner.
We conclude this section by enumerating the limitations of the new version of the MG5 aMC code whose release is associated with the present paper. Firstly, the implementation of eqs. (2.3) and (2.4) implies that the automated calculation of NLO QCD+EW effects is always carried out in an additive scheme. If one is interested in a multiplicative approach (where QCD and EW terms are assumed to factorise), one must compute the leading and second-leading contributions separately, and manually combine them afterwards 6 . Secondly, the implementation of the convolution of short-distance cross sections with fragmentation functions is not made publicly available (thus, observables with tagged photons and/or leptons cannot be constructed). This stems from phenomenological considerations, given that the perturbative terms for which those functions are a necessity are strongly suppressed, and the functions are not easy to extract from current datasets 7 . Thirdly, although the implementation of the MC counterterms for QCD-like QED showers 5 Obviously, differences remain. In particular, the different scales associated with the evolution of the photon-to-photon fragmentation function and the photon-in-photon PDF imply that the issues (see e.g. ref. [71]) relevant to photon-initiated processes are simply not present here. 6 We point out that all user-selected cross section contributions are separately available during the course of a single run. Therefore, their combination can be straightforwardly achieved at the analysis level. 7 Having said that, the possibility of using purely-theoretical, leading-behaviour photon and lepton fragmentation functions will be made available in the MG5 aMC code in the near future. This will allow one to recover standard results for processes where such particles are tagged.
(which are available in both Pythia8 [73] and Herwig7 [74]) has been completed, some further validation will be necessary in order to bring the MC@NLO matching in the QED sector on par with the QCD one. Thus, in such a sector the present release of MG5 aMC is restricted to performing fixed-order computations. Finally, ISR and beamstrahlung effects are not yet implemented in MG5 aMC, and therefore we refrain from presenting e + e − -collisions results in this paper.

FKS subtraction in mixed-coupling expansions
In this section we introduce the basic mechanisms and the notation that underpin the implementation of a mixed-coupling scenario in MadFKS; in the process, we shall thus extend the FKS subtraction to cover such a case. In keeping with what has been done so far and in order to be definite, we shall deal explicitly with the QCD+QED case, which has the further advantage of being essentially a worst-case scenario as far as IR subtractions are concerned 8 . For backward compatibility, the notation is modified in a manner as minimal as is possible w.r.t. the one introduced in ref. [50]; eq. (n.m) of that paper will be denoted by eq. (I.n.m) here.

Matrix elements
What is done in sect. 3.1 of ref. [50], and in particular eq. (I.3.1)-eq. (I.3.18) is unchanged, bar for the following amendments: • n (B/R) L is the total number of light quarks, gluons, massless charged leptons, and photons, at the Born (B) and real-emission (R) level; • n H is the number of massive particles that are either strongly interacting, or electrically charged, or both; • n ∅ is the number of particles that do not belong to either of the two previous categories; • n I = 1 regardless of the type of incoming particles (see eq. (I.3.11)-eq. (I.3.13)); • the symbol ⊕ as is used in eq. (I.3.18) understands either a QCD or a QED vertex.
The colour operators are defined in eq. (I. 3.26); that definition must be supplemented with: In a mixed-coupling scenario there will be a QED analogue of eq. (3.6), i.e. the charge-linked Born's that stem from the insertion of a soft-photon line. We have: M (n,0) QED(p,q)kl (r) = − 1 2s 2 − δ kl ω(I 1 )ω(I 2 ) colour spin p 1 p 2 q 1 q 2 δ p 1 +p 2 ,p δ q 1 +q 2 ,q (3.8) where e(I) is the electric charge of particle I in unit of the positron charge, and: s(I) = 2 outgoing (anti)particle , 1 incoming (anti)particle . This result is identical to that which one obtains directly from the Altarelli-Parisi QED kernels, as it must as a consequence of the commutation of the soft and collinear limits (see e.g. ref. [75]). Such kernels are reported in appendix A, together with the results for all of the QED-charge factors that are relevant to FKS subtraction, and that will appear in sect. 3.2. By construction, a given colour-linked (charge-linked) Born will be non-null only if both k and l correspond to strongly-interacting (electrically charged) lines. This implies that, for any given l and T, several terms in the sum on the l.h.s. of eq. (3.14) will be identically equal to zero. Finally, the one-loop matrix elements that generalise those of eq. (I.3.25) are:

Short-distance cross sections
The short-distance cross sections emerge from applying the FKS subtraction procedure to all possible IR-divergent configurations, be them QCD-or QED-induced. The key observation is that such configurations are "squared" in nature (i.e. they never lead to relative O(g S e) corrections). In other words, one either inserts a soft-gluon line or a softphoton line; likewise, a collinear splitting is of either QCD or QED type. Hence, provided that the S functions take into account all IR-divergent sectors regardless of their origins, the subtraction proceeds as before. In turn, this is equivalent (see sect. 5.2 of ref. [50]) to constructing the P FKS set so that all singularities are represented there. This can be done by following the procedure advocated in sect. 2.4.1 of ref. [17], which has been implemented in that paper for the QCD case, and extended in the course of the present work to cover mixedcoupling scenarios. Basically, such a procedure simply requires one to consider, within the adopted theory model, all possible 1 → 2 branchings with at least one outgoing massless particle. In so doing, one treats all massless particles in a democratic manner (which is consistent with item 4 in sect. 2) which, among other things, implies that any two-particle combination in the initial state is allowed. This is the main reason that motivates the amendments listed at the beginning of sect. 3.1.
The hadron-level result is obtained by extending eq. (I.C.1) (for simplicity, we start by setting the renormalisation and factorisation scales equal to each other): where the sums are constrained by eqs. (3.2) and (3.3), which are always understood. The sums over incoming-parton species are also left implicit. The symbol in eq. (3.18) denotes the usual initial-state convolution of short-distance cross sections and PDFs, i.e. for a given initial-state leg associated with hadron momentum P : With the caveat on the S functions mentioned above, the (n+1)-body contribution has the same form as in eq. (I.4.29), with the matrix elements taken from eq. (3.4). The degenerate (n + 1)-body contributions can be derived from eq. (I.4.41) and eq. (I.4.42). We adopt a short-hand notation to re-write the sum of those contributions for a given initial-state parton configuration (a, b) as follows: where The result for the mixed-coupling case then reads as follows: where we have admitted the possibility of different PDF schemes relevant to QCD and QED evolution. It is immediate to see that the first (last) two terms on the r.h.s. of eq. (3.22) correspond to the blue right-to-left (red left-to-right) arrows of fig. 1 of ref. [17]. Thus, they will identically vanish when (p, q) are associated with the pure QED (QCD) NLO term. The n-body cross section in eq. (3.18) can be read from eq. (I.4.3): We point out that, for a given (p, q), either the Born or the other three terms on the r.h.s. vanish, owing to eqs. (3.2) and (3.3); the final result is correct because p and q are summed over in eq. (3.18). The Born term differs from that of eq. (I.4.4) only by notation: The term dσ (C,n) (p,q) collects the Born-like remainders of the final-and initial-state collinear subtractions. Thus, in full analogy with eq. (3.22) and by taking eq. (I.4.5) and eq. (I.4.6) into account, we have: with (T = QCD, QED): The colour and charge factors that appear in eq. (3.27) are reported in appendix A. Analogously, from eq. (I.4.12) one obtains the soft term: Note that the integrated eikonal has the same expression in the QCD-and QED-induced terms, being of kinematical origin. Finally, the virtual contribution is analogous to that of eq. (I.4.14): The finite part of the virtual contribution that appears in this expression depends on what is included in the divergent part of the one-loop matrix elements. We use the same conventions as in ref. [50], and write the latter as in eq. (I.B.1) and eq. (I.B.2): where (T = QCD, QED): Finally, the separate dependence on the renormalisation and factorisation scales can be reinstated by following the procedure outlined in appendix C of ref. [50]. We report here in appendix B its extension to the mixed-coupling case. We remind the reader that MG5 aMC allows the evaluation of the hard-scale and PDF uncertainties at no extra computational costs, by following the procedure introduced in ref. [76]. This feature is also supported in mixed-coupling applications. Note, however, that the scale dependence of α is ignored, and that, in keeping with what is done in current PDF fits, the QCD-and QED-factorisation scales are set equal to each other.

FKS subtraction with fragmentation
As was discussed in sect. 2, the computation of sufficiently subleading NLO contributions might require that some physical objects be defined through fragmentation functions (FFs henceforth). This motivates the extension of the FKS subtraction procedure to cover the cases where FFs would also be present -in other words, measurable quantities may be obtained by means of FFs, or with a jet-finding algorithm, or as stable-particle taggable massive objects (e.g. a Higgs boson); any combination of the objects thus constructed is allowed. This must be done by respecting the chief FKS subtraction characteristics, namely that it be universal, and observable-and process-independent.
In order to proceed, we observe that in NLO computations at most two final-state partons are not well separated. Since hard and isolated partons can be fragmented in a trivial manner, this implies that the only non-trivial case we shall have to deal with is that which features a single FF. Furthermore, sect. 3 shows that the FKS subtraction in a mixed-coupling scenario can readily be obtained from its QCD counterpart, and from the knowledge of the relevant kernels and colour or charge factors. Therefore, in order to simplify the notation of the formal proof that follows, we limit ourselves to presenting explicitly only the QCD case; the mixed-coupling results can then be readily obtained by following the same procedure as in sect. 3. For the same reason, we assume that all final-state particles are massless.
Following what has been done in the original paper [46], such a proof is essentially that of the cancellation of the IR divergences that arise in the intermediate steps of the calculation of a generic final-state observable. Its by-products are the main results we are interested in, namely the IR-finite short-distance cross sections that can be numerically integrated.

Fragmentation in perturbative QCD
The case we are considering, that of a single FF, is by definition equivalent to the computation of a single-hadron cross section in QCD. That is typically written as follows (see e.g. refs. [77,78]): Here, dσ is the short-distance partonic cross section (for a given partonic process), already convoluted with the PDFs and integrated over the degrees of freedom not explicitly indicated in eq. (4.1); D (ap) H (ζ) is the FF of parton p (with momentum k p and flavour a p ) which fragments into hadron H (with momentum K H ); the latter carries a fraction ζ of the three-momentum of the former (in leading-twist QCD factorisation, fragmentation is strictly collinear); the index p runs over all final-state partons. Equation (4.1) is not suited to the definition of a parton-level generator through the FKS (or any equivalent) procedure; to that end, a fully differential form is required. Such a form is: where in the r.h.s. one understands that the momenta of hadron H and parton p are related by K H = ζ k p .
Equation (4.2) has an unusual form, because the variables are seemingly separated (i.e. at variance with eq. (4.1) there is no dependence on ζ in the partonic cross section dσ on the r.h.s.). However, this is purely an artifact that stems from the fully differential nature of eq. (4.2). One must keep in mind that formulae of this kind (including those of FKS) are never meaningful if they do not understand the inclusion of an IR-safe observable. In other words, eq. (4.2) formally expresses the computation of the expectation value: In order to see that eqs. (4.1) and (4.2) are mutually consistent, we demonstrate that the former equation can be derived from the latter one. By writing the partonic cross section in terms of a matrix element and the phase space: with N a normalisation factor, one obtains: and keeping in mind that parton p is massless (k 0 p = | k p |) we have: whence: since by construction for anyk p . It is apparent that, after summing over p, eq. (4.8) is identical to eq. (4.1), if one assumes that hadron H is massless, i.e. K 0 H = | K H |. This is of course not the case for any physical hadron, but it appears to be a reasonable approximation in the present context; note, also, that by forcing the fragmented parton to have a non-zero mass one is obliged to perform an arbitrary kinematic reshuffling. Therefore, in what follows hadrons will always be treated as massless. Finally, note that by inserting the identity d 3 K H δ( K H − ζ k p ) into the r.h.s. of eq. (4.3), and by using eq. (4.1), one readily obtains: The rightmost side of eq. (4.10) is by definition the expectation value of O, which shows again the correctness of eq. (4.2), and the way it must be understood.
In summary, we shall use eq. (4.2) as our master formula for the fully-differential implementation of single-FF cross sections. We remark that, at variance with eq. (4.1) whose l.h.s. is a function of a given four-momentum K H , eq. (4.2) is naturally associated with n different n-body final-state kinematic configurations, each of which corresponds to a different parton p undergoing fragmentation.

FKS: notation
While the MG5 aMC implementation of the FKS subtraction is done according to ref. [50], the notation used in that paper is not ideal to carry out formal manipulations such as those required to prove the cancellation of IR singularities, and certain aspects of the original FKS paper, ref. [46], have to be preferred. In particular, ref. [46] employs partonic cross sections that are highly symmetric in the kinematic and flavour spaces, which is a key feature we want to use here. Conversely, the S functions of ref. [46] are cumbersomethey include the definition of the observable jets and use it to partition the phase space, which was shown later in ref. [47] to be unnecessary; thus, we shall use here the S functions as defined in ref. [50]. The merging of the two notations is straightforward and we shall not explicitly spell it out here. We shall denote by eq. (II.n.m) eq. (n.m) of ref. [46]. We remind the reader (see sect. 3) that by eq. (I.x.y) we denote eq. (x.y) of ref. [50].
When one of the final-state partons is fragmented into a hadron, the short-distance cross sections have to be modified as discussed in sect. 4.1. Hence, we write the unsubtracted (n−1)-jet (m+2)-parton L-loop partonic cross section with one parton fragmented into a hadron as follows: Here, F H denotes a purely kinematical factor, responsible for enforcing final-state cuts on hadron H. While there will be no need to specify the functional form of this factor, its key property is that it must vanish when the transverse momentum of hadron H is below a given threshold. This implies that: (4.12) The factor J reconstructs n − 1 jets (the momentum of parton p being excluded from the reconstruction); as the notation suggests, an isolation condition (of the jets w.r.t. hadron H) is enforced as well. These two features are not mandatory, and are given here merely as examples; in particular, it is possible to include hadron H in the jet-finding procedure, and to obtain n final-state jets. The hadron-level contribution induced by eq. (4.11) is: where the sum over flavours runs over all possible QCD partons (i.e. the gluon and all the light quarks and antiquarks). We stress that the flavour of the fragmenting parton p is also summed over; this is obviously right from a physics viewpoint, and amends the notation used in sect. 4.1, where such a sum has been neglected since that case was relevant to a single partonic process. Loosely speaking, this sum restores the complete symmetry over final-state partons, which is originally broken by the fragmentation of one of them, and is thus ultimately responsible for the fact that the symmetry factor in eq. (4.11) is equal to 1/m!; more details on this point are given in appendix C. Since no confusion is possible, in eq. (4.13) we have employed the symbol to denote the convolution of the short-distance cross section with both the PDFs and the FFs. The former convolution has already been introduced in eq. (3.19), while the latter one is defined as follows: which has to be interpreted as explained in eq. (4.3). Note that by construction, the condition of hadron-level four-momentum conservation is: where the last term on the r.h.s. of this equation is the momentum associated with the fragmentation remnants. As is customary, we shall neglect to indicate explicitly any dependence on this momentum.

Final-state collinear counterterms
By means of the FKS procedure one arrives at subtracted partonic cross sections. In the case dealt with in ref. [46], the combination of all subtracted cross sections is still IR divergent, owing to the fact that the complete integration over initial-state degrees of freedom cannot be performed at the short-distance level (the incoming-parton momenta are constrained by their hadronic counterparts). In order to obtain IR-finite cross sections one must thus add the so-called (initial-state) collinear counterterms, whose forms can be derived from first principles following the procedure outlined at the beginning of sect. 2.1 of ref. [46] (see eq. (II.2.2)-eq. (II.2.5)).
In the fragmentation case one of the final-state momenta is constrained as well, and therefore final-state collinear counterterms (one for each fragmenting parton) must also be envisaged. To work out their forms, we start from eq. (4.2) and proceed by analogy to sect. 2.1 of ref. [46]. This implies that we need to replace hadron H with parton a q , and interpret the cross section on the r.h.s. of that equation as a subtracted one. As shown in eq. (4.13), we also need to consider explicitly the sum over the flavour of the fragmenting parton. Thus, the p-parton contribution to eq. (4.2) becomes: where a shortened notation is used. Now we expand the cross sections in eq. (4.16) by introducing their LO and NLO contributions, as is done in eq. (II.2.3): and we use the analogue of eq. (II.2.2): By replacing eqs. (4.17) and (4.18) into eq. (4.16) one obtains: Eqs. (4.19) and (4.20) are the analogues of eq. (II.2.4) and eq. (II.2.5), respectively. Note, however, that the role of the indices of the Altarelli-Parisi kernels is opposite in eq. (4.20) w.r.t. that in eq. (II.2.5), consistently with the fact that the former is a timelike branching, while the latter is a spacelike one. The second term on the r.h.s. of eq. (4.20) defines the final-state collinear counterterm we were seeking; as was said before, there will be one such counterterm per fragmenting parton.

Born and virtual contributions
As in the original FKS formulation, the Born and virtual cross sections follow trivially from the master equation, eq. (4.11) -we simply need to replace there the values m = n and L = 0 : to obtain the Born cross section, and the values m = n and L = 1 : to obtain the virtual cross section. Equations (4.21) and (4.22) are, as expected, extremely similar to their standard-FKS counterparts -eq. (II.2.6) and eq. (II.2.15) respectively, the only formal differences being due to the fragmenting parton and its associated factor. In both of these equations 3 ≤ p ≤ n + 2, owing to eq. (4.13).
Equation (4.22) implies that the IR divergences of virtual origin in the fragmentation case are given by eq. (II.3.2), times the fragmenting-parton factors. Therefore, the cancellation of the singularities is proven if one shows that all of the contributions of real-emission origin have the same forms as those given in ref. [46], times the fragmentation factors. It is obvious that this is a sufficient but not necessary condition (since the individual realemission contributions might have different forms w.r.t. the standard-FKS ones, with only their sum being equal to the sum of the latter, up to the fragmentation factors); but it will turn out that by and large this is indeed the case, which is a significant simplification.

Real-emission contribution
We shall follow here as closely as possible the procedure outlined in sect. 4 of ref. [46]. We split the real-emission contribution into a soft and a non-soft part, as is done in eq. (II.4.11) and eq. (II. 4.12). For the latter, the initial-and final-state collinear contributions are dealt with separately, by using a decomposition identical to that which leads to eq. (II.4.14) and eq. (II.4.15). By exploiting eq. (4.12), the complete symmetrisation over final states, and the invariance over parton relabeling (see appendix C), the results of ref. [46] can straightforwardly be extended to the fragmentation case for the soft and the initial-state collinear contribution. In particular, the analogue of eq. (II.4.27) reads: The quantity in eq. (4.25) is a genuine (n + 1)-body term, and thus 3 ≤ p ≤ n + 3. On the r.h.s. of that equation factors appear that stem directly from the FKS representation of the (n + 1)-body phase-space in 4 − 2 dimensions (see eq. (II.4.4)), namely: Thus, in eq. (4.25) dφ ( =0) is such that the measure on the r.h.s. times the volume factor V ( =0) is equal to dφ n+1 /ξ i in four dimensions. The initial-state collinear remainders of eqs. (4.26) and (4.27) are (quasi-)n-body terms, and therefore 3 ≤ p ≤ n + 2.
As was anticipated in sect. 4.4, eqs. (4.23), (4.26), and (4.27) have the same IR-singular structure as their counterparts in ref. [46] (bar for the fragmentation factors). By showing that this is the case also for the contributions of final-state collinear origin, we shall achieve the sought proof of the cancellation of the IR divergences in the presence of fragmentation.
In order to do this, we start from the analogue of eq. (II.4.15), which reads as follows: where The first steps are the same as in ref. [46], and in particular one can follow them up to eq. (II.4.63) in order to achieve the splitting into a singular and a non-singular part that appears in eq. (II.4.64). The finite part, eq. (II.4.65), can directly be used also in the fragmentation case, through the usual change of notation and the multiplication by the fragmentation factor. It reads: i dξ j dy j dΩ (2) j dφ ( =0) .
This is a genuine (n + 1)-body term, and thus one has 3 ≤ p ≤ n + 3. By construction, dφ ( =0) is such that the measure on the r.h.s. times the volume factor V 2 ( =0) is equal to dφ n+1 /(ξ i ξ j ) in four dimensions.
As far as the singular part is concerned, it requires a more careful treatment than the finite one, which constitutes the only non-trivial bit of the proof given in this section. To this end, we need to distinguish two cases in eq. (4.31): A. p = i and p = j; In words, either the fragmenting parton is different from both partons associated with the S ij functions (case A.), or it is equal to either of them (case B.). Let us first simply count the number of contributions, and make sure that cases A. and B. exhaust all possibilities. The correct counting emerges by considering both the sum over all possible S functions, and that over the fragmenting partons. Thus, in the present case: (1 − δ ij ) = (n + 1)(n + 1)n , (4.34) where the δ ij term enforces the condition i = j that stems from the S functions. By direct computation: given that a term δ ip δ jp must vanish if i = j. Thus: The first term inside the square brackets on the r.h.s. of eq. (4.36) corresponds to case A., while the sum of the other two terms corresponds to case B. The number of contributions they are associated with, defined in analogy to eq. (4.34), is: The fact that S = A + B follows from eq. (4.36), and can be checked directly from the r.h.s.'s of eqs. (4.34), (4.37), and (4.38). We then go back to the treatment of final-state collinear singularities in standard FKS (sect. 4.4 of ref. [46]). The idea is that partons i and j (associated with S ij ) are combined into what is called, in ref. [46], parton 7 (here it will be n + 4) in the intermediate steps of the computation. Thus, one arrives at eq. (II. 4.83) where, similarly to the cases of soft and initial-state collinear singularities, there is no dependence upon i and j in the reduced cross sections. Therefore, all S ij contributions are identical after relabeling. The number of such contributions is: (1 − δ ij ) = (n + 1)n . (4.39) When combined with the real-emission symmetry factor, this gives: 1 (n + 1)! (n + 1)n = n n! . (4.40) The factor 1/n! is then included in the Born-level cross sections; the factor n at the numerator is cancelled by turning the identical contributions into a fully-symmetric sum over final-state partons implicit in eq. (II.4.86). That sum is what appears as sum over j in eq. (II.4.88). Now consider case A. Owing to its definition (the Kronecker delta's in eq. (4.37)), which implies that no dependence upon i and j enters the fragmentation part, the manipulations of sect. 4.4 of ref. [46] do apply to the present case as well, up to eq. (II. 4.83). At this point, however, the fragmentation factor interferes with the counting of contributions and the relabeling. Eq. (4.39) must be replaced by eq. (4.37) and therefore, rather than eq. (4.40), one has: 1 (n + 1)! (n + 1)n(n − 1) = n(n − 1) n! . (4.41) After including the factor 1/n! into the Born cross sections, one is therefore left with n(n−1) contributions, that must account for the sum over the fragmenting partons and the sum over j on the r.h.s. of eq. (II.4.88). The fact that this is really so stems from observing that the definition of case A. implies that, after relabeling, the index of the fragmenting parton can take any value different from that of the parton associated with the colour factor Z(a) (see appendix A of ref. [46]), which is j in eq. (II.4.88). Furthermore, we must take into account that, after relabeling, the index p of the fragmenting parton is such that 3 ≤ p ≤ n + 2. These two facts amount to saying that the sum over fragmentation contributions and a fully symmetrised final state must be: (1 − δ jp ) = n(n − 1) , (4.42) which is exactly the numerator of eq. (4.41). This proves that in case A. the FKS result of eq. (II.4.88) still holds (with the usual changes of notation), with an extra factor of fragmentation origin: that must be inserted under the two sums over j on the r.h.s. of that equation, for any given fragmenting parton p. Eq. (4.43) explicitly enforces the conditions that all partons except j may fragment. Let us now turn to case B. To this end, consider the singular part that emerges from eq. (4.31) after the decomposition analogous to that of eq. (II.4.64); the quantity of interest will be (see eq. (4.36) and the comment that follows it): (4.44) This can be further written as (see eq. (II.4.66)-eq. (II.4.69)): where D i is given in eq. (4.32) and (4.47) As the notation suggests, in eq. (4.46) the terms A i and A j are proportional to F The indices α and β are such that {α, β} = {i, j}, with the two possible choices denoted as follows: so that: eq. (4.53) coincides with eq. (II.4.80). Furthermore, by using eq. (I.4.17), eq. (II.4.58), eq. (II.4.71), and eq. (II.4.72) (and by neglecting the ∆ term in the latter, owing to its vanishing upon azimuthal integration), we obtain: × P < aαS(aα,a β ) (z, ) M (n,0) (S(a α , a β )) h(1 − z)F By construction, the reduced matrix element M (n,0) features the parton that splits into the pair (i, j), whose flavour is equal to S(a α , a β ), and which is explicitly indicated in eqs. (4.54) and (4.55) -the dependence on the other partons, being irrelevant in what follows, is omitted. Note that some minor abuse of notation concerns the fragmentation factor, since according to the general notation introduced in sect. 4.2 the argument of F H (zk n+4 ). However, given that in what follows only the dependence upon z matters, that on k n+4 is dropped from the notation. Finally, observe that the Altarelli-Parisi kernels in eqs. (4.54) and (4.55) are identical; this is because the collinear limit of the (n + 1)-body matrix elements can be written (in a fully equivalent manner) by factoring either of the two following kernels: We have used the former for the change of variable c i (applied to the term A i ), and the latter for the change of variable c j (applied to the term A j ). Conversely, the arguments of the function h are different, since that function results from the collinear limit of the S ij function, which has the unique form given in eq. (I.4.17).
In the collinear limit enforced by δ(1 − y j ), the expressions of the measures given in eqs. (4.52) and (4.53) simplify as well. Firstly, observe that, thanks to eq. (II.4.71), in such a limit the polar and azimuthal angles of parton α in eq. (4.52) (parton β in eq. (4.53)) coincide with those of the splitting parton (labeled by n + 4 here, and by 7 in ref. [46]). Thus: where we have used eqs. (4.28) and (4.29). The quantity: is the reduced n-body phase space, which originates from the (n + 1)-body one by "combining" (in the collinear limit) momenta k i and k j into the momentum k n+4 of the mother parton; note that eq. (4.58) coincides with eq. (II.4.82). We can now observe that in the collinear limit the dependence on the azimuthal variables of parton β are trivial (see eq. (4.54)); therefore, we can integrate over them and use eq. (II.4.42), whence eq. (4.57) becomes: By performing the same manipulations (and by taking into account that the roles of α and β are swapped) we obtain from eq. (4.53): Let us now first consider the contribution Y j to eq. (4.45); we can use the results of eqs. (4.55) and (4.60), use the Dirac delta to get rid of the y α integration, and obtain: where we have used the fact that: for any function f (). Owing to the fact that the change of variables c j of eq. (4.51) is identical to that of eq. (II.4.70), the quantity (ξ n+4 D β ) that appears in eq. (4.61) is rewritten by using eq. (II.4.76)-eq. (II.4.79): Now we consider the contribution Y i to eq. (4.45); by proceeding analogously to what was done above, using eqs. (4.54) and (4.59), we obtain: It is important to note that the change of variables c i of eq. (4.50) is not the same as that of FKS. Hence, the identity of eq. (4.63) must not be used to handle the term (ξ n+4 D α ) that appears in eq. (4.64). It is clear than an analogous expression exists that can be applied to the present case; however, this is not necessary. The reason is the following: thanks to the presence of h(z)F  H (z)) and at z = 1 (because of h(z)). Therefore, the subtraction terms in D α are identically equal to zero, which implies that D α is a regular function (and not a distribution), whence: By substituting this expression into eq. (4.64) we obtain: Eq. (4.67) is fairly similar to eq. (4.61). One can in fact show that such a similarity is even closer than it appears at a first look. One starts by observing that: Then, one uses again the fact that h(z) damps singularities at z = 1. Thus: for any F (which is in general a (4 − 2 )-dimensional coefficient). The quantity in curly brackets on the r.h.s. of eq. (4.69) has the same form as the r.h.s. of eq. (4.63), and can be made identical to it by suitably choosing the coefficient F (namely, by setting it equal to the coefficient of δ(1 − z) in eq. (4.63)). This implies: Thus, eq. (4.67) becomes: In this way, it is apparent that Y i in eq. (4.71) has exactly the same form as Y j in eq. (4.61) (because of eq. (4.63)), except for the fact that the former features a factor h(z), while the latter features a factor h(1 − z). This implies that the sum Y i + Y j , which is the quantity of interest (see eq. (4.45)), will feature: thanks to eq. (I.4.23). After some algebra, and by using eq. (4.30), one arrives at: Note, for future reference, that: with the pre-factor on the r.h.s. being the one that is conventionally factored out in FKS formulae (see e.g. eq. (II.3.1) and eq. (II.4.88)). Eq. (4.73) is the analogue of eq. (II.4.83), and it is apparent that these two equations are rather similar. In order to proceed, one needs to perform the same operations as in eq. (II.4.86): namely, to sum over all (i, j) pairs, and to relabel the parton indices. We have already gone through this process earlier in this section, when we have dealt with case A. In the case B. we are discussing here, the number of contributions is given by the r.h.s. of eq. (4.38), divided by two. This is because in eq. (4.38) the overall factor of 2 comes from the fact that the results stemming from the δ ip and δ jp terms are identical. However, the quantity Y ij already includes both of them (see eq. (4.44)), whence the necessity of dividing by two the r.h.s. of eq. (4.38) for a correct counting. The bottom line is that the number of contributions is the same as in FKS (eq. (4.39)), and we can therefore proceed as done there as far as the symmetry factors are concerned. Note that, at variance with case A., in case B. the sum over fragmenting partons is no longer present.
There is a further subtlety associated with relabeling that we must mention. The fragmenting parton α and its mother parton n+4 have degenerate momenta in the collinear limit, k α = zk n+4 . In eq. (4.73) there is a sum over the flavour (d) of the mother parton, which of course does not affect the kinematics: k n+4 is the same for all terms in such a sum. By construction (see eq. (4.58)) this momentum is generated in the reduced phase space dφ n . Finally, keep in mind that the rescaled momentum ξ n+4 enters the expressions of D (0) (z) and D (1) (z) (see eq. (4.63)). Upon relabeling, while the flavours of the fragmenting parton and of the mother will need to be kept distinct, one can assign the same label to these two partons, thus resulting in a more compact notation (no confusion being possible) 9 .
At the end of the day, we can therefore write: where one understands that the expressions of D (0) (z) and D (1) (z) are given in eq. (4.63) with the formal replacement there. From this point onwards, k k is the momentum of the mother parton (because is the one that appears in the phase space), with the momentum of the fragmenting parton being equal to zk k by construction.
It is manifest that eq. (4.75) does not have the same singular structure as eq. (II. 4.88). This is what one expects, because we have so far neglected the contribution of the finalstate collinear counterterms. These are given in the second term on the r.h.s. of eq. (4.20), and we simply have to write them in the same form as eq. (4.75), which is straightforward to do. Firstly, we observe that the counterterm of eq. (4.20) is relevant to a single parton that fragments; in keeping with eq. (4.13), we have therefore to sum over all possible fragmenting partons, as well as over their flavours. Secondly, eq. (4.20) originates from eq. (4.16), where the fragmented parton momentum is fixed, and the fragmenting parton momentum is obtained from the former by a 1/ζ rescaling. As discussed above, eq. (4.75) is such that one fixes the parent momentum; thus, we shall use eq. (4.20) by formally replacing k → ζk (we shall also rename ζ as z, for consistency with eq. (4.75)). In this way, the final-state counterterms to be added to eq. (4.75) read as follows: Note that, precisely as in the case of initial-state collinear singularities, the AP kernel in eq. (4.77) is not the same as that in eq. (4.75): the former is in four dimensions but includes the contribution at z = 1, while the latter is in 4 − 2 dimensions but lives in the z < 1 space. Thus, we shall now manipulate eq. (4.75) in order to show that it features a singularity which has the same form as that in eq. (4.77). By expanding eq. (4.75) in , and using eq. (II.4.52) one obtains: By exploiting the identities given in eq. (II.4.48) and eq. (II. 4.49), and using the definition of D (0) (z) (see eq. (4.63)), we have: By multiplying both members of this equation by the factor 1/ that appears in eq. (4.78), and by using the definition of (see eq. (4.74)), one obtains: By replacing the r.h.s. of this equation into eq. (4.78), one sees that the pole term proportional to the AP kernel has exactly the same form, except for an overall sign, as that in eq. (4.77), and therefore drops in the sum. Thus: Let us consider eq. (4.81). The sum over a k is what remains of the sums over a i and a j in eq. (4.44). According to the definitions of parton-and hadron-level cross sections given in eqs. (4.11) and (4.13) respectively, it belongs to the latter; hence, in the parton-level quantity defined by eq. (4.81) one fixes a k and does not sum over it. It is then clear that, apart from the fragmentation factor F (k) H (1), the divergent parts of eq. (4.81) and of eq. (II.4.88) are identical. As far as the fragmentation factor is concerned, we can employ the identity: When the r.h.s. of this equation is substituted into eq. (4.81), the divergent part of the latter is complementary to the divergent part of the contribution of case A. (compare eqs. (4.43) and (4.82)). Thus, regardless of the value of p (the label of the fragmenting parton) the union of case A. and case B. is identical to the divergent part of eq. (II.4.88), up to the fragmentation factor. This is therefore the same structure as that of all of the other NLO contributions; we have therefore proven that the cancellation of the IR singularities is fully achieved, and that it proceeds similarly to what happens in the standard FKS case. In summary, the contribution of final-state origin can be written as follows: The first term on the r.h.s. of eq. (4.83) is the analogue of eq. (II.4.65) and is given by eq. (4.33). The other two terms emerge from the manipulations performed above, and can obtained from case A. and eq. (4.81). Remember that they are n-body terms, and thus 3 ≤ p ≤ n + 2. The quantity dσ  H (1) in eq. (4.81); some are already explicit, but others have to be obtained from the δ(1 − z) terms contained in D (0) (z) and D (1) (z)in order to do so, use has to be made of eq. (II.4.48) and eq. (II.4.49). Finally, the terms not proportional to δ(1 − z) in eq. (4.81) are collected into dσ (n+1,0;p) out,+,z . After some algebra, and by using eq. (4.76), the results read as follows (3 ≤ p ≤ n + 2): and: As was already said, the divergent part of eq. (4.84) is identical, up to the fragmentation factor, to that of eq. (II.4.88). As far as the finite part is concerned, the log Q 2 term is again identical to that of eq. (II.4.88); this must be so, since the Ellis-Sexton-scale dependence of the present contribution must be compensated by that of the soft and virtual ones, which in turn depend on the presence of fragmentation factors only in a trivial manner. The C(a k ) bit is also identical to that of ref. [46], while the γ(a k ) and γ (a k ) ones are not (owing to the conditions enforced by δ kp ). A simple explanation for this fact is that the γ's are "integrated" quantities: in the absence of fragmentation a full integration over the i and j degrees of freedom is possible, but this is not the case when fragmentation is present, and thus there is an amount of talk-to between eqs. There is an extra logarithmic term in eq. (4.85) which features the energy (zE p ) of the fragmenting parton: this is the related to the "missing" γ and γ contributions in eq. (4.84) to which we have alluded above. We shall soon return to comment on this extra logarithmic term. Before doing that, it is convenient to re-cast eq. (4.85) in a slightly different form, so that its similarities with eq. (II.4.54) and eq. (II.4.55) can be made more explicit. We observe that, while in the initial-state case we have used distributions whose subtractions are controlled by the parameter ξ c , the final-state cross section is so far expressed in terms of standard plus distributions. We can amend the situation by means of the identities for any 0 ≤ z c < 1, and having defined: The replacement of eqs. (4.86) and (4.87) into eq. (4.85) generates terms of the same form as those that appear in eq. (4.84) (i.e. that factorize the Born without any z convolution), and are proportional to C(a p ) owing to eq. (II.4.49); thus, they can conveniently be included in the definition of dσ (n+1,0;p) out,+,δ . In order to write them explicitly, we choose: Some trivial algebra shows that dσ (n+1,0;p) out,+,δ now reads as follows: The structure of eq. (4.90) is interesting. Its finite part for k = p is the same as before, i.e. identical to its counterpart in eq. (II.4.88) or eq. (I.4.6). Conversely, the finite contribution associated with the fragmenting parton features a log µ 2 /Q 2 term which is identical to an analogous piece of initial-state origin (see eq. (II.4.54), eq. (II.4.55), and the first line of eq. (I.4.6)). It is the other term that seemingly does not have an initial-state analogue. However, as for the logarithmic term in eq. (4.85) or eq. (4.91) which features the energy zE p , this is an accident due to kinematics, and in particular to the fact that FKS cross sections are written in the rest frame of the incoming partons, where the energies of the latter are E 1 = E 2 = √ S/2. We have in fact explicitly verified that, by repeating the same computations that has led to eq. (II.4.54) and eq. (II.4.55) without assigning a specific value to E i , i = 1, 2, we arrive at the same results as in ref. [46], bar for two extra terms. The first of these reads as follows: and must be inserted in the curly brackets on the r.h.s. eq. (II.4.54) and eq. (II.4.55), thus rendering these two equations fully analogous to their final-state counterpart, eq. (4.91). The second term is such that the finite part of the Born-like term in eq. (II.4.54) and eq. (II.4.55) reads: which is identical (up to the obvious replacements i → p and δ I → δ O ) with the one that appears in the last two lines of eq. (4.90). As they must, eqs. (4.92) and (4.93) vanish by setting E i = √ S/2, thus recovering the results of ref. [46]. We have therefore found that the FKS short-distance cross section of initial-and finalstate collinear origin have a larger degree of similarity in the presence of fragmentation than in the inclusive case. This need not be surprising, since it is only in the former case that both initial-and final-state partons are kinematically constrained by "external" objects (the PDFs and FFs, respectively). Having said that, observe that by construction that scale µ that appears in eqs. (4.84) and (4.85) (or eqs. (4.90) and (4.91)) is the scale that enters the FFs, which in principle can be different from the factorisation scale used in the PDFs.

Summary
In this section, we have revisited the FKS subtraction by applying it to the cases in which one final-state parton is fragmented. By following as closely as possible the procedure of the original paper [46], we have proven the cancellation of the IR singularities that emerge in the intermediate steps of the computations, and obtained in the process the IR-finite short-distance cross sections. The forms of the latter are rather similar to those relevant to the inclusive (i.e. not fragmented) case, with the most significant differences due to the contribution of final-state collinear origin.
The origin of the final results, and the formulae where they appear, are the following: Born (eq. (4.21)), virtual (eq. (4.22)), soft (eq. (4.23)), subtracted initial state (eq. (4.25)), initial-state collinear remainders (eqs. As a check of the correctness of the IR-finite parts of the above formulae, we have applied them to the computation of the initial conditions of the perturbative b-quark fragmentation function, finding full agreement with ref. [72] (note the erratum of that paper).

Introduction
A long-standing difficulty of perturbative computations of scattering amplitudes in Quantum Field Theories is that of the handling of unstable short-lived particles. The situation is complicated by the fact that contributions from such particles can spoil gauge invariance and unitarity, and their treatment possibly necessitates relaxing strict fixed-order accuracy. In particular, when carrying out NLO EW computations in the SM, care is needed when considering processes that involve the intermediate massive vector bosons Z and W ± , the Higgs boson, and the top quark.
A Feynman diagram that features the propagator of an unstable particle P with mass M and virtuality p 2 diverges when p 2 → M 2 . In those cases where p 2 is associated with a physical observable (e.g. the e + e − invariant mass in a Z decay) such a divergence can be avoided by means of final-state cuts, but this is not desirable in general, since it prevents one from studying the kinematically-dominant pole region. Thus, a universal solution entails a regularisation of the propagator. The ones which are best motivated from the physical viewpoint all stem from the Dyson summation of the geometric series that results from the insertion of two-point 1PI graphs; this leads to the propagator 10,11 : with M 0 the bare mass of P , and Σ(p 2 ) equal (possibly up to a factor of i) to the 1PI contribution. Because of the presence of an imaginary part in the denominator, eq. (5.1) achieves the sought regularisation. Unfortunately, one cannot naively Dyson-sum all P propagators in an arbitrarily complicated Feynman diagram. Among other things, by effectively setting the pole mass of the propagator of P to be different from the corresponding Lagrangian mass parameter, one might violate gauge invariance [79][80][81][82]. However, eq. (5.1) does suggest that a better-defined mass (and width) should indeed be associated with the position of the complex polep 2 : In fact, note that the appeal of eq. (5.2) is that it reads the same when expressed in terms of either bare/unrenormalised quantities or renormalised ones: the position of the pole must not change under renormalisation, lest the analiticity properties of the S matrix be changed too (see e.g. chapter 10 of ref. [83]). Equation (5.1) also suggests an alternative way (still potentially gauge-violating) in which the P propagator can be regularised, namely: with M the on-shell mass, Γ the total decay width of P , and where one exploits the fact: which follows directly from the optical theorem. At the level of squared amplitudes, eq. (5.3) will lead to a Breit-Wigner (BW) form: In the timelike region p 2 ≡ s > 0, i.e. when P becomes resonant at s M 2 , the term Γ 2 M 2 prevents BW (s) from diverging. Furthermore, a BW function admits the following expansion in terms of distributions: 10 We simplify the notation by considering only the case of a scalar particle. It should be clear that, for the arguments of this section, this implies no loss of generality. 11 Throughout this section we understand the +i0 prescription of the Feynman propagators.
where by P we denote the principal-value operator. Through eq. (5.6) one can study the impact of off-shell effects (where s = M 2 ), and address gauge violations. For example, by keeping the first term on the r.h.s. of eq. (5.6) one works in the well-known (and gauge invariant) narrow-width approximation. Equation (5.3) is also employed in the context of the so-called pole approximation [79,80,84,85], where a suitable subset of amplitudes (typically associated with NLO corrections), stripped of the propagator of eq. (5.3), are expanded around p 2 M 2 ; the leading term of such an expansion, which can be shown to lead to a gauge-invariant result, is then kept. A systematic generalisation of the pole approximation can be achieved within the unstable-particle effective theories [86][87][88], which at one-loop accuracy are equivalent to the former. Here, one exploits the fact that at high energies there is a natural hierarchy established by the relationship E 2 M Γ, which allows one to separate the production and decay mechanisms, and to achieve a formal expansion in Γ/M in the pole region p 2 M 2 .
Neither the pole approximation nor the unstable-particle EFTs are apt to study the offshell region, where non-resonant contributions might become important, and the transition between the off-and on-shell regions. Furthermore, automated calculations are essentially impossible in the former approach, and would be highly non-trivial in the latter (beginning with the fact that the construction of the underlying theory model is a very involved operation). Conversely, these issues are absent if one works in the complex mass (CM henceforth) scheme [48,49]. Thus, the CM scheme is the strategy of choice in MG5 aMC for dealing with unstable particles, and we shall limit our discussion to it in the present paper.
The core idea of the CM scheme, which stems from the observations related to eq. (5.2), is to modify the renormalisation conditions of the theory, yielding complex-valued renormalised parameters that include the masses of the unstable particles, and also potentially a subset of the coupling constants. By construction, this procedure leaves unaltered all algebraic relations realizing gauge invariance. At the same time, it naturally regularises unstable-particle propagators, which assume the same functional form as in eq. (5.3).
At the LO, such a modification is rather innocuous, since the analytic continuation of LO amplitudes, which involve only rational expressions of kinematic invariants, is unambiguous. At the NLO, however, the extension of the on-shell renormalisation condition [89] presents many subtleties in Quantum Field Theory, such as the proof of perturbative unitarity investigated in ref. [90,91], and its applicability beyond NLO. In what follows, we shall limit ourselves to discussing the more pragmatic concerns of assessing the correctness of the CM implementation in MG5 aMC, as well as listing and providing the necessary ingredients for guaranteeing the formal NLO EW accuracy of our results.

Complex mass scheme formulation
Let us first point out that when referring generically to the CM scheme we actually understand the common properties of a class of schemes, whose members differ in the choice made for the independent input parameters. This implies that, when only a specific member of this class is relevant, we may characterise it more precisely as, for example, CM G µ scheme, or CM α(m Z ) scheme, and so forth.
As was anticipated in sect. 5.1, at the LO the use of the CM scheme simply amounts to redefining the mass of all unstable particle fields, according to eq. (5.2): All derived parameters must then be expressed in terms of the complex masses of eq. (5.7) and of the independent inputs, thereby possibly acquiring an imaginary part. In the SM, this is for example (but not exclusively) the case of G µ (in the α(m Z ) scheme), of α (in the G µ scheme), of the cosine of the Weinberg angle, and of the Yukawa couplings.
At the NLO, one needs to properly define the renormalisation conditions. Given the similarities between the CM scheme and the standard on-shell (OS henceforth) one, we start by recalling those that in the latter are relevant to the self-energies. We start by introducing the mass and wave-function counterterms: where we have denoted by Σ R and Σ U the renormalised and unrenormalised self-energy, respectively, of the unstable particle 12 . Since in the OS scheme the renormalised mass must remain real, the renormalisation conditions need only to take into account the real part of the self-energy 13 : which yield the following expressions for the counterterms: where we have employed the usual shorthand notation: for any q 2 . We point out that the real part that appears in eqs. (5.9)-(5.12) is trivial for stable particles (since for them the self energies are real quantities), but is necessary when dealing with unstable particles, and thus encompasses all situations. Indeed, for unstable particles 12 Where possible, we shall adhere to the notation convention established in eqs. (5.7) and (5.8), within which lowercase (uppercase) symbols are associated with Lagrangian parameters in the CM (OS) scheme. Note that, in general,M = M . 13 See footnote 11, which is especially relevant to the evaluation of self energies on the mass shell.
the one-loop Σ develops an imaginary absorptive part that can be related, through the optical theorem, to the total LO decay width (see eq. (5.4)), and this happens independently of whether the free propagator is assigned a zero or a non-zero width. It is then the zerowidth case that implies that the use of the real part is not academic, because the OS scheme can be employed to carry out NLO calculations 14 for processes that either feature only stable particles, or where unstable particles are present but finite-width effects can be neglected -the typical situation being that of a non-resonant t-channel exchange, for which the renormalisation conditions of eqs. (5.11) and (5.12) are directly relevant (conversely, non-vanishing widths pose the problems already discussed in sect. 5.1, thus rendering the OS scheme an option not viable in general in this case). We shall further discuss the consistency of renormalisation procedures when setting particle widths equal to zero in sect. 5.5. We now turn to the CM scheme, which allows one to address the finite-width scenario. Fundamentally, in such a scheme one imposes on-shell-type renormalisation conditions that involve both the real and the imaginary part of self-energies. The ensuing UV counterterms feature an imaginary component that order by order helps avoid any double counting between the non-zero width and the bubble insertions in internal propagators. We write the analogues of eqs. (5.8) as follows, by taking eq. (5.7) into account: One then generalises eqs. (5.11) and (5.12) as follows: These definitions, together with those relevant to coupling renormalisation (which we do not report here, owing to their being functionally identical to those relevant to the OS scheme), ensure that by working in the CM scheme one can proceed analogously to what is done in other renormalisation schemes. For example, we observe that the renormalised Lagrangian expressed in terms of renormalised parameters is equal to the bare Lagrangian expressed in terms of bare parameters. Furthermore, while all the derived parameters might acquire an imaginary part as was already the case at the LO, they do so without spoiling the relationships among them which are constrained by gauge invariance. Finally, we note that [79,82]: By writing the rightmost equality in eq. (5.17) we have assumedΓ = O(α), which may appear in contradiction with the fact thatΓ is regarded as an independent input parameter.
14 The explicit expressions for the NLO OS counterterms relevant to the EW sector of the SM can be found e.g. in ref. [89].
What we understand here is that such an input is associated with a given value of α (in general, the value of the independent coupling relevant to the chosen scheme). If we were to measureΓ for progressively smaller values of α, we should obtain a series tending linearly to zero. Likewise, this is the meaning we associate with saying that the difference between M 2 andM 2 is of O(α 2 ). Given the actual very slow running of α, these remarks do not play a major role in physics simulations, but must be taken into account when considering some of the mathematical properties of the CM scheme, which we have to exploit in its MG5 aMC implementation, and which we now turn to discussing.
In what follows, we shall need to consider some α → 0 limits. It is therefore convenient to assume to work in the CM α(m Z ) scheme, where α is real. There is no loss of generality in this, since in other CM schemes one would instead consider the limit in the real-valued input coupling from which α is derived. By using eqs. (5.14) and (5.15), one obtains: We can now exploit eq. (5.17), the optical theorem (eq. (5.4)) applied to the one-loop self-energy, and the fact that at this order the imaginary part of such self-energy is finite. By expanding the rightmost term of eq. (5.18) around p 2 = M 2 , and by keeping only the dominant term, we obtain: with Γ (0) the LO total decay width. Eq. (5.19) complements our previous comment about Γ. Namely, in the context of numerical tests not only this quantity must vanish linearly with α, but also it must do so in such a way to guarantee that: We shall show the importance of eq. (5.20) in appendix E.1, where we shall present a technique for some systematic and automated tests of CM-scheme implementations. Equation (5.17) and the "perturbative" vanishing ofΓ are also directly relevant to the proper definition of the mass and wave-function counterterms, eqs. (5.15) and (5.16)here, we shall consider only the former in order to give a definite example. The self-energy Σ typically features branch cuts in the p 2 complex plane; hence, its analytical continuation must be performed with care, and we shall discuss this issue in sect. 5.3. Here, we limit ourselves to point out that, owing to the strict connection between the OS and CM schemes, a natural cross check on the final result is to verify that, in the α → 0 limit, the mass counterterms in the OS scheme (eq. (5.11)) and in the CM scheme (eq. (5.15)) be equal: It is crucial to understand that, when taking this limit, the propertiesM → M andΓ → 0 must be enforced. Eq.  Contrary to what is customary done, we do not circumvent this issue by Taylor-expanding Σ, but instead evaluate the complete self-energies in the appropriate Riemann sheet. We conclude this section by mentioning the fact that, as far as the SM is concerned, all renormalisation counterterms have been derived and implemented by hand in two UFO models (distributed with the MG5 aMC release), relevant to the CM α(m Z ) and CM G µ schemes (see sect. 5.4 for further details), respectively. Thus, as is the case for all NLOgrade UFO models, these include in particular all UV and R 2 counterterms necessary for MadLoop to compute an arbitrary one-loop amplitude 15 .

Definition of mass and wavefunction UV counterterms in the CM scheme
As was just discussed in sect. 5.2, UV-renormalisation conditions in the CM scheme imply, in particular, the necessity of evaluating massive-particle self energies at p 2 =M 2 − iΓM (see e.g. eqs. (5.15) and (5.16)). Thus, owing to the presence on the extra (w.r.t. that of the OS scheme) imaginary part −ΓM in this mass-shell condition, analytical continuation might lead one to compute logarithms (or any other multi-valued function) in Riemann sheets different from the first. The goal of this section is to discuss the strategy put in place in MG5 aMC in order to tackle this problem.
One example relevant to the computation of a logarithm can be readily given by considering the one-loop fermion contribution to the W self energy in the SM, depicted in the left panel of fig. 1, whose unrenormalised transverse part reads as follows on the CM-scheme mass shell 16 : Here, s W is the complex-valued sine of the Weinberg angle, N f c is the number of colours relevant to the SU c (3) representation to which fermions f 1 and f 2 belong, and by µ we denote the mass scale that needs to be introduced in the context of dimensional regularisation. We have also employed the logarithm in the second negative Riemann sheet, defined according to the general formula: where by log 0 z we have denoted the principal-value or first-Riemann-sheet definition of the logarithm with a branch cut on the negative real axis: whence: Since in practice we shall end up using only either first or second Riemann sheet logarithms, we shall adopt the shorthand notation: The use of log − () in eq. (5.22) instead of log() guarantees that the result for Σ f 1 f 2 U,T is obtained consistently with the Feynman-propagator prescription it stems from. Ultimately, together with similarly consistent treatments of the other contributions to the self-energy (which we shall discuss below), this leads to eq. (5.21) being fulfilled.
In general, the situation is complicated by the possible presence of multiple scales and/or unstable particles that circulate in the loops. The solution proposed in ref. [49] is to Taylor-expand in p 2 aroundM 2 any self-energy Σ(p 2 ) relevant to the computation of UV counterterms. Given that the latter are evaluated on the complex-mass pole of eq. (5.7), this would be essentially identical to an expansion inΓ aroundΓ = 0, were it not for the possible dependence onΓ due to sources different from the mass shell of the particle whose self-energy is being computed (we shall give below an example of this). We point out that, in order to be NLO-accurate (assuming that O(Γ/M ) = O(α)), at least the first two terms of the Taylor expansion must be retained, in order to get rid of a contribution ∼ 1/(iΓM ) due to an intermediate propagator being on-shell in the resonant region.
The Taylor-expansion technique has some tricky aspects, in particular due to the possible sensitivity to soft kinematics of loop propagators. In order to illustrate such aspects, we consider another contribution to the transverse part of the W self energy, namely that due to a (γ, W ) loop (see the right panel of fig. 1). For our current purposes, it is sufficient to deal with the B 0 part of such a contribution: We remind the reader that the second and third arguments of the B 0 function correspond to the masses of the two particles that circulate in the loop. By using the explicit expression for B 0 one obtains: Conversely, the first-order Taylor expansion 17 Equation (5.30) explicitly shows one of the features of the Taylor expansion mentioned above, namely that it is not exactly equivalent to expanding inΓ W , owing to the remaining dependence upon this quantity in the third argument of B 0 . Equations (5.28) and (5.30) lead to: The difference in eq. (5.31) is thus of NLO, which implies that the Taylor expansion of eq. (5.29) is not sufficient to obtain a correct result at this perturbative order. The problem stems from the derivative of the logarithm that appears in the third line of eq. (5.28), that induces a contribution proportional to 1/(Γ WMW ) when p 2 =M 2 W , and ultimately from the fact that, in theΓ W → 0 region relevant to the Taylor expansion, one becomes sensitive to the branch point of the logarithm. This is also the reason why the situation does not change if considering higher-order terms in the Taylor expansion. Indeed, one can show that an n th order expansion leads to: The problem exposed above has indeed been pointed out in ref. [49], and a pragmatic solution proposed there is that of adding the "missing term" (i.e. the r.h.s. of eq. (5.32)) back to the expanded counterterms. It is clear that a straightforward Taylor expansion works for all graphs where one does not cross branch cuts of logarithms and other multivalued functions whenΓ → 0. On the other hand, such an expansion requires thatΓ can be viewed as a small parameter (which is not necessarily equivalent to saying that . This is certainly true in the SM, but not necessarily so in arbitrary newphysics theories.
The obvious alternative to Taylor-expanding is to keep the full self-energy expressions, and to figure out the appropriate analytical continuation of the OS results which becomes necessary when imposing CM-scheme renormalisation conditions. This is equivalent to being able to choose the appropriate Riemann sheets where the multi-valued functions that appear in the OS results are to be computed. Such an approach has the advantage of being immediately applicable to any BSM theory, regardless of its mass spectrum and width settings.
In order to pursue this strategy, we start by reminding the reader that for the computation of mass and wavefunction UV counterterms we are concerned only with 1-point and 2-point scalar integrals, and the former do not pose any problems. As far as the latter are concerned, the basic results we need to consider are the following [92] (see also ref. [93]): and where p 2 and µ i (i = 1, 2) are the virtuality of the incoming particle and the masses of the particles that circulate in the loop, respectively. Equations (5.33)-(5.35) are derived by assuming p 2 , µ 1 , and µ 2 to be real-valued parameters; both the square roots and the logarithms are meant to be evaluated in the first Riemann sheet. These results can be thus directly applied to an OS-renormalisation procedure by setting p 2 = M 2 , µ 1 = M 1 , and µ 2 = M 2 , where M , M 1 , and M 2 are the relevant OS masses. The next step is to show that eqs. (5.34) and (5.35) still apply to the case of complexified loop masses 18 (in our notation, this implies setting . This is immediately obvious in the case of eq. (5.34), and requires only a slightly more elaborate proof for eq. (5.35), which we refrain from showing here. We understand the assumptions Therefore, in order to arrive at the results that one needs to use in a CM-scheme computation, the only thing that remains to be done is to show how to analytically continue the expressions given above for real and positive p 2 to complex p 2 values -since we need to evaluate the self energies at p 2 = m 2 ≡M 2 − iΓM . This then entails figuring out whether at p 2 =M 2 − iΓM the square roots and logarithms can still be evaluated in the first Riemann sheet, or other Riemann sheets are also necessary.
We proceed as follows. Let f (T (γ)) (5.39) denote any of the elementary functions that appear in eqs. (5.33)-(5.35), with explicit or implicit dependencies on the following quantities: Therefore, the only non-trivial cases are those where f () is either a logarithm or a square root. Furthermore, the latter case can be derived from the former, since for any complex number z and a given real number a one defines: where Log z is the logarithm function whose codomain is the full Riemann surface. Then: with log k z given in eq. (5.23), and R k the (|k| + 1) th Riemann sheet (positive or negative) defined as the complex plane with a branch cut on the negative real axis. The fundamental property of the logarithm function needed here is the direct result of its analytic continuation along a curve C and of the monodromy theorem, namely: Here, C is an oriented curve that starts from the arbitrary complex number z 0 and arrives at the arbitrary complex number z. By ∆ C Arg z we have denoted the variation of the 18 This fact has already been used in the case of the Taylor expansion, where it is crucial (see footnote 17).
A quick, if not fully rigorous, way to argue that this is the case is to observe that, by giving a negative imaginary part to the masses of the particles that circulate in an one-loop bubble diagram, the signs of the imaginary parts of the +i0-regulated Feynman propagators do not change.
multi-valued Arg z function along C. By construction, Arg z = arg k z if Arg z ∈ R k . One also has the properties: Here, C is a curve with the same endpoints as C that can be obtained by a continuous deformation of C without passing through the origin z = 0; C −1 is the same curve as C, with the opposite orientation. Given eq. (5.39), we call the (infinite) set of complex numbers: We further observe that, given eqs. (5.33)-(5.38), T has a zero winding number 19 around the origin. This implies that the final, analytically-continued form of B 0 will feature either first-or second-(both positive and negative) Riemann sheet logarithms. We remark that multi-valued functions can be nested; an explicit example is given in eq. (5.35), where the arguments of the logarithms feature a square root. In such a case, we proceed in an iterative manner. Namely, we first deal with the inner function (the square root), and determine whether the first or the second Riemann sheet is to be used (for each γ ∈ [0,Γ]); then, we apply the procedure to the logarithms. We point out that, in this way, the trajectories relevant to the logarithms include the information on the Riemann sheet employed to evaluate the square roots which, among other things, guarantees that such trajectories are continuous.
The procedure advocated above 20 implies that, starting fromγ = 0, we follow the trajectory by increasingγ till the endpointγ =Γ is reached, counting the number of times in which we cross the branch cut of f (); for each of them, the direction in which the cut 19 The procedure proposed here is valid also for trajectories with non-zero winding numbers; those simply entail the use of logarithms with values in Riemann sheets different from the first or second ones. 20 To the best of our understanding, this is analogous to the integral-level one of ref. [94], at least for a class of trajectories (loosely identified as "monotonic" in the paper). is crossed needs also to be considered. If we denote by n +− (n −+ ) the number of such crossings from the positive to the negative (negative to positive) imaginary part, then:  fig. 2, where the branch cut on the negative real axis is represented by a shaded region, and to be definite we choose f () = log(). According to eq. (5.46), in the case of trajectory A we shall need to turn log() into log + (). For both trajectory C and D log() will have to be turned into log − (), while in the cases of trajectories B and E one will end up using the principal-value logarithm, although for different reasons (in the case of B, the branch cut is crossed twice, in opposite directions; no crossings occur for E). Note the different behaviour of D and E, in spite of the fact that these trajectories have the same endpoints. Figure 2 can also be employed to sketch a couple of shortcuts to following the complete trajectory (which is a numerically involved procedure). In the endpoint method, only the endpoints of T are considered. If T (0) < 0 and T (Γ) < 0, then either log + () (when T (0) > 0 and T (Γ) < 0) or log − () (when T (0) < 0 and T (Γ) > 0) must be employed. Conversely, if either T (0) > 0 and T (Γ) > 0, or T (0) and T (Γ) have the same sign, then log() is used instead. These criteria give the correct results for trajectories A, B, and D. However, they also imply that there are cases, where both the real and the imaginary parts have opposite signs at the two endpoints, which cannot be addressed with this methodan example is that of trajectory C. Finally, in other cases (such as that of trajectory E) the method just gives an incorrect result. As a refinement of the endpoint method the straighttrajectory method can be used. This entails replacing the trajectory with a segment that connects the two endpoints, and then proceeding as before (which is then simply equivalent to finding the intersection of such a segment with the real axis). Thanks to the leftmost identity in eq. (5.44), this method is guaranteed to give the correct result, provided that one is able to understand whether the continuous deformation from the trajectory to the segment does not cross the origin. This is the case for all of the trajectories in fig. 2, except for E -in the latter case, the straight-trajectory method would lead to an incorrect result. It is important to note that both the endpoint and the straight-trajectory methods cannot self-diagnose a failure; some further information on the complete trajectory is necessary. Thus, they cannot be reliably used in the context of an arbitrary model with arbitrary parameter assignments.
The trajectory approach leads straightforwardly from eqs. (5.33) and (5.34) to the following results: otherwise.
Conversely, eq. (5.35) renders an analytical formulation impractical (although possible), and it is best to resort to numerical methods. In what follows, we present sample results obtained by setting masses and widths as specified in table 1 -note thatΓ =M 21 . In fig. 3, we show the trajectories of: which correspond to the argument of the square root that appears in eq. (5.37) and to the argument of one of the logarithms in eq. (5.35), respectively. The quantities γ 0 , γ 1 , and γ + are defined in eqs. (5.37), (5.38), and (5.36), respectively. As far as the left panel of fig. 3 is concerned, we see that for all configurations bar D the square root ends up being computed in the second Riemann sheet. The endpoint method gives the correct results for trajectories A, B, and D, fails for E, and cannot be applied in the case of C as explained before. The straight-trajectory method gives the correct results in all cases except for E, where it fails. 21 The trajectory associated with any valueΓ <M is simply a subset of the trajectory relevant toΓ =M .
This choice of a very large width value is thus simply a practical way to address all situations of interest within a single study. Table 1: Mass and width settings relevant to illustrative studies of the trajectory method. For all of these cases we also setΓ =M . The trajectories we consider are dimensionless, hence there is no need to specify a reference mass scale. Figure 3: Trajectories of eq. (5.49) for the five configurations of table 1. In the left panel, trajectories relevant to configurations A and E are multiplied by the numerical factors reported in the labels in order for them to fit into the layout, and to improve visibility; note that (T (γ = 0)) > 0 for configuration E.
Conversely, the trajectories of the argument of the logarithm in eq. (5.49), depicted in the right panel, show that in all cases log() must be used, except for D that requires the use of log + (). For all trajectories, both the endpoint and the straight-trajectory methods give the correct results.
We now consider the case of configuration E in table 1 (except forΓ, which is allowed to vary in the range 0 <Γ <M ), and present the corresponding results for the finite part of B 0 in fig. 4 (the left panel displays the real part of B 0 /(iπ 2 ), the right panel the imaginary one) as a function ofΓ/M ; we have also set µ 2 =M 2 1 . Several curves appear in fig. 4, each of them obtained by computing B 0 with different approaches and approximations. The black solid line is the prediction of the trajectory method. The endpoint-method result is shown as a red long-dashed curve. By Taylor expanding as explained earlier in this section we obtain the brown short-dashed and blue dotted curves, that correspond to keeping the first two 22 and the first one hundred terms in the expansion, respectively. Furthermore, the green line overlaid with open boxes is the result of the contour integration proposed in ref. [94], while the yellow curve overlaid with full circles is analogous to the former, but with an alternative definition of the contour (which we dub "revised"). Although strictly necessary only forΓ >Γ √ , withΓ √ 0.37M , we have employed the revised contour also forΓ <Γ √ where, as is apparent from fig. 4, its results are identical to those of the contour of ref. [94]. The interested reader can find more details on the revised contour in appendix E.2. Finally, the vertical lines indicate theΓ values at which the trajectories associated with configuration E cross the negative real axis. The leftmost and rightmost ones (atΓ 0.35M andΓ 0.43M , respectively) correspond to the trajectory of (γ + − 1)/γ + , while the central one (atΓ =Γ √ ) corresponds to the trajectory of γ 2 0 − 4γ 1 . We see from fig. 4 that, for very smallΓ values, all methods give identical results. By increasingΓ, the Taylor expansion limited to two terms departs quickly from the other predictions, which are essentially on top of each other up toΓ Γ√ . AtΓ =Γ √ the endpoint method is unable to figure out correctly the appropriate Riemann sheet for the calculation of the square root (see the left panel of fig. 3), and thus differs significantly from the full-trajectory result for all valuesΓ ≥Γ √ . This is also the region where the Taylor expansion breaks down, irrespective of how many terms are kept, and where the contour integration of ref. [94] leads to a cusp-like behaviour inΓ for both the real and imaginary parts of B 0 . Conversely, by using the revised version of the contour, the resulting predictions are in perfect agreement with those of the full-trajectory approach described

earlier.
We conclude this section with a few remarks. We point out that in the SM trajectories are necessarily "short" (sinceΓ M ), and start from points located close to the real axis (sinceΓ i M i ). This is the reason why, apart from exceptional cases (such as that of eq. (5.34)), the Taylor expansion approach can be shown to be a reasonable one; likewise, the endpoint method is expected to yield the correct analytical continuation (also thanks to the fact that the SM mass spectrum is not close to being degenerate) 23 . Thus, for the CM-scheme SM UFO models shipped with MG5 aMC, we have limited ourselves to implementing the latter approach. We stress, however, that it is crucial to maintain the full flexibility of the trajectory method in view of simulations within models with arbitrary assignments of masses and width, in keeping with the philosophy that underpins MG5 aMC. In fact, it is relatively easy to automate the study of eq. (5.45), e.g. starting with a coarse discretisation of the parameter rangeγ ∈ [0,Γ], which can be refined if necessary (for example when nearing the branch cut). We present a couple of possible approaches in appendix E.2, one of which we have used to obtain the results presented in this section.

5.4
About the phase of α in the complex-mass G µ scheme As was already discussed, derived couplings can potentially become complex when they depend on couplings and/or masses that acquire an imaginary part as a result of the renormalisation conditions in the CM scheme. In particular, this is the case of the EW coupling α in the CM G µ scheme, which reads: (5.50) Here, we have denoted by G (Gµ) µ the renormalized Fermi constant, whose superscript Gµ reminds one that we are working in the CM G µ scheme, where G (Gµ) µ is real-valued. The complex values of the Z and W boson masses, defined according to eq. (5.7), imply that α has a non-zero phase.
In the context of NLO computations this fact is problematic, since it may lead to uncancelled IR singularities in the context of subtraction procedures. In order to illustrate this point with a simple example, we consider selected virtual and real-emission contributions to the NLO 3 corrections (i.e., of O(α 2 S α 2 )) the process qq → q q g -such contributions are shown in fig. 5. At the level of matrix elements, one will need to consider (see sect. 3.1) the real parts of the diagrams of fig. 5. We sketchily write these as follows: Thus, for the process we are considering we can write this integrated kernel symbolically as follows: with 0 ≤ n ≤ 2 a number that depends on the specific Born-level amplitudes B 1 and B 2 . It is therefore apparent that there is a mismatch due to α between eq. (5.54), and eqs. (5.51) and (5.52). Even by bearing in mind that eq. (5.54) has been worked out by assuming that α is real, the different ways in which α enters eqs. (5.51) and (5.52) imply that the generalisation of eq. (5.54) for complex-valued α is involved, if at all possible.
We note that this problem is not specific to the computation of some subleading term such as the NLO 3 of the present example. Suppose that one studies the same process, but in the pure-QED case, by replacing gluons with photons. Equations (5.51)-(5.54) become then: While it appears that eq. (5.58) can be modified so as to have the same singular behaviour as both eq. (5.55) and (5.56), by doing so it would not any longer be equal to the integral of its local counterpart, which is employed in the subtraction of eq. (5.57).
In order to address the issue we have just discussed, what is commonly done in the CM G µ scheme is to use, in place of of eq. (5.50), a real-valued α Re : (5.59) We point out that the choice of an absolute value in eq. (5.59) is arbitrary, and e.g. employing instead the real value of the r.h.s. of eq. (5.50) is possible, as it would serve the same purpose. However, any such solution is not really appealing, in that it might spoil a gauge relation by higher-order terms. In particular, derived parameters other than α can assume different values, depending on whether they are expressed directly in terms of all of the input parameters of the CM G µ scheme, or by using those inputs with the exception of G (Gµ) µ , replaced by α Re . Although counter-intuitive, in view of the problem posed by IR singularities it is the latter approach that one must adopt in order to prevent gauge relations from being broken.
Another way to remedy this situation, equivalent at the NLO to the one we have just discussed, is that of turning G This effectively amounts to working in another scheme of the CM class, which we dub G µ scheme and where one uses G In conclusion, the real-valued α of eq. (5.59) emerges automatically in the G µ scheme, where one reshuffles higher-order terms in order to trade the phase of the coupling α for a suitable phase of G (Gµ) µ . As a result, one retains both order-by-order IR pole cancellations and gauge-invariant expressions of derived SM couplings. The numerical predictions presented in sect. 6 have been obtained in the G µ scheme, i.e. are identical to those of the G µ scheme with α = α Re used everywhere.

About enforcing zero widths for final-state unstable particles
In order to guarantee the unitarity of the S matrix, one must consider stable particles in the final state (see e.g. ref. [95]). Therefore, the widths of potentially unstable final-state particles must be set equal to zero; in the CM scheme, this implies that the corresponding fields are then renormalised with the same conditions as in the OS scheme. On the other hand, it has been already argued in sect. 5.1 that, for processes that feature unstableparticle resonances, it is desirable to resum the relevant 1PI insertions, thereby obtaining natural regulators in the propagators, which prevent the numerical integration from diverging in certain regions of the phase space. Because of this, the situation arises where one must assign non-zero widths to the unstable particles that appear as resonances in a given process, and simultaneously set equal to zero the widths of final-state unstable particles. We shall denote the sets composed of such particles by P res and P fs , respectively. Whenever these two sets P res and P fs overlap, one cannot perform the computation consistently, and the decays of the final-state unstable particles must be included 24 . This is for example the case of the top quark in the process 25 pp → tW −b , or of the Z boson in the process pp → Ze + e − .
In this section we discuss the case in which the sets P res and P fs do not overlap, and consider the question of whether a consistent NLO computation within the CM scheme is possible. In order to be definite, it may be useful to refer to an explicit example, which we shall take to be the NLO 3 corrections to the process pp → ttj (this corresponds to O(α 2 S α 2 )) due to real-emission contributions such as ud → tt(W + * →)ud or uū → tt(Z * →)dd . The latter require non-zero W -and Z-boson widths, while the width of the top quark must be set equal to zero -thus, P fs = {t,t} and P res = {W + , Z}. Given that the top quark or antiquark does not appear as a resonance in this process, one could anticipate that this width assignment is not problematic.
More generally, our conclusions are the following. In the context of a CM-scheme calculation, it is always possible to impose OS-type renormalisation conditions on potentiallyunstable particles, provided that the latter only appear in the final state and not as intermediate resonances that may go on shell.
The previous seemingly trivial statement crucially depends, among other things, on a correct interpretation of the operators˜ and †, as introduced and used in ref. [89]. We 24 This statement relies on the implicit assumption that the pole-mass regions of the resonances are integrated over. 25 It is crucial that EW corrections are driven by the same interactions through which the top decays. As is well known, pure-QCD NLO corrections to tW associated production can be computed also with stable tops, although admittedly in a non-trivial manner.
shall comment on these aspects in what follows. For what concerns˜ , in particular, we point out that in a straightforward application of the CM scheme this operator is simply irrelevant, while in the OS scheme, because of the absence of any complex-valued couplings, the operator is sufficient -it is only in a "mixed" CM-OS setup that˜ becomes necessary. Let us consider the particular contribution 26Γt (W + ,b) (p 2 ) to the top quark renormalised irreducible two-point functionΓ t (p 2 ) due to a (W + , b) bubble graph. We choose it to exemplify a general situation with analytical expressions which are as simple as possiblethe conclusions apply universally. Such a contribution reads as follows: where s W is the sine of the Weinberg angle, M t the (OS) top-quark mass (see sect. 5.2), and w ± = (1 ± γ 5 )/2. The application of OS-renormalisation conditions for a massive fermion, sketched in sect. 5.2 and reported explicitly in ref. [89], leads to the following expression for the top quark wavefunction renormalisation counterterm: According to ref. [89], the notation˜ understands that the real part must be taken of only the loop functions, and not of the factorised couplings. The expressions for the unrenormalised self-energy functions Σ t,S/L/R (W + ,b) (where we assume a diagonal CKM matrix) read as follows: where the coefficients κ(p 2 ) are finite and real. As it has been already noted in sect. 5.2, the optical theorem relates the onshell coefficients κ L/R ,fin (p 2 = M 2 t ) to the perturbative LO width of the top quark, independently of any input value Γ t 27 . This absorptive part of the loop function is precisely what is intended to be removed by the˜ operator in eq. (5.62). 26 At variance with what is done elsewhere and in order to avoid the proliferation of subscripts, in this section we denote renormalised quantities by means of a hat symbol. 27 Note that the finite widths of the particles running in the loop can also contribute an imaginary part to the self-energies of unstable particles. Such contributions are always beyond NLO accuracy and can be finite, as in δ (W + ,b) Zt, but also UV divergent as in δ (W + ,γ) Z W + . We have chosen to define˜ so that it removes the complete imaginary finite part (irrespective of its origin) of the loop functions defined with couplings factored out, but not its singular part; thus retaining the exact IR poles cancellation at fixed order.
By substituting eqs. (5.63)-(5.65) into eq. (5.61), we find the following expressions for the top quark wavefunction counterterms 28 : In the context of our case-study computation of the NLO 3 correction to pp → ttj, the OS renormalisation conditions apply to the top quark only and the factors αs −2 W and α m −2 W s −2 W take the complex values assigned by the CM renormalisation conditions. It is therefore crucial that the real part˜ does not apply to the factorised couplings, so as to guarantee the cancellation of UV divergent terms. However, the finite imaginary absorptive part iκ L/R ,fin (p 2 = M 2 t ) must not be included in the top mass OS counterterm δ (W + ,b) M t , since the latter is the fixed-order counterpart of the width regulator that appears in G (t) R (p 2 ) and that is set to zero there just in force of the OS conditions. In fact, this selective behaviour of the operator˜ is mandatory in order to maintain gauge relations (see eqs. (3.30) and (3.31) in ref. [89]), guaranteeing the fermion-flavour universality of the coupling constant α. We finally point out that the same logic demands that the operator˜ act as the identity on the CKM matrix elements.
We point out that the observations made above in principle apply also to the wavefunction renormalisation constants δZ L/R q for a massless fermions q. However, their contributions to the one-loop matrix-element end up being added incoherently: with A L/R is the Born amplitude for the production of a left-or right-handed massless quark q. Thus, the imaginary part of δZ L/R q is irrelevant in practice for one-loop computations. This is not the case for massive fermions, because they can undergo a chirality flip, whence the imaginary part of δZ L/R t does contribute to the one-loop matrix elements of processes featuring final-state top quarks. It is therefore important to perform any complex-conjugate operation in a correct manner.
In view of what has just been discussed, we conclude this section by remarking that we find the † operator as introduced in ref. [89] to be potentially misleading. To give one example, eq. (3.21) therein reads:  [89]). See also eq. (5.70) and the discussion related to it. and may appear inconsistent since in general δZ ii . For this reason, we stress here that we understand any quantity δZ L/R † ij to be an independent renormalisation constant associated with the antiquarksī andj, and not as a derived parameter obtained by (possibly) transposed complex conjugate of δZ To underscore this fact, we prefer to use the notation: [89]) . (5.70) Thus, δZt(our work) ≡ δZ † 33 (ref. [89]), whence δZ . This makes it more explicit that the quantities δZ L/R ij are independently derived from the renormalisation conditions of the anti-quark fields, identical to the quark ones but with the following substitution of Dirac spinors: u(p) →v(p) andū(p) → v(p).

Results
In this section we present illustrative results relevant to several hadroproduction processes. We start by discussing in sect. 6.1 the setup for the MG5 aMC runs -the model used, the input parameters, and the definition of the basic observables. In sect. 6.2 we give our predictions for total rates and some selected differential distributions, by considering processes that display a significant diversity in their final states 29 . For such processes, we compute that is, the leading LO and second-leading NLO (i.e. NLO EW) contributions. We do not include the HBR cross sections (see sect. 2) in the latter. Finally, in sect. 6.3 we study the production of a tt pair, possibly in association with either a heavy boson (Z, W + , and H) or a light jet, again at both the fully-inclusive and the differential level. In the case of tt(+B, j) production, we shall report the complete LO and NLO results, namely: with Σ LO 4 and Σ NLO 5 being relevant only to ttj production. As for the case of eq. (6.1), HBR contributions have not been included. We stress that the goal of this section is to document the achievement of the full automation in MG5 aMC of mixed-coupling calculations, and the typical usage of the code. Therefore, we did not choose our settings on a process-by-process basis, as one would do if the primary concern were a phenomenology study, but rather preferred to impose the same conditions on all processes, so as to facilitate direct comparisons among them. Note that while the phenomenological impact of HBR cross sections should always be assessed (see e.g. refs [22,25,39,[96][97][98][99][100], and refs. [101][102][103] for shower-type approaches) in the context 29 As was mentioned at the end of sect. 2, features related to FFs are not publicly released yet, and therefore we do not consider processes with tagged photons. Photons are treated on equal footing with QCD partons (see appendix D), and either enter a jet or help define a dressed charged lepton. of a realistic analysis, this paper emphasises the technicalities of genuine NLO automated computations. Therefore HBR results, that are obtained by means of IR-finite tree-level calculations which are straightforward to carry out with LO-type runs in MG5 aMC, are unimportant here.
Likewise, as far as systematics are concerned one of the key aspects for the present paper is that of a good control over the numerical accuracy of the final results. This is nontrivial, owing to the expected strong numerical hierarchy among the different contributions to the LO and NLO cross sections in a mixed-coupling expansion. We shall therefore limit ourselves to reporting the errors associated with the Monte Carlo integrations over the phase spaces. Thus we shall neglect the uncertainties due to renormalisation-scale, factorisation-scale, and PDFs variations. However, we point out that such uncertainties can be obtained with MG5 aMC at no extra computational costs, thanks to the reweighting procedure introduced in ref. [76] which has been extended to the mixed-coupling scenario -see appendix B.

Setup of the calculation
The UFO [104] model we use, loop qcd qed sm Gmu, is included in the standard MG5 aMC package. It contains the UV and R 2 counterterms relevant to NLO QCD and EW corrections, the latter in the G µ scheme (see sect. 5.4). The model features five massless quark flavours, sets the CKM matrix equal to the identity, and is compatible with the usage of both the OS and the CM schemes for all massive particles (W ± , Z, H, and top quark). Prior to process generation, we therefore execute the following commands: MG5 aMC> set complex mass scheme true MG5 aMC> import model loop qcd qed sm Gmu MG5 aMC> define p = g d d~u u~s s~c c~b b~a MG5 aMC> define j = g d d~u u~s s~c c~b b~a The latter two instructions tell the code that the photon (denoted by the symbol a) must be considered part of the proton and of the jets, in keeping with the democratic approach -more details are given in appendix D.
We consider pp collisions at a center of mass energy of 13 TeV (LHC Run II). The values of masses and widths of the two heavy vector bosons are: In the case of the Higgs, given the smallness of its SM width (Γ SM H = 4.07 MeV) the choice Γ H = 0 is always an excellent approximation from a physics viewpoint, that allows one to compute Higgs production in various channels without bothering with its decays 30 . Therefore, we have always assumed Γ H = 0, except for the pp → e + ν e jj and pp → e + e − jj processes, where Γ H = Γ SM H has been employed 31 . We remind the reader that the self-consistency of the calculation (see also sect. 5.5) demands that, in processes where any massive particle (W , Z, top quark, and Higgs) is treated as stable, i.e. is left undecayed, its width has to be set equal to zero, and it must not appear as an intermediate resonance that may go on shell. In the case of the vector bosons, this implies setting Γ BW Z = 0 and Γ BW W = 0 before employing eq. (6.6) to obtain the CM-scheme inputs.
The value of the Fermi constant as extracted from the muon decays, which is an input in the theory model adopted here, is set equal to: The EW coupling α is derived from eq. (6.9) and from the W and Z CM-scheme masses as explained in sect. 5.4 -see in particular eq. (5.59) and the discussion given there. The PDFs are the central ones of the LUXqed plus PDF4LHC15 nnlo 100 set [109,110], extracted from LHAPDF6 [111] with number 82000; these are associated with α S (m Z ) = 0.118 , (6.10) with a three-loop running (which is performed by LHAPDF6). The renormalisation and factorisation scales have been set as follows: where the sum is extended to all final-state particles at the parton level (i.e. prior to possibly combining them into jets or dressed leptons). In order to define final-state objects in an IR-safe way, we apply the following set of minimal selection cuts. 30 This fact is often abused in BSM models with large-width scalars, where a "Higgs cross section" is not necessarily a meaningful concept. 31 These processes receive contributions from a subset of one-loop diagrams that feature an s-channel Higgs boson propagator. At the orders that we are considering, eq. (6.1), such diagrams are interfered with LO diagrams that do not have Higgs propagators, thus leading to integrable singularities even when ΓH = 0. Therefore, we have set ΓH = 0 in these cases solely to improve the behaviour of the numerical integration.
• All photons that are within ∆R f ± γ ≤ 0.1 of a charged fermion (either a quark or a lepton), are recombined with that fermion (by definition, ∆R ij = (∆φ ij ) 2 + (∆η ij ) 2 ). If there is more than one charged fermion candidate, the photon is recombined with the one that yields the smallest ∆R f ± γ . The four-momentum assigned to each dressed fermion is the sum of the momenta of the original fermion and the photon. Hence, after recombination, these dressed fermions possibly feature a non-zero mass.
• Jets are reconstructed from all massless QCD partons (i.e. gluons and quarks, the latter might be a recombined quark-photon pair) and from photons that were not recombined with charged fermions. They are defined with the anti-k T algorithm [112] (as implemented in FastJet [113]) with R = 0.4, and subject to the conditions p T (j) > 30 GeV, |η(j)| < 4.5.
• All jets so obtained are kept in the analysis (in other words, photon-jets [32] are not tagged as such). Leptons and jets are not required to be separated in ∆R.
The results that follow are obtained by integrating simultaneously all of the contributions we have considered, in an iterative manner till a target accuracy is reached or surpassed. By introducing the shorthand keywords LO and NLO as follows: for the processes of sect. 6.2 (see eq. (6.1)), and: for the processes of sect. 6.3 (see eqs. (6.2) and (6.3), with n LO = 3 and n NLO = 4 for tt(+B) production (n LO = 4 and n NLO = 5 for ttj production), the target accuracy is defined to be a relative MC error equal to 0.02% on the total NLO cross section, within cuts (NLO being that of eqs. (6.12) and (6.13)). We remind the reader that MG5 aMC determines automatically the number of phase-space points sampled in each MC iteration to achieve a given accuracy -more details on this item can be found in sect. 2.4.3 of ref. [17]. The overall runtime to compute all of the results presented in this section is a couple of weeks on O(200) CPUs. As was mentioned in sect. 2, MadLoop can choose dynamically which integral-reduction module to employ; this is done in an order that is pre-defined by the user. In the current MadLoop in this order, if necessary). The Ninja and CutTools integral-reduction modules obtain the scalar integrals from OneLoop [114]. For the processes considered in this paper, we have found that, with the accuracy as specified above, an overall (i.e. relevant to all of the processes combined) negligible amount of O(100) phase-space points have required quadruple-precision calculations, all of which were then deemed to be numerically stable.

NLO EW corrections
In this section we present the leading LO and second-leading NLO (i.e. NLO EW) results for a variety of processes, whose complete list can be found in the first column of table 2. The second column of that table reports instead the MG5 aMC commands used to generate those processes. These adhere to the general syntax reported in sect. 2; note in particular the keywords 32 that determine which coupling-constant combinations are considered in the calculations, according to eqs. (2.8) and (2.9).
We start by looking at fully-inclusive rates, obtained with the conditions and acceptance cuts given in sect. 6.1. The third and fourth columns in table 2 report the LO and NLO results, defined according to eq. (6.12). The fifth column displays instead the fractional correction (given in percent) due to NLO EW effects, i.e.: As was anticipated in sect. 6.1, all of the uncertainties reported in the three rightmost columns in table 2 are MC integration errors; as one can see, in absolute value they are almost always well below the per-mille level 33 . Table 2 confirms the well-known fact that NLO EW effects to fairly inclusive observables are mostly negative, and rather small in absolute value (a few percent). Several (but not all) of the triple-boson production processes constitute an exception to the latter rule, the largest correction being that associated with the HHW + final state, where δ EW = −12.8%. For such processes, the four largest δ EW 's are all negative and of O(−10%), while the largest positive correction is δ EW = 6.2% in W + W − W + production.
It should not come as a complete surprise that the corrections are typically larger for processes with large transferred momenta and several final-state bosons, such as in the case of triple heavy boson production. Indeed, two of the most significant contributions to EW corrections, the EW Sudakov logarithms and the initial-state photon-induced terms, are relatively enhanced in this kinematics regime and for these final states. In the differential distributions plotted and discussed below, we shall touch upon these effects in more detail. 32 The keyword [QED] is conventional, and it implies that both electromagnetic and weak effects (i.e. the complete O(α) corrections) are taken into account, since both are included in the loop qcd qed sm Gmu model. Restrictions to the QED-only or weak-only cases can be achieved by adopting a simpler theory model (for those processes for which these restrictions are meaningful). 33 The largest fractional error (still a mere 1.1·10 −3 on the NLO cross section) affects HHW + production.
We have checked that this is dominated by the opening at the NLO of a new t-channel configuration where an initial-state photon couples directly to the W + . This channel is not mapped ideally by our phase-space parametrisation. 1.0613 ± 0.0001 · 10 2

Process
1.0539 ± 0.0001 · 10 2 −0.70 ± 0.02 Table 2: Processes considered in sect. 6.2. The second column reports the MG5 aMC syntax used to generate them. The third and fourth columns display the fully-inclusive results for the quantities defined in eq. (6.12). The fifth column shows the results for the fractional correction defined in eq. (6.14). All uncertainties are due to MC integration errors.  We now turn to presenting a few sample differential distributions. We have collected in a single figure the results relevant to processes characterised by similar final states. All of the figures have the same layout, constituted by a main panel and an inset. In the main panel, LO and NLO results (defined according to eq. (6.12)) are shown as dashed and solid histograms, respectively. The histograms are normalised so that the content of each bin is the total cross section in that bin. The inset displays the predictions for 1 + δ EW , which is therefore equal to the ratio of the NLO over LO results presented in the main panel.

Vector boson (plus jets) production
In left panel of fig. 6 the scalar sum of transverse momenta, H T , is shown for W + * production, possibly in association with jets. The processes we have considered are: pp −→ e + ν e , pp −→ e + ν e j , pp −→ e + ν e jj . (6.15) In order to be definite, we have restricted ourselves to presenting results only for the cases of positively charged leptons. Obviously, the generation of W − * +jets processes is fully analogous to that carried out here. NLO EW corrections to W +jets production have already appeared in the literature in various approximations [21,23,24,[115][116][117][118][119]; we shall not compare those predictions with ours in this paper 34 . The contributions to H T by the final-state particles that appear in eq. (6.15) are due to the positron (possibly dressed with a photon), to the missing transverse momentum (set equal, by using MC truth, to the transverse momentum of the neutrino), and to any reconstructed jets. For the simplest process in eq. (6.15), pp → e + ν e , at the LO and in the narrow width approximation it would be kinematically impossible to have H T > M W . The inclusion of finite-width effects for the W boson opens up this region of phase space, which 34 Comparisons of results obtained with then-private MG5 aMC versions with those computed by means of other tools have been previously reported in ref. [120] (for ttH production) and ref. [1] (for the processes of eq. (6.17)). is thus populated exclusively by configurations where the W boson is off-shell. Therefore, the cross section for H T > M W remains very suppressed at the LO, as is evident from the green dashed curve in the main frame. When NLO corrections are included, real emissions contribute to this region also when the W is on-shell, leading to very large and positive corrections -the K factor, shown in the lower inset, is equal to about 2 for H T 500 GeV. It is again because of the outsize impact of the real corrections that the K factor for the process pp → e + ν e j grows rapidly and assumes large values in the region H T 1000 GeV. There, one is dominated by configurations where there are two hard, back-to-back jets, and the W boson has a small p T relative to either jet. These kinds of enhancement due to kinematic features that become available only at a certain perturbative order beyond the leading one are of course not peculiar to EW corrections, and indeed they have often been studied in QCD as well (see e.g. refs. [118,[121][122][123][124]). Eventually, by increasing the jet multiplicity, the complexity of the final state is sufficiently large already at the Born level, and such effects are pushed towards highly-constrained corners of the phase space. For example, in the case we are considering here, the EW corrections to pp → e + ν e jj behave as one naively expects them to do, becoming negative and growing with H T (as is typical in a regime dominated by EW Sudakovs).
On the other hand, it is in the comparison between processes that differ by final-state jet multiplicities where the QCD and EW cases deviate from each other. This is because in pure QCD the real corrections to pp → e + ν e are identical to the Born of pp → e + ν e j. Hence, in kinematic regions dominated by real emissions (e.g. at large H T in the present case) one expects the NLO result of the former process to be rather close to the LO result of the latter one. The same argument holds when comparing pp → e + ν e j with pp → e + ν e jj. But when one computes EW corrections instead of QCD ones, these considerations are no longer valid: the Born-level jet(s) of pp → e + ν e j(j) is (are) still predominantly emerging from a QCD branching, while the one(s) that contributes at the NLO to pp → e + ν e (j) production has (have) a rate governed by EW effects, and therefore much smaller.
NLO subleading contributions to pp → e + ν e j feature doubly-resonant diagrams relevant to pp → W + * Z * and pp → W + * W − * production, with the Z * and W − * "decaying" hadronically. Among such subleading terms, the one with the largest α S power (NLO 2 ) is due to the interference of these di-boson diagrams with QCD non-resonant ones. This, and the related fact that several of these interference contributions are identically equal to zero owing to their colour structures, implies that their presence in our calculations poses neither numerical nor phenomenological problems. When they are squared, di-boson doubly-resonant diagrams enter only the NLO 3 term (that is of O(α/α S ) relative to the NLO 2 one), which we do not consider here 35 .
We now turn to discussing Z * (+jets) production; as in the case of W + * (+jets) production, no comparisons will be given in this paper with previous NLO EW results [18,24,[125][126][127][128][129]. The plot on the right panel of fig. 6 presents the invariant mass of the charged-lepton 35 In an analogous manner, pp → W + * Z * j and pp → W + * W − * j diagrams enter the LO subleading contributions to pp → e + νejj production, starting at LO2 where they interfere with QCD non-resonant diagrams, and at the squared level in the LO3 term. Because of eq. (6.12), these terms do not concern us. pair in the processes: The invariant mass is computed for dressed leptons (according to the definition given in sect. 6.1). All of the three processes of eq. (6.16) clearly display the expected Z-boson mass peak. The flatness of the histograms to the left of the peak is to a certain extent an artifact of our plotting choice (a logarithmic x axis) -we stress that the contributions due to an s-channel photon exchange (as well as those without any s-channel "decay" to the e + e − pair), which rapidly grow with decreasing invariant masses, are included in our results.
From the inset of the right panel of fig. 6, we see that the NLO EW corrections are fractionally the largest just below the Z-boson peak. Indeed non-soft wide-angle (i.e. out of the dressed-lepton cone) real-emission photon radiation from the final state leptons decreases the reconstructed invariant mass w.r.t. that determined by the bare leptons, thus essentially shifting contributions from the Z-boson peak to the bins to its left 36 . On the other hand, the narrowness of the ∆R cone that defines the dressed leptons renders it much harder to shift the bare-lepton mass to the right of the peak, which would occur by dressing either leptons with a photon emitted by any non-leptonic hard particle.
Furthermore, in the region M (e + e − ) < 90 GeV the EW corrections are larger (and significantly so in the two bins to the left of the Z peak) for pp → e + e − j production than for pp → e + e − jj and pp → e + e − production. There is a hierarchy among these NLO effects in the latter two processes as well, which is visible given our MC uncertainties, but marginal in absolute value. We point out that the relative behaviours of the corrections affecting different processes are in part controlled by the chosen acceptance cuts. More specifically, when one considers pp → e + e − j at the LO, the requirement that the jet have p T (j) > 30 GeV implies p T (e + e − ) > 30 GeV. This renders it very likely that the leptons cuts p T (l ± ) > 10 GeV are trivially satisfied. In turn, this implies that when either lepton emits a real photon (a configuration that, as we have argued before, is responsible for giving the dominant contribution to NLO corrections to the left of the Z peak), lepton cuts will still be easily passed. On the other hand, in the case of pp → e + e − , at the LO p T (e + e − ) = 0. Therefore, lepton cuts imposed on real-emission configurations where the photon does not belong to a dressed lepton will be more stringent than those imposed on the associated LO configurations, and in this way NLO corrections will be relatively more suppressed than for pp → e + e − j production. Finally, the case of pp → e + e − jj in an intermediate one between the other two processes. Note that it is still comparatively easy to have a LO p T (e + e − ) 0 configuration with two hard back-to-back jets, in which case the mechanism advocated for pp → e + e − production plays a role here as well.
Four-lepton production, and VBF Higgs and associated production In the left panel of fig. 7 we present the transverse momentum of the hardest same-flavour lepton pair for the processes:   We have chosen the two processes in eq. (6.17) in order to be definite, as representatives of the class of reactions with four final-state leptons; both have been studied before [26,27,30,35,130,131]. In fact, without any additional complications, MG5 aMC is able to deal with any process that belongs to this class, regardless of the particular flavour and charge combinations.
In detail, the definitions of the p T (ll) (relevant to pp → e + e − µ + µ − ) and p T (lν) (relevant to pp → e + ν e µ −ν µ ) observables are the following. For the former, one uses dressed leptons; the e + e − and µ + µ − pairs transverse momenta are then computed, and the largest of the two is set equal to p T (ll). In the latter case, charged leptons are again dressed first; then, the transverse momenta of the e + ν e and µ −ν µ pairs are computed (by using the MC truth information to find the neutrinos), and the largest of the two is set equal to p T (lν). The NLO EW corrections behave rather differently for the two processes. While for the four charged lepton process they display the typical Sudakov behaviour at high p T , for the other process the corrections are positive and growing for p T 40 GeV, starting to decrease only towards p T 400 GeV. We point out that the two processes have significant differences in their underlying mechanisms. Firstly, although both 2l2ν and 4l production are dominated by di-boson resonant contributions (namely, di-W and di-Z, respectively), it is only the former case that features diagrams with t-channel spin-one exchanges (thus enhanced at large momentum transfers). These appear in γγ-initiated processes, owing to the direct γW + W − coupling. Secondly, partonic processes such as γq → W + * W − * q that give rise to 2l2ν final states may be enhanced at large lepton-pair p T 's owing to quasi-collinear q * → W * q splittings (see e.g. ref. [121]). While a similar mechanism also occurs in 4l production, in that case its effects are balanced by a stronger suppression than in the case of 2l2ν production 37 . Finally, at the NLO 2l2ν production features a 37 The overall impact of quasi-collinear enhancements on observable cross sections ultimately depends on the interplay between their kinematics characteristics, the partonic matrix elements, and PDF effects -see real-emission contribution due to an underlying tW doubly-resonant mechanism, which might induce very large corrections driven by the top-quark on-shell region. In practice, this does not happen -the resonant channels are always associated with an initial-state γb (or γb) pair, and are thus significantly suppressed by parton luminosities, so that the "convergence" of the perturbative series is not spoiled. We note that the same remark applies to the computation of the NLO QCD corrections to W + W − production where the effects, although larger [134] than those relevant to the NLO EW case considered here, are still under control 38 .
The right panel of fig. 7 displays the transverse momentum of the Higgs boson that emerges from the processes: The EW corrections to these processes have appeared in the literature in various approximations [135][136][137][138][139]. Obviously, the study of He −ν e production, which we have chosen to ignore, is identical to that of the rightmost process in eq. (6.18).
The cross sections for pp → He + e − and pp → He + ν e receive (dominant) contributions from the underlying resonant production of Higgs in association with a Z * and a W + * boson, respectively. These mechanisms do also contribute to the leftmost process in eq. (6.18), in which case the vector bosons "decay" hadronically rather than leptonically. By weighting with the different hadronic and leptonic vector-boson branching ratios, with a back-of-the-envelope calculation one can compare the total rates for Hjj production with the sum of those relevant to He + e − plus He + ν e production, and find that the former is about four times larger than the latter sum. This is a first naive indication that pp → Hjj production is dominated by a VBF mechanism, rather than by the associated production ones; clearly, a much better evidence for this comes from studying the topology of the outgoing jets.
The NLO EW corrections behave rather similarly for the three processes in eq. (6.18), and are about O(−5%) at small transverse momentum, p T (H) 200 GeV. They increase in absolute value with p T (H) and remain negative, as is typical for these types of corrections in regions of phase space dominated by EW Sudakov logarithms.

Triple-boson production
In the left panel of fig. 8 we present the transverse momentum of the hardest vector boson in triple vector boson production. This class of processes has been considered previously e.g. refs. [132,133] for discussions on this point. 38 The presence of resonant contributions does not constitute a problem per se; other quantities entering the cross sections, particularly the couplings and the parton luminosities, play a crucial role as well. Probably the best known among the problematic cases is that of the NLO QCD corrections to tW production, which receives an enormous gg → tt-channel enhancement. We point out that it is this underlying tW − tt mechanism that renders the NNLO QCD calculation of the W + W − cross section prone to large corrections. Although solutions can be devised (see footnote 39 for some examples) that eliminate or reduce the impact of resonant contributions, they are always arbitrary to a certain extent, and it is difficult to quantify the uncertainties associated with them. Hence, they are best avoided if not deemed to be strictly necessary.
in refs. [140][141][142][143][144]; we focus on: The computations relevant to final states which are obtained by those in eq. (6.19) by charge conjugation, namely W − W + W − and ZZW − , do not pose any additional complications, and are thus ignored here. Conversely, pp → ZW − W + is more problematic, since it receives contributions from real-emission diagrams with an s-channel top quark (i.e. from an underlying t * W − Z ort * W + Z production mechanism). Thus, while technically this process is doable in our setup by setting the top width equal to its physical value in order to prevent the matrix elements from diverging on the top resonance (see sect. 5.5), potentially it still poses the problems common to all processes which, at the NLO, "interfere" with a top-induced "background" (such as instabilities in the numerical integration caused by extremely large K factors). We have already discussed an example (W + * W − * production, eq. (6.17)) where such an interference in practice does not lead to any issues at the perturbative orders we are interested in. However, the case of ZW − W + production is much more involved, and therefore we prefer to postpone its study to when MG5 aMC will feature an automated treatment of the subtraction or removal of resonant contributions, with procedures analogous to those already considered in the literature in different contexts 39 . Another, simpler, solution is that of performing the computation in a scheme with four flavours. This will not be done here, but it is feasible with the present version of MG5 aMC (we note that a 4FS restriction of the OS model is available, while its CM counterpart has still to be constructed) 40 . From the inset in left panel of fig. 8, we see that ZZZ production exhibits the typical behaviour of NLO EW corrections, which are small at small transverse momentum, and grow in absolute value with p T . The other two processes in eq. (6.19) display a more intricate behaviour, owing to a combination of effects: the virtual Sudakov corrections, which decrease the rates; and the positive enhancement of the cross section, due to the presence of photon in the initial state. The latter are related to the direct coupling of W 's with photons, and entail the opening at the NLO of new t-channel contributions which are absent in pp → ZZZ production.
In the right panel of fig. 8 we consider the transverse momentum of the hardest Higgs boson in the following heavy-boson production processes: The NLO EW corrections for these processes are computed here for the first time. For reasons fully analogous to those that apply to triple vector-boson production (reported in the discussion below eq. (6.19)), the process pp → HW + W − , for which top-quark resonant contributions appear at the NLO, is not studied in what follows. Triple-H production 39 The procedures that are being implemented in MG5 aMC are fully local in the phase-space of finalstate particles, such as those of refs. [145][146][147][148][149][150][151][152][153]. Global [134,[154][155][156] or semi-local [157][158][159][160] approaches are not suited to automated observable-independent short-distance computations. 40 Another possibility in the context of a five-flavour computation is that of adding a dedicated integration channel for each of the new resonant contributions that open at the NLO level. is loop-induced, and therefore is also ignored. As in all of the other cases treated so far, processes obtained by means of charge conjugation from those of eq. (6.20) can be generated without problems by MG5 aMC, but have not been considered here. For inclusive rates (see table 2) the NLO EW corrections are −9% for pp → HZZ, −11% for pp → HHZ, and −13% for pp → HHW + , while for pp → HW + Z they are a positive 1.6%. At the differential level, all of the four processes display the typical behaviour of EW corrections (i.e. negative and growing in absolute value with p T ) at large transverse momenta; however, the p T values for which these effects become dominant do depend on the specific process. In particular, as is the case for the inclusive rates, it is HW + Z production that stands apart, since up to relatively large transverse momenta (p T 200 GeV) the negative contributions due to the EW Sudakovs (which are present in the other three processes as well) are compensated by positive contributions. Among these, the dominant one is driven by a quasi-collinear enhancement stemming from γq → HW + q * (→ Zq), a mechanism fully analogous to that already advocated for the second process in eq. (6.17), and that cannot be present in the other three processes in eq. (6.20). Finally, we notice that (smaller) differences between the triple-boson processes of eq. (6.20) can be induced by virtual corrections, owing to the different ways in which the bosons enter the one-loop diagrams (chiefly, by being directly attached to the heavy-quark loop, or by resulting from the branching of a parent particle that is directly attached to the loop).

Associated top-quark, and jet production
In the left panel of fig. 9 we consider the transverse momentum of the tt pairs in the following processes: These have been studied before in the literature [22,25,31,120,161], also with a then- private version of MG5 aMC. The p T (tt) distributions behave as is typical for EW corrections dominated by EW Sudakov logarithms. Full agreement with our previous results [22,25] is found, which constitutes a further cross check of the full automation of the mixed-coupling expansion achieved in the current version of MG5 aMC. The processes of eq. (6.21) are also considered in sect. 6.3, where we include all of the LO and NLO terms, as anticipated in eqs. (6.2) and (6.3). On the right panel of fig. 9 we show the transverse momentum of the hardest jet in triple jet, single-top, and ttj production: pp −→ tj , pp −→ ttj . (6.22) As the notation suggests, in the single-top process we do not include single anti-top production. NLO EW corrections to triple-jet production are computed here for the first time. As fig. 9 shows, we find them to be small for this observable, but not entirely negligible at the upper end of the considered transverse momentum range (∼ 1 TeV), where they are of O(−10%). Up to small differences, they thus exhibit the same pattern as the EW corrections to the inclusive transverse momentum in dijet production [32,162]. We have verified that similar effects are present in the second-and third-hardest jet p T 's. Conversely, the impact of NLO EW corrections is seen to be essentially negligible on any of the two-and three-jet invariant masses (up to 4 TeV) that can be constructed from the three-hardest jet momenta.
The single-top process of eq. (6.22) includes both t-and s-channel mechanisms. Its NLO EW corrections have been computed before in the context of supersymmetric extensions of the SM [163][164][165] (a soft approximation has been employed in ref. [163] to deal with real-emission contributions). We find that EW corrections follow the typical pattern of cross-section suppression at large transverse momenta, due to Sudakov effects. As far as ttj production is concerned (which will also be discussed in sect. 6.3), the impact of the EW corrections is significantly smaller in this p T range than for single-top.
6.3 Complete NLO corrections to pp −→ tt(+V, H) and pp −→ ttj production In this section we focus on top pair production, possibly in association with either a heavy boson: The syntax of these commands has already been discussed in sect. 2. We point out that in the case of ttj production at these perturbative orders massless leptons must also be included in the definition of both the p and j multiparticles, in keeping with what is explained in appendix D. This can be done by executing the following commands: MG5 aMC> define p = p e+ e-mu+ mu-ta+ ta-MG5 aMC> define j = p immediately after the p and j definitions given at the beginning of sect. 6.1, and before the process-generation command. The computation of ttW − production would not pose any additional problem w.r.t. that of pp → ttW + ; it is not carried out here. The results for all the LO and NLO terms have already been computed with a private version of MG5 aMC for the pp → tt and pp → ttW + processes, and presented in refs. [39,40], respectively (in the latter paper, predictions for pp → tttt are reported as well). Recently, the NLO corrections to ttj production, bar for photon-induced processes, have been computed in ref. [41]. The complete NLO corrections for pp → ttZ and pp → ttH are given here for the first time.
We start by considering total rates, which we report in table 3. The first row displays the LO 1 contributions to the cross sections, given in pb. Rows 2-9 present instead all of the other contributions, as fractions over the LO 1 one, namely:  fig. 9), and are therefore slightly different from the latter (while being statistically compatible with them) -see the discussion immediately before eq. (6.12). As expected, for fully inclusive rates all contributions apart from the LO 1 and NLO 1 ones are small, with the exception of the NLO 3 term (and, to a smaller extent, of the NLO 2 one as well) in ttW + production -this constitutes a +12% correction of the LO 1 cross section, and can be understood as due to the opening of a tW scattering process, as was already suggested in ref. [40,166]. More in details, Σ NLO 3 and Σ NLO 2 are equal to about +7.6% and −2.9%, respectively, of the total NLO cross section for such a process. We now turn to presenting predictions for selected differential distributions in figs. 10-14. These figures have all the same layout. In the main frame there are eight histograms (ten in fig. 14). The solid back one is the sum of all the three (four) LO i and four (five) NLO i contributions (as in sect. 6.2, these are given as cross sections per bin), which has been denoted by "NLO" in eq. (6.13). The three (four) dashed/dot-dashed ones are the Σ LO i terms (green for LO 1 , blue for LO 2 , red for LO 3 , and yellow for LO 4 in fig. 14), while the four (five) solid/dotted ones show the Σ NLO i terms (green for NLO 1 , blue for NLO 2 , red for NLO 3 , yellow for NLO 4 , and light blue for NLO 5 in fig. 14). A dot-dashed or dotted pattern is used when the corresponding result is negative -what is displayed on the figure is then the absolute value of the cross section. In the lower insets of the figures, the ratios are shown of the individual LO i and NLO i contributions over the total NLO result, with the same patterns as those used in the main frames.     The transverse momentum of the top quark in pp → tt production is presented in fig. 10. The LO 1 contribution is dominant in the whole p T (t) ≤ 3 TeV range considered, accounting for at least 60% of the total NLO cross section. The second largest contribution is the NLO 1 one, that constitutes a correction of the LO 1 term equal to about 40% of the latter, with a rather mild dependence on p T (t) -needless to say, this is well know from standard NLO QCD results. All of the other contributions are small and equal to at most 1% of the LO 1 one. The exception is the NLO 2 term, which monotonically decreases with p T (t), becoming negative at around p T (t) 50 GeV, and growing up to −10% when p T (t) → 1 TeV. This example confirms that, when an accuracy at the percent level is required, corrections subleading w.r.t. the dominant QCD ones must be computed. In the particular case of p T (t), and up to 3 TeV, the only contribution that can be safely neglected is the NLO 4 one, which remains everywhere below the per-mille level w.r.t. LO 1 . However, it is important to point out that this conclusion is both process-and observable-specific.
In fig. 11 we show the transverse momentum of the Z boson in pp → ttZ production. The LO 1 contribution is again dominant, at the level of 70% of the total NLO cross section. For p T (Z) 0.5 TeV the NLO QCD corrections (NLO 1 ) essentially account for the  remaining 30% of the rate. However, for larger transverse momenta the NLO 2 contribution decreases very rapidly towards negative values, which can be as large as −25% of the total at p T (Z) 3 TeV. There is thus a significant cancellation between NLO 1 and NLO 2 , since the former also grows (towards larger positive values) with increasing p T 's, but slower than the latter. In general, the pattern of the impact of the subleading terms is an interesting one, in that it systematically violates the hierarchy one would naively expect on the basis of a simple coupling-constant counting. For example, at small p T 's the largest contribution among the subleading ones is that due to LO 3 , that amounts to about 2.5% of the total NLO rate, followed by LO 2 (equal to about −1% of the total). Moving towards larger p T 's the NLO subleading terms become increasingly important. Apart from the case of NLO 2 , which we have already discussed, it is worth noting at p T (Z) 2 TeV we have with all these contributions being relatively close to each other and thus featuring non-negligible cancellations (since Σ NLO 4 < 0).
The transverse momentum of the hard W + boson in pp → ttW + production is presented in fig. 12. As was the case for Z transverse momentum of fig. 11, QCD-induced mechanisms are responsible for the dominant contributions to the cross section, at both  the LO (LO 1 ) and the NLO (NLO 1 ). However, there are also notable differences w.r.t. the case of Z associated production. More specifically, we observe what follows. Firstly, the Σ LO 2 term is identically equal to zero because of colour. Secondly, for p T 's in the TeV region, the Σ NLO 1 contribution is comparable to or larger than the Σ LO 1 one. Thirdly, for p T (W ) 400 GeV the largest subleading term is NLO 3 , and in particular Σ NLO 3 > Σ NLO 2 . This is the manifestation, at the differential level, of what has been already observed in the case of fully inclusive rates in table 3. More details on this process can be found in ref. [40].
We consider ttH production in fig. 13, where we display the transverse momentum of the Higgs boson. Apart from the very large p T (H)'s, in this case subleading contributions do tend to be numerically subleading. All of them, apart from Σ NLO 2 , are well below 1% of the total NLO rate for p T (H) 1 TeV. As was already observed in fig. 9, the NLO EW corrections (NLO 2 ) are positive (3 − 4%) at small p T 's, but become negative at around p T (H) 150 GeV, and approach the −10% level in the TeV range.  order) LO 1 , NLO 1 , and NLO 2 . As is evident from the plot, in particular from the lower inset, the relative impact of the former NLO contribution increases with p T -being equal to about 5% of the total cross section at the threshold, growing significantly immediately afterwards, and reaching a value of about 30% for p T (j 1 ) 150 GeV. As far as Σ NLO 2 is concerned, it is equal to about −2% of the total cross section at the threshold. It decreases slightly up to p T (j 1 ) ∼ 100 GeV, and then increases (in absolute value) significantly, to reach values of O(−5%) at p T (j 1 ) ∼ 1 TeV. As was already observed in several of the cases discussed so far in this section, the hierarchy among the various contributions does not really follow the one based on naive coupling-constant counting (apart from the two dominant contributions), with large violations associated with increasingly subleading terms.

Conclusions
In this paper, we have studied a number of topics relevant to the perturbative computations that are accurate to NLO in the simultaneous expansion in two coupling constants, which we refer to as mixed-coupling scenario, with the final goal of applying our findings  to the automation of such computations in the framework of MadGraph5 aMC@NLO. In order to be definite, and given its importance for the current and future physics collider programmes, we have explicitly discussed the QCD+EW case, which offers the additional advantage of exposing the problems posed by an infrared sector which is by far and large maximally involved. However, the validity of our treatment is not restricted to those theories, and importantly this remark applies to the upgraded public version of the Mad-Graph5 aMC@NLO code, which we release in conjunction with this paper. In particular, mixed-coupling capabilities do not impair one of the guiding principles that underpin our software, namely that the characteristics specific to a given theory are not hardwired in the code, but loaded dynamically into it as part of a model that is fully under the user's control. One example not related to EW corrections has been presented in ref. [167], in the context of charged-Higgs studies.
In order to exemplify the procedure adopted to subtract the infrared singularities in a mixed-coupling expansion, we have explicitly extended the relevant formulae of the FKS method to the QCD+QED case. We have also shown how to re-formulate the FKS subtraction in presence of final-state tagged particles, represented by means of fragmentation functions. In doing so, we have proven the cancellation of the dimensionally-regulated singularities that emerge in the intermediate steps of the computations in an observable-and process-independent way, and expressed the final results in terms of IR-finite short-distance partonic cross sections. These features are thus fully analogous to their counterparts in the inclusive non-tagged case. We have also discussed several aspects of interest, both formal and relevant to their practical implementation, of the complex-mass scheme approach to theories with unstable particles, so that they can be readily applied also to cases other than the Standard Model, characterised by arbitrary mass spectra and possibly large decay widths.
We have computed the total cross sections and selected differential distributions for a large number of processes at the 13-TeV LHC. We have always included NLO EW effects, and in the case of tt production (typically in association with other objects) we have actually given the complete LO and NLO predictions in QCD+EW. Several of these results are presented here for the first time. From a phenomenological viewpoint, one of the most interesting consequences of our study is the fact that the numerical impact of the various subleading terms is more often than not impossible to predict with simple counting arguments based on the hierarchy of the couplings. Needless to say, the extent to which these subleading effects (and the NLO EW corrections, for that matter) will be important in practice will depend on the ability of the experiments to collect data samples with sufficiently large statistics in the relevant corners of the phase space Having the LHC phenomenology in mind, the primary application of this paper in the near future will be that of the computation of QCD+EW corrections in a systematic manner for all hadroproduction processes of interest. It is thus important to bear in mind the current limitations of MadGraph5 aMC@NLO, which will be removed in later versions; we point out that such limitations are due to a lack of either implementation work or phenomenological information, since at the formal level all of the necessary theoretical ingredients have been given either here or in previous publications. Firstly, only fixed-order predictions can be obtained, since NLO+PS capabilities are not publicly released, in spite of the fact that the MC counterterms appropriate for the MC@NLO modified subtraction of QCD-like QED showers are available. This stems from both the necessity of thorough tests against current collider data of the models employed in MCs for QED showers, and some interesting issues posed by the presence of tagged light particles. Secondly, we have not implemented the FKS subtraction with fragmentation functions. We plan to do so starting from theoretically-motivated functions associated with photons and leptons, in the hope that this will help their extraction from actual data. We also point out that the NLO predictions automatically computed by the code are in the additive scheme; this is mandatory, since the calculations are not restricted to NLO QCD and NLO EW effects, but might include further subleading terms. If a multiplicative-scheme result is desired, one must combine manually the Σ NLO 1 and Σ NLO 2 cross sections given in output by MadGraph5 aMC@NLO. Finally, applications relevant to e + e − collisions will require the implementation of ISR and beamstrahlung effects.
MG5 aMC package. SF thanks the CERN TH Division for the hospitality during the course of this work, and P. Nason for discussions on the results of ref. [72]. RF, VH and DP thank the Munich Institute for Astro-and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe" for its support. SF, VH, and HSS are grateful for the support from the organisers of the Les Houches 2017 "Physics at TeV Colliders" workshop, where interesting discussions related to this work took place. HSS thanks A. Denner for discussions about the complex-mass scheme. MZ thanks the LPTHE Paris, where he was employed during the early stages of this paper. RF and DP are supported by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award Project "Event Simulation for the Large Hadron Collider at High Precision". The work of VH has been supported by the ERC grant 694712 "pertQCD". HSS is supported by the ILP Labex (ANR-11-IDEX-0004-02, ANR-10-LABX-63). The work of MZ has been supported by the Netherlands National Organisation for Scientific Research (NWO) and in part by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodovska-Curie grant agreement No 660171.

A Altarelli-Parisi kernels and FKS charge factors in QED
We start by reporting the 4−2 dimensional forms of the unpolarised unregularised (z < 1) one-loop Altarelli-Parisi (AP henceforth) QED kernels: The AP splitting functions in QED must still obey eqs. (4.48) and (4.49) of ref. [46], namely: By using eq. (A.1) in eq. (A.7), one obtains again eq. (3.16). Conversely, eq. (A.6) can be exploited to derive the γ QED factors. In order to do so, one can use the momentum conservation condition and the QED evolution equations for PDFs. Namely: a dx x f a (x) = 1 , with the sum extended over all flavours. Eq. (A.8) leads to: Since this must hold for any set of PDFs, it implies: Eq. (A.10) can be solved explicitly by expanding the AP kernels perturbatively. At the NLO, the following formula suffices: By considering the terms proportional to α in eq. (A.10), one arrives at: dy y P QED γq i (y) + P QED q i q i (y) = 0 , (A.12) dy y P QED γl i (y) + P QED l i l i (y) = 0 , (A.13) where the range of the index i must include both fermions and antifermions. By using the forms of the QED AP kernels given before, explicit computations lead to the following results: Owing to the absence of QED contributions to gluon branchings at this order, we also set γ QED (g) = 0 ; (A. 18) note that C QED (g) = 0 from eq. (3.16). The QCD γ's factors are given in eq. (I.4.7) and eq. (I.4.8) for strongly-interacting particles (bear in mind that they are denoted here by γ QCD (I)). We also need to define such factors for photons and leptons, and obviously: These and eq. (A.18) serve the sole purpose of expressing short-distance cross sections in a compact way.
In the FKS formulae for Born-like final-state remainders, there appear factors denoted by (with the present notation) γ QCD . These are given in eq. (I.4.9) and eq. (I.4.10) for gluons and quarks respectively. As was done above, we also need to set: It is clear that an FKS subtraction of QED singularities will lead to similar quantities. One can identify them with the O( ) non-logarithmic term of Z(I) with a minus sign in frontsee eqs. (A.10) and (A.11) of ref. [46]. Explicit computations of Z(q), Z(l), and Z(γ) with AP QED kernels lead to: and trivially γ QED (g) = 0. The results of this appendix imply, in particular: with the normalization (note the sign of the QED coefficient):

B RGE invariance
In order re-instate in the cross sections the separate dependence on the different hard scales, let µ R be the QCD renormalisation scale, µ α the QED renormalisation scale, and µ F the factorisation scale. Then, the rule given at the beginning of appendix C of ref. [50] is modified as follows: If µ F = µ R or µ F = µ α , all the formulae for the short-distance cross sections given in this paper must be computed with µ = µ F , except for the argument of α S , which must be set equal to µ R , and for the argument of α, which must be set equal to µ α .
Because of this, two extra terms must be added on the r.h.s. of eq. (3.18) at the level of cross sections (that is, convoluted with PDFs and summed over p and q): By neglecting terms beyond NLO, one arrives at the analogue of eq. (I.C.4): By using eq. (A.27), one sees that a sufficient condition for the solution of eq. (B.2) is: for any p ≥ 1 and q. Analogously, the independence of µ α leads to:  4) show that both C (p,q) and D (p,q) are proportional to α p S α q , which justifies the notation adopted. In ref. [50] the RGE for µ F was used to insert the dependence on multiple scales into the virtual corrections computed with a single scale; the same operation can be performed here. By keeping only terms of LO and NLO, from eqs. (3.18) and (B.1) we obtain: The terms that feature the derivative of the PDFs are dealt with the AP equations and the following identity, where we explicitly indicate initial-state flavours: (B.6) By using the perturbative expansion of the AP kernels given in eq. (A.11), each of the first two terms on the r.h.s. of eq. (B.5) generates two contributions, corresponding to the QCD and QED evolution of the relevant incoming leg. From eqs. (3.21) and (3.22), one sees that the derivative w.r.t. µ F of the degenerate (n + 1)-body terms is determined by (T = QCD, QED): By using the identities of eqs. (A.6) and (A.7) in eq. (B.7), one has: The first term on the r.h.s. of this equation cancels the contributions due to PDFs evolution (see eq. (B.6)). The other term cancels the contribution due to dσ (C,n) (p,q) in eq. (B.5), which results from the first line on the r.h.s. of eq. (3.27). This implies that a sufficient condition for eq. (B.5) to be satisfied is: for any (p, q), which generalises eq. (I.C.8). Eq. (B.9) can now be employed in the same way as eq. (I.C.8) in appendix C of ref. [50], starting from the analogue of eq. (I.C.9): while eq. (I.C.10) is still valid. Therefore, if given a one-scale finite virtual contribution v (p,q) (M 2 ) such that: we can obtain the analogue of eq. (I.C.13): For all practical purposes, the scale dependence of α can be neglected. In the formulae above, this is achieved by formally setting β QED 0 = 0, and by replacing α(µ α ) with a constant value appropriate to the EW renormalisation scheme that is being employed.

C Symmetry factors in FKS
The symmetry factor associated with the (m + 2)-body matrix elements in FKS formulae is always equal to m!, regardless of the flavours of the m final-state partons (see eq. (4.11)). To see why this is so, one starts from the identity 41 : By introducing an index j = σ −1 (i), eq. (C.3) can be re-written as follows: which therefore is manifestly identical to a contribution to the r.h.s. of eq. (C.1). The same argument allows one to start from a term on the r.h.s. of eq. (C.1), and associate it with a term on the l.h.s. of that equation, thus concluding the proof Equation (C.1) is essentially the definition of a fully-symmetric final state, which therefore entails a symmetry factor equal to m!. However, this does not yet imply that: is equal to the quantity that enters the definition 42 of the short-distance cross sections, namely: The arguments of the matrix elements in eq. (C.1) would best be written as ordered sets, since we understand that, for any l, parton l with flavour a l has momentumk l on the l.h.s. and momentumk σ(l) on the r.h.s. However, since we are about to prove that in all cases the final state can be regarded as fully symmetric, we thought it unnecessary to introduce an ordered-set notation just for the sake of the present discussion. 42 We assume for the moment that all other factors in the definition of the cross sections are invariant under permutations of final-state partons, as in the standard FKS formulation. Other cases, including that of fragmentation, will be treated shortly.
Here, restricts the sum over all non-redundant flavour combinations (i.e. all of the partonic subprocesses). The symmetry factor is: where we have assumed that the m-body final state is partitioned into k sets, each of which composed of n i identical particles (i = 1, . . . k), so that: In order to show that eq. (C.5) is identical to eq. (C.6), we first re-write the former as follows: ) is the number of identical contributions to the sum of eq. (C.5). Such a number is equal to the number of permutations of m objects that do not leave invariant its k subsets composed of n i identical objects (because all ordered sets left invariant by a permutation that acts non-trivially only on identical elements are counted once in flavour sums such as those that appear in eq. (C.1)). To compute this, start from the total number of permutations of m objects, which is equal to m!. Choose the set generated by one of them, and consider its subset labelled by k = 1; each of the n 1 ! permutations that operates within this subset is equal to the identity. This argument is valid for each of the original permutations chosen, which implies that the number of permutations that do not act as the identity on the k = 1 subset is m!/n 1 !. By repeating this argument for each of the k subsets, one arrives at the number sought: (C.10) By inserting this result into eq. (C.9) and by taking eq. (C.7) into account, one proves that eqs. (C.5) and (C.6) are indeed identical. Let us now consider the expression: for any function g(a) (which can possibly also depend on the four-momentum of parton a) and index 3 ≤ p ≤ m + 2. Examples of eq. (C.11) are the hadron-level fragmentation cross section, eq. (4.13), or cross sections one arrives at in the intermediate steps of the FKS procedure, such as eq. (II. 4.20). By applying to eq. (C.11) the same arguments made before in this appendix one arrives at: for any 3 ≤ p ≤ m + 2, p = p. Here, there is a subtlety which is worth stressing: at variance with eq. (C.1), that is fully local in the momentum space, the fact that eqs. (C.11) and (C.12) are identical understands a phase-space integration (whence the presence of the phase-space dφ m in these equations), owing to the possible momentum dependence of the function g(). This is obviously not restrictive, since all IR-safe observables are obtained by an integration over the phase-space of the relevant short-distance quantities.
The arbitrariness of p in eq. (C.12) thus leads to: (C.13) The factor 1/m on the r.h.s. of eq. (C.13) implies that the symmetry factor associated with quantities such as X defined in eq. (C.11) summed over p is indeed equal to m!, since m! = m(m − 1)!, and (m − 1)! is the symmetry factor of X with p fixed. We point out that the identity in eq. (C.13) is used not only when dealing with fragmentation cross sections, but also in the manipulation of the cross sections that emerge in the intermediate steps of the computation in the standard FKS case -an explicit example will be given later in this appendix.
In the context of fragmentation cross sections, eq. (C.13) can also be understood in a less formal way by using simple physics considerations. Consider an m-gluon final state. A gluon is fragmented, and one is thus left with (m − 1) gluons, with an associated symmetry factor equal to (m − 1)!. It is clear that it does not matter which particular gluon fragments; one the other hand, in order to avoid counting the same contribution more than once, the sum over different fragmenting gluons need not be performed: one is sufficient. However, precisely because all possible fragmentation contributions are identical, one can symbolically write: which is what eq. (C.13) expresses more precisely 43 . Note, finally, that the r.h.s. of eq. (C.14) is also consistent with the idea that the sum over fragmentation contributions is related to final-state multiplicities. By formally thinking of a gluon as a hadron (which can be done by setting the FF equal to a Dirac delta), in an m-body configuration the gluon multiplicity has to be equal to m, and the corresponding symmetry factor equal to m!, which is precisely eq. (C.14).
We conclude this appendix with a remark on the label of the fragmenting parton, and specifically on its summation range. Owing to eq. (4.13), one starts with a realemission contribution (which has an (n + 1)-body final state) where 3 ≤ p ≤ n + 3. After the procedure that leads to the cancellation of the IR singularities, one is left with several (quasi-)n-body contributions (namely, soft (eq. (4.23)), initial-state collinear remainders (eqs. (4.26) and (4.27)), and final-state collinear remainders (eqs. (4.84) and (4.85), or eqs. (4.90) and (4.91))) for which 3 ≤ p ≤ n + 2. It is important to realise that this reduced summation range is not imposed in order to be consistent with the n-bodiness of the corresponding contribution, but rather that it emerges naturally from the computation. In order to show that this is the case, we consider the soft cross section of eq. (4.23) in order to be definite.
Firstly, we note that when dealing with the analogue of eq. (II. 4.11), the fragmentationspecific information is taken into account by using the identity: where we have exploited eq. (I.4.18) and eq. (4.12). The key point of the procedure of sect. 4.2 of ref. [46] is that the dependence upon the soft parton i is restricted to the eikonal factors, and disappears from the reduced n-body cross sections (upon relabeling). Crucially, eq. (C.15) guarantees that this remains true also in the fragmentation case, thanks to the factor 1 − δ ip . Because of this, the flavour a p of the fragmenting parton does not play any special role: it appears also in the reduced matrix elements, where its properties are never used in the manipulations of sect. 4.2 of ref. [46]. Indeed, the only parton flavour which needs a specific treatment is that of the soft one, a i , the sum over which must be carried out explicitly (and turns out to be trivial owing to a factor δ a i gsee the comment immediately below eq. (II.4.18)). The bottom line is that also in the present case one arrives at eq. (II.4.20) (with minor notation changes, and the fragmentation factor F (p) H ). At this point, however, some care is required. In the FKS procedure, one passes from eq. (II. 4.20) to the individual terms in the sum on the r.h.s. of eq. (II.4.25) by means of a relabeling, i.e. of a map R σ : with σ a permutation of n objects. The relabeling formalises the independence of eq. (II.4.20) of i, and allows one to see that the sum on the r.h.s. of eq. (II.4.25) amounts to an overall factor equal to n + 1 (equal to 4 in ref. [46]). It is the arbitrariness of σ implies that one can put eq. (II.4.27) in its fully symmetric form in terms of final-state quantities. One starts by writing the Born cross section times 44 e.g. nC(a 3 ) (and by choosing the colour-linked Born indices e.g. as k = 3 and k = 4), and then uses the procedure that leads to eq. (C.13) to arrive at the sums that appear in eq. (II.4.27).
In the fragmentation case this procedure is unchanged. However, the presence of 1 − δ ip implies that the overall factor that results from the sum over i in eq. (II.4.25) is not equal to n+1, but rather to n. Therefore, when this is combined with the symmetry factor of the real-emission matrix elements (n + 1)! (i.e. the 1/4! of eq. (II.4.20)), one obtains n/(n + 1)! rather than 1/n! (which in FKS is then interpreted as the Born-level symmetry factor, and absorbed into the Born cross sections of eq. (II.4.25)). Therefore, in the fragmentation case we can still arrive at eq. (II.4.27), up to two differences: a fragmentation factor and an overall factor n n + 1 , (C. 18) which arises from: with the 1/n! term then included in the Born cross sections as in FKS. Now observe that by construction (see eq. (C.16)) the index of the fragmenting parton in eq. (C.17) is such that 3 ≤ p ≤ n + 2, while for the original index 3 ≤ p ≤ n + 3. Owing to the symmetry of eq. (II.4.27), the property of eq. (C.12) can be exploited here as well, to relabel p as 1. This implies that the sum over p in eq. (4.13) now contains n + 1 identical terms. By performing that sum, one gets rid of the n + 1 factor in the denominator of eq. (C.18). The remaining factor n there is then cancelled by performing the procedure that leads to eq. (C.13) (by what is there the factor 1/m), procedure that allows one to re-instate a fully symmetric form in terms of the fragmentation factors. This concludes the proof that the analogue of eq. (II.4.27) is given by eq. (4.23).

D Process generation and infrared safety
We remind the reader that the generic expression of a cross section in a mixed-coupling expansion is given by the master equation (2.1). At the LO and NLO, in particular, this reads as in eqs. (2.3) and (2.4), respectively. The integer numbers k 0 , c s (k 0 ), c(k 0 ), and ∆(k 0 ) are process-specific quantities (with k 0 = c s (k 0 ) + c(k 0 ) + ∆(k 0 )), whose Born-level interpretation is apparent in eq. (2.3), namely: • k 0 is the overall power of the coupling constants combined (i.e. α n S α m is such that n + m = k 0 ) at the Born level, which is the same for all the Σ LO i contributions.
• c s (k 0 ) is the power of α S common to all Σ LO i contributions.
• c(k 0 ) is the analogue of c s (k 0 ), relevant to α.
• ∆(k 0 ) + 1 is the number of contributions to the complete LO cross section Σ (LO) .
As was already anticipated in sect. 2, the typical MG5 aMC process-generation command may read as follows: where p i denotes either a particle or a multiparticle label. The integers n max and m max help decide which Σ LO i and (indirectly) Σ NLO i will be included in the predictions, according to the constraints given in eqs. (2.8) and (2.9) for the LO and NLO cross sections, respectively. n max and m max can be freely set by the user. However, we point out that some assignments might not correspond to any physical contribution. In particular, by comparing eq. (2.8) with eq. (2.3), we obtain the following conditions: By imposing that the max and min operators in eq. (D.2) be non-trivial, and that the q range in the sum on the r.h.s. of eq. (2.3) be non-null, one further obtains: which, together with eq. (D.1), give the complete conditions that guarantee that there will be at least one Born-level contribution 45 . We point out that the keywords QCD=n max and QED=m max may both be omitted. When this is the case, MG5 aMC generates the process with the smallest possible number of QED vertices at the Born level (which is equivalent to using QED=m max with the smallest m max compatible with eqs. (D.1) and (D.3)).
In a graphical manner, the generation procedure is shown in the example of fig. 15, where each of the black blobs in the upper (lower) row represents a Σ k,q contribution at the LO (NLO). By setting n max = N + 1 as in fig. 15, one selects the set S 1 of all the LO contributions which lie to the right of the vertical blue dashed line. Conversely, by setting m max = M , the selected set S 2 is that of all the LO contributions which lie to the left of the vertical red dashed line 46 . The LO cross section that will enter the physical predictions is then obtained by summing the contributions that belong to the intersection S 1⊕2 = S 1 ∩ S 2 of the two sets obtained previously -hence, these are the terms of O(α N +1 ). The procedure just described is implemented in MG5 aMC essentially by constraining the powers of the coupling constants. Naively, one might assume that eqs. (2.8) and (2.9) are all that is needed, but in fact care must be exercised at the NLO, where one must relax the strict implementation of such constraints in the case of a mixed-coupling expansion. In order to illustrate the problem, let us consider the process generated by:

MG5 aMC> generate p p > t t~QED=0 QCD=2 [QED]
which thus computes the second-leading ("EW") corrections of O(α 2 S α) to the leading Born terms of O(α 2 S ). Among the real-emission contributions, one finds e.g. the partonic process: γq −→ ttq , (D. 5) with q any massless quark. There are two FKS sectors associated with eq. (D.5), relevant to the collinear configurations in which the outgoing quark is parallel either to the incoming photon or to the incoming quark. In these configurations, the Born processes that factorise are:q q −→ tt , γg −→ tt , (D. 6) respectively. The corresponding matrix elements are of O(α 2 S ) and O(α S α). Therefore, in spite of the fact that Born-level O(α S α) terms must not contribute to the predictions that result from the generation command given above, they must nevertheless be generated in order for the code to be able to construct the counterterms that render the cross section IR finite.
The bottom line is the following: at the Born level, the constraints of eq. (2.8) are applied to the predictions given in output, as requested by the user. Internally, those of eq. (2.9) are used instead. More precisely, the Born matrix elements that are generated for the sole purpose of constructing IR counterterms factorise a coupling-constant combination α n S α m that has 47 : • If QCD corrections are computed: n one unity larger and m one unity smaller than the corresponding exponents associated with the user-selected Born contribution with the largest α S power; • If QED corrections are computed: n one unity smaller and m one unity larger than the corresponding exponents associated with the user-selected Born contribution with the largest α power.
In the example of fig. 15, these prescriptions lead to compute the α N +2 We now turn to commenting on the role of the multiparticles in the context of NLO computations. We remind the reader that in MG5 aMC a multiparticle is a subset of the particles that belong to the loaded physics model, which may be defined by the user at runtime. For example: MG5 aMC> define p = g d d~u u~s sd efines the multiparticle p to be composed of the gluon and the three lightest quarks. Multiparticle names can be freely chosen, but in practice the two most commonly used ones, namely p and j, are conventionally associated with the incoming hadrons and the outgoing "jets", respectively, and serve to define the elementary (i.e. those that are relevant at the level of short-distance cross sections) components of these objects. Because of this, multiparticles can be effectively used at the LO 48 to control which partonic processes can contribute to the calculation of a given observable cross section. For example, with the definition of the proton given above, no processes will be taken into account that have at least one charm or bottom quark in the initial state. The situation is unfortunately more involved at the NLO, as the example of eq. (D.5) shows. If one defines the proton as a multiparticle that does include the photon but not the gluon, one generates (among others) the process of eq. (D.5) and the qq-initiated Born of eq. (D.6), but excludes the γg-initiated Born, which is essential in order to have an IR-finite cross section. We point out that this problem could be avoided by using a different multiparticle definition, to be used exclusively when constructing IR counterterms. We have refrained from implementing this option, because it would lead to IR-finite, but still unphysical predictions 49 . The immediate consequence of what we have just said is the following: • Whenever EW corrections are considered, the photon must be included in the definition of the p and j multiparticles.
An explicit example, which applies to the computations performed in this paper, is given at the beginning of sect. 6.1.
Besides the mandatory inclusion of photons in multiparticles, there might be cases when (massless) leptons must be included there as well. In general, this is the case when one is interested in computing EW corrections to processes that have initial-or finalstate photons at the Born level in the user-selected contributions. For example, leptons need not be part of the multiparticles in the calculation of the EW corrections to the tt hadroproduction considered before (the initial-state photon that appears in eq. (D.6) is relevant to a real-emission, not a Born, contribution). Conversely, they must be included in the computation of the complete NLO corrections to dijet hadroproduction, as was done in ref. [32] where QCD partons, photons, and leptons have been treated democratically.
The discussion given above on IR safety should render it clear that the inclusion of leptons in the p multiparticle serves the purpose of the correct implementation of collinear subtractions, and is thus independent of whether lepton PDFs are non-zero. In fact, although examples of sets with non-zero lepton PDFs exist (see e.g. refs. [169,170]) in practice their contributions are in general safely negligible. Therefore, if we assume f (H) l (x) ≡ 0, we may wonder which Born-level matrix elements with initial-state leptons will be generated that will give a null contribution to the cross section, and can therefore be discarded before generating them in order to increase the efficiency. This is the case when two initial-state leptons are present. In fact, if the corresponding matrix element is regarded as a physical Born contribution, it is equal to zero because of f (H) l (x) ≡ 0. Conversely, if it used to construct an IR counterterm for a real-emission matrix element, the latter will necessarily have one initial-state lepton, which will again be equal to zero because of the vanishing of the lepton PDFs. Processes with a single initial-state lepton (Xl ± → Y ) will have to be kept. In fact, although they vanish as Born contributions, they give the factorised matrix element in real-emission processes of the type Xγ → l ∓ Y , which are in general different from zero.
In order to take the above into account, the shell variable include lepton initiated processes is made available in MG5 aMC, that helps optimise the generation of lepton-initiated processes. In particular, when such a variable is set equal to False (its default value), all processes with initial-state leptons will be discarded, except for those at the Born level that feature a single lepton, and such that the realemission processes associated with them cross that lepton to the final state, and replace it with an initial-state photon. When the variable above is set equal to True, no process will be discarded.
In conclusion, whenever the following NLO corrections are computed: for any k such that 1 ≤ i − k ≤ i, we recommend to define the p and j multiparticles as summarised in table 4. Note, therefore, that the definitions are dictated by the most Processes without jets Processes with jets Physical objects PDF(qg) PDF(qgγ) PDF(qg) PDF(qgγ) i = 1 p = q g p = q g a p = q g p = q g a j(qg), γ, l, ν, j = q g j = q g massive particles i = 2 inconsistent p = q g a inconsistent p = q g a j(qgγ), l, ν, j = q g a massive particles i ≥ 3 inconsistent p = q g a inconsistent p = q g a l j(qgγl), ν, j = q g a l massive particles Table 4: Recommendations for the definitions of the p and j multiparticles in the computations of the NLO corrections given in eq. (D.7). q stands for all of the massless quarks and anti-quarks of the loaded physics model, g is the gluon, a is the photon, and l collects all of the massless charged-leptons (for example, l = e+ e-mu+ mu-ta+ ta-). In the rightmost column, by j, γ, l, and ν we denote jets, photons, charged leptons, and neutrinos, defined as explained in the text; q and g denote light quarks and gluons, respectively.
subleading term among those selected in eq. (D.7). We point out that the definitions in table 4 stem from the fact that the current version of MG5 aMC does not handle tagged photons and leptons, and thus will be modified when this limitation will be lifted. By PDF(qg) we have denoted PDF sets whose contents are limited to light quarks and gluons; conversely, PDF(qgγ) denotes PDF sets that include light quarks, gluons, and photons 50 . The inclusion of charged leptons into p and j summarises the discussion presented before on IR-safety requirements. Note that it holds regardless of whether lepton PDFs are zero or non-zero; however, in the former case the variable include lepton initiated processes can be left to its default value (i.e. False).
In the rightmost column of table 4 we list the final-state objects that can be defined; not surprisingly, their nature depends on the most subleading perturbative contribution one considers. In particular, jets must be defined democratically in quarks and gluons (i = 1), and photons (i = 2), and charged leptons (i ≥ 3); we have symbolically denoted the jets thus obtained by j(qg), j(qgγ), and j(qgγl), respectively, in order to render their particle contents explicit. When i ≥ 2, charged leptons must be defined as dressed, i.e. recombined with nearby photons; photons cannot be tagged, but only found inside jets (note that it is legitimate to have a jet composed of a single photon) 51 . Finally, when i ≥ 3 in general all 50 In the case where a PDF set features charged-lepton distributions, such leptons must always be included in the definition of p. Furthermore, the shell variable include lepton initiated processes must be set equal to True prior to the generation of the process. 51 This implies that photons enter both the jet-finding algorithm and the charged-lepton definition. There are at least two IR-safe procedures to accomplish this. The first one requires to start from the recombination of photons with charged fermions (both leptons and quarks), and then to reconstruct jets using gluons, quarks (possibly dressed), and not-recombined photons. The second procedures starts from the recombi-of the light particles must be considered in the jet-finding procedure; leptons can possibly be defined at the analysis level as jets that feature a non-zero lepton number. The recommendations summarised in table 4 must be seen as minimal. While they will always give consistent results if appropriate final-state cuts are applied, for specific processes some simplifications might be possible for the definitions of both the multiparticles and the physical objects.
E Technicalities of the complex-mass scheme E.1 A systematic test of the CM scheme implementation in the off-shell region Given the relative complexity of the derivation of the UV and R 2 counterterms for mixed QCD+EW corrections in the SM, as well as the intricacies of the expansions in the coupling constants that formally affect the masses, widths, and matrix elements in the CM scheme, it is important to provide validation techniques which are complementary to the check of IR-poles cancellation (performed by MG5 aMC for each generated process), and probe the parts of the loop amplitudes and of the counterterms that are both UV and IR finite.
In this section, we describe one such technique, that we have implemented in MG5 aMC. In essence, it consists in comparing kinematically-local results obtained in the CM and OS schemes, checking that they differ by higher-order terms (i.e. whose accuracy is higher than that one considers, which is either LO or NLO). The rationale is as follows. We have already illustrated in sect. 5.2 the strict similarities between the CM and OS schemes. We have also recalled there that OS computations are possible in the presence of unstable particles, provided that width effects can be neglected. The latter condition does not necessarily imply that all widths are set equal to zero, which is on the other hand the framework we want to work with. Because of this, we formally introduce a zero-width (ZW henceforth) setup, that we define as follows: for tree-level (L = 0) and one-loop (L = 1) matrix elements 52 . The index r numbers the unstable particles, and {p k } is a given kinematic configuration. The idea, then, is that for a suitable off-shell configuration {p k } (i.e. where the virtualities of all unstable particles are not close to the corresponding pole masses) the ZW matrix elements and their CM-scheme counterparts will differ by higher-order terms.
The meaning of "higher orders" requires a clarification in the context of a mixedcoupling scenario. Here, it implicitly understands terms associated with an increasing power in the coupling that governs unstable-particle decays which, in the SM case, coincides with α. In what follows, we thus use α in order to be definite, and we assume to work in a scheme where α is a real number (see sect. 5.4).
nation of photons with charged leptons, and then defines jets using gluons, quarks, and not-recombined photons. In the former procedure, there is no need for a jet-lepton separation, while in the latter one IR-safety demands that such a separation be imposed. 52 At variance with what was done in sects. 3 and 4, in the notation for the matrix elements we have omitted the index that corresponds to the final-state particle multiplicity, since it plays no role here.
Therefore, we shall deal with the following expressions for the ZW and CM-scheme matrix elements: with b the power of α associated with the leading Born term. The dependence on α S , or on any other couplings, is left implicit in the coefficients κ. The CM-scheme matrix elements on the l.h.s. of eq. (E.3) are understood to be computed by usingM andΓ (see sect. 5.2).
The series that appears on the r.h.s. of that equation stems from the perturbative expansion of theΓ r 's; we point out that, for the sake of the present numerical tests, it is crucial that eachΓ r obeys eq. (5.20). Conversely, given eq. (5.17), the fact that M r (at variance with Γ (0) r ) does not vanish when α → 0, and the fact that we are working at either the LO or the NLO, the quantitiesM r may be treated as fixed external parameters in the tests. The off-shell kinematic configuration used for the tests is such that 53 : where Ω r denotes the set of the decay products 54 of the r th resonance in the process considered, and we set ρ min = 10 by default. Both the list of resonances and a random kinematic configuration that satisfies eq. (E.4) are constructed automatically by MG5 aMC. At the tree level, the validation test consists in effectively verifying the following equality: κ Since the actual evaluation of the parameters κ in the CM scheme would involve the analytic Taylor expansion of eq. (E.3), for an automated numerical test is more convenient, and fully equivalent, to construct the following quantities: and verify that that they are both finite real numbers: The use of Mr and Γr in eq. (E.4) in place ofMr andΓr would lead to the same results. 54 Decay products may be resonances themselves, as is the case of a W that emerges from a top-quark decay.
In eqs. (E.7) and (E.8), λ is a real-valued parameter, and α ref denotes an arbitrary userspecified fixed coupling value, which serves as a reference. We stress that the actual value of the coupling constant used in the computations, α = λα ref , must be employed everywhere in the matrix elements, including in the expression ofΓ r (α) chosen for the unstable particle widths. The limits on the r.h.s. of eqs. (E.7) and (E.8) are computed numerically by MG5 aMC by choosing progressively smaller values of λ (equally spaced on a logarithmic scale, starting from λ = 1), and by verifying that the resulting sequences have constant asymptotes for λ → 0 -in other words, the sequences must feature vanishingly small slopes before becoming sensitive to numerical inaccuracies in matrix element evaluations.
In the remainder of this section, we shall give an explicit example of the test described above, by considering the partonic process ud → cs (i.e. one contribution to W + * production) in the α(m Z ) scheme; the quarks are taken to be massless, and a diagonal CKM matrix is assumed. We start by considering the tree-level amplitude A (0) CM , which reads 55 : with V (0)µ CM the qq → W * +µ offshell current. The corresponding ZW amplitude, defined in the same way as for matrix elements (eq. (E.1)), can be obtained from the r.h.s. of eq. (E.10) by means of the formal replacement CM → ZW. We can now expand each term on the r.h.s. of eq. (E.10) in series of α. We begin with the W propagator G µν CM : .

(E.11)
We shall be working in the Feynman gauge (ξ = 1) for the rest of this section, in which case the expression above simplifies to: (E.13) 55 Overall i factors do not play a role in what follows, and are thus ignored.
However, we stress that the conclusions drawn in this section apply independently of the chosen gauge. The offshell current V (0)µ CM differs from its ZW counterpart solely because of the complexified couplings present in the former. We choose 56 the input values for α in the CM and ZW schemes to be real and equal to each other. Thus, the only difference between the currents V (0)µ CM and V (0)µ ZW originates from the different values of the Weinberg angle, that we denote by c W and C W , respectively, in the two computations. We then have: The expansion in α of the derived parameter c W reads: .
where g ∝ √ α appears in the vertex of the current V By moving A ZW from the l.h.s. to the r.h.s. of eq. (E.17), and then taking the absolute value squared of both sides, one sees that the difference of the matrix elements that appears in the numerator of δ (0) (λ) (see eq. (E.7)) is of O(α 3 ), and thus shows that eq. (E.5) is satisfied. Actually, since the O(α 2 ) term in eq. (E.17) is purely imaginary, the O(α 3 ) contribution to the difference of the matrix elements in δ (0) (λ) is identically equal to zero, which implies that in the present example ∆ (0) = 0. This is in fact the case for all 2 → 2 scatterings, which renders the test proposed here less stringent for these processes. Because of this, our actual numerical validation has been based on processes with final-state multiplicities larger than two, as we shall explicitly show in the following.
The terms of O(α 2 ) on the r.h.s. of eq. (E.17) would cause the limit of eq. (E.8) to diverge, were they not canceled by analogous terms emerging from the difference of the one-loop matrix elements M In order to show that this is the case, we start by noting that most amplitudes are such that the difference between the CM and ZW results is at least of O(α 3 ). We denote the sum of such amplitudes by A (1)reg CM ; by definition, this quantity therefore obeys the following relationship: Hence, the contributions of A (1)reg CM and its ZW counterpart to eq. (E.8) cannot possibly lead to a divergent limit, and for this reason we call them "regular". Note that they include both unrenormalised one-loop amplitudes and UV counterterms. As far as the former are concerned, for illustrative purposes we consider here explicitly the self-energy insertion: Equation (E. 19) is based on the fact that the self-energy components Σ T /L,CM of the W boson differ from their ZW counterparts only because of the complexified coupling c W , which evidently only leads to NLO-subleading terms in α. Similar consideration hold for all other unrenormalised one-loop amplitudes, as well as for the α and wave-function UV counterterms. It follows that the only non-regular contributions to the full one-loop amplitude are due to the mass and c W UV counterterms. Thus, we write: As far as the contribution due to the mass counterterm is concerned, its expression can be formally read from the r.h.s. of eq. (E.19) with the formal replacement Σ αβ CM → g αβ δm 2 W which, in a ZW computation, becomes Σ αβ ZW → g αβ δM 2 W . Therefore: Here, we have substituted the mass counterterms with their explicit expressions, given in eqs. (5.11) and (5.15). In eq. (E.21) we have used the fact that the real parts of the mass counterterms in the ZW and CM setups differ by higher-order terms, provided that the analytical continuation in the latter scheme is performed appropriately (see eq. (5.21) and sect. 5.3). As for the imaginary part of Σ(M 2 W − iΓ W M W ), it has been expressed in terms of the LO W -boson total decay width Γ (0) W by using eq. (5.18). This underscores the importance of employing an expression forΓ W which obeys eq. (5.20) when evaluating eq. (E.8) numerically.
We now turn to considering the second non-regular contribution on the r.h.s. of eq. (E.20), namely that stemming from the UV counterterm δc W : where the last equality assumed A (0) ZW . Note that the number of insertions of the counterterm δC W is always equal to the power of C W factorised in A (E.23) We can now proceed analogously to what has been done in eq. (E.17). Namely, one moves the two ZW amplitudes from the l.h.s. to the r.h.s. of eq. (E.23), and then computes the absolute value squared of both sides. In this way, one shows that the linear combination of matrix elements that appears in the numerator of δ (1) (λ) is of O(α 4 ), thus proving that ∆ (1) is a finite number. We conclude this section by presenting the numerical results of the tests advocated here for a couple of representative processes. Such tests can be run directly from the interactive MG5 aMC interface by issuing the command: MG5 aMC> check cms <process definition> <options> Automated CM scheme check for process u d -> e + ν e γ |δ (1)  Automated CM scheme check for process g g > e + ν e µνµ b b - where process definition follows the usual [17] tree-level or loop-level MG5 aMC syntax (depending on whether one wants to evaluate the limit on the r.h.s. of eq. (E.7) or that of eq. (E.8)); further options can be added 57 . An explicit example of this is:

MG5 aMC> check cms u d~> e+ ve a [virt=QED] --tweak=alltweaks
In order to be definite, we limit ourselves to considering here the one-loop test of eq. (E.8), which is obviously more involved than its tree-level counterpart of eq. (E.7). We do so for the two processes ud → e + ν e γ and gg → e + ν e e −ν e bb. As was already said, the ∆ (i) 's of eqs. (E.7) and (E.8) are computed by evaluating the arguments of the limits on the r.h.s. of those equations (denoted there by δ (0) (λ) and δ (1) (λ), respectively) for increasingly smaller values of the scaling parameters λ. The results are then plotted for each value of λ, as is done here in fig. 16; in this way, one can easily see the behaviour of the relevant matrix element combination for λ → 0.
The left and right panels of fig. 16 display the results relevant to the ud → e + ν e γ and gg → e + ν e e −ν e bb processes. In the main frame, δ (1) (λ) is shown, while the inset presents λ δ (1) (λ) , which must tend to zero in the case of a successful test. Both quantities are computed in different scenarios, each of which corresponds to one curve in the main frame and one curve in the inset of fig. 16. More in detail, the three black curves are obtained by setting:Γ for the solid, long-dashed, and short-dashed patterns respectively. Equations (E.24)-(E.26) understand that these settings for the widths apply to all relevant unstable particles (i.e., the W boson for ud → e + ν e γ, and the top quark, the W boson and the Z boson for gg → e + ν e e −ν e bb -we point out that in the latter process all resonant and non-resonant contribution are included). The widths of eqs. (E.24)-(E.26) obey eq. (5.20); fig. 16 clearly shows that the three choices differ in the large λ region, but converge to the same value when λ → 0, and ultimately lead to a successful test. Conversely, by setting: that correspond to the solid red and green curves, respectively, the tests do fail, with δ (1) (λ) departing from a constant behaviour for λ 10 −3 . We point out that this occurs for a mere 1% deviation from eq. (5.20), and demonstrates the high numerical sensitivity of the present checks.
Finally, the blue solid curves in fig. 16, labelled there as log ± → log ∓ , correspond to an alteration of the CM-scheme implementation that consists in changing the selected Riemann sheet (see eq. (5.23)) in the logarithmic part of the UV mass counterterms of unstable particles. This leads to a violation of the consistency relation of eq. (5.21), which results in mis-cancellations much more severe than those induced by a small modification of the LO term in the functional form of the widths. This provides further evidence that any inconsistency in the CM-scheme implementation is very unlikely to be undetected when the present numerical validation procedure is carried out 58 . Having said this, there is a subtle point that is worth making. Namely, that by verifying that eq. (5.21) is fulfilled we establish that the appropriate Riemann sheet is used in the limit in which all widths vanish. This is not necessarily the same sheet as that relevant to the physical (i.e. with non-zero widths) configuration, be it either because (in the language of sect. 5.3) the trajectory is "long" (i.e.Γ M is not fulfilled), or because itsγ = 0 endpoint is not close to an OS-like configuration (i.e.Γ i M i is not fulfilled for some i). This does not happen in the SM, but there is no reason that prevents it from happening in an arbitrary model, where thus we do not expect the test outlined in this appendix to be particularly effective in detecting an incorrect analytical continuation. 58 We remark that the plots of fig. 16 differ from those generated automatically by MG5 aMC only in their layout, improved here for the sake of clarity. In order to get the curves that correspond to the incorrect choice of Riemann sheet for the logarithms and to the various functional forms of the widths, the option --tweak=alltweaks must be specified. The gg → e + νeµ −ν µbb process also necessitates to force quadruple precision (and a single helicity for increased speed) with --CTModeRun=4 --helicity=1.
In conclusion, we have successfully run the automated numerical test discussed in this section for a variety of processes, in both the G µ and the α(m Z ) CM schemes, making sure that those collectively involved all possible SM unstable particle. Together with the customary IR-pole cancellation check, we have thus validated the MG5 aMC implementation of the CM scheme in the SM.

E.2 On the numerical study of trajectories and contour integration
In this section, we briefly sketch two numerical approaches to the study of the trajectory T of eq. (5.45). In particular, our final goal is that of determining the difference n +− − n −+ that appears in eq. (5.46). We shall also briefly comment on the corresponding contour integration.
We start by observing that, irrespective of the winding number of T , the real functions T (γ) and T (γ) of the real variableγ ∈ [0,Γ] are single-valued. Let: 0 < ζ 1 < ζ 2 < . . . < ζ n <Γ , (E. 29) be the n zeros of T (γ) in the interval (0,Γ): Define: Although numerical routines that find the zeros of a real function are fast and efficient, it is desirable to have a less rigorous, but much quicker, alternative to the procedure outlined so far. One possibility stems from approximating the imaginary part of the trajectory with a sequence of segments. In other words, one introduces: with N σ a given integer number. The idea is that of replacing T (γ) for σ i−1 ≤γ ≤ σ i with a straight line that connects the two points T (σ i−1 ) and T (σ i ), i.e. with: One then finds the zero of T (γ): Note that for this solution to be acceptable, we must have σ i−1 < ζ i < σ i , which happens when T (σ i ) and T (σ i−1 ) have opposite signs. This is nothing but eq. (E.35) (for a given i). In fact, as the notation used here suggests, the roles of the σ i and ζ i quantities This procedure gives a correct result if the partition of the range [0,Γ] achieved by eq. (E.44) is sufficiently fine-grained to capture the behaviour of T (γ). This is obviously the case for N σ → ∞, which suggests that a self-diagnostic test (still not fully watertight) is that of repeating the computation of eq. (E.47) by increasing N σ -the same result must be obtained. We point out that this constitutes an automated test. In a case-by-case situation, one can always check visually that the discretisation of the trajectory according to eq. (E.45) represents a continuous deformation of the original curve that does not cross the origin. Indeed, this is what has been done with the trajectories considered in sect. 5.3, which have been dealt with by setting N σ = 100.
We point out that eq. (E.47) is equivalent to the straight-trajectory method when N σ = 1. Because of the leftmost property in eq. (5.44), a difference between the results of eq. (E.47) with N σ = 1 and with a sufficiently large N σ implies the failure of the straighttrajectory approach, which can only be due to the fact that the given trajectory cannot be continuously deformed into a straight line without crossing the origin.
We now turn to discussing the contour integration that has led us to the results shown in fig. 4. In sect. 6.6 of ref. [94], the authors have proposed a strategy for integrating I = 1 0 dx log − (V (x)) with V (x) = ax 2 + bx + c, where the polynomial coefficients a, b, and c can assume arbitrary complex values. This allows one to directly evaluate the most general bubble function B 0 (p 2 , µ 2 1 , µ 2 2 ) with all arguments complex (under certain restrictive conditions -see eq. (81) of that paper, and eq. (37) there for the definition of log − (); note that the latter differs from the log − () function of eq. (5.26)). In the complex plane, the branch cut of log − () relevant to the integral I is the set of all pointsx that satisfy the two conditions V (x) = 0 and V (x) > 0. The former condition results in two branches of an hyperbola, and the latter condition selects a subset of points (one of which possibly empty) in each of these branches; we denote such subsets by H 1 and H 2 . If either (or both) H 1 or (and) H 2 crosses the real axis in the range [0, 1], then I must be computed by deforming the x ∈ [0, 1] range into a contour C that lives in C. The contour must be such that it crosses neither H 1 nor H 2 , i.e. C ∩ H 1 = ∅ and C ∩ H 2 = ∅. The authors of ref. [94] have proposed a general parametrisation for C, given in terms of the coefficients a, b, and c. We have found that this parametrisation works as expected if both of the sets H i ∩ [0, 1], i = 1, 2 contain at most one point. Conversely, if for either i = 1 or i = 2 the set H i ∩ [0, 1] contains two points, then C ∩ H i = ∅, and C is therefore not suited to numerical integration. One actual example of this situation is the case of configuration E of table 1. We have addressed this issue by using different parameters w.r.t. those that emerge from the proposal of ref. [94] 61 , while keeping the same functional parametrisation of C suggested there. With these modified parameters we have obtained the results labelled as "revised contour" in sect. 5.3.

E.3 SMWidth: an SM decay-width calculator at the NLO QCD+EW accuracy
As was documented in sect. E.1, in off-shell regions widths can be computed at the LO without spoiling the overall NLO accuracy of the calculation. However, this is not true in general, and NLO-accurate widths must be used. In the SM, the widths of the top quark [171][172][173][174][175][176][177], W boson [178][179][180][181][182][183][184], Z boson [181,182,[185][186][187], and Higgs boson [188] have been known for a long time at the NLO accuracy in QCD and EW. In spite of this, no single public tool provides one with all of these decay widths at such an accuracy. We have amended this, by writing a self-contained package, dubbed SMWidth, which can be used either standalone or within the MG5 aMC framework. SMWidth computes the top, W , and Z total widths 62 from first principles, and calls HDecay [188,189] to obtain the Higgs width. In what follows, we briefly describe its workings, with the standalone usage presented in sect. E.3.1, and that within MG5 aMC in sect. E.3.2. Sample results are given in sect. E.3.3. Before going into that, we introduce some general features of the code.
SMWidth is written in Fortran90, and is included in the MG5 aMC distribution. It implements the analytically-integrated decay amplitudes of the top quark, W boson, and Z boson, by including NLO QCD and EW corrections. Such amplitudes have been generated by FeynArts [190], and the relevant loop and phase-space integrations have been carried 61 In the notation of eq. (91) of ref. [94], the chosen parameters for dealing with benchmark point E are: α1 = 0, α2 = 0, β1 = −0.05, β2 = 0, and αc = 0.7. 62 Widths relevant to specific decay channels can be computed only by means of calls to functions which are not part of the standard user interface.
out by means of an in-house Mathematica module, whose output has been converted into the Fortran90 code; in this way, only elementary numerical operations 63 are performed by SMWidth, and the program has thus a negligible CPU load. The top-decay width is computed according to the following formula: with n + m = 3, n ≥ 0, m ≥ 1, and: 49) As was done in appendix E.1, the subscripts CM and ZW denote CM-scheme and zerowidth-setup results 64 ; for the definition of the latter, see eq. (E.1). By W and W * we have denoted a final-state on-shell W boson, and its intermediate off-shell counterpart, respectively. D W denotes the set of the fermion pairs into which a W boson can decay. The superscript "NLO" that appears on the r.h.s. of eq. (E.48) imply that the corresponding process includes one-loop and real-emission corrections stemming from QCD or EW interaction vertices. The superscripts "(0)" in eq. (E.49) indicate that the two widths are computed at the tree level. Equation (E.48) can be easily generalised to take into account the decays of the top into down-type, non-bottom quarks.
As far as the Higgs boson is concerned, SMWidth simply acts as a calling interface to HDecay, which is embedded in our package. The user is thus expected to also cite the relevant HDecay papers [188,189] in that case.

E.3.1 Standalone usage
We now turn to describing the standalone usage of SMWidth. The package contains a module, test.f90, which acts as the driver when working in standalone mode. The code is compiled with the shell command: > make -f makefile_test We assume a gfortran compiler with the same characteristics as that employed when using MG5 aMC. The test program embedded in test.f90 can be run by executing: 63 OneLoop [114] is called on-the-fly as a library that returns the value of the C0 function, necessary in the computation of EW corrections to the top and W widths. 64 SMWidth uses OS-type renormalisation conditions and inputs. The computation of the first term on the r.h.s. of eq. (E.49) is then performed by complexifying both the W and Z masses according to eq. (5.7) (by using the internally computed LO values of ΓW and ΓZ ) and the relevant coupling factors. Note that this is done withMi = Mi andΓi = Γi, which is consistent with the fact any difference between CM-and OS-type input parameters would be beyond accuracy in eq. (E.49), owing to eqs. (5.17) and (5.19).

> ./ test
There are two input files, one that sets the EW renormalisation scheme to be adopted, and one that sets the physics parameters 65 . The former file, named scheme.dat, contains a single integer, equal to 1(2) for the α(m Z ) (G µ ) scheme. In the standalone mode, the use of this file can be bypassed by setting the global variable Decay scheme equal to either 1 or 2. As far as the physics inputs are concerned, they must be written in a file phyinputs.dat, which the Fortran subroutine ReadParamCard is tasked to read as follows: CALL ReadParamCard ( ' phyinputs . dat ') As is suggested by the call above, the name phyinputs.dat can be freely changed by the user; what matters is that it contains the basic physics input parameters, namely the masses of the top quark and of the W , Z, and Higgs bosons, and the value of either α(m Z ) −1 or G (Gµ) µ , depending on which EW scheme is chosen. We assume α S (M BW Z ) = 0.1184 (with M BW Z = 91.1876 GeV, see eq. (6.4)), and obtain the value of α S (M t ), α S (M Z ), and α S (M W ) employed in the calculation of the NLO QCD corrections to the top, Z, and W widths by means of a two-loop evolution. Two templates, relevant to the two valid EW-scheme choices, are included in the package, param card MZ.dat and param card Gmu.dat. As these names suggest, they are derived from the param card.dat file of MG5 aMC.
The W -decay width can be computed as follows: width = SMWWidth ( qcdord , qedord , finitemass ) The settings of qcdord and qedord determine whether QCD and EW corrections, respectively, are included -valid inputs are 1 to include the corrections, and 0 to exclude them. By setting the entry finitemass equal to .true. the effects are included due to the nonzero bottom, charm, tau, and muon masses. Otherwise, when finitemass=.false. all of these fermions are treated as massless.
In a manner fully analogous to the case of the W , the Z decay width is computed as follows: width = SMZWidth ( qcdord , qedord , Wrad , finitemass ) The entries are the same as those relevant to the W decays, and have the same meanings. On top of those, the parameter Wrad can be set equal to .true. in order to include in the computation the contributions due to the channels Z → W ± + f 1 +f 2 (where the W is on-shell and the fermions are treated as massless, regardless of the value of finitemass). These contributions are typically very small, owing to phase-space suppression.
The top-quark decay width can be calculated as follows: width = SMtWidth ( qcdord , qedord , finitemass , wwidth ) 65 While such parameters are fully under the user's control, one must not vastly depart from the measured SM values, lest SMWidth give incorrect results. In particular, no new decay channel must open.
The entries qcdord and qedord have the same meanings as before. finitemass is used here in order to include (when set equal to .true.) non-zero bottom mass effects -the bottom is regarded as massless when finitemass=.false.. The other quarks and the leptons are always treated as massless. Finally, by setting wwidth=.true. one is able to include the contribution due to a non-zero W width, which is an essential ingredient in order to attain NLO EW accuracy for a process that features intermediate top quarksthis corresponds to the contribution of eq. (E.49). Finally, the width of Higgs boson is computed as follows: with idummy is an INTEGER-type dummy argument. As was already said, the function SMHWidth simply acts as an interface to HDecay, with the input parameters given in param card.dat written onto the HDecay input file hdecay.in.

E.3.2 Usage in MG5 aMC
The usage of SMWidth from within MG5 aMC can proceed either in an interactive manner, through the command line, or implicitly, through suitable settings in an input file. In the former mode, the syntax is as follows:

MG5 aMC> compute widths <particle(s)> --nlo
where <particle(s)> is the list of the labels (according to the MG5 aMC syntax) for all of the particles whose widths one wants to compute. For example:

MG5 aMC> compute widths w+ z h t --nlo
The option --nlo specifies that the widths have to be computed by including both NLO QCD and EW corrections. In terms of the calls to the internal SMWidth routines described in app. E.3.1, this corresponds to qcdord=1 and qedord=1. In the implicit mode, the physics-input file must be named param card.dat. The correct template for the latter is chosen automatically by MG5 aMC upon loading the physics model (see below). In turn, its contents (in particular, the presence of the numerical value of either α(m Z ) or G (Gµ) µ ) tell the code which EW scheme must be adopted; thus, scheme.dat is irrelevant in this mode of operation. Finally, one must use the keyword Auto@NLO in the decay block of param card.dat, as per the following example: MG5 aMC will then compute the widths for the top quark and the heavy bosons at the NLO QCD+EW accuracy, by calling SMWidth internally. This option is only available after importing one of the NLO-EW-compatible UFO models available in MG5 aMC, namely either loop qcd qed sm (for the α(m Z ) EW scheme) or loop qcd qed sm Gmu (for the G µ EW scheme). For consistency with these models, light-fermion and b-mass effects are always neglected; the contributions of the channels Z → W ± + f 1 +f 2 to the Z width are also ignored. Conversely, finite-W -width effects in top decays are included. In terms of the calls to the internal SMWidth routines described in app. E.3.1, these settings correspond to finitemass=.false., Wrad=.false., and wwidth=.true..

E.3.3 Illustrative numerical results
We conclude this section by reporting sample results obtained with SMWidth. As was already said, in the case of the Higgs boson they are obtained by an implicit call to HDecay.
The values of the input parameters in the α(m Z ) scheme are given in table 5. The lightfermion masses have been set equal to the values reported in table 6. The renormalisation scale has been set equal to the mass of the mother particle, and a two-loop running has been used for α S . The corresponding results are reported in table 7, where we also show a breakdown into the various types of contributions to the top quark, W , and Z widths according to the following definitions: with δ i the fractional corrections due to NLO QCD effects (i = α S ), NLO EW effects (i = α), finite fermion-mass effects (i = m f , included only at the LO), and finite W -width effects (i = Γ W -thus, Γ LO t δ Γ W is equal to the r.h.s. of eq. (E.49)). In the G µ scheme, the same parameters as in tables 5 and 6 are employed, with the exception of α(m Z ), which is replaced by G (Gµ) µ . This in turn leads to a different value for the QED coupling constant, which we denote by α Gµ . The situation is summarised as follows: The corresponding results are presented in table 8.
We have checked numerically that the contribution of the channels Z → W ± f 1f2 to the total Z boson width is negligible owing to phase-space suppression, in spite of being     formally of O(α 2 ). As far as the top decay is concerned, we see that δ α and δ Γ W give comparable contributions (in absolute value), which thus must be both taken into account in order to have a sensible result at the NLO EW accuracy. Conversely, the impact of δ m f is significantly smaller (for the top) or negligible (for the W and the Z), and compatible with being thought as due to power corrections (m f /M t,W,Z ) 2 . Where possible, the results of SMWidth have been cross-checked with those available in the literature [89,177]. We also point out that MG5 aMC is itself capable of computing the W , Z, and top-quark total widths with the same accuracy as SMWidth (obviously, with running times longer than those of the latter); we have verified that the predictions of the two codes are in excellent agreement with each other.