QCD matrix elements + parton showers: The NLO case

We present a process-independent technique to consistently combine next-to-leading order parton-level calculations of varying jet multiplicity and parton showers. Double counting is avoided by means of a modified truncated shower scheme. This method preserves both the fixed-order accuracy of the parton-level result and the logarithmic accuracy of the parton shower. We discuss the renormalisation and factorisation scale dependence of the approach and present results from an automated implementation in the Sherpa event generator using the test case of W-boson production at the Large Hadron Collider. We observe a dramatic reduction of theoretical uncertainties compared to existing methods which underlines the predictive power of our novel technique.


Introduction
Events with large final-state multiplicity, and in particular with many jets, are notorious in recent and present particle physics experiments at the energy frontier. This is especially true for measurements at the Large Hadron Collider, which operates at the highest centre-of-mass energies ever achieved in an accelerator-based experiment. Multi-jet events constitute a significant fraction of both signals and backgrounds. They must therefore be described with a theoretical precision that aims at matching the experimental one, or better. Triggered by this necessity, astounding practical improvements in event simulation have occurred over the past decade, with far-reaching consequences for particle physics phenomenology, including precision physics and searches for new phenomena alike. These developments range from the (near) automation of cross section calculations for multi-particle final states at the next-to leading order [1][2][3][4][5][6][7][8][9][10][11][12][13] to the construction of algorithms that allow to use multi-particle calculations at leading and next-to-leading order to drive event simulation at the particle level. For the latter point, it is necessary to consistently combine fixed-order matrix elements with parton showers. Two alternative ideas have been pursued in this regard: An exact matching of next-to-leading order calculations to parton showers has been worked out in the MC@NLO [14][15][16][17][18][19] and POWHEG [20][21][22][23][24][25] approaches.
For the latter, first steps to also incorporate electroweak corrections have been reported in [26,27]. In general, both these methods are suitable for simulating inclusive processes, but they fail to describe the precise kinematics of multi-jet final states as part of the inclusive reaction. Alternatively, algorithms to merge various leading-order matrix elements with increasing parton multiplicity into one inclusive sample have been proposed [28][29][30][31][32][33][34]. They describe the inclusive cross section at Born level only, but each additional jet is also described by exact leading-order matrix elements, while respecting the overall logarithmic accuracy of the parton shower. This is best achieved by implementing truncated vetoed parton evolution [32,33]. A combination of MC@NLO and MEPS, called MENLOPS, was recently introduced [35,36], which allows to promote the cross section of the inclusive MEPS sample to NLO accuracy. In this publication we propose a new method to merge next-to leading order predictions for the production of multi-jet final states in hadron-hadron collisions with the parton shower, maintaining both the fixed order accuracy for each jet multiplicity and the logarithmic accuracy of the parton shower. This is achieved by combining individual MC@NLO simulations [37] with suitably modified truncated parton showers. The actual implementation takes advantage of the various matrix element and parton shower generators which are implemented in the event generation framework SHERPA [38,39]. The main aim of our work is to reduce renormalisation and factorisation scale dependence of the fixed-order part in the simulation. We have elucidated the new merging technique, MEPS@NLO, for electron-positron annihilation in a parallel publication [40]. The present paper is therefore organised as follows: Section 2 reviews the basic formulae and discusses their extension to hadron-hadron collisions. For a more pedagogical introduction of the new method, based on the case of e − e + annihilations into hadrons, the reader is referred to a parallel paper, [40]. Section 4 presents first applications of the procedure to W -production at the Large Hadron Collider. Section 5 contains our conclusions and an outlook.
2 Multijet merging at next-to leading order -MEPS@NLO We have summarised the status of matrix-element parton shower matching and merging techniques in a parallel publication [40]. In the following, we only explain the basic ingredients and formulae which are needed for the MEPS@NLO merging method.

