$hhjj$ production at the LHC

The search for di-Higgs production at the LHC in order to set limits on Higgs trilinear coupling and constraints on new physics is one of the main motivations for the LHC high luminosity phase. Recent experimental analyses suggest that such analyses will only be successful if information from a range of channels is included. We therefore investigate di-Higgs production in association with two hadronic jets and give a detailed discussion of both the gluon- and weak boson fusion contributions, with a particular emphasis on the phenomenology with modified Higgs trilinear and quartic gauge couplings. We perform a detailed investigation of the full hadronic final state and find that $hhjj$ production should add sensitivity to a di-Higgs search combination at the HL-LHC with 3 ab$^{-1}$. Since the WBF and GF contributions are sensitive to different sources of physics beyond the Standard Model, we devise search strategies to disentangle and isolate these production modes. While gluon fusion remains non-negligible in WBF-type selections, sizeable new physics contributions to the latter can still be constrained. As an example of the latter point we investigate the sensitivity that can be obtained for a measurement of the quartic Higgs-gauge boson couplings.


I. INTRODUCTION
After the Higgs boson discovery in 2012 [1] and subsequent analyses of its properties [2], evidence for physics beyond the Standard Model (BSM) remains elusive. Although consistency with SM Higgs properties is expected in many BSM scenarios, current measurements do not fully constrain the Higgs sector. One coupling which is currently unconstrained and has recently been subject of much interest is the Higgs self-interaction ∼ η, which is responsible for the spontaneous breaking of electroweak gauge symmetry in the SM via the potential with µ 2 < 0, where H = (0, v + h) T / √ 2 in unitary gauge. The Higgs self-coupling manifests itself primarily in a destructive interference in gluon fusion-induced di-Higgs production [3][4][5] through feeding into the trilinear Higgs interaction with strength λ SM = m h η/2 = gm 2 h /(4m W ) in the SM. The latter relation can be altered in BSM scenarios, e.g. the SM coupling pattern can be distorted by the presence of a dimension six operator ∼ (H † H) 3 , and di-Higgs production is the only * Electronic address: mdolan@slac.stanford.edu † Electronic address: christoph.englert@glasgow.ac.uk ‡ Electronic address: nicolas.greiner@desy.de § Electronic address: k.nordstrom.1@research.gla.ac.uk ¶ Electronic address: michael.spannowsky@durham.ac.uk channel with direct sensitivity to this interaction [6]. A modification solely of the Higgs trilinear coupling, which is typically invoked in di-Higgs feasibility studies, is predicted in models of µ 2 -less electroweak symmetry breaking, e.g. [7].
After the Higgs discovery, analyses of the di-Higgs final state at the high-luminosity LHC and beyond have experienced a renaissance, and di-Higgs final states such as the bbγγ [6,[8][9][10], bbτ + τ − [11][12][13], bbW + W − [11,12,14] and bbbb [11,12,15] channels have been studied phenomenologically, often relying on boosted jet substructure techniques [16] (see also an investigation [17] of rare decay channels relevant for a 100 TeV collider). Recent analyses by ATLAS and CMS [18,19] have highlighted the complexity of these analyses and the necessity to explore different production mechanisms to formulate constraints on the Higgs self-interactions in the future. This program has already been initiated by feasibility analyses of the hhj, hhjj and tthh production modes in Refs. [10,12,[20][21][22].
Di-Higgs production in association with two jets is a particularly important channel in this regard since this final state receives contributions from the weak boson fusion (WBF) production mode. The phenomenological appeal of the WBF mode is twofold. Firstly, the weak boson fusion component of pp → hhjj is sensitive to modifications of the gauge-Higgs sector [20,23], which can lead to large cross-section enhancements. Secondly, the QCD uncertainties for the WBF topologies are known and under theoretical control [24,25], such that a search for BSM electroweak-induced deviations is not hampered

General Information
The maximum effective rank in this group is 6.

Group 1 (4-Point)
General Information The maximum effective rank in this group is 4. by QCD systematics. This situation is very different from QCD-induced production [26], and can be attributed to the particular phenomenology of WBF-like processes [27,28] However, an additional source of uncertainty that was neglected until recently [20] is the correct inclusion of the gluon fusion contribution to pp → hhjj analyses. In contrast to single Higgs phenomenology, the correct inclusion of massive fermion thresholds is crucial to a reliable prediction of QCD-induced pp → hhjj [12].
Given that the cross sections in WBF hhjj production are very suppressed compared to WBF hjj production (the WBF hhjj cross section is ∼ 750 times smaller), we have to rely on the dominant hadronic Higgs decay modes to be able to observe this final state. This rules out one of the most crucial single Higgs WBF selection tools -the central jet veto [29]. The observation of WBF-induced pp → hhjj production is further hampered by the top threshold in the QCD-mediated process. Since the top threshold sets the scale of the di-Higgs subsystem, an analysis that tries to retain as many low p T Higgs bosons as possible leads to a QCD contribution that dominates over the WBF component when minimal WBF-like cut requirements are imposed [20].
In this paper we extend the discussion of Ref. [20] in a number of directions. We first perform a detailed comparison of EFT-approaches to QCD-mediated pp → hhjj against a calculation keeping the full mass dependencies of top and bottom quarks in Sec. II. We compare the QCD-induced pp → hhjj phenomenology to the WBF signature in Sec. III before we discuss general approaches to isolate the signal from the dominant top backgrounds in a hadron level analysis in Sec. IV. This sets the stage for a discussion about the prospects to isolate the WBF and GF components in Secs. IV A and IV B, followed by a study on constraining V V hh couplings using the WBF induced signal in Section IV C. We focus on collisions with 14 TeV throughout.

A. Finite top mass effects
It is well known that effective field theory approximations in the m t → ∞ limit cannot be invoked to study di-Higgs final states at colliders in a reliable way due to the effects of top-quark threshold [4,30]. Further, the breakdown of the m t → ∞ approximation is worsened in the presence of additional jet emission [12,31]. Finite  m t effects must therefore be considered for all QCD di-Higgs production channels, which will be required to set the best limits on the Higgs self-coupling or formulate a realistic estimate of the GF contribution in a WBF-like selection.
The computational challenges in QCD-mediated hhjj production are significant, with the gluon-fusion channels particularly time consuming even when using state-ofthe-art techniques. The standard method of simulating a differential cross section from unweighted events is not feasible in this case, and we instead use a reweighting technique that is exploited in higher order calculations and experimental analyses (see e.g. [32]).
We generate GF hhjj events by implementing the relevant higher dimensional operators in the m t → ∞ limit obtained by expanding the low-energy effective theory [33] in MadEvent v5.1 [34] using the FeynRules/Ufo [35] framework. * This allows us to sample a weighted set of events that we subsequently feed into our analysis solely depending on their final state kinematics. If an event passes the selection requirements of a certain search region, we correct for the full mass dependence using a reweighting library based on GoSam package [36] at this stage. The reweighting employs exactly the same matrix elements used for the event generation and the trilinear coupling is steered through a modification of the GoSam matrix element, i.e. variations of the trilinear coupling are part of the reweighting. A selection of Feynman diagrams which contribute to the gluon fusion signal are illustrated in Fig. 2. The GoSam code used for the reweighting is based on a Feynman diagrammatic approach using QGRAF [37] and FORM [38] for the diagram generation, and Spinney [39], Haggies [40] and FORM to write an optimised fortran output. The reduction of the one-loop amplitudes was done using Samurai [41], which uses a d-dimensional integrand level decomposition based on unitarity methods [42]. The remaining scalar integrals have been evaluated using OneLoop [43]. Alternative reduction techniques can be * The effective theory implementation can be modified in the sense that only one effective vertex insertion is allowed. This is gives only a mild ∼ 10% effect in the tail of the distribution, and is not relevant for an order one EFT/full theory rescaling, see below.
used employing Ninja [44] or Golem95 [45]. To validate the reweighting procedure we regenerated the code that has been used in [20] with the improvements that became available within GoSam 2.0, in particular improvements in code optimisation and in the reduction of the amplitudes. For the reduction we used Ninja, which employs an improved reduction algorithm based on an Laurent expansion of the integrand. This leads to substantial improvements in both speed and numerical stability compared to the previous version. We combined the code with a phase space integration provided by MadEvent [46]. Further substantial speed-up has been obtained by Monte Carlo sampling over the helicities rather then performing the helicity sum. This enabled us to perform a full phase space integration and we found full agreement within the statistical uncertainties between the result obtained from reweighting and the result from the full phase space integration.

B. Phenomenology of QCD-mediated hhjj production
Top thresholds are particularly prominent in the di-Higgs invariant mass distribution, which is thus well suited to benchmark the relation of the finite m t limit to the effective theory of Eq. (2). Other observables constructed from the six particle final state are also relevant when performing a targeted phenomenological analysis, and we discuss both these and the phase space dependent parton-level reweighting in detail in the following.
In Figures 1, 3, and 4 we show a selection of hhjj final state observables for inclusive cuts p T,j > 20 GeV and |η j | < 5, no cuts on Higgs bosons are imposed. We label Higgs bosons and jets according to their hardness, i.e. p T,h1 > p T,h2 and p T,j1 > p T,j2 . The cross sections are given in Tab. I. The inclusive gluon fusion cross section is about 2.5 times larger than the WBF cross section approximately independent of the value of the Higgs trilinear coupling.
As previously established in [4,12,20] the di-Higgs system is badly modelled by the effective theory which under-and overshoots the full theory cross section at low and high momenta respectively. For pp → hhjj this is a qualitatively similar behaviour compared the pp → hh(j) production: The m hh distribution is the crucial observable which parametrises the finite top quark mass effects. The EFT describes the low maximum transverse Higgs momenta p T,h1 reasonably well, as shown in Fig. 1(a). The jet emission on the other hand integrates over a considerable range of m hh , and the ratio of full theory vs effective theory smaller than one for the finite m t limit produces a smaller integrated cross section than the m t → ∞ limit for the jet kinematics.
Considering just the dijet system in Fig. 3, we observe that the jet kinematics is not severely impacted by the reweighting procedure upon marginalising over the di-Higgs kinematics. The phase space dependence of the    dijet invariant mass Fig. 3(a) is relatively mild aside from the total rescaling of the inclusive cross sections, and the ratio for the pseudo-rapidity distribution of the jets is nearly flat, Fig. 3(b). This is also true for the azimuthal angle difference ∆φ jj . The angular distributions of the leading members of the jet-Higgs system are relatively mildly impacted by the reweighting too Fig. 4(b). This agrees with the m hh being the observable most sensitive to the top threshold (as in pp → hh(j)), and is also sup-ported by the larger impact of the reweighting of m hh in Fig. 4(a). A reweighting based on m hh to correct for finite top mass effects suggests itself for future analyses as a time-saving approach with reasonable accuracy.

III. THE WEAK BOSON FUSION CONTRIBUTION
The weak boson fusion contribution to pp → hhjj has received considerable attention recently and precise higher-order QCD corrections have been provided in [24,25,28]. Due to the sensitivity of the WBF contribution to both the trilinear coupling and the quartic V V hh (V = W, Z, γ), as shown in the Feynman diagrams in Fig. 5, weak boson fusion to two Higgs bosons can, in principle, provide complementary information about BSM physics which remains uncaptured in pp → hh(j) and pp → tthh [23].
We generate WBF samples with varying λ using MadEvent v4 [47] and normalise the cross section to NLO accuracy [24]. The WBF hhjj contribution shares the QCD properties of WBF hjj production [27] which means it shares the distinctive ∆η(j1, j2) distribution shown in Fig. 6(a): To produce the heavy di-Higgs pair we probe the initial state partons at large momentum fractions. This together with a colour-neutral t-channel exchange of the electroweak bosons [48] (see also [49]) leads to energetic back-to-back jet configurations at large rapidity separation and moderate transverse momenta with a centrally produced Higgs pair. The production of an additional Higgs boson in comparison to single Higgs production via WBF leads to a cross section reduction by three orders of magnitude (see Tab. I) in the SM. Such a small inclusive production cross section highlights the necessity of considering dominant Higgs decay channels such as h → bb and h → τ + τ − and the non-availability of central jet vetos [29] as a means to control the background and GF contribution in a targeted analysis as a consequence.
The gluon fusion contribution is bigger by a factor of 2.5 than the WBF component of hhjj production, however, with increasing invariant di-Higgs mass the WBF  contribution is enhanced relative to GF production as a consequence of the suppression above the 2m t threshold, as shown in Fig. 6(b).
Since we cannot rely on vetoing hadronic activity in the central part of the detector, a potential discrimination of GF from WBF needs to be built on the following strategy, which we will investigate in Sec. IV: • To isolate the di-Higgs (WBF+GF) signal we can exploit the relative hardness of the di-Higgs pair which peaks around ∼ 2m t . Such hard events are less likely to be produced by (ir)reducible backgrounds.  • Focussing on large m hh we can enhance WBF over GF by stringent cuts on the jet rapidity separation. This will also imply a significant decrease of QCDdominated backgrounds.
• By explicitly allowing central jet activity, we can exploit the colour correlation differences in WBF vs GF to further purify our selection. Since colour flow is tantamount to energy flow in the detector, event shapes are particularly well-suited observables for unravelling the colour correlations in the final state once the reconstructed di-Higgs pair has been removed † . This strategy was first proposed for single Higgs analyses in [51] (see also [52]).

IV. TAMING THE BACKGROUND
For our hadron-level analysis we shower our signal samples with Herwig++ [53] and generate backgrounds as follows: ttjj, tth, Zhjj, and ZZjj with Sherpa [54], and ZW W jj with MadEvent v5. We find the dominant backgrounds to be ttjj and tth production, for which next-to-leading order results are available [55,56] and we use inclusive K factors K ttjj ≃ 1 and K tth ≃ 1.5 to estimate the higher order contributions to these backgrounds. Higgs branching ratios are set to the values agreed upon by the Higgs Cross Section Working Group [57]. † A detailed discussion of event shapes at hadron colliders can be found in [50]. We begin the hadron-level analysis implemented in Rivet [58] by recreating a base selections similar to [20]: ‡ 1.) We require two tau leptons using a two tau-trigger based on staggered transverse momentum selection cuts p T ≥ 29, 20 GeV in |η τ | < 2.5 and assume a flat tau tagging efficiency of 70% with no fakes. Jets are constructed by clustering R = 0.4 anti-k T jets using FastJet [59] with p T,j ≥ 25 GeV and |η j | ≤ 4.5.
2.) The two leading jets are b-tagged with an acceptance of 70% and fake rate of 1% [60] in the central part of the detector |η j | < 2.5. We remove events if either of the two leading jets overlaps with a tau. Any additional jets which do not overlap with a tau are considered as potential "tagging jets", of which we require at least two.
3.) As a final step of this base selection we require the b jet and tau pairs to reproduce the Higgs mass of 125 GeV within ±15 and ±25 GeV respectively. § The signal and background cross sections after these cuts are presented in the Base Selection column of Table II. We find that the background contribution of ttjj dominates with tth also providing a larger-than-signal background resulting in S/B ∼ 1/300, making a study based ‡ Our analysis has been validated with two independent implementations. § A high mass resolution is a crucial cornerstone of any successful di-Higgs analysis to assure a minimum pollution of Z boson decay backgrounds [13]. only on these selections extremely challenging. Since we only have ∼ 40 expected gluon fusion and ∼ 10 expected weak boson fusion events at 3 ab −1 luminosity, additional selections must also be careful to retain enough signal cross section to allow statistically meaningful statements to be made with a finite amount of data.
Shape comparisons for the rapidity and dihiggs invariant mass distributions as motivated in the previous section are shown in Fig. 7. Indeed, as expected, cutting on the angular distance of the jets will serve to both purify towards a WBF-only selection at a reduced background rate. The dominant backgrounds are unlikely to produce a large invariant mass m hh . However the WBF contribution, due to the lack of the 2m t threshold peaks at a considerably lower invariant mass, leading to significant decrease of the WBF contribution for a reasonably strong cut on m hh , which is required to observe the hhjj signal at the given low signal yield, even at 3 ab −1 luminosity.

