Giant QCD K-factors beyond NLO

Hadronic observables in Z+jet events can be subject to large NLO corrections at TeV scales, with K-factors that even reach values of order 50 in some cases. We develop a method, LoopSim, by which approximate NNLO predictions can be obtained for such observables, supplementing NLO Z+jet and NLO Z+2-jet results with a unitarity-based approximation for missing higher loop terms. We first test the method against known NNLO results for Drell-Yan lepton pt spectra. We then show our approximate NNLO results for the Z+jet observables. Finally we examine whether the LoopSim method can provide useful information even in cases without giant K-factors, with results for observables in dijet events that can be compared to early LHC data.


Introduction
At CERN's Large Hadron Collider (LHC), it is widely anticipated that signals of new physics, for example supersymmetry, may manifest themselves as large excesses of data compared to expected QCD and electroweak backgrounds at high momentum scales [1,2,3,4,5,6,7]. The estimation of these backgrounds will be one of the elements in ascertaining the presence of any new physics from such signals. Consequently, considerable effort is being invested across the particle physics community in the development of methods to understand and predict backgrounds (some of the issues involved are described nicely in ref. [8]). Given the QCD methods that are available today, some of the best prospects for obtaining systematic, accurate predictions of backgrounds involve next-to-leading order (NLO) QCD calculations. By carrying out a systematic expansion in the strong coupling  [19] for three observables in Z+jet production: the Z transverse momentum (left), the p t of the hardest jet (middle), and the scalar sum of the transverse momenta of all the jets, H T,jets (right). The bands correspond to the uncertainty from a simultaneous variation of µ R = µ F by a factor of two either side of a default µ = √ p 2 t,j1 + m 2 Z . The jet algorithm is anti-k t [20] with R = 0.7 and only events whose hardest jet passes a cut p t > 200 GeV are accepted. The cross sections include the branching ratio Z → e + e − . and obtaining the first two terms (leading order (LO) and NLO) for a given process, one often obtains predictions that are accurate to 10 − 20%. The importance of NLO predictions in the LHC programme has motivated a large calculational effort destined to extend the range of processes known at NLO (for reviews, see refs. [9,10]).
While the majority of NLO calculations show some degree of convergence relative to the LO results, several groups have commented in recent years on the appearance of K factors, ratios of NLO to LO results, that grow dramatically towards high transverse momenta [11,12,13,14,15,16] (similar behaviour is visible also in [17,18]). The problem generally occurs for hadronic observables (jet transverse momenta, etc.) in processes that involve heavy vector bosons or heavy quarks, at scales far above the boson or quark mass. Fig. 1 illustrates this for the pp → Z+jet process at LHC (14 TeV) energies. It shows the distributions of three observables that are non-zero for configurations involving a Z-boson and one or more partons: the transverse-momentum of the Z-boson (p t,Z ), the transverse-momentum of the highest-p t jet (p t,j1 ) and the effective mass (scalar sum of the transverse momenta) of all jets (H T,jets ). At LO, all three distributions are identical. At NLO, the p t,Z observable is rather typical of a QCD observable: its distribution has a NLO K-factor of about 1.5, fairly independently of p t,Z , and its scale dependence is reduced with respect to LO. The p t,j1 distribution is more unusual: at high p t it has a K-factor that grows noticeably with p t,j1 , reaching values of about 4 − 6, which is anomalously large for a QCD correction. The H T,jets observable is even more striking, with K-factors approaching 100.
Given that fig. 1 involves momentum scales where α s ∼ 0.1, one is driven to ask how it is that such "giant" K-factors can arise. As touched on in [13], and discussed in more detail in [14,15] for the p t,j1 case, the answer lies in the appearance of diagrams with new  Figure 2: A) a LO contribution to Z+jet production; B) and C) two contributions that are NLO corrections to Z+jet observables but whose topology is that of a dijet event with additional radiation of a soft or collinear Z-boson either from a final-state quark (B) or an initial-state one (C).
kinematic topologies at NLO. This is illustrated in fig. 2: at LO the only event topology (A) is that of a Z-boson recoiling against a quark or gluon jet. One type of NLO diagram involves gluon radiation from this basic topology, giving modest corrections to all our observables. However, there are also NLO diagrams (B,C) whose topology is that of a dijet event, in which a soft or collinear Z-boson is radiated from outgoing or incoming legs. These diagrams do not contribute significantly to the p t,Z distribution, because the Z-boson carries only a moderate fraction of the total p t . However when examining p t,j1 , it is irrelevant whether the Z boson is soft or not. Contributions B and C then lead to a result that is of order α 2 s α ew ln 2 p t,j1 /m Z , where the double logarithm comes from the integration over soft and collinear divergences for Z emission. The ratio of the NLO to LO results is therefore O α s ln 2 p t,j1 /m Z , 1 rather than just O(α s ), hence the K-factor that grows large with increasing p t . 2 For the H T,jets observable the enhancement is even bigger because the dijet topology leads to H T,jets ∼ 2p t,j1 instead of H T,jets = p t,j1 at LO.
While it is reassuring that we can understand the physical origins of the large Kfactors in fig. 1, we are still left with doubts as to the accuracy of the NLO Z+jet predictions for p t,j1 and H T,jets , since they are dominated by the LO result for the Z+2parton topologies. One way forward would be to calculate the full NNLO corrections for the Z+jet process. However, while work is progressing on NNLO calculations of 2 → 2 processes with QCD final states (see e.g. [22] and references therein), results are not yet available; nor are they likely to become available any time soon for some of the more complex processes where giant K-factors have been observed (e.g. some observables in pp → W bb [13,17]). Alternatively one could simply try to avoid observables like p t,j1 and H T,jets in inclusive event samples. For example, with additional cuts on the vectorboson momentum or a second jet, refs. [13,15] showed that the K-factors are significantly reduced. However, given the many analyses that are foreseen at the LHC, it is likely that at least a few will end up probing regions where giant K-factors are present.
To understand how else one might address the problem of giant K-factors, one can observe that in our Z + jet example, the bottleneck in obtaining a NNLO prediction is the inclusion of the two-loop 2 → Z + 1 parton contributions and proper cancellation of all infrared and collinear divergences. Yet the two-loop (and squared one-loop) contribution will have the topology of diagram A in fig. 2 and should not be responsible for the dominant part of the NNLO correction, which will instead come from diagrams with the topology of B and C, with either an extra QCD emission or a loop. So if one includes treelevel 2 → Z + 3 and 1-loop 2 → Z + 2 diagrams (i.e. Z + 2 jets at NLO) and supplements them with even a crude approximation to the two-loop 2 → Z + 1 result, one that suffices merely to cancel all divergences, then one should have a good approximation to the full NNLO result (a related observation has been exploited to obtain approximate NNLO results for high-p t J/ψ production in [23]).
The purpose of this article is to develop a general method for obtaining such rough estimates of missing loop corrections. Our approach, called LoopSim, will be based on unitarity. After explaining how it works in section 2, and outlining a secondary "referenceobservable" approach for control purposes, we will test the method by comparing its results to full NNLO predictions for lepton-p t spectra in Drell-Yan production in section 4, apply it to our Z+jet observables in section 5 and finally, in section 6, examine whether it can be of use even in the absence of giant K-factors, specifically for a number of dijet observables.