Ingredients of the method
One can formulate this technique in terms of the expectation value of an arbitrary, infrared-safe observable O, evaluated by taking the average over sufficiently many points in the multi-parton phase-space. Let us assume that m next-to-leading order predictions are merged, starting with an n-particle final state, and that p leading-order predictions are added on top. We consider an observable which is sensitive to at least (n + k) partons. This observable is determined as where O j are the contributions from j-parton final states. The respective partons can be generated by either the parton-level generator or the shower algorithm. We will not explicitly mention any parton-shower emissions beyond the fixed order or logarithmic order which we intend to maintain. For simplicity we assume that the (n + k)-particle final state is always described at next-to-leading order in the strong coupling, i.e. k ≤ m. Thus, for the observable O above, we have to trace up to three terms in the case that (n + k + 1)particle final states are generated at next-to-leading order (k < m) and up to two terms in the case that they are generated at leading order (k = m). The latter case was analysed extensively in [36] and shall not be reviewed here. In the following we will focus solely on the scenario that the (n + k + 1)-particle final state (and the corresponding (n+k +2)-particle real-emission contribution) is simulated with an NLO parton-level generator.
It is important to distinguish between partons and parton-level jets. We do not consider jets as defined experimentally, using any existing jet algorithm, but introduce a separate parton-level definition which is needed simply to separate parts of our inclusive event samples. To this end, we introduce the jet criterion Q. It is defined for any pair of final-state partons, labelled i and j, as For initial-state partons, labelled a, it is given by C k a,j = C k (aj), j where (aj) is the virtual t-channel parton originating from a and j. In all cases, the spectator index k runs over the remaining parton indices in the event. Any (n + k)-parton final state can be assigned a unique value of the jet criterion, Q n+k = Q(Φ n+k ), equal to the minimum of all Q ij constructed from the (n + k) partons. This value can be below or above a predefined value Q cut , which is called the merging scale. If it is above, Φ n+k is called an (n + k)-jet configuration.
In the context of next-to-leading order calculations, one also needs to define a criterion for the corresponding real-emission configuration to be of (n + k)-jet type. To this end, we reduce the (n + k + 1)-parton final state to an (n + k)-parton final state by backward clustering according to the kinematics of the parton shower. The pair of partons to be clustered and the corresponding spectator parton are identified to be the ones leading to the smallest jet measure. To simplify our notation, we define the j-jet inclusive and exclusive expectation values They include contributions from all final states with at least j partons, which must form j (but not j + 1 in the case of O excl j ) parton-level jets according to the definition of Q. The j +l-parton configuration is reduced to a j-parton configuration for the purpose of jet identification by an exact inverse parton shower [32]. In order to describe the computation of O j we need to introduce the following additional quantities • Squared leading-order (Born) matrix elements, B n+k (Φ n+k ), for n outgoing particles, summed (averaged) over final state (initial state) spins and colours, including symmetry and flux factors and multiplied by parton luminosities. The corresponding real-emission matrix elements are called R n+k (Φ n+k+1 ). Note that R n+k (Φ n+k+1 ) = B n+k+1 (Φ n+k+1 ).
• Sudakov form factors of the parton shower, which are given by where K n+k is the sum of all splitting kernels for the n-particle final state. The one-emission phasespace element for the splitting, dΦ 1 , is parametrised as where t is the ordering variable, z is the splitting variable and φ is the azimuthal angle.
For truncated showering [20,32], the corresponding Sudakov form factors are cut into two regimes, a hard one (>Q cut ) and a soft one (<Q cut ) such that ∆ . Pictorially speaking these two regions can be identified with emissions that populate the matrix element and the parton shower region, respectively. The latter are needed only in situations, where the parton shower evolution parameter t and the measure of hardness of an emission Q are not equivalent, e.g. for angular-ordered parton showers.
• The resummation scale µ Q , which defines an upper limit of parton evolution for the core process with multiplicity n.
• MC@NLO evolution kernels D n+k . They account for the first (hardest) emission in the MC@NLO generator, which is produced with a modified Sudakov form factor, ∆ (A) n+k (t, t n+k ). In our implementation, the kernels carry full colour information, which is needed to maintain NLO accuracy [37].
• The NLO-weighted Born differential cross sectionsB V n+k is the Born-contracted collinear-subtracted one-loop amplitude, while the sum of integrated subtraction terms is given by I n+k term in the square bracket) and off the intermediate lines. They populate different regions of phase space, as indicated by the Θ-functions in the evolution parameter t n+k+1 of the respective emission. Compared to an MC@NLO simulation, including the third line is necessary because the truncated parton shower can generate emissions in the region Q n+k+1 < Q cut , which alter the real-radiation pattern at O(α s ).
• The hard remainder functions which contains both the standard MC subtraction terms, D (A) n+k , and the subtraction terms for the truncated parton shower.
Note that for our purposes it is more useful to treat the expressioñ as a compound subtraction term, leading to a compound evolution kernel. It defines the Sudakov form factor Correspondingly we define a compound Sudakov form factor for the parton shower (2.10)

