Finite quark-mass effects in Higgs boson production with dijets at large energies

The production of a Higgs boson in association with at least two jets receives contributions both from the fusion of weak vector bosons (VBF) and from QCD processes, especially gluon fusion (GF). The former process is important for measuring the coupling of the Higgs boson to weak bosons, whereas the latter process plays an important role in determining any CP-admixtures in the Higgs sector. In this paper we go beyond the current state-of-the-art for fixed order calculations of the GF process (i.e. one loop H + 2j including full quark mass effects) by including the all-order effects in leading log(ŝ/pt2), together with full quark mass and loop-propagator kinematic effects. We calculate the mass-dependent components and implement the resummation within the framework of High Energy Jets. The high-energy effects suppress the prediction compared to fixed order at large Δy12 and mjj (and therefore within the usual VBF cuts of widely separated jets), just as found in the limit of mt → ∞. The mass dependence is more significant than at fixed order, because the systematic inclusion of the leading logarithms in ŝ/pt2 results in a hardening of the transverse momentum of the Higgs boson, which in turn probes in more detail the loop-structure of the coupling. In particular, the full mass dependence reduces the cross section within VBF cuts by 11% compared to a calculation based just on the infinite top mass limit, but the impact of the bottom quark remains small. This all implies that the gluon-fusion contribution within VBF-cuts is less severe than current estimates suggest.


JHEP04(2019)127
In this paper we calculate the gluon fusion component of Higgs boson production in association with two jets supplementing the fixed order results with the leading logarithmic corrections inŝ/p 2 t to all orders in the coupling and including full quark mass effects. The necessary mass-dependent components for the resummation are derived, and allow for a first calculation of the interference of top and bottom quark mass contributions in H + 2j. Within VBF cuts, NLO predictions exhibit an instability which is cured by the inclusion of the all-order high-energy logarithms. Precise, stable predictions of the gluon fusion component in H + 2j are essential for the analysis of the Higgs boson couplings both to weak bosons and to gluons mediated by heavy quarks.
Let us first review the current status of the calculation of the quark-mass effects in Higgs Boson production. Higgs bosons are most copiously produced at the CERN LHC through heavy-quark mediated gluon fusion (GF), where the Born-level process is at oneloop and at order α 2 s . Even after the discovery [1,2] of the Higgs boson, an accurate prediction of the production mechanism is needed to reveal any possible effects on the production rate from the much sought-after physics beyond the Standard Model. Since the coupling of the Higgs boson to a quark is proportional to the quark mass, the gluon fusion process is dominated by the contribution from a top-quark loop. This process is known to N 3 LO (to order α 5 s ) in the limit of infinite top-mass [3][4][5][6]. Finite top-mass effects in the inclusive cross section can be taken into account at one order lower (to N 2 LO in α s ) by a formal expansion in (m h /m t ). The effects on the total cross section of the finite top-mass are found to be very small indeed [7][8][9]. An explicit calculation of the loop contribution using the full propagator dependence allows the inclusion also for the contribution from bottom-quarks. It is here found that the bottom-top interference effects are of the order of −5% [5].
Higgs boson production in association with one jet obviously forms a subset of the higher order corrections to the calculation of inclusive Higgs-boson production. As such, it is known to N 2 LO in α s (to order α 5 s ) in the limit of infinite top-mass. It was very recently calculated to NLO with full dependence on the heavy-quark propagator [10]. This explicit calculation of the quark loops can be used to check the earlier reported top-bottom interference effects, which were approximated using a small-mass expansion for the amplitudes involving bottom quarks, and the infinite mass limit for the top-quark contribution [11]. For H + j-production at NLO, the effect of the full dependence on the heavy-quark propagator momentum is a 9% increase in the overall cross section over the result obtained in the infinite top-mass limit. The quality of the approximations obtained using the infinite top mass limit in H + j-production is therefore much worse than for inclusive Higgs boson production. Furthermore, there is a strong phase-space dependence: for transverse momenta of the Higgs boson larger than 800 GeV, the effects of the full dependence on the heavy-quark propagators leads to a suppression over the result for an infinite top mass of more than an order of magnitude. For processes with more than just the Higgs boson in the final state, the limit of infinite top-mass loses not only the dependence on the mass of the propagating quarks, but also the full kinematic dependence on the propagators in the loop-diagrams.

Quark mass effects in Higgs boson production with HEJ
The HEJ framework is a perturbative framework for QCD processes which achieves leading logarithmic accuracy in variables which scale as s/t. These are seen to arise at all orders in α s from BFKL [36][37][38][39]. Logarithmic accuracy at the integrated level is obtained through a power series expansion of the relevant matrix elements in s/p 2 ⊥ . This gives the dominant terms in the multi-Regge kinematic (MRK) limit (or high-energy limit), defined for a 2 → n QCD process as s ij → ∞, |p i⊥ | finite i, j = 1, . . . , n. (2.1) In this limit, only a subset of flavour and momentum configurations contribute at leading power. From Regge theory, these are the configurations which permit the maximum number of t-channel gluon exchanges in a ladder chain with particles ordered in rapidity. We will call these "FKL configurations". 2 For example, in pure jet production qQ → qggQ and gg → gg are FKL configurations while qq → gg and qg → qQQg are not (all particles written in rapidity order). Furthermore, it is known that scattering amplitudes in the high-energy limit factorise into products of independent scalar factors and emission vertices, which each depend on a reduced set of momenta [28,40]. This is critical for obtaining simple expressions for large n. In order to achieve this simple form though, many approximations are required which limit the power of the description away from the strict limit, i.e. in physical regions of phase space. Within HEJ we achieved the t-channel factorised structure while making fewer approximations through the use of vector currents in place of scalar factors [29]. This immediately renders the HEJ description of 2 → 2 processes exact and preserves the correct position of the t-channel poles, and significantly improves the description of the matrix elements in the physical region.
The production of a Higgs boson with at least two jets using currents in HEJ has been described in the infinite quark mass limit in [34]. The t-channel factorised structure of the amplitudes lends itself to the inclusion of finite quark mass effects because the diagrams and loops do not become any more complicated than those at leading order for arbitrarily large n. This has been a bottleneck for fixed order calculations where calculations with finite quark masses have stalled at H + 3-jets at leading-order. Finite quark mass effects in H + 2j were studied in the high-energy limit in [28]. In this section, we recap some results and describe the necessary adaptation to incorporate the results in the current-structure of HEJ. We will start our discussion with the simplest leading-order configurations and gradually move on to more complex arrangements. We then describe how we supplement these expressions with the leading-logarithmic real and virtual corrections to all orders in α s which allows us to resum the logarithms in s/t.

