Higgs-boson plus dijets: higher-order matching for high-energy predictions

Several important processes and analyses at the LHC are sensitive to higher-order perturbative corrections beyond what can currently be calculated at fixed order. The formalism of High Energy Jets (HEJ) calculates the corrections systematically enhanced for a large ratio of the centre-of-mass energy to the transverse momentum of the observed jets. These effects are relevant in the analysis of e.g. Higgs-boson production in association with dijets within the cuts devised to enhance the contribution from Vector Boson Fusion (VBF). HEJ obtains an all-order approximation, based on logarithmic corrections which are matched to fixed-order results in the cases where these can be readily evaluated. In this paper we present an improved framework for the matching utilised in HEJ, which for merging of tree-level results is mathematically equivalent to the one used so far. However, by starting from events generated at fixed order and supplementing these with the all-order summation, it is computationally simpler to obtain matching to calculations of high multiplicity. We demonstrate that the impact of the higher-multiplicity matching on predictions is small for the gluon-fusion (GF) contribution of Higgs-boson production in association with dijets in the VBF-region, so perturbative stability against high-multiplicity matching has been achieved within HEJ. We match the improved HEJ prediction to the inclusive next-to-leading order (NLO) cross section and compare to pure NLO in the h → γγ channel with standard VBF cuts.


Introduction
Fixed-order perturbation theory delivers a good description of inclusive rates of colliderprocesses involving jets; however, logarithmic corrections of various origins are important for observables in different regions of phase space. For example, the detailed description of the dependence of the cross-section on jet-sizes R receives systematic logarithmic perturbative corrections of the type α n s ln n 1/R [1,2]. These logarithms are controlled by DGLAP-like evolution equations, which also govern the formalism of parton showers [3][4][5]. While corners of phase space characterised by large ratios of transverse scales are well described by the parton-shower formalism, measurements at D0 at 1.96 TeV [6] and ATLAS at 7 TeV [7] indicate clearly that even when matched with fixed-order matrix elements, the parton showers do not describe well the regions of large invariant mass or large rapidity spans of the jet systems. This region of phase space is of particular interest in the process pp → Hjj, with contributions (at Born level) of α 2 w through weak boson fusion and α 4 s through gluon-fusion (GF). It is reasonable to distinguish the two contributions to the same final state, since the quantum interference is negligible [8][9][10]. The impact of the radiative corrections to each process is rather different though; in particular, the t-channel colour octet exchange in the gluon-fusion process leads to increased jet-activity [11], which allows for a distinction of the production mechanism within the phase space populated by weak boson fusion. The two jets in weak boson fusion are often separated by a large JHEP08(2018)090 JHEP08(2018)090 The matching procedure described here can thus be viewed as merging the results of fixed-order calculations by use of the power-expanded matrix elements of HEJ coupled with the logarithmic virtual corrections, similarly to the CKKW-L-method [20,21] of using the logarithmic accuracy of a shower-algorithm to merge fixed-order cross sections of varying multiplicity. This paper will present a complete reformulation of the procedure for merging and all-order summation. With the same input (as in use of the same matrix elements to the same order), the results are unchanged, but the new procedure for obtaining the allorder results and the merging will allow for merging results beyond tree-level, and will be computationally much more efficient.
In section 2 we describe the original mechanism for matching leading-order samples within HEJ before a detailed discussion of the new formulation. This includes both analytical aspects and practical aspects of implementation. In section 3 we study the results obtained in the new formalism in the context of Higgs boson plus dijets in three studies. Firstly we confirm that when matching to fixed order samples is limited to a maximum of three jets, we find consistent results with the previous formalism. Secondly, we study the impact of increasing the multiplicity in the fixed-order samples, now possible for the first time. Thirdly, we compare the matched all-order results of HEJ with those obtained at next-to-leading order accuracy. In section 4, we conclude with a final discussion.

Matching
In the original formulation, the cross sections within HEJ are calculated by explicitly constructing the all-order result by first generating a 2 → n+l kinematic point for each number of partons n = 2, . . . , N , where N is chosen sufficiently large (in practice around 22), and l describes the non-partonic particles produced, e.g. Z, W, H or their decay products. In order to simplify the notation we will only discuss the purely partonic case. Likewise, we will restrict our discussion to the leading-logarithmic contribution. Note that all our arguments apply equally to the more general scenario. We demonstrate this by showing results for the production of a Higgs boson in association with at least two jets, including recently computed sub-leading corrections [18].
The high-energy limit is dominated by Fadin-Kuraev-Lipatov (FKL) configurations, where two partons scatter in such a way that there is no radiation outside the rapidity range spanned by the scattering partons and only gluons are emitted inside this range. To ensure the high-energy limit applied is valid, it is required that the extremal (in rapidity) partons are perturbative (hard in terms of transverse momentum), and are members of the extremal jets. The transverse momenta of the remaining partons are all generated down to effectively 0 GeV, technically to a very small scale of order 200 MeV (which can be varied), below which there is perfect cancellation between the subtraction terms (used in the organisation of the cancellation of the IR divergences [18]) and the real-emission terms. The matching to LO accuracy for all m-jet rates, m ≤ n, is then obtained by first projecting the kinematics of the generated all-order events into Born kinematics according to the number of hard jets as described in the previous section. The event weight is multiplied with a ratio of the square of the full Born-level matrix element to the HEJ approximation of the same.

