$VHH$ production at the High-Luminosity LHC

We study the phenomenology of associated production of a vector boson with a pair of Higgs bosons ($VHH$) at the High-Luminosity LHC. Despite the low rate of this channel, the scaling of the cross section suggests a measurement could be a useful probe of modifications of the trilinear Higgs boson coupling and anomalous interactions in the gauge-Higgs sector. We focus on both $WHH$ and $ZHH$ production, using the leptonic ($W \to l \nu$, $Z \to ll$, $Z \to \nu \nu$) decay modes of the vector bosons and the $HH \to 4b$ di-Higgs decay mode. We show that top pair backgrounds are problematic for the $W \to l \nu$ and $Z \to \nu \nu$ channels, leaving $Z \to ll$ as the most promising decay mode. However, even for this channel, we find limited sensitivity due to a low signal rate. We discuss some potential avenues for improvement.


I. INTRODUCTION
Since the discovery of a scalar resonance with the properties of a Standard Model-like Higgs boson by the ATLAS and CMS collaborations in 2012 [1,2], the focus of the experiments and phenomenological community has shifted towards precise measurements of the properties of this new particle and its interactions [3]. Of particular interest is the trilinear self-interaction coupling, λ, since this would provide a first model-independent measurement of the shape of the scalar potential, that could be related, e.g., to models of strong first-order phase transition necessary for baryogenesis [4].
At colliders the focus has been on direct di-Higgs boson production as the measurement channel of choice. We will continue this tradition here. We note, however, that single Higgs production channels are also sensitive to the trilinear coupling at one-loop [5][6][7] and these may turn out to be competitive with direct production, especially at future lepton colliders where single Higgs-strahlung can provide the most sensitive channel overall.
The leading gg → HH + X channel of di-Higgs production at the LHC has been studied in a range of final states in the phenomenological literature [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], and by ATLAS and CMS [26][27][28][29][30][31]. Additionally di-Higgs final states have been studied in associated production with a top quark pair [32][33][34] and in associated production with two jets (which includes the leading vector boson fusion contribution) [35][36][37][38]. In this paper we focus on associated production with a weak boson which has previously been studied at hadron colliders [39][40][41][42][43] and lepton colliders [44]. Our objective is to perform a realistic sensitivity analysis of this production channel to modifications of the trilinear Higgs coupling λ and anomalous quartic gauge-Higgs interactions of the form V V HH. Foreshadowing the results of our analysis, we choose to work in a simple "anomalous" coupling framework, where the Standard Model (SM) couplings are modified by simple rescaling: λ = λ SM (1 + c 3 ), where λ SM is the SM value of the triple Higgs boson coupling and c 3 parametrises the modifications coming from new physics. Equivalently, we separately consider modifications of the quartic gauge-Higgs couplings through g V V HH = g V V HH,SM (1 + c V V HH ) for V = {Z, W }.
The diagrams which contribute to this process at tree level are presented in Figure 1. As noted in [41], the interference between the λ contribution and the other diagrams in this channel is such that it could potentially offer a sensitive probe of λ > λ SM . This can be most easily seen if we amputate the quarks to focus on the V µ → V ν HH subdiagrams setting the gauge boson mass to equal the Higgs boson mass, m V = m H , for simplicity, where we find FIG. 1: The diagrams that contribute to pp → V HH at tree level. The third one is accompanied by a 'u channel' one with pH 1 ↔ pH 2 which we do not draw.
at threshold: where λ = (1+c 3 )λ SM . This suggests the interference pattern is such that the cross section is smallest for λ ≈ −2λ SM , in contrast to gg → HH and vector boson fusion qq → qqHH where this occurs for λ > λ SM , as demonstrated in  Another potentially interesting feature of V HH production is that it gives us sensitivity to the W W HH and ZZHH vertices, which could be modified in strongly coupled and extra dimensional models of electroweak symmetry breaking [45]. These are probed very efficiently by vector boson fusion qq → qqHH production [35,36,38] since the contributing diagrams are such that the high di-Higgs invariant mass (m HH ) region is only suppressed by valence quark parton density functions which leads to large sensitivity when using boosted reconstruction techniques. However this channel is only sensitive to a linear combination of the W W HH and ZZHH vertices, whereas the ability to tag the vector boson in V HH production could potentially allow us to constrain these independently as ZHH and W HH production are sensitive to g ZZHH and g W W HH separately at leading order. A similar observation motivated the study of W W H production in [46].
We will work at √ s = 13 TeV proton-proton collision energy throughout the paper. Since the signal processes are qq-initiated at tree-level the cross sections grow more slowly when going to higher collider energies as compared to gluon fusion-initiated di-Higgs production channels; for this reason the qualitative changes that can be expected by considering, e.g., a 100 TeV proton-proton collider will be limited to a potentially larger data samples and detector improvements.