Definition of the method
In terms of these quantities, the exclusive expectation value of Eq. (2.3) in the MEPS@NLO method is determined to O(α s ) as follows: where t c is the infrared cutoff of the parton shower. Note that there is a slight mismatch between the terms D It is now explicit, that the Sudakov form factor on the first line accounts for a veto on emissions with Q > Q cut . This is exactly the same reasoning as in the MEPS method [32]. If the definition of hardness, Q, and the evolution parameter of the parton shower, t, coincide, we obtain in the same manner (2.13) We can write this in a form that makes the O(α s ) correction more explicit (2.14) Note that we have added arbitrary higher-order terms, which allow to include the sum over truncated shower subtractions inB n+k , since the two coincide if t n+k = Q n+k > Q cut and t n+k+1 = Q n+k+1 < Q cut , which is enforced by the two Θ-functions on the last line. The various terms are interpreted easily. The products in square brackets correspond to truncated vetoed parton showers, with their O(α s ) terms partially subtracted. In practice, these expressions can be generated by running a truncated vetoed shower and skipping/reweighting the first veto, depending on B n+k /B (A) n+k . The remainder of the expression corresponds to an ordinary MC@NLO simulation, consisting of S and H events. This scheme is particularly easy to implement in practice, because no emissions need to be generated in the truncated shower. In [40] this version of the method has been introduced and proven to be correct. In fact, it is worth stressing here that within SHERPA, indeed the definitions of t and Q are equivalent.

Fixed-order and logarithmic accuracy
The proof of next-to-leading order and logarithmic accuracy of the MEPS@NLO method in e + e − -collisions was presented in a parallel publication [40], for the simpler case where truncated showering effects can be neglected. First, the proof presented there will be extended to the case where truncated showering effects must be included. Proving fixed order accuracy is a fairly straightforward exercise. For this it is sufficient to expand Eq. (2.11) to the first order in α s . By construction, this yields exactly the same result as is obtained in MC@NLO for Q n+k+1 < Q cut . The effects of Sudakov suppression and modified subtraction cancel to first order in α s , and the correct fixed-order radiation pattern is recovered. The proof of logarithmic accuracy inherent to the parton shower is a bit more cumbersome. Firstly, we will show that below Q cut the method proposed here reproduces the formal accuracy of a corresponding MC@NLO approach for (n + k) jets. This guarantees that the first emission follows the corresponding treelevel matrix element for (n + k + 1) final state particles. Further emissions of course are generated by the parton shower and therefore, by construction, exhibit the correct behaviour. Then one needs to show that at the merging cut the combination of the (n + k)-jet and (n + k + 1)-jet exclusive samples does not generate unwanted terms in the (n + k + 1)-parton ensemble. As already noted before, in the (n + k)-jet contribution the second emission is generated by the parton shower and therefore correct by construction, while in the (n + k + 1)-jet contribution the (n + k + 2)th parton is generated through MC@NLO techniques, which, again by construction, maintain the logarithmic accuracy of the parton shower. Let us start the proof by constructing the inclusive observable where the MC@NLO description of the observable is obtained by dropping the Θ(Q cut − Q n+k+1 ) in (2.11). This contribution is given by (2.16) Before turning to the evaluation of the correction term O corr n+k , we note that O MC@NLO indeed is the MC@NLO result for the n + k-parton final state, up to a product of Sudakov form factors, which describe the implementation of truncated vetoed parton showers on the underlying Born configurations, starting from the n-particle final state. It is also worth noting that the emissions by truncated showering in the second line of this expression are, at O(α s ), compensated by the corresponding term in the hard remainder function. The correction term for the (n + k + 1)-parton configuration above the jet-cut, up to contributions of order α 2 s , and ignoring the effect of emissions below t n+k+1 , thus reads .