The LoopSim method
The main ingredient of the LoopSim method is a procedure for taking a tree-level event with n final state particles and supplementing it with a series of events with n−1 particles (approximate 1-loop events), n − 2 particles (approximate 2-loop events), etc., such that the sum of the weights of the full set of events is zero. This "unitarity" property will ensure that all the soft and collinear divergences of the tree-level matrix elements will cancel against identical divergences in the simulated loop contributions.
An outline of the procedure is given in fig. 3. Given a tree-level input event (a), the first step is to interpret it as a sequence of emissions (as if it had been produced by a parton shower), so that for example (diagram b) one can view particle 2 as having been emitted from particle 1, and particle 4 as emitted from the beam. The attribution of an emission sequence can be performed with the help of a suitable sequential-recombination jet algorithm and will be most meaningful in the limit that emissions are strongly ordered in angle and energy. The next stage is to decide which particles reflect the underlying hard structure of the event. If the event structure at the lowest possible order is that of a 2 → 2 scattering, then one should identify two outgoing "Born" particles. The Born particles will remain present in all the approximate "loop" events that are generated. They are represented as thick red lines in diagram (c). Again this step is most meaningful when all non-Born emissions are soft and collinear.
One then generates a set of simulated "1-loop" events by finding all ways of recombining one emitted particle with its emitter, diagrams (d,e). Each such "1-loop" event comes with a relative weight of −1 compared to the tree-level diagram. Similarly the set of simulated "2-loop" events is obtained by finding all ways of recombining two emitted particles with their emitter(s) (diagram f), each with relative weight +1; and so forth down to events where only Born particles remain (in fig. 3 this is already reached at the two-loop level). Note that the loop-diagrams drawn in fig. 3 are not intended to represent  Figure 3: Sketch of the LoopSim procedure as applied to a tree-level event (a) with 4 outgoing particles (numbered) and the beam (horizontal line); diagram (b) shows the attribution of the emission sequence, (c) the identification of the Born particles (thick red lines), and (d)-(f) the resulting "looped" diagrams. These diagrams are relevant in approximating next-to-next-to-leading corrections to a process whose LO contribution has a 2 → 2 structure. the actual Feynman diagrams that would be relevant at 1 and 2-loop level. Instead they indicate the way in which we have approximated the loop divergences, as the unitarising counterparts of the divergences that appear for each emission in the soft and collinear limits.
Given the above procedure for unitarising tree-level events, we shall see that it is then straightforward to extend it to event sets that also include exact loop diagrams.

The tree-level pure glue case
We start by examining the LoopSim procedure in the simple case of purely gluonic treelevel events. This will suffice to introduce most of the relevant concepts. Section 2.2 will then discuss some of the additional issues that arise for events with quarks and vector bosons, while the handling of events sets that include exact loop diagrams will be left to section 2.3.
It is helpful to introduce some notation: Firstly, b is the number of final-state particles present in the lowest relevant order (i.e. the number of final-state "Born" particles). For instance b = 2 if considering higher-order corrections to dijet events, as in fig. 3. E n represents a generic event with n final state particles. So the starting event of fig. 3 would be labelled E 4 . Finally, U b l will be an operator that acts on an event E n and returns all the events at l loops obtained from E n using the LoopSim method. For instance, fig. 3d,e represents the action of U b=2 l=1 on the input E 4 event (a). The central part of the LoopSim method involves the construction of the operator U b l acting on E n for all l = 0 . . . n − b (l ≤ n − b because the number of real final state particles cannot be smaller than that of the lowest order event).

Attribution of structure to events
Recall that the primary function of the LoopSim method is to cancel the divergences that appear in the soft and collinear limits. In these limits, events can be interpreted as stemming from a sequence of probabilistic (parton-shower) type 1 → 2 splittings of some original hard Born particles. The knowledge of the splitting structure will help us generate loop events to cancel the divergences.
The attribution of a branching sequence is most easily performed using a sequential recombination jet algorithm (and is inspired by the CKKW matching procedure [24]). We will use the Cambridge/Aachen (C/A) algorithm [25,26], which has important advantages over the k t algorithm when dealing with nested collinear divergences (avoiding "junk" jets [25]). 3 As a first step, to each of the i = 1 . . . n particles in the event E n , we assign a unique "identity" index I i ≡ i.
We then run the Cambridge/Aachen (C/A) algorithm on the event. It repeatedly clusters the pair of particles that are closest in angle, i.e. with smallest is the usual squared angular distance in the (y, φ) plane, and R LS is a free parameter, the LoopSim radius. The C/A algorithm continues until all the d ij > 1, at which point the remaining particles are deemed to cluster with the beam.
An ij → k clustering in the C/A algorithm can be reinterpreted as a k → ij splitting. The C/A algorithm does not distinguish in any way between i and j. However, in the soft limit, say p tj ≪ p ti , rather than viewing k as splitting to i and j, it is a better reflection of the divergent structure of the amplitude to view k as having emitted a soft gluon j. Then i is nothing other than particle k with some small fraction of its energy removed. To account for this, in an ij → k clustering, if p tj < p ti , then we declare that the "identity" I k of particle k should be the same as that of particle i, I k = I i . Also we record I i as being a "secondary emitter" and remember that the object with identity I j has been emitted from the object with identity I i . (Exchange i ↔ j if p ti < p tj ). This is represented in fig. 3b by the fact that particle 1 is a straight line, off which particle 2 has been emitted; the identity of the 1 + 2 combination is I 1+2 ≡ I 1 ≡ 1. For an iB clustering, we record I i as having been emitted from the beam. The next step in attributing structure to the event is to decide which event particles should be viewed as the Born particles, i.e. which particles are responsible for the hard structure in the event. Inspired by the original formulation of the Cambridge algorithm [25], for every ij → k recombination we assign a k t algorithm type hardness measure h ij = min(p 2 ti , p 2 tj )∆R 2 ij /R 2 LS [29,30]. 4 For every beam recombination, we assign a hardness h iB = p 2 ti . We then work through the recombinations in order of decreasing hardness. For an ij → k recombination (or k → ij splitting), assuming i is harder than j, we mark I k ≡ I i as a Born particle. If fewer than b particles have already been marked as Born particle, we also mark I j as a Born particle. For an iB recombination, we mark I i as a Born particle. This is repeated until b particles have been marked as Born (a particle may be marked more than once; in such a case its marking counts only once). As an example, in fig. 3, the hardest recombination will be between particle 3 and the beam, so particle 3 is marked as a Born particle. The next hardest recombination will that of (1 + 2) with the beam. Therefore we mark I (1+2) = I 1 = 1 as a Born particle. This exhausts the number (b = 2) of Born particles that need to be marked.
At the end of the above procedure, every particle will have been marked as emitted either from the beam or from another particle, and some particles will also have been marked as secondary emitters and/or Born particles. Thus in figure 3, particle 1 is labelled as having been emitted from the beam, it is a secondary emitter and a Born particle; particle 2 is labelled as having been emitted from particle 1; particle 3 is a Born particle, emitted from the beam; and particle 4 was emitted from the beam. The structure that we attribute is of course physically unambiguous only in the presence of strong ordering of emission angles and energies. However, as we shall argue in section 2.4, the mistakes that we make for non-ordered configurations should have a small impact for observables with giant K-factors.