Finite mass dependence in qQ → qHQ
We first consider the gluon-fusion production of a Higgs boson with two jets originating from two quarks with different flavours q and Q. At leading order, only two diagrams 2 Named after Fadin, Kuraev and Lipatov.
with all momenta left-to-right. The quark-loop insertions in the two diagrams are symmetric under charge conjugation. They can be written as a colour-diagonal vertex of the form where m is the quark mass, α s = g 2 s 4π is the strong coupling constant and v ≈ 246 GeV is the Higgs vacuum expectation value. For convenience, we list the form factors T 1 and T 2 in appendix A. The expressions are given there for a single propagating quark flavour; in practice, we sum the contributions from a top quark and a bottom quark propagating in the loop at the amplitude level. We therefore automatically have contributions from both flavours and the interference between them when we square the amplitude. The colour and helicity summed and averaged square of the matrix element then takes the factorised form where the invariants t 1 = q 2 1 , t 2 = q 2 2 are defined in terms of the t-channel momenta q 1 = p a − p 1 , q 2 = p 2 − p b , and S m qQ→qHQ is the current contraction We use the usual conventions for the QCD colour factors N C = 3, C F = 4 3 , C A = 3. Throughout this section we use double-bar notation to indicate summing over spins.
The amplitude has the same structure as in the limit of an infinite quark mass [29,41,42]. Indeed, the only difference is in the expression for V µν H , which in the limit is given by The factorisation of eq. (2. 3) into a current contraction and a product of t-channel propagators is exactly what enables us to perform HEJ resummation, as we will discuss in more detail in section 2.5. We stress that up to this point we retain the exact expression for the amplitude without having to resort to approximations valid in the high-energy limit.

JHEP04(2019)127
Here we have referred to initial quarks; the treatment of initial antiquarks is completely analogous. The only qualitatively different amplitude is qq → qHq, which receives contributions from both t-channel gluon exchange as in figure 1 and two annihilation diagrams with s-channel gluon exchange. Both sets of diagrams are individually gauge independent, but the annihilation diagrams are subdominant for a large invariant mass between the two jets, so that the leading contribution to the amplitude is the same as the pure-quark amplitude eq. (2.3). One may also consider the process qq → gHg, but this is not an FKL configuration and hence will not contribute a leading power in s/t to the matrix-elementsquared.

Finite mass dependence in gq → gHq
It was shown in [30] that in pure jet production the matrix elements for gluon-initiated 2 → 2 processes can also be described exactly as a contraction of currents where the current for the equivalent quark process (j ± µ (p o , p i ) in eq. (2.4)) is multiplied by a scalar factor K g /C F . For a backward-moving incoming gluon for example, K g is given by Using this instead of a scalar gluon impact factor or a current multiplied by C A /C F (the limiting value also referred to in [34]) improves the HEJ description of jet processes with incoming gluons for finite rapidity differences. The t-channel factorisation of an amplitude implies not only that each factor is independent of the momenta of the rest of the process, but that it is also independent of the particle-content of the rest of the process. This description of incoming gluons in inclusive dijet production is therefore also valid in H + 2j production and we will use it in what follows. At leading order, the gq → gHq in the initial state is significantly more involved than the process with two incoming quarks. Of the 20 diagrams contributing to the leadingorder amplitude, 10 can be obtained from charge conjugation. The remaining diagrams are depicted in figure 2.
The amplitude with full quark-mass dependence is known for general kinematics [20,21]; it does not have a t-channel factorised form. In the following subsections, we discuss the different hierarchies which can exist between invariants in the process and the expressions we will use to describe this process in the corresponding regions of phase space.

Central Higgs-boson emission
Let us first discuss the case where the Higgs boson is emitted in between the two jets, with a large rapidity separation from each jet. More concretely, we consider the momentum assignment g(p a )q(p b ) → g(p 1 )H(p H )q(p 2 ) with the hierarchy where s ij = (p i + p j ) 2 are invariant masses of the outgoing particles, Assuming the gluon to be emitted backwards, this hierarchy implies the rapidity ordering y 1 y H y 2 . The forward emission of the gluon is of course completely analogous.
In [28], it was shown that the amplitude in this limit assumes a similar factorised form as in the pure quark case. This is also true within the HEJ formalism [34]. As an example, the colour summed and averaged square of the helicity-conserving amplitude for a positive-helicity gluon and a negative-helicity quark can therefore be written as was given in eq. (2.6). The current contraction is given by as in eq. (2.4) for qQ → qHQ.

Peripheral Higgs-boson emission
We now consider the case where we drop the strong-ordering requirement between the Higgs boson and one of the jets. For the case where the Higgs boson is not strongly separated from a quark, we again invoke t-channel factorisation to treat this as in the qQ → qHQ -7 -process. There the result with V µν H is exact wherever the Higgs boson is emitted and hence in this case we again use eq. (2.8).
We now consider Higgs boson rapidities which are not strongly ordered with respect to the rapidity of the gluon in gq → Hgq. We will still require a separation from the quark, i.e. y 1 , y H y 2 or We denote such configurations as gq → Hgq, reserving the notation gq → gHq for the central Higgs-boson emission discussed in section 2.2.1. The t-channel factorisation of the amplitude is only guaranteed where there is a large rapidity separation between outgoing particles, and hence in this reduced limit we should not expect to recover a form with two t-channel poles as in eq. (2.8). However, as there is still a large rapidity separation to the quark line, we expect to find a factorised form about the pole in t 2 as follows: where the remainder of the amplitude has been written as an effective current j µ H . This current, dependent on the reduced set of momenta (p1, p H , p a ), is derived in appendix B. The derivation follows closely the approach in [30] for calculating the effective gluon current. In contrast to that case, or indeed that of a central Higgs boson, amplitudes flipping the gluon helicity also contribute. Explicit expressions are given in appendix B.

