Sensitivity to New Physics in final states with multiple gauge and Higgs bosons

We analyse the sensitivity to beyond-the-Standard-Model effects of hadron-collider processes involving the interaction of two electroweak and two Higgs bosons, VVHH, with V being either a W or a Z boson. We examine current experimental results by the CMS collaboration in the context of a dimension-8 extension of the Standard Model in an effective-field-theory formalism. We show that constraints from vector-boson-fusion Higgs-pair production on operators that modify the Standard Model VVHH interactions are already comparable with or more stringent than those quoted in the analysis of vector-boson-scattering final states. We study the modifications of such constraints when introducing unitarity bounds, and investigate the potential of new experimental final states, such as ZHH associated production. Finally, we show perspectives for the high-luminosity phase of the LHC.


Introduction
Scattering processes involving the production of multiple electro-weak (EW) as well as Higgs (H) bosons are central to the entire physics programme of the Large Hadron Collider (LHC). Such reactions typically feature rich and complex signatures through which it is possible to probe several complementary properties of the underlying dynamics at the same time, allowing a detailed scrutiny of the gauge interactions they stem from. This can be exploited both for a precise extraction of fundamental Standard Model (SM) parameters, as well as to cast stringent bounds on potential beyond-the-Standard-Model (BSM) scenarios.
In the SM, the underlying SU (2) L ⊗ U (1) Y gauge symmetry fully dictates the pattern and relative intensity of the triple VVV and quartic VVV V gauge interactions (with V ( ) = W ± , Z), as well as of the VVH and VVHH Higgs-gauge interactions. Experimental measurement of such trilinear and quartic pure-gauge and Higgs-gauge couplings, which we collectively denote with TCs and QCs respectively, is thus especially suited to gather insight on the nature of the EW symmetry-breaking mechanisms, and to unveil potential deviations from the SM couplings in the EW gauge sector [1][2][3][4][5][6][7].
BSM scenarios can be conveniently analysed in a model-independent manner by resorting to an effective-field-theory (EFT) framework, in which indirect New Physics effects are parametrised by including in the Lagrangian higher-dimension operators that modify the production pattern of SM particles. In the Standard Model Effective Field Theory (SMEFT) [8,9], operators affecting both TCs and QCs start at dimension 6 [8,10,11], and their effects are commonly constrained at the LHC by measuring the production of two EW bosons [12][13][14][15][16][17][18], or the vector-boson-fusion production of EW bosons [18][19][20][21], even if weaker constraints can also be obtained from more complex processes such as vector-boson scattering (VBS) [22]. In particular, in a dimension-6 linear EFT extension of the SM VVH and VVHH vertex factors are totally correlated. Modifications of QCs with no TC contamination are instead more subtle as, assuming EW symmetry is linearly realised, they only arise in operators of dimension 8 or higher, which are dubbed for this reason genuine anomalous quartic operators [23,24]. As a consequence, they mainly affect the production of three bosons or the VBS di-boson channels. Effects of Higgs compositeness, encoded at low energy in an extension of the SM like the Higgs EFT (HEFT) [25,26], do also induce genuine QC variations, and their relation with dimension 8 or higher SMEFT operators has been studied in Refs. [27,28].
In this article we concentrate on the analysis of genuine SMEFT anomalous quartic operators, and set new constraints on their effects by means of a thorough investigation of the VVHH interaction 1 . The production processes we consider are Higgs-pair production, both in vector-boson fusion mode (VBF-HH) and in association with a Z boson (ZHH), as well as the loop-induced production of a single Higgs in association with two Z bosons (ggZZH). We study these reactions by means of a simplified version of the analysis that the CMS experiment employs to constrain genuine QC operators in VBS. After validating this procedure against public CMS QC constraints [29][30][31][32] in absence of unitarity bounds, we show that constraints from VBF-HH on genuine QC operators are comparable with or more stringent than those quoted in the analysis of the VBS signature. Similar investigations are performed on the ZHH and ggZZH channels, both at current LHC luminosity as well as for the high-luminosity-LHC (HL-LHC) phase. We do not consider in this study the qqinitiated ZZH process. At the LHC its cross section, known at NLO in the SM [33,34], is larger than that of the gluon-fusion process, however it is not within current data sensitivity and so no results from LHC collaborations exist. As our study constrains operators that specifically modify QCs without altering TCs, we only focus on processes exhibiting quartic vertices at the lowest order.
We choose the subset of possible dimension-8 operators which modify VVHH vertex strengths following Refs. [23,24]. Among these operators, the ones containing no EW field strength tensor are categorised as scalar, and denoted with a subscript 'S', while mixed operators (subscript 'M') feature two covariant derivatives of the Higgs doublet and two field strengths. We do not consider transverse ('T') operators in our analysis, as their effect on the VVHH vertex is found to be negligible, compatibly with the results of Ref. [35].
By analysing data at a centre-of-mass energy of 13 TeV, the CMS collaboration has determined limits on the Wilson coefficients f X /Λ 4 (Λ being the New Physics scale) of the above operators from several VBS final states, which are identified by the presence of two vector bosons (V or γ) in the final state and two jets with a large invariant mass and rapidity difference. The ATLAS collaboration has performed several VBS measurements at 13 TeV as well, but it has not extracted limits on dimension-8 operators so far. A full review of VBS measurements can be found in Refs. [36,37] and a compilation of 1 We neglect all contributions from SMEFT operators of dimension 6, in order to directly compare our results with the procedure adopted by CMS. A more refined analysis should consistently include constraints on dimension-6 operators, which is beyond the scope of our work.
constraints on dimension-8 operators is available in [38]. The most stringent limits on the coefficients of the dimension-8 operators f S0 , f S1 , f M0 , f M1 , and f M7 are derived from the simultaneous study of leptonic WZ and same-sign WW VBS [29], and from semileptonic WV VBS searches, where the hadronic decays of W and Z are not resolved [30]. The latter work also contains semileptonic ZV searches, which have weaker sensitivities to all operators. The most stringent limits quoted by the CMS collaboration on the remaining operators f M2−5 result from VBS Wγ and Zγ analyses [31,32]. However, we will show that such constraints can be improved by investigating these operators in the aforementioned VV analyses. The CMS collaboration does not quote any limit on the coefficient f S2 , while limits on f M6 do not add further information, as the corresponding operator was found to be redundant [39]. Concerning multi-Higgs-boson final states, both the ATLAS and CMS collaborations have searched for Higgs-pair production in the VBF-HH mode, in the 2b2b2j, bbτ + τ − 2j, and bb2γ2j channels [40][41][42][43]. In these studies, BSM effects are described in a different parameterisation which is equivalent to inserting an effective VVHH vertex modifier κ 2V where the SM corresponds to κ 2V = 1 [44]. For this study, we consider the constraint −0.1 < κ 2V < 2.2 at the 95% confidence level (@95% CL) set by the CMS 2b2b2j analysis [40]. 2 For other tri-boson final states that we will consider in this work no experimental studies exist. In this work, approximate estimates of the sensitivities are performed for the ZZH and ZHH final states only, since tri-boson processes containing simultaneously W ± bosons and H → bb decays are contaminated by large-rate tt(+X) backgrounds and cannot therefore be estimated correctly in a simplified analysis.
The typical effect of QC modifications is to distort differential spectra, primarily the invariant mass of the produced multi-boson system (e.g. m V V for di-bosons), with respect to the SM baseline, leading to enhanced rates in the high-energy tails. By increasing QCs, scattering amplitudes grow up to the point in which they violate unitarity bounds, signalling the breakdown of a truncated EFT approach, and the necessity of including explicit New Physics degrees of freedom in the Lagrangian. Unitarity constraints based on a partial-wave decomposition of matrix elements [24] are then important to guarantee the reliability of the exclusion bounds cast by the experimental searches, although CMS EFT-based analyses do not follow a uniform prescription to enforce them. In particular in [30] unitarity effects are neglected, and the limits are set using events with arbitrarily large values of m V V . In [31,32] no upper m V V bound or signal hypotheses are applied on the data, but the operator-dependent unitarity limits of [24] are computed a posteriori based on the measured limits on Wilson coefficients: since the signal model extends however beyond this limit, this may lead to a circular argument. Finally, in [29] physical limits are obtained by cutting the EFT contribution off at the unitarity bound, and removing all but the expected SM contribution for higher values of m V V ('clipping' method [46]). This approach is also questionable, as the signal model features an unphysical sharp drop at the unitarity bound, which depends on the value of the Wilson coefficient under examination. Also, more importantly, experimental data beyond the bound are still used to set the limit, stretching the EFT approach outside its validity domain. Studies on the uncertainty related to the unitarisation methods are also available [47].
In this article we adopt a dedicated technique to incorporate unitarity effects on EFT operators in a physically consistent manner, similarly to what done in Refs. [48][49][50]. First, limits on a given EFT operator are set as functions of m V V neglecting unitarity, by using the simulated cumulative m V V distributions (as opposed to total cross sections within fiducial cuts) for the analysed channels. The obtained m V V -dependent exclusion is then compared to the unitarity constraints for that operator, which are naturally in the form |f X | < f max,X (m V V ), and the best attainable limits are defined by considering the experimental data in the maximum m V V range in which the effects of the operator do not violate unitarity.
The structure of the paper is as follows: in section 2 we introduce the employed inputs from experimental analyses and the setup of our simulations, with particular reference to the inclusion of unitarity bounds; sections 3 and 4 show the validation against the CMS analysis of VBS; in sections 5 and 6 results are presented relevant for the VBF-HH, ZHH, and ggZZH channels at the LHC, with projections for HL-LHC detailed in section 7; we finally conclude in section 8.