JHEP08(2018)090
The cross section (and kinematic distributions) is then obtained through the formula:- where |M reg HEJ ({p i })| 2 is the square of the regularised all-order matrix element within HEJ for the 2 → n phase space point (see ref. [18] for further details), and is the ratio between the square of the matrix element evaluated at full tree-level accuracy and within HEJ for the state projected to tree-level 2 → m kinematics described by the jet is the exclusive m-jet measure applied to the generated event kinematics. T y indicates rapidity ordering. The limits of the integral over the transverse momentum of the extremal partons combined with the two-jet measure is set to guarantee that the extremal partons carry the dominant momentum of the extremal jets. We choose a cut-off p ⊥,min corresponding to 90% of the transverse momentum of the respective extremal jet. 1 We use here the phrase 'kinematics of the generated allorder event' to mean the n-parton kinematic point of the resummation event sampled in eq. (2.1). In order to match each m-jet rate to tree-level accuracy, each generated event in the all-order phase-space is mapped to a m-jet tree-level kinematic point, and requires an evaluation of the full m-jet matrix element.
The scale-variation of the normalisation of the cross sections is determined by the treelevel matrix elements, and mostly unchanged by the leading logarithmic (LL) high-energy resummation implemented in HEJ. This could be reduced by extending the reweighting factor w m−jet to next-to-leading order accuracy. However, in order to do this, one would have to integrate over all m+1 parton real emission phase space resulting in a specific m-jet Born level kinematics. This would be prohibitively time-consuming. An opposite approach is to begin with fixed order samples of exclusive jet rates and then merge these using HEJ to generate all-order results. We demonstrate how to do this in the next subsection and find significant benefits already at tree-level accuracy. In particular, each phase space point used for the tree-level matrix element maps into all the relevant resummation phase space points which leads to fewer evaluations of the tree-level matrix elements. This in turn allows for matching to higher multiplicity for a given CPU envelope.

Supplementing fixed order samples with HEJ resummation
The reformulation of the resummation and matching should reproduce the results of eq. (2.1). Starting from this equation, we introduce a δ-functional and an integration over the Born level kinematics of the on-shell, reshuffled jets {j i B } reconstructed from the resummed kinematics. Eq. (2.1) is rewritten to .
The first two lines are now the phase space integration over the LO matrix element, which can be represented in terms of (potentially weighted) tree-level events. Obviously, the Born-level partonic momenta are identical with the Born-level jet momenta, i.e.
only depends on the Born-level HEJ approximation to the matrix element. Lines 4-6 are the integration of the HEJ matrix elements over all of the resummation phase space, which map onto the given fixed-order kinematics. Finally, line 7 removes the factors introduced in the first line of eq. (2.3) compared to eq. (2.1) in order to write the matching in terms of a standard phase space integration over fixed-order PDFs and matrix elements. The δ-functionals of the fifth line in eq. (2.3) connect the reconstructed Born-level kinematics with the kinematics of the jets arising from the resummation. The algorithm devised for projecting the jet momenta of the resummation onto Born-level kinematic gives [19]  dσ/dp B ⊥,min [ab/GeV] Figure 1. The distribution of the minimum transverse momentum of jets used in the matching, p B ⊥,min , for Higgs boson plus dijet production of transverse momenta larger than 30 GeV. The fact that the distribution falls off quickly below the jet analysis scale ensures the resummation phase space with a given minimum jet transverse momentum is covered by a fixed-order generation with a slightly smaller requirement on the jet transverse momentum. For example, a generation of fixed-order events with a minimum jet transverse momentum of 20 GeV is sufficient for an analysis requiring a transverse momentum of at least 30 GeV. plus the constraint that the rapidities of the jets are kept fixed. Here p B J l is the momentum of the fixed-order, matching level jet, q ⊥ is the sum of the transverse momenta of partons outside jets after resummation, which equals minus the transverse momentum of the jets after resummation. P ⊥ is the scalar sum of the jet transverse momenta after resummation. This algorithm can be straightforwardly applied when the resummation event has been constructed, and had a jet-clustering applied. If, however, we want to start from fixed-order generated events, the algorithm needs to be inverted, such that all resummation-momenta on the right-hand side of eq. (2.5) are explored for a given Born-level kinematic point.
While eq. (2.3) is mathematically equivalent to eq. (2.1), it does not prove that the approach is viable. The first challenge is to ensure that in fact, the integration over the matching, or fixed-order phase space, in the first line of eq. (2.3) does not actually extend to zero transverse momentum of the matching jets. This would lead to a divergence in the fixed-order cross section and invalidate the starting point. In figure 1 we investigate the minimum transverse momentum of jets used in the matching for the evaluation of fixed-order matrix elements. The plot shows dσ/dp B ⊥,min , where p B ⊥,min is the minimum jet transverse momentum used in the merging with matrix elements (i.e. the minimum transverse momentum in the resulting on-shell Born level kinematics after reshuffling) for Higgs-boson production in association with at least two jets with transverse momentum of at least 30 GeV. One sees that the matrix element sample needs to include events with a JHEP08(2018)090 minimum jet transverse momentum below the final analysis scale -but not too far below. It is observed that this distribution gets broader, and that the weight for small p B ⊥,min is relatively more important, both for larger rapidity spans, and if more hard jets are required (obviously these two requirements are linked).
The next challenge is to generate all resummation kinematics corresponding to a specific fixed-order or matching kinematics. This is not an obvious switch to make: substituting a requirement on (N)LO kinematics to result in a given Born level jet configuration with that of the full resummation event resulting in a given Born level jet kinematics. However, the formalism will offer a number of benefits. Statistical convergence can be controlled at a more fine-grained level. Stability can be ensured first at the fixed-order stage before attempting resummation, and each jet multiplicity can be considered separately. We are free to choose whichever generators we find most suitable for producing fixed-order events. A further improvement is due to the fact that the fixed-order matrix element is evaluated only once for each fixed-order kinematic point, so we expect a significantly enhanced computational efficiency, especially for high jet multiplicities.