Finite mass dependence in gg → gHg
As noted above, the t-channel factorisation which arises in the limit of large rapidity separations implies that the building block corresponding to each end of the chain is independent of the rest of the process. We can therefore describe the gg-initiated state by taking the expressions for gq → gHq and adding the necessary change for an incoming gluon derived from pure jets. We find for central Higgs boson emission from eq. (2.8): (2.14) Likewise, for a Higgs boson emitted backward of both gluons where only the subset of hierarchies applies (eq. (2.10)), we find from eq. (2.11): Note that the gluon closest to the Higgs boson will always have the more complicated treatment derived in the previous subsection (i.e. it is the momentum of the gluon closest in rapidity to the Higgs boson which will enter j µ H ).

The first set of next-to-leading logarithmic corrections
We have now derived the HEJ description of the scattering amplitudes for all 2 → H + 2 processes which will contribute at leading power for Higgs-plus-dijets. At the start of this section, we identified the necessary subprocesses as the FKL configurations of flavour and momenta. In [34], we extended the HEJ framework to also describe Born processes where we have relaxed the requirement of rapidity ordering on exactly one gluon, allowing it to be emitted outside of the rapidity range defined by an outgoing quark, e.g.: The addition of the colour-neutral Higgs boson does not affect the following argument, so for a moment we will consider only the coloured particles. When constructed as a rapidity-ordered ladder diagram, see figure 3, the ordering above contains one t-channel quark propagator and one t-channel gluon propagator as opposed to the two t-channel gluon propagators one would find if y 1 and y 2 were reversed. It is therefore a non-FKL configuration, and as it is missing one gluon propagator it is suppressed by one power of s at matrix-element-squared level compared to the FKL configuration. This is formally therefore a next-to-leading logarithmic contribution to the dijet cross section; however one can still construct the leading logarithmic contributions to each particular subprocess. We denote these "unordered" configurations. This particular class of processes was chosen as it had been observed after matching to leading-order that they contributed significantly in regions of phase space with large transverse momentum for example. Their inclusion therefore allows HEJ to reduce its dependence on fixed-order matching [34] (see section 3 for a full discussion of all matching included in the current study). From eq. (2.17), these subprocesses at Born level have just one strong rapidity-ordering between the coloured particles and one therefore constructs an effective current to describe the q(p a ) → g(p 1 )q(p 2 )g * end of the chain, denoted j uno µ (p 2 , p 1 , p a ) (illustrated in figure 4). This must now carry two colour indices as it consists of terms with different colour flow. The matrix element for the case of central Higgs boson emission (where the Higgs boson is emitted between the quarks) is then given by where K uno = −1/2 and we are using figure 4). It is clear therefore that the arguments of S uno qQ→gqHQ are not independent; we have chosen to display the implicit dependence on q 1 , q 2 in anticipation of processes with additional gluons, where the dependence is explicit. The expression for j uno µ (p 2 , p 1 , p a ) is given in appendix C. Figure 3. Rapidity-ordered ladder diagrams for qQ → gqQ. The ordering in eq. (2.17) (left) contains just one t-channel gluon propagator, whereas the FKL configuration (right) contains the maximum number, two. Figure 4. For the rapidity ordering in eq. (2.17) where there is only strong ordering between y 2 and y 3 , we should only expect to find factorisation about the t-channel pole between these particles so the structure of the amplitude should be as shown. We find this in eq. (2.18).
So far in this section, we have discussed how to construct the HEJ approximation to Born-level matrix elements. In the rest of this section, we outline how to supplement these with the dominant corrections at all orders in α s .

HEJ resummation
In the previous subsection, we have outlined the HEJ approximation to the Born-level matrix element for f 1 f 2 → f 1 Hf 2 , f 1 f 2 → Hf 1 f 2 , qf 2 → gqHf 2 and their symmetric equivalents. These skeletons will provide all leading-logarithmic contributions and an important class of next-to-leading contributions to H+ ≥ 2-jet production. In order to achieve an all-order resummation we must add both the dominant real and virtual corrections. Our method was described in great detail for the case of an infinite quark mass in [34]. The presence of a finite quark mass does not affect the resummation procedure as the mass dependence enters only in the Higgs boson vertex V µν H and the effective current j µ H , and so the method remains unchanged. In the rest of this subsection we therefore briefly summarise the real corrections, the virtual corrections and the organisation of the cancellation of the poles to give an all-order finite result.

Real emissions
The dominant real corrections in the MRK limit arise in FKL flavour and momentum configurations. Where there are more than two outgoing coloured particles, the unique FKL configuration for a given incoming state is the process which has additional gluons emitted with rapidities in between two outgoing particles of identical flavour to the incoming particles. The high-energy limit, eq. (2.1), implies that these must be well separated from all other outgoing partons. Therefore, taking central Higgs-boson emission as an example, for a total of n emitted partons with momenta p 1 , . . . , p n , and Higgs boson emitted between j and j + 1, we have hierarchies of the form In this limit, it has been shown that the emission of the ith gluon can be described by a Lipatov vertex [43], V L , given by [29] where q i = p a −p 1 −. . .−p i is the incoming t-channel momentum and q i+1 = q i −p i+1 is the outgoing t-channel momentum (the q i will also have p H subtracted if the emission is forward of the Higgs boson). This vertex is derived from the five possible tree-level diagrams for qQ → qgQ, and then employing t-channel factorisation. At the matrix-element-squared level after summing over polarisations each contributes a factor of which multiplies the spinor string function, S (see below). This description of real corrections is the same as for the case of pure tree-level multijet production. A priori, additional gluons could be emitted off the massive quark loop coupling to the Higgs boson (see figure 5). However, in comparison to the emissions described before, such corrections are suppressed for a large rapidity separation between the Higgs boson and the gluons and will be neglected in the following. The absence of t-channel enhancement is obvious in the limit of a large quark mass, where the quark loop is absorbed into an effective local interaction.

JHEP04(2019)127
Defining the notation ·g· to mean an arbitrary number of gluons (including zero), we may therefore compactly write the HEJ approximation to the Born-level matrix element for where we define K q = Kq = C F for any flavour of quark/antiquark. Analogous expressions exist for the case of the Higgs boson emitted outside of the coloured particles in rapidity.