Simulation setup and observables
The simulation of the final states under examination is performed using version 2.7.3 of the MadGraph5_aMC@NLO event generator [51,52] at leading order (LO) in QCD 3 at the parton level. LHC conditions are reproduced using symmetric proton-proton collisions at √ s = 13 TeV. We consider the following processes which are sensitive to dimension-8 EFT effects: 1. EW production of two same-sign W bosons and two quark jets (W ± W ± VBS); 2. EW production of a W and a Z boson, associated with two quark jets (W ± Z VBS); 3. EW production of two opposite-sign W bosons and two quark jets (W + W − VBS); 4. EW production of two Higgs bosons and two quark jets (VBF-HH); 5. EW production of two Higgs and a Z boson (ZHH); 6. Loop-induced (LI) production of two Z and a Higgs boson from a gluon-gluon initial state (ggZZH). Figure 1 shows typical lowest-order diagrams for the processes outlined above, with (left) and without (right) EFT insertions. Since signal process 5 has not been searched for at the LHC so far, we also simulate the main expected background for this process when the dominant H → bb decay channel is chosen to reconstruct Higgs bosons: 7. QCD-mediated production of a Z boson with two b and two anti-b quarks.   Table 1. LHC processes of interest studied in the present work, numbered according to the main text. The number of active flavours considered in protons and jets is also listed. σ[m min , m max ] is the cross section in the interval m min ≤ m ≤ m max , with m being the invariant mass of the produced di-or tri-boson final state. Each cross-section value is evaluated after cuts on final-state jets described in the text and reports in brackets the integration uncertainty on the last digit.
simulations share the same choice for SM EW parameters: The parton distribution functions used in the simulations are taken from the NNPDF2.3 LO set [57] with α s (m Z ) = 0.130. EFT effects are simulated by loading in MadGraph5_aMC@NLO suitable UFO models [58], containing the relevant modifications to SM vertex structures. For all signal LO simulations, the latest models available at [59] were used. It must be noticed that these differ in several aspects from those used in the CMS analyses [29,30], which date back to 2013. Taking into consideration observations in recent papers such as [39], the list and the content of EFT operators were revised by removing redundant operators (e.g. M6), redefining others (e.g. M5) and adding new ones (e.g. S2), in order to obtain a complete basis.
For the loop-induced ggZZH simulation, a dedicated UFO model has been constructed, since the original one did not support loop diagrams. Starting from a NLO QCD description of the SM, new Lorentz structures corresponding to the additional operators were inserted. This was possible since none of these operators involves coloured particles. While in this paper such a model has been used only for the loop-induced process ggZZH, it could be exploited to refine our analysis including NLO QCD corrections for all processes, which is left for future work. In the VBF-HH samples, modifications in terms of κ 2V were instead performed by changing the WWHH and ZZHH couplings by the same scaling factor. Our simulations include: • for processes 1-3 (VBS) and 5-6 (ZHH, ggZZH), sets of samples with one Wilson • for process 4 (VBF-HH), sets of samples with coefficient values as above, and additionally samples with an effective vertex factor κ 2V = 0, 1, ±2, ±5, ±10; • for process 7 (background), just the SM generation.
Since the EFT-sensitive region is located at high scattering energies, the decays of vector and Higgs bosons are expected to produce secondary particles with large transverse momentum (p T m H /2, m V /2). For this reason we do not apply parton showering to parton-level events, assuming modifications to the transverse momentum p T to be moderate. Following the same reasoning, no selections are applied to the decay products of the gauge and Higgs bosons, since we assume that at large p T the detector acceptance be constant over the considered phase space (and similar for signals and physical backgrounds): since the experimental inputs already have acceptance effects folded in, these factorise in our approach. The only exception is for processes 5 and 7, where no experimental analysis exists and the approach followed in this case is discussed in Section 6.
In the VBS and VBF processes (1-4), on the other hand, typical experimental selections on the additional jets are applied: In the background process (7) a cut on pairs of b-jets is applied to emulate a possible experimental selection. It is required that the invariant mass of at least two of the possible bb pairs, which can be obtained from combining the four b-jets, satisfies the requirement 115 GeV < m bb < 135 GeV. Examples of simulated di-boson or tri-boson invariant-mass spectra are shown in Figs. 2 and 3. We notice that, while final states with Higgs bosons have SM cross sections which are a few orders of magnitude smaller with respect to VBS, EFT effects which affect treelevel amplitudes induce very large enhancements at high invariant mass, especially in the case of mixed M-type operators. In all cases the ratios BSM/SM in the lower panels of Figs. 2 and 3 tend to unity at low invariant mass: this is an indirect confirmation that the modifications induced by dimension-8 operators on SM structures, like EW boson propagators, are numerically negligible for the purpose of our analysis.
CMS uses Monte Carlo templates to perform binned data analyses and extract limits on the Wilson coefficients of dimension-8 operators. Typical data distributions used with this aim are the di-boson invariant mass, or the transverse mass if the final state contains one or more undetected neutrinos whose transverse momenta are inferred from the event imbalance in the transverse plane. This procedure requires a full description of the SM backgrounds as well as a precise simulation of detector effects. In order to obtain a simplified estimate of the relative EFT sensitivities for different final states, we use as a physical observable the cross section (or equivalently the expected number of observed events at given integrated luminosity) in the interval m min ≤ m ≤ m max , with m being the invariant mass of the produced di-or tri-boson states in the various channels; we denote such a cross section as σ[m min , m max ], or simply as σ.
Unitarity bounds for all operators are computed as functions of f X /Λ 4 using Ref. [24], resulting in more stringent bounds than those provided by the VBFNLO program [60], which only considers s-wave scattering. In all analyses presented in the following, m min is fixed to a default value of 1.1 TeV. The choice is a compromise between defining an EFT-enriched region where possible SM backgrounds are negligible, and having a m min value below the unitarity bound in the whole range of Wilson coefficients considered in the analysis. The values of σ using m min = 1.1 TeV and m max = √ s (namely with no upper bound) are shown in the rightmost column of Table 1. Repeating the analysis with m min shifted by ±100 GeV modifies the outcome of the validation performed in Section 3 by less than 10%, for all operators.

