Automation of NLO QCD and EW corrections with Sherpa and Recola

This publication presents the combination of the one-loop matrix-element generator Recola with the multipurpose Monte Carlo program Sherpa. Since both programs are highly automated, the resulting Sherpa+Recola framework allows for the computation of -in principle- any Standard Model process at both NLO QCD and EW accuracy. To illustrate this, three representative LHC processes have been computed at NLO QCD and EW: vector-boson production in association with jets, off-shell Z-boson pair production, and the production of a top-quark pair in association with a Higgs boson. In addition to fixed-order computations, when considering QCD corrections, all functionalities of Sherpa, i.e. particle decays, QCD parton showers, hadronisation, underlying events, etc. can be used in combination with Recola. This is demonstrated by the merging and matching of one-loop QCD matrix elements for Drell-Yan production in association with jets to the parton shower. The implementation is fully automatised, thus making it a perfect tool for both experimentalists and theorists who want to use state-of-the-art predictions at NLO accuracy.

the parton shower. Section 4 then begins with the validation of the NLO EW computations of these benchmark processes before presenting combined predictions for NLO QCD and EW corrections to all of the process classes considered. After the conclusions in Section 5, details on the installation procedures and run-card commands are given in the Appendix.

Details of the implementation
This section presents some details of the techniques and algorithms used in both Sherpa and Recola, as well as information on the interface between them. For a detailed description of the methods employed by these programs, the reader is invited to explore the references given below.

The Sherpa framework
Sherpa [12,15] is a multipurpose event generator, which aims to simulate the entirety of exclusive scattering events in high-energy particle collisions. In order to achieve this, Sherpa contains several modules and algorithms to cope with the many different physics challenges of collider physics. These include methods for the generation and integration of hard-scattering matrix elements, QCD parton-shower simulations, and the modelling of the parton-to-hadron fragmentation process and the underlying-event activity.
Leading-order matrix elements are provided by the built-in generators AMEGIC [33] and COMIX [34]. For virtual matrix elements contributing at one-loop order, Sherpa relies on interfaces to dedicated tools. In particular, for the evaluation of QCD corrections, Sherpa has interfaces to BLACKHAT [2], GOSAM [7], NJET [6], OPENLOOPS [3] as well as the BLHA interface [35]. Infrared divergences appearing in the QCD virtual and real-emission amplitudes are treated by the Catani-Seymour dipole-subtraction method [16,17] that has been automated in the Sherpa framework [36]. The default parton-shower algorithm in Sherpa is based on Catani-Seymour factorisation [37,38].
Sherpa can also correctly combine the matrix elements for multijet-production processes at both LO and NLO QCD with the parton shower. There are several techniques available on the market for this task. The method employed in Sherpa for the merging of varying parton-multiplicity tree-level processes, called MEPS@LO, is based on shower truncation and is presented in Ref. [39]. Its generalisation to NLO QCD matrix elements, dubbed MEPS@NLO, is discussed in Ref. [40]. The latter is based on the MC@NLO-style matching of NLO QCD matrix elements to the Catani-Seymour dipole shower [41]. In Ref. [42] the inclusion of finitemass effects in particular for bottom-quark-initiated processes has been discussed. To efficiently provide matrix-element and parton-shower computations with theoretical uncertainties, usually estimated by renormalisation and factorisation scale variations, or modified input parameters such as the strong coupling α s or the parton-density functions (PDFs), an event-wise reweighting approach has recently been implemented [43].
The present article contributes to the recent efforts to extend the Sherpa framework to the automated computation of NLO EW corrections, thus further improving the perturbative accuracy of the predictions. To this end, based on the results presented in Refs. [18,19], the implementation of the Catani-Seymour dipole-subtraction formalism has been extended to account for QED corrections [29,44]. Furthermore, the interfaces to the available one-loop amplitude generators had to be generalised to account for the variable orders in the strong and EW coupling parameters. In Refs. [29,30] this has been presented for the OPENLOOPS generator and applied to the production of on-and off-shell gauge-boson production in association with jets. Here, the corresponding developments for the Recola program are presented, along with the validation of the NLO QCD and EW corrections for a variety of phenomenologically important processes.

The matrix-element generator Recola
Recola [8,14] is a public matrix-element generator that is able to compute all tree and one-loop contributions to matrix elements squared in the Standard Model. For a detailed description of the code, the interested reader is referred to the manual of the program [8] as only the features relevant for the interface are discussed here. All amplitude computations are performed in the 't Hooft-Feynman gauge, and ultraviolet (UV) and infrared (IR) divergences are treated by dimensional regularisation. All tree and one-loop matrix elements squared are summed over spin and colour. In Recola everything is generated dynamically, and no process-specific libraries are needed to compute arbitrary processes at NLO QCD and EW accuracy in the Standard Model. Several schemes for the renormalisation of the electromagnetic coupling are available and are briefly discussed below. Furthermore, Recola allows for a flexible and generic treatment of the flavour scheme used for the running of the strong coupling constant as detailed in the following. An important feature of the code is its general implementation of the complex-mass scheme [45][46][47], which allows for the consistent computation of processes featuring resonant massive particles.
The computation of both the tree and one-loop amplitudes is performed in a recursive fashion [14]. For the tree-level amplitudes, the algorithm employed is inspired by the Dyson-Schwinger equation. At one-loop level, the amplitude can be written as the sum of tensor integrals where A CT and A R2 denote the counter terms and the rational terms [48], respectively. The former cancel the UV divergences present in the tensor integrals, the latter provide additional finite terms originating from dimensional regularisation. The core algorithm [14] used for the recursive computation of the tensor coefficients numerically is based on an idea of van Hameren [49]. To numerically evaluate the one-loop scalar [50][51][52][53] and tensor integrals [54][55][56], Recola relies on the Collier library [57,58]. Finally, we note that the implementation of the interface is valid for Recola-1.2 and subsequent releases.