II. SETUP OF CALCULATION AND ANALYSIS
We use the MadGraph5 aMC@NLO framework [47] to generate both signal and background events at leading or next-to-leading QCD order [48][49][50][51][52][53][54] depending on the number of legs of the process. We use MadSpin [55] to decay particles with the correct spin correlations (reweighting the branching ratio to the HXSWG recommendation for m H = 125.0 GeV in [56]) and Herwig 7 [57][58][59][60][61][62] to shower and hadronise the parton-level events. We employ the Rivet framework [63] to analyse the hadron level events.
We use the default b-tagging implementation in Rivet which "ghost-associates" [64,65] b and c mesons to the jets to define b and c jets, which are then assigned b-tags with an efficiency of 77% for b jets, falling to (100/6)% for c jets and (100/134)% for light jets, corresponding to a standard operating point for the ATLAS MV2c10 algorithm [66].

III. ZHH PRODUCTION
ZHH production does not have an a priori obviously superior candidate between the invisible (Z → νν) and and leptonic (Z → ll) decay channels of the Z due to the lower branching ratio to charged leptons. The cross section before branching ratios at NNLO QCD * is ∼0.37 fb [43,67]. Previous phenomenological studies of ZHH production at parton level have employed the invisible decay channel [41], † and it also forms the most sensitive ZH channel in the recent H → bb analysis by ATLAS [68]. These exploit the fairly soft missing energy spectrum of top pair backgrounds, leaving Z+heavy flavour backgrounds dominant while gaining statistics from the relatively large Z → νν branching ratio. For this study we have investigated both of these Z decay channels.

A. Z → νν
The Z → νν channel is attractive due to the relatively large branching ratio: This gives a ZHH (Z → νν, HH → 4b) NNLO QCD cross section of ∼0.025 fb. Additionally, aggressive cuts on the missing transverse energy |E miss T | and vetoes on identified leptons allows multijet and top pair backgrounds to be controlled enough for this channel to be the most sensitive in the ZH (H → bb) context [68]. We do not take into account the |E miss T | trigger efficiency in our analysis: according to current ATLAS performance [69], this would lower our cross sections by a factor of ∼ 2 due to the cut of 100 GeV. However, one can expect this to change in the future HL-LHC runs. Unfortunately we find that top quark pair (associated) production is a very challenging background for the Z → νν channel of ZHH (HH → 4b) production. Following the strategy of the ZH analysis, we consider further avenues to control it: • Vetoing on prompt leptons requires them to be hard enough to be identified and within the detector volume.
Defining veto-able leptons as those with p T > 5 GeV within the inner detector |η| < 2.5 (similar to the lepton veto used in monojet analyses, see, e.g. [70]) with perfect identification efficiency (even for τ leptons) improves the signal-to-background ratio considerably. However some 9% of semi-leptonically decaying top pairs in tt production (the fraction varies only slightly for ttbb and ttH production) still pass and in fact form the dominant background.
• The missing energy spectrum is softer for the top backgrounds than the signal. However due to the already low signal cross section we keep the missing transverse energy cut at a relatively low value, |E miss T | > 100 GeV. * The relatively large NNLO/NLO K-factor of 1.2 is caused by the introduction of gg → ZHH contributions at NNLO. Due to the low total cross section of the pp → ZHH process, we only approximately consider this channel through the K-factor and do not include it explicitly in our Monte Carlo simulations. † As far as we understand, no parton showering was employed in the study of [41].  • We require the kinematics of the four b-tagged jets i, j, i , j in the event to resemble those of a pair of Higgs decays by minimising • We further attempt to reduce top backgrounds by finding the three jet permutation t * which most closely reconstructs the three jet mass m t = 172.5 GeV with a sub-permutation W * which reconstructs the two jet mass m W = 80.4 GeV at the end of the analysis by minimising However we ultimately find that it is difficult to use this to improve the sensitivity since the scales in the signal (with 4 jets from two Higgs decays) make it easy to fake a top candidate, in particular when the cross section is dominated by diagrams containing the trilinear coupling, as shown in Figure 3b.