Constructing virtual (loop) events
Once every particle is labelled in an event E n , one can compute the result of U b l (E n ), which is a set of events E n−l . For an event E n with respectively b Born particles and n s non-Born secondary emitters, we define to be the maximum number of particles that will be allowed to become virtual in a given event. It is obvious that Born particles will not become virtual. Additionally, secondary emitters will also not become virtual. To understand why, consider the event (2) in which particle 3 is a secondary emitter, since it emitted 4. There is a divergence for 4 to be collinear to 3 only if 3 is a final-state particle. If instead 3 is made virtual, then the divergence for emitting 4 no longer exists (there is no divergence for emission from internal lines in a diagram), so that the weight for the diagram in which 3 is virtual would be the weight of the tree-level diagram times a small coefficient ε ≪ 1. This simplest way of accounting for this is to approximate ε = 0 and thus not generate events in which secondary emitters are made virtual (a more detailed discussion is given in appendix C).
Having understood which particles can be made virtual, the operator U b l , when applied on an event E n , generates all the v l diagrams in which l particles become virtual. For the virtual events to cancel the infrared and collinear divergences that appear in the treelevel diagram, we need an infrared and collinear (IRC) safe procedure to make particles virtual. For instance, the divergent weight of an event with two collinear partons i and j has to be cancelled by that of corresponding virtual event (j makes a loop over i) when computing the distribution of any IRC safe observable; and two collinear partons, if not virtualised, have to remain collinear when another particle becomes virtual.
There are two ways for a particle j to make a loop: • If it is labelled as clustering with particle i, then one has to spread the momentum of particle j over i and all the particles that are labelled as clustering with it but which were emitted after j according to the C/A clustering sequence (i.e. at smaller angle). The exact procedure is explained in detail in appendix A, and is designed to ensure that the recombination maintains any collinearity properties of non-looped particles and is invariant under longitudinal boosts. When j is the only particle that clusters with i, then the procedure becomes equivalent to adding the momenta of particles i and j, p k = p i + p j , and then rescaling the momentum p k such that its mass is set to 0, while leaving its transverse components p x , p y and its rapidity unchanged.
• If particle j is labelled as clustering with the beam, then when it is "looped" it is simply removed from the event. Note that looping particles with the beam is less trivial than it may seem at first sight, because of an interplay with factorisation and the PDFs. Nevertheless it can be shown, appendix B, that for particle types that are included in the PDFs it does make sense to loop them. A p t imbalance will result from the looping of particles with the beam, and so after all loops have been made, we apply a transverse boosts to all remaining event particles, conserving their rapidities, so as to bring the total transverse momentum to zero (again, see appendix A).
There is some arbitrariness to our procedures for producing physical kinematics in the looped events. One avenue for future work would be to examine the impact of making different choices. The operator U b l has the following properties If w n is the weight of event E n , then each of the events generated by the U b l (E n ) operator has a weight w n−l = (−1) l w n .
Once all the U b l (E n ) have been calculated for l = 0 . . . n − b, one has to combine them in order to subtract all the soft and collinear divergences that appear in the calculation of E n and the virtual diagrams generated from it. This is done by the operator U b ∀ , which is defined as It generates all the necessary looped configurations that have the same order in α s as the original tree-level diagram. It is straightforward to see that the total weight of the diagrams obtained from the U b ∀ operator is 0. Indeed, if we apply it to an event E n whose maximum number of virtual particles is v, we get for v > 0. We note that the above procedure for approximating loop diagrams does not generate the finite terms needed to cancel the scale-dependence of lower-order diagrams. While it would be straightforward to include such terms, we believe that in the absence of full loop calculations, not including them helps ensure that the standard procedure of variation of renormalisation and factorisation scales is more likely to provide some form of reasonable estimate of the uncertainties on our results.

Some examples
In order to illustrate the action of the operator U b l , we give below some simple examples in the pure glue case. In each of these examples, only the Born particles are labelled with numbers Eq. (7a) gives an example of singly-looped configurations ("1-loop diagrams") generated by LoopSim when studying the 2 → 4 contributions to QCD dijet production. Eq. (7b) shows the "2-loop diagrams" generated from the same event. The next equation shows what happens if we add one more particle to the final state. If, eq. (7d), we now set the number of Born particles for the same event to be 3, we obtain only one 2-loop diagram instead of three, as represented in eq. (7d). 5 Finally, the last two examples of eq. (7) give a case with a splitting: the emitter is not looped, even if it is not a Born particle. We also give a few examples of the action of the U b ∀ operator: (8c) In the last case, only one particle can become virtual because there are two secondary emitters which cannot be looped.

Treatment of flavour within LoopSim
Let us now examine some of the issues that arise if we are to extend the LoopSim method to processes with quarks and vector bosons. We start with quarks and consider the situation depicted in fig. 4. In this case, applying the C/A algorithm as in the previous section will lead to the recombination of the two quarks q 1 and q 2 , which is clearly not physical. If flavour information is available 1 q q 2 Figure 4: Example of an event where two quarks q 1 and q 2 may get recombined by the C/A algorithm.
for events, then one can veto on such a clustering, for instance by defining the clustering distance d qq between two quarks to be infinite. As discussed in [31] such a modification alone is not sufficient to systematically guarantee sensible treatment of flavour in jet clustering. Refs. [31,32] have both discussed the further modifications needed in the case of the k t algorithm. A proper handling of flavour within LoopSim might seek to extend those modifications to the C/A algorithm. However, neither of the NLO programs that we use, MCFM and NLOJet++, provide information on particle flavours, so we defer such modifications to future work and just maintain the d ij = ∆R 2 ij /R 2 LS distance for all partons. For observables that are not flavour-sensitive this should not be a major drawback, given the observation [31] that divergences associated with the mistreatment of flavour are strongly subleading. Were we to be interested in heavy tagged quarks, more careful treatment might well be needed. Note that there are also subtleties related to flavour and PDFs, discussed in appendix B.
What about non-QCD particles, specifically vector bosons? Let us examine the case of Z bosons. A Z can be emitted from quarks or antiquarks and we would like this to be reflected when establishing the approximate emission sequence, because if the Z has been emitted from a quark, then that quark is a secondary emitter and should not be looped. In other cases a Z boson may be the hardest isolated object in an event. Then it is to be considered a Born particle. On the other hand we won't necessarily wish to consider diagrams where a Z boson is looped, because they would represent electroweak corrections, not QCD corrections.
One issue in dealing with electroweak particles is that they are not emitted from gluons. If one could distinguish between quarks and gluons, then this could be accounted for during the C/A clustering, by defining the distance d gZ between a Z and a gluon to be infinite. Since we will not know which partons are quarks or gluons, we adapt Frixione's isolation procedure [33] to decide if a Z boson relatively close in angle to a parton i is likely to have been emitted from i. More precisely, if then we define d iZ = ∆R 2 iZ /R 2 LS , otherwise d iZ = ∞. When recombining i and Z into a particle k, then the identity index I k is set equal to I i (a quark and a Z give a quark). Our procedure means that a Z that is very collinear to a parton is always considered to be emitted from that parton -this makes sense because such configurations are much more likely to occur when the parton is a quark. In contrast a soft parton in the general vicinity of a Z is not clustered with the Z, which is sensible given that most soft partons tend to be gluons. Finally, for a recombination between a parton i and Z, we define the hardness of the branching h iZ as 6 while a recombination of a Z with the beam has a hardness The latter means that for an event with just a parton and a recoiling Z boson, the parton's beam hardness will always be lower than the Z's, implying that for b = 1 it is the Z-boson that will be the single Born particle, as should be the case, at least when the parton has p t ≪ m Z , i.e. in the kinematic regime that dominates the total cross section for the Z.
Once a structure has been assigned to an event with a Z boson, the next question is that of the looping procedure. When looping partons it remains identical to before, with just a small extension of the recoil procedure in order to deal with decay products of the Z boson (see appendix A). In the situations where the Z is not a Born particle (it is never an emitter), straightforwardly following the procedure of section 2.1.2, one would deduce that one should loop the Z as well: (straight lines are partons, either quarks or gluons). The rightmost diagram, with the looped Z, is not, however, a QCD loop diagram: it is an electroweak loop correction to a multijet event. The LoopSim procedure does not aim to reproduce electroweak loop corrections (though in this case it might be a reasonable approximation). Furthermore, in any analysis that tags on Z bosons, such a diagram would not be tagged and so would not contribute. Thus, although the LoopSim procedure naturally generates events with looped Z bosons, events like the rightmost diagram of eq. (12) are simply to be discarded. For events with W ± bosons, the same procedure can be used as for Z's. Note, however, that while the "looped" Z events may give a reasonable approximation to actual electroweak loop diagrams, looped W ± events will not. This is because W-boson emission changes quark flavour: consider a tree-level diagram ud → bbW + , with the W + emitted collinearly off the incomingd, converting it into aū. The LoopSim procedure would give a "loop" diagram ud → bb, with the W + looped. However no such loop diagram exists and the correct loop diagram would instead involve uū → bb. 7 This is closely related to the phenomenon of Bloch-Nordsieck violation [34] that is found when considering electroweak double logarithms. Since we in any case discard events in which electroweak bosons are looped, this should not be a problem for the practical use of LoopSim in events with W ± bosons.

