Z/gamma* plus Multiple Hard Jets in High Energy Collisions

We present a description of the production of a di-lepton pair (through Z boson and virtual photon) in association with at least two jets. This calculation adds to the fixed-order accuracy the dominant logarithms in the limit of large partonic centre-of-mass energy to all orders in the strong coupling. This is achieved within the framework of High Energy Jets. This calculation is made possible by extending the high energy treatment to take into account the multiple t-channel exchanges arising from Z and gamma*-emissions off several quark lines. The correct description of the interference effects from the various t-channel exchanges requires an extension of the subtraction terms in the all-order calculation. We describe this construction and compare the resulting predictions to a number of recent analyses of LHC data. The description of a wide range of observables is good, and, as expected, stands out from other approaches in particular in the regions of large dijet invariant mass and large dijet rapidity spans.


Introduction
The Large Hadron Collider (LHC) sheds ever more light on Standard Model processes at higher energies as it continues into Run II. One "standard candle" process for the validation of the Standard Model description in this new energy regime is the production of a dilepton pair through an intermediate Z boson or photon, in association with (at least) two jets [1][2][3][4][5][6][7]. This final state can be entirely reconstructed from visible particles (in contrast to pp → dijets plus(W →)eν) making it a particularly clean channel for studying QCD radiation in the presence of a boson. Experimentally, this process is indistinguishable from the production of a virtual photon which has decayed into the same products, and we will consider both throughout.
W and Z/γ * -production are excellent benchmark processes for investigating QCD corrections, since the mass of the boson provides a perturbative scale, while the event rates allow for jet selection criteria similar to those applied in Higgs boson studies. W, Z/γ * -production in association with dijets is of particular interest, since in many respects it behaves like a dijet production emitting a weak boson (i.e. electroweak corrections to a QCD process rather than QCD corrections to a weak process). This observation means that a study of W, Z/γ * -production in association with dijets is relevant for understanding Higgs-boson production in association with dijets (which in the gluon-fusion channel can be viewed as a Higgs-boson correction to dijet production). This process is interesting (e.g. for CP -studies) in the region of phase space with large dijet invariant mass, where the coefficients in the perturbative series have logarithmically large contributions to all orders. As an example of the increasing importance of the higher orders, it is noted that the experimental measurement of the (N + 1)/N -jet rate in Z/γ * +jets increases from 0.2 to 0.3 after application of very modest VBF-style selection cuts even at 7 TeV [1,2,4].
The current state-of-the-art for fixed-order calculations for this process is the next-to-leading order calculation of Z/γ * plus 4 jets by the BlackHat collaboration [8]. While it has become standard to merge next-to-leading order QCD calculations with parton showers [9][10][11][12][13][14], results for jet production in association with Z/γ * bosons have so far only appeared with up to two jets [15,16] (corresponding results for a W boson with up to three jets were given in [17], following those for a W boson plus two jets in [16,18]). Indeed, W/Z + 0−, 1− and 2−jet NLO samples have been merged with higher-multiplicity tree-level matrix elements and parton shower formulations [19,20]. Beyond the matching, the parton shower cannot be expected to accurately provide a description of the large-invariant mass limit, from its resummation of the (soft and collinear) logarithms which are enhanced in the region of small invariant mass. An alternative method to describe the higher-order corrections is instead to sum the logarithmic corrections which are enhanced at large invariant mass between the particles. This is the approach pioneered by the High Energy Jets (HEJ) framework [21,22]. Here, the hard-scattering matrix elements for a given process are supplemented with the leading-logarithmic corrections (in s/t) at all orders in α s . This approach has been seen to give a good description of dijet and W plus dijet data at both the TeVatron [23] and the LHC [24][25][26][27][28]. In particular, these logarithmic corrections ensure a good description of W plus dijet-production in the region of large invariant mass between the two leading jets [28] and in large invariant mass regions in a recent 4-jet ATLAS study [29]. It is not surprising that standard methods struggle in the region of large invariant mass, since the perturbative coefficients receive large logarithmic corrections to all orders, and perturbative stability is guaranteed only once these are systematically summed.
The purpose of this paper is to develop the treatment of such large QCD perturbative corrections within High Energy Jets to include the process of Z/γ * plus dijets. While this process has many features in common with the W plus dijets process, one major difference is the importance of interference terms, both between different diagrams within the same subprocess (e.g. qQ → qQ(Z →)e + e − with emissions off either the q or Q line) and between Z and γ * processes of the same partonic configuration. For processes with two quark lines, the possibility to emit the Z/γ * from both of these leads to profound differences to the formalism, since the t-channel momentum exchanged between the two quark lines obviously differs depending on whether the boson emission is off line q or Q. Furthermore, the interference between the two resulting amplitudes necessitates a treatment at the amplitude-level. High Energy Jets is formulated at the amplitude-level, which, together with the matching to high-multiplicity matrix-elements, sets it apart in the field of high energy logarithms [30][31][32][33][34][35][36][37][38]. The added complication over the earlier High Energy Jets-formalism (and indeed in any BFKL-related study) by the interfering t-channels introduces a new structure of divergences in both real and virtual corrections, and therefore a new set of subtraction terms are needed, in order to organise the cancellation of these divergences. The matching to full highmultiplicity matrix elements puts the final result much closer to those of fixed order samples merged according to the shower formalism [15,16,19,20] -although of course the logarithms systematically controlled with High Energy Jets are different to those controlled in the parton shower formalism. In particular, High Energy Jets remains a partonic generator, i.e. although it is an all-order calculation (like a parton shower), it is not interfaced to a hadronisation model. Initial steps in combining the formalism of High Energy Jets and that of a parton shower (and hadronisation) were performed in Ref. [39].
We begin the main body of this article by outlining the construction of a High Energy Jets amplitude and its implementation in a fully flexible parton level Monte Carlo in the next section. In section 3 we derive the new subtraction terms which allows us to fully account for interference between the amplitudes. The subtraction terms allow for the construction of the all-order contribution to the process as an explicit phase-space integral over any number of emissions.
Specifically, the main result for the all-order summation is formulated in Eq. (28): where σ is the sought-after cross section, and the rest of the equation is discussed in the relevant section. Section 3 also discusses the necessary modifications in order to include fixed-order matching. In section 4 we show and discuss the comparisons between the new predictions obtained with High Energy Jets and LHC data. We conclude and present the outlook in section 5.

