Production of two, three, and four Higgs bosons: where SMEFT and HEFT depart

: In this article we study the phenomenological implications of multiple Higgs boson production from longitudinal vector boson scattering in the context of effective field theories. We find compact representations for effective tree-level amplitudes with up to four final state Higgs bosons. Total cross sections are then computed for scenarios relevant at the LHC in which we find the general Higgs Effective Theory (HEFT) prediction avoids the heavy suppression observed in Standard Model Effective Field Theory (SMEFT).


Introduction
The Higgs boson was first observed in 2012 at the Large Hadron Collider (LHC) [1][2][3], confirming the existence of a field that endows some of the other elementary particles with mass, playing a fundamental role in our understanding of the universe.However, many questions remain open regarding its origin and its interplay with the rest of the Standard Model (SM) fields.
In recent years, effective field theories (EFTs) have proven to be a useful tool to parametrize the low energy effects of the underlying, and yet unknown, ultraviolet physics.Among these EFTs, the Standard Model Effective Field Theory (SMEFT) and Higgs Effective Field Theory (HEFT) have emerged as powerful frameworks for describing particle interactions at various energy scales.In this paper, we explore the theoretical foundations and predictive capabilities of SMEFT and HEFT, with a particular focus on their applicability to the observed phenomena at the LHC.By comparing the predictions of these theories with experimental data obtained from the LHC, we aim to gain a deeper understanding of their strengths and limitations in describing the behaviour of subatomic particles.
When trying to compare EFT predictions to experimental observations, two approaches are the standard: either through a global fit of cross sections and differential distributions or through the definition of realistic (pseudo-)observables, such as fiducial cross sections, ratios, and asymmetries.Several groups have applied these techniques to the available data from LHC Runs I and II, for example, see [4][5][6][7][8][9][10][11][12][13][14][15].In the spirit of the latter approach, we propose a systematic study of longitudinal vector boson scattering cross sections.
Here, we focus on the electroweak sector of the Standard Model, allowing us to explore further the connection between the Higgs and gauge bosons.In this sector, by only looking at a handful of EFT operators (initially only one or two), we can already extract a great deal of information on the Higgs-Goldstone interactions and hence, on the electroweak symmetry breaking (EWSB) mechanism.
At the energies currently reached by accelerators, both HEFT and SMEFT seem to fit the available observations.Recently, several works have addressed the distinction between the two EFTs, see for example [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20].Moreover, in previous works [21][22][23], we showed that, from the SMEFT point of view, there is essentially one operator of the Warsaw basis that can, alone, shed light on the open questions of the symmetry breaking mechanism (this operator coincidentally being one of the least constrained by current LHC fits [24][25][26][27] due to the narrow width approximation used by the experiments).Furthermore, we showed that by comparing different scattering amplitudes of longitudinal gauge bosons one can draw a clear distinction between a SMEFT-like and a HEFT-like scenario.While current bounds on the ωωh and ωωhh interaction vertices are still loose, we found that this method could be useful in the future to eventually falsify the SMEFT approach (and hence the traditional EWSB realization à la SM).
In this work, we go one step further by calculating analytical and numerical predictions for different cross sections.We focus on the electroweak production of 2, 3, and 4 Higgs bosons, which we study in various benchmark scenarios.In particular, we study the longitudinal vector boson scattering (VBS) into multiple Higgs bosons, which has been subject of study of several recent works, for example [18,[28][29][30][31][32][33][34].
For the comparison of the SMEFT scenarios and the HEFT one, it will be fundamental to introduce the flare function, F(h).This function takes the role of parametrizing the Higgs-gauge interaction and it has the advantage that can be defined for both the SMEFT and HEFT theories 1 , allowing us to map individual measurements to both scenarios.The flare function will be discussed in detail below (see also ref. [21]) and provides the effective coupling a n of a pair of EW gauge bosons, W W , to an arbitrary number n of Higgs bosons: It is customary that the first coefficients are also denoted by a ≡ a 1 /2 and b ≡ a 2 .In order to compute the production of Higgs bosons from the scattering of longitudinally polarized weak bosons, W L W L → n × h, we will employ the equivalence theorem (EqTh) [35][36][37][38][39][40][41], where the EW gauge boson scattering W L W L → n × h can be approximated by the EW Goldstone scattering ωω → n × h for s ≫ m 2 W . Furthermore, both the EW gauge boson and Higgs boson masses will be neglected with respect to the center of mass energy of the process √ s.Hence, we emphasize that the present article will be fully focused in the computation of these scattering amplitudes, ωω → n × h, in the massless limit.More precisely, we will explore the production of two, three and four Higgs bosons, which we consider that are representative enough of the issue at hand.Thus, this function F(h) is going to provide the irreducible tree-level vertices for ωω → n × h.However, in general, several couplings a j are going to contribute to the tree-level amplitude ωω → n × h.All the amplitude and cross section calculations in this article will use this standard HEFT Lagrangian [42][43][44][45][46][47][48][49] and the associated flare function F(h) [21][22][23].
We display compact analytical expressions for the processes at hand.These are fundamental to understand patterns and simplifications that happen in the New Physics scenarios.For example, we find a factorization of the s-dependence in our analytical VBS expressions, such that only two numerical integrations are required to evaluate the cross sections at any energy, vastly reducing the computational cost.When looking at even higher Higgs multiplicities and full collider simulations this feature may be crucial.
Still, a superficial observation of these combinations of a j effective couplings may fail to reveal a pattern in the amplitudes.In fact, it is possible to eliminate the hωω vertex by means of appropriate field redefinitions (see appendix C).Although this distorts the effective Lagrangian and makes its symmetries less evident, this simplification makes the massless ωω → n × h computations much more transparent: the number of tree-level topologies is greatly reduced, with the diagrams now provided by an effective flare function with the form F(h) = 1 + â2 5 ), where â2 = a 2 − a 2  1 /4, â3 = a 3 − 2 3 a 1 a 2 − a 2 1 /4 and â4 = a 4 − 3 4 a 1 a 3 + 5 12 a 2 1 a 2 − a 2 1 /4 .Further effective combinations are of no interest for this work although they can be extracted in a similar way.We have employed these simplified interactions to double-check all the computations.The technical details of this additional test have been relegated to appendix C.
By studying longitudinal VBS into multiple Higgs bosons, this article explores the delicate cancellations one finds in these scattering amplitudes in the case when the lowenergy range can be described by SMEFT.These cancellations are absent in more general theories that do not accept a low-energy EFT in terms of a Higgs doublet, leading to cross sections order of magnitude larger than those one obtain in SMEFT analyses, as we will see later in the main text.
In section 2, we compute the scattering amplitudes and cross sections for ωω into two, three and four Higgs bosons.Section 3 particularizes these outcomes to the case when the underlying theory accepts a SMEFT description at low energies.Section 4 contains a phenomenological study of the obtained effective theory amplitudes.We use some benchmarks to compare both effective approaches, HEFT and SMEFT, showing how SMEFT models in general provide much more suppressed cross sections than theories which does not accept a SMEFT description.We show in section 5 how stringent this suppression must be, providing cross section exclusion plots for scenarios that admit SMEFT.Our final conclusions are provided in section 6.Some technical aspects have been relegated to the appendices.All the results in this article have passed several checks: we have used an specific Mathematica code for the generation of arbitrary ωω → n × h amplitudes [50]; we have used a HEFT model file [51] in FeynRules [52] and FeynCalc [53] for the analytical calculations of the VBS amplitudes; the massless phase-space integration code MaMuPaXS [54] has been developed for the numerical evaluation of the ωω → n × h cross sections in the approximations of this article; finally, we have used the HEFT model file [55] in combination with MadGraph5 aMC [56] for a numerical test of the analytical expressions.