Virtual corrections
The leading-logarithmic terms of the virtual corrections can be obtained via the Lipatov Ansatz [37]. This is a prescription where, given a t-channel factorised matrix element as in eq. (2.22) and a corresponding hierarchy of scales, eq. (2.19), each t-channel pole is replaced: where t i = q 2 i , y i+1 and y i are the (ordered) rapidities of the emissions connected by the propagator. Using dimensional regularisation (D = 4 − 2 ) contains divergences as ε → 0. In the following we describe how to organise the cancellation of these poles with the poles arising from soft, real corrections. This prescription correctly describes the leading and even next-to-leading logarithmic terms, as verified up to N 2 LO [44].

Organisation of cancellation of poles
It is clear that eq. (2.24) will generate poles in ε. The divergences in eq. (2.22) are not immediately obvious; these appear when the matrix-element-squared is integrated over phase space where the momentum of the additional gluons described by V L goes to zero. The two sources of poles cancel exactly. In order to organise this cancellation, we introduce subtraction terms which in real-emission phase space apply for 0 < |p ⊥ | < λ. We use a subtraction term of for p k⊥ < λ in each of the squares of the Lipatov vertices, such that (2.27) for p k⊥ < λ. The subtraction term is integrable in D = 4−2 and the contribution is added to the virtual corrections, such that the finite contribution from the virtual corrections is then described through the prescription where the regularised Regge trajectory ω 0 is It can then be shown analytically (see e.g. [29,42]) that the poles which arise when integrating over the region 0 < |p ⊥ | < λ are equal and opposite to the virtual poles in eq. (2.24). The real radiation is regulated since (2.30) For values of |p ⊥ | 200 MeV, the limit is sufficiently good that the integral of the sum of Lipatov vertices and subtraction term can be ignored for lower values of |p ⊥ |. In that case we define a lower cut-off κ 200 MeV of the phase space integrals. The final results are stable under variation of both κ and λ; the results presented in this study were obtained with λ = κ. The final regularised, resummed expressions for the squared matrix element for the production of a Higgs boson in between partons j and j + 1 with incoming flavours f 1 , f 2 is then 3

JHEP04(2019)127
The same formula holds for a backward (forward) Higgs-boson emission if f 1 (f 2 ) is a quark or antiquark with j = 1 (j = n). For a peripheral emission close to a gluon, there is an equivalent expression but the first two lines instead mirror eq. (2.11) and we have one fewer t-channel pole as there is one fewer hierarchy. Explicitly, the matrix-element-squared is given by . . , n and t i = q 2 i . The regulated matrix element above is valid in the phase space of an arbitrary number of extra real gluon emissions each with |p ⊥ | > κ, provided they are between the extremal partons in rapidity. Note that the extremal partons play a special role and are not allowed to become soft (we do not include the necessary virtual corrections to regulate the fundamental spinor strings). In practice we require the extremal partons to carry a significant fraction of the extremal jet momentum to ensure that they remain perturbative.

Matching to fixed order
This section describes the matching of results within HEJ both to the full leading-order finite quark-mass matrix elements and to the NLO cross sections obtained for infinite top mass. We generate our event weights using the procedure outlined in [35], where we begin from fixed-order samples and supplement these with resummation. The resummation step only applies to the particle and momentum configurations discussed in section 2 (namely FKL configurations or FKL with one unordered gluon). If a given fixed-order event is not one of these configurations, it enters our final event sample with its weight unaltered.
As in [35] we choose the central factorisation and renormalisation scale for all the predictions as µ r = µ f = max(m H , m 12 ), where m 12 is the invariant mass of the system of the two hardest jets. This scale is chosen because the study in [45] indicates a better convergence of the perturbative result for pp → 2j than traditional p t -based scale choices such as e.g. µ r = µ f = H T /2, in particular for large dijet rapidity-spans. Since the formalism of the current study is of particular interest within the VBF-cuts, we choose a central scale which obtains reasonable uncertainty-estimates for the distributions that the cuts are based on. We note that, with this choice, the p t -based observables such as p H⊥ show the same pathological scale variance for large values that the invariant mass-based observable develops for renormalisation and factorisation scales based on the transverse momenta.
We observed in [35] that a central scale choice of µ r = µ f = H T /2 leads to distributions in rapidity and dijet invariant mass with values close to the upper edge of the scale variation -14 -

JHEP04(2019)127
band obtained when µ r and µ f are varied independently by a factor of two around this central scale choice, keeping their ratio between 0.5 and 2. The scale variance obtained with a central scale of µ r = µ f = H T /2 is therefore pathologically decreased for distributions at large dijet invariant mass or large rapidity separations, which are relevant for the VBF studies. While the scale variations obtained at NLO with a central scale choice of µ r = µ f = max(m H , m 12 ) are larger, they are also more reasonable as an indication of the uncertainty due to higher order corrections within the VBF-cuts. Section 3.1 describes the fixed-order samples available which we use as our starting point and the point-by-point matching applied to the resummation events. Section 3.2 then describes the matching performed for all events at the level of the total cross section.

Matching of exclusive amplitudes
HEJ allows for the perturbative series in each n-jet phase space point to be matched to fixed order [35]. This, obviously, is possible only if amplitudes for the n-jet phase space point are readily available. In this study, the fixed-order calculations are performed using Sherpa [26] with the extension of OpenLoops [27] for the evaluation of the pp → H + 2jprocesses with full quark-mass dependence. However, this fixed-order setup includes just the effects from the top-quark, and not also those of the loops of bottom-quarks. The effect of both top and bottom quarks is included in the resummation, and will be discussed later.
Additionally, the H + 3j-processes are not readily evaluated with the full quark-mass dependence even at leading order, and even using the infinite top-mass limit, only the pp → H + 2j-process is available at NLO (and then obviously pp → H + 3j at tree-level).
The limitations in the fixed-order results mean that the matching [35] within HEJ has to use many more different components than usual. We will describe them in the following sections. No point-by-point matching is performed for events with six or more jets. For such multiplicities the fixed-order results are expensive to compute, while typically only contributing to less than a percent for all shown observables.
On top of the matching of exclusive events described in the following, the final predictions for HEJ will be scaled with the ratio of the inclusive cross section for pp → H + 2j calculated at infinite top-mass for NLO and HEJ.
The described procedure obtains top and bottom mass dependence through the allorder results, matching to the full top-mass results for pp → H + 2j, and to pp → H + 3j and pp → H + 4j in the limit of infinite top-mass.