Merging NLO calculations and beyond
Before explaining how we merge exact higher orders calculations, let us mention how we use the LoopSim method in practice on tree-level events at several different orders. We introduce the notation X@n p LO to denote an approximation to the N p LO cross section for producing X, with all loop terms estimated through the LoopSim procedure. It is obtained by applying the U b ∀ operator to all tree-level diagrams that can contribute up to N p LO. For instance, one can write Notice that U 1 ∀ (Z@LO) = Z@LO and U 2 ∀ (Z+j@LO) = Z+j@LO. The terms U 1 ∀ (Z+j@LO) and U 2 ∀ (Z+2j@LO) simulate up to one-loop corrections, and U 2 ∀ (Z+3j@LO) simulates up to two-loop corrections. Now let us see how things work beyond tree-level accuracy. We define E n,l to be a generic event at l loops (exactly calculated) with n particles in the final state. We first consider the case where only one-loop corrections are computed exactly, so that we have tree-level events E n,0 and exact one-loop events E n−1,1 . As before we can apply the unitarisation operator to the tree-level events, U b ∀ (E n,0 ). However, since we now include exact 1-loop contributions, E n−1,1 , we must remove the approximate 1-loop contributions . This alone is not sufficient, because among the extra contributions from the exact 1-loop terms, there will be pieces that are finite for a given (n − 1)-parton configuration, but that can lead to divergences when integrated over the (n − 1)-parton phase space. To cancel these extra divergences, we should introduce additional approximate higher-loop contributions, which can be obtained by applying the unitarisation operator U b ∀ to the difference between the exact and approximate one-loop terms. So, rather than including just events E n−1,1 and subtracting U b where the extra subscript 1 on the U b ∀,1 indicates that it is the form to use when the exact 1-loop result is to be included. The action of U b ∀,1 depends on the number of loops already included in the event on which it operates: we subtract the one-loop contribution returned by LoopSim only in tree-level events. With this notation, one can compute the higher order corrections to eqs. (13) to one-loop accuracy, where the "only" subscript on Z+nj@NLO only means that we take the highest order that contributes, i.e. here α n+1 s α ew , since the LO, α n s α ew , piece of Z+nj@NLO, is already taken into account in the Z+(n − 1)j@NLO contribution. This implies that one should use consistent renormalisation and factorisation scale choices across all different orders of the calculation. Note that in eq. (15) we have introduced the notationn p N q LO to denote an approximation to the N p+q LO result in which the p highest loop contributions have been approximated with LoopSim.
The extension of the procedure beyond one-loop accuracy is simple. For instance, at two-loop accuracy, one has to subtract the approximated two-loop contribution , and the other approximated two-loop contribution U b 1 (E n,1 ) in eq. (14b), giving Therefore, once Z+j@NNLO is calculated, one may compute for instance To be complete, let us mention the generalisation of our procedure to m-loop accuracy We noted at the end of section 2.1.2 that the plain LoopSim procedure does not generate the finite terms needed to cancel residual scale dependence from lower orders. With the introduction of the exact loop contributions, those finite terms do now get included. Thus for a given number of exact plus simulated loops, as we increase the number of exact loops, we should expect to see reductions in scale uncertainties.

Expected precision of the method
Let us briefly explain why the LoopSim method is expected to work in the presence of giant K-factors. We consider an observable A computed respectively at NLO andnLO. We define K and we assume that K nLO gives the exact real part of the NLO calculation as well as the divergent terms of the virtual correction. Therefore where, in writing O(α s σ (A) LO ), we mean that the term missing in thenLO calculation, the finite part of the 1-loop correction, is not especially enhanced. This leads to The relative difference between the approximate and exact NLO calculations is thus suppressed by the inverse K-factor.
Next, considernNLO accuracy. The difference between σ NNLO comes from the parts of the two-loop corrections that are finite and associated with the LO topology, so that they should be free of the enhancements that led to the large NLO K-factor. This implies σ If we define K If K

The reference-observable method
Given the novelty of the LoopSim method, it is useful to have an alternative way of estimating the size of the NNLO contributions that we will approximate with LoopSim.
Here we outline such an alternative method, which, though less flexible than the LoopSim approach, will provide a valuable cross-check and help us build our confidence in results of the LoopSim method. Let us explain it for observables in the Z+j process. Our aim is to estimate σ NNLO for some observable A. 8 We assume that we have a reference observable which is identical to the observable A at LO. For instance, one might consider ref = p t,Z and A = p t,j . We can write the NNLO Z+j prediction for A in terms of the NNLO prediction for the reference observable plus the NLO Z+2j difference between A and the reference cross section The second equality is possible because 2-loop NNLO corrections to Z+j have the topology of Z+j at LO. Therefore, their contributions to the observables A and ref are identical and cancel in the difference in eq. (24a). If we have reason to believe that the perturbative expansion for the reference observable converges well, we can conclude that σ Z+j@NLO is genuinely a small correction. Then i.e. we approximate the NNLO distribution for A in terms of the NLO distribution for the ref observable and a NLO calculation for difference between the A and ref distributions, both of which are exactly calculable. The missing part is suppressed by a relative factor 1/K (A) , as for the LoopSim method. For Z+j, one can see from fig. 1 that p t,Z seems to be an acceptable reference observable for p t,j and H T,jets .
In the sections that follow we shall, for brevity, refer to the RHS of eq. (25) as "ref. nNLO" even though it does not quite adhere to our the meaning ofnNLO as set out in section 2, i.e. in terms of the specific sets of tree-level and loop diagrams that are included exactly.