Calculation in HEFT
In [21] we presented scattering amplitudes for different multiplicities of Higgs bosons in the final state.Here, we will focus on the cancellations occurring for SMEFT-like theories, in contrast to a generic HEFT scenario.We will be using the derivative terms of the leading order (LO) HEFT Lagrangian [42][43][44][45][46][47][48][49], where the unitary matrix L+R coset in terms of the EW Goldstone bosons ω a .A usual representation is provided, e.g., by the exponential form U = exp{iω a σ a /v}.Nevertheless, for the tree-level ωω → n × h scattering discussed in this article, one only needs the O(ω) term in U , which is identical in all representations.Finally, we will neglect next-to-leading order, O(p 4 ), HEFT contributions [49,[57][58][59][60][61][62][63], both at tree and one-loop level.For this work, eq.(2.1) provides the relevant part of the HEFT Lagrangian because we are working under the following set of restrictions: First, we are only considering the derivative part of the lowest order chiral Lagrangian, O(∂ 2 ), as we will be choosing a kinematic regime well over the production threshold (s ≫ m 2 W ∼ m 2 h ); Second, eq.(2.1) incorporates only the scalar sector, as we will be using the equivalence theorem to compute the amplitudes W L W L → n × h at tree-level for s ≫ m 2 W [35][36][37][38][39][40][41].In this article we will focus on the charged Goldstone (ω ± ) scattering, although a similar analysis could be done with the neutral ones (ω 0 ).This exercise can be repeated similarly for the Yukawa sector through the study of t t → n × h processes, where one would also need to introduce the Yukawa flare function, which parametrizes the relevant local vertices (see [12,22]).
Thus, we will compute the ω + ω − → n × h amplitudes, T ωω→n×h , and analyse the corresponding cross section, It is useful to note that, in the massless approximation, one can fully factor s out of the phase-space integration dΠ n = s n−2 d Π n .In this way we are left with a pure angular integration in d Π n .Further details like, e.g., the analytical expression for the total phasespace volume V n = dΠ n , can be found in appendix A.
In the following subsections we describe the amplitudes and cross sections for the most general case, that of HEFT.In section 3 we will be restricting them to the SMEFT case.
The amplitude in eq.(2.3) is a pure J = 0 wave, depending only on the total centerof-mass (CM) energy, not on the scattering angle.Thus, one can readily compute the total cross section, (2.4)

ωω → 3h
For the ωω → 3h amplitude, we have checked our result with the aid of FeynRules and FeynCalc (based on the previous work in [51,55]) and an independent Mathematica implementation by the authors [50], obtaining, This amplitude readily leads us to the corresponding total cross section: (2.6) Previous works have analyzed the W W → 3h scattering in the context of effective theories [69,70], considering modifications in the W W h, W W h 2 , h 3 and h 4 vertex couplings.One of the novelties of the present work is to incorporate the full generality of allowed ωωh n effective vertices relevant for the process (and the general W W h n vertices beyond the EqTh).In particular, we considered the ωωh 3 coupling (a 3 ).Nonetheless, we find that at high energies its contribution to the cross section yields the same energy and angular structure as those from a 1 and a 2 , being the relevant information encapsulated in the combination â3 = a 3 − 2 3 a 1 a 2 − a 2 1 /4 .

ωω → 4h
In this subsection we compute the amplitude of ω in HEFT in the same way as in the previous cases.This amplitude can be written as ) The dimensionless kinematic function B can be written as with where ] and the total fourmomentum q = k 1 + k 2 = p 1 + p 2 + p 3 + p 4 .In the centre-of-mass (CM) rest-frame these relations define: the three-momentum fractions for each outgoing Higgs boson; the angular functions z i = 2 sin 2 (θ i /2) with θ i being the angle between the i-th Higgs boson and the incoming ω + Goldstone boson momenta, ⃗ k 1 (that is, z 1 = 1 − cos θ, z 2 = 1 + cos θ as usual in a two-body problem with t and u channels); z ij = 2 sin 2 (θ ij /2), with θ ij being the angle between the i-th and j-th Higgs bosons.Total 4-momentum conservation allows us to establish some relations between these parameters which are used to simplify the expression for the amplitude.More details are given in appendix B, where we also present B in terms of f 1,2,3 , z 1,2,3 , z 12 and z 13 , reducing the number of free kinematic variables.
The term with B provides the irreducible contribution from cross-propagator exchanges (see appendix C).For the production of 2 and 3 Higgs bosons, even though we had crosschannel exchanges, the propagator poles always cancelled with the vertices when the external legs were set on-shell.This is not fully possible in the 4 Higgs case and the B term remains.
Figure 1 shows the tree-level diagrams contributing to the ωω → 4h process.From the structure of the various possible tree-level diagrams, containing crossed-channel ω exchanges, one might think that the amplitude (2.7) may become singular for some points of . Relevant tree-level topologies for ωω → 4h.Note that in addition one needs to consider all possible permutations for the assignment of the external particles.
the phase-space integration.Indeed, there are diagrams that contain up to three intermediate ω propagators.Nonetheless, a careful observation shows that the momentum-dependent vertex structure exactly cancels the propagators in several cases.
Diagram a) has no propagator.Likewise, when the external legs are set on-shell, the intermediate propagators in diagram b) and c) cancel out -diagram by diagram-and their contribution to the amplitude is simply a polynomial in momenta.However, for this process we find that there are topologies where one of the intermediate propagators survive, d), e), and f ), generating the B denominators which one finds in (2.7) and (2.8).
One might think that the surviving propagator in diagrams d), e) and f ) yields poles in the amplitude at certain energies.However, these potential singularities are always proportional to the structure )/T , with P A and P B the four-momenta of two disjoint pairs of final Higgs bosons (e.g., P A = p 1 + p 2 and A,B and C a constant proportional to v −4 .In the physical region where √ s ≥ M A + M B one has T < 0 when both M A > 0 and M B > 0. Thus, the amplitude never turns singular in the phase-space integral as the singularity is located at non-physical energies 2 .Furthermore, when either M A or M B are zero, the diagram becomes non-singular in the whole complex plane.In summary, as it happened with the two and three Higgs final states, the ωω → 4h cross section is infrared safe and does not contain pole singularities.
This amplitude provides the total cross section for ωω → 4h: , where Tmin ≤ T ≤ Tmax.For MA > 0 and MB > 0, ∆s ≥ 2MAMB > 0 in the physical region, leading to strictly negative values of T max (min) and T .
with the normalized phase-space integrals, and the total 4-particle phase-space volume V 4 = dΠ 4 = s 2 24(4π) 5 −1 , rendering the χ n dimensionless and s-independent (see appendix A for further details).These two constants are dimensionless and independent of the effective Lagrangian parameters.We evaluate them numerically through our phase-space integration code (MaMuPaXS) [54], obtaining: 10) , χ 2 = 0.0193760 ( 16) . (2.12) The angular function B is given in eq.(2.8) and its value is independent of the Higgs parameters.Notice that χ 1 and χ 2 are defined in the massless limit, therefore s-independent.Due to the normalization factor V −1 4 in (2.11), the χ j are pure angular integrations.The total phase-space volume is given by −1 (see appendix A for further details).Finally note that both type of diagrams, either with â4 or â2 2 , have a similar weight in the cross section; the â2 2 terms appear multiplied by (B − 1) in the amplitude and hence lead to the integrals V −1 4 dΠ 4 (B − 1) n ∼ O(1) in the cross section.For convenience, we have explicitly separated the B n integrals, which need to be computed numerically.
We have cross-checked these analytical results for ωω → 2h, ωω → 3h and ωω → 4h with MadGraph5 aMC, finding a complete agreement with its numerical determinations.In particular, we have also used the latter to check the numerical values for χ 1 and χ 2 obtained with MaMuPaXS [54].

SMEFT treatment: suppression of multi-Higgs production
In a previous work [21] we discussed that, neglecting custodial symmetry breaking terms, there is only one operator that modifies the flare function F(h) at each perturbative order in the SMEFT expansion.As it is also done for dimension D = 6 in the Warsaw basis [71], one can use equations of motion and field redefinitions to remove SMEFT operators proportional to □H (with H the Higgs doublet).We can choose an EFT operator basis that optimizes the present analysis so only one operator is relevant for the flare function at each D dimension.Hence, for theories that admit a SMEFT description one has, up to O(1/Λ 4 ), explicit expressions for the coefficients of F(h), in terms of the SMEFT Wilson coefficients and scale Λ: with, where c H□ is the dimensionless Wilson coefficient for the operator |H| D−4 □|H| 2 /Λ D−4 .Note that the expressions in (3.1) are truncated at O(d 2 ) = O(v 4 /Λ 4 ), and higher powers of d have been neglected.Higher coefficients a n with n ≥ 7 vanish at this order in SMEFT.
From (3.1) and section 2 results we can read the relevant combinations of a j couplings that determine the ωω scattering into two, three and four Higgs bosons in SMEFT: We note that these combinations are also found after considering the field redefinitions of appendix C. It is not difficult to observe that for d ̸ = 0 these expressions obey the O(d) ∼ O(1/Λ 2 ) SMEFT relations [22,23], with ∆a 1 ≡ a 1 − 2 and ∆a 2 ≡ a 2 − 1.These SMEFT relations have been refined and positively tested up to O(d 2 ) ∼ O(1/Λ 4 ).On the other hand, it is worth mentioning here that some UV-completions, such as the 2-Higgs Doublet Model, lead to a SMEFT low-energy theory but do not contribute to H□ and d vanish).Their first contribution to the flare function appears at O(1/Λ 4 ) through a non-zero c (8) H□ Wilson coefficient (this is, through a non-zero d 2 ρ).In this case, the SMEFT relations in [22,23] get modified [73] into the O(1/Λ 4 ) constraints, with ∆a 1 = 2∆a = 2c H□ v 4 /Λ 4 = d 2 ρ ̸ = 0.

