Constraining the Higgs self couplings at $e^+e^-$ colliders

We study the sensitivity to the shape of the Higgs potential of single, double, and triple Higgs production at future $e^+e^-$ colliders. Physics beyond the Standard Model is parameterised through the inclusion of higher-dimensional operators $(\Phi^\dagger \Phi- v^2/2)^n/\Lambda^{(2n-4)}$ with $n=3,4$, which allows a consistent treatment of independent deviations of the cubic and quartic self couplings beyond the tree level. We calculate the effects induced by a modified potential up to one loop in single and double Higgs production and at the tree level in triple Higgs production, for both $Z$ boson associated and $W$ boson fusion production mechanisms. We consider two different scenarios. First, the dimension six operator provides the dominant contribution (as expected, for instance, in a linear effective-field-theory(EFT)); we find in this case that the corresponding Wilson coefficient can be determined at $\mathcal{O}(10\%)$ accuracy by just combining accurate measurements of single Higgs cross sections at $\sqrt{\hat s}=$240-250 GeV and double Higgs production in $W$ boson fusion at higher energies. Second, both operators of dimension six and eight can give effects of similar order, i.e., independent quartic self coupling deviations are present. Constraints on Wilson coefficients can be best tested by combining measurements from single, double and triple Higgs production. Given that the sensitivity of single Higgs production to the dimension eight operator is presently unknown, we consider double and triple Higgs production and show that combining their information colliders at higher energies will provide first coarse constraints on the corresponding Wilson coefficient.