Validation on VBS without unitarity constraints
As a first step we show that, although experimental constraints result from template analyses of the observable final states, constraints purely based on σ[m min , m max ] are a reliable proxy of the actual EFT sensitivity. Since most CMS results do not use unitarity regularisation, in this validation phase we consider m max = √ s, i.e. no upper bound on the system invariant mass.
The σ values are computed for each process and simulated sample. The corresponding cross sections for a given experimental final state are obtained by multiplying σ by the branching fractions (BF ) of the gauge and Higgs bosons involved. For instance, the cross section for the 2b2b2j final state considered in [40] is obtained as σ(VBF-HH) × BF 2 (H→ bb). For semileptonic WV VBS [30], the BF -weighted cross sections of VBS W + W − , VBS W ± Z and VBS W ± W ± are summed, as they all contribute to the observed final state. For the ggZZH and ZHH processes, the H→ bb mode is chosen, while only Z decays to electrons and muons are considered, since hadronic decays are more challenging experimentally and beyond the scope of this analysis.
Cross sections for each channel are computed as functions of f X /Λ 4 , and quadratic fits are performed on the obtained results. Examples of such parabolic fits are shown in Fig. 4 for two of the processes under consideration. The validation proceeds through the steps below.
1. A CMS experimental analysis is selected, and the published 95% CL exclusion limit on a randomly-chosen operator coefficient is considered.
2. The exclusion limit is imposed on the parabola corresponding to the chosen operator, yielding a 95% CL exclusion limit on σ.
3. Upon applying such a σ-based exclusion on the parabolae corresponding to all other operators, exclusion limits on the corresponding coefficients are in turn determined. 4. As a closure test, the limits obtained with this procedure are systematically compared with the CMS published ones.