ωω → 3h
For n = 3 Higgs bosons in the final state, we have the scattering amplitude and cross section, (3.9) We are providing the three-Higgs amplitude and cross section at the lowest non-trivial order in the SMEFT expansion, O(d 2 ) ∼ O(1/Λ 4 ) and O(d 4 ) ∼ O(1/Λ 8 ), respectively.

ωω → 4h
For n = 4 Higgs bosons in the final state, we have the scattering amplitude and cross section, with the angle-dependent function B from eq. (2.8) and the numerical evaluation in χ 1 and χ 2 provided in eq.(2.11) through our phase-space integration code MaMuPaXs [54].Again, we have computed the amplitude and cross section at their lowest non-trivial order in the SMEFT expansion, O(d . At lowest non-vanishing order in the SMEFT expansion, for ωω → n × h, we find the SMEFT amplitude: Up to O(1/Λ 2 ), the SMEFT Lagrangian for the scalar sector is given by the operators [21], where the dots stand for non-derivative operators, and other SMEFT operators of higher canonical dimension.The first term in the r.h.s.comes from the dimension D = 4 part of Lagrangian and the second one from L SMEFT (D=6) .At high enough energies over the production thresholds (m ≪ √ s) one may neglect the non-derivative terms and higher dimension operators will be subdominant for energies much lower than the effective theory cutoff ( √ s ≪ Λ).Although this article has been focused on custodial preserving operators, custodial breaking terms would show a similar structure at D = 6, with all the conclusions drawn in this subsection remaining unchanged if they were also incorporated.
For the reasoning in this subsection it will be better to handle the Higgs doublet in Cartesian coordinates, in such a way that 2|H| 2 = (v + h) 2 + a (ω a ) 2 .Also, notice that we are considering a linear representation of H, different to the non-linear one used in the previous sections.Abusing notation, we have denoted both linear and non-linear Goldstone and Higgs fields as h and ω a .Hence, the fields h and ω a in each representation should not be confused.Nonetheless, both representations of the Goldstone and Higgs fields lead to the same on-shell S-matrix elements.
With this consideration in mind, the first term on the r.h.s of eq. ( 3.13) provides the kinetic term of the Higgs h and EW Goldstones ω a .
The second term on the r.h.s. of eq.(3.13) yields three types of contribution: vertices with 4 scalars, vertices with 3 scalars and a correction to the kinetic term (which will need to be rescaled to provide the canonically normalized propagator).Note that all these operators contain two derivatives.
Feynman diagrams built by such Lagrangian will, in general, contain N V 3 vertices of 3 legs, N V 4 vertices of 4 legs, N I internal propagators, N E external lines and L loops.Here we will set the number of loops to zero and the amplitudes ωω → n × h will have where we eliminated N I by means of the relations 3N Tree-level diagrams will show a SMEFT 1/Λ 2 scaling of the form, ) at tree-level.Thus, for a given number of final Higgses n, we find that the leading diagrams scales like, Further clarifications on the suppression of the multi-Higgs amplitudes in SMEFT are provided in appendix D. Amplitudes with a larger number of Higgs bosons in the final state require a higher number of internal propagators N I .For this reason we can conclude that, in general, the production of multi-Higgs states is strongly suppressed in SMEFT, unlike in generic HEFT theories.The first non-trivial contributions to each process scale as: Note that higher SMEFT operators, with dimension D > 6, are never providing an amplitude that dominates over the diagrams coming from the leading SMEFT Lagrangian (3.13).Thus, the suppression (3.17) is a general feature of SMEFT.This result fully agrees with our findings in the framework of the non-linear HEFT realization if one imposes the SMEFT correlations in the flare function F(h) [21].A consistent SMEFT approach tells us that, e.g., the 1/Λ 4 amplitudes are incomplete unless we incorporate the L SMEFT (D=8) interactions to the analysis, as they also produce O(1/Λ 4 ) contributions to the scattering amplitude, as we have shown in previous subsections.
In contrast, non-SMEFT scenarios that allow the flare-function to deviate from those constraints can avoid this strong suppression in the multi-Higgs production cross sections.We will illustrate this with some numerical simulations in the next sections.