The High Energy Limit of QCD and Real Corrections
Fadin and Lipatov observed [30,31] that QCD scattering amplitudes at large invariant mass (compared to the transverse momenta involved) exhibit the scaling expected from Regge-theory.
In particular, this means that for a given configuration of the transverse momenta in a 2 → nscattering, the limiting behaviour of the scattering amplitude as the invariant mass between each pair of partons increases is dictated by the maximum spin of any particle, which could be exchanged in what is termed the t-channel between partons neighbouring in rapidity. This is found by ordering both initial and final state particles according to rapidity (or light-cone momenta in the case of incoming particles), and drawing all possible colour connections between these. If a colour octet connection is allowed between pairs of particles, this corresponds to the possibility of a spin-1 gluon exchange, whereas colour-singlet exchange is identified as a spin-1/2 quark exchange.
The contribution to the cross section from a given momentum configuration of the jets (as opposed to partons) from the different flavour assignments will have a different limiting behaviour, since the large invariant-mass scaling is different e.g. in the process of qg → qg, if the rapidity ordering of the final state q and g is swapped. Considering a specific transverse momentum configuration of the jets in a simple 2 → 2-process, the full amplitude (which will then be squared in the calculation of the cross section) will scale as s ω , where s is the invariant mass of the final jets and ω is the spin of the particle which would be exchanged in the t-channel. Some cases, e.g. gg → gg, always allow for a gluon to be exchanged, and hence the amplitude scales as s 1 for large s. In other cases, e.g. qg → qg, the t-channel particle exchanged is either a quark or a gluon depending on the rapidity order of the flavour assignment, and hence the amplitude scales as s 1/2 or s 1 for large s. However, in this case, it is clear that in the limit of large s the contribution to the resulting jet momentum configuration will be dominated by the process with the gluon exchange. This discussion is illustrated further in Fig. 1. This argument may be further generalised to the case of more than two outgoing partons, where now a 2 → n amplitude scales as where the outgoing particles are ordered in rapidity, s ij is the invariant mass of particles i and j and ω i is the spin of the particle exchanged in the t-channel of neighbouring particles. Γ({t i }) Figure 1: The two lines above illustrate the two possible rapidity orders for the process qg → qg.
In the first case, where the rapidity of the gluon is greater than the quark, the allowed colour connection is a singlet corresponding to a quark exchange in the t-channel. This leads to a contribution to the amplitude which scales as s 1/2 . In the second case, the allowed colour connection is an octet which corresponds to a gluon exchange in the t-channel and a scaling of s 1 . The latter will clearly be the dominant configuration in the limit of large s.
depends only on the square of the t-channel momenta (which in the limit corresponds to minus the square of their transverse components).
We have thus identified the flavour-assignments of partons which will yield the dominant contribution in the limit of large invariant mass between the jets, for any given configuration of the transverse momenta: the dominant contribution is obtained in the flavour configurations which allow for colour-octet (gluon) exchanges between all neighbouring particles. We call these "FKL configurations". Within High Energy Jets we concentrate on describing to all orders in the strong coupling these scattering amplitudes, which contribute to the leading power behaviour of the cross section.
These scaling arguments are unaffected by the additional emission of an electroweak boson and specifically here we discuss the description with an additional Z boson or virtual photon. The emission of an electroweak boson is viewed merely as an electroweak correction to the underlying QCD dijet production.
We begin by considering qg-initiated processes where the quark is the backward-moving incoming parton and take the leptonic decay of the Z/γ * . The ordering described above motivates a unique definition of t-channel momenta, namely if p a is the momentum of the backward quark, p b is the momentum of the forward gluon and y 1 y 2 ... y n , one then defines t i = q 2 i , where q 1 = p a − p 1 − p + − p − and q i = q i−1 − p i for 2 ≤ i ≤ n. Furthermore, the leading contribution, which satisfies the requirement of maximal t-channel gluon exchanges, arises purely from the outgoing state where all of the intermediate particles in rapidity (those labelled 2 to n − 1) must be gluons. As discussed later, the factorisation property of amplitudes in the high-energy limit then allows us to describe the emission of each of these gluons with an independent effective emission vertex, a generalised Lipatov vertex V µ [21], multiplying the corresponding expression for the equivalent 2 → 2 process, qg → qg(Z/γ * →) − + (see Fig. 2). At matrix-element-squared Figure 2: The schematic structure of the high-energy description of the matrix element for qg → qg...g(z →) − + , given in eq. (2). In that specific case particles a and 1 are quarks and particles b, 2,...,n are gluons. level this gives where The lowest order expression on the right-hand-side of Eq. (2), |M HE qg→Z/γ * qg | 2 , is the high-energy process, which will be described in full detail in section 2.3. While p a +p b = p − +p + +p 1 +p n for n > 2, the expression is built of two independent factorised pieces, so this is not a problem. Care needs to be taken with the expression for the t-channel pole, which must be taken symmetrically as 1/t 2 = 1/(t 1 t n−1 ). If the quark is instead the forward moving incoming parton, the expression is identical except for the definition of q 1 where the lepton momenta is removed.
For other initial states contributing to Z/γ * plus dijets, however, the situation is more complicated. In particular for qQ-initiated processes, as the Z/γ * may be emitted from either quark line, and there is interference from the two possibilities of exchanged t-channel momenta. The effective emission vertex remains valid, but we must now work at amplitude level to take into account this interference, both here and for the virtual corrections as described in section 3. In the remainder of this section we will develop the equivalent of eq. (2) for all channels of Z/γ * plus dijets. We begin this in the next subsection, by describing our method of constructing