Validation: comparison to DY at NNLO
The cross section for the Drell-Yan process is known with exclusive final states up to NNLO accuracy [35,36]. Above a certain value of lepton transverse momentum, one finds giant corrections to the lepton p t spectra when going from LO to NLO and large ones from NLO to NNLO. This gives us an opportunity to directly test the performance of the LoopSim method by comparing itsnNLO results to exact NNLO spectra for lepton pair production.
Before examiningnNLO results, it is useful to comparenLO with NLO. If they are in reasonable agreement for some observable, then that serves as a first indication that the LoopSim estimate of missing loop corrections is sensible for that observable. Fig. 5 gives the comparison of thenLO, NLO and LO results for the production of an e + e − pair within the mass window of 66 < m e + e − < 116 GeV at a proton-proton centre of mass energy of 14 TeV. The left-hand plot shows the cross section differential in the transverse momentum of the harder of the two leptons. The right-hand plot gives the corresponding K factor with respect to LO. The results were obtained with MCFM 5.3 [37,19], with its default set of electroweak parameters and NNLO MSTW2008 parton distribution functions. The uncertainty bands in Fig. 5 correspond to varying the renormalisation and factorisation scales µ r = µ f by a factor of 1 2 and 2 around a default choice of m Z . In thenLO result we fixed the value of the LoopSim radius parameter to be R LS = 1, which naturally places interparticle and particle-beam clustering on the same footing (though thenLO result here is actually independent of R LS , because there is at most one isolated QCD parton in the final state).
There are three relevant regions of transverse momentum in fig. 5. For p t,max 1 2 m Z (low p t ) the distribution is dominated by on-shell Z-bosons and its shape is governed by the angular distribution of the Z decays in their centre-of-mass frame. The peak close , the LO distribution comes from Z-bosons that are off shell, which allows the p t of the lepton to be larger than 1 2 m Z . The narrow width of the Z causes the distribution to fall very steeply. The 58 GeV upper edge of this region is a consequence of our cut on m e + e − < 116 GeV. Above 58 GeV (high p t ) the LO distribution is zero.
In the low p t region, the NLO correction is moderate and negative. There is no strong reason to believe that the LoopSim method should work here, but it turns out that thenLO result reproduces the structure of the correction, even if its scale dependence remains much larger than that of the NLO result (this is because the LoopSim procedure does not include the finite terms that would partially cancel the LO scale dependence). In the intermediate p t region, we see a "giant" NLO K-factor. It comes about because initial-state radiation can give a boost to the Z-boson, causing one of the leptons to shift to higher p t (it becomes the "max" lepton). The spectrum of QCD radiation falls much less steeply than the Z-boson lineshape, so this NLO correction dominates over the LO result. In this region the exact loop correction, proportional to the LO result, becomes almost irrelevant and we see near perfect agreement betweennLO and NLO. In the high-p t region only the real emission diagrams of Z@NLO contribute andnLO becomes identical to NLO (both correspond to the Z+j@LO result). Similar results hold for the p t,e ± distribution, while the p t,min lacks the giant K-factor in the intermediate region.
A similar comparison betweennNLO and NNLO spectra is shown in fig. 6. The NNLO results were obtained with DYNNLO 1.0 [36,38,39], used with a set of electroweak parameters compatible with that of MCFM. 9 In the low-p t region we find quite good agreement between thenNLO and NNLO results (with somewhat larger uncertainty bands fornNLO). Such a result was not guaranteed a priori, even if it is not entirely surprising given the reasonable agreement that we saw betweennLO and NLO. In the intermediate p t region, where the NNLO/NLO corrections are substantial, the agreement is excellent. This was expected. At high p t the agreement should be exact, and does seem to be, within statistical fluctuations. The dependence on R LS (shown in the right-hand plot) has been estimated by varying its value from 0.5 to 1.5. The effects are small.
Finally, we note that similar features and a similar level of agreement betweennNLO and NNLO are to be found in the p t,min and p t,e ± distributions.

Results for the Z+jet process
In the previous section, we studied the Z production process and showed that our procedure correctly reproduces the p t distribution of the hardest lepton at NNLO, even, unexpectedly, in regions where the K-factor is not large. In this section we study the Z+j process, whose NNLO cross-section is not known yet, but which leads to giant Kthe Z boson. The cut is applied to both real and virtual terms and its impact should vanish as it is taken towards zero. It is, however, required to be non-zero for the numerical stability of the MCFM Z+j NLO calculation that is among the components of DYNNLO. We set the cut equal to 0.1 GeV in the O(α s ) term and to 1 GeV in the O α 2 s term. A related 1 GeV cut was placed on the O α 2 s piece of thenNLO result (while none was used at O(α s )). The impact of the 1 GeV cut is small but not entirely negligible close to the peak (where, physically, NNLO should in any case be supplemented with a resummation). factors at NLO for some observables as explained in the introduction. Therefore, their NNLO contributions are expected to be accurately described by the LoopSim method. Throughout this section we use MCFM 5.7, including the Z+2j process at NLO [40], with the NLO CTEQ6M PDFs. We will take three different values for the renormalisation and factorisation scales: µ r = µ f = 1 2 µ 0 , µ 0 and 2µ 0 , with where p t,j1 is the transverse momentum of the hardest jet. At high p t , this scale choice should be quite similar to that used in [41] and has the same p t scaling as those in [14,15]. The R LS uncertainty is measured at µ r = µ f = µ 0 using three different values for it: R LS = 0.5, 1, 1.5. In addition to the 3 observables shown in the introduction, p t,Z , p t,j1 and H T,jets = ∞ i=1 p t,ji , we will also consider We only include events for which p t,j1 > 200 GeV.

Validation atnLO
As a first investigation of the performance of the LoopSim method, let us examine how thenLO approximation compares to the full NLO result. Fig. 7 shows the K-factors for thenLO and NLO predictions, with uncertainty bands from scale and R LS variations.
In the upper-left plot, one sees that thenLO prediction for the p t,Z distribution gives a somewhat smaller K-factor than the NLO result. We interpret this as being because certain genuine loop effects are not taken into account by the LoopSim method, for example those related to threshold logarithms, which depend crucially on the factorisation scheme of the parton distribution functions. ThenLO result does, however, reproduce the p t dependence of the K-factor, i.e. the dip towards p t = 200 GeV. This dip arises because of the requirement in our event selection that there should be at least one jet with p t > 200 GeV. At LO this induces a step-function in the p t,Z distribution at 200 GeV. At NLO, soft and collinear emissions smoothen out that threshold and thenLO calculation correctly reproduces the resulting interplay between real and virtual terms.
In the three remaining plots of fig. 7, for p t,j1 , H T,jets and H T,tot , all of which have giant K-factors, one sees good agreement between thenLO and NLO results. This is because the dominant NLO contribution comes from events in the B and C-type configurations of fig. 2, for which there is no corresponding QCD loop correction. The LoopSim method merely serves to cancel the divergences that arise from soft and collinear emissions off A-type configurations and these are not dominant overall.
The R LS dependence, also shown on these four plots, only comes from 1-loop events generated by LoopSim. Therefore, for an observable A studied in Z+j@nLO with two different values R 0 and R 1 for R LS , one can write: . This means that the absolute uncertainty due to R LS is the same for A and p t,Z . Therefore, the relative uncertainty due to R LS is expected to be roughly inversely proportional to the K-factor for A, in analogy with the discussion of sec. 2.4. This explains why the R LS dependence (solid cyan band) looks significantly smaller for p t,j1 , H T,jets and H T,tot than it does for p t,Z plot.