Cross section phenomenology
In this section, we perform a phenomenological analysis of the previous results.We study the multi-Higgs cross sections in a scenario where the SMEFT provides a proper lowenergy description and compare it with cases where it is not applicable, and only a HEFT description is appropriate.We will find an interesting suppression in SMEFT-compatible scenarios.
In general, we will show that the SMEFT cross sections are much smaller than those in the pure HEFT cases.This characteristic does not necessarily result from a large deviation from the SM couplings.Instead, we find that even slight deviations in the a j 's correlations, typical of the SMEFT, can lead to a significant enhancement of the cross section, even when the deviations from the SM couplings are small.
For this study we will consider a series of models and benchmark points (BP).In particular, for SMEFT we will be using, Here, d parametrizes the D = 6 correction and it is related to the hW W coupling a = 1 + d/2 at tree level.The chosen value d = 0.1 is then related to the SM deviation ∆a = a − 1 = 0.05, within the range of the most precise experimental determinations up to date from ATLAS [74], a = κ V = 1.035 ± 0.031, and CMS [75,76], H□ ) 2 ] parametrizes the relation between D = 8 and D = 6 SMEFT Wilson coefficients, which has been taken as ρ = 1 for all our BP.We will see that the precise values of these SMEFT parameters are not particularly relevant for the features observed in our HEFT-vs-SMEFT comparison.
We remark that, although the precise value of ρ is not really relevant for the qualitative behaviour of our results (as far as ρ ∼ O(1)), the choice of d is crucial for the rate of convergence of the SMEFT power series.We could have chosen a BP with a much smaller d (e.g., d = 0.01).However this choice would have made the following examples less illustrative: the large differences we show in this section between SMEFT and non-SMEFT models would be even larger by several orders of magnitude; the convergence of the SMEFT power series would improve and D = 8 corrections would be much more suppressed w.r.t. the D = 6 ones; and the scanning of the cross sections in terms of the a 3 and a 4 couplings we later perform would be more difficult to visualize.
The rationale behind the analysis presented here is as follows: we will consider benchmark points for different models and values of couplings which are very close to the SM and SMEFT in their first flare-function coefficients.The coupling a 1 = 2a is relatively well known from Higgs decay studies at the LHC.It is found to be close to the SM, up to O(10%) uncertainties.Current experimental determinations of a 2 = b from hh production at the LHC [75,[77][78][79][80][81][82][83][84] yield a very wide marginalized allowed range a 2 ∼ [ 0 , 2] (see appendix G for a detailed compilation of a 2 measurements).However, the combination of these results with the determinations for a 1 from Higgs decays will potentially reduce the uncertainty on a 2 to the O(10%) level.This can be observed in figure 2, where the hh production data is showing a fair enough determination of â2 = a 2 − a 2 1 /4, with |â 2 | ∼ O(10%) in the best cases [77] (this feature is easy to understand with the structure of our amplitudes and field redefinitions).Thus, even if the marginalized determinations of a 2 in those works carry a large dispersion, this is no longer true if one superimposes the experimental information on a 1 .
Recent phenomenological LHC studies [67] have shown how these banana shapes in the (a 1 , a 2 )-plane are expected to generally appear in future analyses, proving that we should eventually obtain a moderately good knowledge on â2 .Thus, potentially stringent experimental determinations of a 1 and â2 can serve to discern and discard new physics models, even if the marginal uncertainty on a 2 remains sizeable in the near future.For this reason, we will consider two types of BPs: one with values of a 1 within the experimentally allowed range while keeping a 2 free, and a second type in which both a 1 and a 2 are constrained to have small deviations from the Standard Model.
In the literature one can find several works that discuss specific UV-completions for the Higgs sector.Some of them accept a SMEFT description and some don't.However, to distinguish between them, it is not enough to know whether the model shows a non-linear HEFT structure, where the Higgs boson is a Goldstone from some BSM symmetry, or whether the EFT derives from some strongly-interacting underlying theory.
A paradigmatic example is that of the minimally composite Higgs Model [85,86], in which the EW Goldstones and the Higgs boson have a composite structure and evident non-linear transformation properties, and nontheless still be described by a SMEFT-like Lagrangian for a compositeness scale f ≫ v [7][8][9].
Another example is that of the dilaton model [87][88][89][90][91], where the Higgs is a (pseudo) Goldstone from the spontaneous conformal symmetry breaking.This model does not have a SMEFT description, since the theory does not obey the required flare function regularity constraints [21,92].Still, in this model, the ωω → n × h cross sections vanish in the EqTh approximation for any number of final Higgs particles (see appendix C for details), and becomes of little interest for our purposes.
Alternatively, there are UV-complete models with perturbative coupling expansions which give place to pure-HEFT low-energy theories, and cannot be described by SMEFT.Singlet-scalar extensions of the SM can lead to a pure-HEFT picture if the singlet's mass and its vacuum expectation value get a large enough separation [6,93,94].Also, the 2-Higgs Doublet Model can produce a non-SMEFT low-energy theory for certain choices of its parameters [95].These studies ensure the preservation of unitarity and perturbativity; nonetheless, the integration of the heavy degrees of freedom results in an EFT expansion that cannot be straightforwardly organized in terms of the canonical dimension of the operators.
The structure of these particular models is indeed quite involved and an exhaustive analysis of the various non-SMEFT cases therein is far beyond the scope of the present work.In this article, we will illustrate our phenomenological studies with the simplified flare function models considered in the next subsection.

Models and benchmark points
Here we present the analytical expressions for each a j in terms of the specified numerical inputs for every BP.To illustrate, we provide approximate numerical values for these a j couplings, although their precise values can always be readily determined using the BP analytical relations.We will examine the following cases: • SMEFT: for this BP we will be using the previously selected inputs d = 0.1 , ρ = 1 in (3.1).
• Non-SMEFT BP1: exponential flare function In this scenario, we provide a flare function, F(h), with no real zeroes (given that f (h) is a polynomial).It describes models which do not accept a SMEFT description [6,9,21,22], but fulfil the positivity requirements pointed out in [21].Furthermore, the general term a n in F(h) is going to have a behaviour of the form a n ∼ 1 n! , decreasing and vanishing at large n.
For a fair comparison with the SMEFT BP, we will match the first coefficients to those in the D = 6 expressions in (4.2).Nevertheless, we will see that it is not really crucial for the results to consider the D = 6 approximation as similar results are obtained if one considers either the D = 8 expressions (4.3) or the SM values a 2 = a 1 /2 = 1, a k≥3 = 0.
-BP1 (a 1 ) : We consider the simplest exponential flare function that matches the D = 6 SMEFT prediction in (4.2) for a 1 : where we are using the input a = a 1 2 = 1.05 from (4.2). -BP1 (a 1 , a 2 ) : In a second approach we consider the simplest exponential flare function that matches the D = 6 SMEFT prediction in (4.2) for a 1 and a 2 : where we are using the input a = a 1 2 = 1.05 and b = a 2 = 1.20 from (4.2).
• Non-SMEFT BP2: rational flare function In this scenario we provide, in terms of a polynomial f (h), a flare function which also has no real zeroes and, hence, the corresponding model does not accept a SMEFT description.By construction, this type of models also fulfill the positivity requirements [21].In this case, however, contrary to the exponential case, we will be dealing with general terms a n in F(h) that turn large as n grows, with a behaviour a n ∼ n.
Again, for an appropriate comparison with the SMEFT BP, we will match the first coefficients to those in the D = 6 expressions in (4.2).As it happened with the BP1 models, we will see that the precise value of a 1 and a 2 is not critical for the conclusions drawn here, as long as these couplings are close to their SM limit.
-BP2 (a 1 ) : We consider the simplest rational flare function that coincides with the D = 6 SMEFT prediction in (4.2) for the value of a 1 : where we are using the input a = a 1 2 = 1.05 from (4.2). -BP2 (a 1 , a 2 ) : In a second approach we consider the simplest rational flare function that matches the D = 6 SMEFT prediction in (4.2) for the values of a 1 and a 2 : where we are using the input a = a 1 2 = 1.05 and b = a 2 = 1.20 from (4.2).

BP study for ωω → 2h
Figure 3 shows the expected ωω → 2h cross sections for SMEFT.We are using the D = 6 coupling values (4.2).Dimension D = 8 corrections are further suppressed and lead to a similar outcome.We show the (D = 6)-vs-(D = 8) comparison in this plot for illustration but we will no further discuss the SMEFT (D=8) in what follows, only the D = 6 predictions.SMEFT is then compared in figure 3 with the two non-SMEFT BP's discussed in this article.We plot the predictions for BP1 (a 1 ) and BP2 (a 1 ) .The refined BP's, BP1 (a 1 ,a 2 ) and BP2 (a 1 ,a 2 ) , are not discussed in the plot.Since by construction a 1 and a 2 are set to the SMEFT (D=6) values in (4.2), they yield the same ωω → 2h cross section as SMEFT (D=6) .
Our conclusion is that, as expected from previous sections, the SMEFT cross section is suppressed by several orders of magnitude with respect to the non-SMEFT ones: in general, we will see that the fractional flare functions (BP2) tend to yield higher cross section than the exponential ones (BP1), both of them being well over the SMEFT determinations.[fb] BP1 (a 1 ) BP2 (a 1 ) Figure 3.Comparison of the ωω → 2h cross section predictions of SMEFT (D=6) , BP1 (a1) and BP2 (a1) .For illustration, we also provide the results for SMEFT (D=8) , which refines the previous order SMEFT determination.The comparison of non-SMEFT BP's with both SMEFT (D=6) and SMEFT (D=8) casts similar conclusions: SMEFT cross sections are suppressed with respect to non-SMEFT models.
It is important to remark that the precise input a 1 = 2a = 2.1 is not critical for the non-SMEFT models BP1 and BP2, obtaining essentially the same cross sections if one considers the SM limit a 1 = 2a = 2.This is not true for SMEFT, which relies on very fine-tuned relations between the a j effective couplings [21,22].In particular, imposing the SM a 1 value in SMEFT (D=6) leads to â2 = 0 and a total suppression of the cross section.

BP study for ωω → 3h
Figure 4 shows our predictions for the ωω → 3h cross sections for SMEFT, where we are using the D = 6 coupling values given in (4.2).This prediction is then compared in that plot with the pure-HEFT scenarios, the two non-SMEFT BP's discussed in this article.We plot the predictions for BP1 (a 1 ) and BP2 (a 1 ) .
In this case, we are also showing the predictions for the refined BP's, BP1 (a 1 ,a 2 ) and BP2 (a 1 ,a 2 ) , as the ωω → 3h cross section now also depends on a 3 .In the case of the BP2 fractional flare functions we find a very similar result and a cross section 5 orders of magnitude larger than the SMEFT one.
The BP1 exponential flare functions show an interesting peculiar behaviour.The BP1 (a 1 ) model, with F(h) = exp {a 1 h/v} leads to a ωω → 3h amplitude that vanishes for any value of a 1 .Thus, it has not been represented in the logarithmic plot of figure 4. On the other hand, the cross sections for other (n ̸ = 3) multiplicities lead to a large cross section with respect to the SMEFT (see Figs. 3 and 6).Although this may resemble the dilaton-Higgs case (a 2 = b = a = a 1 /2 and a k≥3 = 0) [87][88][89][90][91], such a model leads to a vanishing ωω → n × h cross section for any number n of final Higgs bosons.This pathology is cured in the refined exponential model BP1 (a 1 ,a 2 ) , which now provides again a sizable cross section, orders of magnitude larger than the SMEFT one.
The former discussion shows indeed the full generality of the HEFT approach.This implies that, in most cases, the presence of 'arbitrary' O(1) a j couplings leads to large cross sections in comparison with those obtained in SMEFT calculations.In the latter, the a 1 and a 2 couplings are O(1) but there are very fine-tuned correlations among all the a j 's [21,22], leading to a very suppressed SMEFT cross section for multi-Higgs production.
An interesting fact of the HEFT models is that, although they typically lead to large cross sections, this is not always necessarily true.There are certain cases with cross sections even more suppressed than the SMEFT one, for instance, the dilaton model.
Our conclusions are that, as expected, the SMEFT cross section is suppressed by several orders of magnitude with respect to the non-SMEFT ones: in general, we will see that the fractional flare functions (BP2) tend to yield higher cross sections than the exponential ones (BP1), being both of them well over the SMEFT determinations.
It is important to remark that the precise input a 1 = 2a = 2.1 is not critical for the non-SMEFT models BP1 and BP2, obtaining essentially the same cross sections if one considers the SM limit a 1 = 2a = 2.This is not true for SMEFT, which relies on very precise relations between the a j effective couplings [21,22].In particular, imposing the SM a 1 value in SMEFT (D=6) leads to â2 = 0 and a total suppression of the cross section.
To end this subsection, figure 5 presents a scan of the ωω → 3h cross-section in terms of a 3 , analogous to that performed in ref. [11] for ωω → 2h with a 2 .For the numerical analysis we set a 1 and a 2 to the SMEFT (D=6) values in (4.2).One can observe that the cross section has a zero for a 3 = 2 3 a 1 a 2 − 1 4 a 2 1 ≈ 0.14 (for these numerical inputs).Close by, one finds the SMEFT (D=6) BP.But one can easily observe how fine-tuned this value is: a ±20% variation in a 3 w.r.t. the SMEFT (D=6) value leads to a rising in the cross section of orders of magnitude, even if a 1 and a 2 are left unchanged.Obviously, O(1) values for a 3 (such as, e.g., those in BP1 (a 1 ,a 2 ) in (4.5) and BP2 (a 1 ,a 2 ) in (4.7) produce even larger cross sections, as one can see in figure 4.

BP study to ωω → 4h
Figure 6 shows the comparison between the SMEFT (D=6) BP and the non-SMEFT BP's, BP1 and BP2.One can observe that in this case the results from matching just the a 1 coupling in (4.2) (BP1 (a 1 ) and BP2 (a 1 ) ) or matching a 1 and a 2 (BP1 (a 1 ,a 2 ) and BP2 (a 1 ,a 2 ) ) are close to each other, respectively.It is also important to notice that on the contrary to want happens with the SMEFT BP, which works on very fine-tuned cancellations, the results for the non-SMEFT BP's are very stable: one could set the SM values a 1 = 2 and a 2 = 1 and the predictions for BP1 and BP2 remain essentially unchanged.
In general one can observe that, as expected, the non-SMEFT cross sections are orders of magnitude larger than the SMEFT one.As mentioned above, the a j couplings are highly correlated in SMEFT, so the outcome is extremely sensitive to the precise value of each of them.Small variations on any of the relevant a j from the values prescribed by SMEFT lead to orders of magnitude modifications in the ωω → 4h cross section.We illustrate BP2 (a 1 ) BP2 (a 1 , a 2 ) Figure 4. Comparison of the ωω → 3h cross section predictions of SMEFT (D=6) , BP1 (a1,a2) , BP2 (a1) and BP2 (a1,a2) .We note that the ωω → 3h cross section is identically zero for the first exponential model, BP1 (a1) , regardless of the value considered for the a 1 input and, therefore, it has not been plotted. .We note that, in between, σ ωω→3h vanishes at a 3 = 2 3 a 1 a 2 − 1 4 a 2 1 = 0.1365.[fb] BP2 (a 1 ) BP2 (a 1 , a 2 ) Figure 6.Comparison of the ωω → 4h cross section predictions of SMEFT (D=6) , BP1 (a1) , BP1 (a1,a2) , BP2 (a1) and BP2 (a1,a2) .this through the a 4 scan in figure 7, analogous to the a 2 scan for ωω → hh in [11].If one considers for the SMEFT (D=6) BP in (4.2) (a 1 = 2.1, a 2 = 1.2, a 3 = 0.1 Û 3, a 4 = 0.0 Û 3) one obtains a very suppressed crossed section σ ωω→4h ∼ 10 −2 fb.We will now leave a 1,2,3 fixed to the SMEFT (D=6) values and vary a 4 .A slight increase in a 4 from its SMEFT (D=6) value leads to a decreasing of the cross section by an order of magnitude.However, in this case its minimum is not zero and it is found at a 4 = 3 4 a 1 a 3 − 5 12 a 2 1 â2 + 1 3 â2 2 (1 − χ 1 ) ≈ 0.0344.From that point on, one can see that increasing or lowering a 4 leads to a rising of the cross section by several orders of magnitude.Figure 7 shows that a ±20% variation in a 4 w.r.t. a

SMEFT(D=6) 4
increases the cross section by two orders of magnitude.

Exclusion plots for SMEFT scenarios
In this section we will assume that the underlying UV theory accepts a SMEFT description in the IR.We plan to provide some illustrative estimates of how large the ωω → n × h cross sections can be according to the present experimental observations for a 1 [74,75] and a 2 [77][78][79].More precisely, current LHC hh production analyses have actually shown an important sensitivity to â2 = a 2 − a 2 1 /4, as one can see in figure 2 (see also [67]).This article is focused on the physics of the hard subprocesses ωω → n × h and a full simulation for LHC and future colliders is left for a future work.For this reason, we will just consider some illustrative rough bounds from h → W W [74,75] and double-Higgs production at the LHC (figure 2): .The cross section's minimum is not zero this time and it is found at In the study in this section we will consider these inequalities under the SMEFT perspective and use them to provide predictions for multi-Higgs production cross sections.For this, we will consider the contribution to the amplitudes from SMEFT operators at the corresponding lowest non-trivial order (1/Λ 2 for 2h, 1/Λ 4 for 3h, etc.).Higher orders introduce corrections which might be relevant for high precision physics but are not going to teach us much here.
In first place, under SMEFT, an experimental restriction on the h → W W coupling implies also a limitation on the allowed range for a 2 and, more specifically, on â2 .At lowest order in the SMEFT expansion, this implies: One can readily see that this implies a bound on the SMEFT D = 6 Wilson coefficient, parametrized by H□ /Λ 2 in (3.2).Nothing is known about the coupling a 3 (or a 4 , a 5 , etc.), although triple-Higgs production analyses such as [69] can be used to assess the sensitivity of current and future colliders to that parameter, a 3 .Indeed, the relevant parameters for 3h and 4h production are not actually a 3 and a 4 but rather the combinations â3 = H□ ) 2 -or the latter is absent [72]-leads to a very different phenomenology driven by the 1/Λ 4 corrections in the a j 's, where they now follow different correlations.This scenario is discussed in appendix F but will not be further analyzed here.For this reason, for the illustration in this section, we will consider for our numerical study the simple bounds: We now consider the allowed values of the parameters d and ρ and maximize the SMEFT 2h, 3h and 4h cross sections at their lowest non-trivial order in the EFT expansion: (5.4) ) For each energy the cross section maxima are reached for d = d max and ρ = ρ max .
Here we have considered the parameters to be energy independent.In this region, the maximum cross section are where the maximum is reached at any ρ ∈ [−ρ max , ρ max ] for 2h and at ρ = ρ max for 3h and 4h (for 1 + χ 1 ≥ 0). Figure 8 shows these maximized cross sections as a function of the center-of-mass energy √ s.At low energies there is a hierarchy in the multi-Higgs production cross sections, where 2h production is larger than 3h, and this larger than 4h.However, as the energy increases there is a moment in which all three cross sections become roughly of the same order and then the hierarchy gets inverted being, according to the theoretical expression, more likely to generate 4h.This feature only highlights that from that point on, the EFT should not be trusted anymore, as we have reached the cut-off of the theory.At that point corrections from higher orders are as big as the contributions already included in (3.7), (3.9) and (3.11).Beyond that point, one should not further rely on the information provided by this plot.Nevertheless, the √ s < 10 1 TeV range in figure 8 shows the exclusion plot for the multi-Higgs production cross section one can expect for the allowed SMEFT range in (5.3).
Figure 8 may be a little bit unsatisfactory, as the high energy part of the plot lies off the EFT applicability region.Thus, one might want to study the EFT prediction for all  possible energies whenever the EFT is still applicable.For this, at a given energy, one needs to establish a practical criterion on what "EFT-valid" means.For theories where one has D = 6 contributions it seems fair to assume that one has an EFT expansion in powers of Hence, for that given √ s, we will allow all values for d that obey, c for some EFT-expansion tolerance parameter ϵ that controls the rate of convergence.This implies that for each value of the energy we are going to consider a range for d that varies with s in the form, .11)This choice of the allowed range for d ensures that the EFT analysis is valid at all considered energies.If now substitute this maximum value d max = d max (s) in the previous cross section bounds (5.7)-(5.9),we obtain the modified expressions, ) which ensure the EFT validity at all explored energies.Figure 9 shows the exclusion regions for 2, 3h and 4h cross sections for the bound ρ max = 1 and the tolerance ϵ = 0.1.