5.
Steps 1-4 are then repeated for different choices of the initial input operator, in order to check that similar bounds are obtained.
The results of the validation phase are presented in Table 2. We first remark that CMS published results are rather incomplete for both the leptonic W ± Z and semileptonic WV analyses, where several operators for which sensitivity is significant were not examined. As anticipated, this sensitivity is in fact comparable to or better than those quoted in Vγ VBS analyses [31,32], which therefore we do not try to reproduce here.
Secondly, we observe that validation is everywhere successful except for the S0-2 operators in the WV semileptonic analysis where validation results are about 60% more pessimistic with respect to the published values. Several tests were performed by changing the relative weight of the processes (especially of W + W − which has a larger background in the experimental analysis) but in none of them closure was found. We notice, however, that in the CMS paper details of the simulation are not reported, in particular it is not mentioned whether the latter was performed including b jets or not. We point out that, when including b jets, other processes are generated by the syntax 3 in Table 1 which are not VBS processes, most notably electroweak top production via qq → γ * → tt → W + W − bb, and enhance the cross section by a factor of ∼ 4. While these processes can be suppressed experimentally via b-jet vetoes, there can be a double-counting effect which artificially leads to lower expected (and observed) limits if tt backgrounds are simulated separately and subtracted. In conclusion, we do not consider our validation successful in this case and cautiously use CMS published limits for scalar operators.
We point out that, as mentioned in item 5 of the above list, the validation procedure has been repeated using different operators as input, and all limits obtained upon input variation have been verified to be essentially identical to those quoted in Table 2. On one hand this confirms the robustness of the whole validation procedure, and on the other hand it allows us to quote limits for the input operators of Table 2