(2.17)
The relevant terms to consider are the ones in the curly bracket. The first as well as the second one consist of one factor directly responsible for the emission of an extra particle,D (A) n+k and B n+k+1 , respectively, which will eventually yield a contribution of O(α s L 2 ). Analysing the factors multiplying these emission terms, reveals that each of them is at most of O(α s L). However, these logarithms, if present, are associated with a factor of 1 N C , they arise from the difference between the evolution kernels in D (A) and D (PS) . Clearly, the combined virtual and real contributions toB n+k do not exhibit any logarithms that could upset the accuracy of the parton shower, because the phase space integrals over the real terms in (2.6) are unrestricted. Taken together, this shows that the correction term does not upset the formal logarithmic accuracy of the parton shower.

Renormalisation and factorisation scale uncertainties
The key aim of the MEPS@NLO approach presented here is to reduce the dependence of the merged prediction on the renormalisation scale µ R and the factorisation scale µ F , which are employed in the computation of hard matrix elements. These scales have not been made explicit so far. Note that only the dependence on renormalisation and factorisation scales are reduced compared to the MEPS method, while the dependence on the resummation scale, µ Q , remains identical. This is a direct consequence of the fact that the parton-shower evolution is unchanged in our prescription. The resummation scale dependence was analysed in detail in [37]. Following the MEPS strategy, the renormalisation scale should be determined by analogy of the leadingorder matrix element with the respective parton shower branching history [32]. In next-to-leading order calculations, however, one needs a definition which is independent of the parton multiplicity. The same scale should be used in Born matrix elements and real-emission matrix elements if they have similar kinematics, and in particular when the additional parton of the real-emission correction becomes soft or collinear. This can be achieved if we define the renormalisation scale for a process of O(α n s ) as Here, µ 2 i are the respective scales defined by analogy of the Born configuration 1 with a parton-shower branching history. The renormalisation scale uncertainty in the MEPS@NLO approach is then determined by varying µ R →μ R , while simultaneously correcting for the one-loop effects induced by a redefinition in Eq. (2.18). That is, the Born matrix element is multiplied by to generate the one-loop counterterm, while higher-order contributions remain the same. Similar reasoning holds for the collinear mass factorisation counterterms. Given µ F as determined by the MEPS algorithm, a different factorisation scaleμ F can be chosen, which leads to the O(α s ) counterterm where a and b denote the parton flavours of the initial state. The sums run over all parton flavours and P (z) denote the regularised Altarelli-Parisi splitting kernels. In the case of only one hadronic initial state, one of the sums would be absent.