A collider estimate illustration: e + e − CLIC at 3 TeV
In this last subsection we want to provide an illustration of the expected cross sections in an actual collider, not just for the hard subprocess ωω → n×h.We do not intend, however, to perform a full detailed simulation of the expected result in a fully realistic collider.Its complexity relegates that analysis for a future work.In this article we are going to focus on the rough estimate provided by the effective W approximation (EWA) [96][97][98][99][100], where the longitudinal gauge bosons are collinearly radiated from the incoming particles.This kinematical region is expected to dominate the full cross section and provide a reasonable estimate.
Regarding the type of collider, for the sake of simplicity, we will consider a lepton collider such as CLIC (e + e − ) with √ s tot = 3 TeV [101].In Stage 3 scenario the integrated luminosity for that center-of-mass total energy may reach up to 5000 fb −1 [102].This roughly means that cross sections with σ < 0.2 • 10 −3 fb are not expected to be observed.In practice, the lowest observable bound is actually even higher, as one needs to take into account the efficiency of the detectors in the reconstruction of the Higgs decays.
the total cross section, σ tot , is provided by the hard subprocess cross section times an appropriate W L -luminosity function of the form [100], This expression provides the dependence of the differential total cross section on the W W invariant mass, √ s, which has been plotted in Fig. 10 for the production of 2, 3 and 4 final Higgs bosons.We are providing the predictions for the SMEFT (D=6) benchmark (dashed lines) discussed in this and previous sections, considering d = 0.1 and ρ = 1.It is possible to observe the cross-section suppression for SMEFT scenarios as we increase the number of final Higgs bosons.For illustration we also provide the differential cross section for a pure-HEFT benchmark, BP2 (a 1 ) .The comparison of the two scenarios shows again that, in general, pure-HEFT models lead to cross sections which are orders of magnitude larger than those in theories that admit a low-energy SMEFT description.We note that all EWA differential cross-sections vanish at the end-point s = s tot .
To present these results in a more transparent form, in Fig. 11, we provide the integrated cross section over bins of 0.5 TeV.For the pure-HEFT scenario BP2 (a 1 ) one finds cross sections ranging in the order from 10 −2 to 10 2 fb, which could be potentially observed with a 5000 fb −1 integrated luminosity.In the SMEFT (D=6) benchmark this only happens for the two-Higgs production: the production of 3 or more Higgs bosons shows cross sections well below the attobarn by several orders of magnitude.In this sense, the observation of multi-Higgs states in that referred collider would indicate that the underlying UV theory could not support a SMEFT-like low-energy description.It should then be described through the more general HEFT formalism.
For illustration, we have selected a specific lepton collider, providing a clean and clear result.An analogous analysis can be performed for other lepton colliders, as well as for current or future hadron accelerators.However, extending this analysis to hadron colliders requires further discussion on quark parton distribution functions and, to make these predictions useful, a good estimation of detector, reconstruction, and analysis effects.Although very interesting for the proposals of future hadron colliders, such a study is beyond the scope of this work.