Implementation of unitarity regularisation in VBS
In order to consider unitarity bounds in the limits, we adopt a different implementation of the clipping method used in the CMS analysis [29], following a similar approach as in Refs. [48][49][50]. In our approach, we evaluate σ[m min , m max ] for several values of m max , varying between m min + 100 GeV and the maximum kinematically allowed mass. In an experimental analysis, this approach would be equivalent to not taking into consideration data or simulated events above m max in the measurement. For each of these values of σ, the procedure used for validation and described in Sec. 3 is followed to obtain m maxdependent limits on operator coefficients. Since only part of the experimental data fall into the selected [m min , m max ] invariant-mass intervals, the 95% CL exclusion limits on the cross section determined at step 2 are rescaled in each test, considering that measurement uncertainties at high mass are statistically dominated, and therefore poissonian errors can be assumed.    Table 2   become less stringent at smaller m max because of the reduced data statistics as well as of the weaker dependence of σ on the operator coefficient in the restricted intervals. In Fig. 5 the unitarity bounds derived from [24] are also shown. The intersection of the experimental and the unitarity-bound curves represents the maximum invariant mass that can be used to set experimental limits which can unambiguously be interpreted in an EFT expansion without violating unitarity. Hence, the corresponding limit on f X /Λ 4 is the best attainable limit which can be set using a specific amount of data for the process under consideration. We notice that in all cases these limits are significantly less stringent than the ones obtained by neglecting unitarity effects; in some specific cases, not shown, the two curves do not cross at all, meaning that the available data are not sufficient yet to set more stringent limits than those imposed by unitarity alone. In Table 3 we present the best limits obtained with the outlined procedure for each of the considered processes and operators. Corresponding CMS limits using the modified clipping method are quoted in Ref. [29] for the fully-leptonic results only. Considering the substantial difference in the unitarity regularisation procedure and in the source of unitarity bounds (which is the VBFNLO program [60] in the case of Ref. [29]), it is not surprising that our results display discrepancies with respect to the limits quoted in [29], although we note that the constraint-reducing factors when imposing unitarity are of similar magnitude in the two analyses. In particular, [29] quotes some very loose constraints on specific operators, which cannot be investigated with our choice of m min , as they are weaker than the unitarity bounds. These are anyway of modest interest, as coefficient values in such loose ranges would produce too large EFT corrections compared to the SM.