Writing Matrix Elements in Terms of Currents
Traditionally, amplitudes in the HE limit are described as a product of two scalar "impact factors", one for each end of the t-channel chain. Instead, in HEJ, we describe the core 2 → X +2 processes in terms of a contraction of two independent currents. This is inspired by the structure of the exact tree-level amplitudes, where each quark line automatically generates a current. Effectively, helicity currents allow for the distinction of the kinematic invariants s and u, which is lost in the standard high-energy factorisation at the cross-section level. This distinction proves necessary in retaining accuracy in the approximations. This can already be illustrated in the simple example of qQ → qQ. For all negative helicities for example, one can immediately write: where we have employed the spinor-helicity notation for the quark spinors, where i|µ|j is shorthand forū − (p i ) γ µ u − (p j ). The repeated colour index d is summed over and the lower colour indices refer to their respective particle.
We will work in lightcone coordinates p ± = E ± p z and further define p ⊥ = p x + ip y and e iφ = p ⊥ /|p ⊥ |. In components, we get (using the spinors parametrised as in Ref. [21]) Let us first discuss the approach traditionally taken: in order to write this in the desired factorised form of a product of scalars, C(p a , p 1 ) × C(p b , p 2 ), it is necessary to use the limits p + This correctly captures the leading behaviour in s/t and gives a factorised expression.
However, by using helicity-currents, it is possible to achieve a form of factorisation without relying on kinematic approximations. Returning to eq. (4), it may immediately be written as a contraction of two factorised four-vectors: V (p a , p 1 ).V (p b , p 2 ), where the vectors depend on the same momenta as the factorised vertices in the traditional approach, but now the vectors (up to constants) are just standard currents j −µ (p i , p j ) = i|µ|j : Each helicity current has two independent components and this extra degree of freedom compared to the impact factors of the traditional approach is precisely what is required in order to keep the first term in eq. (5) and therefore describe the amplitude exactly.
This illustration is clearly for a very simple process, but the same conclusion applies more generally. One can exactly describe qg → qg as the contraction of a standard quark current and a gluon current j g µ , consisting of a product of a standard quark current and colour factors depending on the gluon momenta only [22]. This holds even though the qg-scattering process has s, t and u-singularities. The same holds for gg → gg as long as the helicities of the two incoming (and outgoing) gluons differ, such that one can define the s, t, u-channels. One can also go beyond pure QCD and describe qQ → W q Q, qQ → Z/γ * qQ and qQ → qQH exactly as the contraction of two currents [21]. In the next subsection we describe the new current for Z/γ * plus jets, and the construction of the resulting amplitude.

A Current for Z/γ * plus Jets
In this section, we will construct a current to describe the emission of a Z/γ * boson and exchange of a t-channel gluon from a quark or antiquark line. We can write the current for the Z emission (only), j µ Z , as a sum of the contributions from the two possible emission sites: one where the Z is emitted before the t-channel gluon and another where the gluon is radiated first, shown diagramatically in figure 3. For definiteness, we could then consider the decay Z → e + e − . We have where M Z is the mass of the Z boson, Γ Z is its width, C Zx is the coupling of the Z to x, x = e, q, ν e , . . . and µ is the Lorentz index for the t-channel gluon propagator. Expanding the quark and lepton momenta using their completeness relations we can fix the helicity of the incoming quark, h in , and the outgoing quark, h out , to be identical, and we are left with a current which only has four possible helicity configurations depending on h q = h in = h out and the = p in p out Figure 3: The current used to describe the quark line with the emission of a Z or γ * is the sum of the contributions arising from the two possible emission sites for the Z/γ * . electron helicity, h e : For the charged lepton channels for Z-decays, we must also include the contribution arising from the exchange of an off-shell photon, γ * . The expression for the current for the off-shell photon has the same form to that shown in eq. (9) with the Z propagator replaced with that of the photon and the couplings modified. Our final current is then the sum of the two:

All-Order Real Corrections for Z/γ * Plus Dijets
With the current derived in the previous subsection, we have the required building blocks to describe the dominant contribution to the real emission in the HE limit, in the manner of eq. (2). We first construct the lowest order description, |M HE qg→Zqg | 2 . Our current, j , is already the sum of diagrams with a mediating Z and diagrams with a mediating γ * . For the quark-gluon initiated processes, this is then all we need for the complete amplitude and we write: The interference term between the Z and γ * processes is immediately included in this construction through squaring the sum of Eq. (10). The equivalent expressions for the gq-initial state and for qg and gq-initial states all have the same simple form. This can then be substituted into eq. (2) to give the real corrections up to any order in α s .
We now turn our attention to the case of two incoming quark lines (or a mix of quark and anti-quarks). Here, it is possible for the Z to be emitted from either quark line, and it turns out that the interference effects are sizeable, see Fig. 4. We must include both possibilities and allow for the interference term. Our high-energy description of the matrix elements relies on the correct description of the t-channel momenta, and this obviously depends on which of the quark lines the Z or γ * was emitted from. We therefore need to modify the simple framework outlined above. We will use the subscript a (b) to label the current at the lowest (highest) end of the rapidity chain. We then define t a (t b ) to be the t-channel momentum exchanged when the bosons are emitted at the lowest (highest) end of the rapidity chain. Then the full amplitude squared for qQ → qQ(Z/γ * →)e + e − is given by: where j a,b are the pure quark currents defined above eq. (7). The coupling constants of the Z to the relevant quarks and leptons are contained within j Z/γ * (h q , h e ), as in eq. (8). Fig. 4 shows the value of this matrix element squared divided by the squared partonic centre-of-mass energy for increasing rapidity separation of the two jets. The result is compared with that obtained from the full, tree-level matrix elements from MadGraph5_aMC@NLO [14]. The slice through phase space here is given by: The matrix element squared divided byŝ 2 tends to a constant when the rapidity separation of the two outgoing partons grows large. This is as expected from BFKL and Regge theory. Fig. 4 also shows the separate contributions to the total matrix element squared coming from the Z/γ * emission from the forward moving quark line (black, dashed) and emission from the backward moving quark line (green, dotted). In this phase space slice, the leptons also have an increasing positive rapidity and so the forward emission matrix element describes the full matrix element most closely, with the contribution from backward-emission falling at large values of ∆y. The sum of the forward and backward emission matrix elements neglecting interference (magenta, dotted) significantly overestimates the final result. Once the (destructive) interference effects have been taken into account, the full sum (red, solid) correctly reproduces the LO matrix element (blue, thick solid). It is therefore clear that at low rapidities the inclusion of the interference effect plays an important role in the accuracy of the matrix element. Neither this effect nor the interference between the Z and γ * channels is included when electroweak corrections are included in a parton shower [41][42][43].
One can also investigate the importance of the virtual photon contributions we include and their interference with the pure Z process. The inclusion of the virtual photon terms is particularly important when studying a combined lepton invariant mass, (p e + + p e − ) 2 , far from the Z mass peak. This can be seen in Fig. 5, where slices through phase space are shown similarly to Fig. 4, but now for an (a) lower and (b) higher value of the dilepton mass. In both cases, the contribution of the virtual photon processes is above 25%.
Having established our description of the 2 → Z/γ * +2 parton process, we now turn our attention to adding the all-order real corrections. Our all-order expression will take the form of a sum of terms like eq. (2) for each of the three terms in eq. (12), such that the squared matrix element for qQ → (Z/γ * →)e + e − q(n − 2)gQ is: Leading Order |M HE | 2 forward emission |M HE | 2 backward emission |M HE | 2 sum (no interference) |M HE | 2 full (with interference) Figure 4: The matrix-element squared divided by the square of the partonic centre-of-mass energy for qQ → ZqQ with the Z decaying to an electron-positron pair for the phase space slice described in eq. (13). Increasing values of ∆ represent increasing rapidity separation between the jets. The different lines show the contributions from different terms in the calculation: only emission from the forward or the backward quark line (black, dashed and green, dotted), their sum without the interference term (magenta, dotted) and their sum including interference (red, solid) which is seen to agree exactly with the LO result (blue, thick solid). (p e + +p e − ) 2 =(120 GeV) 2 Tree Level -MadGraph Tree Level Z/γ * -HEJ Tree Level Z only -HEJ (b) Figure 5: The matrix-element squared divided by the square of the partonic centre-of-mass energy for qQ → Z/γ * qQ with the Z/γ * decaying to an electron-positron pair. The O(α 2 s α W ) tree-level contribution as described in HEJ (red, dashed) exactly matches that of Madgraph (blue, solid). The terms corresponding to the production of a Z boson only (green, dotted) significantly undershoots the full result. The virtual photon terms are, therefore, clearly an important contribution to the matrix element away from the Z Breit-Wigner peak.
In the case of n = 2, this reduces back to eq. (12). If either a or b is an incoming gluon, there is once again a unique set of t-channel momenta and one can set the relevant j to zero in the formula above. This then gives eq. (2) up to a factor of C A /C F which corrects the colour factor.
We therefore have a compact expression for the real-emission contribution to a given process at any order in α s . All real corrections can then be added by summing over n ≥ 2, provided that each contribution is finite. We will organise the cancellation of singularities using a phase-space slicing method which we describe in the next section.