Introduction
In the Standard Model (SM), the breaking of the electroweak symmetry, SU (2) L ×U (1) Y → U (1) QED , is induced by the potential: where Φ is the Higgs doublet and the parameters µ and λ depend on the vacuum expectation value of the Higgs field v (or equivalently, the Fermi constant G F ) and the Higgs boson mass m H , i.e., µ 2 = m 2 H /2 and λ = m 2 H /(2v 2 ). The form of eq. (1.1) is dictated by the symmetries of the SM and the requirement of renormalisability. It is therefore a firm prediction of the SM that, once m H is known, the Higgs boson (H) self interactions are uniquely determined; λ SM 3 = λ SM 4 = λ, where λ SM 3 (λ SM 4 ) is the factor in front of the vH 3 (H 4 /4) interaction in the SM Lagrangian after ElectroWeak Symmetry Breaking (EWSB). Since its discovery in 2012 [1,2] a wealth of information has been accumulated on the scalar particle at 125 GeV of mass. Its couplings to vector bosons and third generation fermions [3] are so far all compatible with the SM expectations. However, no confirmation on the SM form of the Higgs potential is yet available from collider experiments. The reason of this lies in the intrinsic difficulty of accessing the relevant information experimentally.
Determining the form of the Higgs potential necessarily implies measuring the strength of the Higgs three-and four-(and possibly higher) point self-couplings. This is a challenging task at colliders for several reasons. As mentioned above, the self-couplings are proportional to λ, which in the SM is about 1/8 and therefore rather weak, i.e. of the same order of the Higgs couplings to the vector bosons and significantly smaller than the Yukawa coupling to the top quark. In addition, for a direct sensitivity on λ 3 (λ 4 ), processes featuring at least three(four) Higgs bosons need to be considered, namely double(triple) Higgs production As a result, effects associated to the self-couplings in the range of the SM values are in general very small. This simple fact has two immediate implications. First, one will need considerable statistics to be collected at the LHC Run II and III and in future electronpositron colliders before reaching the sensitivity to values close the SM predictions. Second, the precision of the experimental determinations of the self-couplings will critically depend on that of the other Higgs boson couplings entering the same process (or more in general the observable) under consideration. For example, at the LHC the largest production rate is due to gluon fusion processes via a top-quark loop. While the leading contribution from the top-quark Yukawa coupling y t scales as y 4 t , the leading contribution from the Higgs selfcoupling scales as (λ 3 ) 2 and is kinematically suppressed at large m(HH) invariant-mass values.
Many studies have been performed for the LHC at √ŝ = 13 TeV aiming to directly access λ 3 from double Higgs production measurements [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. However, due to the complexity of the corresponding realistic experimental setups, it is still unclear what is the final precision that could be achieved on the determination of λ 3 . At the moment the strongest experimental bounds on λ 3 (from non-resonant double-Higgs production) have been obtained in the CMS analysis for the bbγγ signature [22]. Exclusion limits on λ 3 have been found to strongly depend on the value of the top Yukawa coupling y t . In particular, in the case of an SM y t value, order λ 3 < −9 λ SM 3 and λ 3 > 15 λ SM 3 limits are obtained. According to optimist experimental projection studies [23], even with the high luminosity (HL) option of 3000 fb −1 it may be possible to exclude values only in the range λ 3 < −1.3 λ SM 3 and λ 3 > 8.7 λ SM 3 . Concerning the quadrilinear coupling λ 4 , it is instead incontrovertibly clear that the possibility of constraining λ 4 via the measurement of triple Higgs production at the LHC is quite bleak [24][25][26][27]. Even at a future 100 TeV proton-proton collider only loose bounds may be obtained with a considerable amount of integrated luminosity [28][29][30].
Additional and complementary strategies for the determination of λ 3 and λ 4 are therefore desirable not only at the moment but also for the (near) future. To this purpose, a lot of new results have recently appeared aiming to access λ 3 via indirect loop-induced effects. This idea has been pioneered by McCullough in the context of e + e − colliders [31], where loop-induced effects in single-Higgs production have been investigated for ZH associated production [32][33][34]. A first evaluation of analogous loop effects at the LHC has been pre-sented in ref. [35] for gg → H → γγ. At the same time, the complete set of (one-and two-) loop computations for all relevant single-Higgs observables at the LHC together with the proposal of combining inclusive and differential observables, has been put forward in [36]. Since then several studies have appeared: the computation of the factorisable QCD corrections to the single-Higgs EW production at the LHC [37], two-loop effects in precision EW observables [38,39] and, more recently, further investigations on the impact of the differential information and the relevance of SM electroweak corrections [40]. Furthermore, global analyses in the context of an SMEFT (SM-EFT) have also been presented for present and future measurements at the LHC [41] and even for the case of future e + e − colliders [42,43]. On the other hand, in these works, effects of λ 4 have been either ignored, being irrelevant for the calculation considered, or assumed to be determined in turn by the λ 3 value.
In the present work we investigate for the first time the (combined) sensitivity to both the λ 3 and λ 4 self-couplings in (multi-)Higgs production at future e + e − colliders. We consider H, HH, and HHH production both in association with a Z boson or via W -boson fusion (WBF) [44]. These processes are listed in Tab. 1, where we have also specified at Process λ 3 λ 4 ZH, ν eνe H (WBF) one-loop two-loop ZHH, ν eνe HH (WBF) tree one-loop ZHHH, ν eνe HHH (WBF) tree tree Table 1: Processes considered in this work and the order at which the λ 3 and λ 4 dependence appears. We do not calculate two-loop effects, but we do calculate one-loop effects for both single and double Higgs production.
which level in perturbation theory the λ 3 and λ 4 dependence appears (we do not calculate two-loop effects in this work). In particular, we perform the computation of one-loop effects in single and (for the first time) double Higgs production. The former pose no theoretical challenge, confirm the results of [31,43] (and mutatis mutandis, of [36,37]); they are presented here for completeness and are also used in our analysis. On the other hand, one-loop effects in double Higgs production can be computed only within a complete and consistent EFT approach, where UV renormalisation can be performed. To this purpose, we work in a theoretical and computational framework where the cubic and quartic couplings can independently deviate from the SM predictions and loop computations can be consistently performed. Specifically, we add the two higher-dimensional operators c 2n (Φ † Φ − v 2 /2) n /Λ (2n−4) with n = 3, 4 to the SM Lagrangian, where the presence of the "−v 2 " term considerably simplifies the technical steps of the one-loop calculation in double Higgs production. On the other hand, Wilson coefficients in this basis or in the standard c 2n (Φ † Φ) n /Λ (2n−4) parameterisation can be easily related at any perturbative order and also after the running to a different scale. While the c 2n coefficient are more suitable for the matching to a UV-complete model, the c 2n ones feature simple relations to Higgs self-couplings, and are more convenient for phenomenological predictions such as those performed here for (multi-)Higgs production. Independently from the choice of the basis, it will be clear in the text that the Wilson coefficients c 6 and c 8 , or c 6 and c 8 , are the relevant parameter to be considered and not directly the λ 3 and λ 4 couplings. By comparing and combining the direct and indirect sensitivities on the Higgs selfcouplings that could be obtained at future e + e − colliders (CEPC [45], FCC-ee [46], ILC [47] and CLIC [48,49]) running at different energies and luminosities, we explore the final reach of such colliders to constrain the Higgs potential. In general, we assume that BSM effects due to the Higgs interactions with the other SM particles are negligible w.r.t. those induced by self interactions. In practice, we work under the same assumptions of the first calculations of one-loop λ 3 effects in single Higgs production at e + e − [31] or proton-proton collisions [35,36], which represented an unavoidable input for the analyses considering a more general class of BSM scenarios [41][42][43]. On the other hand, precisely one of these recent global analyses, ref. [43], has shown that in high-energy e + e − collisions, where ZHH and WBF HH production are kinematically available, working with our assumption or allowing for additional BSM effects does not affect the constraints that can be obtained for the trilinear Higgs self-coupling, thus justifying our working assumption. In this work we will investigate the precision that can be achieved on the trilinear Higgs self-coupling, not only when it is close to its SM value, i.e., we will explore also regimes where its effects entering via loop corrections can be relevant also in double-Higgs production. Moreover, these loop corrections, similarly to ZHHH and WBF HHH production at the Born level, involve effects due to the quadrilinear Higgs self-coupling. We will investigate the constraints that can be set on this coupling under two different assumptions. In the first, we consider the case of a well behaving EFT expansion, where dimension-eight operators induce effects smaller than dimension-six ones. In other words, the value of the trilinear coupling automatically sets also the value of the quadrilinear coupling. In the second, we lift this assumptions and we allow for independent trilinear and quadrilinear couplings, namely, we allow for similar effects from (Φ † Φ) 3 and (Φ † Φ) 4 operators.
The paper is organised a follows. In section 2 we introduce the notation and we discuss the EFT framework used in our calculation. The details concerning the definition of the renormalisation scheme at one loop and all the necessary counterterms for the calculation performed here are given in Appendix A. In section 3 we provide the predictions for cross sections of single, double and triple Higgs production at different energies, discussing their dependence on the c 6 and c 8 parameters. The one-loop calculations in single and double Higgs production are performed via one-loop form factors, the explicit results for which are provided in Appendix B. In section 4 we determine the reach of several experimental setups at future e + e − colliders for constraining the cubic and quartic couplings; both individual and combined results from single, double and triple Higgs production are scrutinised. The maximum c 6 and c 8 values beyond which perturbative convergence cannot be trusted are derived in Appendix C. We draw our conclusions in section 5.

Notation and parametrisation of New Physics effects
In this work we are interested to the effect induced by the modification where the New Physics (NP) modifications of the potential are all included in V NP and the symbol Φ denotes the Higgs doublet. The term V SM has already been defined in eq. (1.1). Following the convention of ref. [50], the most general form of V NP that is invariant under SU (2) symmetry can be written as It is important to specify from the beginning why for our calculation it is convenient to parametrise the NP contributions as done in eq. (2.2) and not using the standard EFT parameterisation 3) The advantages of the parametrisation in eq. (2.2) w.r.t the one in eq. (2.3) are due to the fact that after EWSB any Φ † Φ n originates H i terms with 1 ≤ i ≤ 2n, while any In other words, at tree-level, the trilinear Higgs self-coupling receives modifications only from c 6 and the quadrilinear only from c 6 and c 8 . Needless to say, when they are summed to V SM , equations (2.2) and (2.3) not only refer to the same quantity parametrised in a different way (V SM + V NP std = V SM + V NP ), but they are also fully equivalent for any truncation of the series at a given order n.
Writing V SM (Φ) + V NP (Φ) after EWSB as allows to define the self-couplings λ n , which can be parametrised by the quantities 1 where λ SM 3 (λ SM 4 ) is the value of λ 3 (λ 4 ) in the SM and reads In other words,c 6 ,c 8 andc 10 are c 6 , c 8 and c 10 normalised in such a way that can be easily related to κ 3 , κ 4 and κ 5 . In particular Eqs. (2.5) and (2.6) make two important points manifest. First, while the trilinear coupling only depends on c 6 , the quadrilinear depends on both c 6 and c 8 . Thus, in a well-behaved EFT, where the effects of higher dimensional operators are systematically suppressed by a large scale, one expects deviations both in κ 3 and κ 4 and such that (κ 4 − 1) 6(κ 3 − 1), see also eq. (2.10). Second, κ 3 and κ 4 do not depend on any c 2n coefficient with n > 4, without any assumption on the c 2n values with n > 4. In other words, for the study of trilinear and quadrilinear Higgs self-couplings at the tree level, only the c 6 and c 8 Wilson coefficients are relevant. On the other hand, at one-loop level also c 10 is in principle needed. As already mentioned in the introduction, in this paper we will calculate one-loop corrections to double Higgs production, taking into account both c 6 and c 8 effects. In general, when loop corrections are calculated and c 2n coefficients themselves are renormalised, the parametrisation of eq. (2.2) is convenient also for other reasons that are particularly relevant when the Wilson coefficients are renormalised. At variance with eq. (2.3), the values of the coefficients c 2n influence only the value of the Higgs self-couplings; they do not alter the SM relations among m H , v, µ and λ This is convenient because the physical quantities are m H and v and not µ and λ, while using eq. (2.3) one has to determine before µ and λ and then the self-couplings (see Appendix A of [36]). Especially, and this is the main motivation for this work, thanks to eqs. (2.11) and (2.12) the modification V SM (Φ) → V SM (Φ) + V NP (Φ) allows to keep the SM relations between the renormalisation constants and the definition of the renormalisation counterterms as done in ref. [51]. On the other hand, the explicit insertion of v 2 in eq. (2.2) deserves a particular attention for the renormalisation procedure and leads to additional counterterms. All the necessary ingredients for the calculations performed here are provided in Appendix A.
It is important to note that the coefficients c 2n in eq. (2.2) and c 2n in eq. (2.3) are connected by very simple relations and can be converted ones into the others at any step of the calculation. Thus, our calculation is fully equivalent to using the c 2n coefficients and renormalising them in the MS scheme (see Appendix A). Via the simple tree-level relations among c 2n and c 2n coefficients one can convert results at any order from one convention to the other, including the renormalisation-group equations. As already said, while c 2n coefficients are more suitable for the matching to a UV-complete model, the c 2n coefficients are more convenient for one-loop calculations of hard-scattering matrix elements such as those performed here for double Higgs production. Independently from the basis choice, the NP effects should be parametrised via c 6 and c 8 , or c 6 and c 8 , rather than directly through the λ 3 and λ 4 couplings. The reason is that the individual Wilson coefficients are the quantities entering the renormalisation procedure, and c 6 , or equivalently c 6 , is affecting both λ 3 and λ 4 values. In the following we will mainly parametrise NP via thē c 2n coefficients, which are simply related to both c 2n and κ n (see eqs.(2.5)-(2.10)).
3 Single, double and triple Higgs production:c 6 andc 8 dependence In this section we describe the calculations for single, double and triple Higgs production via WBF, the e + e − → ν eνe H(H(H)) processes, or in association with a Z boson, the e + e − → ZH(H(H)) processes. We will denote the latter also as ZH n . Besides the ZH n and WBF production modes, at e + e − colliders there are other possibly relevant processes such as ttH n , Z-boson-fusion (ZBF) or loop-induced H n production via photon fusion [50] from the initial-state radiation. However, these processes have considerably smaller cross sections than WBF or ZH n production modes, so we do not consider them in our analysis. Part of them have been considered in ref. [43] and their impact has been indeed found to be negligible.
While triple Higgs production processes are calculated at the Born level, for both single and double Higgs production we take into account also one-loop corrections involving the additionalc 6 andc 8 dependence. The sensitivity onc 6 andc 8 , and in turn on κ 3 and κ 4 , depends on the multiplicity of Higgs bosons in the final state and the numbers of loop corrections considered, as summarised in Tab. 1. We expect complementary information from ZH n and WBF, when the different collider energies of the possible future e + e − colliders are considered. While the cross section of ZH n is maximal for energies slightly larger than its production threshold, the WBF cross section grows with the energy. Moreover, based on results in refs. [31,36,37,40,43], in ZH n production we expect a strong dependence of the Higgs self-coupling effects from loop corrections, with larger effects at lower energies. On the contrary, in WBF this energy-dependence is expected to be much smaller.
It important to note that ZH n and WBF cross sections, and in turn their sensitivity onc 6 andc 8 , depend on beam polarisations, which can be tuned at linear colliders. 2 First of all, WBF contributes only via the LR polarisations, 3 since W -boson couples only to the left-chirality fermions. Conversely, the ZH n processes can also originate from RL polarisations (right-handed e − , left-handed e + ), also denoted as P (e − , e + ) = (1.0, −1.0).
On the other hand, results for RL polarisations can be easily obtained from those with LR via the relation In all our calculations we use following input parameters [53] G µ = 1.166 378 7 × 10 −5 GeV −2 , m W = 80.385 GeV , m Z = 91.1876 GeV We assumec 6 andc 8 measured at the scale µ r = 2m H , which therefore we will also use as MS renormalisation scale forc 6 in the double Higgs computation. The WBF production modes feature the same final states of ZH n production with Z → ν eνe decays. The latter are not considered as part of the WBF contribution in our calculation. Although NLO EW corrections in the SM would jeopardise the gauge invariance of this classification at the amplitude level [54], this is not the case for the one-loop corrections induced by additional c 2n interactions, which is the kind of effects we consider in this work on top of LO effects, as better specified later in eqs. (3.5) and (3.11). Moreover, the interference of ZH n -type diagrams with "genuine" WBF configurations is negligible