Phase space generation
In order to perform the resummation, we are tasked with the numerical evaluation of the last four lines of eq. (2.3). In principle, we have to integrate over the phase space of arbitrarily many further real emissions. This is made feasible by the fact that for a given fixed-order configuration with finite rapidity span, only a limited number of additional gluons actually lead to a non-negligible contribution in the resummation. Still, the typical multiplicities in the interesting region of large rapidity separations will be quite high and we are required to inspect the corresponding high-dimensional phase space carefully for an efficient integration. In the following, we discuss how to construct an efficient importance sampling.

Gluon multiplicity
The typical number of extra emissions depends strongly on the rapidity span of the underlying fixed-order event. Let us, for example, consider a fixed-order FKL-type multi-jet configuration with rapidities y j f , y j b of the most forward and backward jets, respectively. By construction of the matching algorithm of ref. [19], the jet multiplicity and the rapidity of each jet are conserved when adding resummation. This implies that additional hard radiation is restricted to rapidities y within a region y j b y y j f . Within HEJ, we require the most forward and most backward emissions to be hard in order to avoid divergences [19], so this constraint in fact applies to all additional radiation.
To simplify the remaining discussion, let us remove the FKL rapidity ordering where all rapidity integrals now cover a region which is approximately bounded by y j b and y j f . Each of the m jets has to contain at least one parton; selecting random emissions we JHEP08(2018)090 can rewrite the phase space integrals as [dp i ] (2.7) with jet selection functions and n g ≡ n − m. Here and in the following we use the short-hand notation [dp i ] to denote the phase-space measure for parton i. As is evident from eq. (2.7), adding an extra emission n g + 1 introduces a suppression factor 1 ng+1 . However, the additional phase space integral also results in an enhancement proportional to ∆y j f j b = y j f − y j b . This is a result of the rapidity-independence of the MRK limit of the integrand, consisting of the matrix elements divided by the flux factor. Indeed, we observe that the typical number of gluon emissions is to a good approximation proportional to the rapidity separation and the phase space integral is dominated by events with n g ≈ ∆y j f j b .
For the actual phase space sampling, we assume a Poisson distribution and extract the mean number of gluon emissions in different rapidity bins and fit the results to a linear function in ∆y j f j b , finding a coefficient of 0.975 for the inclusive production of a Higgs boson with two jets. In figures 2 and 3 we compare the fit with the actual outcome.

Number of gluons inside jets
For each of the n g gluon emissions we can split the phase-space integral into a (disconnected) region inside the jets and a remainder: We choose an importance sampling which is flat in the plane spanned by the azimuthal angle φ and the rapidity y. This is observed in BFKL and valid in the limit of Multi-Regge-Kinematics (MRK). Furthermore, we assume anti-k t jets, which cover an area of πR 2 [24]. In principle, the total accessible area in the y-φ plane is given by 2π∆y f b , where ∆y f b ≥ ∆y j f j b is the a priori unknown rapidity separation between the most forward and backward partons. In most cases the extremal jets consist of single partons, so that ∆y f b = ∆y j f j b . For the less common case of two partons forming a jet we observe a maximum distance of R between the constituents and the jet centre. In rare cases jets have more than two constituents. Empirically, they are always within a distance of 5 3 R to the centre of the jet [25], so ∆y f b ≤ ∆y j f j b + 10 3 R. In practice, the extremal partons are required to carry a large fraction of the jet transverse momentum (cf. section 2) and will therefore be much closer to the jet axis.
In summary, for sufficiently large rapidity separations we can use the approximation If there is no overlap between jets, the probability p J ,> for an extra gluon to end up inside a jet is then given by (cf. figure 4) For a very small rapidity separation, eq. (2.10) obviously overestimates the true probability. The maximum phase space covered by jets in the limit of a vanishing rapidity distance between all partons is 2mR∆y f b . We therefore estimate the probability for a parton to  end up inside a jet as (2.11) In figure 5 we compare this estimate with the actually observed fraction of additional emissions into jets. We observe good agreement over the whole rapidity range and for different jet multiplicities.

