Testing anomalous $H-W$ couplings and Higgs self-couplings via double and triple Higgs production at $e^+e^-$ colliders

In the present work we study the implications at the future $e^+e^-$ colliders of the modified interaction vertices $WWH$, $WWHH$, $HHH$ and $HHHH$ within the context of the non-linear effective field theory given by the Electroweak Chiral Lagrangian. These vertices are given by four parameters, $a$, $b$, $\kappa_3$ and $\kappa_4$, respectively, that are independent and without any constraint from symmetry considerations in this non-linear effective Lagrangian context, given the fact the Higgs field is a singlet. This is in contrast to the Standard Model, where the vertices are related by $V_{WWH}^{\rm SM}=v V_{WWHH}^{\rm SM}$ and $V_{HHH}^{\rm SM}=v V_{HHHH}^{\rm SM}$, with $v=246$ GeV. We investigate the implications of the absence of these relations in the Electroweak Chiral Lagrangian case, which is a consequence of $a$, $b$, $\kappa_3$ and $\kappa_4$ being generically different than 1, i.e., their particular value in the Standard Model case. We explore the sensitivity to these Higgs anomalous couplings in the two main channels at these colliders: double and triple Higgs production (plus neutrinos). Concretely, we study the access to $a$ and $b$ in $e^+e^- \to HH \nu \bar{\nu}$ and the access to $\kappa_3$ and $\kappa_4$ in $e^+e^- \to HHH \nu \bar{\nu}$. Our study of the beyond the Standard Model couplings via triple Higgs boson production at $e^+e^-$ colliders is novel and shows for the first time the possible accessibility to the quartic Higgs self-coupling.