Single Higgs production
In this section we briefly (re)-describe the calculation of loop-induced effects fromc 6 (κ 3 ) in ZH and single-Higgs WBF production at e + e − colliders (representative diagrams are shown in Fig. 1). We introduce the notation that will be generalised to the case of double Higgs and triple Higgs production and we show how it is related to the previous calculations [31,36,40,43].
For both WBF and ZH production channels no Higgs self-coupling contributes at the tree level. On the other hand, one-loop corrections depends on the trilinear Higgs selfcoupling, but not on the quartic. Thus, while the LO cross section σ LO (H) is only of SM origin, NLO predictions includes also effects fromc 6 : proportional toc 8 or any other c 2n coefficient in this expansion, meaning that eq. (3.5) is actually exact; no other terms can enter at all even for higher orders in the (v/Λ) expansion. Furthermore, we remind that, at variance with the case of double Higgs production, in single Higgs production at one-loop the anomalous coupling approach (κ 3 ) is fully equivalent to the calculation in the EFT (c 6 ). For our phenomenological study we ignore the SM NLO EW corrections [54,55]. Our main focus is not the precise determination ofc 6 , but the study of its impact via its leading contributions. As discussed in detail in ref. [40], SM NLO EW corrections have a tiny impact on the extraction of the value ofc 6 and do not affect the accuracy of the determination ofc 6 . Therefore, we considerc 6 effects at one loop via the following approximation With this approximation, the sensitivity to the trilinear coupling can be expressed via the ratio where we have expressed the σ i /σ LO ratios directly 5 using the symbols C 1 and C 2 introduced in ref. [36]. C 1 denotes the one-loop virtual contribution involving one triple Higgs vertex, while C 2 originates from the Higgs wave-function renormalisation constant (see eqs. (A.2),(A.14) and (A.17)), which is the only source ofc 2 6 and thus κ 2 3 dependence at one loop level. Both C 1 and C 2 are independently UV-finite and, for simplicity, we choose not to resum higher-orders contributions to the wave function, at variance with ref. [36]. Indeed, given the results already presented in ref. [31], we expect to bound κ 3 close to the SM (κ 3 = 1) and in this scenario such a resummation would not make a noticeable difference anyway. Moreover, even considering κ 3 in the range |κ 3 | < 6 from ref. [56], the difference between the formula in eq. (3.7) and including the resummed higher-order contributions to Z H is below 1% (see also ref. [40]). Considering C 2 in eq. (3.8), the difference w.r.t. the definition in ref. [36] is only due to this choice, however, in the limitc 6 → 0(κ 3 → 1) the two different definitions are equivalent as can be seen from the value of C 2 : Moreover, in the limitc 6 → 0, a linear expansion of eq. (3.7) for ZH would lead to the result in ref. [31]. As explained in ref. [36] for hadronic processes, C 1 parametrises contributions that are process and kinematic dependent. In Fig. 2, we show σ LO (left plot) and C 1 (right plot) for ZH (red) and WBF (green) production as function of the energy of the collider √ŝ . As expected, while C 1 strongly depends on √ŝ for ZH, it does very mildly for WBF H. In particular, for ZH, when increasing the energy, C 1 decreases at the beginning, then changes its sign around √ŝ = 550 GeV and remains small. On the other hand, the total cross section for ZH production peaks at around √ŝ = 240 GeV and decreases as √ŝ increases, while the cross section for WBF H production increases with √ŝ . Thus, while for the range 200 − 500 GeV the ZH production is expected to be more sensitive than WBF onc 6 (κ 3 ), at higher energies the situation is reversed. The information from collisions at different energies, or even at different colliders, increases the sensitivity on κ 3 , as it has been discussed in ref. [43]. We will show analogous results in sec. 4. We have also looked at the differential distribution for the transverse momentum of the Higgs boson, but we have not seen any strong dependence on C 1 . Hence, for single Higgs production at e + e − colliders differential distribution cannot increase the sensitivity on κ 3 , at variance with the case of hadron colliders [36,37,40,41] and of double-Higgs production [57].
The range of validity of this calculation in κ 3 and in turnc 6 is mainly dictated by the effects from δZ NP H , as discussed in ref. [36], from which the bound |κ 3 | < 20 can be also straightforwardly applied here. A more cautious and conservative condition can be derived by requiring perturbative unitarity for the HH → HH scattering amplitude and/or perturbativity for the loop corrections to the HHH vertex in any kinematic configuration. This bound has been derived in ref. [56] and leads to the requirement |κ 3 | < 6, independently from the value of κ 4 . However, the kinematic configuration leading to this bounds are those involving two Higgses on-shell and the virtuality of the third Higgs close to 2m H , which is not relevant for the trilinear interaction entering in single-Higgs production. We independently re-investigated this bound onc 6 (and analogous ones onc 8 ) in Appendix C, where its derivation is discussed in detail.

