An automated subtraction of NLO EW infrared divergences

In this paper a generalisation of the Catani–Seymour dipole subtraction method to next-to-leading order electroweak calculations is presented. All singularities due to photon and gluon radiation off both massless and massive partons in the presence of both massless and massive spectators are accounted for. Particular attention is paid to the simultaneous subtraction of singularities of both QCD and electroweak origin which are present in the next-to-leading order corrections to processes with more than one perturbative order contributing at Born level. Similarly, embedding non-dipole-like photon splittings in the dipole subtraction scheme discussed. The implementation of the formulated subtraction scheme in the framework of the Sherpa Monte-Carlo event generator, including the restriction of the dipole phase space through the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α-parameters and expanding its existing subtraction for NLO QCD calculations, is detailed and numerous internal consistency checks validating the obtained results are presented.


Introduction
As Run-I of the LHC has been successfully completed, culminating in the celebrated experimental confirmation of the existence of the Higgs boson, Run-II proceeds its data-taking at the unprecedented centre-of-mass energy of 13 TeV. As the much anticipated discovery of signals of beyond-the-Standard-Model physics is still lacking, precision tests scrutinising the Standard Model are of prime importance, now and in the foreseeable future. At the same time, new physics searches are looking for increasingly small signals demanding more precise estimates of the Standard Model backgrounds. This expansion of sensitivity of both precision measurements and new physics searches in the multi-TeV region demand an immense improvement of theoretical predictions.
This precision can be achieved by the inclusion of next-to and next-to-next-to-leading order (NLO and NNLO) corrections in the strong coupling and next-to-leading order electroweak (EW) corrections. Here it should be noted that both NNLO QCD and NLO EW corrections are expected to be of a similar magnitude for inclusive observables as numerically α 2 s ≈ α. On selected differential distributions, however, electroweak corrections can grow much larger. They are dominated by photon emissions in the distributions of final state leptons, for example. In invariant mass spectra of lepton pairs below a resonance, for example, O(1) correc-tions can be present, in which case a proper resummation should be included [1]. Similarly, looking at the (multi-)TeV regime, the NLO EW corrections quickly grow considerably, reducing cross sections by a few tens of percent, due to the emergence of large electroweak Sudakov corrections arising as the scattering energies Q 2 m 2 W [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. In this regime they are larger than even the NLO QCD corrections in many cases and their omission becomes the dominant uncertainty in experimental studies and searches.
To this end, benefiting from the well-established techniques developed for the automation of NLO QCD corrections many NLO EW corrections have been calculated recently . To fully automate these computations at NLO EW accuracy in a Monte-Carlo framework all infrared divergences need to be regulated, where various incarnations of subtraction methods have proven to be the methods of choice for practical implementations [39][40][41][42][43]. Similar subtractions have also been published for NLO EW calculations using a dipole picture [1,[44][45][46]. Contrary to the QCD case, only the implementation of [46], restricted to photon emissions of fermions, is publicly available though.
Besides the generalisation to all divergent splittings at O(α), including photon splittings and photon emissions off massive scalars and vector bosons, this publication addresses the issue of automatically detecting simultaneously occurring QCD and QED singularities and subtracting them consistently. These occur as soon as the Born process is defined at multiple orders O(α n s α N −n ). In this case the NLO EW correction to the O(α n s α N −n ) process, being of O(α n s α N −n+1 ), is at the same time the NLO QCD correction to the O(α n−1 s α N −n+1 ) process and will in general exhibit the corresponding infrared singularities. Further, matters of the organisation of the contributing partonic processes and their mapping to reduce the computational complexity along with the provision of infrared safe phase space cuts are discussed. The algorithm is implemented in the AMEGIC [47] matrix element generator which is part of the SHERPA [48] Monte-Carlo event generator framework. It bases on the automated subtraction of massless NLO QCD divergences therein [49,50]. The implementation presented in this publication has, in various preliminary forms, already been used to calculate electroweak corrections to a multitude of important signal and background processes [19,21,25,31,32,34,37,[51][52][53], highlighting its versatility.
The present paper is structured as follows: first, in Sect. 2 the Catani-Seymour dipole subtraction method is reviewed and the general modifications to its differential and integrated subtraction terms are discussed. Section 3 then details its automation in SHERPA's matrix element generator AMEGIC, highlighting the necessary changes and improvements with respect to [49,50]. This section also discusses various options implemented for the incorporation of photon splittings, general infrared safe fiducial phase space definitions and flavour scheme conversions. Essential cross checks validating the presented implementation are then presented in Sect. 4 before concluding in Sect. 5. Explicit formulae for all differential and integrated dipoles are given in Appendix A-C.

Catani-Seymour subtraction at NLO EW
In order to be applicable to NLO EW calculations the wellknown Catani-Seymour dipole subtraction [40,41] needs to be recast in a suitable form. To highlight the changes from the original formulation for NLO QCD calculations the complete structure of the formalism is reviewed. This subtraction formalism starts from the the expectation value of any infrared safe observable O described at NLO accuracy through Therein, the Born term B consist of the squared matrix element |M m | 2 = m s 1 , . . . , s m |s 1 , . . . , s m m with helicity states s n and further includes all parton densities, parton fluxes, symmetry and averaging factors. The virtual and real corrections, V and R, as well as the collinear counter term, C, are defined analogously. When regulating their respective divergences through dimensional regularisation, they have to be evaluated consistently in d = 4 − 2 dimensions for all singularities to cancel. Only after their summation can the limit → 0 be taken. d (4) n and d (d) n are the four and d dimensional phase space element. As V and C on the one hand side and R on the other are defined on phase spaces of different parton multiplicity, Eq. (2.1) cannot be used for numerical evaluation straight forwardly and it is rewritten as collinear counterterm on the other side are cancelled separately by the integral of the subtraction term over the oneparticle phase space, rendering the integrals of the second line finite in four dimensions as well. Section 2.1 now describes the construction of the differential subtraction term, D, used to subtract all divergences from the real emission correction, while Sect. 2.2 presents the integrated subtraction terms, I D = d 1 D, that subtracts all divergences from the virtual corrections and the collinear counterterm in its explicit Laurent expansion after being analytically integrated over the factorised one-particle phase space. All infrared divergences that occur at NLO EW are of QED origin. No subtractions of potentially large, but finite, corrections involving the emissions of real and virtual massive electroweak gauge bosons will be considered. The practical implementation described in Sect. 3 follows the general lines of [49,50].

Differential subtraction terms
To describe all singular limits of a given real emission matrix element R, related to the real emission term of Eq. (2.1) in a similar fashion as the Born term, it is decomposed into a sum over dipoles D [40,41] as Therein, i is the emitter in the final state, j is the emittee, k is the spectator in the final state, a and b are the initial state partons. Each dipole thus encodes the singularity structure caused by the emission of j in the presence of the charge of the spectator. While the divergence associated with the soft emission of j off the dipole ıj −k is partially fractioned into a piece associated with a splitting ıj → i + j in the presence of spectator k and a piece where i and k swap their roles, the divergence associated with the collinear emission of j off ıj is recovered through charge conservation once all dipoles having ıj as emitter are summed over. All four dipole types are depicted in Fig. 1. The individual dipoles take the form [40,41,54,55] (2.5) The Q ıj and Qk are the charges of the emitter and the spectator and their θ ıj/k are 1(−1), if they are in the final (initial) state. Of course, Qk = Q k and θk = θ k . In the case of photon splittings no soft divergence is present. Thus, these splittings have no dipole character. To include them in the dipole-formalism nonetheless and to distribute the recoil in splittings away from the collinear limit, spectator partons need to be assigned. As their only role is to absorb transverse momentum of the splitting process, any other particle may be considered as spectator. Each thus assigned recoil partner may be assigned a weight κ ıjk , only their sum is constraint by Eq. (2.5) in order to add up to the correct collinear limit. Various options to assign recoil partners are implemented, they are detailed in Sect. 3.2.
Since the QED charges are real numbers, the chargecorrelator simply multiplies the matrix element and only leaves the spin-correlators V i j,k , V a i j , V a j,k and V a,b j as insertions in the spin-correlated underlying Born matrix elements. The spin-correlators directly correspond to their QCD counterparts and are detailed in Appendix A. It also defines the initial state momentum rescaling parameters x i j,a , x aj,k and x aj,b , as well as the splitting variables y i j,k , u j and v j . The {α dip } = {α FF , α FI , α IF , α II } parameters serve to restrict the phase space where the individual dipole terms are non-zero and therefore need to be evaluated [54,55]. They are constructed such that for every α ıjk > 0 the singularity is fully subtracted. The introduction of a parameter κ in dipoles where a final state photon splits into a massive fermion pair in the presence of a final state spectator similarly allows a redistribution of finite terms, cf. Appendix A.

Integrated subtraction terms
By the above construction, the integral of the subtraction terms D over the one-particle phase space possesses all the necessary poles to render the second line in Eq. (2.2) finite as → 0. This section now summarises the formulae and findings of [40,41], translated to the QED case, and discusses their important features. The integrated subtraction terms, together with the collinear counterterm, are commonly reorganised into I, K and P operators through the following identification The I operator contains the necessary infrared poles to cancel all divergences of the virtual correction such that, after summing both, the limit → 0 can be taken and the integral can be evaluated in four dimensions. The K and P operators are infrared finite by construction and, thus, evaluation in four dimension is unproblematic as well. Please note, that due to the fact that, contrary to the colour correlator in QCD, the charge correlator in QED is a simple real number and thus the convolution of the Born matrix element with the respective insertion operators in QCD becomes a trivial product of the Born matrix element and the operators in QED. The spincorrelation that was still present in the real subtraction terms has been integrated out in full analogy to the QCD case. In the following, the structure of all three operators in the QED case is discussed. The I operator. The I operator contains all flavour-diagonal endpoint contributions and cancels all divergences present in the one-loop matrix elements. It takes the general form and contains single and, in case of massless emitters, double poles in . Due to the presence of such poles a dependence on the regularisation scale μ 2 enters. It is commonly identified with the renormalisation scale μ 2 R . The I operator is further dependent on the choice of the {α dip } and κ parameters that effect non-singular terms only. It is further decomposed into dipoles [41] wherein i labels the emitter and k the spectator. The charge insertion operator, which is a trivial real number in the QED case, is defined in Eq. (2.5). The full crossing invariance of the QCD I operator may be broken in the QED case in the presence of photon splittings as their recoil partner assignment is arbitrary and may involve information on initial or final state particles. Different possible choices are discussed in Sect. 3.2, some of which may break this crossing invariance. The divergences of the I operator are encoded in the functions V ik and i . While the former contains all soft-(quasi-)collinear divergences the latter contains the pure (quasi-) collinear ones. They do not only differentiate whether the emitter is a photon or not, but also between different spins of the emitter. Their precise form as well as the the flavour constants γ i and K i are given in Appendix B. A I ik encodes the dependence on the phase space restriction of the individual dipoles {α dip }. Finite terms originating in dipoles involving initial state legs, however, can be pushed into the K operator. Thus, A I ik by convention only depends on α FF . Its precise form is given in Appendix C. The K and P operators. The K and P operators collect all pieces of the integrated dipole terms that are not collected in the I operator and combines them with the collinear counterterms C to give a finite result as → 0. By construction they contain only remainders of splittings where either the emitter or the spectator is in the initial state. Thus, they are comprised of terms arising due to the change of the flavour or the partonic momentum fraction x of an initial state due to a splitting.
The K-operator is given by [41] (2.9) It depends on the partonic x, and the flavour change from the Born process initial state flavour a to a of the convolution Eq. (2.6). Therein, the K collect universal terms present in all splitting involving an initial state as either emitter or spectator. Then, while K contains solely remaining terms from final state splittings in the presence of initial state spectators, the K t are their counterparts for initial state splittings in the presence of a final state spectator, i and k running over all final state partons in each case.K contains solely related correlations between both initial states, arising from dipoles where both the emitter and the spectator are in the initial state. The A K terms collect all finite terms arising when any of α FI , α IF or α II is different from unity, thus restricting the phase space of the respective dipoles. Again, the Q 2 ik are the charge correlators of Eq. (2.5). Finally, K FS contains the factorisation scheme dependence. Currently, both only the MS schemes is supported, setting these terms identically zero.
The P-operator now collects the remaining initial state collinear singularity from all dipoles involving initial states either as emitters or as spectators and cancels them against the collinear counterterm. Through this counterterm a dependence on the factorisation scale enters. The P operator is given by [40,41] (2.10) Only initial state splittings are present, either in the presence of a spectator in the final state, which is encoded in the sum of k, or with the opposite initial state b acting as the spectator. It otherwise only depends on the Alterelli-Parisi splitting function detailed in Appendix B.

Implementation
The implementation of the QED generalisation of the Catani-Seymour dipole subtraction scheme in SHERPA's matrix element generator AMEGIC proceeds along the lines of [49]. As in general real and virtual corrections of O(α n s α m ) contain divergences of both QCD and QED origin, both cases are included in this section. In the following, the general structure of the implementation is reviewed.

Identification of dipoles
The starting point to construct the involved subtraction terms in the Catani-Seymour subtraction formalism is a given flavour configuration in the Born or the real emission phase space and the perturbative order O(α n s α m ) in accordance with the respective virtual or real correction to be computed. For all parts, on-the-fly variations of both the factorisation scale μ F and the renormalisation scale μ R are available through an extension of the algorithm detailed in [56]. Differential subtraction terms. A given real emission configuration {ab} → {1, .., m + 1} at order O(α n s α m ) can in general exhibit both QCD and QED divergences simultaneously. The following therefore describes the identification of both types of dipoles. Thus, all triplets {i, j, k} that can be built from the external particles of the process are tested for the presence of an infrared divergence, QCD or QED, by checking for the existence of a dipole subtraction term. In these triplets i and k may be in the initial or final state while j may be in the final state only. Likewise, i = j, i = k, j = k and triplets that only differ in a permutation of i and j are considered identical. Then, the following steps are executed. If the above steps do not lead to a valid dipole subtraction term, no divergence can be present. The real emission, a conventional tree-level process, is grouped together with all its subtraction terms into one computational unit and their respective cross sections and observable values, O( m+1 ) and all O( m i ), are treated as correlated. Integrated subtraction terms. Similar to the above discussed real emission corrections, the virtual correction configurations {ab} → {1, . . . , m} at order O(α n s α m ) in general exhibits poles due to both QCD and QED origins. To subtract them, both QCD and QED integrated subtraction need to be included. In fact, they naturally arise as counterparts to differential subtraction terms constructed for the corresponding real emission correction, as guaranteed by Bloch and Nordsieck [57] or Kinoshita-Lee-Nauenberg [58,59] theorem. Consequently, QCD and QED I, K and P operators are constructed. While their QCD version are described in detail in [49,50], their QED version of Eq. (2.6) are discussed below.
The I operator, on the one hand side, has the same initial state flavours and momentum fractions as the virtual correction. The K and P operators on the other hand, resulting from the combination of the integrated subtraction terms and the collinear counterterms, involve a summation over possible initial state flavours and comprise the following general structure in their dependence on the additional x integration variable Therein, both h(x) and k(x) are regular functions in x, while the plus distribution of g(x) is defined as Hence, the potentially computationally intensive matrix elements in B have to be calculated twice for every phase space point in addition to the flavour summation. To remedy this, the original expression of Eq. (2.6) is reformulated. Dropping the dependence on the {α dip } parameters and explicitly stating the dependence of the underlying Born term on the initial state momentum fractions and parton densities, B ab = f a f b B ab , the K and P operators it can be recast to with f a (η) being the parton density of flavour a and momentum fraction η in the proton, otherwise absorbed in the Born term B. All PDFs are evaluated at the same scale μ F .
are the analytically computed divergence free parts of the integral of the functions under the plus distribution. Its divergence at x = 1 is cancelled numerically on the second last line. Effectively, this reformulation results in a redefinition of the PDF for incoming parton a of B ab . As the integrand of the remaining integral over x a is a simple one-dimensional function without a pronounced structure, its numerical evaluation is very stable and can be accomplished by a single point for each summand per phase space point d (4) m . The K and P operators for the second incoming parton b are recast similarly.
The thus transformed form of the K and P operators require only a single evaluation of the potentially costly matrix elements in B while retaining the number of compu-tations of PDFs needed, speeding up the computation considerably for involved processes. This allows to generate and evaluate the underlying Born matrix element for both the I operators and the K and P operators at the same time. In fact, due to these three operators being simple multiplicative scalars, the common underlying Born matrix element is identical to the standard Born matrix element, allowing for their simultaneous calculation at no extra cost. The operators themselves are then built from dipoles constructed from all doublets {i, k}, i = k, of external partons available in the partonic process for which the charge-correlator Q 2 ik does not vanish.

External photons
External photons can play different roles in a calculation: they can either be resolved or unresolved. According to this distinction they should be treated differently at NLO EW [31,60]. Initial state photons are always unresolved at a hadron collider. They thus should be treated in a short-distance scheme, allowing them to split into massless fermions, necessitating a proper subtraction of infrared divergences. Final state photons on the other hand can play both roles. If they are considered resolved, they should be treated in an on-shell scheme and no explicit photon splitting is allowed. Concerning the dipole subtraction discussed in this paper they are thus neutral particles and do not form part of a dipole, except as possible recoil partner of another unresolved photon. A final state unresolved photon, on the other hand, again must be treated in a short-distance scheme, necessitating their splittings to be subtracted.
As discussed in Sect. 2.1, the dipole picture is not necessary to capture the divergences of photons splitting in pairs of massless fermions, leptons and quarks, due to the absence of soft divergences that necessitate the correlation with the emissions off other partons of the event. Nonetheless, it offers a possibility to assign one or more recoil partners to absorb the transverse momentum of the splitting, thus fitting these purely collinear splittings into the dipole picture. The choice of spectator is essentially arbitrary, and all other partons of the event offer being viable spectators. Only the following condition has to hold for every splitting photon ıj , cf. The κ ıjk are set to the phase space point independent value −n −1 spec , with n spec the number of assigned spectators, thus trivially fulfilling Eq. (3.4).

Process mappings
In general, physical cross sections include multiple different partonic channels. However, many of these partonic channels share identical squared matrix elements, potentially differing by constant factors. They, thus, do not need to be recomputed for every flavour channel but can be reused. As a typical example, consider the production of a lepton pair in association with two jets. The process g d → e + e − g d shares a common squared matrix element with g s → e + e − g s and g b → e + e − g b on Born level at O(α 2 s α 2 ). Hence, the squared matrix elements of the latter two partonic processes are mapped on the first, reusing its computed value. In this way, of the 95 partonic channels of this process, only 30 have to be computed. Further mappings of individual graphs and subgraphs are implemented but are not discussed in the following, see [47,48].
In SHERPA's matrix element generators various forms of process mappings are implemented to reduce the computational complexity and memory footprint, both for the virtual and real emission corrections, cf. [47][48][49][50]. However, while for NLO QCD corrections to the leading order Born process it is true that if the Born process is mappable onto another existing process, then also both the virtual correction and the insertion-operator-augmented colour-correlated underlying Born process of the integrated subtraction term are mappable to the same process. This is no longer true when considering the NLO EW or NLO QCD corrections to subleading Born processes.
In the present implementation of a full NLO QCD and NLO EW subtraction in the matrix element generator AMEGIC the following process mapping strategy is followed for the real subtracted contributions: • The real emission process and its associated dipole subtraction terms are grouped in one computational unit. A given partonic channel of the real emission process can be mapped onto another already existing one if both processes consist of the same diagrams and all involved (internal and external) particles have the same masses and widths and the same underlying interaction (coupling factors may differ by a constant). Is this the case the whole computational unit can be mapped and the result of the mapped-to process can simply be reused.
• Individual dipoles can be mapped if the emitter, emittee and spectator indices are identical, and the underlying Born process can be mapped according to the above rules. In this case, the result of the mapped-to process can simply be reused. • Underlying Born matrix elements can be mapped if the Born-level emitter ıj carries the same indices and the underlying Born process itself can be mapped. The spin correlation insertion operator, needed if parton ıj is either a gluon or a photon and described in [49], is encoded in the calculational routines and necessitates the above restriction. Can the underlying Born process be mapped, only the calculational routines can be shared, reducing the memory footprint, but due to the potentially differing underlying Born momenta its result has to be recomputed.
The strategy for the virtual subtracted contributions reads as follows: • The virtual correction process, interfaced from an external virtual correction provider, as well as its associated integrated dipole subtraction terms and the collinear counterterm, taken together and reformulated according to Sect. 3.1, are grouped into one computational unit. If the underlying Born processes of the integrated subtraction contribution (in all orders required) can be mapped according to the above rules and the virtual correction provider confirms that the virtual correction can be mapped onto the same process that SHERPA's tree-level matrix element generator maps the underlying Born processes onto, the whole computational unit is mapped. In this case, the result of the mapped-to process can simply be reused. The K and P operators, whose internal PDF factors depend on the initial state flavour, still have to be recomputed, but their matrix element coefficients are cached. • If the virtual correction provider cannot confirm the mapping performed for the underlying Born process, only the underlying Born process is mapped and virtual correction is recomputed. Again, the K and P have to be recomputed, but their matrix element coefficients are cached.
Here, efficiency is lost if the virtual correction provider uses less efficient process mappings.

Fiducial phase space definition
Phase space restrictions are an essential part of the implementation of a framework for automated NLO calculcations. These cuts, however, need to be applied in an infrared safe way. At NLO, they must not discriminate between a massless parton before and after its collinear splitting or before and after a soft gluon or photon emission. Thus, if QCD singularities are present, massless QCD charged particles must be clustered into jets before any further cuts are applied. Similarly, in case of the presence of QED divergences, massless charged particles must either be dressed with the surrounding photons or be included in the jet algorithm. Subtleties arise in the presence of both QCD and QED singularities simultaneously. Here, usually, only a fully democratic jet finding can consistently treat all singularities, although specialised solutions exist for simplified situations. Massive QCD and QED charged particles may be treated as bare as their mass shields the collinear singularity, but can also be included into jet finding and dressing algorithms. An intermediate scheme which includes only the logarithms of the parton mass needed to regulate the collinear divergences [45], but are otherwise treated massless in the calculation, is not implemented.
To this end, the implementation in SHERPA is equipped with a range of algorithms to define infrared safe quantities on which further restrictions can be applied. Multiple such selectors can be nested.
DressedParticleSelector This selector takes a choice of dressing algorithm (cone or sequential recombination) and (flavour-dependent) dressing parameters (cone radius or radial parameter and exponent). All charged particles of the process are then dressed with all photons using the specified algorithm with the given (flavour-dependent) parameters. Their four momenta are added such that four momentum is conserved. The dressed charged particles may no longer be on-shell and the photons used to dress the charged particles are removed from the list of particles. The resulting list of particles and their momenta are then passed to all subselectors. Jet_Selector This selector uses FASTJET [61] to build jets from a given list of input particles. It takes a list of flavours that are considered as jet finding input particles, the jet finding algorithm and its parameters including phase space boundaries in p ⊥ , η or y, as well as a minimal and maximal number of jets to be found as arguments. Additionally, clustered jets can be tagged or anti-tagged based on their flavour content, including relative and absolute constituent momentum requirements (e.g. b tagging a jet if one of its constituents is a b-quark or anti-γ -tagging a jet if its photon constitents carry in excess of z thr fraction of the total jet momentum). The clustered jets as well as all particles not used as jet finding input are passed on to all subselectors. Isolation_Selector This selector uses the smooth cone isolation of [62] to isolate particles of a given flavour against particles of another flavour. It takes both the isolation flavour and rejection flavour list as well as the algorithm parameters as input. It further can be specified how many isolated particles of the given flavour should minimally and maximally be found in a given phase space volume (bounded by p ⊥ and η or y ranges), e.g. exactly two isolated photons with p ⊥ > 30 GeV and |η| < 2.35. The list of found isolated flavours and all other particles except for those that should be isolated, but are not, are passed to all subselectors.
Once an infrared safe definition of particles is found, the standard single-or multi-particle selectors such as PT (implementing transverse momentum requirements on a given flavour), PTmis (missing transverse momentum build from all neutrinos), or DR (angular distance between the two given particles) can be used. Further, SHERPA is equipped with a user hook system, providing the possibility for users to implement arbitrary routines for phase space cuts and dynamically load them at run time, without the need to modify their SHERPA installation.
Identified particles should in principle be defined using fragmentation functions, denoted D j i (z) for finding parton i in parton j at momentum fraction z. As, however, all D i i (z) have a δ(1 − z) as leading term in on-shell renormalisation schemes (and only differ by ratios of couplings in other schemes) simplified schemes exist that are applicable to many practical situations. Hence, no fragmentation function is implemented yet.

Flavour scheme conversion
All publicly available PDF sets containing QED effects in their evolution are fitted with five (MRSTqed [63], CT14qed [64], LUXqed [65,66], HKR16 [67], NNPDF3.0qed [68], NNPDF3.1luxqed [69]) or six (NNPDF2.3qed [70]) light flavours. Thus, for next-to-leading order calculations with four or less light flavours the following scheme-conversion terms need to be added for consistency [71] and NLO QCD/EW is the expectation value of an arbitrary observable O computed at NLO QCD or NLO EW, respectively, with n-flavour scheme parton densities and the corresponding strong coupling summed over all initial state contributions. B (n) ab is the Born term, as defined in Sect. 2, in the ab channel. The sums run over all (n f − n lf ) flavours of mass m i that are part of the n f -flavour PDF parametrisation but not the n lf -flavour scheme of the calculation and all combinations of partonic channels ab occurring in the n lf -flavour scheme, respectively. p is the power of the strong coupling in the Born process and gg ab = δ ga + δ gb , i.e. it takes the values 2 in the gg channel, 1 in all qg andqg channels, and zero otherwise. γ γ ab is defined analogously. Each logarithm of course only contributes if the scale is larger than the respective mass. As the electroweak coupling is not taken running in the common renormalisation schemes it is independent of the number of massless flavours in the calculation. For MSlike renormalisation schemes a similar term proportional to its power in the Born process is to be added.

Checks of the implementation
The implementation of the formalism described in the previous sections in SHERPA needs to be validated. To this end, both its independence of its internal free parameters, {α dip } = {α FF , α FI , α IF , α II }, κ and the choice of spectator in photon splittings, and its agreement with independent implementations for fixed values of these parameters need to be tested. While the latter were carried out in [19,21,25,31,32] against the private implementations in Munich [74], in [51] against MadGraph5 [75], and in [51,53] against Recola [76] and found full agreement, the former represent a powerful check of internal consistency. Further, as the {α dip } parameters regulate the phase space coverage of the differential subtraction terms, they can be used to lower the average number of contributing dipoles for a given real emission contribution and, thus, reduce the computational costs of the real-subtracted contribution.
There are now two independent aspects of the calculation that need to be checked: of spectator for photon splittings. It thus needs to be verified that the sum of the contributions of the 0 coefficient of the integrated subtraction terms and the differential subtraction terms is independent of these parameters and choices. In order to arrive at a finite result in d = 4, cf. Eq. (2.2), the real emission correction and the collinear counterterms are added. The corresponding quantity is defined as and will be evaluated for processes containing all available dipole configurations in the following.
In the following, the inclusive or fiducial cross section contribution σ IRD is evaluated for a range of different processes, testing all possible dipole configurations. Throughout the input parameters of Table 1 are used in the G μ scheme, although several other EW input parameter schemes are available in general. Similarly, all unstable particles are treated in the complex-mass scheme [77]. Further, the CT14nlo [78] and CT14qed [64] 1 PDF sets with five active and massless flavours and their corresponding α s parametrisations with α s (m Z ) = 0.118 are used. While the use of the CT14nlo PDF set for NLO EW calculations is, strictly speaking, inconsistent due its missing QED evolution of the initial state quarks, it offers the quantification of the importance of photon initiated processes in these technical comparisons. To the same end, the CT14qed PDF set is not employed using its the best fit value for the intrinsic inelastic photon momentum fraction at the reference scale of Q = 1.295 GeV p γ 0 = 0.05%, but rather evaluate it for the extremes of its 1σ uncertainty, p γ 0 = 0 and p γ 0 = 0.14%, where applicable. In processes where jets need to be constructed to define a fiducial phase space volume, anti-k t jets with R = 0.4 and p ⊥ > 30 GeV are used [80] and both partons and photons are considered as constituents. Similarly, if leptons need to be defined for the same purpose, they are dressed with all photons in a cone of R = 0.1. The number of phase space points in the computation of σ IRD ({α dip }, κ, c γ k ) is kept constant for each process, such that the indicated statistical Table 1 Numerical values of all input parameters. The gauge boson masses are taken from [72], while their widths are obtained from stateof the art calculations. The Higgs mass and width are taken from [73]. The top quark mass is taken from [72] while its width has been calculated at NLO QCD. In calculations where a massive particle is present as an external state, its width is set to zero  13 TeV, σ IRD is computed for both the production of at least two jets at a hypothetical ν e −ν e collider at a centre-of-mass energy of 1 TeV and inclusive single jet production in equally hypothetical ν e p deep inelastic scattering (DIS) with the same centre-of-mass energy are calculated. Figure 2 now details σ IRD for all three setups. ν eνe → j j production, detailed in the left hand side plot, comprises only FF dipoles and, thus, only depends on α FF . Varying its value over four orders of magnitude leads leaves the value of σ IRD unchanged within the statistical uncertainties. The black line is placed at the central value of the computation with the smallest statistical uncertainty to guide the eye. As is evident, lowering the α parameter too much, i.e. restricting the subtraction to act only on configurations very close to the divergence, results in large cancellations between the real-subtracted and the integrated dipole contributions of Eq. (4.1), degrading the statistical prowess of the calculation. Similarly, pp → ν en u e production, detailed on the right hand side, comprises only II dipoles and, thus, depends on α II only. Showing the results for all three choices of photon content, contributing directly through the respective γ q/γq channels as well as indirectly through the impact on the quark PDFs and, through momentum conservation, on the gluon PDF, a similar picture as for ν eνe → j j emerges. The three coloured lines now indicate the central value of the calculation with the smallest statistical uncertainty for each PDF choice. In each case, σ IRD is found to be stable under the variation of α II . Finally, the centre plot shows the correction Fig. 2 {α dip }-dependence of the sum of integrated subtraction term and differentially subtracted real emission for ν eνe → j j, ν e p → ν e j and pp → ν eνe contribution for the hypothetical DIS process ν e p → ν e j requiring at least one jet in the lab frame. As FI and IF dipoles always occur in pairs the α FI and α IF dependence is evaluated together here. The resulting σ IRD (α IF = α FI ) are also found to be stable when varying both parameters simultaneously over four orders of magnitude.
Moving away from the simplest configurations, Fig. 3 displays the results for electron-positron pair production. The left hand side plot again displays their production at the hypothetical ν e -ν e collider used before, finding very similar results and their independence of α FF . The centre plot now, however, displays the production of an electron-positron pair at the LHC. At leading order, this process proceeds through qq → e − e + and γ γ → e − e + at O(α 2 ). Consequently, the qq channel exhibits six dipoles of all four types. In the γ γ channel, the number and types of dipoles present depends on the choice of possible photon splitting spectators c γ k . To regulate all LO singularities the fiducial phase space is defined by requiring the dressed electrons to have a transverse momentum of at least 20 GeV and the pair to have an invariant mass of at least 60 GeV. As σ IRD now potentially depends on the full set {α dip } no continuous variation over four orders of magnitude is performed. Instead, each of the four parameters is varied independently to 0.001 keeping all others at their default value of 1. These four variations are completed by setting α FF = α FI = α IF = α II = 1 and 0.001. The resulting correction contributions are found to be independent of {α dip }.
Finally, the right hand side plot of Fig. 3 displays the dependence of σ IRD on the choice of spectators in photon splittings. Only the γ q/γq channel is considered as both the γ γ and qq channels are independent of this choice in this process. The γ q/γq channel, however, still receives contributions from photon radiation off quarks in addition to the sought after photon splittings into quark-antiquark pairs. Hence, the minimum invariant mass is raised to 2 TeV to increase the relative importance of the photon PDF, enhanc-ing the photon splitting contribution. The resulting correction contribution at O(α 3 ) is found to be independent of all five choices of photon splitting spectators available. Please note that for this process scheme 0 and 3 as well as scheme 1 and 2 lead to the same set of allowed spectators, respectively, and therefore to identical results. Massive dipoles. Massive particles are so far only allowed in the final state. 2 Dipoles involving massive partons, either as emitter or spectator, comprise only three types: FF, FI and IF. Emittees are always considered massless, otherwise no singularity would be present. In Fig. 4 again top-anti-top pair production at a hypothetical ν e -ν e collider and single top production at a hypothetical ν ep collider is considered in order to study the individual dipoles separately. In the left plot, the α FF (in)dependence of the O(α 3 ) correction contribution σ IRD , containing only massive FF dipoles, is shown. It exhibits the familiar picture of decreasing statistical prowess of the calculation with too small α FF , but otherwise consistent results. The right plot details the α IF and α FI dependent massive dipoles in the hypothetical DIS scenario. As before, α IF and α FI are varied simultaneously for this purpose and the independence of the corrections contribution on both parameters is observed.
The left plot of Fig. 5 now investigates the dependence of σ IRD on the κ parameter. In only arises in FF dipoles of gluons splitting into gluons or massless quarks or photons splitting into massless fermions in the presence of a massive spectator. Consequently, to restrict the number of additional contributions, top-anti-top pair production in association with one jet at the hypothetical ν e -ν e collider is con-     Fig. 6 {α dip }-dependence of the sum of integrated subtraction term and differentially subtracted real emission for ν eνe → tt and pp → tt neither gluon nor photon splittings. Due to the relative size of g → gg and g → qq splittings in relation to gluon radiation off the top quarks the κ dependence of the O(α 2 s α 2 ) contribution is much more pronounced than at O(α 4 ), where photon radiation off the top quarks overwhelms the photon splitting contribution. Nonetheless, at both orders an independence of the σ IRD of κ is found. In addition, the right plot investigates the influence on the choice of spectators for the final state photon splitting ocurring at O(α 4 ). No dependence on this choice is observed. Please note that for this process scheme 0 and 3 as well as scheme 1 and 2 lead to the same set of allowed spectators, respectively, and therefore to identical results. And indeed, the correction contribution σ IRD is found independent of {α dip } for each such order. External W bosons. Lastly, it may become necessary to also consider external W bosons (or other massive charged particle with spin > 1 2 in BSM theories) as stable final state particles, e.g. to reduce the computational complexity for high final state multiplicity processes where off-shell effects and effects in the decays can be ignored or recovered through other means [83]. In this case the literature does not provide expressions for the respective massive dipole functions. As their mass, however, can be assumed to be large enough to suppress collinear radiation well enough, this only leaves the spin-independent soft photon emission limit. 3 Here, both the expressions for the radiation of a photon off a massive scalar or a massive fermion can be used.  Figure 7 details the production of a W + W − pair in the hypothetical ν e −ν e collider, separating the FF dipoles and their α FF dependence, on the left hand side. As before, σ IRD is found to be independent within the statistical accuracy and also independent of the whether the massive scalar or massive fermion subtraction terms are used. The right hand side focusses on their production at the LHC. Again, all dipoles contribute at O(α 3 ), leading again to the afore described sixpoint variation. Also in this case, the result is independent of {α dip } and the choice of scalar or fermionic subtraction term. The FI and IF dipoles cannot be investigated separately, as in the ν e p → ν e t case, due to charge conservation.

Conclusions
This paper detailed the construction and implementation of the adaptation of the Catani-Seymour subtraction formalism for NLO EW calculations. Besides the translation of the QCD dipole functions to the QED case, several other issues have been addressed. They include the special role photon splittings play in the formalism, embedding extermal massive emitters of spin > 1 2 into the formalism and the interplay of QCD and QED subtractions for processes exhibiting both kinds of divergences. The resulting general subtraction for NLO EW calculations has been implemented in the SHERPA Monte-Carlo event generator framework. Interfaces to OPENLOOPS, GOSAM and RECOLA to access the needed virtual corrections exist and are fully functional.
In addition to the checks against independent implementations on the level of partial and total cross sections performed in previous publications, numerous internal cross checks for independence of technical parameter choices, {α dip } = {α FF , α FI , α IF , α II }, κ and the choice of spectator in photon splittings c γ k , have been presented here. This imple-mentation will become publically available in the near future with the next major SHERPA release and an extension to the COMIX matrix element generator [84] is foreseen.

A Differential splitting functions
This appendix details the complete functional form of the dipoles introduced in Eq. (2.3), D i j,k , D a i j , D a j,k and D a,b j . The functional forms of the spin-dependent insertions V i j,k , V a i j , V a j,k and V a,b j and the kinematic maps of the differential splitting functions are taken from the original QCD case [40,41]. They are summarised in the following.

A.1 Final-final dipoles
Dipoles with both emitter and spectator in the final state take the form Therein, the charge-correlator is defined in Eq. (2.5). All momenta of the dipole are on-shell and the total four momentum flowing through it is given by It is thus invariant under the emission. The momenta of the parton before and after the splitting are connected through the map wherein mk = m k and the Kallen function is defined as λ (a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc. The splitting variables read y i j,k = p i p j p i p j + p i p k + p j p k , which are defined on the intervals y i j,k ∈ [y − , y + ] and z i ∈ [z − , z + ]. The boundaries are given by with q 2 i j,k = q 2 −m 2 i −m 2 j −m 2 k . As can be seen, a divergence, residing at y i j,k = 0, is only present if either i or j are massless. The relative velocities between ıj andk, i + j and k, and i + j and i are given by They are introduced to facilitate the analytic integration and only take effect away from the singular limit.  The charge-correlator is defined in Eq. (2.5). All momenta of the dipole are on-shell .22) and the total four momentum flowing through it is given by 23) before and after the emission, respectively. Although q acquires transverse momentum after the emission, q 2 =q 2 holds, i.e. the mass of the dipole is invariant. The momenta of the partos before and after the splitting are connected through the map While the initial state cannot absorb recoil transverse to the beam access, it is transferred to the final state through a collective boost of all final state partons. The splitting variables read As the splitting partons are all in the initial state and therefore massless, the singularity at v j = 0 is present in any case. Adopting the above convention for labeling the spin of the emitter aj the dipole functions are defined as As only massless particles are considered as initial states, no splitting functions involving massive scalars are needed.

B Integrated splitting functions
This appendix summarises the complete functional form of the integrated dipole terms cast in the form of the I, K and P operators [40,41].

B.1 Terms of the I perator
Each dipole contribution to the I operator reads cf. Eq. (2.8). Therein, V ik takes the general form ik contains all infrared singularities (double poles for massless fermion splittings, single pole for massive splittings), and a non-singular piece V (NS) ik that incorporates all finite dependences on the masses of the emitter and the spectator. The latter vanishes in the massless case, see below.
In the case of photon splittings, only the non-singular part survives, again being non-zero only in the case of massive spectators. The singular part vanishes in case of photon splittings.
The singular term V (S) ik is symmetric in i and k, it is given by and the velocity factors parametrise the mass dependence. Therein, q 2 ik = ( p i + p k ) 2 = s ik + m 2 i + m 2 k is the dipole invariant mass. In d = 4 − 2 dimensions, > 0, the massless limit is approached smoothly. Subsequently taking the limit → 0, however, results in −2 poles only if at least one of both particles making up the dipole is massless. Completely massive dipoles only result in single infrared poles.
The non-singular V (NS) ik contributions encode the additional finite m i/k > 0 contributions. They vanish in the limit that both m i → 0 and m k → 0 and their precise form can be found in Eq. (6.21)-(6.26) of [41]. Of particular interest is the term for a photon emitter in the presence of a massive spectator, with |q ik | = q 2 ik . Thus, choosing κ = 2 3 somewhat simplifies the integrated subtraction term. Further, In case of the photon, the sum runs over all massless charged flavours present in the model. The flavour constants γ i and K i are given by and are independent of the particle masses. For the photon, the sums run over all charged massless flavours of the model, N C, f being their respective colour multiplicity. Lastly, the dependence of the I operator on the technical

B.2 Terms of the K and P operators
The K operator of Eq. (2.9) reads (B.9) Therein, the K collect universal terms present in all splitting involving an initial state as either emitter or spectator. They are defined as As can be seen, they are independent of the final state and thus independent of any final state particle masses. Here, the regular part of the Altarelli-Parisi splitting functions are given by 11) with N C,a the multiplicity in the QCD representation of parton a (3 for quarks, 1 for leptons). Similarly, the -dependent part of the splitting functions is defined aŝ Thus, for a = a = γ the K term takes the simple form and comprises an end-point term only. The factorisation scheme dependent terms K FS vanish in the MS scheme, and the K i,aa , containing remnants from final state splittings where the initial state a forms the spectator. As it also depends on the final state flavour i, its forms are given as follows. The contributions, listed separately for scalars, fermions and photons, read and The missing function J a γ Q (x, y) is given in Eq. (5.58) of [41]. In the massless limit this reduces to Again, for a = a = γ the end-point contribution cannot arise as no corresponding splitting is integrated over, K t γ γ,k = 0. Finally, theK terms, containing correlations between both initial state partons, read

(B.19)
They thus vanish for a = a = γ . Finally, the P operator is given by The explicit dependence of the I operator on the {α dip } parameters, or more precisely only on α FF , is given by [54,55,[86][87][88] if i is a massless fermion − 2 ln 2 1 + y + − x 1 + y + + 4 ln 1 + y + 2 ln 1 − y + + x 1 − y + + ln 1 + y + 2y + ln 1 − y 2 + + 2x + y + 1 − y 2 + + Li 2 1 − y + 1 + y + − Li 2 1 − y + 2 or if i a massive fermion or if i is a photon or a massive scalar a + x a ln a(a + x)(a + x + ) 2 m k = 0 ( C . 4 ) As can be seen, the massless case can be recovered from the massive one in the limit that m k → 0. The A I ik term could in principle also be shifted into the K operator to leave the I operator invariant. To be applicable to pure final-state calculations, though, it is left where it is. The various abbreviations above are defined as The dependence of the K operator on the {α dip } parameters, or more precisely only on α FI , α IF and α II , is given in [54,55,86,87,89]. It has the general form The superscripts indicate the general functional form of the contribution. In case of the α FI and α IF depedence the existing plus distributions are modified in the K t terms as follows (C.8) As the singular point sits at x = 1 this modification amouts to a finite contribution. The α IF on the other hand leads to finite regular functional contribution in the A K t ab term given in [87,89].