Gluons outside jets and observed jet momenta
Using our estimate for the probability of a gluon to be a jet constituent, we arrive at a number n g,J of gluons inside jets. Before integrating over their remaining phase space, we first have to determine the momenta p J i of the observed (resummation) jets from eq. (2.5).

JHEP08(2018)090
To this end, we have to determine the total transverse momentum q ⊥ of the gluons outside jets. After generating soft transverse momenta for these n g − n g,J gluons, we solve the nonlinear system eq. (2.5) using GSL routines [26]. Note that we have to postpone the rapidity integration, since at this point the rapidity span in the phase space integral is not yet known. The most forward and backward partons have to be part of the extremal jets. Therefore, their momenta will only be determined in the next step.

Gluons inside jets
Recall that after the first step in the phase space parametrisation, eq. (2.7), each jet has exactly one constituent. We now assign each of the n g,J gluons to a random jet. For jets with a single constituent, the parton momentum is fixed completely by the constraints in eq. (2.3). In the case of two constituents, we observe that the partons are always inside the jet cone with radius R and often very close to the jet centre. This allows an efficient integration by choosing a distance to the jet centre and an azimuthal angle with respect to the jet axis for one of the partons, which determines all momentum components of both constituents. As is evident from figure 5, jets with three or more constituents are rare and an efficient phase-space sampling is less important. For such jets, we exploit the observation that partons with a distance larger than R max = 5 3 R [25] to the jet centre are never clustered into the jet. Assuming N constituents, we choose distances, angles, and transverse momenta for N − 1 of them and determine the momentum of the last constituent from the requirement that the constituent momenta have to add up to the jet momentum. Since this last momentum may lie outside the jet cone, it is mandatory to check explicitly whether all candidates are actually clustered into the considered jet. This is to ensure the correct coverage of phase space.
After constructing the resummation jets, we are now in the position to evaluate the rapidity integrals for the partons outside the jets. Finally, we use fastjet [27] to recluster all emitted partons into jets again to check whether the reshuffling conditions imposed by eq. (2.3) are fulfilled. We also ensure that all partons are assigned as intended, i.e. the n g,J designated jet constituents are indeed part of their respective jet and all remaining partons end up outside jets.
We have now outlined the practical steps necessary to implement the rewritten formula of eq. (2.3). In the following section we discuss the results obtained with the new formalism in the key process of Higgs boson production in association with at least two jets. Firstly we confirm that if we limit ourselves to matching with fixed order samples with up to three jets that we reproduce the results obtained with the previous formalism, but now with a much higher efficiency. We will then show and discuss the impact of now being able to increase the multiplicity in the fixed order samples and also compare our results to fixed-order next-to-leading order predictions.

Results
The matching procedure described in this work is significantly more efficient and flexible than the approach used in previous versions of HEJ. To illustrate this, we present new JHEP08(2018)090 results for the production of a Higgs boson in association with at least two jets matched to leading-order events with up to four jets. Previously, matching of HEJ to just three jets was achieved for this process, while using significantly more CPU resources than necessary with the current approach. In its new formulation, the matching is in practice only limited by the capabilities of the underlying fixed-order generator. For instance, the generation of one set of 1000 unweighted leading-order events for the production of a Higgs boson with four jets typically took a few CPU days using MadGraph5 aMC@NLO [28]. It is then just a few additional CPU seconds to generate 100 weighted resummation events from each of the fixed-order 4-jet events, so 100 000 weighted trial resummation configurations in total. In the previous matching approach, generating 100 000 resummation configurations would require the same number of computationally expensive fixed-order matrix element evaluations. 2 Since the resummation is not followed by any computationally intensive steps, we only consider the generation of weighted events here. Nonetheless, we also observe a marked improvement in a short test simulation with unweighted events.
This section will present the results obtained with the new procedure for matching and resummation. Section 3.1 describes the cuts and analysis used. Section 3.2 compares new results with matching up to three jets with those obtained previously, and demonstrates that the two methods yield equivalent results. Section 3.3 investigates the stability of the results obtained by investigating the impact of increasing the order to which matching is achieved. In general, the matching to higher multiplicities should have little impact for configurations where the four-jet contribution is insignificant or the approximation within HEJ already provides a good description. Conversely, the corrections from matching to successive multiplicities can serve to indicate the stability of the HEJ predictions for observables sensitive to additional hard radiation. Finally, in section 3.4 we match the inclusive Hjj-cross section to NLO accuracy, thus obtaining the most precise predictions for Hjj-production, including the effects of VBF cuts and central jet vetos. These results are compared to those obtained at fixed next-to-leading order accuracy.