Double Higgs production
We now consider double Higgs production. The cross sections for the production of two Higgs bosons in association with a Z bosons (e + e − →ZHH) and via WBF (e + e − →ν eνe HH) do depend on the trilinear Higgs self-coupling already at the tree level (see diagrams in Fig. 3). Moreover, for both processes, one-loop corrections depend on both the trilinear and quartic Higgs self-couplings. At leading order ZHH and WBF HH cross sections can be written as 6 : (3.10) 6 The σi terms entering eq. (3.10) are not the same quantities appearing in eq. (3.5). where σ 0 is the SM result, σ 1 represents the leading contribution in the EFT expansion (order (v/Λ) 2 ), while σ 2 is the squared EFT term of order (v/Λ) 4 . Note that within our choice of operators there is no contribution proportional toc 8 in this expansion. Actually, no c 2n coefficient with n > 3 enters at the tree level. The NLO corrections involve several different contributions. First we classify all of them and then we specify those relevant for our study. Using a notation that is analogous to eq. (3.10), the cross section at NLO accuracy can be parametrised as (3.12) (3.13) (3.14) where the σ ij quantities refer to the one-loop terms that factorisec i 6c j 8 contributions and the σ i0j to those proportional toc i 6c j 10 . Some comments on the terms in (3.12), (3.13), (3.14) and (3.15) are in order.
The terms in (3.12) are the NLO EW corrections to the contributions that appear already at LO. The quantity σ 00 , for instance, corresponds to the NLO EW corrections in the SM and has been calculated in ref. [58]. The terms σ 10 and σ 20 are O(α) corrections to σ 1 and σ 2 , respectively, and thus they are always subdominant. They should be included for precise determination ofc 6 values, yet being subdominant, we neglect them together with σ 00 in this first analysis, similarly to the case of single Higgs production.
The terms in (3.13) collect contributions that appear at NLO for the first time. For small valuesc 6 1, these terms are suppressed w.r.t. σ 1 and σ 2 in (3.10), and may be neglected. However, at variance with σ 10 and σ 20 , for large values ofc 6 , they are not subdominant. Thus, we keep them in order to study thec 6 and in turn κ 3 dependence beyond the linear approximation, which as explained is not sufficient for large values of c 6 . Moreover, this allows also to better quantify the range of validity of our perturbative calculation (see Appendix C). These contributions originate from the left diagram of Fig. 4, which shows also the other possible one-loop corrections to the HHH vertex. This diagram induces a (c 6 ) 3 dependence in the amplitude and in turn a (c 6 ) 4 dependence in σ 1−loop (HH) via the interference with Born diagrams. As it has been discussed in ref. [56], its contribution can be large. Also, the presence of (c 6 ) 3 effects indicates that terms up to the order (v/Λ) 6 have to be taken into account in the one-loop amplitudes and thus in the renormalisation constants. Schematically, each order in the (v/Λ) expansion implies that the following terms can be in principle present Thus, the full dependence on λ 3 and λ 4 of the diagrams appearing in Fig. 4 is taken into account. 7 On the other hand, (v/Λ) 6 terms include c 10 contributions, which we reparametrised in term ofc 10 ≡ (c 10 v 6 )/(λΛ 6 ); they lead to an independent value also for λ 5 , the factor in front of the H 5 /v term appearing in V NP (Φ) after EWSB. The origin of the terms in (3.14) and (3.15) can be now understood on the base of Fig. 4 and eqs. (3.16)- (3.18) and are commented in the following. The terms in (3.14) are the contributions that depend onc 8 . Thus, σ 01 , σ 11 and σ 21 are the most relevant quantities in our one-loop study of double Higgs production, as they provide the sensitivity toc 8 and therefore to the deviation from the quadrilinear that one expects on top of the one determined byc 6 only. Although σ 11 and σ 21 would be suppressed for small c 6 we keep them to study the validity of our calculation in the (c 6 ,c 8 ) plane, or equivalently (κ 3 , κ 4 ) plane.
Finally, the last term (3.15), is related toc 10 -dependent contributions. These contributions arise from the diagram with the H 5 interactions in Fig. 4 and the corresponding term in the renormalisation constant δc 6 (see eq. (A.16) for the explicit δc 6 formula), and can be expressed as At one-loop in ZHH or WBF production their sum can be written as a kinematically independent shift toc 6 , In practice we can only constrain a linear combination ofc 6 andc 10 that is in eq. (3.20). In the following we work in the assumptions thatc 10 effects are negligible and we setc 10 = 0, however, for not too large values ofc 10 , i.e., where the linear expansion inc 10 is reliable, results ofc 6 can be translated into a linear combination ofc 6 andc 10 via eq. (3.20). 8 In order to be directly sensitive toc 10 one would need to consider one-loop effects in triple Higgs production, or evaluate quadruple Higgs production at the tree level.
In conclusion, in our phenomenological analysis, we evaluatec 6 andc 8 effects at one loop via the following approximation σ pheno NLO (HH) = σ LO (HH) + ∆σc 6 (HH) + ∆σc 8 (HH) , ∆σc 6 (HH) =c 3 6 σ 30 + σ 40c6 , The analytical results for the form factors used for the calculation of ∆σc 6 (HH) and ∆σc 8 (HH) are given in Appendix B. We show now the impact ofc 6 andc 8 in the σ pheno NLO predictions at different energies. First of all, in Fig. 5 we show the LO cross section σ LO of ZHH (left) and WBF (right) production as function of √ŝ for different values ofc 6 . In ZHH production, the LO cross section peaks around √ŝ = 500 GeV, which is the optimal energy for measuring this processes, while WBF HH cross section grows with energy. As can be seen by comparing the left and right plot, the dependence onc 6 is different in ZHH and WBF HH production. Especially, at variance with ZHH, WBF HH cross sections in general increase when c 6 = 0. This feature is even more clear in the top-left plot of Fig. 6, where we show the dependence of σ LO onc 6 for the different phenomenologically relevant configurations that will be analysed in sec. 4, namely, ZHH at 500 GeV collisions and WBF HH at 1, 1.4 and 3 TeV collisions.
Using a similar layout, in Fig. 6 we display three other plots, which show the dependence of σ LO , ∆σc 6 (HH) and ∆σc 8 (HH)/c 8 onc 6 for different processes and energies. Specifically, in the upper-right plot we show the case of ZHH at 500 GeV, while in the lower plots we show WBF HH at 1 TeV (left) and 3 TeV (right). In these three plots we display σ LO , which has also been shown in the top-left plot, as a black line and ∆σc 6 (HH) and ∆σc 8 (HH)/c 8 as a blue and red line, respectively. Thus the blue line directly shows thec 8 -independent part of σ pheno NLO , while the red one corresponds to the coefficient in front of thec 8 -dependent part ∆σc 8 (HH), which in turn depends onc 6 . For both cases, a shortdashed line is used when ∆σc 6 (HH) or ∆σc 8 (HH)/c 8 are negative. From Fig. 6 we can see that not only for the LO prediction (top-left plot) but also for one-loop effects thec 6 (as wellc 8 ) dependence is very different in ZHH (top-right plot) and WBF HH (lower plots) production. On the other hand, as can be seen in the lower plots, besides a global rescaling factor, WBF HH results are not strongly affected by the energy of e + e − collisions. 9 In the case of ZHH production at 500 GeV, the minimum of the LO cross section is at c 6 ∼ −3, while for WBF HH it is atc 6 ∼ 0.5. This minimum is given by cancellations induced by the interference of diagrams featuring or not the HHH vertex. Such pattern of cancellations is different in the ∆σc 6 (HH) one-loop contribution, which in absolute value is instead minimal atc 6 = 0 and very large at large values ofc 6 . For this reason, e.g., for c 6 < −3 the ∆σc 6 (HH) one-loop contribution is larger than the LO cross section. This does not signal the breaking of the perturbative convergence, rather it is due to the large cancellations that are present in this region only in the LO cross section; as already said, the perturbative limits, which are derived in Appendix C, require |c 6 | < 5 and correspond to the range of the plot. In the case of WBF HH production ∆σc 6 (HH) is always smaller than σ LO , being negative forc 6 > 0 and positive forc 6 < 0.
Regarding the ∆σc 8 (HH) contribution, which we display in the red lines normalised with 1/c 8 , the effect is very different in ZHH and WBF HH production. In the case of ZHH production ∆σc 8 (HH) is always negative and the minimum in absolute value is very close to the minimum of the LO prediction. In the case of WBF HH production ∆σc 8 (HH) change sign atc 6 ∼ −2 andc 6 ∼ 0.5, being positive between these two values and negative outside them. In general, in absolute value, the ratio ∆σc 8 (HH)/σ LO is always below c 8 · 2% value. Still, given the allowed perturbative range |c 8 | < 31 (see Appendix C), effects from large values ofc 8 can be in principle probed.

Triple Higgs production
In triple Higgs production cubic and quartic self-couplings are present already at the treelevel and therefore both the leading dependences onc 6 andc 8 are already present at LO (see diagrams in Fig. 7). Following the same notation used for double Higgs production, the cross section used for our phenomenological predictions can be written as where the σ 00 term corresponds to the LO SM prediction. Similarly to the case of double Higgs production at one loop, terms up to the eighth power in the (v/Λ) expansion are present at the cross section level, although in this case only the fourth power is present at the amplitude level. The upper bounds onc 6 andc 8 mentioned in the previous section and discussed in Appendix C have to be considered also in this case. It is important to note that although for large values ofc 6 andc 8 loop corrections may be sizeable, at variance with double Higgs production,c 6 andc 8 are both entering at LO. Thus, when limits onc 6 andc 8 are extracted, loop corrections may slightly affect them, but only for largec 6 andc 8 values. In Tab. 2 we give all the σ ij /σ 00 ratios, so that the size of all the relative effects from the different NP contributions can be easily inferred. 10 In Fig. 8, we show σ LO at different energies for representative values ofc 6 andc 8 , including the SM case (c 6 = 0,c 8 = 0) where σ LO = σ 00 . There, we also explicitly show the value of the σ 02 component, which factorises 10 There are large cancellations among the different contributions; more digits than those shown here have to be taken into account in order to obtain a reliable result. the (c 8 ) 2 dependence. We can see that for ZHHH production (left) the sensitivity toc 8 is rather weak. The σ 02 component is just around 1% of σ 00 , which means that even for large values ofc 8 the total cross section would not be large enough to be measurable at the future colliders considered in this study (see discussion in sec. 4). On the other hand, the total cross section of WBF HHH increases with the energy, as for single and double Higgs production. Especially, the σ 02 component is much larger; it is of the same order of the SM σ 00 component. As an example, assumingc 8 = 1(c 8 = −1) andc 6 = 0, σ LO at 3 TeV is 1.7 (2.2) times larger than σ 00 . For largec 8 values, σ LO ≈c 2 8 σ 02 ≈c 2 8 σ 00 . As can be seen in Tab. 2, WBF is also very sensitive onc 6 ; for large values ofc 6 indeed σ LO ≈c 4 6 σ 40 and in particularc 4 6 σ 40 ≈c 4 6 σ 00 at 3 TeV. All these effects are even larger at lower energies.