A. Prospects to isolate gluon fusion
We can extend the analysis outlined in Sec. IV with the aim to purify the selection towards the GF component. ¶ We make use of the hard Higgs candidates to greatly reduce the backgrounds by requiring m hh ≥ 500 GeV and additionally require ∆η(j1, j2) ≤ 5 to minimise the weak boson fusion contribution. The signal and background cross sections after these cuts are applied are presented in the 'GF Selection' column of Table II. ¶ Following the analysis of [61], we can expect negligible interference between WBF and GF and which allows us to make this distinction.
The total background is reduced by a factor of ∼ 100 while the gluon fusion contribution only is reduced by a factor of ∼ 2.5 which allows for an encouraging S/ √ B ∼ 1.7 with 3 ab −1 of data. The weak boson fusion contribution is also suppressed compared to GF which allows for a clean probe of the physics accessible in the gluon fusion contribution.

B. Prospects to isolate weak boson fusion
Similarly we can extend the analysis towards isolating the WBF component. Since it has slightly softer Higgs candidates we require m hh ≥ 400 GeV and ∆η(j1, j2) ≥ 5 to reduce both the gluon fusion and background contributions. The signal and background cross sections after these cuts are applied are presented in the 'WBF Selection' column of Table II. The total background is reduced by a factor of ∼ 300 while three times more of the weak boson fusion contribution is retained compared to the GF selection, resulting in S/ √ B ∼ 0.8 with 3 ab −1 of data due to the large reduction in the gluon fusion contribution. However even so the WBF selection is composed of one-to-three parts GF to WBF, which means measurements of physics that only enters the weak boson fusion contribution will need to take this gluon fusion "pollution" into account.