Introduction
One of the not yet fully explored pieces of the Standard Model (SM) of Particle Physics is the Higgs boson potential, involving the triple and quartic Higgs self-couplings, which are both given in terms of the λ parameter of the potential by λ SM HHH = λ SM HHHH = λ, and whose strength is related to the Higgs boson mass by m 2 H = 2λv 2 , with v = 246 GeV. The implementation of all the Higgs boson interactions within the SM is done usually by setting the Higgs boson as a component of a doublet, the simplest linear realisation of the SU (2) × U (1) electroweak (EW) symmetry. This linear realisation leads to some correlations among the Higgs interaction vertices within the SM. In particular, the triple and quartic self-interaction vertices are related by V SM HHH = vV SM HHHH . Similarly, there are also correlations among other SM couplings in the bosonic sector, like the Higgs interaction vertices with W gauge bosons, which are also related by V SM W W H = vV SM W W HH . Deviations from these correlations may lead to new signals of beyond the Standard Model (BSM) physics. Therefore, an important task at future colliders will be the improvement in the measurement of all these bosonic couplings and the test of their correlations (for a recent review on measurements and constraints on Higgs boson couplings at present and future colliders, see for instance [1][2][3][4][5][6][7][8][9]).
In the present paper we study the sensitivity to these particular four BSM Higgs interactions, W W H, W W HH, HHH and HHHH, and focus on future e + e − colliders. To explore the BSM predictions we work within the context of an Effective Field Theory (EFT) for EW interactions (for a recent review on EFTs for the SM see, for instance, [10]). In particular, we choose the non-linear effective field theory given by the Electroweak Chiral Lagrangian (EChL). This EChL, also called Higgs Effective Field Theory (HEFT) in the literature, was prosed long ago [11][12][13][14][15][16][17][18][19][20][21] to describe the BSM, SU (2) × U (1) gauge invariant interactions without the explicit Higgs field and using a non-linear realisation of the EW chiral symmetry. In the last years, after the discovery of the Higgs boson particle, the EChL has been renewed with the incorporation of the Higgs field. Consequently, the new version of the EChL contains more effective operators that now include also the Higgs boson [22][23][24][25][26][27][28][29][30][31][32][33][34].
In the EChL, the Higgs field has the peculiarity of being a singlet under SU (2) × U (1), in clear contrast with the SM case and other EFTs like the SM effective field theory (SMEFT). As a consequence, the previously mentioned SM correlations between couplings involving two and three Higgs bosons are not present, leading to singular signals in multiple Higgs boson production. Thus, double and triple Higgs productions may exhibit very different rates with respect to the SM. In particular, we focus our studies here on two production channels at e + e − colliders, e + e − → HHνν and e + e − → HHHνν, where the subprocesses initiated by W − W + scattering (WWS, also called Vector Boson Fusion (VBF) in the literature) play a very important role. Indeed, at high energies, in the TeV range, these WWS subprocesses are known to be the dominant ones (for a recent review on VBF, see for instance [35,36]. Within the context of the EChL, the BSM Higgs couplings to bosons we are interested in are parametrised by four independent EChL coefficients, a, b, κ 3 and κ 4 (also called anomalous couplings), one for each of the four BSM interactions: W W H, W W HH, HHH and HHHH. Some of the phenomenological consequences of these couplings via WWS have been examined. The effects of a and b in double Higgs production via WWS have been studied mostly for the case of proton-proton collisions, focusing mainly on the LHC and its future upgrades, both in energy and luminosity [37], although the effects of b at e + e − future colliders have also been explored in [38]. The BSM effects of the triple Higgs coupling, λ HHH = κ 3 λ, at LHC via double Higgs production from WWS have been studied in [39]. The highest sensitivity to κ 3 in hadron colliders is reached, however, via gluon gluon fusion (for a review, see for instance [4]). The quartic Higgs self-coupling, λ HHHH = κ 4 λ, and the consequences of a BSM κ 4 , have barely been explored in the literature. Its effects via radiative corrections have been studied in [40,41], whereas the sensitivity to κ 4 in multi-TeV muon colliders has been considered in [42]. In the present study we explore the sensitivity to the four EChL coefficients a, b, κ 3 and κ 4 via these two e + e − → HH(H)νν processes. In particular, we present a novel study of triple Higgs production at e + e − colliders, that we compare with the better known double Higgs production channel which gives no access to the quartic Higgs couplings, and introduce a new proposal to test both the triple and quartic Higgs self-couplings at once via the e + e − → HHHνν process.  But the measurement of the Higgs self-couplings is not an easy task, as the processes that depend on them involve multiple Higgs production, which is difficult to study and generally suffers of low rates. In the current LHC, with proton-proton collisions at a center of mass (CM) energy of 13 TeV, the SM production cross section of two or more Higgs bosons lies at the picobarn (pb) scale or below [4], and the luminosity is not high enough to produce good statistics. Nevertheless, the triple Higgs self-coupling has been constrained by the ATLAS Collaboration [3] to −2.3 < λ HHH /λ SM HHH < 10.3 at the 95% CL, under the assumption that the possible new physics affects only the triple Higgs coupling. The quartic Higgs self-coupling λ HHHH remains experimentally unaccessible at the moment. A more precise determination of the triple Higgs self-coupling via double Higgs production is one of the aims of several future projects, such as the High Luminosity LHC (HL-LHC) and its high energy upgrade (HE-LHC) [5], the International Linear Collider (ILC) [6] and the Compact Linear Collider (CLIC) [7]. Much research on this topic has already been carried out. The study of triple Higgs production as a way to test new couplings of the Higgs boson to SM particles has already been proposed for future proton-proton colliders, such as the Future Circular Collider (FCC) [8,9]. In this work we will focus on the two electron-positron colliders in the list: the ILC and CLIC. They will operate in several stages with different values of the CM energy and integrated luminosities, as it can be seen in Tab. 1. This work is organized as follows: in Sect. 2 we review the relevant interactions within the EChL for the present work (W W H, W W HH, HHH and HHHH) and write them in terms of the four relevant EChL coefficients a, b, κ 3 and κ 4 . In Sect. 3 we analyze the role of WWS in multiple Higgs production at e + e − colliders. In Sect. 4 we present the computation and results of the cross sections for the two relevant WWS subprocesses, W − W + → HH(H), in the BSM case. We also compare these BSM results with the SM predictions, and check how sensitive they are to variations of the Higgs boson anomalous couplings. In Sect. 5, the corresponding results for e + e − collisions are presented for the various colliders. Finally, Sect. 6 explores the sensitivities to BSM couplings in multiple b-jet events that are expected considering the final Higgs decays. The parameter space region accessible at the future e + e − colliders is also derived. The conclusions are summarized in Sect. 7.

BSM Higgs couplings to bosons within the EChL
To explore the BSM predictions in this work we use the non-linear EFT given by EChL, where the Higgs field is a singlet under both the EW gauge (local) symmetry, SU (2) L × U (1) Y , and the EW chiral (global) symmetry, SU (2) L × SU (2) R . This is in contrast to the usual linear realisations, like the SM itself or other EFTs like the SMEFT, where the Higgs field is one component of an SU (2) doublet. The Higgs singlet is usually introduced in the EChL via polynomial functions, generically leading to uncorrelated BSM Higgs couplings involving one or more Higgs particles, while the three EW Goldstone bosons (GB) that arise from the EW symmetry breaking are usually included in an exponential representation and transform non-linearly under the EW symmetry. We will assume here that the dynamics of the Higgs, Goldstone and EW gauge bosons are given by the EChL, while those describing fermions will be the same as in the SM. The EChL follows a pattern and counting rules [31] similar to those in the chiral Lagrangian for low energy QCD [43]. It contains a tower of effective operators ordered by their chiral dimension, providing predictions ordered in even powers of momentum, p n . It also provides well defined predictions of observables beyond the tree level approximation, i.e. including loop corrections which scale as (p/(4πv)) n , such that the convergence of this momentum expansion is ensured for low energies, say below 4πv ∼ 3 TeV. As in any other EFT of EW interactions, the predictions from the EChL are model independent. The information on the particular ultraviolet theory and the corresponding cut-off that lead to this EChL at low energies is encoded in the coefficients of the effective operators. Therefore, to access the BSM physics via the EChL, the challenge is to find sensitivity to the EChL coefficients at present and future experiments.
For our purposes, in this study of the BSM Higgs couplings to bosons we work at the tree level approximation, so it is enough to use the leading order EChL, which reads: Here v = 246 GeV, g and g are the EW gauge couplings, H is the Higgs boson field, and U is the 2 × 2 matrix that contains the three GBs, ω 1,2,3 : with τ a being the Pauli matrices. The field strengths and covariant derivative are defined as: whereB µ = g B µ τ 3 /2 andŴ µ = gW a µ τ a /2. The rotation from the interaction basis to the physical gauge bosons given by the weak angle θ W is the same as in the SM, so the usual tree level relations among gauge couplings and gauge boson masses still hold. Similarly, the tree level relation m 2 H = 2λv 2 also holds here. The gauge fixing Lagrangian, L GF , and the Fadeev-Popov Lagrangian, L FP , do not enter in our computation of cross-sections, since we work at the tree level and choosing the unitary gauge.
In summary, the relevant interactions in this EChL for the present study, working at the tree level and in the unitary gauge, are given by the four vertices V W W H , V W W HH , V HHH and V HHHH collected in Fig. 1. We see in this figure that the coefficients a and b (yellow and green dots respectively) control the so-called anomalous couplings of the Higgs boson to W and Z bosons, while the coefficients κ 3 and κ 4 (red and blue dots respectively) control the so-called anomalous triple and quartic Higgs self-couplings, respectively. The SM vertices are recovered when a = b = κ 3 = κ 4 = 1. It is worth noticing that the correlations among these vertices within the SM, V SM W W H = vV SM W W HH and V SM HHH = vV SM HHHH , are not present within the EChL. Thus, exploring deviations with respect to these correlations (which occur if any of the EChL coefficients differs from 1) is an interesting way to access the BSM physics at colliders. Regarding the experimental bounds on these EChL coefficients from present colliders, there are available constraints from the LHC for a, b and κ 3 , but not for κ 4 . The ATLAS Collaboration has set the following 95% C.L. bounds on each of the Higgs couplings to gauge bosons using single and double Higgs production: The value of κ 3 has also been constrained by the ATLAS Collaboration, using combined information from single and double Higgs production, but the sensitivity to this parameter is lower in comparison, which leads to a much less restrictive bound than the ones for a and b: Again, we would like to emphasize that currently there are no bounds for κ 4 . For the forthcoming numerical analysis in the present work we will explore the effects of BSM Higgs couplings to gauge bosons and Higgs self-couplings by varying a, b, κ 3 and κ 4 within the following theoretical intervals: Notice that although some of these values are outside the experimental bounds, we believe that it is illustrative to study the behaviour with these parameters when they vary in a wider range. As we will show in the next section, there are indeed further constraints on some of these parameters from the potential violation of unitarity that may occur at high energies for too large values of these coefficients. Concretely, a and b will be further restricted by unitarity, but κ 3 and κ 4 will not. In the last part of this work, all these constraints are taken into account.
3 The role of WWS in multiple Higgs production at e + e − colliders These BSM Higgs couplings modify the predictions of double and triple Higgs production at future e + e − colliders with respect to the SM predictions and can lead to testable departures. A first indication of the size of such departures in the particular e + e − → HH(H)νν channels is provided by the study of the subprocesses W − W + → HH(H) for the case of two (three) Higgs bosons, respectively. These WWS channels are known to be the dominant subprocesses at e + e − colliders with energies at the TeV scale. In this section we will illustrate the WWS dominance in the SM by performing cross sections computations in MadGraph5 (MG5) [44]. We will then compare these results to those predicted by the simplified effective W approximation (EWA). This approximation is very useful, but as we will see it does not provide a sufficiently accurate prediction in all cases.

The dominance of WWS within the SM
Within the SM and neglecting the Yukawa couplings, there are 8 diagrams that mediate the e + e − → HHνν process: 4 of them are of the WWS type, while in the other 4 the neutrino pair is produced via the Z boson. In the e + e − → HHHνν process there are 50 diagrams in total, 25 of which are WWS and 25 Z-mediated. We do not draw all these diagrams here, for shortness, but the subset of diagrams that are mediated by WWS are easily extracted from those of the corresponding subprocess, which have been collected in the appendices. MG5 takes into account all these diagrams, allowing to distinguish the different contributions when computing total cross sections. In Fig. 2 we show these cross sections, taking into account all diagrams (dark blue lines) and only the Z-mediated ones (light blue lines). From Fig. 2 we learn that the dominant contribution to σ(e + e − → HH(H)νν) at low energies comes from the associated production of two (three) H's with a Z boson which then goes to νν, denoted here as "σ(e + e − → HH(H)νν) via ZHH(H)". The cross section for this particular production mechanism can be computed as the cross section of the process e + e − → ZHH(H) times the branching ratio BR(Z → invisible) = 20%. As the energy increases, the associated production of two (three) H's with a Z boson becomes subdominant, which suggests that HH(H)νν production at the TeV scale is dominated by the other diagrams (WWS). However, to estimate the size of the WWS contributions one cannot naively extract these diagrams, since they generically do not form a gauge invariant subset. Thus, a different method must be used. Fig. 3, shows that the enhancement in σ(e + e − → HH(H)νν) at high energies actually comes from σ(e + e − → HH(H)ν eνe ), that is, processes with electron neutrinos in the final state. This fact becomes clearer when we also represent the contribution of σ(e + e − → HH(H)ν µνµ ), which decreases at high energies. 1 . The main difference between both processes is that  whereas e + e − → HH(H)ν eνe includes both types of diagrams, WWS and ZHH(H), this is not the case in e + e − → HH(H)ν µνµ where all diagrams are of ZHH(H) type. Then, the contribution to σ(e + e − → HH(H)ν eνe ) coming from WWS can be isolated using the following "theoretical observable": Therefore, this difference between the channels with electron and muon neutrinos quantifies the weight of the WWS diagrams in the full process without the need of separating diagrams which, as mentioned above, is not the proper way if one wishes to preserve gauge invariance. We also define the related quantity R WWS , This is an adimensional quantity that clearly determines how large is the contribution of WWS in the e + e − → HH(H)ν eνe process. If R WWS is close to 1, the cross section will be mostly dominated by WWS. The isolation of these WWS contributions is of relevance for the present study, since, as we will see next, they are indeed the most sensitive ones to the anomalous Higgs couplings. The importance of the WWS contribution in the two processes, e + e − → HHν eνe and e + e − → HHHν eνe , is clearly illustrated in Fig. 4, where we have included the predictions for σ(e + e − → HH(H)νν), σ(e + e − → HH(H)ν eνe ) and σ WWS as defined in Eq. (7). Above the O(1 TeV) energies, σ WWS clearly approaches σ(e + e − → HH(H)ν eνe ), which in turn approaches the total σ(e + e − → HH(H)νν). Thus, we conclude on the dominance of WWS in σ(e + e − → HH(H)νν), via the particular channel with electron neutrinos.
It is also important to notice that ZZ scattering (ZZS) contributes to the final state HH(H)e + e − , which has a much lower production cross section (the probability of producing two or three Higgs bosons via ZZS is approximately ten times lower). Therefore, this ZZS will not be relevant in this work.
Although the previous figures show separately the various contributions from the different neutrino species, σ(e + e − → HH(H)ν iνi ), one must realize that in the experiment it is only possible to measure their sum, σ(e + e − → HH(H)νν). However, the key point on studying such theoretical observables is to show, in a gauge invariant way, that WWS is indeed the largest contribution at high energies and. As a consequence, the predictions for σ(e + e − → HH(H)νν) will approximately follow the same behaviour at these large energies, specifically that of σ(e + e − → HH(H)ν eνe ).
Comparing triple with double Higgs production in Figs. 2, 3 and 4, we find that they are roughly similar in shape, being the WWS enhancement in triple Higgs production displaced to higher energies, since one extra particle is produced. This causes ZHHH to be more relevant compared to HHHνν than ZHH compared to HHνν, especially at energies below 2000 GeV. It is also important to note that the SM cross sections for triple Higgs production at the TeV energy scale are typically three orders of magnitude below those for double Higgs production, which is why it is not expected to measure SM-like HHH production in future linear colliders. For this reason, we will focus our attention on BSM scenarios.

The effective W approximation versus full MG5 simulations
Before exploring the effects of BSM Higgs couplings, it is instructive to first evaluate the validity of the so-called effective W approximation (EWA) [45] by comparing its predictions to the full MG5 simulations. The EWA is a generalisation of the well known effective photon approximation, and treats vector bosons, W 's and Z's, that are radiated from the initial fermions, as if they were partons inside the fermions. This allows to define distribution functions for these vector bosons in a similar way as the PDFs of quarks and gluons inside a proton. The EWA also assumes that the vector bosons are radiated colinearly by the fermions, and then they scatter on-shell in the subprocess. This allows to employ factorisation, obtaining the total cross section of the whole process by convoluting that of the subprocess with the distribution functions of the vector bosons, providing a more analytical and intuitive approach to the computation than with MG5. The key point then is that the EWA assumes that the WWS is the dominant subprocess and, therefore, the validity of this approximation will depend on how much accurate this assumption is. This procedure is much simpler than computing the full cross section by a Monte Carlo method like MG5. In particular, the EWA approximation has been used in the literature very often to simplify the estimates of BSM physics in colliders, both for e + e − and pp collisions, via WWS. For the present case of e + e − colliders the generic representation of the WWS participating in the collision and producing multiple (two or three in our case) Higgs bosons is drawn in Fig. 5. The simple formula displaying the mentioned factorisation is: Here, σ(s) = σ(e + e − → HH(H)ν eνe ) is the total cross section of the process of interest at a CM energy of √ s, andσ ij (ŝ) =σ(W i W j → HH(H)) is the cross section of the WWS subprocess at a CM energy of √ŝ . x 1 and x 2 are the momentum fractions carried by each W boson and define the CM energy of the subprocess byŝ = x 1 x 2 s. The subindices i, j refer to the polarization of the W bosons (longitudinal or transverse). Different polarizations must be taken into account separately, as the probability of radiating a W boson depends on whether it is longitudinally or transversely polarized. Consequently, each polarized cross section is convoluted with the corresponding combination of distribution functions f i (x). Note that this formula assumes that WWS is the dominant contribution to σ(e + e − → HH(H)ν eνe ), so it is expected to work better at high energies. To compute this cross section we writê σ ij (ŝ) in terms of the polarized amplitudes M ij , which we generate using FeynArts-3.10 [46] and FormCalc-9.6 [47], and then perform the integration using VEGAS [48] and a private PYTHON code. We do it for both the HHν eνe and HHHν eνe channels, with the corresponding phase space factors, which we omit here for shortness. The analytical expressions we use for the W distribution functions are taken from Ref. [45] and correspond to the so-called improved EWA, that keeps corrections of order m 2 W /E 2 , with E being the energy of the parent fermion radiating the W . This improved EWA works better than the most frequently used Leading Log Approximation (LLA) EWA, which is only valid in the very high energy limit, E m W . The formulas for the LLA-EWA can also be found in Ref. [45]. The Feynman diagrams contributing to the scattering amplitudes of the W − W + → HH(H) subprocesses are collected in the appendices. We omit the corresponding analytical expressions for shortness. The numerical results of our estimates of the cross sections with the improved EWA and their comparison with the MG5 results are presented next. In Fig. 6 we show the total cross section of the two processes, e + e − → HHν eνe (left) and e + e − → HHHν eνe (right,) as a function of the e + e − CM energy. We include the predictions from the improved EWA and from MG5. We also include the predictions of the contributions to the total cross section from WWS, σ W W S , as defined in Eq. (7). We see in this figure that at high energies the EWA predicts the HH production cross section with good accuracy, while it fails in the predicted cross section in more than a factor of 2 for the HHH case. The energies that we are considering in this work seem to be too low for the EWA to be a good approximation of triple Higgs production. We have also checked that the dominant contribution to the total cross section comes from longitudinally polarized W bosons in both cases.
Another condition for the EWA to be a useful approximation is that it can reproduce the differential cross section distributions. As an example, Fig. 7 shows the differential cross section with respect to the invariant mass of the final state Higgs bosons, M HH and M HHH respectively, which is equivalent to √ŝ due to 4-momentum conservation. The results for the SM case are again accurate for double Higgs production, but not so satisfactory in the case of three Higgs bosons. The good agreement of the EWA with the MG5 computation in the double Higgs production channel also occurs in the case of the BSM predictions. Indeed, the EWA works even better for BSM than for the SM. This will be shown in the next section.
From this brief analysis we conclude that the EWA is an interesting alternative to compute our results in the double Higgs case, but it is not accurate enough to study triple Higgs production in the range of energies available at the e + e − colliders under consideration. For the final analysis at e + e − colliders in sections 5 and 6, we will not use the EWA anymore, and our results on total and differential cross sections will be extracted from the full Monte Carlo simulation with MG5.

Deviations from BSM Higgs couplings in WWS
In this section we study the effects of the anomalous Higgs couplings in double and triple Higgs production via WWS. In particular, we consider the W W H, W W HH, HHH and HHHH interactions, given by the EChL coefficients a, b, κ 3 and κ 4 , and explore the departures with respect to the SM predictions both in total cross sections and in some differential distributions. First, we study the effects of a and b in W − W + → HH and next we analyse the effects of κ 3 and κ 4 in W − W + → HH(H). To simplify the computations, we take κ 3,4 = 1 when exploring the sensitivity to a and b, and a, b = 1 when studying the effects of κ 3 and κ 4 . The computations presented in this section have been performed using FeynArts-3.10 for the generation of diagrams in the unitary gauge and FormCalc-9.6 to perform the analytical calculation of the corresponding amplitudes. In addition, we use VEGAS to integrate numerically over the corresponding phase space in order to obtain the corresponding cross sections. All these computations have been additionally checked with MG5.

Effects of a and b in W − W + → HH
Here we present the effects in the W − W + → HH subprocess of the two anomalous Higgs couplings parametrised by a and b, and compare them with the SM case, i.e, with a = b = 1. Fig. 8 shows how the total cross section of this subprocess depends on the CM energy when varying the EChL parameters a and b. We restrict our study to positive values of a but consider both positive and negative values of b. Although some values are outside the experimental bounds, we believe that at this intial stage it is illustrative to show the effects of these parameters when they vary in a wider range. Some first conclusions can be extracted from Fig. 8. It is clear that the behaviour of the cross section with the energy in the EChL is very different than that of the SM prediction. While in the SM case (dashed yellow line in the left plot), the behaviour is flat with the energy at the TeV region and above, in the rest of the cases, with a and/or b different from 1, the cross section grows very steeply. This growth with the energy occurs typically in predictions from EFTs, due to the ultraviolet incompleteness of these theories. In particular, the EChL provides predictions, as previously said, that generically grow with powers of the external momenta.
In the SM case, the flatness with the energy occurs because there is a strong cancellation of the terms that grow with the energy among the contributions from the various diagrams. If we consider the dominant contribution to this SM cross section, that comes from the longitudinally polarized W gauge bosons, it turns out that, at high energies (for √ŝ m W , m H ), the contributions from the contact, t and u channels to the scattering amplitude M(W L W L → HH) are proportional toŝ, whereas the contribution from the s channel has a constant behaviour withŝ [39]. In this case, i.e for a = b = 1, there is an exact cancellation of the terms growing linearly withŝ among the contact, t and u channels, and what remain are just the terms having a constant behaviour withŝ. However, this strong cancellation among diagrams does not happen in the case of the EChL for arbitrary a and b values. More concretely, this cancellation of the terms growing with energy occurs in the SM because of the relation mentioned previously among the two vertices, V W W H = vV W W HH , which is not fulfilled in the EChL case. To get a similar cancellation of the terms growing linearly withŝ in the EChL amplitude one should restrict the parameters to the particular setting given by b = a 2 , which is obviously not necessarily true in the general case. Notice also that this restriction on a and b is not required by any symmetry argument (remember that the EChL is gauge and chiral invariant for all a and b). Within the SM, the relation V W W H = vV W W HH occurs as a consequence of the symmetry being linearly realised, i.e. because the Higgs field is placed into a doublet. But in the EChL, the Higgs field is a singlet under the gauge and chiral symmetries and a and b can be arbitrary. The two plots in Fig. 8 also show that, the bigger the deviations of a and b from 1, the bigger the cross section. Thus, some possible BSM physics, given by the EChL parameters departing from their SM values, would yield clear experimental evidence, with much larger cross sections than those predicted by the SM. For instance, at √ŝ = 3 TeV, and for the considered values of a and b in the left panel, the cross section can be several orders of magnitude larger than in the SM. For negative values of b (right plot), the departure with respect to the SM prediction can be even larger. For instance, for a = 1, b = −1 the cross section is about a factor 1000 larger than that for a = b = 1. This is a remarkable enhancement, produced by just changing the sign of b. Notice that the variations chosen for b with respect to 1 are much bigger than for a, simply because the experimental bounds existing on a are more restrictive, while the constraints are looser for b. In these plots, a varies at most a 10% with respect to 1, and b does so in up to 300%; these very different ranges are responsible for the apparent dominance of b.
This large increase of the cross sections when a and/or b are not equal to 1 might indicate that unitarity is not guaranteed in the EChL predictions. In order to check if unitarity is preserved, we plot, in Fig. 9, the energy dependence of the longitudinally polarized s-wave partial amplitude for the same values of a and b as those chosen in the previous figure.
These results show that unitarity is, in fact, violated for some parameter (a, b) values at high energies. Except in the SM case (which, as it is well known, perfectly respects unitarity), the rest of the chosen values of the parameters render partial wave amplitudes larger than 1 at energies of several TeVs, which could be reachable at colliders. As it could be anticipated, the further the departure from the SM, the sooner unitarity is violated.
Focusing on one particular energy, it is also very illustrative to study the failure of unitarity as a function of the parameters a and b. Fig. 10 (left panel) shows the predictions at √ŝ = 3 TeV for the lowest partial wave amplitude as a function of b and for several values of the parameter a. It can be seen that, at this energy, when a varies from its SM value approximately a 10% (which, roughly speaking, coincides with experimental bounds), b can depart from 1 a 50% at most, in order to obtain predictions which preserve unitarity. Thus, for the case of CLIC with the largest planned energy, we learn from this figure that by keeping a within the experimentally allowed interval and by restricting the b parameter in a rather conservative interval of b ∈ [0.5, 1.5], the predictions will respect unitarity. With this conservative choice, unitarity will consequently be also preserved at other colliders with lower energies.
Finally, to illustrate the effect of the EChL parameters on the issue of unitarity violation when restricted by the relation b = a 2 , we show in Fig. 10 (right panel) the predictions for the lowest partial wave amplitude as a function of the energy, for various values of a (b is set here to a 2 ). It is clear from this plot that, at high energies ( √ŝ m W , m H ), a constant behaviour with the energy, similar to the one in the SM, is obtained in the EChL when the a and b parameters are restricted by b = a 2 . We also conclude from this figure that all the displayed predictions respect the unitary bound, expect for the extremely large value of a = 10, where the constant prediction excedes the unitary limit. Such large values of a are totally out of the experimentally allowed interval, and thus they will not be considered anymore in this work.
All in all, the results in this subsection clearly show that the EChL predictions are way different than those of the SM when varying the parameters a and b and provide quite large cross sections for the W − W + → HH subprocess. We will see in the following that this can lead to measurable departures from the SM predictions at future e + e − colliders. And this is true even after including the restrictions on the parameters a and b such that unitarity is respected in all the predictions.

Effects of κ 3 and κ 4 in W − W + → HH(H)
In this section we analyse the WWS subprocesses and the BSM deviations with respect to the SM predictions arising from the anomalous triple and quartic Higgs self-couplings, parametrised by κ 3 and κ 4 respectively. For this study we explore channels of both double and triple Higgs production and, for simplicity, we set a, b = 1. Since double Higgs production is not sensitive at all to κ 4 , we focus here mainly on the triple Higgs production case, which is the novel one, and use the double Higgs case rather as a reference case (better known) regarding the sensitivity to κ 3 . For a more devoted study on the effect of κ 3 in double Higgs production see, for instance, Ref. [39].
The results of our calculations are collected in Figs. 11 to 13, in which we show the behaviour of the cross section with the energy for different values of κ 3 and κ 4 , and Figs. 14 to 16, in which we display the cross section variation with κ 3 and κ 4 , independently, at fixed energies. Notice again that since double Higgs production does not involve the quartic coupling (at least at tree level), the value of κ 4 is not relevant for the predictions in that channel. We can extract some first conclusions by looking at these plots. Starting from Fig. 11, that shows the double Higgs production case, we observe that modifying κ 3 in the interval [-10,10] leads to a notable enhancement in the total cross section with respect to the SM prediction, here represented by the dashed line (κ 3 = 1). The maximum deviation occurs slightly above the threshold energy 2m H , and it is larger for negative values of κ 3 . In the most extreme case, κ 3 = −10, the BSM prediction can deviate up to two orders of magnitude with respect to the SM value in this close to threshold energy region. Notice also that there are no regions disallowed by unitarity in this figure. It can be checked from Ref. [39] that for κ 3 ∈ [−10, 10] and for the energies considered here, √ŝ ≤ 3 TeV, unitarity is preserved in all the predictions for double Higgs production. Fig. 12 shows the case of triple Higgs production. Comparing this figure with the previous Fig. 11, we can see that the consequences of modifying κ 3 , with κ 4 fixed to 1, are similar in triple and double Higgs production, but now the maximum appears slightly above the new threshold energy, 3m H . The deviations with respect to the SM prediction are larger again for negative values of κ 3 , and can reach values of up to five orders of magnitude larger than the SM in that close to threshold energy region. On the other hand, we learn from Fig. 13 that varying the value of κ 4 within the interval [-10,10], with κ 3 fixed to 1, does not modify the shape of the cross section significantly. However, it can also increase its value by one or even two orders of magnitude with respect to the SM prediction in the most extreme cases. The dependence on the sign of κ 4 in triple Higgs production is very mild, in contrast with the  strong dependence on the sign of κ 3 that has been found in both channels.
Coming back to W − W + → HH, Fig. 14 shows the variation of the cross section with κ 3 at a fixed CM energy. Here we see that there is a minimum in the region κ 3 ∈ [2, 5]. The  greater the energy, the larger the value of κ 3 that minimizes the cross section. Deviations from the SM can reach two orders of magnitude in the most extreme case, that is, for κ 3 = −10 at the lowest energy. The sensitivity to variations of the κ 3 parameter decreases with increasing energy; note that far from the region of the minimum, the behaviour of the cross sections is inverted, being higher at lower energies. This is the reason for the appearance of the peak near the threshold in Fig. 11, which does not appear in the SM case. When reproducing this same plot for W − W + → HHH as a function of κ 3 (see Fig. 15) we find a very different picture. First of all, in the left plot in Fig. 15 we see that the maxima at the extreme values of κ 3 = ±10 are displaced in energy and they are now reached around √ŝ = 500 GeV. This is in contrast to the double Higgs production case, where these maxima are produced at the minimum energy of √ŝ = 380 GeV. Second, the shape of the curve, especially at low energies, exhibits two minima , in contrast to just one minimum in HH.
One of them is near the SM value κ 3 = 1, and the other one appears at a positive κ 3 and is displaced to higher values of this parameter as energy grows. The maximum deviations with respect to the SM are large in any case, varying from two orders of magnitude at √ŝ = 3000 GeV to even five orders of magnitude at √ŝ = 380 GeV. One of the most interesting findings is that the cross section σ(W − W + → HHH) is significantly more sensitive to variations in the κ 3 parameter than σ(W − W + → HH), mainly at low energies. However, triple Higgs production has the disadvantage that it leads to smaller cross sections than double Higgs production, due to the obvious phase space suppression.
Finally, Fig. 16 shows the dependence with the κ 4 parameter, fixing κ 3 = 1, for various values of the CM energy. We see that, in this case, the predictions of the cross sections only exhibit one minimum for each energy. Besides, the maximum deviations from the SM vary from one order of magnitude at √ŝ = 3000 GeV up to three orders of magnitude at √ŝ = 380 GeV.
Since triple Higgs production depends on both parameters, κ 3 and κ 4 , it is also interesting to check what happens if we vary the two of them at the same time. In Fig. 17 we represent the W − W + → HHH cross section at three fixed energies, √ŝ = 500, 1000 and 3000 GeV, in the (κ 3 , κ 4 ) plane. The additional information that we can extract looking at these plots is that variations in the cross section are strongly dominated by deviations along the κ 3 direction, and the higher values are reached when κ 3 approaches −10. The dependence with κ 4 is softer, and the maximum cross section can be reached either at negative or positive values of κ 4 . For example, in the negative κ 3 region, the cross section tends to be higher around κ 4 = 10, while in the positive κ 3 region the trend can be this same or the opposite one, depending on the energy. Note that at high energies the combined effect of modifying both parameters at the same time can lead to a notable increase of the cross section, especially if they have opposite signs. For instance, at √ŝ = 3000 GeV and for (κ 3 , κ 4 ) near the corner (−10, 10), the cross section exhibits large values of O(100) pb, to be compared with the SM prediction of 3.7 × 10 −2 pb. We also notice that the lowest values (darker region) are approximately arranged along a line which does not coincide with the κ 3 = κ 4 direction. This means that modifying the Higgs self-couplings without altering their ratio, that is, λ HHH /λ HHHH = λ SM HHH /λ SM HHHH = 1, can also produce an enhancement in the cross section. The size of this darker region is related to the depth of the minimum, and is smaller for lower energies. Note that the SM prediction is contained in this region, which is why we do not expect to measure this triple Higgs production process if the self-couplings are close to their SM values. In any case, what we are studying here is the subprocess, which cannot be seen in a real experiment, so to complete our analysis, in the next section, we will study the full e + e − process to understand how sensitive it is to deviations in these parameters.
Regarding unitarity, we find that the W − W + → HHH subprocess does not violate unitarity for the energies and parameter values considered here, i.e. κ 3,4 ∈ [−10, 10] and √ŝ < 3 TeV. We also confirm the good unitarity behaviour in double Higgs production with respect to κ 3 , for these same settings. It is worth mentioning that for triple Higgs production, the unitarity check has been performed with a different criterion than for double Higgs production. This difference is because to study the unitarity bounds for a 2 → n process, one cannot employ the usual partial waves expansion, since the amplitude depends on a higher number of variables, other than s, t and u. Following for instance Ref. [9], one way of checking unitarity in a 2 → n process is to insert a complete set of intermediate states into its left-hand side, separating the elastic and the inelastic parts: Here Π n refers to the n-body phase space. In this case, the following bound for σ inel (2 → n) is obtained after introducing the partial wave expansion for M el (2 → 2): This is the constraint that we have imposed to check the unitarity of our predictions for σ(W − W + → HHH). According to this criterion, all our predictions for BSM physics coming from deviations in κ 3 and κ 4 are fully unitary for all the energies considered in this work. Concretely, in the worst-case scenario, i.e. for the largest energies considered here ( √ŝ = 3000 GeV), preservation of unitarity requires that: which is not reached even in the most extreme cases in Fig. 17. Thus, in contrast to a and b, which suffer of strong restrictions from unitarity, the parameters κ 3 and κ 4 do not.

Deviations from BSM Higgs couplings in e + e −
In this section we present the results for the effects of the BSM Higgs couplings within the EChL parametrized by a, b, κ 3 and κ 4 in double and triple Higgs production at the real scattering processes e + e − → HHν eνe and e + e − → HHHν eνe , respectively. We will explore these effects at the planned e + e − colliders, and compare the BSM behaviour with the SM predictions. Again, to simplify this study, we take κ 3,4 = 1 when exploring the sensitivity to a and b, and viceversa. As in the SM case, the full computation of the BSM rates is done with MG5 (we neglect the Yukawa couplings to the electrons) and includes 8 contributing diagrams in e + e − → HHν eνe , 4 mediated by WWS and 4 by Z, and 50 contributing diagrams in e + e − → HHHν eνe , 25 mediated by WWS and 25 by Z. We do not draw all of them here, for shortness, but the ones mediated by WWS are easily extracted from those of the corresponding subprocess, which have been collected in the appendices.  Table 2: Predictions of the cross sections σ(e + e − → HHν eνe ) (in pb) within the EChL for different values of the parameters a and b and e + e − CM energies. The SM case (a = b = 1) is also included for comparison. The predictions from both MG5 and the EWA are displayed, as well as the ratio R WWS , defined in Eq. (8), which quantifies the relevance of the WWS within the whole e + e − process.
see from the predictions with MG5 that the distortions introduced by a and b with respect to the SM can be very large, leading to large enhacements in the cross sections. If we set one of the parameters to 1 and vary the other one in a certain amount, the enhancement will be larger when the variation is applied on a rather than on b (notice the difference between the cases a = 1, b = 1.5 and a = 1.5, b = 1). This could be expected since the a parameter enters quadratically in the amplitude whereas b enters linearly. We also see in this table that the size of the enhancements obtained for a given ∆a and ∆b depends on the energy and grows in general with it. The maximum departures from the SM are found for the largest energies and ∆'s considered here, and can lead to ratios of BSM cross section over SM cross section of up to O(100). Second, from this table we also learn that the EWA gives quite good predictions for high energies at and above 1 TeV, not only for the SM rates as we have already seen in Sect. 3.2, but also for the BSM rates. Indeed, we see that the accuracy of the EWA is even better for the BSM predictions than for the SM ones, and this occurs because the dominance of the WWS subprocess is more pronounced in the BSM scenarios than in the SM case. This is illustrated by the ratio R WWS , defined in Eq. (8), which clearly approaches one at the highest energies, indicating that WWS largely dominates the full BSM cross sections. This  WWS dominance can also be seen in the predictions of the cross sections distributions with the invariant mass M HH , as it is shown in Fig. 18, where the collider energy has been set to 3 TeV. We also see in this figure that the prediction from the EWA is very close to the prediction from MG5 for M HH above 1 TeV. All in all, we learn from Tab. 2 and Fig. 18 that the main features found in the predicted cross sections for e + e − → HHν eνe with respect to variations in the a and b parameters follow a similar pattern as those found previously at the W W → HH subprocess level in Sect. 4.1. The reason for this similarity is again the dominance of the WWS subprocess, which is even more pronounced in the BSM scenarios than in the SM case. Finally, we present the results by scanning both parameters a and b together in the (a, b) plane and for various collider energies. Fig. 19 shows the cross section contour lines, obtained with MG5, for e + e − collider energies of 500 GeV, 1 TeV and 3 TeV. The range for the b parameter in these plots, b ∈ [0.5, 1.5], is chosen such that unitarity is guaranteed for the three considered centre-of-mass energies. Although the range considered for a is larger than the one it is actually constrained to, we wanted to display the same ranges for both parameters in order to obtain a global view. Notice that in these plots we have also marked the region with a > 1, where the so-called positivity constraint applies. According to Refs. [49] and [50], some values in this region could yield a potential problem due to violation of causality in this EFT. We include in these plots a dotted region, that is theoretically disfavoured by positivity, and and orange band that displays the experimentally allowed region for a, within a 95% confidence level. We have checked that within this experimentally allowed region for a unitarity is fully preserved.
The following plots finally provide solid conclusions on the sensitivity of the e + e − → HHν eνe process to variations of the a and b EChL parameters. As both a and b are modified simultaneously and in the same range, it is possible to compare their role on the behaviour of the cross section. Notice that the plots are not symmetrical with respect to the (a, b) = (1, 1) point, which means that increasing any of the parameters with respect to 1 in a certain amount is not equivalent to diminishing it in the same amount. This happens for both parameters: the plots do not exhibit any symmetry with respect to the a = 1 or b = 1 axes. This means that the cross section is sensitive to the sign of a − 1 and b − 1. This was already observed in Fig. 18, and is now confirmed. Also, the plots are not symmetrical with respect to the a = b line: equal variations in a and b are not equivalent. Clearly a is the dominant parameter: cross sections grow faster when keeping b constant and varying a than viceversa, which could be expected due to the quadratic dependence on a in the amplitude.
The most relevant conclusion that we learn from Fig. 19 is that the BSM rates are considerably larger than the SM ones, especially for the higher energy colliders, even if we limit ourselves to the region of the (a, b) parameter space allowed by unitarity, positivity and present experimental constraints. Thus, it is interesting to perform a more detailed analysis in a collider framework, taking into account the Higgs bosons decays. We will carry out such an analysis in Sect. 6.

Deviations from κ 3 and κ 4 in e + e − → HH(H)ν eνe
As we have said, the cross section of e + e − → HHν eνe is only sensitive to κ 3 , whereas that of e + e − → HHHν eνe is sensitive to both parameters κ 3 and κ 4 . The behaviour of the cross section with the energy for e + e − → HHν eνe and for different values of κ 3 is shown in Fig. 20. The behaviour of the cross section of e + e − → HHHν eνe for different values of κ 3 and κ 4 is shown in Fig. 21 and Fig. 22, respectively. In Fig. 23 the dependence with κ 3 of σ(e + e − → HHν eνe ) at various fixed energies is represented. The dependence with κ 3 and κ 4 of σ(e + e − → HHHν eνe ) at various fixed energies is displayed in Fig. 24 and Fig. 25, respectively. In all cases we include the SM prediction, corresponding to κ 3 = κ 4 = 1, for comparison. In a first look to all these plots we see clearly that the main features found at the WWS subprocess level in Sect. 4.2, for W − W + → HH and W − W + → HHH, are again found here at the collider level, for e + e − → HHνν and e + e − → HHHνν, respectively. In the following we comment in more detail the results in each of these figures. First, we start with the discussion of HH production in Fig. 20. In these plots we can see how, in general, deviating from κ 3 = 1 causes an enhancement of the cross section that is approximately constant with energy. The strongest deviation occurs when κ 3 = −10, and it differs from the SM prediction by two orders of magnitude. Another thing that can be seen in HH production and will be more significant in HHH is that the bump near the threshold (2m H in HH and 3m H in HHH), where the associated Z subprocess dominates, disappears when we deviate from the SM (dashed lines).
Second, comparing the behaviour of HH to the HHH case, which is presented in Fig. 21, we observe something that we had already noticed in the previous section: triple Higgs production is extremely sensitive to variations in the κ 3 parameter, much more than double Higgs production, reaching BSM deviations of even five orders of magnitude with respect to the SM prediction in the most extreme case (κ 3 = −10). We also learn from this Fig. 21 that in the HHH channel, similarly to the HH channel, the difference with respect to the SM prediction is approximately constant with energy. The bump that is dominated by the associated Z production subprocess disappears as we separate from κ 3 = 1, showing once more that the BSM deviations are clearly dominated by WWS. Fig. 22 is the equivalent to Fig. 21, this time fixing κ 3 to 1 and varying the κ 4 parameter. The profile of the deviations due to κ 4 is very similar to those correspoding to κ 3 , but softer. The maximum deviation now occurs for κ 4 = −10, this time yielding rates around two orders of magnitude above the SM prediction. As we noted previously, BSM deviations mediated by the ZHHH subprocess due to κ 4 = 1 are much smaller than the ones coming from WWS.
Third, we can now look at the dependence with the κ 3 and κ 4 parameters at a fixed CM energy. Starting with the HH case shown in Fig. 23, the main difference that we see when comparing this plot with the results for the subprocess in Sect. 4 (see Fig. 14) is that in the e + e − case the highest cross section for all κ 3 values is always achieved at the highest collider energy. The second observation is that there is not a large difference in the sensitivity to κ 3 depending on the energy, confirming, as stated above, that the deviations with respect to the SM do not depend appreciably on the energy. All the curves in Fig. 23 for the various energies experiment a variation between one and two orders of magnitude with respect to the SM, having the maximum at the extreme value of κ 3 = −10, and they all have just one minimum at the region κ 3 ∈ [0, 2] (the higher the energy, the larger the value of κ 3 at the dip).
In Fig. 24 we show the corresponding plot for HHH production, where we see the dependence of the cross section with κ 3 , setting κ 4 = 1, at various CM energies. Here  we find again, similarly to the HH case, that the highest cross section for all κ 3 values is always achieved at the highest collider energy. Also the deviations with respect to the SM prediction are practically insensitive to the energy, and the largest BSM cross section is obtained again at κ 3 = −10, reaching values of up to more than three orders of magnitud higher than the SM value. In contrast to the previous HH case, the triple Higgs production shows a new feature, namely, the appearance of two minima instead of one. These two minima are obviously correlated with the two minima already observed in Fig. 15 at the W − W + → HHH subprocess level, and they are more clearly visible at lower collider energies.
In particular, for √ s = 380 GeV, there is one minimum around κ 3 = 1 and the other one is between κ 3 = 2 and κ 3 = 3. We also observe a deformation in the √ s = 1000, 1500 and 3000 GeV curves, in the region κ 3 ∈ [2, 6], apart from the minimum around κ 3 = 1. This deformation does not appear (at least so clearly) in the curve √ s = 500 GeV, which only has one minimum around κ 3 = 1. Finally, Fig. 25 shows the dependence of σ(e + e − → HHHν eνe ) with κ 4 , setting κ 3 = 1, for various fixed energies. As we already commented, the deviations with respect to the SM are softer in this case. In contrast with the previous plot, in this one the curves exhibit only one minimum, which is around κ 4 = 1. The only curve in which this does not occur is the one corresponding to √ s = 500 GeV. The cross section at this energy is with difference the least sensible to variations in the κ 4 parameter. Now that we have seen the consequences of varying each of the κ 3 and κ 4 parameters separately, we study next the combined effect in triple Higgs production of varying both parameters at the same time. For reasons that we will motivate later, we will restrict ourselves in this study to the particular case of √ s = 3000 GeV. The results of the contour lines for σ(e + e − → HHHν eνe ) in the (κ 3 , κ 4 ) plane are shown in the left plot in Fig. 26. These results are consistent with the observations made in the previous plots, and confirm clearly that the deviations in the cross section are much stronger in the κ 3 direction than in the κ 4 one. In addition, the sensitivity to κ 4 is larger at the central region of the explored κ 3 interval, and the largest rates in this plot are obtained in the upper left corner, i.e. for the extreme values of (κ 3 , κ 4 ) = (−10, 10). Furthermore, if we look back to the subprocess plots in Fig. 17, we notice that, qualitatively, the variations in the cross section of the e + e − process at √ s = 3000 GeV behave similarly to those of the W − W + → HHH subprocess in the region around  √ŝ = 1000 GeV. This suggests that the 'effective energy' for WWS in e + e − collisions at √ s = 3000 GeV is approximately around √ŝ ∼ 1000 GeV. Finally, we show in the right plot in Fig. 26 the corresponding predictions for the expected number of events in the (κ 3 , κ 4 ) plane in the most favourable case, namely, at the last stage of CLIC with √ s = 3000 GeV and L = 5 ab −1 . As it can be seen in this plot, the expected number of HHHνν events from BSM Higgs couplings can be sizeable in a large region of the (κ 3 , κ 4 ) plane. Indeed, it is much higher than in the SM case, which as we have seen produces negligible event rates. This is true even if we stay within the experimental limits, κ 3 ∈ [−2.3, 10.3] at a 95% CL (remember again that there are no present constraints on κ 4 ). Now that we have characterized both the W − W + → HH(H) subprocesses and the whole e + e − → HH(H)νν processes, we are close to being able to explore the final sensitivity to the EChL parameters. However, we still need to study these processes in the framework of a collider, studying the real final state particles, namely, after the Higgs bosons decays. In the following section we will perform such an analysis, motivating why 3 TeV is the optimal energy to obtain the best BSM signals in double and especially triple Higgs production.
6 Sensitivity to BSM couplings in multiple b-jet events We will focus our forthcoming analysis in the two future linear colliders that are currently under study, the ILC and CLIC (see Tab. 1). As we already mentioned in the introduction, they are both e + e − colliders, and will operate at energies between a few hundreds of GeV and 3 TeV. Each of these energy stages serves a different purpose, being the higher energy configurations the ones oriented to measuring the SM triple Higgs self-coupling via HH production. In principle, none of them is expected to yield measurable signals of HHH production, assuming SM rates. Therefore, our studies of BSM signals via HHH production are very singular in the sense that they will not have the SM as a competitor, since it produces negligible rates (less than 1 event in all cases). On the other hand, to obtain testable results in these e + e − colliders via the processes of our interest, e + e − → HH(H)νν, it is necessary to analyze final states where the Higgs bosons have decayed. Here we will choose the Higgs main decay channel, H → bb, yielding final states with multiple b-jets, arising from the final b quarks, plus missing energy, associated to the final neutrinos. Thus, we explore the sensitivity to the EChL parameters is this type of multiple b-jet events: 1) a and b in events with missing transverse energy and four b-jets from the decays of the two final H's; and 2) κ 3 and κ 4 in events with missing transverse energy and six b-jets from the decays of the three final H's. For the present study we ignore potential backgrounds to multiple b-jets production accompanied by missing energy, which a priori are expected to be negligible in this e + e − context. A more refined analysis, including realistic backgrounds and considering the peculiarities of the planned detectors, is clearly beyond the scope of the present work and is left for future studies.