Bounds on the Higgs self-couplings
In this section we study how thec 6 andc 8 parameters can be constrained at future lepton colliders via the analysis of single, double, and triple Higgs production. We consider four future e + e − colliders, CEPC [45], FCC-ee [46], ILC [47], and CLIC [48,49], with different operations modes 11 that are summarised in Tab. 3. In the following, we will refer to the different scenarios as "collider-√ŝ " like, e.g., CLIC-3000. Although higher integrated luminosities can be attained at the CEPC and FCC-ee, energies as high as at the ILC and CLIC cannot be reached, since they are circular colliders. As a result, only single Higgs production can be measured at the CEPC and FCC-ee, and therefore only indirect constraints via loop corrections can be set onc 6 . Instead, at the ILC and CLIC double Higgs production can be measured. With this process, bothc 6 andc 8 can be constrained, the former via the direct dependence at the Born level and the latter via the indirect dependence 11 At the ILC also an operation mode at √ŝ ∼ 350 GeV is expected, but studies mainly focused on the scan of the tt production threshold, ignoring Higgs physics. At CLIC also a slightly different scenario at 380 GeV instead of 350 GeV may be possible.  through loop corrections. Moreover, even triple Higgs production is kinematically allowed at the ILC and CLIC, allowing to set direct constraints onc 8 . In our analysis we consider the following two scenarios 12 : 1. As expected from a well-behaving EFT expansion, the contribution fromc 8 is suppressed and we can safely setc 8 = 0. We explore how well we can measurec 6 , not only assumingc 6 ∼ 0, i.e., an SM-like configuration, but also allowing for large BSM effects viac 6 = 0.
2. The value ofc 8 can be different from zero and leads to non-negligible effects. We explore how well we can constrainc 8 and how muchc 8 can affect the measurement ofc 6 .
First, we study the sensitivity of ZH n and WBF processes at the various colliders considered. Then we show combined results for the ILC and CLIC. It is important to note, however, that single Higgs production depends onc 8 only via two-loop effects, which we did not calculate in this work (see Tab. 1). Thus, we cannot directly combine single Higgs with double Higgs and triple Higgs in the case of Scenario 2. Nevertheless, we discuss the limit that can be obtained in single Higgs production under the assumption that thē c 8 -dependent two-loop effects are negligible.