The interface
The purpose of the interface is to provide arbitrary Standard Model one-loop matrix elements, generated with Recola, to Sherpa. This section outlines the choices made for the implementation of Recola in Sherpa in order to compute the virtual piece of NLO QCD and EW corrections. The other tools and ingredients needed for NLO computations, such as tree-level matrix elements or the subtraction method, are not addressed here. Indeed, the tree-level, integrated subtraction term and the real-subtracted parts are all computed by Sherpa without the help of Recola. The details of these implementations can be found in the related publications cited above. In Recola, the generation of the matrix elements is done on the fly, and once Recola is installed and linked to Sherpa, any processes can be computed readily. 3 In the initialisation phase of a Sherpa integration run, all necessary partonic one-loop processes are registered in Recola and automatically generated as soon as the actual integration starts. No extra processspecific libraries are needed. The combination Sherpa+Recola allows thus to compute-in principle-any tree-level-induced Standard Model process at NLO QCD and EW accuracy and any one-loop-induced process at LO. As for any computation in Sherpa, the processes, event selection, parameters, etc. are defined in a run card. The masses and widths issued in the run card are by default the on-shell masses, apart from the Z and W masses which are assumed to be the pole masses. Switching to on-shell masses also for the latter is possible through specific keywords (see Appendix B).
For the renormalisation of the electromagnetic coupling constant, three schemes are available [8]. In the G µ scheme, the electromagnetic coupling α is derived from the Fermi constant G µ . In the α(0) scheme, α is fixed from the value measured in Thomson scattering at p 2 = 0, while in the α(M Z ) scheme, α is renormalised at the Z pole and thus implicitly takes into account the running from p 2 = 0 to p 2 = M 2 Z . Due to the ambiguity in defining a real value of α in the complex-mass scheme, the numerical value of α computed by Sherpa is taken as an input by Recola to ensure compatibility.
Concerning the strong coupling, α s , the computation of the corresponding counter term is done consistently according to the PDF set used. The strong coupling constant is extracted from the PDF set as well as the flavour scheme and the quark masses used for the PDF evolution. These parameters are then used for the computation of the strong-coupling counter term δZ gs which reads [8] δZ gs = − α s Q 2 4π where N l is the number of light quark flavours and the index q runs over the heavy flavours. The parameter ∆ UV contains the poles in D − 4 as described in Ref. [8]. In variable-flavour schemes, all quarks lighter than the scale Q are considered light while the remaining ones are treated as heavy.
In fixed N f -flavour schemes, the N f lightest quarks are considered light while the others are treated as heavy. As emphasised above, the flavour schemes and quarks masses used to compute δZ gs are set consistently in Recola according to the PDF set. Nevertheless, specific commands described in Appendix B allow the user to choose all possible flavour schemes in combination with arbitrary quark masses. Finally, note that the quark masses used to compute δZ gs can, in principle, differ from the ones used in the rest of the matrix element. Ensuring consistent mass values between the run card and the PDF set used is left to the user. Finally, the IR and UV renormalisation scales are fixed to 100 GeV by default in the Sherpa +Recola interface. Thus in general, the virtual part calculated with Sherpa+Recola cannot be directly compared to the virtual part from a Sherpa+OPENLOOPS [29] calculation. However, the sum of the virtual and integrated subtraction-term contributions is independent of the regulators. It is also possible to set the IR scale in Recola equal to a fixed renormalisation scale, in which case a direct comparison of the virtual corrections of OPENLOOPS and Recola is possible (see Appendix B).

NLO QCD validation
This section is devoted to the validation of the Sherpa +Recola interface for NLO QCD computations. This is accomplished by direct comparisons to results obtained from public 4 and private codes as well as published work at three different levels. First, for a broad range of processes, squared matrix elements for individual phase-space points are compared with Sherpa +OPENLOOPS. Next, a comparison of cross sections and differential distributions with NLO QCD accuracy is presented where the results are compared with those obtained from public codes, private codes, and the literature. Finally, to illustrate the applicability of the Sherpa +Recola framework for NLO QCD calculations matched to parton-shower simulations, results for Drell-Yan lepton-pair production in association with jets at MEPS@NLO QCD accuracy are presented. Furthermore, information on the memory consumption and run times are provided.

Phase-space point comparison
As a first validation, the sum of the virtual and integrated dipole part of the NLO QCD corrections to a wide variety of processes, calculated with both Sherpa+OPENLOOPS and Sherpa +Recola, are compared at the level of individual phase-space points. Using identical set-ups, this provides a stringent test of the implementation of the two one-loop generators in Sherpa. Sherpa's Python interface [59] has been extended for this purpose. For each process, i.e. individual partonic channel, 1000 randomly chosen phase-space points are considered. These correspond to the parton momenta in proton-proton collisions with a hadronic centre-of-mass energy of 13 TeV. Both for the factorisation and renormalisation scales, we choose µ R = µ F = √ŝ . We employ the NNPDF-3.0 NNLO set [60,61], featuring α s (M Z ) = 0.118 with a variable number of active flavours up to N F = 5 and two-loop running of α s . Electroweak input parameters are defined in the G µ scheme. For all considered phase-space points all final-state particles have to pass the following set of cuts: where i, j are any particles apart from a photon. (version 3.6.1) in this comparison. In Table 1 the logarithmically averaged relative deviation, ∆ VI , of the sum of the virtual corrections and integrated dipoles between Sherpa+Recola and Sherpa+OPENLOOPS is presented for 62 partonic processes. The logarithmically averaged relative deviation ∆ VI is defined as where dσ i VI is the sum of virtual and integrated-dipole contributions at the phase-space point i, and N p is the number of phase-space points. In Table 2 the squared one-loop amplitudes for 13 loop-induced processes are compared. The logarithmically averaged relative deviation ∆ LI is defined similarly as in Eq. (4) with dσ i VI replaced by the loop-induced cross section dσ i LI . For all processes considered, good agreement is found between the results of Sherpa+Recola and those of the public Sherpa+OPENLOOPS. For most processes the average relative deviation lies between O(10 −9 ) and O(10 −12 ), corresponding to an agreement to 9-12 digits on average.
As one can expect, the agreement decreases for processes with higher final-state particle multiplicity as well as for the loop-induced processes. This originates from the increase in complexity with the number of external particles and the fact that one-loop amplitudes appear squared in loop-induced processes. In addition, the presence of external gauge bosons, in particular gluons, worsens the agreement due to the additional spin and colour degrees of freedom.
In addition to the one-loop results presented here, the squared tree-level matrix elements of Recola have further been compared against the ones provided by Sherpa, through the matrix-element generators AMEGIC [33] and COMIX [34]. For all the processes listed in Table 1 the logarithmically averaged relative differences per phase-space point are well below O(10 −11 ). The number of phase-space points with differences above O(10 −8 ) amounts to at most a few per cent in the worst channels and for the vast majority of processes all 1000 phase-space points have differences below O(10 −8 ).
Regarding the time spend on the evaluation of the one-loop matrix elements, the performances of the generators are overall very similar. We have compared Recola against OPENLOOPS in the default configuration with CutTools, as well as against OPENLOOPS in combination with the tensor reduction of the Collier library (albeit using version 1.0 instead of version 1.1 which is used by Recola). While OPENLOOPS with CutTools is slightly slower, especially for the loopinduced processes, the performance of both generators with the Collier library is comparable. Recola features a longer initialisation time as it does not rely on pre-compiled process libraries as OPENLOOPS does. However, this initialisation phase becomes negligible compared to the overall run time in realistic applications. Overall, no significant differences in the run times have been observed. In Section 3.3, the performance and memory usage for both one-loop providers for a full-fledged matrix-element plus parton-shower simulation is presented.