Conclusions
In this article we have performed an EFT analysis for longitudinal VBS into several Higgs bosons.In order to clarify the relevant interactions at energies well over the production threshold, we have employed the EqTh and considered a series of well-motivated approximations, such as assuming the mass corrections to be subdominant with respect to momentum dependent interactions.These approximations have allowed us to relate the VBS with the EW Goldstone scattering process ωω → n × h.
In section 2, we have presented the calculations for 2h, 3h and 4h production from ωω in the HEFT approach at its lowest order.This means we have considered only tree-level processes and contributions from HEFT operators with up to two derivatives, O(p 2 ).Loops or operators with a higher number of derivatives lead to further subleading corrections, neglected in this article.
We provide analytical expressions for the scattering amplitudes and cross sections, the only exception being the ωω → 4h cross section which is provided through an analytical expression in terms of two coefficients χ 1,2 , which must be computed numerically.A code (MaMuPaXS [54]) has been specifically developed for the phase-space integration with the massless particle approximation assumed here.This code greatly reduces the computational resources required to sample the cross section at different energies.
We want to note that we have identified the relevant combinations â2,3,4 that rule the ωω → 2h, 3h, 4h scattering amplitudes.These are the relevant parameters one should aim for in future experimental analyses, and not the W W h, W W h 2 and W W h 3 couplings a 2,3,4 , respectively.Note that, trying to measure such couplings experimentally leads to unwanted strong correlations in the data analyses.We have shown that through scalar field redefinitions in the effective action (changes of coordinates in the scalar manifold) one can derive these precise relevant combinations âj (see appendix C).Indeed, it is easy to extract, e.g., the relevant combinations â5 and â6 for 5h and 6h production even before starting the computation of this amplitude.
In section 3 we restricted those results to the case when the EFT accepts a SMEFT representation, this is, when one can write down the effective Lagrangian in terms of a Higgs doublet and sort out the operators according to their canonical dimension.It is obviously implicit that, in addition to the latter consideration, this SMEFT power counting in canonical dimension gives a convergent EFT expansion for the observables, which does not always need to be true for an arbitrary EFT.As before, we focus on a regime well above threshold where momentum dependent interactions dominate over non-derivative ones.We compute the contributions to the tree level ωω → 2h, 3h, 4h amplitudes from operators of canonical dimensions D = 6 and D = 8.We provide the results for the amplitudes and total cross sections in compact analytical expressions.Remarkably, we find that under our approximations the multi-Higgs productions is highly suppressed in the SMEFT, since producing a large number of Higgs bosons requires the insertion of an increasing number of D = 6 operators of terms of higher dimension (D = 8, 10, etc.).
In section 4 we have considered a series of benchmark points to compare scenarios where the effective theory accepts a SMEFT description and pure-HEFT scenarios, where the theory cannot be written in the SMEFT form.The study of VBS into 2h, 3h and 4h shows that in general the non-SMEFT cross sections are orders of magnitude larger than the SMEFT predictions.This happens even if we set the W W h and W W h 2 couplings identical to the SMEFT ones.Moreover, in the non-SMEFT scenarios the cross-section values do not dramatically depend on the precise values chosen for the a j and they are stable under small modifications of these parameters.By contrast, SMEFT highly rely on very fine-tuned correlations between the W W h n couplings, which lead to very precise cancellations in the cross sections that drop their value by several orders of magnitude.
Although in some cases, the amplitude for a specific number of Higgs bosons can be abnormally small in a particular model (e.g., ωω → 2h for BP1 (a 1 ) ), in general for any other number of final Higgs bosons non-SMEFT cross sections are orders of magnitude larger than the SMEFT predictions.Thus, a comparative analysis of VBS into channels with different numbers of Higgs bosons can tell us whether the (non-resonant) experimental data can be described through SMEFT.Indeed, section 5 shows the exclusion cross section plots for ωω → 2h, 3h, 4h if one assumes SMEFT: as it does not allow large cross sections in multi-Higgs production one finds that the maximum allowed SMEFT cross section is much more suppressed than the pure-HEFT one.For our BP's we found SMEFT yields cross sections below the attobarn for √ s > 2 TeV for VBS into more than two Higgs bosons (see figure 9).For illustration, we have studied the multi-Higgs production cross section in a 3 TeV e + e − collider (CLIC) within the EWA, where the longitudinal W + L W − L are radiated collinearly from the incoming e + e − .The cross section estimates in Fig. 11 lead to the same conclusions found in previous sections: multi-Higgs production is strongly suppressed in SMEFT scenarios in comparison with pure-HEFT cases.
In summary, it is complicated to produce large deviations from the SM in multi-Higgs VBS within the SMEFT framework.An eventual experimental observation of cross sections that exceed the bounds prescribed by SMEFT (Figs. 8 and 9) would be an indication that a broader approach is needed.Thus, the joint study of the 2h, 3h and 4h production channels can serve as a smoking gun indicating that (non-resonant) new physics, if present, requires an interpretation beyond SMEFT, requesting the use of the full HEFT machinery.We are confident that this work can serve as a motivation for studies in future colliders, illustrating what possible (non-resonant) new physics could we find.
There are several directions that can be followed as next steps.First, one should go beyond the EqTh approximation and consider mass corrections.This would allow us to also study the region close to threshold and improve the accuracy of our high energy EFT predictions.In addition, the cross sections studied here referred to the hard subprocesses W W → n × h while in colliders, either leptonic or hadronic, the W 's need to be radiated from the colliding particles, diminishing the total n × h production cross section for the VBS category.Likewise, simplifications as the effective W approximation may not be fully satisfactory and a full MonteCarlo simulation may be required for a precise LHC estimate.One may find then that the signal studied in this article may appear surrounded by additional reducible and irreducible background that needs to be disentangled.But nevertheless, the signal is there and it should give some definite prediction for pp or e + e − colliders if the experimental data and MonteCarlo simulations are improved enough.A final obvious improvement is the consideration of next-to-leading order and loop corrections in the EFT [103][104][105], which would further refine the predictions.
Galileo Galilei Institute for hospitality and support during the scientific program on "Theory Challenges in the Precision Era of the Large Hadron Collider," where part of this work was carried out.

A Phase-space relations
The Lorentz invariant phase-space for n final particles is given by with E j = ⃗ p 2 j + m 2 j .In the massless case (m j = 0) one has It is then possible to provide an analytical expression for the full phase-space hypervolume [92]: 2) It will be convenient in this massless case, for the formal understanding of the energy dependence and the numerical simulations, to factorize out the global s m dependence of the phase-space integral through the change of variables p i = p i / √ s and q = q/ √ s, defining: with the delta imposing the energy conservation relation i ‹ E i = 1 and, for the center-ofmass frame, the three-momentum conservation i ⃗ p i = 0.For the simplest case, n = 2, this expression recovers the standard result V 2 = V 2 = dΠ 2 = V 2 = (8π) −1 .The advantage of this s-independent phase-space is that the integration of an arbitrary phase-space of n particles will always handle O(1) numbers and trigonometric functions.Thus, if one can factorize out the s dependence of the scattering amplitude, one does not need to repeat the dΠ n (⋄) phase-space integration for every energy.