C. Constraining the quartic V V hh contribution
As mentioned in Section III there is a contribution from quartic V V hh vertices to the WBF induced signal, and modifications of the corresponding g V V hh couplings away from their SM values using the Higgs Cross Section Working Group κ framework [57] will greatly enhance the signal cross section. This allows us to constrain ζ defined by g V V hh = ζ ×g SM V V hh . To achieve this we have generated events with varying ζ using MadEvent v5 and applied the WBF selections described in Section IV B to estimate the enhancement of the signal, which is compared to expected cross section limits on the signal with 3 ab −1 of data in the WBF selection under the assumptions of no systematic uncertainties and 20% total systematic uncertainties for comparison. The results are presented in Figure 8. We find that in the more realistic scenario of 20% systematic uncertainties the expected constraint on the g V V hh couplings is 0.55 < ζ < 1.65 at 95% confidence level. A measurement of pp → hhjj is therefore crucial to constrain new physics which enters predominantly through enhancements to g V V hh .

D. Event shapes of the tagging jets system
The analysis strategies outlined so far have mainly relied on exploiting correlations in the di-Higgs system, with only ∆η(j1, j2) carrying information about the tagging jets. Following similar applications in the context of single Higgs production [51], we investigate a range of event shapes in the tagging jets system in the following, which could offer additional discriminating power through capturing colour correlations in the different signal contributions beyond angular dependencies. More specifically, we will focus on N -jettiness [62,63] and thrust major which provided the best results.
We calculate N -jettiness by minimising where C is a normalisation which cancels when taking the ratio of two τ s, the sum is taken over all visible momenta which do not belong to one of the identified Higgs candidates within |η| < 5, and ∆R k,n is the distance in the η−φ plane between the k-th momentum and the n-th reference vector. τ 3/2 is then explicitly given by τ 3 /τ 2 . Thrust major is defined by where n T is the normalised thrust vector Again the sums run over all visible momenta which do not belong to one of the identified Higgs candidates within |η| < 5. We find τ 3/2 and T maj show promise for improving the WBF selection, but the signal cross section is already too low for us to be able to make meaningful use of this insight. The τ 3/2 and T maj distributions after the GF and WBF selections have been applied are presented in Fig. 9. Cutting, e.g., on T maj < 0.05, the gluon fusion contribution is reduced by 80%, while the WBF contribution is reduced by only 55% amounting to a total of 2 expected WBF and 0.3 expected GF events, with backgrounds very strongly suppressed. This means that WBF can in principle be observed at a small rate that can be used to set constraints on new physics in an almost GFfree selection with greatly reduced backgrounds.
The event shape distributions can also be used to greatly reduce the background in the GF selection, Fig. 9(c). It should be noted that these improvements of GF vs WBF vs background ultimately depend on underlying event and pile up conditions and have to be taken with a grain of salt at this stage early in run 2. However the clear separation that can be achieved with these observables indicate that an analysis employing MVA techniques could, at least in theory, significantly improve the results presented here. These techniques may also prove useful at a 100 TeV collider where the dihiggs production cross-section is substantially higher [10].