NLO QCD fixed-order calculations
Next, we present results for cross sections and differential distributions at NLO QCD accuracy. Using the Sherpa +Recola framework, four process classes are considered, namely off-shell Z-boson production in association with up to two additional jets (a full MEPS@NLO set-up of this process is presented in Section 3.3), off-shell Z-boson pair production, on-shell top-quark pair production in association with a Higgs boson, and Higgs-boson production in association with an on-shell Z boson.

Z-boson production in association with jets
Input parameters: Proton-proton collisions at a centre-of-mass energy of 13 TeV are considered, and the Rivet analysis package [64] is used to analyse events. As before, the NNPDF-3.0 NNLO PDF set with up to five active flavours α s (M Z ) = 0.118 and two-loop running is employed. EW input parameters are defined in the G µ scheme.
In the Drell-Yan plus jets processes, the QCD jets are reconstructed by means of the anti-k T jet algorithm [65,66] with radius parameter R = 0.4, p T,j > 25 GeV, and |η j | < 3.5. The Z boson is assumed to decay into an electron-positron pair with a dilepton invariant mass in the range  Table 1: Average relative deviations ∆ VI between Sherpa+Recola and Sherpa+OPENLOOPS for the sum of the NLO QCD virtual and integrated-dipole contributions evaluated for 1000 phase-space points. The average relative deviations are calculated by averaging the logarithms of the relative deviations between the two results for the considered phase-space points. process 3.541 · 10 −10 gg → HZg 3.013 · 10 −07 gg → HH 4.023 · 10 −10 gg → HHg 1.420 · 10 −09 gg → ZZ 5.576 · 10 −08 gg → W + W − 6.350 · 10 −08 Table 2: Average relative deviations ∆ LI between Sherpa+Recola and Sherpa+OPENLOOPS for the matrix element squared evaluated for 1000 phase-space points for QCD loop-induced processes. The average relative deviations are calculated by averaging the logarithms of the relative deviations between the two results for the considered phase-space points. 66 GeV < M e − e + < 116 GeV. The renormalisation and factorisation scales are event-wise set to M e − e + .

NLO QCD validation:
The Drell-Yan-plus-jets processes are considered at the orders O α 2 , O α s α 2 , and O α 2 s α 2 for the LO cross sections for 0, 1, and 2 jets, respectively. The corresponding NLO QCD cross sections contribute at the orders O α s α 2 , O α 2 s α 2 , and O α 3 s α 2 . The total cross sections calculated with Sherpa+Recola and Sherpa+OPENLOOPS presented in Table 3 turn out to be identical within the given accuracy as they have been obtained using the very same phase-space points and OPENLOOPS and Recola agree to more than 9 digits.
In Fig. 1, the resulting transverse-momentum distribution for the Drell-Yan pair is shown, considering the production processes with zero, one and two additional jets evaluated at NLO QCD accuracy. Evidently, the low-p T region is dominated by the pure Drell-Yan process, i.e. without additional final-state jets at the Born level, where the finite recoil originates solely from the real radiative corrections. On the other hand, the tail of the transverse-momentum distribution is dominated by Drell-Yan-plus-one-jet processes, where the Z transverse momentum results from the recoil against the Born-level jet and the real radiation.
Here only the results obtained with Sherpa+Recola are displayed. All results have been dσ/dp T,Z [fb/GeV] Figure 1: The transverse-momentum distribution of the reconstructed Z boson in Drell-Yan plus zero, one and two-jet production at NLO QCD.
cross-checked against Sherpa+OPENLOOPS, using the identical phase-space points, and no significant deviations have been observed. In fact, for each observable bin the relative difference of the two predictions is well below 10 −5 .

Z-boson pair production
Input parameters: The set-up employed is the one of Ref. [27] which is, for completeness, repeated in the following. The predictions are for the LHC operating at a centre-of-mass energy of √ s = 13 TeV. The NNPDF-2.3QED NLO PDF set [61] with a variable flavour-number scheme and α s (M Z ) = 0.118 has been used for all computations both at LO and NLO. Furthermore, a fixed factorisation and renormalisation scale at the Z-boson pole mass µ R = µ F = M Z is employed. The strong coupling α s is extracted from the PDF set at the renormalisation scale µ R , and the electromagnetic coupling α is calculated in the G µ scheme according to denoting the Fermi constant. The on-shell (OS) values for the masses and widths of the massive vector bosons read They are converted into pole masses and pole widths according to Ref. [67], Since the top quark and the Higgs boson do not appear as internal resonances, their widths are set equal to zero. The corresponding masses read  and NLO QCD with the Sherpa+Recola interface, compared against the independent private multi-channel Monte Carlo program that was employed for the computations in Ref. [27]. The difference is expressed in standard deviations.  The masses and widths of all other quarks and leptons have been neglected. As in Ref. [27], the following acceptance cuts are imposed on the charged leptons ± : The jets from real QCD radiation are treated fully inclusively.