Virtual Corrections and the Cancellation of Divergences
In the previous section, we derived a description for the dominant real emission corrections in the HE limit for a given process contributing to Z/γ * plus jets. Here we describe the corresponding virtual corrections and the organisation of the cancellation of divergences.
For a general QCD amplitude, the Lipatov Ansatz gives an elegant prescription for the leading logarithmic and next-to-leading logarithmic terms of the virtual corrections in the HE limit [30]. Each t-channel pole is supplemented with the following exponential factor: where q i⊥ is the transverse components of the relevant t-channel momentum and we have used dimensional regularisation with d = 4 + 2ε. Given the different 't's which enter the different terms of eq. (14), it is clear we must now also calculate the virtual corrections in three separate terms. We define ∆y i = y i+1 −y i and then incorporate the all-order virtual corrections as follows: To find the physical result (cross section, distributions, etc.), we now need to integrate over n-particle phase space and then sum over all n ≥ 2. However, before it is possible to do that, we must first organise the cancellation of divergences. There are two sources of divergences in eq. (16): the poles in ε within the virtual corrections and, upon integration over all phase space, the divergences which arise from any of the parton momenta going to zero. We do not have collinear singularities in our description, because by construction the particles are assumed to be well-separated.
We will use a phase space slicing method in which we divide the available phasespace into two regions by the introduction of a cut-off scale λ cut on p 2 ⊥ . Above the cut-off, we consider the emissions 'hard' and below the cut-off, we consider them to be 'soft'.
The divergence arising from the emission of a soft gluon can be seen directly from the effective vertex given in eq. (3). In the limit p 2 i⊥ → 0, we find Therefore, the effect of the i th emitted parton becoming soft at the level of the matrix element squared is: where the matrix element squared on the right-hand side is the corresponding one for the momentum configuration of the matrix element on the left-hand side after p i has been set to zero. The relation is identical if either q or Q is replaced by a gluon.
The integration over the soft phase space for the ith parton gives: 13 where we have used a change of variables from p z to rapidity. We will eventually go on to integrate over the momenta of all other particles, but the cancellation occurs already at the integrand level so we will not do so at this point. We have therefore found that the first-order correction to the qQ → Z/γ * q(n − 3)gQ process from this soft real emission is The corresponding first-order virtual correction is found by expanding the exponentials in eq. (16). We find We can now go through term-by-term to show the divergences cancel and find the resulting finite contribution to the matrix element squared. For the backward line Z/γ * emission squared terms, we have the following terms: Performing the expansion in of the final bracket yields: The poles in and the γ E terms have identically cancelled and we are left with a finite logarithm. This is a similar form to that found in [21,44]. The procedure for the forward line Z/γ * emission squared terms is identical and we find The cancellation for the interference terms is also similar and here we find as the finite remainder from the cancellation. These results are valid for any emission between the outer quarks/gluons which becomes soft. If either of the outer quarks/gluons becomes soft, this will also produce a divergence. To remain within the perturbative framework, we require that the outer particles are constituents of the jets and that their transverse momentum is above a minimum value.
It is clear that this result can be iterated order by order in α s . We would then form our final regulated all-order result as where we have defined One can easily check by expansion that this correctly reproduces the results in eqs. (23)- (25). However, the limit we have used from eq. (17) is a limit and not an exact identity. We therefore have to account for the difference between −V 2 (q i−1 , q i )/(t i−1 t i ) and its strict limit of 4/p 2 i⊥ for values of p i⊥ below λ cut . In practice, we include this correction for c cut < |p ⊥ | < λ cut with c cut = 0.2 GeV and find stable results around this value. We demonstrate that our numerical results are also insensitive to the precise value of λ cut in appendix A.
A total (differential) cross section can then be obtained by summing over all values of n and integrating over the full n-particle phase space, using an efficient Monte Carlo sampling algorithm [44,45]: where x a,b are the momentum fractions of the incoming partons and f f k (x k , Q k ) are the corresponding parton density functions for beam (k) and flavour f k . The factor ofŝ 2 is the usual phase space factor. The function Θ cut imposes any desired cuts on the final state. The minimum requirement is that the final state momenta cluster into at least two jets for the desired algorithm 2 .
In the regions of phase space where all final state particles are well separated in rapidity, this gives the dominant terms in QCD at all orders in α s (the leading logarithmic terms in s/t). However, in other areas of phase space, the differences due to the approximations used in |M HEJ−reg qQ→Z/γ * q(n−2)gQ | 2 will become more significant. We can therefore further improve upon eq. (28) by matching our results to fixed order results. Here, we match to high-multiplicity treelevel results obtained from Madgraph5_aMC@NLO [14] in two different ways. This amounts to merging tree-level samples of different orders according to the logarithmic prescription of HEJ.

Matching for FKL configurations
As described in section 2, these are the particle assignments and momentum configurations which contain the dominant leading-logarithmic terms in s/t. The first step of the HEJ description was to develop an approximation to the matrix element for these processes which was later supplemented with the finite correction which remained after cancelling the real and virtual divergences: |M HE qg→Zqg | 2 (eq. (11)) or |M HE qQ→ZqQ | 2 (eq. (14)). The approximation is necessary to allow us to describe the matrix element for any (and in particular, large) n and for including both the leading real and virtual corrections. However, if the parton momenta cluster into four or fewer jets 3 , the full tree-level matrix element remains calculable. In these cases, we perform the matching multiplicatively, so we multiply the integrand of eq. (28) by Here, {j i } are the jet momenta after a small amount of reshuffling. This is necessary because the evaluation of the tree-level matrix elements assumes that the jet momenta are both on-shell and have transverse momenta which sum to zero, neither of which is true in general for our events due to the presence of extra emissions. Our reshuffling algorithm [47] redistributes this extra transverse momentum in proportion to the size of the transverse momentum of each jet. The plus and minus light-cone components are then adjusted such that the jet is put on-shell and the rapidity remains unaltered. This last feature ensures that after reshuffling the event is still in an FKL configuration.
After this multiplicative matching factor has been included, the regularisation then proceeds as before.

Matching for non-FKL configurations
Away from regions in phase space where the quarks and gluons are well-separated, the non-FKL configurations will play a more significant rôle. These have so far not been accounted for at all, and hence we add three exclusive samples of leading-order two-jet, three-jet and four-jet leading-order events to our resummed events. The distinction between the samples is made following the choice of jet algorithm and parameters.
These two matching schemes complete our description of the production of Z/γ * with at least two jets, including the leading high-energy logarithms at all orders in α s . In the next two sections, we compare the predictions from this formalism to LHC data.