Single Higgs production
In this section we discuss the constraints that can be obtained onc 6 via the single-Higgs production modes. As said, since the effects ofc 8 are unknown, we restrict our study to the case where it can be ignored, i.e., Scenario 1. We start by considering the case in which we assume that the Higgs potential is like in the SM (c 6 = 0) and then we consider the BSM case withc 6 = 0. 12 One may be tempted to explore the regimec6 = 0 andc8 = 0, too. However, this condition is neither motivated by an EFT expansion nor protected by any symmetry. As can be seen from eq. (A. 16 Table 4: Expected precision for the measurements of single Higgs production modes and the expected 1σ and 2σ constraints onc 6 , assuming an SM measurement, are listed. The value of for the CEPC has been taken or obtained via a luminosity rescaling from ref. [45], for the FCC-ee from ref. [46], for the ILC from refs. [47,59] and for the CLIC from ref. [49]. In Tab. 4 we show 1σ and 2σ constraints onc 6 that can be obtained via ZH and WBF H at different energies and colliders, using eq. (3.7). We show also the value of C 1 and the accuracy that can be achieved in any experimental setup, as provided in [45-47, 49, 59] or obtained from them via a luminosity rescaling. 13 In general in this work, unless differently specified, we assume Gaussian distributions for the errors and no correlations among them, and the errors are rescaled according to cross section in BSM cases. In the results of Tab. 4 we did not take into account effects due toc 6 in the Higgs decay, since, at variance with the LHC case, they can be in principle neglected at e + e − colliders. Indeed, the total cross section of e + e − → ZH production can be measured via the recoiling mass method [47], without selecting a particular H decay channel. Using the same method, the branching ratio of any (visible) decay channel can be precisely measured and used as input in the WBF H analysis, so that also in this case effects due toc 6 in the Higgs decay can be neglected. Nevertheless, we explicitly checked that taking into accountc 6 effects in the decay for the H → bb channel, which will be the one most precisely measured, results in Tab. 4 are almost unchanged.
As can be seen in eq. (3.7), not only a linearlyc 6 dependent term is present, but also ac 2 6 one. Since C 2 is negative and C 1 is positive for both ZH and WBF H, the SM cross section value is degenerate inc 6 ; besides the SM casec 6 = 0 also a second differentc 6 = 0 condition is giving the same value of the cross section. While for the WBF H this second solution is close toc 6 = 2, in ZH at 240-250 GeV this is aroundc 6 = 9, depending on the energy. As a result, the two solutions being close to each other, in WBF H the 1σ and 2σ intervals are always broad, while in ZH at 240-250 GeV we see two narrow intervals: one aroundc 6 = 0 and one aroundc 6 = 9. Note that for CLIC-350 also ZH is yielding a broad interval as a constraint, since is larger and C 1 is smaller. Via the combined measurement 13 In the case of WBF H at ILC, e.g., only the H → bb has been considered for obtaining the value of .
Thus, smaller values of may be also achieved. in single Higgs production in the Scenario 1 described in the text. Two representative cases are considered: ZH at CEPC-250 and WBF H at CLIC-1400.
of ZH and WBF H processes, or including LHC results in a global fit, thec 6 region around c 6 = 9 can be excluded. In conclusion, assuming no other BSM effects, the best constraints onc 6 via single Higgs production can be obtained at low energy and high luminosity.
We now consider the situation in whichc 6 has a value different from zero, which will denote asc true 6 , and we explore the constraints that can be set onc 6 , by varying the value ofc true 6 . In Fig. 9 we consider ZH at CEPC-250 and WBF H at CLIC-1400 as examples. 14 The bands in the plot show which constraints onc 6 (y-axis) can be set, depending on the value ofc true 6 (x-axis). We considered only the −5 <c 6 ,c true 6 < 5 range, so that results can be directly compared with the analogous analysis performed in the next section for double Higgs production, where this range cannot be extended without violating perturbativity (see Appendix C). The "X" shapes of the ZH and WBF H bands can be understood as follows. In the limit of zero uncertainties two solutions can be obtained from the equation For ZH production, due to the large value of C 1 , only one branch is present in the −5 <c 6 ,c true 6 < 5 region. Instead, for WBF H, since C 1 is small, SM-like scenariosc true 6 ∼ 0 14 The case of CLIC-350 can be directly seen in Fig. 14, which is described later in the text. In this case, since the value of is larger, weaker constraints can be set w.r.t. CEPC-250. 15 If we consider ZH at FCC-ee-240 we obtain P = (4.9, 4.9).  ∼ 4, WBF H constraints are stronger. We remind the reader that it is not obvious that the LHC, even after accumulating 3000 fb −1 of luminosity, will be able to exclude a valuec 6 ∼ 4. Still, with a single measurement forc true 6 ∼ 4 both the intervals aroundc 6 ∼ 4 andc 6 ∼ −3 are allowed, but the latter may be probed also at the LHC. As shown in Tab. 4, also for ZH andc true 6 ∼ 0 there is a second interval in the constraints, but it is outside the range of the plot.

Double Higgs production
We now turn to the case of double Higgs production. The expected precisions for the measurements considered in our analysis 16 are listed in Tab. 5. Although double Higgs production cannot be measured as precise as single Higgs production, it depends onc 6 at LO and therefore the sensitivity on this parameter is much higher.
We start our analysis considering Scenario 1, where we setc 8 = 0. As can be seen in sec. 3.2, the WBF HH dependence onc 6 is similar for different energies. For this reason, for Scenario 1, we show WBF HH only for CLIC-1400, together with ZHH at ILC-500. Similarly to Fig. 9, which concerns the case of single Higgs production, in Fig. 10 we plot the constraints that can be set onc 6 , by varying the value ofc true 6 . Also in σ LO (HH) both a linear and quadratic dependence onc 6 are present, leading to "X"-shape bands. The "X"-shape is slightly asymmetric due to the one-loop σ 30 and σ 40 contributions that are present in σ pheno NLO (HH), see eq. (3.21), which we always use in our study. The central points of the "X" bands are around (c true 6 ,c 6 ) = (−2.5, −2.5) for ZHH at ILC-500, and around (c true 6 ,c 6 ) = (0.5, 0.5) for WBF HH at CLIC-1400. For this reason, although the WBF HH band is narrower due to a largerc 6 dependence, for valuesc true 6 ∼ 0, ZHH at ILC-500 is giving better constraints. On the other hand, for valuesc true 6 = 0 and especiallȳ c true 6 ∼ −2.5, WBF HH at CLIC-1400 is leading to better constraints. It interesting to note that the central points of the "X" bands in WBF H and WBF HH are very close, while for ZH and ZHH they are different. This implies that the combination of the information 16 Note that the value of listed in ref. [49,60] are for a different luminosities than those considered in Tab. 3. Since the statistical uncertainty is the dominant one, the values of in Tab. 5 have been obtained by rescaling those of ref. [49,60] proportionally to the square root of the luminosity. in double Higgs production in the Scenario 1 described in the text. Two representative cases are considered: ZHH at ILC-500 and WBF HH at CLIC-1400.
from WBF single and double Higgs production would not exclude any of the branches of the "X" shape. Thus, the information from ZH or ZHH is necessary for this purpose. We will comment again this point in sec. 4.4.
We now consider Scenario 2. Specifically, we assume that the true value forc 6 isc true 6 and that the measured cross section for double Higgs production is σ measured = σ pheno NLO (c 6 = c true 6 ,c 8 = 0) and we show which value of (c 6 ,c 8 ) can be constrained via the prediction of σ pheno NLO (c 6 ,c 8 ). Starting with the SM case, we show results for σ measured = σ pheno NLO (c 6 = 0,c 8 = 0) in Fig. 11. We consider the range |c 6 | < 5 and |c 8 | < 31, because as explained in Appendix C for larger values the perturbative calculations cannot be trusted. The plot on the left shows the constraints for ZHH and WBF HH at the ILC-500, while the one on the right those for WBF HH at CLIC-1400 and CLIC-3000. First of all we can notice that the constraints onc 6 are weaker than in Scenario 1. Also, no constraints onc 8 independently fromc 6 can be set. On the other hand, the largest part of the (c 6 ,c 8 ) plane can be excluded and the shape of the band depends on the process. It is important to note that this results depend on the choice of the renormalisation scale µ r and therefore the scale at which c 6 (µ r ) andc 8 (µ r ) are measured. Our results refers to µ r = 2m H , which corresponds to the production threshold for the HH pair and therefore to the phase-space region associated to the bulk of the cross section. While the region close to the SM (c 6 ∼ 0,c 8 ∼ 0) is very mildly affected by this choice, we warn the reader that the border of the plane |c 6 | ∼ 5 and |c 8 | ∼ 31 can be strongly affected.
We then consider how the constraints in the (c 6 ,c 8 ) plane depend on the value of σ measured . We consider BSM configurations σ measured = σ pheno NLO (c 6 =c true 6 ,c 8 = 0) with is given. For these plots, only results for ZHH at ILC-500 and WBF HH at ILC-1000 are displayed. Similarly to the SM case, given a value ofc true 6 , the constraints onc 6 independent from c 8 are weaker than those in Scenario 1. However, also in these cases, the largest part of the (c 6 ,c 8 ) plane can be excluded and the shapes of the bands strongly depend both on the process and the value ofc true 6 . In all cases, ZHH and WBF HH sensitivities are complementary; as we will see in sec. 4.4, their combination improves the constraints in the (c 6 ,c 8 ) plane. This is a clear advantage for the ILC, where both ZHH and WBF HH can be precisely measured.
The shapes of the green and red bands can be qualitatively explained as follow. Withoutc 8 effects the green and red bands would simply consist of either two separate (narrow) bands or a single large band, consistently with the results that could be obtained by vertically slicing the bands in Fig. 10. Thec 8 effects bend the bands, leading to the shapes that can be observed in Fig. 12. It is interesting to note that the improvement from CLIC-1400 to CLIC-3000 is rather mild. The main reason is that the increment of the WBF HH cross section is compensated by the decrement of its dependence onc 6 , which can be directly observed in the top-left plot of Fig. 6.

Triple Higgs production
We now consider the case of triple Higgs production. In the SM ZHHH and WBF HHH production processes have a too small cross section for being observed. As an example, if we consider LR-polarised beams at 1 TeV and the dominant decay into a bb pair for the three Higgs bosons and into jets for the Z boson, about 6 ab −1 of integrated luminosity would be necessary for one signal event in the SM. As can be seen in Fig. 8, with WBF HHH the cross section is even smaller in the SM, on the other hand this process has a strong 17 As the total cross section depends onc8 mildly, we do not expect that the constraints depend onc true 0). Indeed, the number of events expected is close to zero and a Gaussian fit cannot be performed. Rather, we have to assume events are zero and compare them with the expected value of events for a given (c true 6 ,c true 8 ) performing a Poissonian analysis. 18 We assume that the other SM backgrounds are giving zero events and we estimate the signal efficiency ε HHH by rescaling the one known for WBF HH production ε HH . In practice, for both WBF HHH and ZHHH production we estimate the signal efficiency to be ε HHH = ε 3 2 HH = 4.7%, where ε HH has been taken from ref. [57]. In Fig. 13, we show the 2σ bounds in the (c 6 ,c 8 ) plane. The plot on the left shows the constraints for ZHHH and WBF HHH at ILC-1000, while the one on the right those for WBF HHH at CLIC-1400 and CLIC-3000. As can be seen, at ILC-1000 almost all the (c 6 ,c 8 ) plane is compatible with a zero event condition, both for ZHHH and WBF HHH production. On the other hand, at CLIC-1400 and especially at CLIC-3000 a vast area of the plane can be excluded via the study of WBF HHH production. In particular, at CLIC-3000, the constraint onc 8 are comparable to those obtainable at a future 100 TeV hadron collider [28,29]. The constraints onc 6 are instead worse than in the double Higgs production case.

