Off-shell Higgs production at the LHC as a probe of the trilinear Higgs coupling

In the context of the Standard Model effective field theory (SMEFT) we examine the constraints on the trilinear Higgs coupling that originate from off-shell Higgs production in proton-proton collisions. Our calculation of the $gg \to h^\ast \to Z Z \to 4 \ell$ process includes two-loop corrections to gluon-gluon-fusion Higgs production and one-loop corrections to the Higgs propagator and its decay. Employing a matrix-element based kinematic discriminant we determine the reach of LHC Run 3 and the high-luminosity option of the LHC in constraining the relevant SMEFT Wilson coefficients. We present constraints that are not only competitive with but also complementary to the projected indirect limits that one expects to obtain from inclusive measurements of single-Higgs production processes at future LHC runs.


Introduction
In the pursuit of observing or constraining deviations from the Standard Model (SM) picture, measuring the couplings of the Higgs boson to the other bosons, fermions and itself lies at the heart of the physics programme at the Large Hadron Collider (LHC). The couplings of the scalar resonance to electroweak (EW) gauge bosons and third-generation fermions have been tested extensively, showing agreement with those of the SM Higgs boson at the level of (10 − 20)% [1,2]. The strength of these couplings will be further scrutinised at the high-luminosity option of the LHC (HL-LHC) and possible future colliders. In contrast to the couplings of the Higgs boson to EW gauge bosons and third-generation fermions, the Higgs couplings to light matter fields and its self-interactions are largely unexplored at present (cf. for example [3] for an overview).
Within the SM, the trilinear and quartic couplings of the physical Higgs field h are parametrised by the scalar potential where m h 125 GeV and v 246 GeV are the Higgs mass and vacuum expectation value (VEV), respectively, and the trilinear (λ) and quartic (κ) Higgs couplings obey the relation The SM potential is thus fully determined by only two parameters, namely v and λ. However, many beyond the SM (BSM) scenarios allow for deviations of the Higgs self-couplings with respect to their SM predictions (a comprehensive collection of such theories can be found in the white paper [4]) and, consequently, could imply a departure from the SM relation (1.2). Measuring or constraining the Higgs self-couplings is therefore essential to our understanding of the mechanism of electroweak symmetry breaking (EWSB) and furthermore provides a way to probe the existence of new physics. The aforementioned relatively loose constraints on the Higgs self-interactions are due to the smallness of the cross sections of double-Higgs and triple-Higgs production, the go-to observables for testing the Higgs potential because of their direct sensitivity to λ and κ. As the triple-Higgs cross section is O(0.1 fb) for proton-proton (pp) collisions at a centre-of-mass energy of √ s = 14 TeV even at the end of HL-LHC, having an integrated luminosity of 3 ab −1 at hand, the quartic Higgs coupling will remain unexplored (cf. [5][6][7][8][9][10][11][12][13][14] for the prospects to determine the quartic Higgs coupling at future colliders). The HL-LHC prospects for double-Higgs production are considerably better but still remain challenging as the cross section is O(33 fb) in pp collisions at √ s = 14 TeV. As a result, only O(1) determinations of the trilinear Higgs coupling from double-Higgs production seem to be possible at the LHC -see [15,16] for the latest prospect studies by ATLAS and CMS. With this in mind, other methods of constraining the trilinear Higgs coupling have been proposed in recent years. The indirect approach first outlined in [17], where the sensitivity to λ arises from loop corrections to the process e + e − → Zh, was later extended to the Higgs production and decay modes relevant at the LHC [18][19][20][21][22][23], at lepton colliders [24,25] and to EW precision observables [26,27].
Analyses of the constraints on the trilinear Higgs coupling including indirect probes have been presented in [21,[28][29][30] and recently also by ATLAS and CMS [1,[31][32][33]. Generally, the direct constraints on λ obtained through double-Higgs production were shown to furnish the most stringent bounds, but indirect constraints from single-Higgs production processes have the potential to be competitive or could fulfil at least a complementary role. As already pointed out in the works [3,19,20], including measurements of differential distributions of single-Higgs processes could turn out to be crucial due to their non-trivial dependence on λ. This point was further investigated in the article [21], where it was found that the associated production of the Higgs together with a EW gauge boson (V h) or a top-antitop pair (tth) provide additional sensitivity to λ at the differential level.
In this work, we investigate the LHC Run 3 and the HL-LHC sensitivity to modifications of the trilinear Higgs coupling from off-shell Higgs production in the gluon-gluon fusion (ggF) channel, considering decays to four charged leptons, i.e. h * → ZZ → 4 . Within the context of the SM effective field theory (SMEFT), we study the effects of a modified trilinear Higgs coupling through a shape analysis of differential distributions. In particular, the use of a matrix-element based kinematic discriminant to improve the sensitivity in regions where the Higgs boson is off-shell will allow us to put constraints on the trilinear Higgs coupling that are not only competitive with but also complementary to constraints from inclusive single-Higgs production cross-section measurements at future LHC runs.
This article is organised as follows. In Section 2, we the introduce the effective interactions that are relevant for the computations performed in our paper. The calculations of the loop corrections to the production, propagator and decay of the Higgs boson in the process gg → h * → ZZ → 4 are presented in Section 3. Our numerical analysis is performed in Section 4, where we also present our LHC Run 3 and HL-LHC projections for the constraints on the trilinear Higgs coupling and compare them to previous results. We conclude and present an outlook in Section 5. Additional material is relegated to Appendix A and Appendix B.

Parametrisation of BSM effects
To allow for a model-independent analysis, we work in the framework of the SMEFT [34][35][36] to parametrise the possible BSM physics entering the trilinear Higgs coupling. While efforts are being made to obtain a general fit of the many dimension-six SMEFT operators to the available experimental data (see for example [37,38] for some recent results), requiring absolute model-independence can in practice impede the extraction of meaningful results due to the large number of free parameters in the fit. In this article, we therefore consider only the subset of dimension-six CP-even operators in the so-called Strongly Interacting Light Higgs (SILH) basis [39] that are build purely from SM Higgs doublets H and derivatives: (2.1) Here we have used the short-hand notation H † We thus implicitly assume that the couplings of the Higgs boson to EW gauge boson and fermions are SM-like and focus our attention on the self-interactions contained in (1.1). Note that the redundant operators O R = |H| 2 |D µ H| 2 and O D = |D 2 H| 2 do not appear in (2.1) as they can be removed in the full SILH basis via an appropriate redefinition of the Higgs field or equivalently its equations of motion [40]. Furthermore, as the operator O T does not modify the trilinear Higgs coupling and is moreover severely constrained through measurements of the ρ parameter describing the degree of custodial symmetry violation [41], it is irrelevant for the purpose of this work.
Within these restrictions, it is sufficient to supplement the SM Lagrangian with the effective operators O 6 and O H only, resulting in the following effective Lagrangian 3 Description of the gg → h * → ZZ calculation In this section we briefly describe the calculation of the O(λ) corrections to the process gg → h * → ZZ that arise in the context of (2.2). The corrections associated to insertions of the operator O 6 are illustrated in Figure 1. They fall into three classes: (i) corrections to ggF Higgs production, (ii) corrections to the Higgs propagator and (iii) corrections to the Higgs decay. In the following, we will discuss separately each of these three ingredients as well as their implementation in our Monte Carlo (MC) code.

Higgs production
The O(λ) corrections to ggF Higgs production receive contributions from both two-loop topologies (see the upper Feynman diagram in Figure 1) as well as from the wave-function renormalisation of the Higgs boson field. The relevant renormalised vertex that describes the process g(p 1 ) + g(p 2 ) → h(p 1 + p 2 ) can be written as [22] Γ µν ggh (p 1 , where α s = g 2 s /(4π) denotes the strong coupling constant, a 1 and a 2 are colour indices and η µν = diag (1, −1, −1, −1) is the Minkowski metric. In addition is the one-loop correction to the Higgs boson wave function associated to insertions of the operator O 6 [18,19] and represents the one-loop triangle diagram with internal top quarks and m t is the top-quark mass. Hereŝ = 2p 1 · p 2 and the C 0 function denotes a three-point Passarino-Veltman scalar integral for which our definition follows the conventions used in the LoopTools package [42]. The non-factorisable two-loop form factor F 2 has been calculated analytically in [19,22] using the method of asymptotic expansions, which in this case is valid up to energies √ŝ m t . To cover the full off-shell range of interest up to √ŝ = 1 TeV, we employ the numerical results for the non-factorisable two-loop form factor F 2 presented in [11,12] in our numerical analysis performed in Section 4.

Higgs propagator
The Higgs propagator also receives corrections from insertions of the operator O 6 (see the centre graph in Figure 1). The resulting renormalised contribution to the self-energy of the Higgs takes the formΣ where the one-loop corrections to the Higgs wave function has already been given in (3.2) and the bare Higgs self-energy and the mass counterterm in the on-shell scheme are given by the following expressions Here the B 0 functions are two-point Passarino-Veltman scalar integrals defined as in [42].

Higgs decay
The full O(λ) correction to the Higgs decay h → ZZ receives a two-loop contribution (see the lower diagram in Figure 1) as well as a counterterm contribution involving Higgs wave-function renormalisation. In the notation of [20], the relevant renormalised vertex describing the h( where m Z denotes the Z-boson mass. The O(λ) corrections to the one-loop form factors G 1 and G 2 are given by and the tensor coefficients C 1 , C 00 and C 11 of the three-point Passarino-Veltman integrals are defined as in the publications [20,42].

MC implementation
The three different types of O(λ) corrections affect not only the overall size of the gg → h * → ZZ → 4 cross section, but also modify the shape of kinematic distributions such as the four-lepton invariant mass m 4 . In order to be able to predict gg → h * → ZZ → 4 we have implemented thec 6 corrections arising from ggF Higgs production, the Higgs propagator and the Higgs decay to Z bosons into version 8.0 of MCFM [43]. Our implementation includes all contributions up to O(λ 2 ) that arise from squaring the full gg → h * → ZZ matrix element (ME) which comprises both the BSM graphs depicted in Figure 1 as well as the LO SM Feynman diagram shown on the left-hand side of Figure 2. We also note that the contribution to the wave-function renormalisation constant (3.2) coming from the propagator corrections at O(λ) exactly cancels against those of the vertices when combined

Numerical analysis
In this section, we study the impact of the O(λ) corrections to the gg → h * → ZZ → 4 process on the m 4 spectrum and a ME kinematic discriminant to be defined later. Following a brief discussion of the impact of QCD corrections, these results are then used to perform a sensitivity study of the two-dimensional constraints on the Wilson coefficientsc 6 andc Hsee (2.2) -that can be obtained at LHC Run 3 and the HL-LHC from measurements of pp → ZZ → 4 . The bounds from Higgs off-shell measurements are finally compared to the limits that are expected to arise from a combination of inclusive single-Higgs measurements at the end of LHC Run 3 as well as the HL-LHC era.

Modifications of differential distributions
In Figure 3 we show m 4 distributions for the Higgs channel alone (left panel) and for the Higgs channel, the gluon continuum background and their interference combined (right panel) at LO in QCD. An example of a one-loop Feynman diagram that contributes to the SM gg → ZZ background is shown in the centre of Figure 2. We restrict ourselves to the off-shell region by considering a mass window of 220 GeV < m 4 < 1000 GeV.
The leptons ( = e, µ) are required to be measured in the pseudorapidity range |η | < 2.5 and the lepton with the highest transverse momentum (p T ) must satisfy p T, 1 > 20 GeV while the second, third and fourth lepton in p T order is required to obey p T, 2 > 15 GeV, p T, 3 > 10 GeV and p T, 4 > 6 GeV, respectively. The lepton pair with the mass closest to the Z-boson mass is referred to as the leading dilepton pair and its invariant mass is required to be within 50 GeV < m 12 < 106 GeV, while the subleading lepton pair must be in the range of 50 GeV < m 34 < 115 GeV. Similar cuts are employed in the ATLAS and CMS analyses [44][45][46][47][48][49][50]. The input parameters used throughout our work are given by G F = 1/( √ 2v 2 ) = 1.16639 · 10 −5 GeV −2 , m Z = 91.1876 GeV, m h = 125 GeV and m t = 173 GeV. The shown spectra assume pp collisions at √ s = 14 TeV and employ NNPDF40_nlo_as_01180 parton distribution functions (PDFs) [51] with the renormalisation and factorisation scales µ R and µ F dynamically, i.e. for each event, set to m 4 . Our predictions include both the different-flavour (i.e. e + e − µ + µ − ) and the same-flavour (i.e. 2e + 2e − and 2µ + 2µ − ) decay channels of the two Z bosons. Two features of the distributions shown on the left-hand side of Figure 3 deserve a further discussion. First, below the two-Higgs production threshold at m 4 = 2m h the BSM gg → h * → ZZ → 4 spectra are both visibly smaller that the SM prediction. This feature can be understood by noting that for sufficiently large values of the Wilson coefficientc 6 the terms (3.5) provide the dominant O(λ) corrections. In fact, the Higgs propagator corrections necessarily reduce the real part of gg → h * → ZZ → 4 amplitude for m 4 < 2m h , and this destructive interference can be so pronounced that the BSM contribution almost exactly cancels the SM contribution. Second, above the two-Higgs production threshold the bare Higgs self-energy Σ(ŝ) develops an imaginary part because both internal Higgs lines in the propagator correction (see the centre graph in Figure 1) can go on their mass shell. As a result the BSM gg → h * → ZZ → 4 spectra can be larger than the SM prediction for m 4 > 2m h . Notice that the two aforementioned features are clearly visible in the case of the BSM distribution corresponding toc 6 = 20, while forc 6 = 10 the spectrum is not enhanced above the two-Higgs production threshold because the O(λ) corrections to Higgs production and decay (see the upper and lower graph in Figure 1) are numerically relevant in this case and tend to cancel the effect of the Higgs propagator.
Our results for the gg → ZZ → 4 distributions including the Higgs signal, the continuum background and their interference are displayed in the right panel of Figure 3. In the vicinity of the two-Higgs production threshold m 4 = 2m h one observes a plateau-like structure in both BSM spectra. This feature arise from the combination of the modified Higgs signal and the interference of the BSM signal with the continuum SM background. This atypical shape change provides a genuine probe of loop corrections to the Higgs propagator involving light virtual particles. Such corrections arise in the case at hand from the insertions of the operator O 6 but they can also appear in ultraviolet (UV) complete models of BSM physics (see for instance [52][53][54][55]). Both BSM spectra also show an enhancement in the m 4 tail. Notice that the observed shape change is qualitatively different from the relative modifications that arise from tree-level insertions of dimension-six SMEFT operators which typically show a roughly quadratic growth with m 4 (see for example [53,56]). Similar statements apply to the case when the SM Higgs boson width is rescaled in such a way that the pp → h * → ZZ → 4 cross section close to the Higgs peak is unchanged [57][58][59]. See Appendix A for a detailed discussion.
We have seen that the inclusion of the O(λ) corrections to the pp → h * → ZZ → 4 amplitude associated to insertions of the SMEFT operator O 6 lead to phenomenologically relevant kinematic features in the m 4 distribution. The analysis sensitivity to the Higgs channel, especially in the off-peak region, has been shown to benefit considerably from the use of ME-based kinematic discriminants (see for instance [44][45][46][47][48][49][50][59][60][61][62]) to separate the gg → h * → ZZ → 4 signal from the main SM background coming from ZZ production in qq-annihilation. A relevant Feynman diagram that contributes to the qq → ZZ background at LO in QCD is displayed on the right in Figure 2. Being sensitive not only to m 4 but to another seven variables such as the invariant masses of the two opposite-sign lepton pairs (for details consult [60][61][62]), the ME-based discriminants fully exploit the event kinematics. In practice, the ME-based discriminants are often embedded in a multivariate discriminant based on a boosted decision tree (BDT) algorithm, but as it turns out in the case of the four-lepton final state the sensitivity of the BDT analysis improves only very little with respect to the ME-based discriminant alone (for example in the case of the analysis [45] the improvement amounts to a mere 2%). In the following, we restrict ourselves to an approach with only a ME-based discriminant, which we define as follows [45,46,48,50] Here P h denotes the squared ME for the gg → h * → ZZ → 4 process, P gg is the squared ME for all gg-initiated channels (including the Higgs channel, the continuum background and their interference) and P qq is the squared ME for the qq → ZZ → 4 process. Following the publications [45,46,48,50] the constant c is set to 0.1 to balance the qq-and gg-initiated contributions. We add that in the SM more than 99% of the pp → ZZ → 4 cross section fall into the range of −4.5 < D S < 0.5 [45]. The kinematic discriminant (4.1) thus presents a null test for BSM models that lead to events with D S < −4.5 or D S > 0.5.
To illustrate the discriminating power of the variable D S we show in the left panel of Figure 4 the normalised SM distributions at LO in QCD for the three contributions corresponding to the MEs that enter (4.1). The discriminant D S , which is calculated for every event in the simulation, is implemented in MCFM which uses the gg-initiated MEs provided in [59]. One observes that the distribution corresponding to qq → ZZ → 4 peaks at D S −3, while the gg → h * → ZZ → 4 spectrum is shifted to higher D S featuring a maximum at around D S −0.5. An enhancement of the gg → h * → ZZ → 4 amplitude in the off-shell region will hence lead to a D S distribution of the full pp → ZZ → 4 process that is shifted to the right compared to the SM spectrum. For BSM scenarios that predict an enhancement in the tail of the m 4 distribution, one thus expects to find an excess of events for D S −1. As discussed in Appendix A, this is precisely what happens if the Higgs boson couplings and its total decay width are modified according to (A.1).
In the case of the O(λ) corrections to off-shell Higgs production resulting from the insertions of the SMEFT operator O 6 , the modifications of the normalised D S distributions are more intricate. For the full gg → ZZ → 4 contribution this is shown on the right-hand side of Figure 4. One observes that in the case ofc 6 = 10 the BSM spectrum is shifted to lower values compared to the SM result. This is a simple consequence of the fact that for c 6 = 10 the gg → h * → ZZ → 4 amplitude and therefore P h in (4.1) is reduced for m 4 values below the two-Higgs production threshold (cf. the left panel in Figure 3). Forc 6 = 20 one instead sees that the D S distribution has two maxima: one at D S −5 and another one at D S −3. The peak at D S −5 (D S −3) is associated to kinematic configurations that lead to a reduction (an enhancement) of P h for m 4 < 2m h (m 4 > 2m h ) -see again the left plot in Figure 3. Finally, forc 6 = 30 the BSM contributions always enhance P h and hence the D S distribution is shifted to larger values. Notice that the sharp cut-off of the SM distribution at D S −3.5, makes (4.1) more sensitive to kinematic configurations that reduce P h rather than enhance it. In fact, this feature turns out to be the key element that allows to set powerful constraints on BSM scenarios withc 6 = 0 by means of the ME-based discriminant D S .

Impact of QCD corrections
In the following we discuss the possible impact of higher-order QCD corrections to the pp → ZZ → 4 cross section differential in the ME-based kinematic distribution introduced in (4.1). In fact, as the ggF contribution to pp → ZZ → 4 is loop-induced (see the centre graph in Figure 2), it enters the ZZ production cross section at O(α 2 s ), i.e. at the next-to-next-to-leading order (NNLO) in QCD. State-of-the-art predictions for fourlepton production at the LHC, obtained at NNLO in QCD [63][64][65][66] and matched to parton shower [67,68], are reaching an impressive accuracy of O(2%) for inclusive cross sections and O(5%) in the case of differential distributions.
While a precision phenomenological study of ZZ production a la [63][64][65][66][67][68][69][70] including the O(λ) corrections associated to the insertions of O 6 is beyond the scope of this work, we wish to assess at least approximately the impact of higher-order QCD corrections on our analysis, in particular the effects on the D S spectra. To this purpose we proceed in the following way. We first calculate for each production channel the so-called K-factor defined as the ratio between the fiducial cross section at a given order in QCD and the corresponding LO QCD prediction. In the case of the gg-initiated contribution we employ the results of the recent work [68], where one of us reported NLO QCD corrections to the corresponding four-lepton invariant mass spectrum (see Figure 2 of that article). One observes that the ratio between the NLO and LO ggF predictions is essentially flat in the region of interest for this work, i.e. for values of the invariant mass within 220 GeV < m 4 < 1000 GeV. By averaging over the ratio of the NLO and LO m 4 spectra within the aforementioned m 4 window we find K NLO gg = 1.83, which is in line with the previous works [66,69,70]. In the case of the qq-initiated contribution we utilise the LO and NNLO QCD results obtained in [66] using MATRIX [71]. The relevant K-factor again turns out to be basically flat in m 4 , with a central value of K NNLO qq = 1.55, which is in accordance with [63].  for each channel and QCD order. The numbers for the gg-initiated channel are obtained utilising results presented in [68], while the results for the qq-initiated channel are taken from [66]. Consult the main text for further details.
The aforementioned K-factors are listed in Table 1, along with the scale uncertainties in each production channel at the relevant order in QCD. They are then used to obtain a QCD-improved prediction for the D S distributions, in the following way: Here the label LO indicates that both the cross sections as well as the ME-based discriminant D S are calculated at LO in QCD. Clearly, (4.2) only captures part of the higher-order QCD corrections to the dσ/dD S spectra, namely those that are associated to the differential cross sections dσ, but ignores beyond-LO effects to dD S . To improve upon this approximation one would need to extend the calculation of the ME-based discriminant (4.1) to the NLO in QCD. While achieving NLO accuracy is in principle possible for ME-based discriminants [72][73][74], the actual calculations are in practice complicated by the fact that they require the use of modified jet algorithms to map resolved and unresolved parton configurations onto the proper MEs -see also [75] for a discussion of this point. Since NLO and NNLO QCD corrections to the shapes of kinematic distributions in pp → ZZ → 4 are small and often indistinguishable when compared to the associated theoretical uncertainties, it is expected that the LO discriminant D S as defined in (4.1) maintains its discriminating power beyond the well-defined LO. This renders the procedure (4.2) a pragmatic and simple approach to incorporate higher-order QCD effects.
In the left plot of Figure 5 we compare the D S spectrum for pp → ZZ → 4 in the SM obtained at LO to the QCD-improved D S spectrum (4.2). One observes a close to flat K-factor of around 1.6 between the LO and the improved prediction of the D S spectrum for pp → ZZ → 4 . The uncertainty bands represent the scale uncertainties obtained by applying the numerical values listed in Table 1 to the relevant the LO and QCD-improved spectra. We observe that the inclusion of higher-order QCD corrections reduces the scale uncertainties by a factor of about 3 from roughly (7 − 8)% to about (2 − 3)%, in line with recent precision SM phenomenological studies. However, the fact that the central value of  Table 1. Right: D S spectra for pp → ZZ → 4 production in the SM (dashed green) and forc 6 = 10 (solid red) andc 6 = 15 (solid blue). The thin green band corresponds to the scale uncertainties of the QCD-improved SM predictions for D S , while the wide green band represents half of the relative difference between the QCD-improved and the LO SM predictions for ME-based discriminant. The lower panel shows the ratio of the BSM and SM predictions including the aforementioned uncertainty bands. All distributions assume pp collisions at √ s = 14 TeV and the improved SM as well as the BSM predictions have been obtained by means of (4.2). Additional explanations can be found in the main text. the improved prediction lies well outside the LO uncertainty bands in all bins demonstrates that at LO, scale variations do not provide a reliable way to estimate the size of higher-order QCD effects. In fact, similar issues are known to occur also beyond LO, for example in the case of inclusive ggF Higgs production (see [76][77][78] for the corresponding state-of-the-art SM calculations) or even for four-lepton production [68], where the loop-induced gg contribution entering at NNLO is unaccounted for by the NLO scale uncertainties. With this in mind and given that the kinematic discriminant D S is only LO accurate, we wish to take a rather conservative approach and estimate the theoretical uncertainties to be half of the relative difference between the QCD-improved and the LO predictions, which corresponds to a relative uncertainty of about 30%.
On the right-hand side of Figure 5 we compare two BSM D S distributions to the SM prediction, where the scale uncertainties quoted in Table 1 and the aforementioned conservative estimate are shown in the lower panel as the solid green and dot-dashed green bands around the SM spectrum, respectively. The depicted spectra have all been obtained using (4.2). In accordance with the general discussion in Section 4.1, one observes that the BSM spectra deviate most significantly from the SM distribution for D S −3.5, while for D S −3.5 the deviations are generically small. In fact, for the two choices ofc 6 shown in the figure the resulting modifications of the D S spectrum for D S −3.5 lie well within our conservative theoretical uncertainty band. Furthermore, since the deviations for D S −3.5 are associated to the modification of the off-shell Higgs production cross section in the region m 4 < 2m h , they will lead to a detectable change in events for a sufficiently large |c 6 |. For instance, forc 6 = 10 (c 6 = 15) the shift in the pp → ZZ → 4 cross section restricted to D S < −3.5 amounts to around 0.4 fb (1.2 fb) compared to the SM. This corresponds to around 1200 (3700) additional pp → ZZ → 4 events at the HL-LHC. Notice finally that the relative modifications of the D S spectrum due to insertions of O 6 are much larger than the shifts seen in the m 4 distribution (cf. the right panel in Figure 3). The ME-based discriminant (4.1) therefore provides a significantly better sensitivity to BSM models withc 6 = 0 than the m 4 spectrum.

Constraints on Wilson coefficientsc 6 andc H
Below we determine the constraints on the Wilson coefficientsc 6 andc H that future LHC runs may be able to set. In the case of the constraints arising from off-shell Higgs production, we utilise the QCD-improved D S predictions obtained by (4.2), assuming a detection efficiency of 99% (95%) for muons (electrons) that satisfy the event selections described at the beginning of Section 4.1. These efficiencies correspond to those reported in the latest ATLAS analysis of off-shell Higgs production [50]. The statistical uncertainties of the computed D S distributions are determined per bin assuming Poisson statistics, i.e. taking the statistical error to be √ N i with N i the number of events in a given bin i. The largest systematic uncertainties in our analysis arise from the theoretical uncertainties on the gg → h * → ZZ → 4 signal process, the gg → ZZ → 4 and qq → ZZ → 4 background processes and the interference between the gg-initiated signal and background. For the theoretical uncertainties on the improved D S prediction (4.2) we take the conservative estimate discussed in Section 4.2, in which we assume also PDF uncertainties at the level of ±5% are included. In our LHC Run 3 analysis we will thus use a total theoretical uncertainty of ±30% when determining the bounds onc 6 andc H . Anticipating theoretical advances in the case of the HL-LHC we assume that the theoretical uncertainties due to scale variations and PDFs are reduced to ±15% which does no seem unrealistic. In fact, similar assumptions are made in the HL-LHC studies [46,47]. Compared to the theoretical uncertainties the experimental uncertainties of systematic origin are close to negligible as they amount to O(2%) (see for example [79]). We will thus ignore the experimental systematics in what follows.
The total number of events in bin i depends on the Wilson coefficientsc 6 andc H in the following way where N i (c 6 ) denotes the number of events which we calculate using MCFM correcting them by QCD effects (4.2) and lepton efficiencies. Notice that N i (0) corresponds to the SM expectation of events. The significance Z i is calculated as a Poisson ratio of likelihoods modified to incorporate systematic uncertainties on the background using the Asimov approximation [80,81]: Here s i (b i ) represents the expected number of signal (background) events in bin i and σ b i denotes the standard deviation that characterises the systematic uncertainties of the associated background. To set bounds onc 6 andc H we assume that the central values of a future measurements of the D S distribution will line up with the SM predictions. We hence employ where ∆ i denotes the relative total systematic uncertainty in bin i. We will employ binindependent systematic uncertainties of ∆ i = 0.3 and ∆ i = 0.15 at LHC Run 3 and HL-LHC, respectively. The total significance Z is obtained by adding the individual Z i values in quadrature. Parameter regions with a total significance of Z > √ 2 erf −1 (CL) are said to be excluded at a given confidence level CL. Here erf −1 (z) denotes the inverse error function. In our numerical analysis, we include 29 bins of equal size of 0.2 that cover the range −5.1 < D S < 0.5.
Before deriving the projected bounds on the Wilson coefficientc 6 from off-shell Higgs production, we recall that the currently best LHC 95% CL limit reads [33]   See also [29,30]. The bound (4.6) has been obtained from a combination of ten double-Higgs and single-Higgs production measurements performed by the ATLAS collaboration. Assumingc H = 0 and employing the fit strategy described above, we find the following 95% CL bounds from future hypothetical measurements of off-shell Higgs production:  minor role. This feature allows to set limits like (4.7) that are competitive with (4.6) using solely off-shell Higgs production in the ggF channel.
On the right in Figure 6 we furthermore show the 68% CL (dashed) and 95% CL (solid) constraints in thec 6 -c H plane that derive from our fit employing LHC Run 3 (blue) and HL-LHC (red) data. One observes that the bounds onc H are essentially independent of the precise value ofc 6 if the latter Wilson coefficient is sufficiently small. In the case ofc 6 = 0, we obtain at 95% CL for instancē We emphasise that non-zero values ofc H do not change the shape of the D S spectrum but only its normalisation cf. (4.3) . This feature explains the approximatec 6 -independence of the exclusion contours in thec 6 -c H plane for smallc 6 . Notice finally that the bounds onc H are significantly more stringent than those onc 6 . This is expected because the SMEFT operator O H changes the prediction for pp → h * → ZZ → 4 already at LO (i.e. one loop) while the corrections due to O 6 start at NLO (i.e. two loops).

Comparison to bounds from inclusive single-Higgs production
To further demonstrate the benefits of off-shell Higgs production in setting bounds on the Higgs trilinear coupling, we compare the results obtained in the previous section to the projected constraints one expects to obtain from inclusive single-Higgs measurements at future LHC runs. In the inclusive case the O(λ) corrections to the various Higgs production and decay channels can be written in terms ofc 6 andc H as 9) where N h has been defined already in (3.2) and C Γ h 1 = 0.23·10 −2 [19,20,22]. Notice that the Wilson coefficientc H leads to a universal correction to all Higgs decay channels. Therefore it leaves the Higgs branching ratios unchanged. The calculations needed to obtain the process-dependent factors C σ i 1 and C Γ f 1 have been performed in [18][19][20][21][22][23]. In our numerical analysis of the inclusive single-Higgs observables we include ggF, W h, Zh, vector-boson fusion (VBF) and tth production and consider the Higgs-boson branching ratios to pairs of photons (γγ), EW bosons (W + W − , ZZ), bottom quarks (bb) and tau leptons (τ + τ − ). The associated C σ i 1 and C Γ f 1 coefficients are collected in Table 2. In terms of (4.9), keeping only terms linear in λ, the Higgs signal strengths for production in channel i and decay in channel f can be written as which we use to build the following χ 2 function: (4.11) Here we have assumed that the central values of the future measurements of the Higgs signal strengths will coincide with the corresponding SM predictions. The variables ∆ f i encode the relative total uncertainties obtained by combining the theoretical and statistical uncertainties associated to µ f i . We collect the values of the ∆ f i used in our LHC Run 3 and HL-LHC analyses in Table 3. Notice that the LHC Run 3 numbers are obtained by combing the current theoretical and the statistical uncertainties in quadrature, while the HL-LHC numbers assume that all theory uncertainties are halved with respect to our current understanding of the relevant signals and backgrounds. The latter assumption corresponds to the scenario S2 in the ATLAS paper [82]. The allowed CL regions are then obtained by minimising (4.11) and determining the solutions to ∆χ 2 = χ 2 − χ 2 min < 2Q −1 (1/2, 1 − CL) with Q −1 (a, z) the regularised incomplete gamma function.
In the left (right) panel of Figure 7 we show the projected 68% and 95% CL constraints in thec 6 -c H plane for the LHC Run 3 (HL-LHC). One observes that the LHC Run 3 fit to the 15 inclusive single-Higgs observable listed in Table 3 shows a pronounced flat direction in thec 6 -c H plane. To understand this feature one first has to realise that only the processindependent coefficients C σ i 1 and C   Table 3: Relative total uncertainties ∆ f i on the Higgs signal strengths defined in (4.10). The quoted LHC Run 3 and HL-LHC numbers are taken from [83] and [82], respectively. Further information can be found in the main text.
The relatively large coefficient C σ tth 1 (cf. Table 2) plays a particularly important role in this respect. Given the large total uncertainties of tth, h → γγ, ZZ at LHC Run 3, the constraining power of the tth channels and thus the impact of C σ tth 1 is however limited. As a result, the inclusive LHC Run 3 exclusions are mainly determined by the ggF channels that have a flat direction forc 6 andc H satisfying µ f ggF δσ ggF (c 6 ,c H ) 0. The situation is visibly improved in the inclusive HL-LHC fit, mostly because the total uncertainties of the tth channels are expected to be significantly reduced. From both panels in Figure 7  is however also evident that the flat direction in the inclusive fit is strongly broken by the constraints arising from off-shell Higgs production.
From the above it should be clear that inclusive single-Higgs and off-shell Higgs measurements should therefore be combined if one wants to exploit the full potential of the LHC in constraining the trilinear Higgs coupling through indirect probes. Performing such a combined analysis, we find forc H = 0 the following 95% CL limits We add that the bounds (4.12) and (4.13) depend in a non-negligible way on the assumed total uncertainties. In this respect one should remember that in the case of the constraints arising from off-shell Higgs production in ggF production we have assumed total systematic uncertainties of ±30% and ±15% in our LHC Run 3 and HL-LHC fit, respectively. We believe that these are conservative uncertainties -results for two additional more aggressive assumptions about the systematic uncertainties entering the HL-LHC off-shell Higgs analysis can be found in Appendix B. In fact, given the steady progress in the calculation of massive higher-loop corrections to pp → ZZ → 4 (see [84,85] for the latest theoretical developments) and in view of the fact that it is theoretically known of how to achieve NLO accuracy for ME-based discriminants [72][73][74], it should be possible to put our naive estimates of theoretical uncertainties on more solid grounds. As can be seen from the left panel in Figure 6, any improvement of our theoretical understanding of the D S distribution for D S −3.5 will notably increase the sensitive of off-shell Higgs measurements to modifications of the trilinear Higgs coupling.

Conclusions and outlook
In this article, we have studied the constraints on the trilinear Higgs coupling that originate from off-shell Higgs production in pp collisions at future LHC runs. To keep the discussion model-independent we have worked in the context of the SMEFT in which the renormalisable SM interactions are augmented by the dimension-six operators O 6 and O H cf. (2.1) . Our computation of the gg → h * → ZZ → 4 process includes two-loop corrections to ggF Higgs production and one-loop corrections to the Higgs propagator as well as the decay h * → ZZ. The resulting scattering amplitudes have been implemented into MCFM where they can be combined with the SM MEs for gg → h * → ZZ → 4 , gg → ZZ → 4 and qq → ZZ → 4 to obtain differential distributions for the full pp → ZZ → 4 process, including the corrections due to insertions of the SMEFT operator O 6 .
Using our MC implementation, we have then studied the shape differences in the fourlepton invariant mass m 4 distribution and the ME-based kinematic discriminant D S defined in (4.1) that are induced by modifications of the trilinear Higgs coupling. We found that the inclusion of BSM effects leads to phenomenologically relevant kinematic features in both spectra. In fact, the discriminant D S turns out to provide particularly powerful constraints on the Wilson coefficientc 6 of the pure-Higgs operator O 6 . The stringent constraints onc 6 arise because BSM scenarios withc 6 = 0 can lead to D S < −4.5, which provides a null test since 99% of the pp → ZZ → 4 events in the SM fall into the range −4.5 < D S < 0.5. We have also assessed the possible impact of higher-order QCD correction to D S , arguing that (4.1) maintains its discriminating power beyond the well-defined LO.
To demonstrate the benefits of off-shell Higgs production in setting bounds on the Higgs trilinear coupling, we have determined the constraints on the Wilson coefficientsc 6 andc H that the LHC with 300 fb −1 and 3 ab −1 of integrated luminosity at √ s = 14 TeV may be able to set. We have then compared the obtained LHC Run 3 and HL-LHC bounds to the projected constraints that a combination of inclusive single-Higgs measurements is expected to provide. Our analysis shows that ggF off-shell Higgs production allow to put constraints on the trilinear Higgs coupling that are not only competitive with but also complementary to the exclusion limits obtained from inclusive single-Higgs production. Specifically, we found that future studies of the D S distribution in pp → ZZ → 4 production should help to remove flat directions in thec 6 -c H plane that remain unresolved in fits that incorporate only inclusive single-Higgs measurements. By combining all single-Higgs boson measurement we find that at the LHC Run 3 (HL-LHC) it might be possible to constrain modifications of the trilinear Higgs coupling as parameterised by (2.4) to the 95% CL range c 3 ∈ [−4.0, 6.1] (c 3 ∈ [−1.7, 5.7]). Additional HL-LHC projections that employ two different assumptions about the systematic uncertainties entering our off-shell Higgs analysis can be found in Appendix B.
The studies performed in this article can be extended in several ways. First, to strengthen the constraints on the trilinear Higgs coupling derived in our work, one should also include in the projections measurements of double-Higgs production as well as EW precision observables. See [29,30,33] for such global analyses based on LHC Run 2 data. Second, the ME-based discriminant D S might also be a powerful tool to constrain BSM effects in pp → ZZ → 4 that arise in the context of Higgs portal models [53][54][55]. Given the limited direct LHC reach in such models [86,87], investigating the indirect sensitivity that off-shell Higgs production can provide seems like a worthwhile and timely exercise. Third, since any improvement of our theoretical understanding of the D S distribution (4.1) will have a tangible impact on the sensitivity of off-shell Higgs measurements to modifications of the trilinear Higgs coupling, we believe that theoretical activity in this direction is important. Since many ingredients (if not all) are already available in the literature on pp → ZZ → 4 production [63-70, 72-74, 84, 85] achieving such a goal is certainly possible.

Acknowledgments
We thank Darren Scott for useful discussions. The Feynman diagrams shown in this work were drawn with TikZ-Feynman [88].

A Higgs width effects
In this appendix we illustrates how rescalings of the form with g SM hXX and Γ SM h denoting the couplings and total decay width of the SM Higgs boson, respectively, modify the kinematic distributions in off-shell ggF Higgs production. Notice that (A.1) leaves the total Higgs production cross sections in all channels unchanged compared to their SM values. This is however not true for the off-shell Higgs cross sections that are essentially independent of Γ SM h and are thus modified if the Higgs couplings g SM hXX are rescaled as in (A.1). By measuring the total number of off-shell Higgs events one can therefore place indirect limits on the total width of the Higgs boson [44][45][46][47][48][49][50][57][58][59].
In Figure 8 we show our results for the m 4 distributions in the gg → ZZ → 4 channel (left) and the D S spectrum of pp → ZZ → 4 (right) for two different rescalings (A.1). The choice ξ = 3 and ξ = 1.5 thereby corresponds approximately to the present LHC Run 2 [48,49] and the projected HL-LHC [46,47] sensitivity, respectively. From the left plot one sees that compared to the SM the BSM predictions have larger off-shell Higgs cross sections with the relative difference between the spectra growing roughly linearly with m 4 . Notice that the observed shape changes are qualitatively different from the relative modifications that occur in the case of the O(λ) corrections associated to insertions of the operator O 6 as shown on the right-hand side in Figure 3. From the right plot in Figure 8 one furthermore observes that compared to the SM the BSM distributions of the ME-based discriminant are enhanced for D S −1. Since they do not feature the enhancements for D S −3.5, the shown D S spectra are hence distinct from the distributions that are displayed on the right in Figure 5, which correspond to the spectra resulting from insertions of the SMEFT operator O 6 . Notice that in contrast to the O(λ) corrections, the effects of (A.1) lead solely to enhancements in the tail of the m 4 distribution. In this case extra care is required in estimating the systematic uncertainties, because the NLO QCD corrections to the gg-induced channel included approximately by means of (4.2) implicitly assume an asymptotic expansion in the top-quark mass (cf. [67,68]) of the relevant two-loop gg → ZZ amplitudes. This expansion fails above the top-quark threshold, i.e. for four-lepton invariant masses m 4 > 2m t , which introduces compared to the case discussed in Section 4.2 an additional systematic uncertainty. To account for this issue, we instead of ±15% assume an enlarged total theoretical uncertainty of ±25%. Employing this uncertainty estimate and performing a shape-fit to the D S spectrum following the procedure outlined in Section 4.3, we obtain for 3 ab −1 of HL-LHC data the 95% CL bound Γ h < 1.49 · Γ SM h . This finding is in line with the limits reported in [46,47] which validates the used fitting approach.

B Additional HL-LHC projections
A crucial ingredient in the shape fit to the D S distribution described in Section 4.3 are the systematic uncertainties σ b i on the background as parametrised by the parameters ∆ i in (4.5). In this appendix we present results for two additional more aggressive assumptions about the systematic uncertainties entering the HL-LHC off-shell Higgs analysis. Specifically, we will employ the two different choices ∆ i = 0.08 and ∆ i = 0.04 of bin-independent systematic uncertainties. These choices can be motivated by recalling that the systematic uncertainties that ATLAS quotes in the HL-LHC study [82] for the on-shell gg → h → ZZ signal strength amount to 5.0% and 3.9% in the baseline scenario S1 and S2 for the expected total systematic uncertainties. The corresponding systematic uncertainties quoted in the CMS work [47] are 7.3% and 4.1%. Since the O(λ) corrections to D S considered in this work are associated to kinematic configurations with m 4 not far above 2m Z , it seems not unreasonable that theoretical predictions of the D S spectra can reach an accuracy that is very similar to the systematics that is expected to be achievable at the HL-LHC in the case of on-shell gg → h → ZZ production.
In the left (right) panel of Figure 9 we show the projected 68% and 95% CL HL-LHC constraints in thec 6 -c H plane assuming ∆ i = 0.08 (∆ i = 0.04). The constraints from inclusive single-Higgs probes (blue regions) are compared to the off-shell Higgs constraints (orange regions). Their combinations (red contours) are also displayed. From a combined analysis of inclusive single-Higgs and off-shell Higgs probes, we find forc H = 0 the following 95% CL limits By comparing the HL-LHC limits given in (4.12) and (4.13) to the above results, one observes that a reduction of systematic uncertainties has only a minor impact in the case ofc 6 , while it has a noticeable impact on the resulting bounds onc H . The limits (B.1) and (B.2) can also be translated into constraints on the modifications of the trilinear Higgs coupling as parameterised by (2.4). The corresponding 95% CL ranges read c 3 ∈ [−1.4, 5.6] and c 3 ∈ [−1.2, 5.4], respectively, which again represent only minor improvements compared to the HL-LHC bound derived in the main body of this article.