Results atnNLO
Results atnNLO are given in fig. 8. In the case of p t,Z the result is similar to the NLO result, and the scale uncertainties remain largely unchanged. In other words, since Z+2j topologies do not dominate the high-p t,Z distribution, adding NLO corrections to them (i.e.nNLO Z+j) makes no difference either to the result or to the uncertainties. We have also shown the dependence on the choice of R in the LoopSim procedure. It is smaller  Figure 9: Comparison between the approximate NNLO/LO K-factor calculated using respectively the LoopSim and the "reference-observable" method for p t,j1 and H T,jets . As a reference observable we have used the differential cross section for p t,Z .
than the scale dependence.
The p t,j1 distribution gets a correction that is just within the NLO uncertainty band, withnNLO uncertainties that are about half the size of the NLO band. Adding in thē nNLO term has made a real difference. This is precisely what we expect: the observable is dominated by Z+2-parton configurations, and these were only present at tree-level in the NLO Z+j calculation. Our use ofnNLO provides the additional 1-loop Z+2-parton and tree-level Z+3-parton configurations that come with NLO Z+2j accuracy.
Given the improvement in scale uncertainty, we need to ask whether the uncertainty due to R LS variation might somehow eliminate part of this benefit. It is, however, small. The reasons are similar to those given around eq. (28).
The H T,jets and H T,tot distributions get significantnNLO corrections, withnNLO/NLO K-factors of about 1.7 − 2. Absolute scale uncertainties increase slightly compared to NLO, but because of the large K-factor, relative scale uncertainties diminish. At first sight, it is somewhat disturbing that thenNLO and NLO uncertainty bands don't overlap. Given the novelty of the LoopSim method, one should therefore ask whether this is reasonable and whether there is any way of cross-checking the result.
A first observation is that sincenNLO Z+j is really NLO of the dominant Z+2j component, the largenNLO corrections that we see are comparable to an O(2) K-factor for going from LO to NLO in the Z+2j prediction. There are many contexts where NLO and LO results are not compatible within scale uncertainties, and so it is not unreasonable that the same should be seen here.
Still, we would like to have some more quantitative cross checks that our results are sensible. One option is to consider the alternative "reference-observable" method presented in section 3, which only makes use of standard NLO calculations to compute the approximate NNLO corrections. The comparison between the two methods is shown in fig. 9 for H T,jets and p t,j1 , where we have taken p t,Z as the reference observable. One notices near perfect agreement for H T,jets and very good agreement for p t,j1 . This gives us some degree of confidence that thenNLO LoopSim results provide an accurate description of the NNLO behaviour for these observables.
A second option for cross-checking the largenNLO effects for H T,jets and H T,tot , is to examine whether H T type observables might generally be "difficult". To do so we look at them in the case of QCD jet events.

QCD jet events as a testing ground
We have seen that thenNLO K-factors for the two effective-mass variables, H T,tot and H T,jets , in Z+jet(s) events are about a factor of two above the NLO K-factor.
Since NLO is the first order at which we see the dominant "dijet" topology for the H T variables in Z+jet(s), fig. 2B,C, it might be instructive to establish a correspondence with a simpler process, QCD dijet production. Having a NLO Z+j prediction is analogous to a LO dijet prediction; and thenNLO Z+j predictions should be analogous to NLO dijet predictions. NLO cross sections for dijet observables can be calculated exactly and therefore we can check whether NLO K-factors of order 2 appear for effective-mass observables in pure QCD events.
We will consider several effective-mass observables: an H T,n variable, which sums over the n hardest jets above some threshold (p t,min = 40 GeV; such a cut is often imposed experimentally 10 ) where p t,i is the transverse momentum of the i th hardest jet. Upper limits on the number of jets included in the effective mass are common in SUSY searches [3,4]. We also define an effective mass for all jets above the p t,min threshold, which is similar to the H T,jets and H T,tot observables of section 5. Finally, for completeness we will consider the distributions of p t,j1 , p t,j2 and the inclusive jet spectrum. All our results in this section will be for a centre-of-mass energy of 7 TeV, to allow comparison to results in the current run of the LHC. At LO, the distributions of 1 2 H T,n (n ≥ 2), 1 2 H T , p t,j1 , and p t,j2 will all be identical. The inclusive jet spectrum will have a distribution that is twice as large (because each of the two jets contributes). Note that we do not impose any rapidity acceptance limits on the jets: though such a cut would have been trivial to include in the LoopSim procedure, it would have complicated somewhat the reference-observable approach that we will consider at the end of the section. LoopSim results with a rapidity cuts of |y| < 2 on the jets are available from the authors on request. Figure 10(left) shows the distributions for two observables, 1 2 H T and p t,2 at LO (where they are identical) and at NLO, as determined using NLOJet++ [42,43] with CTEQ6M PDFs. A first comment is that H T receives a NLO K-factor of order 2, just like thenNLO −< µ − p t,j1 <2 Figure 10: Left: differential cross sections for the p t,j2 and 1 2 H T observables, at LO, where they are identical, and at NLO where they have substantially different K-factors. Right: the NLO K-factors for the 400 < V / GeV < 500 bin for each choice of variable V among the following: the inclusive jet spectrum, the p t distribution of the hardest (p t,j1 ) and second hardest (p t,j2 ) jets, (half) the effective mass of the two hardest jets (H T,2 ), three hardest jets (H T,3 ) and of all jets above 40 GeV (H T ). Also shown on the right are thenLO results for the K-factors. The NLO andnLO (µ) widths correspond to the uncertainty due to simultaneous renormalisation and factorisation scale variation by a factor of two around a central value µ = p t,j1 . ThenLO(R LS ) width shows the uncertainty from a variation of R LS in the range 0.5 < R LS < 1.5. enhancements in the Z+j case. This provides supporting evidence as to their legitimacy. A second comment is that the cross sections are large: these observables will be easily accessible with a few pb −1 of integrated luminosity at a 7 TeV LHC, allowing for an early experimental verification of the large K-factor for H T .
The other observable in the left-hand plot of fig. 10, p t,j2 , has a very different K-factor, somewhat below 1. The right-hand plot shows the NLO K-factors for our full range of observables, focusing on a single bin of the left-hand one, from 400 − 500 GeV. The pattern that we see here allows us to make some deductions. Firstly, the H T,2 variable, which sums the p t 's of the two leading jets, is free of large NLO enhancements. It is the addition of the third jet in H T, 3 and H T that brings about the enhancement. A natural interpretation is the following: it is common for a third, soft jet to be present due to initial state radiation. This third jet shifts the H T distribution to slightly larger values, and because the distribution falls very steeply, that leads to a non-negligible enhancement. This suggests that if, in section 5, we had used effective mass observables with at most two objects in the sum, then thenNLO/NLO ratios would have been close to 1. We have verified that this is indeed the case.
The pattern for p t,1 and p t,2 in fig. 10 can also be explained in similar terms: a soft ISR emission boosts the hard dijet system, breaking the degeneracy between the p t 's of the two hardest jets. It is jet 1 that shifts to larger p t (giving a K-factor > 1), while jet 2 shifts to lower p t and so it gets a K factor below 1. For the inclusive jet spectrum, and for H T,2 , this effect balances out. In addition, final-state radiation from one of the jets can cause it to shift to lower p t (becoming the 2nd jet), further reducing the K-factor for the distribution of p t,j2 . Of the different variables, it is only the inclusive jet p t and H T,2 for which there is a clear reduction in scale uncertainty in going from LO to NLO. Figure 10(right) also shows thenLO results (including uncertainties both from scale variation and from the LoopSim parameter R LS ). Despite the fact that none of the Kfactors is parametrically large (except arguably for H T,3 and H T ), thenLO results are remarkably effective at reproducing the pattern of NLO K-factors, albeit with a small systematic shift and generally larger scale uncertainties. One can also verify that, to within 10 − 20%, the p t dependence of the NLO K-factors is reproduced atnLO.
Given this success ofnLO, and the observed limited convergence of some of the observables at NLO, it is interesting to examine what happens atnNLO, where the additional 3j@NLO contribution that we require is again obtained using NLOJet++. Results are shown in fig. 11.
For the inclusive jet spectrum and H T,2 , which already saw large reductions in scaledependence at NLO, thenNLO corrections have essentially no meaningful effect: they neither significantly affect the central values, nor reduce the scale uncertainties. For these observables, NLO already converged well, and adding a subset of the NNLO corrections without the 2-loop part cannot improve the result.
For the other effective mass observables, the situation is quite different. With H T,3 , thenNLO result is close to the NLO result and the scale uncertainty is much reduced, i.e. this observable seems to come under control atnNLO. In contrast, H T is subject to quite a large further correction, with the central value atnNLO lying outside the NLO uncertainty band, and thenNLO uncertainty band (dominated by scale variation) only marginally smaller than at NLO. Why is this? Perhaps we are seeing the effect of a second ISR emission, which shifts the H T distribution to even higher values? Given that H T,3 converges and H T does not, such an explanation is not unattractive. It is also consistent with the decrease in K-factor at low H T , where the 40 GeV p t cutoff on the jets contributing to the H T sum will eliminate the ISR enhancement. A definitive conclusion would however probably require further study.
For the remaining two observables, p t,1 and p t,2 , thenNLO contribution goes in the opposite direction from the NLO correction and at low p t it seems that the series fails to converge. This is, we believe, closely related to observations of insufficiencies of NLO predictions for dijet cross sections in DIS and photoproduction when identical p t cuts are imposed on both jets [44,45,46,47,48] (equivalent to integrating the p t2 distribution above that cut). The worse convergence at low p t is probably due to the larger fraction of subprocesses that involve gluons in the underlying 2 → 2 scattering, so that perturbative corrections tend to go as (C A α s /π) n rather than as (C F α s /π) n at higher p t .
Considering that we do not have giant NLO K-factors for the jet processes shown here, one may question the validity of the information obtained from the LoopSim procedure. An important cross check comes from a comparison with the reference-observable technique. Examining fig. 10 (right), one sees two natural reference observables: the inclusive jet spectrum and H T,2 , both of which show "perturbative" K-factors and small scale dependence at NLO. Here we will use (half) the inclusive jet spectrum as the reference observable (results with H T,2 would be almost identical). Figure 12 provides a comparison of the LoopSimnNLO results (showing the envelope of the scale and R LS uncertainties) with the reference-observablenNLO results. The comparison is given for all observables except the reference observable itself. The agreement between the two methods is striking, with the reference-observable method giving just a small shift of the K-factors relative to the LoopSim results. The shift is identical for all the observables, as it has to be: it is simply equal to the difference between the NLO andnNLO results for the reference observable. Insofar as we believe the scale dependence to be representative of the true NLO uncertainty on the inclusive jet spectrum, 11 the results for the other observables should therefore be good approximations to the full NNLO results.