V. SUMMARY AND CONCLUSIONS
After discovering single Higgs production at the Large Hadron Collider, new analysis strategies need to be explored to further constrain the presence of new physics beyond the Standard Model. Higgs pair production is pivotal in this regard as constraints from multi-Higgs production contain complementary information, in particular with respect to the Higgs boson's self-interaction. Cross sections for di-Higgs production are generically small at the LHC, which highlights the necessity to explore other viable channels than pp → hh to enhance sensitivity in a combined fit at high luminosity. To this end, we have investigated pp → hhjj production in detail in this paper. Keeping the full top and bottom mass dependencies, we find sensitivity of pp → hhjj searches at the LHC for production in the SM and beyond. The gluon fusion contribution remains important at high invariant di-Higgs masses where the dominant backgrounds can be suppressed to facilitate a reasonable signal vs background discrimination. Unfortunately, the gluon fusion contribution remains large even for selections that enhance the weak boson fusion fraction of pp → hhjj events. This "pollution" is important when such selections are employed to set constraints on new physics effects that enter in the WBF contribution exclusively. Large new physics effects in the WBF contribution can still be constrained, which we have illustrated through an investigation of the constraints that can be set on deviations of the quartic V V hh couplings from their SM values with the HL-LHC, demonstrating that a measurement of pp → hhjj will provide a powerful probe of these. Employing observables which are intrinsically sensitive to the different colour correlation of WBF compared to GF, the discrimination between GF, WBF, and background can be further improved. However, the signal cross section is typically already too small to use such a strategy to constrain the presence of new physics if those effects are only a small deviation around the SM. If new physics effects are sizable, such an approach will remain a well-adapted strategy to minimise GF towards a pure WBF selection.