Sensitivity to a and b in events with 4 b-jets and missing transverse energy
We explore here the sensitivity to a and b in e + e − → HHνν → bbbbνν events. For this study we perform a Monte Carlo simulation employing MG5: the same computation that was already worked out, but including now the decays of the Higgs bosons. Regarding the particular collider project, we consider here three of them: 1) The ILC with √ s = 500 GeV and L = 4 ab −1 ; 2) The ILC with √ s = 1 TeV and L = 8 ab −1 ; and 3) The CLIC with √ s = 3 TeV and L = 5 ab −1 .
In order to characterise the BSM signal arising from the EChL anomalous couplings (a, b), we have first to define the final state with particles that can be detected in the experiments.
The events that we study here contain four b-jets from the hadronization of the final b quarks, which are produced in the Higgs decays, and missing energy corresponding to the final ν and ν escaping detection. This characteristic missing energy from the neutrinos and antineutrinos could allow to differentiate the BSM signal from potential SM backgrounds, like those of QCD origin. These channels exhibit multiple jets in the final state due to the hadronization of quarks and gluons, but typically show no relevant missing energy. Therefore, to explore the BSM signal we will implement some cuts on the relevant variables of the final particles.
In particular, the b-jets will be required to have a minimum transverse momentum in order to be detected, and also the missing transverse energy will be required to be above a minimum value. Besides, a cut in the polar angle θ, or equivalently, a maximum pseudorapidity η j ≡ − log tan θ 2 for the jets, is required. Finally, in order to be detectable, the jets need to exhibit a certain angular separation among them. This is equivalent to stablishing a minimum value for the variable ∆R jj ≡ (∆η jj ) 2 + (∆φ jj ) 2 , where ∆η jj and ∆φ jj are the separations in pseudorapidity and azimuthal angle of the two jets jj , respectively. The values chosen for all these cuts in this work, which are similar to those taken in references [38] and [51], are summarised as follows: The decay of the two Higgs bosons will lead to a reduction factor in the event rates of 0.58 2 , due to the branching ratios (BRs) of the decays to b quarks. Once the cuts on the final state are implemented, the signal event rates will obviously suffer a further reduction. To estimate the effects of these cuts in Eq. (13) we have evaluated the acceptance, A, given by the ratio of the predicted bbbbνν rates after applying these cuts divided by the rates before cuts (the latter are the ones already shown in Fig. 19). We have evaluated A in the (a, b) plane, considering the region delimited by a and b in the interval [0.5,1.5], as in Fig. 19, and we have done it for the three energies chosen here, of 500 GeV, 1000 GeV and 3000 GeV. We have found that A diminishes when increasing the collider energy. For 500 GeV it varies approximately in the range 0.65-0.72, for 1 TeV in the range 0.55-0.65, and for 3 TeV in the range 0.3-0.5. There is a qualitative difference between the case of 500 GeV and the other two collider stages, with 1 TeV and 3 TeV respectively. The acceptance for 500 GeV is maximal, above 0.7, at the region close to the SM point (a, b) = (1, 1), whereas in the 1 TeV and 3 TeV cases it is minimal in this region, with values around 0.55 and 0.3 respectively. This qualitative difference can be understood from the fact that the kinematical configuration of the final state particles in two cases of 1 and 3 TeV is driven from the dominant WWS subprocess, which is not the case for the 500 GeV collider, where the Z-mediated diagrams dominate.
In order to conclude on the potential accesibility to these a and b parameters at future e + e − colliders, we also have to take into account the efficiency in the detection of the four b-jets. Thus, the computed rates must be reduced by an extra factor of 4 , where is the b-tagging efficiency factor. We will assume here of an 80%, which is the value commonly used in the literature, see for instance [37] and [38]. Taking into account all the cuts shown above, the BRs of the decays, and the efficiency reduction factors, we have finally computed the bbbbνν event rates of our BSM signal for each choice of the anomalous couplings a and b and for each selected collider setup. We display, in Fig. 27  predicted number of events in the (a, b) plane for: ILC(500 GeV) with 4 ab −1 (left plot), ILC(1TeV) with 8 ab −1 (middle plot), and CLIC(3TeV) with 5 ab −1 (right plot). Notice that we have adjusted the interval displayed for the paramater a so as to coincide with its presently experimentall allowed interval, whereas for the b parameter we have displayed the full interval [0.5,1.5], as in the previous figures. From this Fig. 27 we can already extract some conclusions regarding the BSM event rates compared to the SM ones. In the three plots shown we find areas (depicted in the lighter colors), where the ratio BSM/SM is considerably larger than 1. Regarding the relative events statistics, the lowest collider energy provides the smallest rates, and the highest collider energy the largest ones, as expected. For ILC(500 GeV) with 4ab −1 we only find a few events, varying from less than 10 in the lower-left part of the plot to around 22 in the lower-right corner. These should be compared with the 9 events predicted in the SM. For ILC(1TeV) with 8ab −1 we find larger rates, around 50-70 events in the region close and around b = 1 with a < 1, and up to more than 750 events in the lower-right corner. The SM predicts 100 events in this case. For CLIC(3TeV) with 5ab −1 the largest rates are found, ranging from around 100 events in the region very close and around b = 1 with a < 1 to more than 8000 events in the lower-right corner. Therefore, CLIC offers better rates than the other two explored projects and seems to be the best option in order to be sensitive to the (a, b) parameters. Finally, to provide a more concrete conclusion on the accessibilty to these two anomalous couplings, we have evaluated the ratio R ≡ (N BSM −N SM )/ √ N SM , where N (B)SM is the number of events in the (B)SM case. This quantity is a way to measure the size of the deviation of the BSM signal with respect to the SM prediction. Since there are very low statistics in the 500 GeV case, both for BSM and the SM, we focus on the other two options. We show in and as R > 10 (purple region bounded by dashed contourline) in the more conservative case. It is clear from this plot that both options, ILC(1TeV) and CLIC(3TeV), will offer a good accessibility to measure a and b beyond their present experimental constraints from LHC. The best option will be clearly CLIC (3TeV) with 5ab −1 , where the unaccessible area (in white) shrinks around the SM point, especially in the b direction, which is the worst explored at present. The sizes of these unaccesible areas provide an approximate estimate of the expected improvements on the constraints on these anomalous couplings. Of course, in order to determine more accurately the sensitivity to variations in a and b with respect to their SM values, it would be necessary to analyze the possible backgrounds as well. Nevertheless, this is beyond the intention of this work and is left for future research.
6.2 Sensitivity to κ 3 and κ 4 in events with 6 b-jets and missing transverse energy In this section we study the sensitivity to the anomalous couplings κ 3 and κ 4 in e + e − → HHHνν → bbbbbbνν events. Therefore, we consider the dominant decays of the final Higgs bosons leading to b-jets and take into account the missing energy left by the final neutrinos and antineutrinos. We perform our analysis in the most favourable scenario, which is the last stage of CLIC, at a CM energy of 3000 GeV with an integrated luminosity of 5 ab −1 . In the previous sections we learnt that it is at that high energies and luminosities where the channel with three Higgs bosons and neutrinos can be most sensitive to variations of the Higgs self-couplings. As we have said, it is just in triple, and not double, Higgs production where there is a unique sensitivity to the quartic Higgs coupling κ 4 . Indeed, considering at the same time deviations in κ 3 and κ 4 in triple Higgs production plus neutrinos can boost the number of events from not more than one in the SM (which does not yield a detectable signal) to tens, hundreds or even thousands in the most extreme cases of BSM scenarios. Thus, in this section we will focus on characterising these BSM signals from the triple Higgs channel, with mainly 6 b-jets and missing transverse energy, at the last stage of CLIC. First, we investigate the main features of the kinematical configuration in these multiparticle events for our BSM signal with (κ 3 , κ 4 ) = (1, 1) and compare it with the SM case with (κ 3 , κ 4 ) = (1, 1) . For this characterization, we will generate a set of samples for different values of κ 3 and κ 4 and study some interesting event distributions. All the events will be generated using MG5 and analyzed using ROOT 6 [52]. Since the hadronization of the b quarks in the final state is a very expensive task, in order to define the b-jets we will use a resolution criterion instead. Following the reasoning in [38], we will consider an energy resolution of ∆E/E = 5% and assume that two quarks with a small separation of ∆R qq < 0.4 cannot be resolved individually. This condition will be applied recursively until we converge to a final list of quarks that we will identify as the b-jets. The plots in Figs. 29 to 34 show the distributions with respect to several kinematic variables for different values of (κ 3 , κ 4 ). Concretely, we chose: 1) the invariant mass of the three Higgs bosons M HHH , equivalent to that of the six b-jets (Fig. 29); 2) the missing transverse energy E / T , equivalent to that of the final νν pair (Fig. 30); 3) the realistic number of b-jets in the events, N b−jet , due to the resolution criterion commented above (Fig. 31); 4) the angular separation between two b-jets, ∆R bb (Fig. 32); 5) the transverse momentum of the b-jets, p b T (Fig. 33), and 6) the pseudorapidity of the b-jets, η b (Fig. 34). For simplicity we have only plotted deviations of either κ 3 or κ 4 , with the other one set to 1. Notice also that all the histograms are normalized to unity (and not to the real number of events) since we are interested now in comparing just their shape. We will next comment some general features learnt from these plots. First, we note that the invariant mass distribution of the three Higgs bosons, which is equivalent to the energy of the subprocess, √ŝ , confirms the observation we made when comparing Fig. 26 with the corresponding plots for the W − W + → HHH subprocess: we suggested that the 'effective energy' for WWS in e + e − collisions at √ s = 3 TeV is approximately around √ŝ = 1 TeV. In fact, in the distribution for the SM case we see in Fig. 29 a maximum centered around M HHH = 1000 GeV, while for other values of κ 3 and κ 4 the maximum is generally displaced to lower energies, and the peaks exhibit different shapes. From the E / T and p b T plots, we also see that both the neutrino-antineutrino pairs and the b-jets tend to be produced with higher transverse momentum in BSM scenarios than in the SM, which is consistent with a smaller value of the pseudorapidity, as we see in the plot of η b . As a final remark, we also notice that the realistic number of b-jets appears to decrease as we deviate from the SM. This is due to the b quarks being produced with smaller relative angles, as can be seen in the distribution with respect to the variable ∆R bb . The area of this distribution in Fig. 32 with larger rates is displaced to lower values of ∆R bb in the BSM cases than in the SM, meaning that the b-jets may not be always identified individually, thus yielding in some cases an apparently lower number of produced b-jets. Finally, for the jet analysis, we need to take also into account the b-jet identification efficiency. As in the previous section, we will adopt here a b-tagging efficiency of 80% [38]. This corresponds to a misidentification efficiency of 10% for c-jets and 1% for light flavour jets.    Overall, considering the features learnt from all these distributions and the needed requirement of not loosing too much signal, we will impose the following set of cuts, which are similar to those in [38], but also requiring a minimum value of the missing transverse energy: Some comments are in order. First, regarding the cut in pseudorapidity we use here the value reported in [53] of |η| max = 2.72 for the particular case of CLIC. Notice that this is slightly less restrictive than the one chosen in the previous section of |η j | < 2. Since the rates for triple Higgs production are much lower than for double Higgs production, we relax this cut to avoid the loss of too many signal events. Second, we require our events to include at least six jets. We also impose that these jets have a minimum transverse momentum of 20 GeV, since very low p T jets can be difficult to detect. Next, we want at least five jets to be identified as b-jets. The reason why we do not require 6 tagged b-jets lies on the efficiency. Assuming that all jets are equal, i.e., that they are not sorted or classified in any way, the probability of identifying five out of the six as b-jets is: while if we tag the six of them: so allowing five tagged jets is way more efficient. It is also important to note that since all b-jets come from on-shell Higgs bosons, the three pairs should reconstruct the Higgs invariant mass. Although the resolution will not be optimal since we are treating with jets, this could be an additional cut to reject backgrounds. We have computed the value of the acceptance, A, for the BSM events after applying all these cuts, and found that for the studied values of κ 3 and κ 4 it is between 0.44 and 0.51. In the SM case, with κ 3 = κ 4 = 1, the acceptance is lower and drops to 0.34. Thus, in summary, to finally compute the number of predicted BSM events, we use the following reduction factors, considering the BRs to b quarks, the b-tagging efficiencies, and the acceptance, A, after applying all the cuts in Eq. (14):  As it can be clearly seen in this Fig. 35, the bbbbbb νν event rates are very low in the region of the (κ 3 , κ 4 ) plane close to the SM point, leading to less than one event, therefore yielding unobservable signals. Separating from this area, the BSM event rates increase reaching values above 10, 100 and even 1000 in the extreme cases with κ 3 near -10. We believe that having this signal statistics as a starting point is very promising and motivates a more complete study including backgrounds, which we are neglecting here. As we already stated, a detailed analysis taking into account all the backgrounds and the characteristics of the particular detectors at CLIC is beyond the scope of this work.
Finally, we wish to conclude with our estimate of the accessible region to the BSM anomalous couplings in the (κ 3 , κ 4 ) plane. This is summarised in Fig. 36, where we have set the 'naive' criterion for accessibilily to a given (κ 3 , κ 4 ) (in absence of background) by requiring a minimum value of 10 events. Thus, we find that the purple area in this Fig. 36 leads to more than 10 events (reaching large values of hundreds and even thousands at the extreme values) and will give access to values of the anomalous self-couplings which are at present unconstrained. In particular, smaller values of |κ 3 | inside the presently allowed interval κ 3 ∈ [−2.3, 10.3] (marked by the red arrows in this plot) could be tested at CLIC. Regarding the κ 4 parameter, it is clear from Fig. 36 that this HHHνν channel opens the possibility to test the yet unexplored BSM quartic Higgs self-coupling at CLIC. Concretely, if we assume BSM scenarios that are outside the unreachable white region at the center of this plot (which contains the SM point), there seems to be access to κ 4 ∈ [−10, 10]. Figure 36: Accessible region in the (κ 3 , κ 4 ) plane to BSM scenarios via e + e − → HHHνν → bbbbbbνν events at CLIC, with √ s = 3000 GeV and L int = 5 ab −1 . The criterion of accesibility assumed here is requiring more than 10 signal events (purple area). The red arrows mark the limits of the present bound from ATLAS [3] for κ 3 , given in Eq. (5).