Setup
To facilitate the comparison with previous results we will adopt the cuts of the experimental analysis of ref. [29], and the parameters of our analysis in ref. [18]. To recapitulate, we consider the gluon-fusion-induced production of a Higgs boson together with at least two anti-k t jets with transverse momenta p ⊥,j > 30 GeV, rapidities |y j | < 4.4, and radii R = 0.4 at the 13 TeV LHC. While it is obviously irrelevant for the considerations of the QCD corrections considered in this paper, we consider the Higgs boson decay into two photons with |y γ | < 2.37, 105 GeV < m γ 1 γ 2 < 160 GeV, Γ H = 4.165 MeV and a branching fraction of 0.236% for the decay into two photons. We use the CT14nlo PDF set [30] as provided by LHAPDF6 [31]. In addition to inclusive quantities with the basic cuts listed above, we also consider additional VBF-selection cuts applied to the hardest jets as in [29]: In the first step, we generate leading-order events with two, three and four jets. With our new matching procedure we are free to use an arbitrary fixed-order event generator for this purpose. For the present analysis we employ version 2.5.5 of MadGraph5 aMC@NLO [28]. For each jet multiplicity we produce about 2000 sets of unweighted events, each comprising 10 000 events for the sets with two or three jets and 1000 events for sets with four jets.
As the transverse momenta of the jets are modified during resummation (cf. eq. (2.5)), we have to generate at least a fraction of events with Born-jet momenta below the threshold of 30 GeV required from the observed jets. As already shown in figure 1 the contribution after the resummation from such tree-level configurations in the matching drops off very rapidly below the jet transverse momentum analysis scale of 30 GeV. Passing this information to the underlying fixed-order generator, such that only a small fraction of events are generated below the nominal transverse momentum threshold could improve the sampling considerably. Having such an option would therefore be highly desirable. For the time being, we manually generate 200 additional sets of Born-level events with transverse momenta down to 20 GeV for each jet multiplicity.
Events with more exclusive jets than can be reasonably evaluated in MadGraph5 aMC@NLO are unmatched and generated with a custom built Monte Carlo generator based on tree-level HEJ matrix elements instead of full leading-order ones. In this way, we can supplement the fixed-order input with events including up to ten jets obtained within the HEJ approximation. These events are simply passed through the same matching mechanism based on eq. (2.3) just as the lower-multiplicity events obtained using MadGraph5 aMC@NLO. The maximum multiplicity of ten is an arbitrary cut-off, based on an explicit check that the impact on observables at this multiplicity is negligible.
Since the final kinematics required for a kinematic scale setting are not known at the point of generating fixed-order events, we use a fixed renormalisation and factorisation scale of µ r = µ f = m H during the fixed-order generation. After resummation the events are rescaled to a central scale of µ r = µ f = H T /2. In order to assess the scale dependence, we independently vary both the renormalisation and factorisation scales by factors {1/2, 1/ √ 2, 1, √ 2, 2} and discard combinations with µ r /µ f < 1/2 or µ r /µ f > 2. In the effective Higgs-gluon coupling, we keep the renormalisation scale at the Higgs-boson mass and apply the limit of an infinite top-quark mass. These scale settings and even the use of the infinite top-mass limit are however not inherent to the use of the high-energy resummation of HEJ, but can be included with modified components of the amplitudes, similar to ref. [32]. After generating the tree-level input events, we apply resummation as presented in the previous sections. Recent progress described in [18] allows us to apply resummation not just for FKL-ordered matching-events, but also the sub-leading contribution from events with three jets or more, where the rapidity-ordering of the two most forward or most backward jets is flipped compared to FKL ordering. This corresponds to a gluon emission outside a rapidity-interval delimited by quark jets. For each resummation-type tree-level event, we generate 100 weighted trial configurations in the resummation phase space. For the remaining sub-leading events we cannot add resummation and simply adjust the factorisation and renormalisation scales as described above.

Comparison to previous results
In order to demonstrate the validity of the new approach we compare here first our results with leading-order matching up to three jets to those obtained in our previous work [18]. We find good agreement within the statistical errors. As examples, we show the transverse momentum distributions for the Higgs boson for inclusive and VBF cuts in figure 6. The previous and new method of organising the calculation are equivalent. For the comparison, we have adjusted our settings to match those in [18] as closely as possible. Apart from restricting the fixed-order matching to configurations with at most three jets, this also means that the extremal partons are required to have a fixed minimum transverse momentum of 27 GeV instead of a fraction of the corresponding jet momentum, as discussed in section 2.