B 4-momentum conservation: kinematical relations
Given a process ω(k 1 ) ω(k 2 ) → h(p 1 )...h(p n ).In the rest frame, the center of mass squared energy, We define the rest-frame three-momentum fractions as for each outgoing Higgs boson; the functions where θ i is the angle between the i-th Higgs boson and the first ω Goldstone boson momentum, ⃗ k 1 (that is, z 1 = 1 − cos θ, z 2 = 1 + cos θ as usual in a two-body problem with t and u channels); and the functions where θ ij is the angle between the i-th and j-th Higgs bosons.
Taking the masses of both the Golstones and Higgses as zero, we have the 4-momentum products: where no implicit summation is assumed in repeated indices.These expressions allow us to write down these functions in terms of Lorentz invariant products: with q ≡ k 1 +k 2 .If one wants the amplitudes in the main text in terms of Lorentz invariant scalar products one just needs to perform the replacements above.Using 4-momentum conservation, we can obtain some relations between the newly defined functions.Firstly, we have Also, from conservation of the total energy in the rest frame, we get for all j ∈ {1, ..., n} and (B.9), we obtain As an example, for three Higgs bosons in the final state, the relevant relations are Furthermore, we can use this relations to write B of eq.(2.8)only in terms of f 1,2,3 , z 1,2,3 , z 12 and z 13 as This expression for B has been checked analytically with Maple and numerically with the integration code in ref. [54].
C HEFT simplifications through field redefinitions: further checks and understanding of ωω → n × h: As stated in the main text, we have used the standard form of the HEFT Lagrangian (2.1) to compute every possible tree-level Feynman diagram contributing to the processes discussed in this article.We have performed the calculation within the equivalence theorem and neglected non-derivative interactions, keeping only the O(E 2 ) contributions to the amplitudes.
After some algebra, we have been able to simplify the amplitudes, finding that both for 2h and 3h production the amplitudes are pure s-waves.The crossed-channel Goldstone exchanges cancel out and we are effectively left with a contact interaction (in spite of the large variety of crossed channel diagrams).For the case with 4h in the final state, we also find strong cancellations and the amplitude is effectively reduced to a contact ωω → 4h interaction plus a contribution with one Goldstone exchange in the crossed channel.
In this appendix we proceed to clarify the origin of such cancellations, which lead to the extremely simple structure of our results.The part of the HEFT Lagrangian (2.1) relevant for the ωω → n × h scattering has the form, with the flare function As an alternative to this canonical HEFT Lagrangian (2.1), we realized that these processes can be simplified by performing a non-linear field redefinition [106][107][108][109][110] in (2.1) before starting the amplitude computation.The aim of these redefinitions is to simplify and eliminate the derivative three-particle vertex hωω [111,112].In this way, the number of possible diagrams for a given process is greatly reduced 3 .
For this simplification we have considered transformations of the form, with the free dimensionless real constant N and the O(h) function g(h) determined by the flare function F(h) through the relation ã .
(C.3) Interestingly, the application of this transformation to the Lagrangian (C.1) leads to a new Lagrangian with exactly the same structure in (C.1), but with a new function F(h) determined by: F(h) = F(h) 1 + g(h) 2 . (C.4) In particular, we will consider here a very specific N normalization: This choice of the scalar manifold coordinates transforms the original Lagrangian (C.1) into but now in terms of the new function, with coefficients provided by the ones of F(h) through the combinations Note that F(h) has now its first non-trivial contribution at O(h 2 ), in comparison with F(h), where its first non-trivial term starts at O(h).We remind that both the O(h 5 ) and O(ω 4 ) operators are irrelevant for the tree-level ωω scattering producing two, three and four Higgs bosons in the final state.Nonetheless, it is very simple to extract the relevant parameters for any ωω → n × h scattering: one just needs to extract g(h) from (C.3), taking N = a/2 and including the power series terms up to order h n , the latter included.One then substitutes this g(h) in (C.4) to obtain F(h) up to that order.Although this work does not consider production of a higher number of Higgs bosons, the procedure above allows us to extract the relevant combination of the flare function coefficients, a j , for a generic amplitude ωω → n×h: 1st) compute g(h) up to O(h n ) by plugging F(h) in (C.3) up to that order with N = a 1 /4; and 2nd) expand F(h) = F(h) (1 + g(h)) 2 up to h n .The corresponding coefficients âj will be the relevant combinations for that process.Data studies that do not take this into account and use directly the a j instead of the âj essentially introduce unnecessary correlations that make the experimental analysis much more involved.As an example, though irrelevant for this work, we show that one can easily extract the next âj coefficients: â5 = a 5 − 1 120 a 1 15a 1 â3 − 16â 2 2 + 96â 4 , â6 = a 6 + 1 180 a 1 7a 1 â2 2 − 27a 1 â4 + 45â 2 â3 − 150â 5 , etc.It is also interesting to note that in the dilatonic model [87][88][89][90][91], where F(h) = (1 + ah/v) 2 , one finds for the transformation above that 1 + g(h) = (1 + ah/v) −2 , leading to F(h) = 1.So that, remarkably, all the âj 's couplings are zero and all ωω → n × h amplitudes vanish at tree level (within the EqTh limit considered in this article).This exact same conclusions apply to the SM, where a = 1.
Given that the generating functional of the quantum field theory is invariant under field redefinitions, both (2.1) and (C.6) Lagrangians produce the same on-shell scattering amplitudes [106][107][108][109][110]. Nonetheless, we have explicitly checked that the Lagrangian in (C.6) exactly reproduces the results in the text for 2, 3 and 4 Higgs production.
In this work we are neglecting the potential V (h) (since E ≫ m h ) and normalizing any pure-Higgs derivative operator of the form f (h) (∂h) 2 into the canonical kinetic form 1 2 (∂h) 2 .Thus, there is no Higgs self-interaction and all the vertices relevant for the processes discussed in this article contain two derivatives and are of the type ωωh n .The critical gain from this field redefinition approach -using the simplified Lagrangian (C.6)-is that the number of diagram topologies is greatly reduced: there is only 1 diagram for ωω → 2h and ωω → 3h, and it is reduced to 2 diagram topologies for ωω → 4h process (up to n! permutations in the labeling of the outgoing Higgs particles in the diagrams).Figure 12 shows all the relevant diagrams generated by (C.6) for ωω → 2h, 3h, 4h, where in every ωωh n vertex one now has the ân effective coupling.c-d) Only two diagrams contributing to the process ωω → 4h.We have used the simplified Lagrangian (C.6) to generate these amplitudes, so every ωωh n vertex carries an ân effective coupling.Note that, in addition, one needs to consider all possible permutations for the assignment of the external particles.
A closer observation of the transformation above shows that, in some scenarios, one can also conveniently choose the normalization N to remove a higher order term, a n h n , from F(h) instead of the first.For instance, provided a 2 < a 2 1 /4, the choice N = a 1 ± » a 2 1 − 4a 2 /4 removes the a 2 h 2 term in F(h), passing this information to the terms of order h 1 and h 3 , h 4 , etc.
Another example is provided by the normalization N = 3 8 a 3 a 2 −a 2 1 /4 , which removes the h 3 term in F(h) and encodes its information in the factors âj now multiplying h 1 , h 2 and h 4 , h 5 , etc.This detail can be important for a proper interpretation of W W → 3h computations: at high energies, in the EqTh, it is possible to describe the ωω → 3h scattering without an ωωh 3 vertex (â 3 = 0) [69], understanding that the ωωh 2 and ωωh couplings are not the original ones (a 2 and a 1 ) but the effective ones in (C.6) (â 2 and â1 ).
The main inconvenience of this approach is that the SU (2) L ×SU (2) R chiral invariance of the action is no longer explicit as in (2.1).The symmetry transformations are more involved in this form.For this reason, the correlations between the a 1 and a 2 couplings in F(h) one finds for SMEFT-type theories with D = 6 contributions (this is, a 2 = 2a 1 − 3 [21,22]) are no longer applicable for F(h) (SMEFT with D = 6 contributions does not fulfill â2 = 2â 1 − 3, as one can see using, e.g., (3.3) and â1 = 0); the chiral operator structures considered in [21,22] are deformed here for the terms of order ω 4 and higher, so the conclusions therein are no longer applicable for the present simplified Lagrangian with F(h).
Aside from these considerations, this type of simplifications may prove useful for the computation of EW processes with Higgs and Goldstone bosons in the equivalence theorem, either at tree or loop levels.

D Multi-Higgs production suppression in SMEFT: further clarifications
In order to obtain the valid tree-level configurations for a fixed number n of final Higgses in ωω → n × h, one has to vary N V 3 from 0 to N max It is illustrative to observe the valid diagrams for the cases with n = 2, 3, 4, 5 Higgses in the final state: • n = 2 allows only diagrams with: N ).