Two-jet matching with finite quark mass
The exclusive two-jet events are matched to full leading order, with finite quark mass effects. However, as our fixed-order setup allows for just the top-quark diagrams, technically the matching is performed by multiplying the all-order results containing both top and bottom mass effects with the ratio of the square of the full Born-level matrix element evaluated with just the top quark and the corresponding approximation within HEJ (using just the propagating top-quark, with no bottom-quark effects). The final event weights are therefore where M LO is the leading-order matrix element, M HEJ the all-order HEJ matrix element, and M HEJ, LO its truncation to leading-order. The superscript indicates the quark masses that are taken into account.

Three-, four-and five-jet matching with infinite quark mass
With the method of [35], the resummation could be constructed starting from event files from the calculation of Born-level Higgs-boson production in association with three jets including full momentum and mass dependence reported in [25]. However, since these are not available, the three-, four-and five-jet events will be matched in the infinite top-mass limit, which can be readily evaluated using Sherpa [26] and Comix [46]. Technically then, the reweighting of the event is performed with the ratio of the Born-level evaluation of the HEJ-approximation in the infinite top-mass and the full Born-level expression in the same limit, while the resummation is performed using the full expressions developed in section 2.5 (top and bottom included). The contribution from the matrix elements to the event weights is then where the "eff" superscript refers to the limit of an infinite top-quark mass. In this approximation, the interaction between the Higgs boson and gluons is described by an operator of dimension five, so that matrix elements exhibit unphysical scaling in the limit of large momenta. Since we choose not to include finite top-mass corrections in the truncated HEJ matrix element this effect cancels out in the ratio in eq. (3.2). The emission of quarks and gluons should resolve the dependence on the loop momenta only for large energies of the emission (compared to m t ). Since the bulk of each jet multiplicity consists of jet transverse momenta close to the defined jet threshold, the quark-mass effects should have only a small effect on the inclusive cross section. The quality of the approximation can be checked by applying a similar strategy of reweighting in pp → H + 2j, where the full result is known and can be checked against. The result of starting with the HEJ-approximation of the matrix element truncated at Born level with full top-mass dependence, multiplied by the ratio of the full Born-level result to the HEJ approximation, both evaluated in the infinite top-mass limit, is shown in figure 6, and compared to the Born-level result with both the full top-mass dependence and in the effective theory of infinite top-mass. The result obtained from m t → ∞ undershoots the full result by 5% for transverse momenta of the Higgs boson up to p H⊥ ≈ m t , at which point it diverges from, and overestimates, the correct cross section. The net result is that for the tree-level results, the infinite top-mass limit gives a good approximation to the integrated cross section obtained with the full top-mass dependence, as observed in [25,28], even if it is clear that the agreement in cross sections is accidental and will depend on the transverse cuts used.   While the corrections are relatively small and uniform for the differential cross section with respect to the azimuthal angle between the two jets and the rapidity of the Higgs boson, there are systematically increasing corrections to the distribution with respect to the invariant mass between the two jets, growing to more than 10% for m j 1 j 2 > 700 GeV.
We now turn the attention to studying the level of approximation to the full massdependent tree-level result by using the full mass-dependence in the HEJ-approximation to the full tree-level result followed by matching of the matrix elements in the infinite topmass limit (as we have to do for 3, 4 and 5 jets). The results are checked in the case of just two jets. If the reweighting factor was also evaluated with finite quark masses, the black and orange lines would be identical, and therefore the difference gives a measure of the quality of the approximation. We see that for rapidity-distributions, exemplified by that of -17 -

JHEP04(2019)127
the Higgs-boson, the level of accuracy obtained is roughly 5%. The accuracy is better than 12% in the distribution of the azimuthal angle between the two jets, and similar for the invariant mass between the jets, although here, and for the transverse mass-distribution of the Higgs boson, the corrections increase with increasing scale.
We emphasise that the seemingly good agreement in the distribution of y H between the results using the infinite top-mass and full top-mass dependence is completely accidental, and that the results presented with the dashed (orange) line and obtained using eq. (3.2) are more accurate.
We conclude that by using finite quark-masses in the simplified HEJ amplitudes, and applying matching in the infinite top-mass limits we can expect the result with finite top (and bottom) quark mass to be well approximated for distributions even of 3, 4 and 5 jets.

Matching of leading-order results to NLO in the infinite quark mass limit
The results of the resummation and matching procedure described so far will be compared to the best possible fixed-order result we can obtain. This consists of Born-level for full top-mass (but not including the small effect of the bottom-mass), reweighted bin-by-bin by the differential NLO K-factor calculated for infinite top-mass. The LO and NLO results for the distribution of (left) the rapidity separation of the hardest two jets and (right) the maximum rapidity-difference between any two hard jets, ∆y f b in pp → H + 2j with infinite top-mass are shown in figure 7. The NLO K-factor is particularly interesting: it has a linear growth in both cases and is large at large ∆y. Where it is plotted as function of ∆y f b , it goes to 1 for ∆y f b = 0. This represents a region of phase space dominated by exclusive 2-jet events. For the rapidity separation of the hardest two jets the K-factor reaches a factor of 3 for rapidity differences of ∆y 12 = 8, and for the most forward/backward jets reaches a factor of 6 at ∆y f b = 8. This obviously brings into question the validity of NLO-calculations at such rapidity-differences. The source of the apparent perturbative instability in the fixed-order result is treated systematically within HEJ. It is here worth mentioning that in the MRK-limit the all-order HEJ-results for pp → H +2j and pp → H + 3j will contain the same effect from the virtual corrections (and soft emissions) of a suppressing factor ∝ exp(ω(k 2 ⊥ )∆y f b ) ∝ (log(ŝ/p 2 t )) ω(k 2 ⊥ ) with ω(k 2 ⊥ ) < 0. However, when the perturbative series is terminated at NLO-accuracy, the effect of the expansion of the exponential suppression is included only in the events with Born-level kinematics. The suppression is missing at NLO in the corrections from real emissions because of the fixed-order termination of the perturbative series. At large rapidity-spans ∆y f b , this will inflate the NLO-result compared to the all-order result of HEJ, irrespectively of the choice of renormalisation and factorisation scale.
The balance between a suppression for H + 2j at NLO at large ∆y of the two-parton contribution and enhancement of the three-parton contribution discussed above is obviously influenced by the value of α s and therefore the scale choices. Indeed, the effect of choosing instead a central scale choice of µ r = µ f = H T /2 is illustrated in figure 8. As seen in the right plot, the K-factor tends to unity for ∆y f b → 0, rises to 1.5 at ∆y f b = 4, stabilises and then starts decreasing at ∆y f b ∼ 7. As a function of ∆y 12 , the K-factor starts at 1.3 for ∆y 12 = 0, and then decreases to 0.7 at ∆y 12  observed for the central scale choice of µ r = µ f = H T /2 may seem more appealing than the behaviour observed in figure 7; however, the variation obtained around this central scale will certainly underestimate the uncertainty from uncalculated higher orders, since the central scale choice leads to a result close to the edge of the results obtained by the variation. Furthermore, the scale variation band for NLO in figure 8 (left) increases with ∆y 12 , reaching −70% in the last bin, above ∆y 12 > 7. This indicates an instability of the NLO calculations for µ r = µ f = H T /2 at large rapidity differences. All the results presented in the following with µ r = µ f = max(m H , m 12 ) are also presented in appendix D for µ r = µ f = H T /2. Just as for fixed-order predictions, other processes like W +jets could be used in order to verify which of the scale choices obtains the best description of data.