Results with VBF-HH
As mentioned before, experimental results for VBF-HH are given in terms of an effective VVHH vertex factor, κ 2V , rather than as limits on EFT Wilson coefficients. We employ an analogous procedure as the one used for the validation, namely: 1. The most stringent published 95% CL exclusion limit on κ 2V is considered.
2. We exploit the VBF-HH simulation as a function of κ 2V to impose an exclusion limit on the corresponding parabola, yielding a 95% CL exclusion limit 4 on σ.
3. Upon applying such a σ-based exclusion on the parabolae corresponding to all EFT operators, exclusion limits on the corresponding coefficients are determined.
4. As a validation test, by repeating the procedure using one of the f X coefficients as input, we check that the CMS limits on κ 2V are reproduced. Figure 6 presents σ fits as a function of κ 2V and of one EFT Wilson coefficient. In Table 4 we present the related upper limits without unitarity regularisation, compared to the best limits found using the VBS analysis in Section 3. In spite of the much smaller yields, we find that, already with the available LHC data, VBF-HH estimated limits supersede those obtained with VBS for f M0 , f M2 and f M3 , and are comparable for most of the Mtype operators. In an experimental interpretation of the results, limits on M-type Wilson coefficients could be combined to yield more stringent limits than those obtained by any single analysis. An analogous procedure is followed to obtain limits taking unitarity into account. Figure 7 presents m max -dependent limits along with theoretical unitarity bounds. In Table 4 we also present the related upper limits with unitarity regularisation, compared to the best limits found using the VBS analysis in Section 4. The conclusions of the previous paragraph remain valid with unitarity as well.