Impact of four-jet matching on distributions
The HEJ approximation is exact in the limit of Multi-Regge kinematics, i.e. for large rapidity separation between hard jets. An equivalent characterisation is to demand the centre-ofmass energy and the invariant masses between all final-state jets to be much larger than the typical transverse momenta of these. If these conditions are fulfilled, we expect a good HEJ prediction and hence small matching corrections. In order to assess the perturbative stability of the final predictions, we will here study the impact on the resummed and matched cross section of scale variations and of successive matching to two-jet, three-jet and four-jet tree-level events.
One of the main goals of HEJ is to improve the prediction of the gluon-fusion background to Higgs-boson production in weak-boson fusion. Standard VBF cuts project out a kinematic region with a large invariant mass between the hardest jets, where the gluon fusion receives significant contributions from higher jet multiplicities. Figure 7(a) displays the relative contribution of the exclusive two-, three-and four-jet component to the distribution on the invariant mass between the two hardest (in transverse momentum) jets. The relative contribution from exclusive three-and four-jet-events increases with increasing m j 1 j 2 . Figure 7(b) displays the impact of matching of successive multiplicity on the distribution of the invariant mass between the two hardest (in transverse momentum) jets. The effect of the four-jet matching is small but non-zero even at large m j 1 j 2 . This is because even in this limit a large separation between all jets is not guaranteed.
The contribution from jet multiplicities of more than or equal to 5 is less than 5% for an invariant mass of at least 1 TeV. We conclude that the uncertainty on the distribution of m j 1 j 2 from terminating the matching at the four-jet contribution is insignificant, and well within the quoted scale variation.
A central prediction of BFKL, which arises also within HEJ, is a linear increase in the number of jets for a growing rapidity span between the most backward and forward jets. 3 This behaviour is demonstrated in figure 8, which also investigates the impact of matching to tree-level of successive multiplicities. Although the contribution from higher jet multiplicities increases with the rapidity separation, the effect of fixed-order matching on this observable actually decreases. This confirms our expectation that the HEJ approximation JHEP08(2018)090   Figure 8. Average number of jets for fixed-order matching up to two, three, and four jets. In (a) we show the average total number of jets vs. the maximum rapidity-separation. In (b) we show the number of jets in the rapidity region of the two hardest jets.
works well for large ∆y j f ,j b . It is this linear increase in the average number of jets versus increasing rapidity span which can be exploited to suppress the gluon-fusion contribution with a central jet veto.
In contrast to this, if the two hardest jets are tagged, and only jets in-between these are counted as a function of the rapidity difference between the hardest jets, then the initial linear growth stalls at an average number of jets of around 2.3. The difference in behaviour to the VBF contribution is therefore less pronounced by tagging the hardest jets, rather than the most forward and backward hard jets. This was investigated further in ref. [33]. Also, the impact of the matching corrections remains sizeable for all rapidity separations.
In observables which are neither dominated by higher jet multiplicities nor completely described by the HEJ approximation we observe that the matching corrections are converging, but the corrections from four-jet matching are non-negligible. In figure 9 we show the distribution of the Higgs-boson transverse momentum with inclusive and with VBF cuts. While there is a notable difference between the matching to fixed-order predictions up to two and three jets, the effect of four-jet matching is much smaller. In all cases the matching corrections are well inside the scale variation.
The azimuthal angle between jets is of particular interest for the extraction of the CP -properties of the effective coupling between the Higgs boson and gluons. Figure 10 shows the effects of fixed-order matching on the distribution of the angle between the two hardest jets. Similar to the transverse momentum distribution in figure 9 the corrections from four-jet matching are uniformly moderate.
In order to achieve a greater reduction of the gluon-fusion background to weak-boson fusion within the VBF-cuts, a veto on further jets can be applied. This has the added benefit of reducing the contribution from higher jet multiplicities, which is harder to predict in perturbation theory. The effectiveness of such a cut relies on the difference in the quantum corrections to the processes of VBF and GF [11]. Since this difference is due to the t-channel colour-octet exchange of the GF process, we will apply a central jet veto only  in the regions away from the tagging jets, since the collinear regions have similar emissions in VBF and GF. This is a slight improvement on the normal central jet veto cuts, and is inspired by the Zeppenfeld variable [34]. Here, we consider a veto of events with jets within a rapidity distance y c to the rapidity centre of either (a) most forward and backward jets or (b) the hardest jets (see also [33,34]). In case (b), we only consider vetoing on further jets which are in between the two hardest jets. The results are shown in figure 11. As expected, the cross section in case (a) converges for large y c to the exclusive prediction for the production of a Higgs boson with exactly two jets, irrespective of the fixed-order matching to higher multiplicities. When applying the jet veto between the hardest jets in case (b), the overall reduction in the GF component is considerably smaller.  Figure 11. Effect of a jet veto between (a) the most forward and backward jets and (b) the two hardest jets. Events with additional jets within a distance of y c to the rapidity centre are discarded.