Results of finite quark masses and all-order resummation
This section will first present results for a separate investigation of the higher-order effects included with HEJ compared to the fixed-order approaches. As mentioned before, we employ Sherpa in combination with OpenLoops to obtain the fixed-order predictions.
To evaluate the finite quark mass corrections within HEJ, we make use of QCDLoop [47]. We adopt our input parameters and cuts from [34,35], following the experimental analysis of [48]. Explicitly, we consider the gluon-fusion-induced production of a Higgs boson decaying into two photons, 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. For the photons, we require |y γ | < 2.37, 105 GeV < m γ 1 γ 2 < 160 GeV, the decay into two photons to 0.236%. We use the CT14nlo PDF set [49] as provided by LHAPDF6 [50]. 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 [48]: A discussion of the values chosen for the quark masses is in order. In the gluon-fusion production of a Higgs boson together with light-flavour jets the heavy quarks only appear in internal loops and are off shell. We therefore do not use on-shell masses, but instead prefer the MS mass-scheme. The scale µ m associated with the MS mass is a priori independent of the renormalisation scale used for the running coupling. It should be set to a scale characteristic for the heavy quark loop.
For the bottom quark, the mass is negligible compared to all other scales in the loop. Since the observables considered in this work depend only mildly on the bottom-quark mass, the exact scale choice has little impact on the prediction. To be definite, we use µ m b = m H and m b (m H ) = 2.8 GeV, which can be obtained from input values of m b (m b ) = 4.18 GeV [51], α s (m Z ) = 0.118 via renormalisation group evolution at two loops. The effect of higher orders in the evolution is negligible.
The effect of the top-quark mass is much more important. While there are ongoing efforts [52][53][54][55][56][57] to relate the very precise values reported by the LHC and Tevatron experiments [58][59][60] to a well-defined short-distance scheme, the top-quark MS mass is not known very precisely at the moment. For this project, the values chosen are µ mt = m H with m t (m H ) = 163 GeV, in line with direct determinations of the MS mass [61,62] and compatible with a pole mass of 173 GeV [63] within the uncertainties quoted in [62].

JHEP04(2019)127
Since the fixed-order setup can take into account the effects only of the top-quark, all results in section 4.1 are for finite top-mass only (no effects from the bottom-quark included). Section 4.2 investigates the effects on the results of HEJ from the finite top and bottom mass compared to the results obtained for infinite top-mass. Finally, section 4.3 compares the most precise predictions from HEJ, including both top and bottom mass effects, and the matching to fixed order discussed in section 3, to that of the fixed-order finite top-mass results, scaled to NLO accuracy, as described in section 3.2. Figure 9 compares the results obtained with finite top-mass at LO, the LO rescaled to NLO accuracy in the limit of infinite top-mass, and in all-order HEJ (using just finite top-mass but no contribution from the bottom quark). Comparing the results of p H⊥ in figure 9a for LO and the rescaling using the bin-by-bin K-factor calculated in the limit of infinite top-mass, one sees that the NLO K-factor (the ratio between the lines in blue and in green, indicated by the blue band in the lower plots) varies locally between 0.8 and 2 within ranges of the distributions checked. The NLO K-factor is decreasing for increasing transverse momentum p H⊥ , crossing unity at p H⊥ = 340 GeV.

Effects of higher perturbative orders
The NLO K-factors for the distributions in the invariant mass between the two hardest jets m 12 (figure 9b) and the rapidity-difference between the two hardest jets (figure 9c) have the same systematic behaviour of increasing K-factor as observed for ∆y f b in figure 7 and discussed there. The NLO K-factor for m 12 increases from 1.5 to 2.2 at m 12 = 1 TeV, and for ∆y 12 the NLO K-factor increases in a straight line from 1.5 to 3 at ∆y 12 = 8. This obviously then induces a large K-factor when a large rapidity-separation and invariant mass is required in the VBF-cuts, as illustrated for the φ 12 within these plots seen in figure 9d. It can also be seen in figure 9c that the ratio between HEJ and LO decreases linearly as a function of ∆y 12 ; this is an illustration of the logarithmic suppression of events with exactly two jets where ∆y = log(s/t) for large s.
The fixed-order matching bin-by-bin (as opposed to phase-space point by phase-space point employed with HEJ) does not ensure the same value for the integrated cross section. The effect of the matching will depend on the binning width etc. The size of the variation in the cross sections from the various distributions is one measure of the residual room for improvement in the matching. The integrated cross sections obtained from various distributions using the method of differential K-factors are listed in table 1. There is found to be very little variation in the integrated cross section of just 0.1 fb, well within the scale variation on the NLO-rescaled cross section of 6.2 +1.1 −1.2 fb, and an overall K-factor of 1.6. Within the VBF-cuts, the overall NLO K-factor is 2.2, and the NLO-rescaled cross section is found to be 0.53 +0. 15 −0.13 fb. Table 1 also contains the result for HEJ, rescaling the all-order results with finite top-quark mass with the ratio between the results obtained for infinite top-quark mass at NLO and at all orders with HEJ. Thus the matching is different to that applied at lowest order and normalises to the NLO cross section in the infinite top-quark mass limit. The inclusive cross section for HEJ matched as described is found to be 5.    Table 1. Cross sections obtained at LO, LO scaled with bin-by-bin K-factor for various distributions, and the HEJ with the inclusive cross sections scaled to NLO.