Monte-Carlo implementation
In this section we describe the Monte Carlo implementation of the merging formula Eq. (2.11) in SHERPA. The techniques needed to combine leading-order matrix elements with parton showers are given elsewhere [32,34]. The algorithm reads as follows: • Draw an event according to the total cross section of the inclusive NLO-expression, effectively a sum overB (A) and H n+k .
• According to the absolute value of the relative contributions, select a standard (S) or hard (H) event. Select flavours and momenta accordingly.
• Reconstruct a parton shower history over tree-level like configurations.
• Start the parton shower, which will work out differently for S and H events.
For the latter, perform a truncated and vetoed shower in the spirit of the LO MEPS method, cf. [32,34].
If the matrix element level configuration could not be reached through standard parton showering (like, e.g. in the case of uū → W + sc) skip this step.
For the former, perform an MC@NLO step, i.e. generate an extra emission t n+k+1 through the D (A) n+k terms. Perform a truncated and NLO-vetoed shower between t n and t n+k . In contrast to a regular vetoed shower, in this NLO-vetoed shower, the first trial emission per cluster step, which is above Q cut will not lead to vetoing the event, but will be ignored, depending on the weight B n+k . This corresponds to the finite correction, i.e. the sum over the splitting kernels in (2.7). In [40], this corresponds to the correction term in the second line of Eq. (3.5). There, we also described in more detail some of the implementation details.
For both types of events continue with a vetoed parton shower, as in the MEPS case.

Results
In this section we present results generated with the previously described merging method. We employ the leading-order matrix element generators AMEGIC++ [43] and COMIX [44] in conjunction with the automated dipole subtraction provided in SHERPA [45] and the implementation of the Binoth-Les Houches interface [46] to obtain parton-level events at next-to-leading order. Virtual matrix elements for W + n jets are provided by the BLACKHAT library [5,8,47,48]. We employ a parton shower based on Catani-Seymour dipole factorisation [49] and the related MC@NLO generator [37] to generate events at the parton shower level. In contrast to the other MC@NLO implementations, no leading colour approximation is made. Our generator therefore recovers the full next-to-leading order accuracy of the fixed-order result throughout the phase space, to all orders in the colour expansion. We compare our predictions to a recent measurement [50] of W +jets events by the ATLAS collaboration. The analysis is used as implemented in the Rivet [51] framework and selects events containing a lepton within |η| < 2.5 with p ⊥ > 20 GeV and E miss T > 25 GeV. From the lepton and neutrino a transverse mass variable is calculated and required to fulfil m W T > 40 GeV. Jets are clustered using the anti-k t algorithm with R = 0.4. Effects from hadronisation and multiple parton interactions are not taken into account in the scope of this publication but can easily be enabled on top of the parton shower level studies pursued here. This is also justified by the fact that the two effects happen to compensate each other to a large extent in this analysis. At this point we would like to stress, again, that in our implementation the definitions of t and Q are equivalent, essentially transverse momenta of the respective splitting, such that truncated showering effects can be neglected. SHERPA predictions are made in three different approaches:

MC@NLO
NLO+PS matched sample for the W + 0-jet process using the MC@NLO-like implementation described in [37].

MENLOPS
The MENLOPS method described in [36] is used to merge an NLO+PS sample for the W +0-jet process with tree-level matrix elements for the higher multiplicity W + 1, 2, 3, 4-jet processes. Here we use this method on top of the W + 0-jet MC@NLO sample, as described in more detail in a parallel publication.