NLO QCD validation:
The Sherpa+Recola interface has been cross-checked for this process against the independent private Monte Carlo program that had been used for the computations in Refs. [23,24,26,27]. The comparison for the total cross section is shown in Table 4. Agreement to four digits is found between the two independent calculations within the statistical uncertainty. Figure 2, shows a comparison at the level of differential distributions between the two Monte Carlo integrations, once for the distribution in the di-muon mass m µ + µ − (left plot) and once for the four-lepton invariant mass m 4 (right plot). The upper panels plot the absolute predictions of the two integrations on top of each other, the differences being almost invisible. The lower panels show the bin-wise ratio between the two calculations. The fluctuations in the four-lepton invariant mass are almost everywhere below 0.5%, from the far off-shell region over the on-shell production threshold at m 4 = 2M Z up to 1 TeV. A similar pattern is observed for the di-muon mass in the left plot with slightly larger fluctuations.

Higgs production in association with a top-quark pair
Input parameters: In Ref. [68], a comparison between MadGraph5_aMC@NLO and Sherpa +OPENLOOPS for various cross sections at LO, NLO QCD, and NLO EW has been presented for the process pp → ttH. In particular, five different cross sections have been reported. The first one is fully inclusive, and no event selections are applied. Two of them are computed when applying a cut on the transverse momentum of the three massive final-state particles at 200 GeV and 400 GeV, respectively. Furthermore, a cross section with a higher transverse-momentum cut (500 GeV) on the Higgs boson only is presented and finally the cross section obtained by excluding events with a top quark with absolute rapidity lower than 2.5. The computations are done for a centre-of-mass energy of √ s = 13 TeV at the LHC. The masses of the involved particles read and the corresponding widths are all set to zero. The bottom quark is considered massless and its PDF contribution is included. The NNPDF-2.3QED PDF set with a variable flavour-number scheme and α s (M Z ) = 0.118 [61,69,70] has been used. The renormalisation as well as the factorisation scales are set to a common scaleĤ T /2, defined aŝ where the index i runs over all the final-state particles. The considered LO production cross section is at the order O α 2 s α .
NLO QCD validation: The obtained LO cross sections listed in Table 5 show a generally good agreement. The corresponding NLO QCD predictions of order O α 3 s α are reproduced in Table 6. The agreement is reasonable but a definite statement is not possible as the predictions made by MadGraph5_aMC@NLO and Sherpa+OPENLOOPS do not have statistical errors. In addition we have cross-checked the fully inclusive set-up at the level of distributions against the public Sherpa+OPENLOOPS implementation, finding perfect agreement in each bin. The distributions for Sherpa+Recola are shown in the plots of Fig. 3. 5 6.840(7) · 10 −3 36.3(2) 37.5 36.9  Figure 3: Differential distributions at a centre-of-mass energy √ s = 13 TeV for pp → ttH at the LHC in a fully inclusive set-up: transverse momentum distribution of the top quark (left) and rapidity distribution of the top quark (right).

Higgs-boson production in association with a Z boson
As last example, Higgs-boson production in association with an on-shell Z boson is considered. This channel is particularly interesting, as it receives sizeable contributions from the loop-induced gg → HZ channel.
Input parameters: The set-up used is similar to the one of the Drell-Yan example described in Section 3.2.1. However, here we use the partonic centre-of-mass energy √ŝ as the renormalisation and factorisation scale. For the Higgs boson we assume the decay into a pair of bottom quarks, that is, however, fully factorised from the production process. The mass of the bottom quark is thereby taken to be m b = 4.8 GeV.

NLO QCD validation:
The LO contribution appears at order O α 2 for the quark-initiated channel, hence the NLO QCD cross section is of order O α s α 2 . The loop-induced gg → HZ process contributes at order O α 2 s α 2 but is enhanced by the gluon PDF. The comparison for the total cross sections calculated with Sherpa +Recola and Sherpa +OPENLOOPS can be found in Table 7. Again, perfect agreement is found, originating from the fact that both cross sections have been evaluated using the same phase-space points such that there is no statistical difference, and the deviations per point are below 10 −9 (cf. Section 3.1).
In Fig. 4 the transverse-momentum distribution of the reconstructed Higgs boson is presented, where the quark-and the gluon-initiated channels are shown separately. Though much smaller than the quark-induced contribution, the loop-induced gluon-initiated channel still contributes around 10% to the total cross section for Higgs transverse momenta in the range 100-200 GeV. This nicely illustrates the importance of precise predictions of this phenomenologically important Higgs-production process. Considering even larger values of p T,H the gluon contribution falls off more steeply, due to the different slope of the quark and gluon parton-density functions. Let us note that the results for the differential distributions have explicitly been checked against Sherpa+OPENLOOPS, and using identical phase-space points perfect agreement has been found.