Matching and comparison to fixed next-to-leading order
The complete reformulation of the formalism for matching and all-order summation described in section 2 has allowed for matching to higher jet multiplicities in HEJ. The impact of the four-jet matching on the studied distributions is small. The method presented in the earlier sections has been concerned with a point-by-point matching of the resummation to full high-multiplicity tree-level accuracy. As extensively demonstrated in section 3.3, this achieves perturbatively stable results for the shapes of distributions. In order to reduce the scale variation and benefit from the full NLO results for Hjj-production, we will now rescale the results for HEJ within the inclusive cuts of eq. (3.1) to the NLO cross section for each choice of factorisation and renormalisation scale. Thereby, full NLO accuracy is obtained for all dijet observables, LO accuracy for trijet observables, and the impact on the shape of distributions from four-jet contributions is accounted for at LO. This method was applied also in ref. [33]. While this approach does not change the shape of distributions, the scale variation is reduced to the level of NLO predictions. We will here compare these predictions to those obtained at fixed NLO using MCFM [35,36] and SHERPA [5]. Figure 12 compares the predictions for the distribution on the invariant mass between the two hardest jets. The scale variation on the HEJ results is vastly reduced to that of figure 7, as generally expected by the inclusion of the full NLO corrections. The distribution obtained with HEJ for the invariant mass between the two hardest (in transverse momentum) jets is still significantly steeper than that at pure NLO, as a result of the possibility of significantly higher jet multiplicity, and the fact that hard central jets have a slightly smaller PDF-suppression than hard forward jets, and therefore the two hardest jets tend to also be central. This means that the prediction for the cross section within the VBF cuts is significantly smaller with HEJ than for NLO, and indeed lies outside the scale-variation band obtained at NLO. In numbers, the cross sections obtained (at NLO) for pp → h(→ γγ)jj for inclusive cuts and with a central scale choice of µ r = µ f = H T /2 is 6.58 +0.08 −0.57 fb. This is obviously the same as that obtained with HEJ

JHEP08(2018)090
for inclusive cuts, once the cross sections are normalised to NLO accuracy. For the VBF cuts, the NLO cross section is 0.872 +0.024 −0.090 fb, and that obtained for HEJ is 0.561 +0.031 −0.067 fb. Even though the inclusive cross section for HEJ is normalised to that obtained at NLO, a sizeable difference in the cross section within the VBF cuts arises due to a difference in the slope of distribution in m j 1 j 2 and the requirement of m j 1 j 2 > 400 GeV for the VBF cuts. The VBF cuts cause a similar reduction in the cross section to 13.2% (NLO) and 8.5% (HEJ) of the inclusive cross section respectively.
Comparing the results of figure 12 and figure 13 we observe that a choice of a central scale for the NLO calculation of µ r = µ f = H T /2 leads to a suspiciously small scale variation -and indeed the central scale choice gives results close to the extremum obtained with the variations, despite the scales being varied either side of the central choice of µ r = µ f = H T /2. Such a behaviour of the scale variation often indicates that the NLO scale variation obtained with this scale choice is underestimating the theoretical uncertainty [37].
Indeed, ref. [37] investigated the distribution in m j 1 j 2 for dijet production at NNLO at the LHC, and found that at large m j 1 j 2 this scale choice is favoured over p T based on perturbative convergence. The invariant mass between the two hardest jets obviously is not a stable perturbative scale choice for all bins in the distribution, which extends to very low values of m j 1 j 2 . With a central scale choice of µ f = µ r = max(m h , m j 1 j 2 ), the central scale choice leads to predictions in the centre of the variation band. The scale variation bands obtained with NLO and HEJ also overlaps in each bin of the distribution. With this central scale choice, the cross sections obtained at NLO for inclusive cuts is 6.23 +1.11 −1.22 fb. For the VBF cuts, the NLO cross section is 0.542 +0.156 −0.125 fb, and that obtained for HEJ is 0.359 +0.045 −0.061 fb. The VBF cuts cause a similar reduction in the cross section to 8.7% (NLO) and 5.8% (HEJ) of the inclusive cross section respectively.
It is worth noting that (ignoring the mass of each jet) since m 2 j 1 j 2 = 2p ⊥j 1 p ⊥j 2 (cosh(y j 1 − y j 2 )−cos(φ j 1 −φ j 2 ))a central scale choice of µ r = m j 1 j 2 systematically runs α s such that α s ∆y j 1 j 2 tends to a constant for large ∆y j 1 j 2 . This would seem to spoil the standard argument of BFKL noting large and systematic leading logarithmic corrections of the form (α s ∆y j f j b ) k at large ∆y j f j b , at least for ∆y j 1 j 2 sufficiently large that m j 1 j 2 is close to the hadronic collision energy that only two jets exists, because hard radiation beyond the two jets required is suppressed. For events with more than two jets, there is no direct correlation between ∆y j 1 j 2 and ∆y j f j b . The results for the scale choice of µ r = µ f = m j 1 j 2 are discussed further in appendix A. Here we just note that the apparent convergence of the perturbative series (i.e. a comparison of the LO and NLO results and scale variation) is not significantly different for the two scale choices. The scale variation around µ r = µ f = H T /2 is accidentally small, since the central scale choice leads to the maximum cross section within the variation.