New experimental final states
The ggZZH and ZHH processes, not yet studied in dedicated experimental analyses, are considered in this paper. The purpose is to investigate the potential sensitivity of such processes by performing an exploratory feasibility study. Table 4. Comparison between VBF-HH constraints and those obtained from semileptonic VBS with or without unitarity constraints, with limits extracted from current CMS Run2 analyses. All entries represent 95% CL lower and upper limits in units of TeV −4 . Bars represent processes for which there is no sensitivity to the corresponding operator, or cases where the theoretical unitarity bound is more stringent than the experimental one for all m max cut-off values.  For ggZZH, we generated samples as described in Section 2, considering EFT modifications of the SM cross section. The cross section computed for this process is very small in the SM (see Table 1): in Figure 8 we show the cross section for different values of the EFT coefficients considered. We find that, even for very significant variations (±40 TeV −4 ) of the EFT Wilson coefficients, the predicted cross section remains extremely small, thus making this process not sensitive enough to be observed at the LHC, not even in the highluminosity phase. Therefore we refrain from furthering the study of this channel in the present paper. We stress however that it is possible to simulate this process (as well as other loop-induced or even NLO-accurate processes) in a consistent manner, employing the NLO UFO model we have constructed including the dimension-8 operators from Ref. [23]. 40 20 0   For the ZZH process, the computed cross section is found to be larger than that of ggZZH in the SM (Table 1), making it worth to estimate the sensitivity that can be obtained from a benchmark analysis of this process. In this analysis, we consider the LHC Run2 luminosity and we estimate the number of detectable events as N = σ · L · ε · A, where A denotes the experimental acceptance for the Z and H decay products and ε is the corresponding total detection efficiency. We compute values of σ for the signal and backgrounds described in Section 2, by assuming the Higgs boson decaying into a bb pair, and the Z bosons decaying into leptons (electrons and muons). We perform Z and H decays using Pythia8 [61] and define typical LHC acceptance requirements as: for all decay products, simulating a realistic experimental analysis. The acceptance is computed to be 93% for signal (in the SM, conservatively taken as equal for all EFT scenarios) and 71% for background. The overall efficiency of identification and selection of electrons, muons, and b jets is estimated from experimental papers [40,62] to be 70%. From the derived estimates of the signal and background yields, we compute 95% CL upper limits on the cross section of the ZHH process using the Feldman-Cousins approach [63], and then we use this result to evaluate the limits on EFT operators with a procedure analogous to those presented in the previous sections, taking into account unitarity regularisation. In Table 5 we present the values of the limits obtained with such procedure. Since no valid limits are found in presence of unitarity with Run2 luminosity, we report as reference the limits obtained without unitarity requirements. We show in the  next section that the ZHH process becomes sensitive to EFT variations when considering an increased integrated luminosity.