Matching to parton shower
To illustrate and validate the use of virtual QCD matrix elements from Recola for full particlelevel simulations, a state-of-the-art QCD calculation is presented where NLO QCD matrix elements of varying final-state parton multiplicity get matched to the QCD parton shower of Sherpa [37] and subsequently hadronised. In particular, the Drell-Yan process is investigated, i.e. pp → γ * /Z * → e + e − /µ + µ − , in association with jets in the MEPS@NLO scheme [40]. To this end, NLO QCD matrix elements are considered for up to two additional jets and the tree-level contribution for three final-state partons. The phase-space slicing parameter of the merging scheme is set to Q cut = 20 GeV. The NNPDF-3.0 NNLO PDF set is employed with α s (M Z ) = 0.118 and five active flavours. For the electromagnetic coupling constant α, the α(M Z ) scheme is used with a numerical value of α(M Z ) = 0.007764. In the MEPS@NLO approach the renormalisation and factorisation scales are chosen dynamically, based on the determination of an event-wise 2 → 2 core process and an associated sequence of nodal splitting scales, obtained through a clustering of the matrix-element partons that effectively inverts the Sherpa parton shower [38,39,71]. To estimate the scale uncertainties, a 7-point scale variation is considered for µ R and µ F , computed by reweighting the central prediction [43]. Both scales are varied independently by factors of 1/2 and 2, thereby omitting the variations with ratios of 4 between the two scales. The corresponding uncertainty is then taken as the envelope of all variations considered.
The Each jet candidate needs to be separated from the reconstructed leptons by ∆R j ≥ 0.5. For the study presented here, the public Rivet implementations of the two analyses were used, allowing for a direct comparison to particle-level results. Further details on the selections can be found in the respective publications. Figures 5 and 6 provide examples for the comparison of particle-level MEPS@NLO simulations based on Sherpa +Recola against data. In the left panel of Fig. 5 the inclusive transverse-momentum distribution of the reconstructed Z bosons is compared to the ATLAS data from Ref. [73]. The inclusion of the Z + j and Z + jj matrix elements at NLO QCD accuracy provides a good description of events with sizeable recoil of the Drell-Yan pair. The partonshower component, on the other hand, dominates the low-p T,Z region. Notably, the theoretical scale uncertainties are of order ±20% and well overlap with the experimental uncertainty band, indicated in yellow. The right panel of Fig. 5 presents the comparison of the jet-multiplicity distribution against the data from CMS [72]. By construction of the MEPS@NLO algorithm, in the given set-up, the first two bins have NLO QCD accuracy, while the third has LO QCD precision. All higher jet multiplicities solely originate from the QCD parton-shower component. The central theoretical predictions agree well with the data though there is a tendency to overestimate higher jet counts. However, taking into account theoretical and experimental uncertainties, simulation and data agree well.
Finally, Fig. 6 presents the transverse-momentum distributions of the two leading jets. The set-up provides NLO QCD accuracy for both observables. The predictions are compared to the respective LHC data from CMS [72]. Very good agreement of data and theoretical prediction is found. For the latter the theoretical uncertainty estimates from variations around the central scale choice increase for larger jet transverse momenta, exceeding the ±20% range seen in the more inclusive boson transverse momentum distribution.
The results for the jet-associated Drell-Yan process based on Sherpa+Recola presented here confirm related comparisons of simulations based on NLO QCD matrix elements matched with QCD parton showers with data, see e.g. Refs. [74,75]. This nicely illustrates the usability of the Recola loop-amplitude generator in these highly non-trivial multi-scale calculations. To this end, the same set-up has been explicitly compared, but using the OPENLOOPS generator for the one-loop amplitudes. For all distributions considered, embracing many more than presented here, full agreement between Sherpa+Recola and Sherpa+OPENLOOPS has been observed. Besides, in this "real-life" application no significant run-time differences have been observed. However, when using Recola instead of OPENLOOPS the allocated memory increases by about 50% for the process set-up considered.

NLO EW validation and combined predictions
To illustrate the capabilities of the combination of Sherpa with Recola, three processes have been computed at both NLO QCD and EW accuracy. These comprise off-shell vector-boson production in association with jets, the production of two off-shell Z bosons, and the on-shell production of a top-quark pair in association with a Higgs boson. These three channels represent phenomenologically very important LHC processes. Furthermore, they are highly non-trivial, and each process features different technical challenges regarding their evaluation to NLO QCD and EW accuracy, thus demonstrating the generality of our implementation. Since up to now there is no public code that allows the computation of arbitrary processes at NLO EW accuracy, the number of processes that have been checked is smaller than for QCD.
For each process, a short introduction is given followed by the description of the calculational  set-up and the actual NLO EW validation. Finally, the cross sections and differential distributions for combined NLO QCD and EW accuracy are presented. The cross sections including NLO QCD or EW corrections read respectively. The additive combination of the two types of corrections is straight-forward, while a multiplicative combination can be defined as The difference between these two ways of combining NLO QCD and EW corrections provides an estimate of the missing higher orders resulting from mixed QCD-EW contributions. The NLO QCD×EW combination can be understood as an improved prediction when the typical scales of the QCD and EW corrections are well separated.