Combined bounds
We now investigate the constraints that can be obtained via the combination of the information from single, double and triple Higgs production. We consider both Scenarios 1 and 2 and, as already mentioned, in the case of Scenario 2 we combine only results from double and triple Higgs production. We show in parallel the limits onc 6 from single Higgs production by assuming that thec 8 -dependent two-loop effects are small.
We start discussing the Scenario-1 analysis, separately considering the ILC and CLIC. For both colliders we progressively include results at higher energies in three stages. In for the ILC (left) and the CLIC (right) in the Scenario 1 described in the text. the case of the ILC, we start with ZH at ILC-250, in a second step we include ZHH and WBF H results from ILC-500 and finally ZHHH and WBF H(H(H)) from ILC-1000. Instead, in the case of the CLIC, we start with ZH at CLIC-350, in a second step we include WBF H(H(H)) and ZHHH results from CLIC-1400 and finally WBF H(H(H)) results from CLIC-3000. In the case of triple Higgs production we assume that we observe as many events as predicted by σ LO (HHH) in eq. (3.22), withc 8 = 0.
In Fig. 14, we show the combined results for the ILC (left) and CLIC (right) assuming Scenario 1. In the first stage, both ILC-250 and CLIC-350 constraints are worse than those of CEPC-250 shown in Fig. 9. This is due to a lower precision in the measurements ( ) and for CLIC-350 also a smaller value of C 1 . However, in the second stage, including results at higher energies, for both colliders constraints are much stronger, since double Higgs production becomes available. Especially, combining single and double Higgs production the "X" shape disappears and only the band around the linec 6 =c true 6 remains. 19 In the case of the CLIC, bumps are still present atc 6 ∼ 1, which originate from the centre of the "X"-shape band for WBF H(H) at CLIC-1400, see Fig. 9 and Fig. 10. For the same reason, also for the ILC the band is slightly larger aroundc 6 ∼ 1. In the third stage, constraints are improved both for the ILC and CLIC. Still, the weaker bounds can be set for ∼ 0 <c true 6 < 1, where the center of the "X"-shape band for WBF HH is located. In this region, constraints are better at the ILC thanks to the ZHH contribution at 500 GeV, which helps to resolve this region.
We now consider Scenario 2. As done in the case of double Higgs production we assume that the true value forc 6 isc true 6 and that σ measured = σ pheno NLO (c 6 =c true 6 ,c true 8 = 0) while that for triple Higgs production we observe as many events as predicted by σ LO (HHH) in eq. (3.22), withc 8 = 0. In the case of the ILC, we consider ZHH at ILC-500 and its combination with ILC-1000 results from ZHHH and WBF HH(H) production. In the case of CLIC, we consider ZHHH and WBF HH(H) production at CLIC-1400 and its combination with WBF HH(H) at CLIC-3000. Thus, while ILC-500 is not a combined result, being simply obtained for ZHH production, all the others include information from both double and triple Higgs production. As already said, single Higgs production cannot be directly included in the combination, since itsc 8 dependence starts at two-loop level.
In Fig. 15 we show results for the SM case (c true 6 = 0,c true 8 = 0) as green bands. There we also show as red bands the limits onc 6 extracted from single Higgs measurements at the ILC and CLIC 20 , assuming that the two-loopc 8 dependence is negligible. Due to the available higher energies, combined double and triple Higgs constraints at the CLIC are better than at the ILC. Indeed the WBF HH(H) production cross section increases with the energy. On the other hand, single Higgs production can be better measured at the ILC and therefore the corresponding constraints onc 6 are better than at the CLIC. We notice that the only case where single Higgs results may be relevant in a further combination with those from double and triple Higgs production is the case of ILC-500, which is actually coming from only ZHH production. Indeed, the combination of ZH at ILC-250 and WBF H at ILC-500 would help in removing the band aroundc 6 = −4, and shrinking the possible region for the band around SM value. On the contrary, at higher energies the WBF HH production is more relevant in constrainingc 6 . Thus, with the exception of ILC-500, single Higgs production could be helpful in constraining the (c 6 ,c 8 ) plane only if the dependence onc 8 at two-loop is larger than what we assumed or if low-energy runs at higher luminosity, such as those at circular colliders, are considered.
In Fig. 16 we show the constraints from the combination of double and triple Higgs for BSM casesc true 6 = −4, −2, −1, 1, 2, 4. As already discussed for the SM case, constraints from single Higgs production are negligible for high energy e + e − colliders in this scenario under our assumptions and for this reason they are not shown. We display in each plot both CLIC and ILC bounds. As we can see, both in the SM and in all BSM cases considered, the combination of results from double and triple Higgs production is always strongly improving the bounds. Also, with higher energies, stronger constraints can be set; the best results can be obtained combining results at CLIC-1400 with those at CLIC-3000, especially forc true 6 = 0 since a non-zero number of events can be observed. It is interesting to note that CLIC bounds around (c true 6 ,c true 8 ) are less sensitive than at the ILC on the value ofc true 6 , featuring vertical elongated contours in the (c 6 ,c 8 ) plane. The reason is that at CLIC bounds mainly comes from WBF HHH, while at the ILC mainly from double Higgs production, both ZHH and WBF HH.
In conclusion we observed that low-and high-energy runs are useful for constraining the shape of the Higgs potential. Under the assumption of Scenario 1, we have shown the complementarity of ZH production at low energy with WBF HH information at higher energies. Under the Scenario 2, we have shown that the combination of the information from double and triple Higgs production, which is possible only at high energy, improves the constraints in the (c 6 ,c 8 ) plane (cf. fig. 12 with fig. 16). 20 More specifically, for the ILC, the single Higgs limit are combined results from ZH at ILC-250, WBF H at ILC-500, and WBF H at ILC-1000, while for the CLIC, the single Higgs limit are combined results from ZH at CLIC-350, WBF H at CLIC-1400, and WBF H at CLIC-3000.

Conclusions
Determining whether the scalar potential for the Higgs boson is the minimal one predicted by the SM is among the main targets of the current and future colliders. In this work, we have investigated the possibility of setting constraints on the shape of the Higgs potential via the measurements of single, double and triple Higgs production at future e + e − colliders, considering the two dominant channels, i.e., Z boson associate production (ZH n ) and W boson fusion WBF. In order to leave the possibility for the trilinear and quadrilinear couplings to vary independently, we have added to the SM potential two EFT operators 2 v 2 4 and calculated the tree-level and one-loop dependence on c 6 and c 8 for single and double Higgs production as well as tree-level results for triple Higgs production (see also Tab. 1 in sec.1).
One-loop corrections to single Higgs production, which depends only on λ 3 and thus c 6 , have already been calculated and studied in the literature and we have confirmed previous results. On the other hand, the one-loop dependence on λ 4 and therefore on c 6 and c 8 of double Higgs production has been calculated for the first time here. At variance with the case of single Higgs production, the EFT parametrisation is in this case compulsory and an anomalous coupling approach cannot be consistently used; the c 6 parameter is itself renormalised and receives corrections from both c 6 and c 8 . We have provided all the necessary renormalisation constants and counterterms and expressed the finite one-loop results via analytical form factors that can be directly used in phenomenological applications. We have also motivated the inclusion of the "− 1 2 v 2 " term in the EFT parametrisation, which simplifies the renormalisation procedure by preserving the relations among the SM counterterms. Nevertheless, results can always be easily translated to the In our phenomenological analyses we have considered several experimental setups at future e + e − colliders (CEPC, FCC-ee, ILC and CLIC) and have analysed the constraints that can be set onc 6 ≡ c 6 v 2 λΛ 2 andc 8 ≡ 4c 8 v 4 λΛ 4 . To this purpose we have considered two scenarios: • Scenario 1: the effects ofc 8 are negligible and we analyse the constraints that can be set onc 6 , both for the SM potential and in the casec true 6 = 0.
• Scenario 2: the effects ofc 8 are not assumed to be negligible and we analyse the constraints that can be set on the (c 6 ,c 8 ) plane, both for the SM potential and in the casec true 6 = 0.
In Scenario 1 the value of λ 4 directly depends on λ 3 , while in Scenario 2 they are independent. We verified that requiring perturbative convergence sets upper bounds on the absolute values ofc 6 andc 8 , i.e., |c 6 | < 5 and |c 6 | < 31. Thus, we have analysed the constraints that can be set in this region of the (c 6 ,c 8 ) plane. In Scenario 1, the best constraints onc 6 can be obtained from the combination of ZH results from low-energy high-luminosity runs and results from high-energy runs for ZHH and WBF (HH, HHH) production. On the other hand, in BSM casesc true 6 = 0, WBF H gives stronger constraints than ZH production and similarly WBF HH production can be more sensitive than ZHH production.
In Scenario 2, since two-loopc 8 effects for single Higgs production are not available, we combine only double and triple Higgs production, and show the single Higgs bounds under the assumption that two-loop effects are negligible. The combination of high-energy results from double and triple Higgs production gives the best constraints and in both cases the WBF channel is in general the most relevant. Single Higgs production is only relevant for low-energy machines, and almost negligible once WBF HH is available. For this reason, the higher is the energy, the stronger are the constraints that can be obtained in the (c 6 , c 8 ) plane, both for the SM case and the BSM configurations withc true 6 = 0. In both Scenario 1 and Scenario 2, although WBF HH constraints alone are stronger than those for ZHH, the two production processes are in fact complementary and lead to improved results when they are combined. At high-energy e + e − colliders triple Higgs production is not measurable in the SM, but its cross section strongly depends on the value ofc 8 . In particular, at CLIC-3000, the constraint that can obtained onc 8 via WBF HHH production are comparable to those obtainable at a future 100 TeV hadron collider.
In conclusion, we have demonstrated that the analysis of single, double and triple Higgs production at e + e − colliders can be exploited for constraining the trilinear and quartic coupling via direct and loop-induced indirect effects. In this first sensitivity study we have assumed that BSM effects on the couplings of the Higgs boson with other particles can be neglected. This assumption has already been shown to be reasonable for the SM case in Scenario 1 in ref. [43], yet further studies will be necessary for the other configurations considered in this work. Also, as already mentioned, another possible sensitivity on thec 8 parameter may be obtained from the high-precision measurements of single Higgs production at future e + e − colliders. To establish what kind of constraints could be reached onc 8 in this case, a two-loop computation of e + e − → ZH will be needed. 4.4517.08. This work has received funding from the European Union's Horizon 2020 research and innovation programme as part of the Marie Sk lodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104) and by F.R.S.-FNRS under the "Excellence of Science -EOS" -be.h project n. 30820817. The work of D.P. is supported by the Alexander von Humboldt Foundation, in the framework of the Sofja Kovalevskaja Award Project "Event Simulation for the Large Hadron Collider at High Precision". where Z H is the Higgs wave function and T is the tadpole contribution, which we cancel via the δt counterterm so that the physical value of v does not get shifted. All the other quantities do not receive additional one-loop contributions on top of the SM ones, including δv, which is completely of SM origin. Thus, for our calculation the necessary ingredients for the renormalisation of the virtual corrections are: All the quantities with "SM" as apex are the SM contributions and can be found in [51], those with "NP", which indeed stands for new physics, are the new contributions from c 6 , c 8 and c 10 . Besides c 6 , which is renormalised in the MS scheme, all the other EW input parameters are assumed to be renormalised on-shell, with exception of fine structure α, which we renormalise in the G µ -scheme. This is relevant for our calculation since in the SM the renormalisation of v is related to the charge renormalisation, δZ e , The appearance of the extra quantity −6 c 6 Λ 2 v 3 δv in eq. (A.3) is due to the presence of v in the parametrisation of eq. (2.2), which as we said has an impact in the renormalisation procedure. Before giving the explicit formulas for δZ NP H , δ(m 2 H ) NP , δt NP and the counterterm for the H 3 vertex and H propagator, we briefly discuss this technical aspect.
The explicit term v used in the parametrisation of eq. (2.2) is a subtle quantity. In a tree-level analysis it can be trivially identified with the location of the minimum of V (Φ), which defines the ground state |0 of the Higgs field Figure 18: The structure of one-loop effects in the HHV V amplitude expressed via form factors.
is the contribution from the trilinear Higgs self-coupling to δZ H in the SM and A 0 and B 0 are the standard scalar loop integrals and ∆ is the UV divergence ∆ ≡ 1/ − γ + log(4π) in D = 4 − 2 dimensions. As discussed in sec. 3.2 terms up to the order (v/Λ) 6 have to be in general considered. However, note that no terms beyond (v/Λ) 2 are present in δt, or beyond (v/Λ) 4 in δZ H and δm 2 H , while c 6 is appearing at order (v/Λ) 2 , so terms up to (v/Λ) 6 are in fact present in δc 6 .
We want to stress that all these contributions have to be taken into account in order to obtain gauge invariance for the final finite result of double-Higgs production at one loop. We kept the explicit dependence on the ξ parameter for a generic R ξ -gauge in order to verify that renormalised amplitudes do not depend on ξ. With this calculation setup, results are equivalent to those of a standard calculations based on the parameterisation of eq. (2.3) and c 2n coefficients renormalised in the MS scheme.