JHEP04(2019)127
the infinite top-quark mass. While the results for the inclusive cross sections are similar at NLO and in HEJ, the distributions differ significantly. As is evident from figure 9, the differential distribution from HEJ is harder in p H⊥ compared to the scaled LO-result, while the spectrum is decreasing significantly faster for both m 12 and ∆y 12 . This means that even though the total cross section for HEJ is matched to NLO (in the infinite topmass limit) with a scale-dependent K-factor of 1.4 +0.4 −0.4 , within the VBF-cuts the result of 0.23 +0.02 −0.03 fb happens to be closer (but with a reduced scale dependence) to the LO cross section of 0.24 +0.12 −0.08 . It is just a numerical coincidence of the cuts applied that the cross sections agree. As seen already in the discussion of the NLO corrections, the perturbative corrections are large in the VBF region. There is no reason to believe the perturbative series has converged already at NLO.

Effects of the finite top mass
The impact of the full top-quark mass-dependence on the Born-level result for pp → H +2j was already investigated in figure 6. While the effect on the integrated cross section is very small, the effect on the differential distribution in p ⊥H is enormous. The infinite top-mass approximation undershoots the full-top mass result by 5% for p ⊥H up to 200 GeV and then increasingly overshoots for increasing transverse momentum, reaching 40% error already at p ⊥H = 340 GeV. Similarly, for the invariant mass between the two hardest jets, the distribution for the infinite top-mass result starts off undershooting the true result by 5%, crossing at m 12 = 150 GeV and increasing to 16% by m 12 = 1 TeV. The error due to the infinite top-quark mass approximation is very small and uniform in the rapidity of the Higgs boson.
We now turn our attention to the impact of both the finite top-quark and bottom-quark mass on the results of HEJ. First, we list in table 2 the result for the cross section with inclusive-and the VBF-cuts for infinite top-quark mass and finite top-quark mass for fixed order (LO scaled with NLO in the limit of infinite top-quark mass). For HEJ we also list the result using both finite masses for the top-quark and bottom-quark. The finite top-quark mass has a much larger impact on the results of HEJ than at fixed order, which might at first seem surprising, since the results of HEJ are matched to the fixed-order results. The larger impact of the top-mass effects are therefore not a result of the approximations in HEJ. Instead, as is evident in the distributions of figure 9, the higher-order corrections of HEJ emphasise the distribution at larger p H⊥ , where the corrections from the finite quarkmass are large. Therefore, the top-quark mass corrections of the HEJ-results amount to a 9% reduction within the inclusive and 11% within the VBF-cuts. We do not observe any effect of the non-zero bottom mass beyond 1% for any of the observables studied. The impact obviously increases, if the bottom mass is chosen larger [25]. Figure 10 compares the results obtained with HEJ using the three different descriptions of quark masses, namely infinite top-quark mass, finite m t but m b = 0, and finite both m t and m b . Evidently, the effect of the finite m b is negligibly small and uniform in all the distributions. As seen already in figures 6 and 9, the approximation of infinite top-mass fails for transverse momenta significantly larger than the top-mass (illustrated here by a plot of the distribution in the transverse momentum of the Higgs boson in figure 10a).  Table 2. Cross sections obtained in fixed order perturbation theory (either full NLO for the results using infinite top-quark mass or LO scaled bin-by-bin with the K-factor obtained in the infinite top-quark mass limit) and in HEJ for pp → H + 2j with inclusive and VBF-cuts. See text for further comments.

JHEP04(2019)127
Similarly, the results using an infinite top-mass overshoot the result of finite top-mass by 20% at an invariant mass between the two hardest jets of 400 GeV, increasing to 40% at 1 TeV. This is relevant for the description of the contribution from the QCD process within the VBF-studies of pp → H + 2j. The corrections from finite quark masses to the distributions in ∆y 12 (figure 10c) or ∆φ 12 with additional VBF cuts (figure 10d) reach just 10%.