ATLAS -Z+Jets Measurements
We now compare the results of the formalism described in the previous sections to data. We begin with a recent ATLAS analysis of Z-plus-jets events from 7 TeV collisions [4]. We summarise the cuts in table 1. Any jet which failed the jet-lepton isolation cut was removed from the Lepton Cuts p T > 20 GeV, |η | < 2.5 ∆R + − > 0.2, 66 GeV ≤ m + − ≤ 116 GeV Jet Cuts (anti-k T , 0.4) p T j > 30 GeV, |y j | < 4.4 ∆R j > 0.5 event, but the event itself is kept provided there are a sufficient number of other jets present. Throughout, the central value of the HEJ predictions has been calculated with factorisation and renormalisation scales set to µ F = µ R = H T /2, and the theoretical uncertainty band has been determined by varying these independently by up to a factor of 2 in each direction (removing the corners where the relative ratio is greater than two). Also shown in the plots taken from the ATLAS paper are theory predictions from Alpgen [48], Sherpa [19,49], MC@NLO [9] and BlackHat+Sherpa [8,50]. We will also comment on the recent theory description of Ref. [20].
In Fig. 6 we begin this set of comparisons with predictions and measurements of the inclusive jet rates. HEJ and most of the other theory frameworks give a reasonable description of these rates. The MC@NLO prediction drops below the data because it only contains the hard-scattering matrix element for Z/γ * production and relies on a parton shower for additional emissions beyond the one hard jet. The HEJ predictions have a larger uncertainty band which largely arises from the use of leading-order results in the matching procedures.
We will now discuss a number of the differential distributions. In Ref. [4] these are displayed as distributions normalised to the inclusive Z/γ * -rate. However, given the excellent agreement between the HEJ-prediction and data for the inclusive 2-jet cross section, we prefer to compare to data directly the prediction obtained with HEJ for the distributions.
The size of the scale variation of the HEJ predictions is largely dictated by the matching to leading order accuracy. The smaller scale variation in the results of e.g. Blackhat+Sherpa is therefore a reflection of the benefit of going to NLO. The choice of not normalising the HEJ predictions further increases the size of the scale variation bands, as there is no cancellation in scale dependence in numerator and denominator. We find, though, that our scale dependence tends to lead to a change in overall normalisation rather than in shape. We demonstrate this by plotting (1/σ((Z/γ * → e + e − )+ ≥ 2j)) dσ/dX for various variables X in appendix B. Including such a normalisation factor significantly reduces the size of the scale uncertainty band, down to less than ±10% in both cases. The quality of agreement with the central line is unchanged.
The first differential distribution we consider here is the distribution of the invariant mass between the two hardest jets, Fig. 7. The region of large invariant mass is particularly important because this is a critical region for studies of vector boson fusion (VBF) processes in Higgs-plus-dijets, and as previously discussed, the corrections arising from QCD are similar in both processes: The radiation patterns are largely universal between these processes, so one can test the quality of theoretical descriptions in Z/γ * -plus-dijets and use these to inform the Hjj-analyses. It is also a distribution which will be studied to try to detect subtle signs of new physics. In this study, HEJ and the other approaches all give a good description of this variable out to 1 TeV. It will be interesting to see if the very good agreement between HEJ and the central data points will survive, once larger data sets lead to a reduction in the experimental uncertainty. The merged sample of Ref. [20] (Fig. 9 in that paper) combined with the Pythia8 parton shower performs reasonably well throughout the range with a few deviations of more than 20%, while that combined with Herwig++ deviates badly. In a recent ATLAS analysis of W -plus-dijet events [28], the equivalent distribution was extended out to 2 TeV and almost all of the theoretical predictions deviated significantly while the HEJ prediction remained flat. This is one region where the high-energy logarithms, included only in HEJ, are expected to become large.
In Fig. 8, we show the comparison of various theoretical predictions to the distribution of the absolute rapidity difference between the two leading jets. It is clear in the left plot that HEJ gives an excellent description of this distribution. This is to some extent expected as high-energy logarithms are associated with rapidity separations. However, this variable is only the rapidity separation between the two hardest jets which is often not representative of the total rapidity 'length' of events with more than two hard jets, since the hardest jets tend to be central in rapidity. Nonetheless, the HEJ description also performs well in this restricted scenario. The next-to-leading order (NLO) calculation of Blackhat+Sherpa also describes the distribution quite well while the other merged, fixed-order samples deviate from the data at larger values. The merged samples of Ref. [20] (Fig. 8 in that paper) describe this distribution well for small values of this variable up to about 3 units when combined with Herwig++ and for most of the range when combined with the Pythia8 parton shower, only deviating above 5 units.
The final distribution in this section is that of the ratio of the transverse momentum of the second hardest jet to the hardest jet. The perturbative description of HEJ does not contain   The scaling pattern is also investigated for a preselection selection of particles produced via vector boson fusion (VBF). Fi cross section as a function of the exclusive jet multiplicity and R jets with m jj > 350 GeV and | y jj | > 3.0, in the following referr The data are consistent with the BlackHat+SHERPA predicti multiplicity well whereas ALPGEN overestimates R 3/2 . Figure 6: These plots show the inclusive jet rates from (a) HEJ and (b) other theory descriptions and data [4]. HEJ events all contain at least two jets and do not contain matching for 5 jets and above, so these bins are not shown.
any systematic evolution of transverse momentum and this can be seen where its prediction undershoots the data at low values of p T 2 /p T 1 . However, for values of p T 2 0.5p T 1 , the ratio of the HEJ prediction to data is extremely close to 1. The fixed-order based predictions shown in Fig. 9 are all fairly flat above about 0.2, but the ratio to the data differs by about 10% for the Blackhat+Sherpa and Sherpa predictions. Clearly the theoretical uncertainties for the fixed-order based predictions for values of p ⊥2 /p ⊥1 close to 1 are very small. Comparing to the normalised distribution in appendix B, this is a region where the theoretical uncertainties in HEJ also become very small when normalisation is taken into account. constructed from opposite-sign leptons with p T > 20 GeV, |⌘| < 2.5, R``> 0.2 and 66 GeV  m`` 116 GeV and for jets with p jet T > 30 GeV, |y jet | < 4.4 and R`j > 0.5. Cross sections as a function of the inclusive and exclusive jet multiplicities and their ratios have been compared, as well as differential cross sections as a function of transverse momenta and rapidity of the jets, angular separation between the leading jets and the inclusive variables H T and S T . Compared with previous publications, the sensitivity has been extended to regimes with larger jet multiplicities and larger jet transverse momenta. In addition, the sample has been compared to theory in specific kinematic regions governed by large logarithmic corrections.
In general, the predictions of the matrix element plus parton shower generators and the fixed-order calculations are consistent with the measured values over a large kinematic range. MC@NLO fails to model not only higher jet multiplicities but also the transverse momentum of the leading jet. The transition from staircase to Poisson scaling of the exclusive jet multiplicity ratio, expected from theory when introducing a large scale difference, (b) Figure 7: These plots show the invariant mass between the leading and second-leading jet in p T . As in Fig. 6, predictions are shown from (a) HEJ and (b) other theory descriptions and data [4]. These studies will inform Higgs plus dijets analyses, where cuts are usually applied to select events with large m 12 .  Figure 11. (a) Measured cross section for Z (!``) + jets as a functio | y jj |, between the two leading jets and (b) as a function of the inva jets, m jj , for events with at least two jets with p jet T > 30 GeV and The cross sections are normalized to the inclusive Z (!``) cross sec in Figure 2. constructed from opposite-sign leptons with p T > 20 GeV, |⌘ 66 GeV  m`` 116 GeV and for jets with p jet T > 30 GeV, |y jet Cross sections as a function of the inclusive and exclusive ratios have been compared, as well as differential cross sections momenta and rapidity of the jets, angular separation betwee inclusive variables H T and S T . Compared with previous publ been extended to regimes with larger jet multiplicities and larg In addition, the sample has been compared to theory in specific by large logarithmic corrections.
In general, the predictions of the matrix element plus par the fixed-order calculations are consistent with the measured va range. MC@NLO fails to model not only higher jet multiplici momentum of the leading jet. The transition from staircase to P (b) Figure 8: The comparison of (a) HEJ and (b) other theoretical descriptions and data [4] to the distribution of the absolute rapidity different between the two leading jets. HEJ and Black-hat+Sherpa give the best description. These results will inform analyses of Higgs plus dijets, where cuts are usually applied to select events with large rapidity separation of jets. ized to the inclusive Z (!``) cross section. The azimuthal distance is well modelled by ALPGEN and by BlackHat+SHERPA. The tendencies observed in the modelling of the distance in and in rapidity are reflected in the measurement of the R spectrum of the leading jets. SHERPA models a too flat spectrum for both and R. The offset of 15% of the SHERPA prediction from the observed cross section in the bulk of the data in figures 11 and 12 is consistent with the results presented in figure 2(a) for the inclusive Z (!``) + 2 jets cross section.

Distributions after VBF preselection
A veto on a third jet is used to reject Z + jets background in selections of Higgs boson candidates produced by VBF. Figure 13 shows the transverse momentum and rapidity dis-

CMS -Z + Jets Measurements
We now compare to data from a CMS analysis of events with a Z/γ * boson produced in association with jets [5]. We show, for comparison, the plots from that analysis which contain theoretical predictions from Sherpa [19,49], Powheg [51] and MadGraph+Pythia [14]. The cuts used for this analysis are summarised in table 2. Table 2: Cuts applied to theory simulations in the CMS Z-plus-jets analysis results shown in Figs. 10-12. As in the previous section, any jet which failed the final jet-lepton isolation cut was removed from the event, but the event itself is kept provided there are a sufficient number of other jets present. The main difference to these cuts and those of ATLAS in the previous section is that the jets are required to be more central; |y| < 2.4 as opposed to |y| < 4.4. This allows less room for evolution in rapidity; however, as we will see, HEJ predictions are still relevant in this scenario. Once again, the central values are given by µ F = µ R = H T /2 with theoretical uncertainty bands determined by varying these independently by factors of two around this value. Once again, the theoretical uncertainty bands on the HEJ predictions are large (we note that they are not displayed in the MadGraph+Pythia6 predictions). The size is dictated by matching to leadingorder. As illustrated in appendix B, the scale variation effects are largely an overall normalisation and not a change in shape and are significantly reduced in normalised distributions. Therefore the agreement between the central predictions and data is more significant than the variation bands initially suggest. HEJ events always contain a minimum of two jets and therefore here we only compare to the distributions for an event sample with at least two jets or above.

Lepton Cuts
We begin in Fig. 10 by showing the inclusive jet rates for these cuts. The HEJ predictions give a good description, especially for the 2-and 3-jet inclusive rates in this narrower phase space. In Figs. 11-12, we show the transverse momentum distributions for the second and third jet respectively (the leading jet distribution was not given for inclusive dijet events). Beginning with the second jet in Fig. 11, we see that the HEJ predictions overshoot the data at large transverse momentum. In this region, the non-FKL matched components of the HEJ description become more important and these are not controlled by the high-energy resummation. The HEJ predictions are broadly similar to Powheg's Z-plus-one-jet NLO calculation matched with the Pythia parton shower. In contrast, Sherpa's central value significantly undershoots the data at large transverse momentum although it is within their scale variation band. Here the Madgraph+Pythia central prediction gives the best description of the data; their scale variation band is not shown.            Figure 12: The transverse momentum distribution of the third hardest jet in inclusive dijet events in [5], compared to (a) the predictions from HEJ and (b) the predictions from other theory descriptions.

Comparisons for the W ± +Jets/Z+Jets Ratio
In this section we briefly comment on the all-order predictions from HEJ for the ratio of W ± plus jets to Z/γ * plus jets events. We compare to data from a recent study undertaken by the ATLAS collaboration [6]. The cuts for both final states are summarised in table 3. Table 3: Cuts applied to theory simulations in the analysis of the ATLAS W ± +jets/Z+jets ratio predictions shown in Tables 4 and 5.

Lepton Cuts
Tables 4 and 5 show the measured values of the ratio between W -plus-jets and Z-plus-jets events, R jet , separated into inclusive and exclusive samples of 2, 3 and 4 jets. Also shown are the corresponding values from HEJ and the ratio of the two. We see extremely good agreement for the 2-jet ratios and the 3-and 4-jet ratios agree at the 10% level. This is comparable with the other theoretical predictions used in the study (BlackHat+SHERPA [8,50,52,53], ALPGEN [48] and SHERPA [19,49]) as can be seen in Fig. 1 of [6].

Conclusions
In this paper we have discussed augmenting the theoretical description of inclusive Z/γ * -plusdijets processes with the dominant logarithms in the High Energy limit at all orders in α s . In particular, the description constructed here is accurate to leading logarithm inŝ/t. This is achieved within the High Energy Jets (HEJ) framework. We began in section 2 by motivating and describing the construction of an approximation to the hard-scattering matrix element for an arbitrary number of gluons in the final state. This uses factorised currents for electroweak boson emission and outer jet production combined with a series of (gauge-invariant) effective vertices for extra QCD real emissions.
In contrast to previous HEJ constructions (for pure jets, W -plus-jets and Higgs boson-plusjets), the complete description of the interference contributions between Z and γ * processes and between forward and backward emissions required a new regularisation procedure. This is described in section 3 where we showed explicitly the cancellation of real and virtual divergences by using the Lipatov ansatz to include the dominant contributions in the High Energy limit of the all-order virtual contributions. The method by which we match our matrix element to the leading order matrix elements was also outlined here. In this way we achieve the formal accuracy of our Monte Carlo predictions to Leading Logarithmic in (ŝ/t) and merge Leading Order predictions in α s for the production of two, three or four jets.
In section 4, we compared the predictions of our construction to Z/γ * -plus-jets data collected at the ATLAS and CMS experiments during Run I. We see excellent agreement for a wide range of observables and can be seen to describe regions of phase space well where some other fixedorder-based predictions do not fare as well. Discrepancies which occur only do so in regions where we do not expect this description to perform as well, for example where there is a large ratio between p T 1 and p T 2 . We also discuss properties of other available theoretical descriptions. This all-order description of Z/γ * -plus-dijets allows predictions for the ratio of W ± +dijets to Z/γ * +dijets at all-orders in α s for the first time. This is an extremely important analysis as many theoretical and experimental uncertainties cancel in this ratio and in section 4.3, we show that we correctly reproduce the ratios of the total cross sections.
Just as for previous analyses of LHC data, it is found that the high-energy logarithms contained in HEJ are necessary for a satisfactory description of data in key regions of phases space, e.g. at large values of jet invariant mass. Such regions of phase space are crucial for the analysis of Higgs boson production in association with dijets. The impact of the high-energy logarithms will only be more pronounced at the larger centre-of-mass energy of LHC Run II, and beyond at a possible future circular collider. The HEJ framework and Monte Carlo is the unique flexible event generator to contain these corrections and will provide important theoretical input for the study of important processes at LHC Run II and beyond. on λ cut and would therefore dilute any dependence in the full sample. We begin in table 6 where we show the value of the cross section for different values of λ cut for exclusive 2-, 3-and 4-jet samples. The cuts applied are the same as in section 4.1. It is clear that the cross section does not display a large dependence on the value of λ cut . λ cut (GeV) σ(2j) (pb) σ(3j) (pb) σ(4j) (pb) 0. 2 5.03 ± 0.02 0.70 ± 0.02 0.13 ± 0.03 0. 5 5.05 ± 0.01 0.70 ± 0.01 0.13 ± 0.01 1.0 5.09 ± 0.01 0.71 ± 0.01 0.13 ± 0.01 2.0 5.16 ± 0.04 0.72 ± 0.01 0.13 ± 0.01 Table 6: The FKL-only cross sections for the 2-, 3-and 4-jet exclusive rates with associated statistical errors shown for different values of the regularisation parameter λ cut . The scale choice was half the sum over all transverse scales in the event, H T /2. Figure 13 shows the effect of the same variation in λ cut on the differential distribution in both the rapidity gap between the two leading jets in p ⊥ , ∆y j1,j2 , (a)-(c), and the rapidity gap between the two extremal jets in rapidity, ∆y jf,jb , (d)-(f). Results are shown for exclusive 2-, 3-and 4-jet samples in each case, once again the cuts applied are the same as in section 4.1. Again the scale choice for the central line was µ F = µ R = H T /2. The variation bands have been determined by varying these two scales independently by up to a factor of two in either direction with the extremal points removed where the relative difference between µ F and µ R is greater than a factor of 2. The distributions also show a very weak dependence on the choice of λ cut .
In practice, our default chosen value for λ cut is 0.2.

B Normalisation Effects on Scale Uncertainties in Z/γ * +Jets
Here we discuss the effect of normalising the predictions shown in section 4.1 to the total crosssection. We see from Fig. 6a that we describe the experimentally observed inclusive two jet rate very well and, as such, do not require normalisation to agree with the data. However, applying a normalisation procedure which consistently applies scale variation simultaneously in numerator and denominator significantly reduces the size of the scale uncertainty bands for High Energy Jets (or any theoretical prediction).
In Figs. 14a, 14b and 14c we show the results from Figs. 7a, 8a and 9a where we have normalised to the total cross-section calculated for each scale combination. We see that, as expected, the central value of HEJ still describes the data well in the regions discussed in section 4.1 and now the size of the theoretical uncertainty band is significantly reduced (by as much as a factor of 16 in the last bin of the p ⊥2 /p ⊥1 -distribution for example, and more typically by a factor of about 4). This illustrates that varying the renormalisation and factorisation scales leads to a change in overall normalisation but not to any significant change in shape. Therefore, it is still valuable to discuss the quality of agreement of the central line, despite their apparently large accompanying uncertainty bands in the unnormalised predictions.  Figure 13: (a)-(c) The effect of varying λ cut on the differential distribution in the rapidity gap between the two leading jets in p ⊥ , ∆y j1,j2 , with the N jet = 2, 3, 4 exclusive selections shown from left to right, and (d)-(f) for the rapidity gap between the most extremal jets in rapidity, ∆y jf,jb , with the N jet = 2, 3, 4 exclusive selections shown from left to right. The different colours represent λ cut = 0.2 (red), 0.5 (blue), 1.0 (green) and 2.0 (purple) and the bands represent the scale variation described in the text. (c) Figure 14: The predictions of section 4.1 normalised to the total cross-section, with scale variation consistently applied to numerator and denominator.