Vector-boson production in association with jets
As a first process, the production of a vector boson plus jets is considered. This process has already been computed for the LHC running at a centre-of-mass energy of 13 TeV for both onand off-shell vector bosons at NLO QCD and EW [22,29,30], and therefore allows a comparison to existing results. From the technical point of view, the process is a mixture of QCD and EW contributions. In particular, for vector-boson production in association with more than one jet at NLO EW, interferences appear between QCD and EW production channels. This makes vector-boson-plus-jets production a good testing ground of the interface, although it is also an important process in its own right.
Input parameters: The input parameters for this process class are taken from Ref. [76] to allow a tuned comparison. For completeness, these parameters as well as the analysis cuts are detailed here. The renormalisation and factorisation scale for off-shell W-boson production iŝ H T,W /2, defined viaĤ The scales used in off-shell Z-boson production,Ĥ T,Z /2, are similarly given bŷ For on-shell weak-boson production, a slightly different scale,Ĥ T /2, is used, defined viâ where V =W, Z for W+jets and Z+jets production, respectively. No scale variations have been considered for this validation, although a variation of a factor of 1/2 and 2 on the central scale was considered in the original article.
For on-shell vector-boson production, the bosons are decayed to leptons in a factorised approach, thereby preserving the spin correlations [59]. For W-boson production, the leptonic decay channels W ± → e ±( ν ) e and W ± → µ ±( ν ) µ are considered. Similarly, for Z-boson production, the allowed leptonic decay channels are Z → e + e − and Z → µ + µ − . For all processes in this section, the masses for the Z boson, W ± bosons, H boson and top quark read The Fermi constant is taken to be G µ = 1.1667 × 10 −5 GeV −2 , and the G µ scheme is used to consistently define the EW parameters. The NNPDF-2.3QED NLO PDF set with a variable flavour-number scheme, QED corrections and α s (M Z ) = 0.118 [61,69,70] has been used for both LO and NLO calculations. For on-shell vector-boson production, the widths of the external bosons are set to zero in general. An exception to this rule is in the QCD-EW interference term introduced in the EW real-subtracted contribution which includes electroweakly produced jets leading to a poorly converging phase-space integration. Because this enters only as an interference term, it does not give rise to a true resonance and maintaining a width of zero is theoretically acceptable. Following Ref. [29], a small, artificial width for the W and Z bosons is introduced in this case in order to control the phase-space integration. In this publication we use 0.3 GeV.
For the off-shell vector-boson production processes, physical values of the vector-boson widths are used, as well as for other unstable particles such as the Higgs boson and the top quark, Furthermore, the complex-mass scheme is employed for the unstable particles in this case, and a unit CKM matrix is assumed. Photons within a rapidity-azimuthal-angle distance of R γq/ = 0.1 from a quark or lepton are recombined with a simple cone-like algorithm with the closest charged particle. Jets are defined with the anti-k T algorithm using R = 0.4 and cut variable  Any jet with more than 50% of its energy originating from a photonic contribution is removed from the jet list. For the cross section validation for on-shell W + j production, two phase-space regions, other than the inclusive cross section subject to the cuts (20), are considered, defined by the additional cuts p T,W > 1 TeV, p T,j > 1 TeV.
Distributions have been analysed with Rivet following Ref. [30] and using the cuts shown in Table 8.
NLO EW validation: Next, we present the validation of the Sherpa +Recola interface against published cross sections obtained with Sherpa+OPENLOOPS [29,76,30]. First, on-shell W + j production at a centre-of-mass energy of 13 TeV at the LHC is considered. Table 9 shows the relative difference between total cross-section calculations, with different phase-space cuts, for the Sherpa+Recola interface against the published numbers from Sherpa+OPENLOOPS. The errors quoted on the Sherpa+Recola cross sections are statistical in origin. While for the Sherpa+OPENLOOPS results scale uncertainties were presented, the comparably negligible statistical uncertainties were not listed. However, for the validation of the codes, it is necessary to demonstrate close statistical agreement with the published numbers for the central scale choice. Taking the statistical errors of Sherpa+Recola as benchmark, good agreement between the two calculations for the NLO QCD and EW total cross sections is observed in all cases.
Besides the integrated cross sections, Ref. [29] also provides distributions, which can be used for a more qualitative validation over a large phase space. Figure 7 displays the p T distribution of the hardest jet in on-shell W + + j production. The left-hand side of Fig. 7 Figure 7: Differential distributions at 13 TeV for pp → W + j at the LHC, both at LO and NLO EW. The left-hand figure shows the leading-jet p T with minimal cuts and the right-hand plot the leading-jet p T with an additional cut ∆Φ(j, j) < 3π/4. The lower panels display the ratio of the NLO EW calculation to the corresponding LO result.
prediction, and the right-hand side the effect of a phase-space cut ∆Φ(j, j) < 3π/4. This effect is very large, in line with the findings in Ref. [29]. The Sudakov behaviour in the large-p T region is clearly recovered once this ∆Φ(j, j) cut is included, because it removes the contributions from dijet-like structures with a soft W boson emitted. At NLO EW, these types of contributions can be better viewed as a real EW correction to dijet production.
Secondly, Z+jets production is considered. Distributions have been published for off-shell Z+jets processes in Ref. [30]. Although this is not such a rigorous validation as the total cross section or phase-space point comparisons for other processes, it allows a qualitative assessment of the agreement between the calculations across a large phase space. From the wide range of distributions presented in Ref. [30], we select here the distributions in the transverse momenta of the leading lepton, p T, 1 , and leading jet, p T,j 1 , (according to p T ordering) (Fig. 8). The pp → + − jj process provides a particularly good test of the Sherpa+Recola interface, because it includes interference terms between EW and QCD produced jets contributing at NLO EW. This makes the process a lot more challenging than the + − j or νj final state. For both plots in Fig. 8, the behaviour across the full p T range is in agreement with the observations in Ref. [30].
Combined predictions: Table 10 presents combined predictions for NLO QCD and EW corrections to the integrated cross sections for Z/ + − + jets processes. These cross sections include all of the cuts from the analyses, and correspond to the distributions presented in this section. The on-shell calculation of Z + jets includes the branching ratios to e + e − /µ + µ − . The large NLO corrections are dominated, in all cases, by the NLO QCD contribution.
Beginning again with W + j production, distributions at a centre-of-mass energy of 13 TeV at the LHC are presented for both on-shell and off-shell W + j production. Figure 9 shows the transverse momentum p T,W of the reconstructed W boson for both additive and multiplicative combination of NLO QCD and EW effects. The prescription used to combine the NLO corrections clearly has a large effect in the high-p T region of the plot, which is indicative of large higher-order  2.34 · 10 2 21.6 2.31 · 10 2 19.7 pp → + − jj 5.53 · 10 1 6.93 · 10 1 0.04 6.88 · 10 1 0.03 Table 10: Integrated cross sections for pp → Z + jets for a centre-of-mass energy of √ s = 13 TeV calculated with Sherpa+Recola. The cross sections at LO as well as for the additive and multiplicative combinations of NLO QCD and EW corrections are given. The cross sections are expressed in pb while the relative corrections are given in percent.
corrections. The right-hand plot shows the effect of imposing a cut ∆Φ(j, j) < 3π/4 on the two jets in the NLO case. This removes the contribution from the production of two hard jets and a soft W boson. At NLO EW, such cuts are useful in order to differentiate processes which should more correctly be considered as an EW real correction to dijet production.
Similarly, Fig. 10 shows the distributions for on-shell W+j production. The decays to leptons are treated in a factorised approach, and the NLO QCD and NLO EW corrections are applied only to the on-shell W + j final state. The NLO QCD and NLO EW corrections show very similar behaviour, whether on-shell or off-shell W-boson production is considered.
Results are also presented for Z+jets production in Figs. 11 and 12, but this time considering both the 1-jet and 2-jet channels as well as the on-shell vs. off-shell effects. The off-shell process naturally includes the effects from γ * → + − interference. As was mentioned in the validation section for this process, the Z + 2 jets final state introduces interference terms between NLO σ/σ NLO QCD Figure 11: Differential distributions at 13 TeV for Z + j production at the LHC. The left-hand plot shows the p T of the Z boson for on-shell Z production and the right-hand plot the same observable for the reconstructed Z boson for off-shell Z production. The lower panels display the ratio of the distributions to the NLO QCD result.
QCD and NLO EW, which are taken into account automatically. These plots display the (reconstructed) p T of the Z boson. The distribution for + − jj shown in Fig. 12, was also presented in Ref. [30], and can therefore be additionally viewed as a further validation of the Sherpa +Recola interface. Figure 11 presents the Z-boson p T in Z + j production at a 13 TeV LHC. Here, there is a larger impact from the NLO EW corrections in the low-p T region for + − j production than for on-shell Z + j production. This effect is smaller for the + − jj final state shown in Fig. 12, where little difference is observed between on-shell and off-shell Z-boson production across the entire phase space. In all cases, the Sudakov behaviour in the large-p T region is clearly observed. Also, the NLO QCD + EW and NLO QCD × EW curves for Z + jj in Fig. 12 show a very good agreement with each other, indicating that the higher-order corrections in this case are a lot smaller than for the distributions of Z + j production, where a large difference is observed in the high-p T Sudakov region. It is reassuring that these differences are largely removed for both the on-shell and off-shell Z+jets processes once higher jet multiplicities are included. This implies that a merged NLO QCD and EW sample would give an accurate picture of both NLO QCD and NLO EW corrections to this process.