Require |E miss
3. Require at least 4 jets with p T > 40 GeV.
4. Require the 4 leading jets to be b-tagged using the definition and efficiencies detailed in Section II.

5.
Require that these b-tagged jets have the kinematics of a HH pair decay, χ HH < 1.6.

B. Z → ll
The Z → ll (l = e, µ) channel, when compared to Z → νν above, suffers from a smaller branching ratio giving an NNLO QCD parton level cross section for ZHH (Z → ll, HH → 4b) of ∼0.008 fb. However it allows for top backgrounds to be controlled through a combination of aggressive cuts on m ll which can be justified due to the excellent lepton momentum resolution of the LHC experiments, and requiring |E miss T | < 50 GeV. Since we find that Z + heavy flavour backgrounds are dominant and are able to reconstruct the Z boson completely we have investigated angular observables in the Z, H 1 , and H 2 systems (where H 1,2 are the leading and sub-leading reconstructed Higgs candidate in p T , respectively) and find that ∆η(Z, H 1 ) can be used to significantly reduce this background at little signal cost, see Figure 4a. Other angular observables may also carry some additional information, however due to the low signal rates a multivariate approach would be required to make use of this (while also balancing the non-negligible top pair backgrounds).
The Z → ll analysis steps are as follows: 1. Require exactly two same flavour opposite charge electrons or muons inside |η| < 2.5 with p T > 25 GeV.
2. Require these leptons to have an invariant mass compatible with that originating from a Z boson decay, |m ll − m Z | < 5 GeV.
3. Require |E miss 4. Require at least 4 jets with p T > 40 GeV. Veto event if any of these overlap with a lepton.
5. Require the 4 leading jets to be b-tagged using the definition and efficiencies detailed in Section II.
6. Require that these b-tagged jets have the kinematics of a Higgs boson pair pair decay, χ HH < 1.6.
7. Require that the leading Higgs candidate H 1 and the reconstructed Z boson are not too far separated in pseudorapidity, ∆η(Z, H 1 ) < 2.
The visible cross sections after these selections are applied for the signal and backgrounds are presented in Table I. Using our selections the two channels end up being competitive with each other: the Z → νν channel has higher signal statistics but a low signal-to-background ratio S/B ∼ 1/660, while the Z → ll channel has lower statistics but a higher S/B ∼ 1/85.
Techniques making use of boosted topologies have been shown to offer sensitivity improvements in the HH → 4b final state in HH [20] and HHjj [38] production, by making use of lower background yields in the high-m HH tail of the distribution and the ability to separate signal from background using boosted reconstruction techniques. In particular the sensitivity of HHjj measurements to modifications of c V V HH is greatly enhanced in the VBF topology. We have checked that the m HH distribution after all other cuts in the Z → ll analysis ‡ itself only contains new information in the low m HH region compared to the dominant Zbbbb background when taking the kinematic dependence of c 3 and c V V HH into account, see Figure 4b. We note that the behaviour when the process is dominated by the two different types of vertices we consider can be explained with reference to the contributions detailed in Figure 1: at tree level the c V V HH and c 2 V V H contributions scale as ∼ v 2 /m 2 HH at large m HH , while the λ contribution scales as ∼ v 4 /m 4 HH . This explains the similarity in scaling in the differential cross section between the Standard Model and when it is dominated by the c V V HH term, and the faster falloff when it is dominated by the λ term. This would allow a multivariate analysis to improve our results for large modifications to λ, but is difficult to make use of in a cut-based analysis.