Perspectives for HL-LHC
In spite of the moderate increase in centre-of-mass energy, it is interesting to explore projections for the HL-LHC dataset, as the increased number of events has a non-trivial impact on operator limits when unitarity bounds are accounted for. In absence of unitarity regularisation, and assuming poissonian scaling in statisticallydominated high-mass regions, 95% CL exclusion limits on σ are expected to improve as L −1/2 , with L being the integrated luminosity. Since the dependence of the cross section on the Wilson coefficients is quadratic, the related f X /Λ 4 limits scale as L −1/4 , yielding an improvement of a factor about 2.2 when passing from the 140 fb −1 of the LHC Run2 to the 3000 fb −1 expected from the HL-LHC 5 .
In Table 6 we present the expected limits in a HL-LHC scenario, with and without unitarity regularisation, for VBF-HH compared to the best limits found using the VBS analysis of Section 3. Limits in absence of unitarity are obtained by simply rescaling the excluded cross sections by the square-root ratio of the integrated luminosities, hence the sensitivity hierarchy in this case is the same as in the Run2 study. Conversely, limits in presence of unitarity have a significant gain from an increased dataset since, while 95% CL limits become tighter, the m max value for which the result remains EFT-interpretable correspondingly moves to larger values, allowing more data to be included in the sensitivity estimate. This is due to the fact that a higher integrated luminosity enhances the sensitivity to deviations from the SM that are generated by smaller operator coefficients, resulting in a wider invariant-mass window allowed by unitarity constraints. This effect would lead to Table 6. Comparison of VBF-HH constraints with those obtained from VBS when applying or not unitarity constraints, in a projection with full HL-LHC luminosity. All entries represent 95% CL lower and upper limits in units of TeV −4 . Bars represent processes for which there is no sensitivity to the corresponding operator, or cases where the theoretical unitarity bound is more stringent than the experimental one for all m max cut-off values.
the first physical limit on f S1 , as well as to improvements on other limits by factors of up to 4-5, as shown for f M3 in Figure 9 (left).   [24] (red lines). The intersection of the experimental and the unitarity-bound curves represents the best limits which can unambiguously be interpreted in an EFT expansion without violating unitarity. While all individual results scale roughly as L −1/4 , the improvement on the EFT-interpretable limit is a factor of 4.6. Table 7 reports the 95% CL limits obtained for the ZHH process in the HL-LHC scenario, with and without unitarity regularisation. As anticipated in Section 6, this process becomes sensitive to dimension-8 operators when increasing the luminosity. In this scenario, in fact, it is possible to set limits on some M-type operators even in presence of unitarity  bounds, as shown in Figure 9 (right). We note that such limits are not quite competitive, due to the small cross section of the ZHH process with respect to its background process: for future analyses, it will be crucial to develop strategies to enhance the signal and reduce the number of background events. As a simple test, by halving the b-jet pair invariantmass window for the background, described in Section 2, we observe an improvement of about 20% on the limits both with and without unitarity regularisation, stemming from the decrease in the background contribution.

Summary
In this article we have studied the sensitivity to BSM effects of the interaction between two electroweak gauge bosons (V) and two Higgs bosons (H). To this aim, we have focused on hadronic reactions that feature the VVHH vertex at lowest order, in particular VBF Higgs-pair production, ZHH associated production, and loop-induced gluon-fusion ZZH production. We have considered modifications to these processes induced by dimension-8 EFT operators affecting the VVHH interaction, and set up a simplified framework to cast stringent bounds on the Wilson coefficients of these operators using such reactions. After validating our framework against published experimental limits, we have performed a systematic comparison of our results with current bounds from the CMS collaboration, based on the VBS signature. This has shown that, as far as dimension-8 VVHHmodifiers are concerned, the VBF-HH channel is able to yield limits which are comparable to, or even more stringent than, those stemming from VBS, thus proving viable in view of a full-fledged experimental analysis, as well as of an experimental combination with VBS limits. The ZHH and ggZZH channel have been instead shown to have more limited constraining power on the same operators with Run2 luminosity.
We have then considered the effect of unitarity constraints on the extraction of exclusion limits in an EFT-based analysis. As this subject is typically treated in a non-systematic manner in current extractions, we have employed a dedicated prescription to consistently take unitarity bounds into account. We find that, as expected, unitarity constraints systematically weaken the experimental limits on EFT Wilson coefficients, including the ones set by current CMS analyses of the VBS channel. However, even in presence of unitarity, VBF-HH limits are equally competitive with VBS-based ones.
We have finally investigated the perspectives of the exclusion limits for the highluminosity phase of the LHC, showing that limits in presence of unitarity can improve at the HL-LHC by a factor of up to 5 with respect to Run2, which is significantly more than what is expected from pure luminosity-scaling considerations. In this scenario, the ZHH final states also can contribute in a combined exclusion of some of the relevant EFT coefficients.