Z-boson pair production
As second process we consider the production of two off-shell Z bosons with subsequent decays into pairs of different-flavour charged leptons. At leading order (LO), this is a purely EW process implying that interferences in different orders of the strong and EW coupling first occur at NNLO. As a consequence, the NLO QCD corrections at O(α s ) and the NLO EW corrections at O(α) may be computed independently. The NLO QCD corrections are known [77][78][79][80], and the complete NLO EW computations have recently been published [24,27]. σ/σ NLO QCD Figure 12: Differential distributions at 13 TeV for Z + jj production at the LHC. The left-hand plot shows the p T of the reconstructed Z boson for on-shell Z production and the right-hand plot the same observable for off-shell Z production. The lower panels display the ratio of the distributions to the NLO QCD result.  Table 11: Total cross sections calculated for pp → µ + µ − e + e − at LO and NLO EW with the Sherpa +Recola interface, compared against the benchmark numbers from Ref. [27]. The difference is expressed in standard deviations.
Input parameters: We use the same set-up as for the QCD validation of Z-boson pair production in Section 3.2. In addition, real photons from QED bremsstrahlung and charged leptons are recombined to dressed leptons if their separation in the rapidity-azimuthal-angle plane fulfils ∆R γ < 0.2, following the prescription of Ref. [27]. In contrast to Ref. [27], we refrain from including photon-induced contributions.
NLO EW validation: In Table 11, a comparison of the NLO EW total cross section is shown once obtained from Sherpa+Recola and once by combining Recola with an independent private multi-channel Monte Carlo integrator, i.e. with the results from Ref. [27]. Perfect agreement is found between the two calculations within the statistical uncertainty at (sub-)permille level. In Fig. 13, a comparison of the two independent calculations at the level of NLO EW differential distributions via the di-muon mass m µ + µ − and the four-lepton invariant mass m 4 is presented. Like in the corresponding QCD validation in Section 3. agreement. This is a highly non-trivial check, since the benchmark calculation of Ref. [27] has been cross-checked internally by two independent calculations both at the level of the employed matrix elements and at the level of the phase-space integration. Furthermore, the results from Ref. [27] were generated in mass regularisation and slightly different conventions of the dipole subtraction terms.
Combined predictions: The combined predictions for the total cross section including the QCD corrections from Section 3.2.2 and the EW corrections from this section are stated in Tab. 12. In Fig. 14, different NLO predictions are presented for distributions in the di-muon and four-lepton invariant mass. In the four-lepton invariant-mass distribution, we observe the typical pattern of large negative EW corrections of around −20% in the high-energy regime at around 1 TeV. The radiative tail below the pair-production threshold at m 4 = 2M Z with corrections around +30% is due to the fact that resonant contributions are shifted to lower values by real photon radiation. Since the LO cross section is falling off steeply in this region, the photonic corrections become large. A similar radiative tail is observed also in the di-muon prediction that amounts to positive corrections of up to +60%. While the QCD corrections are positive over the whole range of both the di-muon and the four-lepton invariant-mass distribution at the  order of +30%, the EW corrections exhibit a non-trivial sign change at the Z-boson resonance m µ + µ − = M Z and at the pair-production threshold m 4 = 2M Z , respectively. Further discussion of this issue can be found in Ref. [27].

Higgs production in association with a top-quark pair
As last example, we consider the on-shell production of a pair of top quarks in association with a Higgs boson. The LO cross section is of order O α 2 s α . Hence the NLO QCD and EW corrections contribute at order O α 3 s α and O α 2 s α 2 , respectively. At NLO EW, this computation features QCD-EW interferences. In addition to computing the EW corrections to the QCD-mediated production of the top-quark pairs, one must also consider the QCD corrections to the interference of the QCD and electroweakly produced top-quark pairs. Moreover, the final state consists exclusively of massive particles which is not the case for any of the previously presented processes. This process constitutes thus a non-trivial validation of the implementation. Concerning on-shell top quarks, the process has already been computed at NLO QCD [81][82][83][84][85] and at NLO EW [86][87][88]. It has also been matched to a parton shower [89][90][91]. On the other hand, for off-shell top quarks, the NLO QCD [21] and EW [28] corrections have been computed only recently.
Input parameters: We use the set-up of Ref. [68] which has been described in Section 3.2.3. Concerning the electromagnetic coupling α, the α(M Z ) scheme is employed. Note that contributions originating from initial-state photons are neglected in order to match one of the set-ups of Ref. [68].
NLO EW validation: As for the QCD validation, we compare five NLO EW cross sections that have been computed by MadGraph5_aMC@NLO and Sherpa+OPENLOOPS in Ref. [68]. The obtained cross sections are reported in Table 13. Generally good agreement has been found with the results presented in Ref. [68]. Nonetheless, as no statistical errors are stated in the aforementioned reference, an exact comparison has not been possible.   Combined predictions: Next, combined NLO QCD and EW predictions for the inclusive cross section as well as for a few distributions are presented. No event selection is applied to the final state, meaning that we consider the fully-inclusive production process. The total cross sections at LO and NLO for an additive and multiplicative combination of NLO QCD and EW corrections as defined in Eqs. (13) and (14), respectively, are listed in Table 14. As the EW corrections are moderate, there are no big differences between the two combinations. This seems to indicate that the missing higher orders of mixed QCD-EW type are small in this case. In Fig. 15, the transverse momentum as well as the rapidity distribution of the Higgs boson are displayed. The effects of the NLO QCD corrections are dominant over the whole transversemomentum range and are typically of the order of 25%. The EW corrections vary from about 1.2% at zero transverse momentum to −8.9% at 700 GeV. This behaviour is characteristic for Sudakov logarithms that grow large when all invariants involved in the process become large. In Fig. 16, the distribution of the transverse momentum as well as the rapidity of the top quark are shown. Again, the QCD corrections are large over the whole transverse momentum range and amount to at 35% at low transverse momentum to go down to 22% at p T,t = 600 GeV. The relative EW corrections also decrease from about 2.6% to reach −6.7% at 600 GeV.