IV. W HH PRODUCTION
The production of a pair of Higgs bosons in association with a W boson, W HH, has the advantage of a larger cross section compared to ZHH: σ(W HH → lν l + 4b) 0.04 fb at NNLO QCD (l = e, µ) [42]. The main challenge is to distinguish this channel from tt + X backgrounds, where X can be jets (including bb), or a Higgs boson decaying to bb.
We follow a similar analysis strategy to [41], which consists of the following steps: 1. Require exactly one lepton inside |η| < 2.5 with p T > 25 GeV. 5. Require the 4 leading jets to be b-tagged using the definition and efficiencies detailed in Section II. 6. Require that these b-tagged jets have the kinematics of a HH pair decay, χ HH < 1.6.
Similar to the Z → νν analysis we also attempt to reduce the top quark backgrounds using the observable χ t defined above. However, since the V HH kinematics are very similar between V HH channels, we find the same situation as in Figure 3b and, consequently, we do not employ χ t here either.

V. PROJECTED LIMITS
The visible cross sections after each selection in all three analyses are presented in Table I. The projected 95% confidence level limits on the trilinear Higgs coupling λ with the full HL-LHC data set using the three analyses are presented in Figures 5 and 7a: these have been derived by calculating the visible V HH cross section dependence on variations of λ (taking the kinematic dependence correctly into account) and fitting this to a polynomial, while the projected cross section constraints are evaluated using Poissonian likelihoods. Results with 20% systematic uncertainty approximated as a Gaussian uncertainty on the background are also presented to provide an estimate of the sensitivity of these results to realistic experimental conditions. Due to the weak sensitivity of the W HH analysis we do not present results including this systematic uncertainty for this channel. ‡ The kinematics of the Higgs boson system in all three channels are similar, so this conclusion also applies to the Z → νν and W → lν analyses. depend on the analysis to allow for the possibility of leptons escaping detection, and to ensure there are two leptons for the Z → ll analysis. The ZHH, W HH, Ztt, Ztt, W tt, ttH, and tt samples are generated at NLO QCD. The tt sample is further reweighted to the NNLO+NNLL QCD cross section [71]. The other samples are generated at leading order and the ttbb sample is reweighted to the NLO QCD cross section [72,73].  Our projected 95% confidence level (C.L.) limits on c 3 , λ = λ SM (1 + c 3 ), calculated using the CLs method [74], for the three analyses are presented in Table II. We also use the same strategy to project constraints on the quartic interaction V V HH and present the results in Figures 6 and 7b, with the 95% confidence level limits on c V V HH , g V V HH = g V V HH,SM (1 + c V V HH ), in Table III. Compared to the projected limits on g V V HH in the HHjj study in [36], and in particular to the projected limits set using boosted reconstruction techniques in [38], we see that the sensitivity of the V HH production channels is quite limited. Nevertheless, the ability to separate the ZZHH and W W HH contributions is as noted above unique to the V HH channels and these could therefore, at least in theory, provide complementary information for this measurement. The weak sensitivity here however suggests we are constraining values which violate perturbative unitarity at the scales the LHC probes: recasting the limits derived on dimension-6 operators in [75] we find that values for c V V HH ∼ O(5) violate perturbative unitarity for √ s ∼ TeV. All of the results presented here assume an idealised detector with perfect performance except for the b-tagging (false) rates given in Section II. To estimate how sensitive our results are to the uncertainties introduced by a realistic detector we have also performed the same analyses using the fast simulation machinery included in Rivet, which allows for lepton efficiencies and jet, lepton, and missing transverse energy smearing according to reported values by ATLAS to be taken into account [76,77]. The full results are presented in Table IV in the Appendix. The most significant effect we find is that top quark backgrounds become even more problematic for the Z → νν analysis. This reduces S/B by another factor of about 10.