Conclusions
Several cases of LHC observables with giant NLO K-factors have come to light in recent years. They are characterised by the presence at NLO of new partonic scattering topologies that have large enhancements over the LO topologies. In these cases, NLO calculations, while important in highlighting the presence of the large K-factors, cannot on their own provide accurate predictions.
In this article we have examined how to address this problem by combining NLO results for different multiplicities, for example Z+j@NLO with Z+2j@NLO. Our main, most flexible method, LoopSim, makes use of unitarity to cancel the infrared and collinear divergences that appear when one tries, say, to apply Z+2j@NLO calculations to observables that are non-zero starting from Z+1-parton. We referred to the result as Z+j@nNLO, where the "n" indicates that the highest loop contribution to the NNLO result (the two-loop part) has been estimated with LoopSim.
In introducing a new approximate method for estimating NNLO corrections, significant evidence needs to be provided that the method is meaningful. Firstly, we gave reasons why, in cases with giant K-factors associated with new NLO topologies, we expectnNLO results to be a good approximation to NNLO results. As a next step, we carried out studies comparing Z/γ * @nNLO (DY) to NNLO predictions for the pp → Z/γ * + X → e + e − + X process. In comparing the DY lepton p tn NLO distributions to NNLO we found near-perfect agreement in a region of giant K-factors, p t − 1 2 m Z Γ Z . Interestingly, even in the region where the NLO K-factor was not large, p t 1 2 m Z , thē nNLO results provided a significantly better approximation to NNLO than did the plain NLO result. This need not always be the case, but is, we believe, connected to the ob-servation that ournLO results reproduced much of the structure seen at NLO (recall, Z@nLO means combining Z@LO with Z+j@LO).
For Z+j production, the first step of our validation procedure was to comparenLO and NLO results. All observables with giant K-factors showed good agreement between the two (one with a moderately large K-factor did not). For those observables,nNLO always appeared to provide extra information: either suggesting a convergence of the perturbative series, with reduced scale uncertainties (for p t,j1 ), or an indication of substantial further higher order corrections (for the effective-mass type observables H T,jets and H T,tot ). Almost identical results were seen with our alternative "reference-observable" estimate of the NNLO contribution.
The largenNLO corrections that we saw for effective mass observables led us to examine a range of effective-mass and jet observables in the simpler context of pure jet events (with the expectation that Z+j@NNLO might be similar to 2j@NLO). There we saw a significant NLO K-factor for all effective mass variables except one, H T,2 , which summed over just the two leading jets. In the Z+j case we had summed over all jets and hence it is not surprising that we should have observed substantialnNLO/NLO ratios.
Even though the observables in the pure jets case did not display giant K-factors, the pattern of NLO results was remarkably well reproduced atnLO. This encouraged us to studynNLO predictions, which provided substantial extra information for several of the observables, with the reference-observable method again giving important cross checks. Since the cross sections for the jet observables are large, these results could easily be tested with early LHC data.
We close this article with a few lines on the relation between LoopSim and other predictive methods. There is a close connection betweennLO (ornnLO) and CKKW and MLM [24,49] matching, since they also both provide ways of combining tree-level results with different multiplicities. Of course CKKW and MLM matching provide an interface with parton showers too, which the LoopSim method does not. On the other hand it is significantly easier to include multiple loop orders into the LoopSim method than it is within matrix-element/parton-showering matching procedures (though work is ongoing in this direction see e.g. [50]).
An interesting cross-check of the LoopSim method will come with the completion of the NNLO calculations for the Z+j and dijet processes. At that point the method could also, for example, be used to merge Z@NNLO with Z+j@NNLO, so as to provide annNNLO prediction for quantities like the Drell-Yan lepton p t spectrum. The value of the LoopSim method also goes hand-in-hand with progress on 1-loop calculations, especially with the prospect of automated NLO calculations now on the horizon (for example [51,52,53]).
Note that the LoopSim code, which will be made public in due course, can currently only deal with hadron-collider processes involving any number of light partons and up to one vector boson. It would benefit from further work to appropriately include heavy quarks and additional bosons.

A Recoil procedure
In this appendix, we provide further details on how we perform the recoil of an event when a particle becomes virtual, including the treatment of the decay products of the Z boson. We first examine the simpler case of a particle that makes a loop with the beam, then we show how to deal with a particle that makes a loop with another particle. A.1 A particle recombines with the beam Let us assume that particle i 0 makes a loop with the beam. To balance the event moment, we follow the following procedure: 1. For each particle i = i 0 , store its rapidity y i .
2. Perform a separate longitudinal boost on each particle so as to bring its rapidity to 0 (i.e. get a purely transverse event).

Compute
where E i is the energy of particle i in the purely transverse event.