E Massless Multiparticle Cross Section (MaMuPaXS)
E.1 General case with 4 or more particles in the final state In order to numerically integrate the dimensionless function B appearing in ωω → 4h in eq.(2.8), we are using the publicly available code MaMuPaXS [54].The code is based on Fortran, Python and the VEGAS library [113,114] and can be accessed through a GitHub repository [54].
For the integration of B n from (2.11) our goal was keeping the computation as simple as possible.In order to achieve this, the change of variables described below (the MaMuPaXS package) is a straightforward procedure for integrating the phase space if all the particles can be considered massless.Indeed, the only external library we use is VEGAS [113,114].However, there is the option of the RAMBO algorithm [115], a classic Monte-Carlo algorithm for generating phase space points and scattering events for multi-particle processes in the non-massless case.The question arising here is whether, in the massless case, MaMuPaXS (based on solving a 4 × 4 system of linear equations and computing a 4 × 4 determinant -the Jacobian) outperforms RAMBO.But, since we were not bounded by CPU time, this was a collateral issue.Hence, we have currently implemented up to 2 → 5 processes on MaMuPaXS (though this article only studies up to 2 → 4 processes).However, an extension to arbitrary 2 → n processes is feasible and, indeed, should be made in order to compare with RAMBO's performance.This is a desirable goal for a continuation of the MaMuPaXS development, but at this stage we consider that this is beyond the scope of this work.
We will consider a process of the form φ The code implements a change of variables for the 4-dimensional delta function in the phasespace integration (A.1), which can be analytically integrated.The initial state will have total four-momentum q = k 1 + k 2 = ( √ s, 0, 0, 0) in the center-of-mass rest frame.In spherical coordinates, each massless final state particle 4-momenta have the form: The argument of the 4-dimensional delta function is equal to the null four-vector if where n is the number of particles in the final state and the new integration variables are the spherical angles {θ i , ϕ i } and the energies E i , equal to the three-momentum modulus |⃗ p i | in the massless case studied here.
In this work, we consider two different cases.The first one, n ≥ 4, where the previous system of equations is actually a linear system of 4 equations that can be solved analytically, This system has a unique solution for E 1 , E 2 , E 3 , E 4 .After using the delta-function in the phase-space integral to fix E 1 , E 2 , E 3 , E 4 , we are just left with an integration over the spherical angles {θ i , ϕ i } and the remaining energies {E 5 , E 6 , . . .} (if any).Furthermore, because of the symmetry of the problem (two colliding massless scalars along the z-axis, producing scalar particles in the final state), the scattering is independent of the absolute azimuthal orientation: We can fix one of the ϕ i angles in all the cases (n = 3, 4, 5, . . . ) without any loss of generality.The integration of this azimuth yields a 2π factor in the phase-space integration.
Hence, we will change the integration variables of eq.(A.1) to spherical coordinates, where the total four-momentum q = ( √ s, 0, 0, 0) T and the p j = (E j , ⃗ p j ) T must be understood as column vectors.F : R4 → R 4 is a linear function on E 1 , E Hence, the 4-dimensional delta can be written as The 4-dimensional delta can be computed via solving E 1 , E 2 , E 3 and E 4 from equation (E.4), and multiplying the integrand by the inverse of the Jacobian, (det A) −1 .

E.2 Case with n = 3
If n = 3, the system of 4 linear equations may be, in principle, overdetermined.Hence, one of the angles will be a function of the other angles.Actually, this additional restriction is the well-known fact that the 3-momenta of the 3 particles in the final state must be coplanar in the center-of-mass rest frame.This coplanarity compatibility requirement can be written as: cos θ 1 cos θ 2 cos θ 3    = 0 .(E.9) This equation implies an additional restriction that fixes one of the angular variables as a function of the others.The energies E 1 , E 2 , E 3 are fixed in the usual way by means of (E.4).
In the particular case n = 3, the F : R 4 → R 4 function that enters on eq.(E.6) is defined on E 1 , E 2 , E 3 and on the angle θ 3 .The roots of F (E 1 , E 2 , E 3 , θ 3 ) − ( √ s, 0, 0, 0) T can still be analytically found.The Jacobian of F can also be computed, so that we still integrate on spherical coordinates.The only difference is that one of the spherical angles will be fixed by the Dirac delta.

E.3 Code repository
All the Jacobians for n = 3, 4 and 5, as well as the required change of variables and the computations of E 1 , E 2 , E 3 and E 4 (for n ≥ 4) or θ 3 (for n = 3) is coded in Fortran on [54].This Fortran code is called from Python via F2PY 4 , and then integrated via VEGAS [113,114].
For the current version, the amplitudes must be hard-coded inside certain FORTRAN files specified in the repository documentation.Then, the FORTRAN code has to be compiled to a Python module.Afterwards, the VEGAS library for Python can be used for the numerical integration.However, since we think our code may be useful for the high energy physics community, we are planning in a future version to build a standalone FORTRAN library so that it can be used without recompiling, either within Python or linking to a Fortran program.The details about the usage of the code are explained on the GitHub repository.
F SMEFT starting at D = 8 In some particular scenarios we may find that the leading corrections (D = 6) to the flarefunction are zero and the first non-vanishing corrections related to the flare function stem from D = 8 SMEFT operators [72].For theories that have such a SMEFT description starting at O(1/Λ 4 ), one finds, up to higher order corrections: H□ in the result from ref. [21,22].This is, in these scenarios one sets d = d 2 = 0 while keeping d 2 ρ ̸ = 0.This type of additional SMEFT suppression is found for these couplings, for instance, in theoretical frameworks such as the 2-Higgs Doublet Model [72].
) up to corrections with higher powers in s.The C n (z; c j ) are O(1) dimensionless combination of angular functions the z = {z i , z jk } trigonometric functions previously defined (see appendix B) and the O(1) SMEFT Wilson coefficients c j (where Λ 4−D has been factored out from the SMEFT Lagrangian term as it is customary).The exponent γ n indicates how many insertions of the dimension-6 SMEFT operators one needs to introduce in the diagram to obtain the first non-vanishing contribution to ωω → n × h.It follows the pattern γ 2 = 1, γ 3 = γ 4 = 2, γ 5 = γ 6 = 3...This can be summarized in the general expression γ n = ⌈n/2⌉ (see below and appendix D).

Figure 2
Figure 2. a) CMS experimental confidence regions for the hW W coupling κ V = a = a 1 /2 and for the hhW W coupling κ 2V = b = a 2 , from a non-resonant hh production search with each Higgs boson decaying into a highly boosted b b pair [77] (white lines and colour map in figure 11 from the additional material in CMS-B2G-22-003).b) CMS confidence regions from processes with one Higgs boson decaying to b b and the other one to W + W − [78] (black lines and colour map).c) ATLAS experimental confidence regions, from the combination of hh → b bb b, b bτ + τ − , b bγγ decay channels [79] (red lines).d) CMS searches for non-resonant production of Higgs boson pairs in b bτ + τ − final states [80] (black lines and colour map confidence regions).We have superimposed the correlation a 2 = 2a 1 − 3 in SMEFT at O(1/Λ 2 ) (red line in a), b) and d); blue line in c)).In addition, we have also plotted the parabolas â2 = a 2 − a 2 1 /4 = 0 (solid black in a) and c); solid green in b) and d)) and â2 = ±0.2(dashed black in a) and c); dashed green in b) and d)), which determine the longitudinal W W → hh scattering in the naive equivalence theorem.Reproduced under the Creative Commons BY 4.0 license.

4 3 d 2 ( 1 +) 2 ó 2 (d 2 8 ) 6 )
ρ) + O(d 3 ) and â4 = 1 3 d 2 (1 + ρ) + O(d 3 ) provided in eq.(3.3).One can identify the two types of SMEFT contributions to ωω → 3h and ωω → 4h: a single insertion of one D = 8 operator proportional to c (8) H□ parametrized by ρ ≡ c (d 2 ρ terms); and a double insertion of D = 6 operators proportional to terms without ρ).In the numerical analysis that follows we will be taking the assumption that the c (H□ terms are similar in size to the (c (H□ ) 2 contributions, or smaller.The case where c

Figure 8 .
Figure 8. SMEFT exclusion plot for the cross sections for 2, 3 and 4 Higgs bosons with |d| ≤ d max = 0.1 and |ρ| ≤ ρ max = 1.The regions above the solid, dashed and dotted lines can be safely excluded if the Wilson coefficients are within the considered range.Notice that the EFT perturbativity condition is not considered in this figure, as the EFT expansion breaks down on the region past the crossing point.

Figure 9 .
Figure 9. Exclusion plot for the maximum value of the cross sections for 2, 3 and 4 Higgs bosons with the constraint |ρ| ≤ ρ max = 1 and EFT-expansion tolerance ϵ = 0.1.

Figure 12
Figure 12. a) Only diagram contributing to the process ωω → 2h.b) Only diagram contributing to the process ωω → 3h.c-d)Only two diagrams contributing to the process ωω → 4h.We have used the simplified Lagrangian (C.6) to generate these amplitudes, so every ωωh n vertex carries an ân effective coupling.Note that, in addition, one needs to consider all possible permutations for the assignment of the external particles.