A. Avenues for improvement
In light of the extremely challenging nature of this analysis that is evident from the results presented in Table I, we discuss some of the improvements that we have attempted. As is evident from Figure 3a, there is additional information to be used in the missing transverse energy distribution for the Z → νν analysis. We have checked explicitly that using a cut of 150 GeV leads to a small but noticeable improvement in the sensitivity. However the improvement vanishes when using a fast detector simulation due to the smearing of the |E miss T |. We have included a shape comparison for the smeared |E miss T | distribution in the Appendix. The plain tt background could also be reduced by using an operating point for the ATLAS MV2c10 b-tagger with higher background rejection, however this is not a "silver bullet" since the ttbb background remains sizeable and we have checked that using the most aggressive point in [66] does not lead to significant change in sensitivity due to the accompanying loss of signal efficiency.
As already discussed for the Z → ll analysis, making use of lower background yields in the high-m HH tail of the distribution does not generate dramatic improvements. The situation for the two other analyses where tt + X backgrounds are dominant could be more promising, however any such information will be difficult to make use of effectively at the HL-LHC due to the tiny cross sections.
A multivariate approach could make use of additional information in angular distributions as discussed in Section III B, potentially combined with observables like m ll and χ t motivated by top-based backgrounds to significantly improve the S/B without losing signal statistics in a way not possible in a simple cut-based analysis as presented here. However even so, the low signal rate after the basic selections we apply here will remain problematic: we have re-calculated the limit on c 3 from the Z → ll analysis presented in Table II assuming backgrounds are reduced by an order of magnitude with no cost in signal efficiency, and find −17.5 < c 3 < 13.5 in the absence of systematics, which illustrates that even large improvements in a multivariate analysis are unlikely to make this channel sensitive to the c 3 range where the cross section scaling favours it over other channels, Figure 2.
Accessing the full phase space information through the Matrix Element method [78][79][80][81][82][83][84] or by using MadMax [85][86][87] or similar would ultimately be necessary to make a definitive statement on the potential sensitivity of an idealised analysis, but due to the complex nature of the relevant backgrounds, this is beyond the scope of our present study. The extremely low signal rates which are evident from our simple sensitivity analysis here suggests that a differential analysis would be challenging even at a 100 TeV collider where signal cross sections are only ∼20 times larger.

VI. CONCLUSIONS
Interference between the trilinear coupling-induced and and other contributions in di-Higgs production could allow subdominant channels to have competitive sensitivity to specific modifications of λ. We have investigated V HH production at the HL-LHC, where this channel shows the strongest scaling with λ for small positive modifications. Additionally, this channel could potentially be used for measuring ZZHH and W W HH vertices separately. Unfortunately top quark pair production backgrounds are difficult to control for the Z → νν and W → lν decay modes, and the signal cross section is small for the Z → ll decay mode, making the ultimate sensitivity limited. We have discussed possible avenues to improve the analysis but concluded that using V HH production as a probe of modifications of the Higgs sector is likely to remain challenging even when employing more advanced techniques.

VII. ACKNOWLEDGEMENTS
We thank Christoph Englert and Adam Falkowski for helpful comments on an earlier draft of the manuscript. We would also like to thank Kalliopi Petraki and Andy Buckley for helpful discussions. KN acknowledges support by the NWO Vidi grant "Self-interacting asymmetric dark matter". AP acknowledges support by the ERC grant ERC-STG-2015-677323.

VIII. APPENDIX
The smeared missing tranverse energy distribution taking detector effects into account in the Z → νν analysis is presented in Figure 8. The cross sections at different stages in the cutflow when using the fast detector simulation are provided in Table IV.  IV: Cross sections in picobarns for ZHH (Z → ll and Z → νν) and W HH (W → lν), and backgrounds after the selections described in the text are applied when using the fast detector simulation included in Rivet. Generation level cuts on the invariant masses of lepton pairs, missing transverse energy, and pT of jets are employed for some of the backgrounds. Top quark branchings depend on the analysis to allow for the possibility of leptons escaping detection, and to ensure there are two leptons for the Z → ll analysis. The ZHH, W HH, Ztt, ttH, W tt and tt samples are generated at NLO QCD. The tt sample is further reweighted to the NNLO+NNLL QCD cross section [71]. The other samples are generated at leading order and the ttbb sample is reweighted to the NLO QCD cross section [72,73].