A global view on the Higgs self-coupling at lepton colliders

We perform a global effective-field-theory analysis to assess the precision on the determination of the Higgs trilinear self-coupling at future lepton colliders. Two main scenarios are considered, depending on whether the center-of-mass energy of the colliders is sufficient or not to access Higgs pair production processes. Low-energy machines allow for ~40% precision on the extraction of the Higgs trilinear coupling through the exploitation of next-to-leading-order effects in single Higgs measurements, provided that runs at both 240/250 GeV and 350 GeV are available with luminosities in the few attobarns range. A global fit, including possible deviations in other SM couplings, is essential in this case to obtain a robust determination of the Higgs self-coupling. High-energy machines can easily achieve a ~20% precision through Higgs pair production processes. In this case, the impact of additional coupling modifications is milder, although not completely negligible.


Introduction
So far, the LHC provided us with a good deal of information about the Higgs boson. The determination of its linear couplings to several Standard Model (SM) particles is nowadays approaching, and in some cases surpassing, the 10% precision, allowing for powerful probes of a broad class of natural beyond-the-SM (BSM) theories. On the contrary, the prospects for measuring the Higgs self interactions, namely its trilinear and quadrilinear self-couplings, are much less promising. At present, the trilinear Higgs coupling is loosely constrained at the O(10) level, and the high-luminosity LHC (HL-LHC) program could only test it with an O(1) accuracy (see for instance the experimental projections in Refs. [1,2]). The prospects for extracting the quadrilinear Higgs selfcoupling are even less promising.
From a theoretical point of view, on the other hand, the determination of the Higgs self-interactions is of primary importance. They characterize the Higgs potential, whose structure could shed some light on the naturalness problem. Moreover, they control the properties of the electroweak phase transition, determining its possible relevance for baryogenesis. Sizable deviations in the Higgs self-couplings are expected in several BSM scenarios, including for instance Higgs portal models or theories with Higgs compositeness. All these considerations motivate the effort spent investigating the achievable precision on the Higgs self-interactions at future collider experiments.
Projections for high-energy hadron machines (100 TeV pp colliders in particular) are already available in the literature [3]. They show that a very good precision on the determination of the trilinear Higgs coupling, of the order of 5%, is possible. High-energy hadron machines, however, might only be constructed in a distant future and could be preceded by lower-energy lepton colliders. It is thus worth studying the impact of future lepton machines on the determination of the Higgs potential. In this work, we perform such an analysis, providing an assessment of the achievable precision on the determination of the Higgs trilinear self-coupling.
We consider a comprehensive set of benchmark scenarios including low-energy lepton machines (such as FCC-ee and CEPC) as well as machines that can also run at higher energies (ILC and CLIC). We will show that low-energy colliders, although not able to access directly the Higgs trilinear coupling in Higgs pair production processes, can still probe it by exploiting loop corrections to single Higgs channels that can be measured to a very high precision. This approach, pioneered in Ref. [4], allows for a good determination of the Higgs trilinear interaction, which can easily surpass the HL-LHC one. In performing this analysis, however, one must cope with the fact that different new-physics effects may affect simultaneously the single Higgs cross sections, see Ref. [5] as well as Refs. [6][7][8][9][10][11][12][13][14]. In such a situation, a robust determination of the Higgs self coupling can only be obtained through a global fit that takes into account possible deviations in other SM couplings. We will show that, within the SM effective field theory (EFT) framework with a mild set of assumptions, the relevant operators correcting single Higgs production can be constrained provided enough channels are taken into account. In this way, a consistent determination of the Higgs self-coupling is possible even without direct access to Higgs pair production.
High-energy machines, on the other hand, are able to directly probe the trilinear coupling via Higgs pair production, through Zhh associated production and W W -fusion. We will see that these two channels provide complementary information about the Higgs self interaction, being more sensitive to positive and negative deviations from the SM value respectively. We will also show, as anticipated in Ref. [15], that a differential analysis of the W W -fusion channel, taking into account the Higgs-pair invariant mass distribution, can be useful to constrain sizable positive deviations in the Higgs trilinear coupling that are hard to probe with an inclusive study.
The paper is organized as follows. In Section 2, we discuss the indirect trilinear Higgs coupling determination through single-Higgs production processes. The impact of pair production is then studied in Section 3. The main results are summarized and discussed in Section 4 for the most relevant benchmark scenarios considered in the analysis. The appendices collect some useful formulae and provide additional results for some secondary benchmark scenarios not included in the main text. Additional numerical results are provided as ancillary files together with the arxiv submission of this paper.

Low-energy lepton machines
In this section, we study the precision reach on the trilinear Higgs coupling through the exploitation of single Higgs production measurements. These are the dominant handles available at future circular lepton colliders, like the CEPC and FCC-ee, which cannot easily deliver high luminosities at center-of-mass energies where the Higgs pair production rate becomes sizable. These machines could run above the e + e − → Zhh threshold, at a 350 GeV center-of-mass energy in particular, but the small cross section (in the attobarn range) and the limited integrated luminosity lead to a negligible sensitivity to this channel. The analysis of single-Higgs production can also be relevant for the ILC. While this machine could eventually reach a center-of-mass energies of 500 GeV (or even of 1 TeV) in a staged development, its initial low-energy runs can have an impact on the determination of the trilinear Higgs coupling that is worth investigating.
According to recent reports [16,17], both CEPC and FCC-ee are planned to collect 5 ab −1 of integrated luminosity at 240 GeV. FCC-ee is also envisioned to collect 1.5 ab −1 at 350 GeV. 1 Although a run at this center-of-mass energy is not officially forecast for the CEPC, it is nevertheless a viable option given its planned tunnel circumference of 100 km. As a general circular collider run scenario, we therefore consider the collection of 5 ab −1 of integrated luminosity at 240 GeV and several benchmark luminosities at 350 GeV, namely 0, 200 fb −1 and 1.5 ab −1 .
The full ILC run plan comprises 2 ab −1 of integrated luminosity at 250 GeV, 200 fb −1 at 350 GeV, and 4 ab −1 at 500 GeV, with these luminosities equally shared between runs with two P (e − , e + ) = (±0.8, ∓0.3) beam polarization configurations [19,20]. Additional results for a 70%/30% repartition of the luminosity between the P (e − , e + ) = (±0.8, ∓0.3) polarizations will be provided in Appendix B. In this section, we focus only on the runs at 240/250 GeV and 350 GeV, and consider a few benchmarks for the integrated luminosity collected at 350 GeV.
To summarize, we focus on the following benchmark scenarios: • Later in this section we also extend these scenarios to cover a continuous range of luminosities at 240 (250) and 350 GeV.