B One-loop amplitudes via form factors
In this section we provide all the form factors that are necessary for the calculations of oneloop amplitudes for ZHH and WBF HH production entering σ pheno NLO (HH) in eq. (3.21). These are the form factors for the • HV V vertex, Figure 19: Feynman diagrams contributing to the V [HV V ] form factor at one loop.
where σ 1 and σ 2 are part of σ LO (HH) in eq. (3.10). Note that σ 30 and σ 40 are written in such a form that can be easily extend to the case in which the δZ NP H contribution from external legs is resummed, as done in ref. [36]. However, considering |c 6 | < 5, resummation is not necessary given thatc 2 6 δZ SM,λ H < 4%.

HVV-vertex
The HV V form factor, which will denote as V [HV V ], enters both the single and double Higgs production calculation and can be written as For our calculation thec 6 -independent part can be ignored, while in a generic gauge V 1 [HV V ] is given by the three diagrams 22 in Fig. 19. Using the convention that the corresponding Feynman rule is iV µ 1 µ 2 [HV V ], as we will do also for the other form factors, we can write V µ 1 µ 2

1
[HV V ] as In particular T µ 1 µ 2 (p 1 , p 2 , m V , m H ) = (−6B 0 − 24m 2 V C 0 + 24C 00 )g µ 1 µ 2 − 24p µ 2 1 p µ 1 2 C 12 , (B.6) where p 1 , p 2 are the (incoming) momenta of the two vector bosons, µ 1 , µ 2 are the corresponding Lorentz indices, m V with V = W, Z is mass of the vector bosons, and B 0 , C 0 , C 00 , C 12 are one-loop scalar/tensor integrals defined according to the notation used, e.g., in ref. [51] and where the following variables are understood: It is important to note that P (HH) does not depend onc 8 . Indeed, although the second diagram, the seagull, depends onc 8 due to the HHHH vertex, it is exactly cancelled by the Higgs-mass counter term. We remind the reader that the −δZ NP H component in the counterterm has been removed from P [HH].

HHH-vertex form factor
The form-factor for the HHH vertex, V [HHH], receives contributions from the diagrams already shown in the main text in Fig. 4  ) the value ofc 6 (c 8 ) such that the one-loop amplitude is as large as the tree-level one, i.e. the value of ofc 6 (c 8 ) from where perturbative convergence cannot be trusted anymore. For the estimation ofc max 6 we take into account the leading contribution from V 30 and P 20 , both yieldingc 3 6 terms. For c max 8 we instead consider as first step the contribution from V 11 , which is the dominant term whenc 6 is large, and we compare it with the linearlyc 6 dependent part of the tree-level vertex. In such a way the value ofc max 8 is independent onc 6 .
The value ofc max 6 (c max 8 ) has a kinematic dependence. In the left plot of Fig. 23 we display the dependence ofc max 6 on m(HH), ranging from 125 GeV to 3 TeV. The equivalent plot forc max 8 , taking leading term inc 6 , is shown on the right. Thus, we explore m(HH) values both below the production threshold and in the tail of the m(HH) distributions. As can be seen in both cases, the most stringent constraints, |c 6 | < 5 and |c 8 | < 31, arise from the threshold condition m(HH) = 2m H , while for different values of m(HH) the bound is weaker.
In the case ofc max 6 the constraint is independent on the value of the renormalisation scale µ r and compatible with the result obtained in ref. [56], where the subdominant contribution of P 20 was not taken into account. Conversely, in the case ofc max 8 the constraint does depend on the value of the renormalisation scale µ r . However, we verified for µ r = m H , 4m H , that the most stringentc max 8 value is anyway arising from the kinematic condition m(HH) = 2m H .
The constraints of eq. (C.1) have been derived via the analysis of the H * → HH amplitude, but also the HV V and HHV V vertexes contribute via loop corrections to the quantity σ pheno NLO (HH), eq. (3.21), that is used for our phenomenological analysis. Thus, it is important to check if they can affect the results of eq. (C.1). To this purpose we directly considered the quantities for ZHHH and WBF HH at different energies. The quantity rc 6 is the ratio between the term with the highest power inc 6 from ∆σc 6 and the one with the highest power in c 6 from σ LO , i.e., the ratio of the dominant contributions at tree and one-loop level for largec 6 values. Similarly, the quantity rc 8 is the ratio between the term with the highest power inc 6 from ∆σc 8 and σ LO . Thus, both of them can be considered as a generalisation of the first step; both HV V and HHV V vertices are taken into account and phase-space integration is performed.
In the left plot of Fig. 24, we show rc 6 for the case of ZHH at 500 GeV and of WBF at 1, 1.4 and 3 TeV, which are the phenomenologically relevant scenarios analysed in sec. 4. Requiring |rc 6 | < 1, we can get |c 6 | < 8 for ZHH at 500 GeV, and |c 6 | < 9, 10, 11 for WBF HH at 1000, 1400 and 3000 GeV, respectively. Thus, as one would expected from Fig. 23 for the H * → HH vertex, at higher energies, far from the production threshold, limits are weaker. In the right plot we show rc 8 for the same energies an process. Also in this case the obtained limits are weaker than in eq. (C.1), |c 8 | <∼ 35 − 40.