JHEP08(2018)090
The results for HEJ are identical to those for 4-jet matching in figure 8 (since just the total cross section has been adjusted to the NLO result for Hjj), but the results are here compared to those obtained using the NLO calculation for Hjj-production. As observed also in previous analyses [38], the results obtained at NLO tends towards 2.5, where the exclusive, hard three-jet cross section is as large as the two-jet cross section. This clearly illustrates the slow convergence of the perturbative series. The results for NLO and for HEJ start diverging already at small ∆y j f j b . It is worth noting that the linear growth in the number of hard jets vs. ∆y j f j b has been experimentally confirmed [6,39] for several processes with colour octet exchanges in the t-channel.
Even though exactly the same events are involved, the breakdown of the convergence is less obvious in figure 14(b). The number of jets in-between the two hardest is obviously smaller, and both the results for HEJ and for NLO appear to asymptote to a value for the average number of jets of 2.2 for NLO and 2.3 for HEJ. Figure 14 also shows the predictions for the Higgs transverse momentum spectrum obtained at NLO and with HEJ both for inclusive (c) and VBF-cuts (d). The distributions are very similar for inclusive cuts, with a peak around 80 GeV, and the spectrum from HEJ slightly harder. For VBF cuts, the prediction for HEJ is lower than that for NLO, as a result of the steeper spectrum in m j 1 j 2 and the requirement of m j 1 j 2 > 400 GeV. The two predictions for the high-p ⊥ tail within the VBF cuts coincide, but in this region the infinite top-mass approximation is certainly not trustworthy.
Finally, figure 14(c) and (d) compares the azimuthal angle between the two hardest jets for (c) inclusive and (d) VBF cuts respectively. In both the distributions, the region of back-to-back jets at φ = ±π is slightly suppressed in HEJ compared to NLO. The valley at φ = 0 for the inclusive cut is due to the jet-algorithm removing the collinear region. The result within the VBF cuts on figure 14(f) show that the reduction in the cross section within the VBF cuts for HEJ compared to NLO predominantly is in the region of back-toback jets and jets in the same azimuthal direction. The region at φ j 1 j 2 = 0 is not collinear within the VBF cuts, and so the structure induced by the jet algorithm within the inclusive cuts of figure 14(e) is not present within the VBF cuts of figure 14(f).

Conclusions
We have presented a reformulation of the matching formalism within HEJ, which recasts the calculation as one of merging fixed-order samples of increasing multiplicity. The merging is performed respecting the resummation of perturbative terms logarithmically enhanced at largeŝ/p 2 t . While the formalism is mathematically equivalent to that previously used, stable results are obtained using orders of magnitudes less CPU time. This allows matching to be performed to higher multiplicity.
The new formalism was used in a study of Higgs-boson production in association with dijets. The impact of the higher-multiplicity merging is minimal on the shape of distributions important for the application of VBF cuts. For a central scale choice of µ f = µ r = H T /2, the VBF cuts reduce the inclusive cross section of 6.58 +0.08 −0.57 fb on the h → γγ-channel to 13% (0.872 +0.024 −0.090 fb) at NLO, or 8.5% (0.561 +0.031 −0.067 fb) once both NLO JHEP08(2018)090 and the HEJ-corrections are accounted for. The further suppression within HEJ is due to a steeper falling spectrum in the invariant mass between the two hardest jets. The NLO scale dependence is estimated by varying the renormalisation and factorisation scale independently by a factor of two. However, the scale variation around H T /2 is artificially small, since the central scale choice achieves a value close to the maximum within the variations. With a scale choice of µ r = µ f = m j 1 j 2 (but bounded from below by m h ), the spectrum is similar at NLO and with the further HEJ-corrections. With the scale choice, the inclusive NLO cross section for h(→ γγ)jj-cross section is 6.23 +1.11 −1.22 fb, and 0.542 +0.156 −0.125 fb within the VBF cuts. The result for HEJ within the VBF cuts is 0.359 +0.045 −0.061 fb. The VBF cuts cause a similar reduction in the cross section to 8.7% (NLO) and 5.8% (HEJ) of the inclusive cross section respectively.
The formalism presented in this report will be instrumental in the further developments, including an account of heavy quark mass effects, and mergings of NLO cross sections.
A Result for the central scale choice of m j 1 j 2 Figure 15 shows the same distributions for NLO and HEJ (with the inclusive cross section scaled to that of NLO for each scale choice) as those investigated on figure 14, but for predictions obtained using a central scale of µ f = µ r = max(m h , m j 1 j 2 ). Panel (a) shows the average number of hard jets vs. the rapidity difference between the forward-backward jet pair. The result for NLO is similar to that obtained with the scale H T /2, but the NLO scale variation is reduced (even though the scale variation on the NLO cross section themselves is increased by using m j 1 j 2 instead of H T /2). The rise in the number of hard jets is slightly stronger for HEJ than with the scale of H T /2. Figure 15(b) shows the average number of jets, counting only additional jets if their rapidity is in-between that of the two hardest jets. The behaviour of the NLO prediction for small ∆y j 1 j 2 is similar to that displayed on figure 14(b) until ∆y j 1 j 2 ∼ 3, after which the average number of jets decreases. This behaviour is a result of the decreasing value of α s for large ∆y j 1 j 2 with this scale choice. The

JHEP08(2018)090
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.