Higher-order corrections to single-Higgs processes
As a first step, we analyze how a modification of the trilinear Higgs coupling affects single-Higgs processes. We parametrize possible new physics effects through the quantity κ λ defined as the ratio between the actual value of the trilinear Higgs coupling λ 3 and its SM expression λ SM 3 (the Higgs vacuum expectation value is normalized to v = 1/( While the trilinear coupling does not enter single-Higgs processes at leading order (LO), it affects both Higgs production and decay at next-to-leading order (NLO). The corresponding diagrams for Higgsstrahlung (e + e − → hZ) and W W -fusion (e + e − → ννh) production processes are shown in Fig. 1. In addition to the vertex corrections, which are linear in κ λ , the trilinear coupling also generates corrections quadratic in κ λ through the wave function renormalization induced by the Higgs self-energy diagram. Such contributions have been computed for electroweak [21][22][23] and single-Higgs observables [4,[24][25][26][27][28]. Figure 1: One-loop diagrams involving the trilinear Higgs coupling contributing to the main single Higgs production processes: e + e − → hZ (top row) and e + e − → ννh (middle row). The Higgs self-energy diagram (bottom) gives a universal modification to all Higgs production processes via wave function renormalization.
Following Ref. [26], we can parametrize the NLO corrections to an observable Σ in a process involving a single external Higgs field as where Σ LO denotes the LO value, C 1 is a process-dependent coefficient that encodes the interference between the NLO amplitudes involving κ λ and the LO ones, while Z H corresponds to the universal resummed wave-function renormalization and is explicitly given by The impact of a deviation δκ λ ≡ κ λ − 1 from the SM value of the trilinear Higgs selfcoupling is therefore up to subleading corrections of higher orders in δZ H and C 1 . 4 The linear approximation in δκ λ is usually accurate enough to describe the deviations in single Higgs processes inside the typical constraint range |δκ λ | 5. We will nevertheless use the unexpanded δΣ expressions throughout this paper to derive numerical results. √ s for the e + e − → hZ and e + e − → ννh single Higgs production processes. Right: The linear dependence of production and decay rates on the δκ λ , δc Z , c ZZ and c Z parameters (see Section 2.2 for details on the meaning of these parameters). For e + e − → ννh, only the W W -fusion contribution is included. The dependence on δκ λ is amplified by a factor of 500.
The value of C 1 in Higgsstrahlung (e + e − → hZ) and W W -fusion (e + e − → ννh) processes are shown in the left panel of Fig. 2 as functions of the center-of-mass energy √ s. Very different energy dependences are observed for the two processes. A quick decrease is seen in Higgsstrahlung, from C 1 0.022 at threshold to about C 1 0.001 at a center-of-mass energy of 500 GeV. On the other hand, a nearly constant value C 1 0.006 is observed for the W W -fusion process over the same range of energy. Further numerical values are provided in Appendix A for both production and decay processes. Beside the inclusive production and decay rates, we also checked the impact of a correction to δκ λ on the angular asymmetries that can be exploited in e + e − → hZ → h + − measurements (see Refs. [29,30]). We found that these effects are almost negligible and have no impact on the fits.
To conclude this section, we show in the right panel of Fig. 2 the linear dependences of a set of production rates and Higgs partial widths on δκ λ and on three EFT parameters that encode deviations in the Z-boson couplings, δc Z , c ZZ and c Z (see Section 2.2 for a detailed discussion of the full set of BSM effects we are considering). Only leadingorder dependences are accounted for, at one loop for δκ λ and at tree level for the other parameters. One can see that the various observables have very different dependences on the EFT parameters. For instance, δc Z affects all the production processes in an energy-independent way. 5 On the contrary, the effects of c ZZ and c Z grow in magnitude for higher center-of-mass energy in both Higgsstrahlung and W W -fusion cross sections. It is apparent that the combination of several measurements can allow us to efficiently disentangle the various BSM effects and obtain robust constraints on δκ λ . From the sensitivities shown in Fig. 2, we can roughly estimate that a set of percent-level measurements in single-Higgs processes has the potential of constraining δκ λ with a precision better than O(1) and the other Higgs EFT parameters to the percent level. We will present a detailed quantitative assessment of the achievable precisions in the following.

Analysis of Higgs data at lepton colliders alone
Having obtained the one-loop contributions of δκ λ to single Higgs observables, we are now ready to determine the precision reach on the Higgs trilinear self-interaction. In order to obtain a robust estimate, we perform here a global fit, taking into account not only deviations in the Higgs self-coupling, but also corrections to the other SM interactions that can affect single-Higgs production processes.
For our analysis, we follow Ref. [5], in which the impact of single-Higgs measurements at lepton colliders on the determination of Higgs and electroweak parameters was investigated. We include in the fit the following processes • Higgsstrahlung production: e + e − → hZ (rates and distributions), • Higgs production through W W -fusion: e + e − → ννh, • weak boson pair production: e + e − → W W (rates and distributions), with Higgs decaying into a gauge boson pair ZZ * , W W * , γγ, Zγ, gg or pairs of fermions bb, cc, τ + τ − , µ + µ − .
New physics effects are parametrized through dimension-six operators within an EFT framework. For definiteness, we express them in the Higgs basis and refer to Ref. [31] for a detailed discussion of the formalism. Since CP-violating effects are strongly constrained experimentally, we exclusively focus on CP-conserving operators. We also ignore dipole operators and work under the assumption of flavor universality. We relax this assumption only to consider independent deviations in the of top, bottom, charm, tau, and muon Yukawa couplings.
To estimate the precision in the measurement of the EFT parameters, we assume that the central value of the experimental results coincides with the SM predictions and we neglect theory uncertainties. For simplicity we compute the SM cross sections at LO, neglecting NLO effects coming from SM interactions. These contributions can be important for the experimental analysis, since the modifications they induce in the SM cross sections can be non negligible compared to the experimental accuracy. For the purpose of estimating the bounds on BSM effects, however, they play a negligible role. We adopt a further simplification regarding electroweak precision observables, treating them as perfectly well measured. Such an assumption can significantly reduce the number of parameters to consider and is straightforward to implement in the Higgs basis which transparently separates the Higgs and electroweak parameters. The potential impact of this assumption will be discussed at the end of Section 4.
Under the above assumptions, we are left with twelve independent dimension-six effective operators that can induce leading-order contributions to single-Higgs and diboson processes. To this set of operators, we add the correction to the Higgs self-coupling parametrized by δκ λ . 6 The full list of parameters included in our fit contains: -corrections to the Higgs couplings to the gauge bosons: δc Z , c ZZ , c Z , c γγ , c Zγ , c gg , -corrections to the Yukawa's: δy t , δy c , δy b , δy τ , δy µ , -corrections to trilinear gauge couplings only: λ Z , -correction to the trilinear Higgs self-coupling: δκ λ .
Since our focus is on the future sensitivity on the trilinear Higgs self-coupling, we present results in terms of δκ λ only, profiling over all other parameters. For a detailed analysis of the sensitivity on the other operators see Appendix B and Refs. [5,14].
In our fit, we only include terms linear in the coefficients of the EFT operators, neglecting higher-order corrections. This approximation can be shown to provide very accurate results for all the parameters entering in our analysis [5]. The only possible exception is δκ λ , which can be tested experimentally with much lower precision than the other parameters. Although we checked that a linear approximation is reliable also for δκ λ , we keep Eq. (2.4) unexpanded in our numerical analyses. For simplicity, cross terms involving δκ λ and other EFT coefficients are however neglected, since the strong constraints on the latter coefficients and the loop factor make these contributions irrelevant.
In order to estimate the precision on Higgs measurements at different luminosities, we use a naive scaling with an irreducible 0.1% systematic error. This systematic error has no impact for the benchmark scenarios we consider, but becomes non-negligible for the largeluminosity projections presented at the end of this section (see Fig. 5). Another important source of uncertainty in our fit comes from the precision on the determination of trilinear gauge couplings (TGCs). In our analysis, we consider a range of possibilities. In the most conservative case, we assume 1% systematic errors in each bin of the e + e − → W W angular distributions used to constrain anomalous TGCs (aTGCs) [5]. In the most optimistic case, we assume that aTGCs are constrained much better than all the other parameters, so that they do not affect our fit. This is equivalent to enforcing the following relations among the EFT parameters: λ Z = 0 . 6 In the notation of Ref. [31] the δκ λ parameter corresponds to δλ 3 /λ. We start our discussion of the fit results by considering the benchmark scenarios for circular colliders. The profiled ∆χ 2 fit as a function of δκ λ is shown in the left panel of Fig. 3. The 68% CL intervals are also reported in Table 1.
The numerical results show that a 240 GeV run alone has a very poor discriminating power on the Higgs trilinear coupling, so that only an O(few) determination is possible (brown dashed lines in the plot). The constraint is also highly sensitive to the precision in the determination of TGCs, as can be inferred from the significantly different bounds in the conservative and optimistic aTGCs scenarios. The inclusion of measurements at 350 GeV drastically improve the results. An integrated luminosity of 200 fb −1 at 350 GeV, is already sufficient to reduce the uncertainty to the level |δκ λ | 1, whereas 1.5 ab −1 leads to a precision |δκ λ | 0.5.
It is interesting to compare the above results with the constraints coming from an exclusive fit in which only corrections to the trilinear Higgs coupling are considered and all the other parameters are set to zero. With 5 ab −1 collected at 240/250 GeV, and irrespectively of the presence of a run at 350 GeV, we find that such a fit gives a precision of approximately 14% in the determination of δκ λ . The strongest constraints come from the measurement of the e + e − → Zh cross section at the 240 GeV run, which is the observable with the largest sensitivity to δκ λ (see discussion in Section 2.2 and left panel of Fig. 2). Other processes at the 240 GeV run and the higher-energy runs have only a marginal impact on the exclusive fit.
The exclusive fit provides a bound much stronger than the global analyses, signaling the presence of a nearly flat direction in the global fits. We found that δκ λ has a strong lepton collider alone lepton collider + HL-LHC  Table 1: One-sigma bounds on δκ λ from single-Higgs measurements at circular lepton colliders (denoted as CC) and the ILC. The first column shows the results for lepton colliders alone, while the second shows the combination with differential measurements of both single and double Higgs processes at the HL-LHC. For each scenario two benchmarks with conservative and optimistic assumptions on the precision on trilinear gauge couplings are listed. The integrated luminosity is assumed equally shared between P (e − , e + ) = (±0.8, ∓0.3) for the ILC.
correlation with δc Z and c gg , while milder correlations are present with c Z and λ Z . 7 This result sheds some light on the origin of the improvement in the global fit coming from the combination of the 240 GeV and 350 GeV runs. The latter runs, although probing processes with a smaller direct sensitivity to δκ λ , are useful to reduce the uncertainty on the other EFT parameters. In particular, the 350 GeV run with 1.5 ab −1 of integrated luminosity allows for a reduction of the uncertainty on δc Z , c gg , c Z and λ Z by a factor of about 4. This in turn helps in lifting the flat direction in the global fit. This effect is clearly visible from the left panel of Fig. 4, which shows the fit on the δκ λ and δc Z parameters obtained with a 240 GeV run only and with the inclusion of a 350 GeV run.

Synergy between measurements at the HL-LHC and lepton-colliders
So far, we only considered the precision reach of lepton colliders on the extraction of the trilinear Higgs self-coupling. Significant information on δκ λ can however also be obtained at the high-luminosity LHC. It is thus interesting to estimate the impact of combining the different sets of measurements. The Higgs trilinear self-coupling can be accessed at the HL-LHC mainly through the exploitation of the Higgs pair production channel pp → hh. An analysis of this channel within the EFT framework has been presented in Ref. [32], in which the most promising channel, namely pp → hh → bbγγ, has been investigated. A fully differential analysis (taking into account the Higgs pair invariant mass distribution) allows to constrain δκ λ to the interval [−1.0, 1.8] at 68% CL. A second minimum is however present in the fit, which allows for sizable positive deviations in δκ λ , namely an additional interval δκ λ ∈ [3.5, 5.1] can not be excluded at 68% CL. Some improvement can be obtained complementing the pair-production channel with information from single Higgs channels, which are affected at NLO by the Higgs self-coupling. In this way, the overall precision becomes δκ λ ∈ [−0.9, 1.2] at 68% CL (with the additional minimum at δκ λ ∼ 5 excluded) and δκ λ ∈ [−1.7, 6.1] at 95% CL [33]. To estimate the impact of HL-LHC, we will use here the results of the combined fit with differential single and pair production (corresponding to the orange solid curve in the right panel of Fig. 3).
The combinations of the HL-LHC fit with our benchmarks for circular lepton colliders are shown in the right panel of Fig. 3. One can see that a 240 GeV run is already sufficient to completely lift the second minimum at δκ λ ∼ 5, thus significantly reducing the 2σ bounds. The precision near the SM point (δκ λ = 0) is however dominated by the HL-LHC measurements, so that the lepton collider data can only marginally improve the 1σ bounds. The situation is reversed for the benchmarks including a 350 GeV run. In this case, the precision achievable at lepton colliders is significantly better than the HL-LHC one. The combination of the LHC and lepton collider data can still allow for a significant improvement in the constraints if limited integrated luminosity can be accumulated in the 350 GeV runs (see Table 1). With 1.5 ab −1 collected at 350 GeV, on the other hand, the lepton collider measurements completely dominate the bounds.
Similar results are obtained for the low-energy ILC benchmarks. In this case, the lower integrated luminosity forecast at 250 GeV (2 ab −1 ) can be compensated through the exploitation of the two different beam polarizations P (e − , e + ) = (±0.8, ∓0.3). The only difference with respect to the circular collider case is the fact that the 250 GeV run fit is more stable under changes in the trilinear gauge couplings precision. This is due to the availability of runs with different polarizations, which provide better constraints on the EFT parameters. Analogously to the circular collider scenarios, the combination of the 250 GeV measurements with the HL-LHC data allows to completely lift the minimum at δκ λ ∼ 5, while a 350 GeV run would easily surpass the LHC precision. We report the results for the ILC benchmarks in Appendix B (see Fig. 14). For completeness, we mention that an exclusive fit on δκ λ at the ILC allows for a precision of approximately 32%, significantly better than the one expected through a global fit. Also in this case a nearly flat direction is present when deviations in all the EFT parameters are simultaneously allowed (see right panel of Fig. 4).
Having observed the significant impact of the combination of measurements at 240/250 GeV and 350 GeV center-of-mass energies, to conclude the discussion, we now explore a continuous range of integrated luminosities accumulated at the various colliders. The onesigma limits as functions of the integrated luminosity are displayed in Fig. 5 for the circular colliders and the ILC. Conservative and optimistic precisions for TGC measurements are respectively assumed to obtain the solid and dashed curves. The combination of runs at these two different energies always brings drastic improvements. The fastest improvements in precision on the δκ λ determination is obtained along the L 350 GeV /L 240 GeV 0.7 and L 350 GeV /L 250 GeV 0.5 lines for circular colliders and the ILC, respectively.  Figure 6: Higgs pair production cross sections at lepton colliders as functions of the center-ofmass energy (based on Fig. 7 of Ref. [36]) and illustrative diagrams. The difference between the two ννhh curves is entirely due to double Higgsstrahlung followed by invisible Z decay.

High-energy lepton machines
Having explored the reach of low-energy lepton colliders in the previous section, we now enlarge our scope to include machines with center-of-mass energies above 350 GeV. They offer the opportunity of probing directly the trilinear Higgs self-coupling through Higgs pair production processes, double Higgsstrahlung e + e − → Zhh and W W -fusion e + e − → ννhh in particular. The precision reach in the determination of δκ λ at ILC and CLIC has already been studied by the experimental collaborations [34,35]. These studies performed an exclusive fit, allowing for new-physics effects only in the trilinear Higgs self-coupling.
In this section, we first review the experimental projections on the extraction of the Higgs self-coupling from double Higgs channels. In this context, we also point out how differential distributions, in particular in the W W -fusion channel, can allow for an enhanced sensitivity to δκ λ . Afterwards, we reconsider Higgs pair production measurements from a global EFT perspective, showing how the determination of δκ λ is modified by performing a simultaneous fit for all EFT parameters. We also evaluate how these results are modified by combining double-Higgs data with single-Higgs measurements from low-energy runs.

Higgs pair production
As already mentioned, Higgs pair production at high-energy lepton machines is accessible mainly through the double Higgsstrahlung e + e − → Zhh and W W -fusion e + e − → ννhh channels. The cross sections for these two production modes as functions of the center-ofmass energy of the collider are shown in Fig. 6. It is interesting to notice their completely different behavior, so that the relevance of the two channels drastically changes at different machines. At energies below approximately 1 TeV, double Higgsstrahlung is dominant whereas, at higher energy, the channel with the larger cross section is W W -fusion. To be more specific, the cross section of double Higgsstrahlung reaches a maximum at √ s 600 GeV before starting to slowly decrease as the s-channel Z boson gets more and more offshell. On the contrary, the e + e − → ννhh cross section initially grows steadily with the center-of-mass energy of the collider and adopts a logarithmic behavior above 10 TeV. Notice that the e + e − → ννhh channel receives non-negligible contributions that are not of W W -fusion type. The largest of them arises from double Higgsstrahlung followed by a Z → νν decay. These contributions can however be efficiently identified at sufficiently high center-of-mass energies since the kinematic of the process is significantly different from that of W W -fusion. Notice, moreover, that both double-Higgs production cross sections are significantly affected by the beam polarization (see Appendix B and Fig. 15).
The e + e − → Zhh process at the ILC with 500 GeV center-of-mass energy has been thoroughly studied in Ref. [34]. A total luminosity of 4 ab −1 , equally split into two beam polarization runs P (e − , e + ) = (±0.8, ∓0.3), allows for a precision of 21.1% on the cross section determination through the exploitation of the hh → bbbb final state. A further improvement can be obtained by also including the hh → bbW W * channel, in which case the precision reaches 16.8%.
The e + e − → ννhh process has also been studied at a 1 TeV center-of-mass energy. A significance of 2.7σ (corresponding to a precision of 37%) could be achieved in the hh → bbbb channel, assuming and integrated luminosity L = 2 ab −1 and P (e − , e + ) = (−0.8, +0.2) beam polarization [37].
Studies of the e + e − → ννhh process at CLIC (both at 1.4 TeV and 3 TeV centerof-mass energy) are available in Ref. [35]. Assuming unpolarized beams and 1.5 ab −1 , the precision on the 1.4 TeV cross section could reach 44%. With 1.5 ab −1 , the 3 TeV cross section could be measured with a 20% precision. Both bbbb and bbW W * channels are included in these analyses, though the sensitivity is mainly driven by the former, as shown in Table 28 in Ref. [35].  The dependence of the Higgs pair production cross sections on δκ λ is shown in Fig. 7 for a set of benchmark scenarios. The SM cross section for each benchmark is provided in the legend. 8 Shaded bands show the precisions on the determination of the SM rates discussed above. Note the experimental collaborations made no forecast for the precision on double Higgsstrahlung at 1 TeV and above.
It is interesting to notice that, around the SM point, the sensitivity of both Higgs pair production channels to δκ λ gets milder at higher center-of-mass energy. On the contrary, the sensitivity to the other EFT parameters tends to increase with energy. Another important feature is the significant impact of terms quadratic in δκ λ on the behavior of the cross section around the SM point, especially for the W W -fusion channel shown in the right panel of Fig. 7. For this reason, a linear approximation is in many cases not sufficient to extract reliable bounds. In Table 2, we list the 68% and 95% CL bounds obtained from the benchmarks ILC and CLIC runs retaining the full dependence of the cross section on δκ λ . From Fig. 7, one can see that the interference between diagrams with and without a trilinear Higgs vertex has opposite sign in double Higgsstrahlung and W W -fusion. These two processes are thus more sensitive to positive and negative values of δκ λ respectively. A combination of double Higgsstrahlung and W W -fusion measurements could hence be used to maximize the precision for both positive and negative values of δκ λ . Such a scenario could be achieved at the ILC through the combination of a 500 GeV and a 1 TeV run. The impact of such combination can be clearly seen from the plot in the left panel of Fig. 8.
Being quadratic functions of δκ λ , inclusive cross sections (for each process and collider energy) can match the SM ones not only for δκ λ = 0, but also for an additional value of the trilinear Higgs self-coupling, resulting in a second minimum in the ∆χ 2 . In W Wfusion, the SM cross section is also obtained for δκ λ 1.08, 1.16 and 1.30 at center-of-mass energies of 1, 1.4 and 3 TeV, respectively. Whereas, for double Higgsstrahlung at 500 GeV, the SM cross section is recovered at δκ λ −5.8. This latter solution poses no practical problem for ILC since it can be excluded by HL-LHC measurements. Alternatively, it can be constrained by Higgs pair production through W W -fusion at 1 TeV, as well as through the indirect sensitivity of single Higgs measurements.
For CLIC, the secondary solutions at δκ λ 1 are more problematic. They can be constrained neither by HL-LHC data, nor by single Higgs measurements which are mostly efficient close to the threshold of the single Higgsstrahlung production. A more promising possibility is to exploit double Higgsstrahlung rate measurements. At center-of-mass energies above 1 TeV, however, they only provide weak handles on δκ λ . The e + e − → Zhh cross section becomes relatively small, being only 0.08 fb at 1.4 TeV with unpolarized beams. Moreover, the sensitivity to the trilinear Higgs self-coupling decreases with energy, as shown in Fig. 7. Since the experimental collaborations did not provide an estimate for the CLIC precision achievable on the SM e + e − → Zhh rate, we estimate it by naively rescaling the ILC 500 GeV projections by the total cross section at CLIC. We find that adding this information to inclusive e + e − → ννhh rates measurements only excludes the second minimum to the 1σ level (dashed orange line in the right panel of Fig. 8).
In addition, we consider the possibility of performing a differential analysis of double Higgs production through W W -fusion, studying whether a fit of the Higgs pair invariant mass distribution M hh can be sufficient to further exclude the δκ λ 1 points. The M hh distribution shows a good sensitivity to the Higgs trilinear, which mainly affects the shape of the distribution close to the kinematic threshold. This can be observed in   We estimate the impact of a differential analysis of the ννhh channel by performing a simple fit of the M hh invariant mass distribution. We consider either two or four bins, whose ranges are listed in Table 3. For simplicity, we work at parton level and assume a universal signal over background ratio across all bins. The right panel of Fig. 8 summarizes the result of the fits. It shows that a differential analysis can be useful in enhancing the precision on δκ λ . In particular, it allows us to exclude the second fit solution δκ λ 1.3 at the 68% CL, and to reduce significantly the 95% CL bounds for positive deviations in the Higgs self-coupling. For instance, the 4-bin fit restricts δκ λ to the range [−0. 18

Global analysis
It is important to verify whether the results discussed in Section 3.1, obtained assuming new physics affects only the triple Higgs coupling, are robust in a global framework once all other EFT parameters are taken into consideration. We therefore perform a global analysis at ILC and CLIC including measurements of both double-Higgs (Higgsstrahlung and W W -fusion) and single-Higgs processes (ννh, Zh, tth and e + e − h) in addition to diboson production.
We adopt the following benchmark scenarios chosen by the experimental collaborations for Higgs measurement estimates: • ILC: we follow the scenario in Ref. [20], assuming ILC can collect 2 ab • CLIC: we follow Ref. [35] and assume 500 fb −1 at 350 GeV, 1.5 ab −1 at 1.4 TeV and 2 ab −1 at 3 TeV can be collected with unpolarized beams. It should be noted that a left-handed beam polarization could increase the ννhh cross section and somewhat improve the reach on δκ λ .
For the global fit, we follow the procedure and assumptions adopted for the single Higgs processes fit at low-energy colliders. We also include the one-loop dependence on δκ λ in single Higgs production and decay processes, as done in Section 2. Such effects are also included in the top-Higgs associated production e + e − → tth and in ZZ-fusion e + e − → e + e − h, although they have a negligible impact. On the other hand, only the tree-level Higgs self-coupling dependence is considered in Higgs pair production processes, since one-loop corrections are numerically insignificant. As already stressed, the quadratic dependence on δκ λ in Higgs pair production processes cannot be neglected. In this case, cross terms between δκ λ and other EFT parameters are also accounted for. The linear approximation is adopted elsewhere. The estimates for the precision of the SM Higgs pair production cross section are taken from Refs. [34,35,37] already discussed in the previous section.
The results of the global fit for the ILC and CLIC benchmark scenarios are shown in Fig. 10. The 68% and 95% CL intervals are also listed in Table 4. It is interesting to compare these results with the ones obtained through the exclusive fit on δκ λ discussed in Section 3.1 (see Fig. 8). The χ 2 curves for ILC (up to 500 GeV or 1 TeV) and CLIC (no binning, 2 bins and 4 bins in M hh ) show very mild differences in the global fit with respect to the exclusive one. This demonstrates that the additional EFT parameters are sufficiently well constrained by single Higgs measurements and therefore have a marginal impact on the global fit. We also analyzed the impact of combining ILC and CLIC measurements with HL-LHC ones. The precision achievable at the LHC is significantly   poorer than the one expected at high-energy lepton colliders, so that the latter dominate the overall fit and only a mild improvement is obtained by combination.
We saw that allowing for other EFT deformations beside δκ λ does not worsen the global fit significantly. This result, however, was by no means guaranteed. To stress this point, we display in Fig. 11 the profiled χ 2 obtained by artificially rescaling the precision in single Higgs measurements. The ILC (up to 500 GeV, left panel) and CLIC (no binning in M hh , right panel) benchmarks are used as examples. For each collider, we show the results of the exclusive δκ λ analysis of the Higgs pair production measurements (solid black curve) and of the global analysis (dashed blue/cyan). The additional dashed curves correspond to global fits in which the precision in single Higgs and diboson measurements is rescaled by factors ranging from 0.5 to 10. It can be seen that the global fit is sizably affected by such a rescaling, in particular the fit precision is significantly degraded if single Higgs measurements become worse. This result shows that a comprehensive global analysis of the single Higgs measurements is crucial for obtaining robust constraints on δκ λ . Notice moreover that an improved precision on single Higgs measurements could have a positive impact on the determination of the Higgs self coupling at the ILC.
The impact of the uncertainty on the EFT parameters measurements on the extraction of the Higgs self-coupling from Higgs pair production was also recently investigated in Ref. [14]. It focused mainly on Higgs pair production through double Higgsstrahlung at ILC 500 GeV and on single-Higgs production in lower-energy runs, taking into account the uncertainties on SM parameters and electroweak precision observables. Loop-level contributions to single-Higgs processes coming from a modified Higgs self-coupling were not included in the fit, and the linear approximation was used to obtain the numerical results. The final fit takes into account runs at 250 and 500 GeV, with 2 and 4 ab −1 respectively equally shared between P (e − , e + ) = (∓0.8, ±0.3) beam polarizations. The estimated precision on the measurement of δκ λ is 30%, which is in good agreement with the constraints we obtained in our ILC benchmark scenario.

fusion.
We performed a global analysis, simultaneously taking into account corrections to the Higgs self-coupling and deviations in EFT parameters affecting Higgs interactions with other SM particles. The results of the analysis are summarized in Fig. 12 for the various benchmark scenarios considered. For each scenario, three sets of bounds are shown. Thin lines with vertical ends show the precision expected from measurements at lepton colliders only. The superimposed thick bars combine them with HL-LHC measurements. Finally, the thin solid and dotted lines are obtained by combining single Higgs measurements only at lepton colliders (1h) with the HL-LHC bounds. As discussed in the main text, unpolarized beams are assumed for the CEPC, FCC-ee and CLIC. For the ILC runs up to 500 GeV, an equal share of the luminosity at the two P (e − , e + ) = (±0.8, ∓0.3) beam polarizations is assumed, whereas a single polarization P (e − , e + ) = (−0.8, +0.2) is adopted at 1 TeV.
We found that a global analysis is essential to derive robust bounds on δκ λ . This is the case, in particular, if only low-energy lepton machines, such as CEPC or FCC-ee, are available. In this scenario, the Higgs self-coupling can be determined with good accuracy, around 40% at the 68% CL, by exploiting single Higgs measurements in the ννh and Zh channels as well as diboson production. In order to achieve this accuracy, it is essential to combine runs at different center-of-mass energy, for instance at 240 GeV and at 350 GeV, both with luminosities in the few attobarns range. Measurements at a single energy, in fact, leave a nearly flat direction unresolved in the global fit and lead to a very poor determination of δκ λ . Runs at two different energies can instead significantly reduce the flat direction by constraining with better accuracy on the other EFT parameters.
The high-energy linear colliders making direct measurements of the triple Higgs selfcoupling through pair production still provide the best constraints. Double Higgsstrahlung and W W -fusion yield complementary information, being more sensitive to positive and negative deviations in the Higgs self-coupling respectively. It is interesting to notice that the dependence of these two processes on δκ λ is stronger at lower center-of-mass energy, as shown in Fig. 7, so that ILC runs at 500 GeV and 1 TeV energy maximize the overall precision allowing for a determination of the trilinear Higgs self-coupling with a 20% uncertainty approximately, at the 68% CL.
High-energy measurements alone, such as the ones available with the 1.4 and 3 TeV CLIC runs, can only rely on ννhh production and have limited sensitivity to positive deviations in δκ λ . In this case, a second minimum in the global fit is present for δκ λ ∼ 1. The additional minimum can be excluded by performing a differential analysis exploiting the Higgs pair invariant mass distribution, whose threshold behavior is strongly sensitive to deviations in the Higgs self-coupling. A differential analysis can provide an order-20% determination of δκ λ at 68% CL, however at 95% CL values δκ λ 1 would still be allowed.
It is interesting to compare the above results with the ones achievable at the HL-LHC and at possible future hadron colliders. The HL-LHC is expected to be sensitive only to   Figure 12: A summary of the bounds on δκ λ from global fits for various future collider scenarios. For the "1h only" scenario, only single Higgs measurements at lepton colliders are included.
deviations of O(1) in the Higgs self-coupling. As one can see from Fig. 12, this precision is comparable to (or better than) the one achievable at low-energy lepton colliders with low integrated luminosity at 350 GeV runs. This is the case for our circular collider benchmarks with 200 fb −1 integrated luminosity at 350 GeV, as well as for the low-energy runs of the ILC. In these scenarios the HL-LHC data will still play a major role in the determination of δκ λ , while lepton colliders always help constraining large positive δκ λ that the HL-LHC fails to exclude beyond the one-sigma level. On the other hand, with 1 ab −1 of luminosity collected at 350 GeV, the lepton collider data starts dominating the combination. The situation is instead different at high-energy hadron colliders which can benefit from a sizable cross section in double Higgs production through gluon fusion. A pp collider with 100 TeV center-of-mass energy is expected to determine δκ λ with a precision of order 5% [3], thus providing a better accuracy than lepton machines. Intermediateenergy hadron machines, such as a high-energy LHC at 27 − 33 TeV could instead provide a precision comparable to that of high-energy lepton colliders. A rough estimate of the δκ λ determination at a 33 TeV pp collider gives a ∼ 30% precision at 68% CL for an integrated luminosity of 10 ab −1 .
To conclude the discussion, let us come back to our assumption of perfectly well measured electroweak precision observables. It seems fully justified if low-energy runs at the Z-pole are performed. This could for instance be the case at the ILC, CEPC, and FCC-ee which could respectively produce 10 9 , 10 10 , and 10 12 Z bosons. A Z-pole run for these machines can provide significant improvements with respect to LEP measurements (2 · 10 7 Z bosons), making electroweak precision observables basically irrelevant for the extraction of the Higgs trilinear self-coupling.
Without a new Z-pole run, evaluating the impact of a limited accuracy on electroweak precision observables might be less straightforward. An analysis of such scenario for the ILC collider has been recently presented in Ref. [14]. This work explicitly includes present constraints on m Z , the A asymmetry at the Z-pole, Γ Z→ll , Γ Z , Γ W and forecasts for improved m W , m H , and Γ W measurements, assuming no new run at the Z-pole. In that scenario, it is argued that Higgs measurements can be used to improve the constraints on the electroweak parameters. The achievable precision is sufficient to ensure that electroweak precision observables do not significantly affect the determination of δκ λ .
The precision necessary to decouple electroweak and Higgs parameters determinations in other benchmark scenarios might deserve further exploration. We think that electroweak precision measurements will have a negligible impact on trilinear Higgs selfcoupling determination at high-energy machines where Higgs pair production is accessible. This conclusion is supported by the results of Section 3 showing that the determination of δκ λ is only mildly affected by the other EFT parameters, once a wide-enough set of single Higgs measurements is considered. The situation for low-energy colliders, in which the Higgs self-coupling can be accessed only indirectly through single Higgs processes, is instead less clear. As we saw in Section 2, the precision on δκ λ obtained through a global fit is significantly lower than the one estimated through an exclusive analysis. Consequently, the precision of the single-Higgs and triple-gauge coupling extractions has a relevant impact on the fit. In principle, electroweak precision parameters could affect the bounds on single Higgs couplings and thus indirectly degrade the δκ λ constraint. This aspect might be worth a more careful investigation, which is however beyond the scope of the present work.  In this appendix we collect the numerical values of the coefficients C 1 , defined in Eq. (2.2), which encode the corrections to single-Higgs processes due to a deformation of the Higgs trilinear coupling. In Table 5 we report the C 1 coefficients for the total crosssection of the main single-Higgs production modes, namely Higgsstrahlung, vector-boson fusion and associated production with top quarks. Several values of the center-of-mass energy √ s are reported in the table, corresponding to the benchmark runs of future lepton colliders considered in main text. The calculation has been performed with the help of the public tools FeynArts, FormCalc, LoopTools, and CUBA [41][42][43].  Figure 13: Value of C 1 as a function of the center-of-mass energy √ s for the e + e − → hZ, e + e − → ννh, e + e − → he + e − and e + e − → htt single Higgs production processes. Notice that the result for Higgs production in association with a top-quark pair has been rescaled by a factor of 0.1.
Notice that the values of C 1 for Higgsstrahlung, W W -boson fusion and ZZ-boson fusion are independent of the beam polarization if we restrict ourselves to diagrams up to one loop, as we did in our analysis. As for e + e − → tth, the Higgs self-coupling gives rise to tiny beam polarization effects. Given the small impact of the latter production mode in our analysis, we can safely neglect such effects. The dependence of the C 1 coefficients on the collider energy is also shown in Fig. 13.
Besides the inclusive rates, we also checked the impact of a modified Higgs trilinear coupling on the angular asymmetries that can be built for the e + e − → hZ → h + − case (see Refs. [29,30]). We found that these effects are almost negligible and have no impact on our analysis.
For completeness, we also report in Table 6 the C 1 coefficients for the Higgs partial widths [26].

B Additional results
In this appendix, we collect some additional numerical results and plots that were not included in the main text.
In Fig. 14, we show the profiled ∆χ 2 as a function of δκ λ for the low-energy ILC benchmark considered in Section 2, including 2 ab −1 of integrated luminosity at 250 GeV and either 200 fb −1 fb or 1.5 ab −1 at 350 GeV with luminosities equally split into P (e − , e + ) =   (±0.8, ∓0.3) beam polarizations. In the left panel, we show the global fit for the ILC alone, while in the right panel we combine these results with the differential single and double Higgs measurements at the high-luminosity LHC. The corresponding 68% CL intervals are listed in Table 1. In Table 7, Table 8 and Table 9, we consider three alternative benchmark scenarios for the low-energy ILC runs. The three scenarios differ from the one considered in the main text by different choices of beam polarizations and luminosity splitting among them. The total integrated luminosities are the same as in the main benchmark, namely 2 ab −1 at 250 GeV, 200 fb −1 fb or 1.5 ab −1 at 350 GeV. In Table 7, we consider P (e − , e + ) = (∓0.8, ±0.3) beam polarizations with luminosity split between them according to a 70%/30% ratio. In Table 8 and Table 9, we consider P (e − , e + ) = (∓0.8, 0) beam polarizations with luminosity split between them with a 50%/50% ratio and a 70%/30% ratio respectively.
If only ILC data are included in the fit, the precision achievable in the case of a  Table 9: One-sigma bounds on δκ λ from single-Higgs measurements at the low-energy ILC. In this table, we consider a benchmark scenario with integrated luminosity split into P (e − , e + ) = (∓0.8, 0) beam polarization with a 70%/30% ratio.
P (e − , e + ) = (∓0.8, ±0.3) polarization with a 70%/30% luminosity split is slightly better than the one of the other scenarios. The impact is however marginal and basically disappears once the ILC data is combined with the high-luminosity LHC one. We find that the differences in the fits are mainly due to the dependence of the pair production cross sections on the beam polarizations. In Fig. 15, we show this dependence for the double Higgsstrahlung and W W -fusion pair production cross sections. These results are obtained with MadGraph5 [38] and do not take into account beam-structure effects. One can see that the largest cross sections are obtained for a P (e − , e + ) = (−0.8, +0.3) beam polarization. The cross sections for P (e − , e + ) = (0, 0) are smaller by a factor ∼ 2, while a much larger suppression is present for P (e − , e + ) = (+0.8, −0.3). 9 As a last result, we show the impact of the inclusion of the δκ λ parameter in the global fit on the EFT operators. For definiteness, we focus on the circular lepton colliders benchmarks. For the fit, we use the 12 EFT parameters considered in the main text, namely δc Z , c ZZ , c Z , c γγ , c Zγ , c gg , δy t , δy c , δy b , δy τ , δy µ , λ Z . As done in Ref. [5], it is convenient to slightly redefine the EFT parameters connected to the Higgs decays into γγ, Zγ and gg. In particular we define First of all, we focus on the fit obtained from low-energy lepton colliders only. In this case, the top Yukawa coupling and the Higgs contact interaction with gluons can not be accessed independently and can only be tested through the Higgs decay into gg. The δy t and c gg parameters thus always appear in combination as shown in Eq. (B.3).
In the global fit we include only thec eff gg parameter and not c gg and δy t separately. The precision on the various EFT parameters with and without the inclusion of δκ λ is shown in the upper panel of Fig. 16. One can see that, if only a 240 GeV run is available, the inclusion of the Higgs self-coupling in the fit significantly degrades the precision on δc Z andc eff gg . In this case, as we already discussed in the text, the precision on δκ λ is very low. The situation changes drastically in the presence of runs at 350 GeV. In this case, the precision onc eff gg is effectively decoupled from the determination of the Higgs trilinear coupling. Some correlation of δκ λ with δc Z is still present with 200 fb −1 of integrated luminosity at 350 GeV, while a much milder effect remains with 1.5 ab −1 of integrated luminosity. In the lower panel of Fig. 16, we show the global fit obtained after combination with high-luminosity LHC measurements. In this case, the top Yukawa and the Higgs contact interaction with gluons can be independently tested. The results of the global fit show that the inclusion of the Higgs trilinear coupling affects only the determination of δc Z . The impact is however much smaller than in the fit with lepton collider data only. Other EFT parameters are affected in a negligible way.