Conclusion
After the very successful Run I, culminating in the discovery of the Higgs boson, the LHC is now operating in the Run II phase. This phase at √ s = 13 TeV might ultimately lead to the  Figure 16: Differential distributions at a centre-of-mass energy √ s = 13 TeV for pp → ttH at the LHC: transverse-momentum distribution (left) and rapidity distribution of the top quark (right). The lower panels show the ratios σ LO /σ NLO QCD , σ NLO QCD+EW /σ NLO QCD and σ NLO QCD×EW /σ NLO QCD . discovery of New Physics or confirm, to an even higher degree of accuracy, the theory of the Standard Model of particle physics. In either case, precise (including at least NLO QCD and electroweak corrections) and appropriate (with event selections reproducing the experimental set-ups) predictions for a plethora of Standard Model processes are needed. Such state-of-theart predictions are, on the one hand, required for precision measurements in the Standard Model.
On the other hand, in New Physics searches, they are compulsory to obtain realistic estimates of the Standard Model expectations in order to discriminate possible New Physics contributions.
To provide such predictions, public Monte Carlo programs are the ideal tools as they can be used by both the experimental collaborations and the theory community. In this publication, the combination of the one-loop matrix-element generator Recola with the multipurpose Monte Carlo program Sherpa has been presented. In particular, a short presentation of both codes as well as the main features of the interface have been given. The Sherpa+Recola framework is designed for Standard Model predictions and offers the possibility to compute-in principle-any process at NLO QCD and electroweak accuracy. A large fraction of this article is devoted to the validation of the implementation of Recola in Sherpa. This entails comparisons of squared matrix elements for individual phase-space points, fixed-order cross sections, and differential distributions, as well as the merging/matching of NLO QCD matrix elements with Sherpa's QCD parton shower. These comparisons are performed against public and private codes as well as results presented in the literature.
Following this validation, predictions have been presented at NLO QCD and electroweak accuracy for three specific processes: both on-and off-shell vector-boson production in association with jets, off-shell Z-boson pair production, and on-shell production of a top-quark pair in association with a Higgs boson. In addition to their distinguished physical relevance, these processes constitute a good testing ground for this fully automatised implementation. They feature both massive and massless final states, as well as strongly and electroweakly interacting final-state objects. In addition to fixed-order computations, all other functionalities of Sherpa (the QCD parton shower, hadronisation etc.) can be used along with Recola. To demonstrate this, NLO QCD matrix elements for Drell-Yan production in association with multiple QCD jets have been merged and matched to the parton shower. For illustrative purposes, some resulting predictions have been compared to actual LHC data.
The Sherpa+Recola combination is readily publicly available for NLO QCD predictions. The required methods to perform NLO electroweak calculations on the Sherpa side will be made public soon. This ultimately makes it an ideal tool for both experimentalists and theorists to obtain NLO QCD and electroweak accurate predictions for Standard Model processes. It opens the possibility to perform systematic studies on the impact of electroweak corrections for a multitude of LHC production processes.

A Installation procedures
Recola-1.2 and subsequent versions are compatible with the interface to Sherpa described above. Recola in association with Collier can be downloaded from http://recola.hepforge.org. Once downloaded, the following command lines have to be issued: tar -zxvf recola-collier-1.2.tar.gz cd recola-collier-1.2 cd build cmake .. make Recola is then installed in association with Collier. Various compilation options can be found in the respective manuals [8,58]. We note that the Recola library has to be compiled dynamically (this is the default setting) to be used with Sherpa.
The Extra configuration options can be found in the manual of Sherpa available on the website. After this installation procedure, NLO computations can be readily performed.

B Specific run-card commands
As mentioned in Section 2.3, some commands allow the user to deviate from the default settings. This appendix is thus devoted to the description of these commands.

On-shell masses for W and Z boson
By default, the input masses for the W and Z bosons are the pole masses. Nonetheless, it is possible to set the on-shell masses instead by including the line RECOLA_ONSHELLZW=1 in the input run card.

Flavour scheme and quark masses
In order to allow all possible combinations of masses and flavour schemes, a few flags for the run card exist. The first one is

RECOLA_FIXED_FLAVS.
The values 4, 5, and 6 correspond to the corresponding fixed-flavour schemes. Setting the flag to 14, 15, 16 allows for a dynamical scheme up to the number of flavours (4, 5 or 6). Finally, the value −2 sets a fixed-flavour scheme according to the masses set in the run card.
It is also possible to set explicitly the masses used in the renormalisation [see Eq.
(2)] of the strong coupling using RECOLA_AS_REN_MASS_C/B/T, corresponding to the charm, bottom, and top-quark mass, respectively. On the other hand, the quark masses used for the hard matrix element are read out only from the run card. This means that the user has to take care of the consistency of her/his computation between the matrix element computed and the PDF set used.

UV and IR scales
The UV and IR scales are both set, by default, to a fixed value of 100 GeV. These technical parameters do not impact physical results. However, the choice of the IR scale does change the individual contributions of the virtual corrections and the integrated subtraction terms. In the Sherpa+Recola interface, it is possible to set these scales in the run card, e.g. to directly compare the virtual and real subtraction contributions separately to independent code. These UV and IR scales are set with the keywords UV_SCALE and IR_SCALE, respectively.