Final results for HEJ
In this section we compare the most accurate results obtained using the methods described in this study in HEJ to those obtained at fixed order. We start by comparing the observables already investigated previously; as such, the red and grey bands in figure 11 are identical to those on figure 10, but are here compared to the results of using Born-level with finite top-quark mass, rescaled bin-by-bin with the NLO K-factor obtained using infinite topquark mass. We see in figure 11a that the fixed-order result is significantly softer in the transverse momentum of the Higgs boson than the result obtained with HEJ. We have already discussed how this leads to a larger impact of the finite quark masses within HEJ than at fixed order. Figure 11b illustrates that the distribution in the invariant mass between the two hardest jets is increasingly suppressed for increasing m 12 in HEJ compared to fixed order. While the results are similar for small m 12 , the ratio of fixed order over HEJ reaches 1.5 at m 12 ≈ 500 GeV. Similarly, as illustrated in figure 11c the results of HEJ are much suppressed compared to NLO at large ∆y 12 . The ratio of fixed order to HEJ found here increases in a straight line finally reaching 5 at ∆y 12 = 8. As discussed around figure 7, this is due to the absence of a logarithmic suppression of the 3-jet component in the NLO prediction. Although the HEJ cross section is matched to the NLO value, this does not change the differences in the shapes of the distributions and the large K-factor at large m 12 therefore persists. Figure 11d shows the distribution with respect to the azimuthal angle between the two hardest jets, measured relative to the positive rapidity direction, thus exploring the full interval from −π to π. VBF cuts have again been applied in addition to the general cuts. These require a significant invariant mass and rapidity separation of the hardest two jets and hence the suppression in figures 11b and 11c translates into a large difference (around a factor of 2) in the cross section between the HEJ and fixed order predictions (as also seen earlier in table 2). The distinctive shape which arises as a result of the CP structure of the ggH vertex [12,13,64] is seen in all the predictions.
We present the results of figures 9-11 for the alternative central scale choice of H T /2 in appendix D. The main conclusions of the plots are unchanged: the impact of the higherorder corrections in HEJ lead to a harder distribution in p H⊥ , which enhances the finite quark mass and loop propagator effects. This in turn leads to a suppression of the prediction at large m 12 , and the predicted impact of a VBF cut is more severe in the all-order calculations of HEJ than that seen in fixed-order predictions.
In the final figure in this section, figure 12, we discuss an alternative to a traditional jet veto [12,13,[64][65][66]. We begin by defining two tagging jets t 1 and t 2 : firstly as the  (d) Figure 11. The results obtained with HEJ compared with fixed order for various key distributions. See text for further details.
hardest two jets in the event t 1,2 = j 1,2 and secondly as the most forward/backward jets, t 1,2 = j f,b . We may then construct y 0 = (y t 1 + y t 2 )/2 for each event. The event will then be vetoed if it contains a further jet with transverse momentum above 30 GeV in-between the two tagging jets which satisfies |y j − y 0 | < y c . This procedure applies to a larger region in rapidity than a traditional jet veto which is only applied to jets in-between the tagging jets. This means that the same level of suppression can be obtained with a higher (and therefore perturbatively safer) transverse momentum cut. In figure 12a, we show the results when we choose the two hardest jets as the tagging jets while in figure 12b the tagging jets are the most forward/backward jets. In both cases, the cross section has reached a plateau by about y c = 2. The difference between the two choices is relatively small but the cross section for a given value of y c is lower for the forward/backward choice for the tagging jets than for the hardest jet choice. As discussed in [67] this type of jet veto has for y c up to 1.5 very little impact on the VBF process itself (since most radiation is produced close in rapidity to the Born-level jets), and is therefore an efficient tool in distinguishing the contribution from the two processes for pp → H + 2j. We saw that the VBF cuts themselves have a relatively larger impact on the cross sections of HEJ than fixed order, because of the steeper fall-off with m 12 and ∆y 12 . Figure 12 shows that a further cut on jet activity will have a yet larger effect on HEJ compared to fixed order. This is all expected since the fixed order results fail to reproduce the rise in jet activity with increasing rapidity separation, which is observed in both data and HEJ [68,69].

Conclusions
We have calculated the gluon fusion contribution to H + 2j including both • leading logarithmic corrections inŝ/p 2 t to all orders in α s , and • full dependence on top and bottom quark masses, including the loop-propagator kinematic effects absent in the m t → ∞ limit.
The components necessary for implementing the full quark-mass dependence within the all-order resummation scheme of High Energy Jets (HEJ) were calculated, such that both the quark mass and the systematic logarithmic corrections within the VBF cuts could be investigated. This goes far beyond the current state-of-the-art fixed order predictions. The results thus obtained have been compared to the fixed-order full top-mass-dependent results evaluated at Born-level, but rescaled bin-by-bin with the NLO K-factor obtained in the limit of m t → ∞. While the fixed-order results obtained with finite m t differ very little from those obtained for m t → ∞ we find a much larger reduction of 9% on the inclusive cross section in HEJ; this is because the transverse momentum of the Higgs boson is found to be harder with HEJ than at LO, and that the finite mass-corrections are larger at large transverse scales. For the first time, our calculation allows the computation of the interference between the top and bottom quark contributions beyond leading order in H + 2j. We find that the interference is extremely small for the running values of m b (m H ) and m t (m H ), less than a percent for all observables.

JHEP04(2019)127
We find that for a scale choice of µ r = µ f = max(m H , m 12 ) the NLO K-factor increases systematically for both m 12 , ∆y 12 and most dramatically for ∆y f b (see figure 7). With a scale choice of µ r = µ f = H T /2 however, the K-factor decreases with increasing m 12 or ∆y 12 . The balance between the large virtual negative corrections and the real positive corrections are clearly scale dependent. The large corrections illustrate a serious perturbative instability of the fixed order expansion within the VBF-cuts. This instability is specifically addressed by HEJ.
At large ∆y 12 and m 12 , the all-order predictions from HEJ are systematically suppressed compared to fixed order. The discussion of scale choice is independent of the discussion of the behaviour at large m jj , and so is the conclusion that a resummation of the leading terms at large m jj leads to a reduction of the cross section within the VBF cuts. Our results show that the gluon-fusion contamination in VBF studies is less severe than the fixed-order estimate would imply. The finite-mass corrections to HEJ within the VBF-cuts lead to a further 11% suppression compared to the results obtained with infinite top-mass.

C The current for a single unordered gluon emission
In section 2.4, we use an effective current, j uno cd µ (p 2 , p 1 , p a ), to describe the emission of an unordered gluon (one additional gluon outside in rapidity of an FKL configuration). The current for q(p a ) → g(p 1 )q(p 2 )g * (q 2 ) was derived in [34] to be: This differs from our other currents as there is no longer a single overall colour factor, and hence colour factors (with free indices c and d) must be included. Upon contracting with another current squaring, this leads to terms with different colour factors. For example, for q(p a )Q(p b ) → g(p 1 )q(p 2 )Q(p 3 ), we find [34] M HEJ tree qQ→gqQ (C. 2) The factor we require in eq. (2.18) is therefore given by where we use the shorthand J µ = V µν H (q 1 , q 2 )j ν (p 3 , p b ).
D Results with µ r = µ f = H T /2 In this appendix we study the effect of using a central scale of µ r = µ f = H T /2 instead of the choice µ r = µ f = max(m H , m 12 ) used in the main text. In  The results obtained at NLO for the two central scale choices µ r = µ f = max(m H , m 12 ) and µ r = µ f = H T /2 are compared in figure 16. It is noteworthy that the difference in the results in figure 16d for the cross section within the VBF-cuts is similar to the difference between the results of NLO and HEJ obtained with the same scale.
Finally, figure 17 compares the results obtained for HEJ with the two central scale choices. The differences in the results for the distributions are larger than indicated by the scale variation. This is not surprising, since the leading logarithmic behaviour at large m jj is unrelated to β 0 -terms from the running of the coupling. As stated earlier, comparisons with data for other processes can determine which of these scale choices obtains the best description. The discussion of scale choice is independent of the discussion of the behaviour at large m jj , and so is the conclusion that a resummation of the leading terms at large m jj leads to a reduction of the cross section within the VBF cuts.   Figure 13. Predictions for various distributions obtained with HEJ, pure leading order, and leading order rescaled with differential K factors for the central scale choice µ r = µ f = H T /2. See figure 9 for the corresponding plots with µ r = µ f = max(m H , m 12 ). 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.