Conclusions
In this work we have explored the sensitivity at future e + e − colliders to BSM Higgs physics induced by the W W H, W W HH, HHH and HHHH interaction vertices within the context of the non-linear effective field theory given by the EChL. These interactions provide the BSM Higgs anomalous couplings of our interest, which depend respectively on the EChL parameters a, b, κ 3 and κ 4 . In this non-linear EFT context, these coefficients are independent and uncorrelated by any symmetry argument, because the Higgs field is a singlet in the EChL. This is in contrast to the Higgs couplings in the SM, where some relations among the interaction vertices, such as V SM W W H = vV SM W W HH and V SM HHH = vV SM HHHH , are derived from the Higgs being a component of a doublet.
We have shown that the double and triple Higgs production channels, in particular e + e − → HHνν and e + e − → HHHνν, provide the proper window to explore efficiently these BSM Higgs couplings. Whereas double Higgs production can access to a, b and κ 3 , the triple Higgs channel is the only one that can access to κ 4 , a parameter which is at present totally unconstrained. Our study here then provides a first try to test the BSM quartic Higgs self-coupling, which is so far evading all experimental searches, at e + e − colliders.
We have also understood why these two particular channels are so efficient in the searches for BSM signals arising from the anomalous Higgs couplings. The main reason is the fact that, at the TeV scale and above, they are dominated by WWS mediated diagrams. It is in these configurations where the largest enhancement due to BSM physics occurs. We also conclude that the main features found in the WWS subprocesses, W − W + → HH(H), are basically reproduced in the corresponding e + e − → HH(H)νν processes. In our study of the dominance of the WWS subprocesses, we also showed that the improved EWA works remarkably well for the double Higgs production case, but it fails in the predictions for triple Higgs production. This is a new result, and convenient to be aware of for future BSM searches. The conclusion in this concern is that the EWA can be safely used for HH but it should not be used for the HHH case.
Our final study of a more realistic experimental scenario, with multiple b-jets and missing transverse energy in the final state, shows that the sensitivity to these parameters will indeed be improved considerably at the future e + e − colliders. The final results in Fig. 27 and Fig. 36 summarise our main conclusion on the accessibility region for these EChL parameters, a and b in Fig. 27, and κ 3 and κ 4 in Fig. 36. These plots also suggest that the best tests of these parameters will be provided by CLIC. In particular, triple Higgs production seems to provide access to κ 4 in the latest stage of this collider, with a luminosity of 5 fb −1 . Of course, in order to obtain a more accurate conclusion, a more complete study should be performed, including realistic backgrounds and taking into account the detector properties.
B Diagrams contributing to W − W + → HHH There are 25 diagrams contributing to the scattering amplitude of the W − W + → HHH subprocess in the unitary gauge. These have been generated using FeynArts-3.10 and are displayed in Fig. 38. The red and blue dots represent the triple and quartic Higgs selfinteractions, respectively. We will omit the analytical result of the corresponding amplitudes for shortness.