Define
and boost all particles into the rest frame of k (so that the total transverse momentum balances).
5. Perform a longitudinal boost on each particle so that it recovers its original rapidity y i .
For the case where two particles, i 0 and i 1 , are looped with the beam, replace i = i 0 with i = i 0 , i 1 and in eq. (32) replace p t,i 0 with p t,i 0 + p t,i 1 , etc. In the case where the Z decays, for instance into 2 leptons, the procedure is identical except that we apply to the  Figure 13: Case where four gluons are emitted from the same quark. Gluon 1 is the last to be clustered with the quark (which roughly corresponds to an early time emission) and gluon 4 is the first to be clustered. In the event where gluon 2 makes a loop over the quark, we spread the gluon 2's momentum over the quark's momentum and the momenta of gluons that were emitted after it, i.e. gluons 3 and 4 (an earlier time emission like gluon 1 cannot be affected).
leptons the same longitudinal boosts as for the Z (the rapidity of the leptons is thus not necessarily 0 when we apply the transverse boost). This conserves the property that the sum of the leptons' momenta is still the Z momentum in the "looped" event.
The logic of the above procedure is that if we had attempted to apply a transverse boost without stages 2 and 5, we would have found that our choice of transverse boost, and the corresponding mapping of high-p t particles' momenta, would be affected by the presence of energetic particles collinear to the beam. This would have made the procedure collinear unsafe.

A.2 A particle recombines with another particle
Let us consider the situation depicted in fig. 13: four gluons are emitted from the same quark, but at different angles: and gluon 2 becomes virtual. The virtualisation of gluon 2 over the quark cannot have an impact on gluon 1, which was emitted earlier in an angular-ordered picture. But it has an impact on gluons 3 and 4. More precisely, let the p i be the momenta in the original event and p ′ i the momenta in the event where gluon 2 is virtual. We define and then set the p ′ i as follows: Subsequently each particle's p ′ i momentum is adjusted such that its mass is 0 (or m Z if gluon 3 is a Z boson rather than a gluon), keeping its transverse components p x , p y and its rapidity unchanged. This can be easily generalised to any number of particles recombining with the same hard one: for each recombined particle i, we spread the looped particle over the hard particle h and over any non-looped emissions from h that are at smaller angle (i.e. earlier in the C/A clustering sequence) than i. In eqs. (34,35a) it is always the original particle momenta that are used to determine the p t,i /p t,tot ratio, so that the result is independent of the order in which we perform the recombinations.
This procedure is designed to ensure collinear safety: if, for instance, gluon 4 is collinear to the quark in the original event, then it remains collinear in the looped event. And if it is gluon 4 (emitted after all the others) that is looped, only the quark momentum is rescaled and its direction barely changes, so that angles between the quark and the other gluons stays the same.
In the case where the Z decays into 2 leptons, one applies the following procedure to each of the leptons: 1. Perform a longitudinal boost of the Z boson respectively in the original event and the looped event such that it has 0 rapidity in each case. Call the momenta obtained p Z,0 = (E 0 , p t,0 , 0) and p Z,1 = (E 1 , p t,1 , 0) respectively.
2. Perform a longitudinal boost of the lepton from the original event into the frame where the initial Z has 0 rapidity.
3. Define a purely transverse vector k such that p Z,0 is transformed to p Z,1 if it is boosted into k's rest frame: with C = ( p t,1 − p t,0 ) 2 (E 1 + E 0 ) 2 .
4. Boost the lepton's momentum into k's rest frame.
5. Apply to the lepton the longitudinal boost that brings p Z,1 to its true rapidity in the looped event.
We are aware of the cumbersome nature of these procedures. A simplification of them that retained the relevant collinear-safety properties would certainly be of interest.

B The LoopSim method and incoming partons
Without going into a full proof, we shall here illustrate why the LoopSim method is sensible even in the presence of incoming hadrons, by considering what happens atnLO. We start with a LO cross section for a process producing n hard objects For compactness of notation, we have dropped the µ r dependence in the differential treelevel partonic cross section dσ ij→n /dΦ n . We have also not yet specified our choice for the factorisation scale µ f . We assume that dσ ij→n /dΦ n contains the necessary constraints to Next, we exchange i ↔ k, replace x ′ a → x a and change the scale µ 2 f in f i/a (x ′ a /z, µ 2 f ) to be Q 2 0 , which is allowed because it corresponds to an O(α 2 s ) change (while here we consider only O(α s )): dz z p ik (z) [C(p 1 , . . . , p n+1 ) − C(p 1 , . . . , p n )] f i/a (x a /z, Q 2 0 ) . (44) Note now that if we take µ 2 f ∼ Q 2 in eq. (40), then the second term in square brackets in eq. (44) cancels the first term in round brackets in the second line of eq. (40). In other words for initial-state radiation, the action of LoopSim is not so much to provide virtual corrections as to cancel the real-emission terms already included implicitly through the PDFs in the leading order cross section. In contrast, the true virtual terms are already included through the PDFs themselves, i.e. through the second term in round brackets in eq. (40).
As an example, consider pp → Z. AtnLO we will have events such as gq → Zq, where the outgoing quark comes from collinear initial-state splitting g → qq, with an underlying hard subprocessqq → Z. From these events LoopSim will generate a configuration in which the outgoing quark is "looped". This will come in with a PDF weight that is the product of a gluon distribution and a quark distribution, so it appears that we have a (negative) gq → Z contribution, which would be unphysical. However in the LO cross section with a factorisation scale µ f ∼ Q, when we writeqq → Z, part of theq PDF comes from g →qq splitting. If we were just to add the real gq → Zq diagram to the LO cross section alone, then in the collinear limit we would be double counting the part already included in the PDF. With the negative "gq → Z" LoopSim contribution, what happens is that we simply remove theq PDF component, generated from g →qq splitting, that was implicitly included at LO with an incorrect final state (i.e. lacking an outgoing quark), since we are now putting it in with the correct final state through the real gq → Zq diagram.
Note that we have not yet worked out the full extension of this discussion to higher orders. The details would depend on the precise higher orders that we have in mind, for examplennLO versusnNLO. However, regardless of these details, the fundamentally unitary nature of the LoopSim procedure is important in ensuring that the simulated "loops" simply bring about an overall consistent set of final states while maintaining the total cross section as calculated with a sensible factorisation scale choice.

C Secondary emitters in LoopSim
In section 2.1.2 we discussed the special treatment needed for "secondary emitters", i.e. non-Born particles that have emitted something. In our procedure, secondary emitters do not get looped: when particle j makes a loop over i, this is justified by the collinear enhancement of the matrix element due to i and j being close in angle. But the emission of the same j from the configuration where i is virtual does not have such a collinear enhancement, so one must not take it into account. Another way to understand it is to consider the example of 2-gluon emission from a qq dipole. The squared matrix element for the emission of 2 real energy-ordered (E 1 ≫ E 2 ) gluons, g 1 , g 2 , can be expressed as [54,55,56] M(k 1 , k 2 ) = (4πα s ) 2 (C 2 with W 1 = 4 (p q .pq) (p q .k 1 )(k 1 .pq) (p q .pq) (p q .k 2 )(k 2 .pq) , (46a) W 2 = 2 (p q .pq) (p q .k 1 )(k 1 .pq) (p q .k 1 ) (p q .k 2 )(k 2 .k 1 ) + (pq.k 1 ) (pq.k 2 )(k 2 .k 1 ) − (p q .pq) (p q .k 2 )(k 2 .pq) . (46b) Since the W 1 term diverges when g 2 is collinear to q orq (unlike the W 2 term), it becomes relevant when g 2 is considered to have been emitted from q orq independently of g 1 . The W 2 term diverges when g 2 is collinear to g 1 (unlike W 1 ), so it becomes relevant when g 2 is considered to have been emitted from g 1 . This is depicted in fig. 14, which also shows the virtual corrections (cf. [56]). One notices that the W 2 term only appears when g 1 is real. The diagrams where g 1 is virtual are taken into account when g 2 is emitted from q orq. Therefore, g 1 cannot become virtual when g 2 makes a loop over it.