MEPS@NLO
The MEPS@NLO method was described in Sec. 2 and is used here for the W + 0, 1, 2-jet processes at NLO. In addition, the W + 3, 4-jet processes are merged using tree-level matrix elements via the MENLOPS technique.
For the two latter approaches we study perturbative uncertainties stemming from variations of the factorisation and renormalisation scale in the matrix elements. The scales are chosen by clustering the 2 → n parton level kinematics onto a core 2 → 2 configuration using a k T -type algorithm with recombination into on-shell particles. The central scale µ F = µ R = µ is then defined as the lowest invariant mass or virtuality in the core process. For core interactions which are pure QCD processes scales are set to the maximum transverse mass squared of the outgoing particles. A variation by a factor of two in each direction as µ F = µ R = µ 2 . . . 2 · µ generates the uncertainty bands of the predictions in the following. We start by looking at the cross section as a function of the inclusive jet multiplicity, Fig. 1. While the "pure" MC@NLO approach is able to describe the inclusive cross section at NLO in good agreement with data, it fails to provide a good description of the number of additional jets. This is to be expected, since the first jet is only described at leading-order accuracy, and any further jet even only in the parton-shower approximation. The MENLOPS method improves that prediction by including higher-multiplicity tree-level matrix elements, which can be seen to lead to a better agreement with data. It becomes obvious though that scale variations within those tree-level matrix elements lead to rather large uncertainties in all jet bins. This is significantly improved in the first two jet bins by the MEPS@NLO method, which has NLO accuracy in these observables and demonstrates the expected reduction in perturbative uncertainties. With the reduced uncertainty comes a near-perfect agreement with experimental data. For the higher jet multiplicities N jet ≥ 3 one recovers the tree-level picture with its larger uncertainties as expected. The picture is very similar when requiring jets with either p ⊥ > 20 GeV or p ⊥ > 30 GeV, also for all other observables we have studied. We thus restrict the plots in the following to the p ⊥ > 30 GeV case. In Fig. 2 the transverse momenta of jets in the different event categories of W + ≥ 1, 2, 3-jets are displayed. The inclusive statements from the previous paragraphs translate directly onto these differential distributions: When sensitive to the W + 1-or W + 2-jet processes only, the perturbative uncertainties become significantly smaller in the MEPS@NLO method leading to a much better agreement with data. A similar observation is made in Fig. 3 for the scalar sum of the transverse momenta of the lepton and jets and E miss T . The inclusiveness of this observable makes the hard tail susceptible to contributions from high jet multiplicities, which are only described at leading-order accuracy and thus cause a larger uncertainty band. In the low-H T region for W + 1-or W + 2-jet events on the other hand one can see the reduced scale band in MEPS@NLO. As a final comparison, the separation between the two leading jets is studied in the ∆R and ∆φ variables. It is clear that the pure MC@NLO approach can not be expected to give a good description of such angular correlations due to its parton shower approximation for the second jet. For the MEPS@NLO method it is impressive to see how the reduced scale uncertainty leads to very good agreement with the data also in these observables. Let us note that here we have analysed perturbative uncertainties stemming from the matrix element part of the calculation. While the inclusion of next-to leading order accuracy definitely reduces the related uncertainties, it is important to note that another source of theory ambiguities has not been discussed here. It is related to the effect of varying the resummation scale, here realised by the introduction of µ Q , in a way similar to the one performed in corresponding analytical calculations. It can be expected, though, that these uncertainties are fairly consistent between the MENLOPS approach and the new MEPS@NLO method presented here. While they certainly deserve further studies, we restricted our discussion here to the effect of theory uncertainties in the matrix element region, where we actually achieved a dramatic improvement over existing methods.

Conclusions
In this publication we have introduced a new method to consistently combine towers of matrix elements, at next-to leading order, with increasing jet multiplicity into one inclusive sample. Our method respects, at the same time, the fixed order accuracy of the matrix elements in their respective section of phase space and the logarithmic accuracy of the parton shower. The analysis of scale dependencies allows for a solid understanding of the corresponding theory uncertainties in the merged samples. Employing next-to leading order matrix elements leads, of course, to a dramatic reduction of the dependence on the renormalisation and factorisation scale and a much improved description of data. The same findings also apply to the case of e − e + annihilations into hadrons, cf. [40]. This allows, for the first time, to use Monte Carlo tools to generate inclusive multijet samples and analyse their uncertainty due to the truncation of the perturbative series in the matrix elements in a systematic and meaningful way.     (right). The predictions are compared to experimental data from ATLAS [50].