Exploring extended scalar sectors with di-Higgs signals: a Higgs EFT perspective

We consider extended scalar sectors of the Standard Model as ultraviolet complete motivations for studying the effective Higgs self-interaction operators of the Standard Model effective field theory. We investigate all motivated heavy scalar models which generate the dimension-six effective operator, |H|6, at tree level and proceed to identify the full set of tree-level dimension-six operators by integrating out the heavy scalars. Of seven models which generate |H|6 at tree level only two, quadruplets of hypercharge Y = 3YH and Y = YH, generate only this operator. Next we perform global fits to constrain relevant Wilson coefficients from the LHC single Higgs measurements as well as the electroweak oblique parameters S and T. We find that the T parameter puts very strong constraints on the Wilson coefficient of the |H|6 operator in the triplet and quadruplet models, while the singlet and doublet models could still have Higgs self-couplings which deviate significantly from the standard model prediction. To determine the extent to which the |H|6 operator could be constrained, we study the di-Higgs signatures at the future 100 TeV collider and explore future sensitivity of this operator. Projected onto the Higgs potential parameters of the extended scalar sectors, with 30 ab−1 luminosity data we will be able to explore the Higgs potential parameters in all seven models.


Introduction
The discovery of the Higgs boson at the Large Hadron Collider (LHC) marked the discovery of the last missing piece of the Standard Model (SM). Precision measurements of the Higgs couplings are a major goal for current and future high energy experiments. Current experimental results provide strong evidence that the nature of the Higgs boson is consistent with the predictions of the SM. The measurement of this behavior is entirely dependent on single Higgs phenomena through precision measurements of the Higgs couplings to the vector bosons and the SM fermions. On the other hand, the Higgs self-interactions, responsible for Electroweak Symmetry Breaking (EWSB), still remain undetermined experimentally.

JHEP05(2018)061
The Higgs self-coupling directly determines the shape of the Higgs potential and therefore measuring possible deviations of the Higgs self-coupling from its SM value is a crucial step in understanding the nature of EWSB, electroweak vacuum stability, and the nature of the electroweak phase transition (EWPT).
In order to investigate the generic features of the trihiggs coupling at the LHC and a future collider we adopt an Effective Field Theory (EFT) approach [1][2][3][4]. In doing so we assume some possible new physics beyond the SM which modifies the Higgs couplings and is heavy with, for example, new physics scales such as Λ NP ∼ TeV. The effects of the new physics are parametrized by higher dimension effective operators, and the dimension-six Q H = (H † H) 3 operator is the leading operator which modifies the momentum independent Higgs self-couplings at low energy. The Q H operator remains the only operator related to the Higgs sector unconstrained by current experiment. In order to motivate this study of the effective operator Q H we consider ultraviolet (UV) complete models which may generate this operator at tree level and therefore with a larger Wilson coefficient. This requirement combined with Lorentz invariance then limits our consideration to extended scalar sectors. 1 Additionally the new scalar must not be charged under SU(3) c as closure of color indices requires Q H be generated at one loop. Such scalar extensions of the SM constitute relatively simple scenarios beyond the SM which are also well-motivated by studies of the electroweak phase transition and baryogenesis [5,6], having dark matter candidates [7][8][9], or mechanisms for neutrino mass generation [10][11][12][13]. The complete list of the scalar extensions which generate a tree-level Q H are real [5,7,8,14] and complex singlets [6], the two Higgs doublet model (2HDM) [15][16][17], real [9,18] and complex [10][11][12][13] triplets, and complex quadruplets. Assuming the new scalars in these models are heavy, we utilize an EFT approach to study their effects on electroweak precision tests, modifications of the single Higgs couplings, and the di-Higgs production process in a model-independent and predictive way.
Many new physics models with SM-compatible single Higgs phenomena could exhibit di-Higgs phenomenology distinct from that of the SM [19,20]. The modifications of the Higgs trilinear couplings can only be directly observed in Higgs boson pair production, therefore the di-Higgs process at the LHC and future colliders is the only direct way to measure the Wilson coefficient of the effective Q H operator. Alternatively the trilinear Higgs coupling can be studied indirectly [21,21,22,[22][23][24][25]. However our paper will focus on the direct constraints, we discuss the indirect constraints briefly at the end of section 2. The di-Higgs production mechanism at hadronic colliders is dominated by the gluon fusion process which includes the triangle and the box contributions from the top quark. Due to destructive interference between these two contributions, the di-Higgs production cross section in the SM is typically small and thus challenging to observe in the near future. However, in the scalar extended models, the di-Higgs cross section may be increased considerably making measurement a possibility at the proposed 100 TeV collider [26,27]. In this paper we investigate the di-Higgs production cross sections in the EFT framework, and study the discovery potential of the Wilson coefficients in the EFT at the proposed 100 TeV collider. 1 Requiring closure of spinor and Lorentz indices implies fermions may only generate the QH operator at one-loop and vectors may only generate dimension-six-operators with two derivatives at tree-level.

JHEP05(2018)061
We proceed by, in section 2, listing the simplest scalar extensions of the SM which generate Q H at tree-level along with their Lagrangians and the corresponding effective Lagrangians after integrating out the new heavy degrees of freedom. Then in section 3 we study the implications of single-Higgs measurements on the corresponding EFTs as well as the implications of these constraints on the UV complete models. In section 4 we study these EFTs' impacts on di-Higgs production at the proposed 100 TeV collider. Finally our conclusions are found in section 5.

The effective Lagrangian
We consider all ultraviolet (UV) complete models which include one additional heavy scalar which generate, after integrating out the new scalar, dimension-six operators affecting the trihiggs vertex at tree level. In order to generate tree level dimension-six operators one needs a term H 2 S or H 3 S, where S is the new heavy scalar, this is a result of all other models having an additional Z 2 symmetry due to the requirements of the gauge symmetry and renormalizability. 2 The relevant theories are then, real and complex scalar singlets, the two-Higgs doublet model, real and complex scalar triplets of SU(2) L with hypercharge Y = 0 and Y = −1 respectively, and finally complex scalar quadruplets of SU(2) L with either Y = 3/2 or Y = 1/2. For each model we write down the Lagrangians for each UV-model along with the corresponding effective field theory (EFT) to dimension-six at tree level, we will only write the new terms in addition to the standard model terms for convenience. In writing the EFTs we will follow the procedure of Henning et al. [28,29]. To clarify our notation and conventions, we write here the general Lagrangian for all UV complete models, neglecting SM fermionic and gauge boson terms, considered: (2.1) Where ∆L contains all terms containing new fields (in the case of the models we consider this is one new scalar multiplet of SU(2) which may or may not have hypercharge). µ 2 becoming negative signals spontaneous symmetry breaking leading to the massive gauge bosons of the SM. After deriving the EFTs we employ the Warsaw basis [2] for the dimension-six operators, translations between the various bases are included throughout much of the recent literature including a package for relating the bases [30]. The operators which are relevant to our analyses are: (2. 2) The fermionic operators should be summed over each generation with an appropriate Wilson coefficient. In general the fermionic operators can have off diagonal components, however for the models considered this is only possible for the two-Higgs doublet model and JHEP05(2018)061 we will employ particular choices of the fermionic matrices in the model to suppress off diagonal components, as is motivated by studies of flavor changing neutral currents, and therefore assume these operators to be diagonal.

Real scalar singlet
The real scalar singlet has Y = 0, it has been studied extensively in the literature both from the UV complete [5,7,8,14] and EFT perspectives [28,31,32]. The Lagrangian, neglecting SM terms, is given by: After integrating out the S field we find the EFT: We note that there are corrections to the renormalizable |H| 4 vertex, which we will find is a common feature of integrating out scalars in our models, as well as the dimensionsix operators Q H and Q H which affect the trihiggs couplings. Additionally the term gg 3 HS /3/M 6 appears to be of the next order in the EFT expansion, we will retain these terms in the text, however in our summary tables 2 and 3 we neglect such corrections.

Complex scalar singlet
For the complex scalar singlet we consider the case of Y = 0. While the complex scalar singlet is technically the same as introducing two real singlets, and therefore doesn't fit our criteria for considered models, we consider it here as it has been studied extensively in the literature. Some examples from the literature which study the complex singlet case and its implications for inflation, the electroweak phase transition, enhancement of the di-Higgs signal, and vacuum stability include [6,[33][34][35][36]. The Lagrangian is then: The M term corrects the dimension-six operator coefficients with terms proportional to M /M which must be small for the validity of the EFT so we neglect them. 3 Integrating 3 M is the parameter which dictates the size of the mass splitting between the components of the complex scalar field. If M were to become large it is possible that the lighter resonances would enter the low energy spectrum and invalidate our EFT approach. Therefore it is a requirement of our EFT approach that this parameter be small. For the same reason we will neglect the effects of Y3 in the 2HDM below. out Φ and Φ † gives the effective Lagrangian:

JHEP05(2018)061
Again we induce corrections to the |H| 4 vertex as well as the effective operators Q H and Q H .
We see that the 2HDM only induces one purely bosonic operator, Q H , at leading order in Y 3 /M 2 , and induces various combinations of rescalings of the Yukawa couplings, i.e. the operators Q eH , Q uH , and Q dH . The only difference between the various realizations of the 2HDM considered are differences in the weight of the fermionic operators, i.e. by tan β or cot β. To make manifest the mass dependence of the Higgs couplings to fermions above we have expanded the fermionic dimension-six operators (in the unitary gauge for convenience) to recast the couplings of H 1 to fermions in terms of their masses, Z 6 , and the mixing angle β. In particular the first line of each expression indicates the shift of the Higgs-fermion couplings relative to the SM prediction, Another unique feature of the 2HDM effective Lagrangians is that they also contain 4-Fermi operators. These are not relevant to our analysis and, as they are weighted by the square of the Yukawa, are unlikely to have large Wilson coefficients except possibly in the case of the top quark which has Y t ∼ 1.

Real scalar triplet
The real scalar triplet model [18,41,42] has been studied in the literature with ambitions of making the electroweak phase transition first order, e.g. in [43], with the possibility of the neutral component being a dark matter candidate [9], as well as from an EFT point of view in [28,44].
The relevant Lagrangian is given by, Integrating out the heavy triplet then gives the effective Lagrangian: 14) It is convenient to make a change of basis here, we may exchange the operator |H| 2 (D µ H) † (D µ H) for the other dimension-six operators at the cost of an error of the next order in the EFT (i.e. O(1/Λ 4 )). While it is frequently simpler to maintain the basis obtained after integrating out the heavy states [45], for the sake of this work which will consider many UV completions and their effective field theories we choose to project onto a common basis. Discussions of the validity of this method including proofs of the invariance of the S-matrix can be found in [46][47][48][49]. We perform the change of basis by using the Higgs equation of motion, scaled up to dimension-six through multiplication by additional Higgs fields, where we have called the renormalized (H † H) 2 coupling, λ R = λ + g 2 /8/M 2 with λ the (H † H) 2 coupling of eq. (2.1), yielding the new form of eq. (2.14): Consistent with our other examples we have again generated the Q H and Q H operators, however interestingly we have also generated the Q HD operator which will have important phenomenological implications which we discuss in section 3.

Complex scalar triplet
Charging the Scalar Triplet under hypercharge, Y = −1, has important uses in the Type II seesaw [10][11][12][13]. The relevant UV complete Lagrangian is then, Integrating out the heavy complex triplet yields the effective Lagrangian, (2.18) which after applying the equation of motion from eq. (2.15) (notice here λ R = λ+|g 2 |/2/M 2 ) gives the final form for the effective Lagrangian: This effective Lagrangian and the effective operators it contains are consistent with our expectations from the other models, particularly the real scalar triplet.

Quadruplet with Y = 3Y H
For the two quadruplet models we follow the notation of [50], the UV Lagrangian is then given by: Integrating out the quadruplet leads to the simple EFT, Note that for a quadruplet we expect a contribution to the T -parameter. This operator does not occur at dimension-six, but does at dimension-eight. Deriving only the dimension-eight operator contributing to the T -parameter yields: Here we have confirmed the sign of [51]. We will see in the case of Y = Y H we obtain a different sign from this work.

JHEP05(2018)061
The UV complete Lagrangian is given by, Again we find a very simple EFT to dimension-six: which we supplement with the dimension-eight T -parameter operator.
This expression agrees with [51] up to a sign. As the sign of the dimension-eight T parameter operators in each quadruplet model come purely from the covariant derivative term of the Lagrangians (other contributions cancel) they should be the same in both eqs. (2.22) and (2.24). Unsurprisingly neither the 2HDM nor the two singlet models generate Q HD , also referred to as the T -parameter operator as they are known not to shift the relation between the Wand Z-masses. It is, however, expected from studies of the dynamics of the triplet models below EWSB that the triplet models considered in this work correct the T -parameter. This is consistent with our findings in Equations (2.16) and (2.19). In the case of the quadruplet we found they were unique in that at dimension-six they generate only one operator, Q H , and that the T -parameter operator was generated at dimension-eight. Additionally, as there are no allowed tree level couplings to Fermions in any of the theories except the 2HDM none of the other theories generate the fermionic operators, however after trading the operator (H † H)(D µ H) † (D µ H) in the triplet models via the EOM we do generate the fermionic operators for the two triplet models. The case of the quadruplets is particularly interesting as studies which indirectly probe the Higgs self coupling, such as [21], only allow the SM coupling λ to vary. Our work indicates that, within the assumptions of our EFT, 4 such a study corresponds to a very specific

JHEP05(2018)061
Theory: UV complete scenario, in the case where one expects the NP to come from dimension-six operators this corresponds to the quadruplets. In the case of the quadruplets the shift in λ due to the effective operators is restricted to be extremely small since the same UV parameter that generates the operator Q H contributes to the strongly constrained T -parameter. This demonstrates that indirect probes of the Higgs self coupling which don't vary other Higgs couplings are incomplete or correspond to specific UV completions which do not satisfy the criterion of the UV complete models considered. Other studies which vary these additional couplings of the Higgs such as [22,23] indicate the bounds on the Higgs self coupling are weakened or even lost without the inclusion of the direct di-Higgs probe. It is useful to project these effective Lagrangians into Lorentz forms relevant to the di-Higgs analysis performed. We do so here, from the perspective of arbitrary Wilson coefficients, when the final analyses are performed we use the expressions for the Wilson coefficients expressed in table 3. We assume that only the heaviest generation for each fermion has a non-negligible contribution to the EFT. Starting from the effective Lagrangian, we can proceed to expand the operators to find the relevant Lorentz forms. Here we have used λ RF to represent the final renormalized coefficient of the (H † H) 2 operator, the expression for λ RF may be found in table 2 in terms of λ of eq. (2.1) and the parameters of each UV-model. This involves finite field renormalizations as the operators Q H and Q HD both alter the Higgs kinetic term below EWSB. Details of this procedure may be found in, for example, [52][53][54]. Below EWSB expanding out the Lorentz forms we find (employing the unitary gauge):

JHEP05(2018)061
Theory: Table 3. Summary of the tree-level effective field theory to dimension-six for the scalar theories considered. "-" indicates the operator is not generated in this theory. The UV operators with normalizations corresponding to each coupling constant should be read directly from the relevant Lagrangians in text. In this table, as mentioned in the text in the Real Scalar singlet discussion, we neglect terms which are of O(g 4 /M 6 ). While the operator Q HD is not generated in the quadruplet models we have entered the contributions to the T parameter in terms of an effective coefficient for this operator into the table.
Here "· · · " indicates the various operators and Lorentz forms which have no impact on our analysis. The coefficients of the terms in the Lagrangian of eq. (2.27) are given by: Note in eq. (2.28) we have introduced m ψ and c ψH as placeholders for the relevant fermion type (i.e. e, u, or d), and in this analysis we only consider couplings to the third generation of each. We have only included g HHu and its corresponding operator as only the top quark h 2ψ ψ operator will have an effect on our analyses as it is proportional to the top-quark Yukawa coupling which is the only large Yukawa in the SM. It is possible to remove the JHEP05(2018)061 HHH operator by a field redefinition of h, however as pointed out in [55] removing this operator by a field redefinition of h (not the full doublet H) requires a nonlinear field redefinition which may prove to make one loop calculations difficult and if done incorrectly gauge dependent. Therefore we retain the g (2) HHH coupling in favor of easier comparison with other works, such as those which study globally the constraints on the h 3 coupling via one loop dependent processes [21,22,24,25,56,57].

Higgs coupling measurements at the LHC
In this section we consider important constraints on our EFTs in section 2. We begin by considering the constraints from electroweak precision data along with a discussion of the loop order at which the S-and T -operators are generated either explicitly via integrating out at the mass scale of the extended scalar sectors or via operator mixing in the EFT while running down to the Higgs mass scale. Next we introduce the effective hγγ coupling in order to add an additional constraint to our global fit to single Higgs processes. We delegate to appendix A unitarity considerations from the EFT perspective, where many amplitudes grow with the square of the center-of-mass energy S, as they do not add additional constraints to our models. Finally with our precision constraints on the EFTs we project these constraints into the UV complete models parameter spaces, this is especially useful in helping to limit the size of the c H coupling which is partially dependent on the same couplings as the hγγ effective coupling.

Electroweak precision measurements
Electroweak precision data (EWPD) provide very strong constraints on the Wilson coefficients of effective operators. We note that the operator Q HD contributes at tree level to the T -parameter, while the operator, contributes to the S-parameter at tree level. However, the only operators contributing to EWPD that are generated at tree-or one-loop level in our theories are Q H and Q HD the operator Q HW B is only generated at two-loop or higher order. From Jenkins et al. [54,58,59] we have the elements of the anomalous dimension matrix for each of these operators: × c HW B + 4g 1 g 2 y h c HB + 4g 1 g 2 y h c HW .

JHEP05(2018)061
Where we have introduced the U(1) Y , SU(2) L , and SU(3) C couplings g 1 , g 2 , and g 3 respectively, n g is the number of active generations at the relevant energy scale, the operators corresponding to the wilson coefficients c W , c HB and c HW are given by, and "· · · " represents other operators not generated at tree-level in our EFTs. The final line of eq. (3.2) is included to indicate that c HW B is not generated at 1-loop by operator mixing and therefore must be generated at two-or higher loop order. However, the T -parameter is generated at tree-level by the triplet models, and one-loop by any theory which induces c H (namely all but the 2HDM). In the quartet models, since the only dimension-six operator is the H 6 operator, there is no contribution to S and T from the H 6 operator. However, the T -parameter can be generated at tree-level by dimension-eight operators.
Including both the one-loop and running effects we have for the S and T parameters (see e.g. [60] and [52]): (3.5) Again we have used · · · to represent operators generated at higher loop order in our theories. For the quadruplet models, the dimension-eight operators generate the following T -parameter where we have defined the Wilson coefficient c T 8 to be the coefficient of the T -parameter operator at dimension-eight. This coefficient c T 8 is then given by, for the Y = Y H and Y = 3Y H quadruplet models respectively. Note that the coefficients c T 8 depend on the same quadruplet parameters as the operator H 6 , and therefore the Wilson coefficient of H 6 is also strongly constrained through this correlation. From GFitter [61] we have the central values of the S and T parameters with correlation matrix ρ as follows, When considering all of the operators discussed above one may perform a sophisticated fit to the EWPD of the many operator coefficients (see e.g. [62]), however for our study we need only consider c H and c HD as discussed above. Therefore performing a simplified chisquare fit relevant to our EFTs, we obtain constraints on the Wilson coefficients (c HD , c H ): We note that c HD is tightly constrained while c H is not as its contribution to S and T is generated at one-loop.

Higgs diphoton rate
In section 2, only the leading tree-level effective operators are written when integrating out the heavy scalars. The leading effective operators which contribute to the Higgs diphoton signature are not included in our framework as they originate from the one-loop contributions. However because of the precision of the H → γγ measurements we will include them in this section. Note that after integrating out the heavy scalars at one loop one may expect contributions to the H → γγ coupling from the following gauge-invariant dimension-six operators, However, since we are only interested in the diphoton rate, and not in corrections to the h → ZZ and h → W W rates we may simplify the calculation of the Wilson coefficients by only considering one effective operator in the broken phase: The general Higgs diphoton Wilson coefficient c γγ for new scalars and fermions at one loop may be found in, e.g. [63]. For the UV complete models considered in section 2 we find the wilson coefficients in table 4. As mentioned in the previous section the Wilson coefficients of the Quadruplet model are all proportional to the parameters contributing to the T -parameter. As such we will not consider the Quadruplet models for the rest of this section.
Finally the diphoton rate relevant to our models is,

JHEP05(2018)061
Where we have defined as the SM part of the h → γγ width taking into account shifts in the couplings of the Higgs to the t-quark and W -bosons due to the effective Lagrangian of eq. (2.27). Here the loop functions A 1/2 (τ ) and A 1 (τ ) are defined in ref. [63].

Higgs global fits
The Run-I Higgs measurements [64][65][66][67] provide constraints on some Wilson coefficients in the effective Lagrangian. For convenience we reproduce our effective Lagrangian below EWSB here: 14) The corresponding Wilson coefficient dependence can be found in eq. (2.28) while the Wilson coefficients for each model can be found in tables 2, 3, and 4. We note that the modified Yukawa coupling of the top-quark also causes a shift the Higgs-digluon effective coupling which we have taken into account in our analyses. These Wilson coefficients contribute to the Higgs signal strengths µ = σ×A× [σ×A× ] SM extracted from the Higgs coupling data, where A × is the product of the acceptance and the efficiency. Since the Higgs discovery global fits to the effective operators relevant to Higgs physics have become an important area of research [52,68,69] and recently they have gone beyond simple inclusion of signal strengths to inclusion of kinematic variables and off-shell measurements [70,71]. They have also been considered in scenarios where EWSB is not linearly realized [72][73][74]. However for the sake of our analyses we require a much smaller set of effective operators, therefore we perform a simplified global fit to the Higgs signal strengths µ i using the program Lilith [75].
In Lilith, all the Run I LHC Higgs measurements [64][65][66][67] are taken into account, and a likelihood statistical procedure is performed to obtain the constraints on the signal strengths. It is based on the assumption that the Higgs measurements are approximately Guassian and thus the likelihood function L(µ) could be simply reconstructed. Under this assumption adapted by Lilith, the −2 log L(µ) follows a χ 2 law for each observable, whereμ i is the theoretical prediction of the measured Higgs signal strengths µ i with Gaussian uncertainty ∆µ i . The full likelihood L(µ) = i L(µ i ) is defined as where C −1 is the inverse of the n × n covariance matrix, with Then the constraints on the signal strengths are recast as bounds on the Wilson coefficients. We perform a global fit on these Wilson coefficients (c HD , c H , c γγ , c iH ) with i = t, b, τ , and then project our results into the sub-space in each scalar model. First we perform the six-parameter fit, and obtain where ρ is the correlation matrix for this global fit. These Wilson coefficients are typically small due to suppression by v 2 M 2 . However from subsection 3.1 we know we must also consider the EWSB constraints. Assuming equal weight and combining with the constraints coming from the S and T parameters, we find that C HD is very tightly constrained: (3.18) We also obtain that v 2 * (c HD − 4c H ) = −0.09256 ± 0.8731, which by eq. (2.28) we see is a very important constraint on both the momentum dependent and momentum independent tri-higgs couplings. In figure 1 we show the v 2 (c HD − 4c H ) × c scalar γγ plane where we have marginalized over the parameters not shown.
We see from figure 1 that the independent constraint on c γγ provides an important constraint in the space of Wilson coefficients which will translate to a constraint on the JHEP05(2018)061 Table 5. The central values and 1σ errors of the model parameters for each UV complete model. We also limit the range of the dimensionless Higgs couplings to be less than ±4π.
various four scalar couplings of the UV models and therefore through their correlation with the Wilson coefficient c H on the affects of the Q H operator. We project these constraints in the EFT framework onto the UV complete model parameters in the next subsection.

Implications for the UV physics
In the global fitting procedure, all the Wilson coefficients are assumed to be independent. We know from section 2 that in the specific scalar extended models some Wilson coefficients are correlated and some Wilson coefficients may be absent altogether. These correlations and absences may be seen in table 3. Therefore, it proves useful to recast the global fit results to obtain constraints on the UV model parameters in each model. We perform the global fit using the Lilith program in each scalar extended model. In figure 2, we show the 1, 2, and 3σ contours on the model parameters in the real and complex singlet, Type-I doublet, and complex/real triplet models. At the same time, we also show the central values and errors for the model parameters in table 5. These plots exhibit similar features. First, the Higgs-Higgs-scalar coupling g/M 2 or Z 6 /M 2 is constrained to be O(0.1 − 1) by the Higgs gauge boson couplings in the singlet and doublet models, while in the triplet models the T -parameter puts tighter constraints on the parameter g/M 2 . Secondly, for the doublet and triplets, the Higgs to diphoton rate puts additional constraints on the couplings which contribute to the c γγ . Converting to the couplings in the UV model, we are not further able to constrain the Higgs-Higgs-scalar-scalar couplings of the triplet models λ HΦ and λ , because the constraints shown in figure 2 and table 5 are very loose. Even the perturbativity constraint, shown as the blue dashed lines in figure 2, is tighter than the constraint from the global fit. So to place constraints on the Wilson coefficients of Q H for the 2HDM and triplet models, we have to rely on di-Higgs collider constraints. Finally, we note that although the global fit cannot constrain the renormalizable Higgs self coupling λ, it is able to constrain the dependence of the h(∂h) 2 effective coupling indirectly. We have neglected to project our global fit into the parameter space of the quadruplet as it is so strongly constrained by the T -parameter and the triplet serves as an example of the affects.
While these indirect constraints on the UV models from the global fit are interesting and useful for our di-Higgs analysis in the following section, stronger constraints may of course be found in UV complete considerations of these models. The ability to loosely constrain numerous models at once from simple Higgs global fits is nonetheless intriguing and (especially in the advent of a significant deviation from the SM expectation) a useful way to direct UV complete searches of greater depth in the future.

Di-Higgs production at the 100 TeV collider
The measurement of the triple Higgs coupling using non-resonant di-Higgs production at both the LHC and future 100 TeV collider has been studied in great detail in the literature which was recently reviewed in [27]. Among all the channels for the Higgs decay final state, the bbγγ channel [76][77][78][79][80][81][82][83] is the most promising due to the combination of large h → bb branching ratio and more accurate reconstruction of photon momentum compared with other channels which helps reduce the backgrounds.
Three different topologies of Feynman diagrams of the pp → hh process via the gluon fusion production are shown in figure 3. Due to the destructive interference between the triangle and box diagram for the di-Higgs production in the gluon fusion channel, it is believed that at 14 TeV LHC with 3 ab −1 luminosity, the triple Higgs coupling −g (1) HHH /λ SM v

JHEP05(2018)061
would be constrained to only [−0.8, 7.7] at 95% CL [84]. In all models considered in this article, the Wilson coefficients of the |H| 6 operator cannot be chosen arbitrarily large. Based on the considerations of the validity of EFT and perturbative constraints, we estimate the value of the modified trilinear Higgs coupling to be within the range (−0.1λ SM , 2λ SM ), and take the cutoff scale to be 2 TeV. The higher the cutoff scale, we expect the narrower range of the trilinear Higgs coupling. On the other hand, at a 100 TeV collider with 30 ab −1 luminosity, the SM value of the triple Higgs coupling can be measured with around 10% uncertainty [27], and even around 4% based on the latest study [85]. Therefore, we expect that 100 TeV collider provides a good opportunity to explore the Wilson coefficients c H in various models we have considered. 5

General formalism on di-Higgs production
In our EFT framework, the effective Lagrangian relevant to the di-Higgs production is where g (1) with the SM vacuum expectation value v ≡ 1 2( √ 2G F ) 1/2 and the SM dimensionless coupling λ SM ≡ √ 2G F m 2 H . From the above Lagrangian, we note that in the Warsaw basis, in addition to the SM trihiggs couplings, we also have derivative triple-Higgs couplings, which may contribute differently to the distribution compared with solely non-derivative couplings.
According to figure 3, the parton amplitude of the di-Higgs production g(p 1 )g(p 2 ) → h(p 3 )h(p 4 ) via the gluon fusion process is HHH vŝ  where the Lorentz structures are and F , F , and G are the form factors for triangle and box diagrams which can be found in ref. [86]. Correspondingly, the differential cross-section for di-Higgs production is given by: where S is the center-of-mass energy squared of the proton-proton system,ŝ = (p 1 + p 2 ) 2 , t = (p 1 − p 3 ) 2 and the parton luminosity function is defined as with f g/p the gluon distribution function, and µ F the factorization scale. As we have previously noted, the triangle diagram and box diagram interfere destructively and the smallest cross section is obtained when g HHH /v ≈ −2.5λ SM assuming no derivative interaction and no corrections to the quark-Higgs couplings. Due to this fact, the variation in the gluon fusion to di-Higgs cross section about the SM value of g (1) HHH decreases, the total cross section decreases, till g (1) HHH reaches −2.5λ SM . Any further decrease in g (1) HHH results in increasing of the cross section with respect to its minimum value at g (1) HHH /v ≈ −2.5λ SM eventually surpassing the SM value for g (1) HHH values lower than −5λ SM . On the other hand as g (1) HHH increases from zero, the total cross section increases. In our case, the situation is more complicated, we now have both an additional vertex and corrections to the quark Higgs couplings.

Di-Higgs cross section
In figure 4 we show the cross section contours of the pp → hh process in the (g (1) HHH /v, g (2) HHH v) plane with three different values of c tH . To evaluate the range of trihiggs couplings g (1) HHH /v and g (2) HHH v, we first use the eq. (2.28) and table 3 to express the two couplings in terms of the parameters in the UV model, then varies the dimensionless parameters in the UV models within the range ±4π, couplings with mass dimension to be in the range ±1 TeV, and the cutoff scale are set to be 2 TeV. These values are chosen such that our EFT matching procedure is valid (dimension-eight operators will not be enhanced by the factor g 2 /M 2 ) and the contribution of the kinematic region larger than cutoff scale to the total rate is negligible due to the suppression of the parton luminosity. After these consideration, we choose relatively loose ranges for the two couplings: g (1) HHH ⊂ (−0.36, 0.07) and g (2) HHH ⊂ (−0.015, 0.015). For c tH = 0, the anomalous Higgs fermion coupling g HHt in eq. (4.5) vanishes and the corrections to the quark Higgs couplings are proportional to c HD − 4c H . In such a case, only the first triangle and box diagrams of figure 3 contribute to the cross section with approximate SM quark Higgs couplings. Hence, one can find that, along the positive vertical direction, given a fixed value of g (2) HHH , the cross section increases. Along the g (2) HHH direction, one can find that a positively increasing value of g (2) HHH will lead to an increase in the total cross-section. This can be understood from eq. (4.6), where we observe that, with a positive g (2) HHH , the second term inside the bracket in front of the F which is induced by the derivative interaction will add destructively with the first term which is induced by the ordinary triple Higgs interaction, such that the effect of destructive interference between the box and triangle diagrams is alleviated.
In the case of c tH = 0.4, the cross section increases significantly when compared with the cross section for c tH = 0, this can also be understood from eq. (4.6) and eq. (4.4): the positive c tH will decrease the magnitude of g tH and also gives a new positive term generated by tthh vertex, which will alleviate the destructive interference. In the case of JHEP05(2018)061 c tH = −0.4, the cross section will reach some minimum value between g (1) HHH /v = −0.1 and −0.15 due to the destructive interference. Below the miminum points, for a fixed g (1) HHH , increasing g (2) HHH will decrease the cross section, because at this point the amplitude from the triangle diagram becomes dominant, increasing g (2) HHH will decrease the magnitude of the term inside the bracket in front of the F , thereby decreasing the cross section.

Monte Carlo simulation and validation
In order to perform our simulations we begin by using FeynRules [87] to generate an UFO model file adding the effects of the dimension-six operators in eq. (4.1). We then modify the model file to include the full triangle and box form factors as computed in [88]. Then we implement MadGraph 5.2.4.3 [89] to generate events. We use Pythia 6 [90] for the parton shower and the FCC card in Delphes 3.4 [91] for simulating the detector. The following analysis is only concerned with statistical uncertainties as the systematical uncertainties are unknown at the moment. When taken into account they will lower the significance levels given in this section.
We refer to the cuts applied while generating the events in MadGraph/Delphes as preselection cuts in the table 6. They are as follows: 6 Important irreducible backgrounds consist of Z(bb)h(γγ), tth(γγ), bbh(γγ), bbγγ production. Apart from these, there are bbjγ, jjγγ, ccγγ and bbjj channel that can potentially have a contribution to the background. Jet fake rates to photons are taken to be 0.012%, while jet and charm mistagging rates to bottom quarks are taken to be 1% and 10% respectively [92]. The backgrounds can be greatly reduced by vetoing extra jets, i.e., by demanding exact two b-tagged jets in each event. This is particularly helpful in reducing the tth background. Applying a Higgs mass window cut of 112.5 < m bb < 137.5 GeV, to the invariant mass of b-jets results in a large reduction in the Zh background due to exclusion of the Z-peak region.
The Higgs mass window cut for the di-photon invariant mass is sharper than that for the invariant mass of b-jets and helps to reduce the background in all the channels. Furthermore, from the normalized distributions for b-jet-pair p T and di-photon p T in figure 5 indicate that the signal is favored for p T values larger than 150 GeV and 140 GeV respectively. Therefore, we further apply these cuts in order to enhance the statistical significance. The resulting efficiencies and cross sections at each stage due to these cuts in our analysis for leading backgrounds and three benchmark (BM) points for the signal are tabulated in table 6.
We first investigate the sensitivity of the trilinear Higgs coupling λ HHH = −g (1) HHH /v in the absence of the derivative Higgs coupling g (2) HHH . In this case, we recover the scenario widely discussed in the literature: how to probe the deviation of the λ HHH from its SM HHH v) plane are as follows: BM1=(0.0225, 0), BM2=(−0.032, 0.0152), and BM3=(−0.141, 0.0152). HHH /v assuming that the derivative Higgs coupling g HHH is zero. The orange and green bands correspond to the 1σ uncertainty in the S/ √ B with assumptions of the theoretical uncertainty for the di-Higgs production cross-section to be 4% and 10% respectively. Right panel: the percentage uncertainties on the measured number of signal events varies with the value of trilinear Higgs coupling. Orange and green lines correspond to theoretical uncertainties of 4% and 10% respectively.
value λ SM at the future collider. Compared with the work in [80], we obtain comparable significance of about 8.25σ for the SM di-Higgs production for luminosity of 3 ab −1 . This corresponds to the significance of ∼ 26σ for 30 ab −1 as can be seen from the black line in the left panel of figure 6, where we plot the S/ √ B for 30 ab −1 and zero derivative interaction. We also estimate the uncertainty in the value of S/ √ B by taking into account the statistical uncertainty for the signal and background as well as the theoretical uncertainty on the di-Higgs production cross-section. It turns out that for a 30 ab −1 luminosity, the statistical uncertainty in the number of signal events due to Poisson fluctuations is around 3%, which is less than the 10% theoretical uncertainty coming from the infinite top mass approximation, the scale, and the PDF uncertainties [27]. The 1σ uncertainty due to this is denoted by the green band in the left panel of figure 6. However, the latest estimation on the theoretical uncertainties places them as low as 4% [85]. Therefore, we also include this case denoted by orange band in the plot shown in the left panel of figure 6.
The right panel of figure 6 represents the percentage uncertainties for the measured number of signal events as a function of the ratio of the triple Higgs coupling to its SM predicted value. Orange and green lines here correspond to the theoretical uncertainty of 4% and 10% respectively. As expected from the above quoted numbers, the theoretical uncertainty dominates except where the ratio of triple Higgs couplings is close to 2.5, where the cross section for di-Higgs production is the lowest leading to enhanced uncertainty due to Poisson fluctuations.
Here we comment on the validity of EFT in our collider analysis. The EFT breaks down when the parton collision center of mass energy approaches the scale of the cutoff scale M = 2 TeV. Therefore, we should in principle add a cut on the kinematic variables JHEP05(2018)061 2.5% 0.38% 0.16% 2 5.1% 1.0% 0.35% -0.9 1.3% 0.26% 0.05% Table 7. The percentage of events with m hh above 1, 1.5 and 2 TeV.
like invariant mass of di-Higgs to only keep the events produced in low energy regime to make our EFT analysis valid. The di-Higgs spectrum is peaked at an invariant mass m hh near the two higgs threshold indicating our EFT approach should be valid (i.e. the processes considered have energy well below our cutoff of 2 TeV). Additionally we have investigated the number of events below 1 TeV, 1.5 TeV, and 2 TeV for three benchmark points: the SM, λ HHH /λ SM = 2, λ HHH /λ SM = −0.9, and finding the results in table 7. As there are only a small number of outlying events with higher energies these numbers support the assertion that the EFT approach is valid in our Monte Carlo simulation. One should note that even if the heavy particles were to be discovered at higher energies that in order to extract the trilinear couplings of the SM Higgs one would still employ an EFT. Such a procedure is analogous to the use of an effective four fermion theory for flavor physics where the heavy W s have been integrated out of the theory in favor of unrenormalizable operators.

Determination of Wilson coefficients
Equation (2.28) and table 3 demonstrate it is necessary to investigate the discovery potential at the 100 TeV collider when both the deviation of the λ HHH coupling from the SM value, and non-zero g HHH exist. Turning on the derivative Higgs coupling g HHH will change the significance of the di-Higgs signatures. In figure 7 we present the reach of the 100 TeV collider with integrated luminosity of 30 ab −1 in the space of g (1) HHH − g (2) HHH in the left panel as well as in the space of Wilson coefficients c H and c HD − 4c H in the right panel, each with c tH = 0. The left and right panels of the figure 7 are not independent. Their values are connected by eq. (2.28), where the contours in the right panel are essentially rotated around the SM values as governed by the eq. (2.28). This represents only a class of models, in which c tH is not important, for example, singlet, triplet and quadruplet models. We plot the statistical significance contours for 2HDMs in c tH − c H space as shown in separate plots of figure 8. c tH = 0 corresponds to tan β → ∞, which is outside the experimental bounds on tan β in 2HDMs. Figure 7 shows the allowed parameter regions in singlet, triplet and quadruplet models, which overlap within the significance contours. In these models, according to table 3, the Wilson coefficients c H and c HD − 4 c H are not independent. More specifically, they are related by linear relations such as c H λ HS(Φ) (c HD − 4 c H ). This linear relation then implies that the boundaries of these regions are governed by the input perturbative limit |λ HS(Φ) | ≤ 4π and are straight lines as can be seen in figure 7. The values of the dimensionless Higgs scalar couplings, such as λ HS , λ HΦ , determine the slopes of the  parameter region in each model. For example, in the real singlet case, along the boundary of the parameter region, the Higgs scalar coupling λ HS should be around ±4π. In the region far from the boundary, the dimensionless Higgs scalar couplings appearing in c H should be small. We choose c tH to be equal to zero in these two plots. This condition is automatically satisfied by singlet and quadruplet models, and also approximately satisfied by triplet models. This is because c tH in triplet models is suppressed by the coupling g 2 which is constrained to be very small by EWPD due to its relation to the T -parameter.
In addition to the allowed region in each model, we also illustrate the region that will generate the expected significance within the 2σ uncertainties around SM value. In principle one should derive the prospective confidence level contour for the parameter space that are consistent with SM prediction in the future experiments. However, this requires the scanning of a fine grid to obtain the selection efficiency of each point, which is beyond the scope of our study. Therefore we simply estimate that this 2σ region roughly gives the region that is hard to differentiate from the SM in the future experiments.  One can observe that, the future di-Higgs experiment is not sensitive to the c H and c HD which have already been strongly constraint by the EWPD. On the other hand, it can constrain the value of c H . Depending on the theoretical uncertainties that can be achieved, it may also be possible to exclude some parameter space of the singlet models, which represents the region outside the 2σ region.
The case of the 2HDM is much more promising for distinguishing between the SM and the NP model as c tH is non-zero. We demonstrate the significance for 100 TeV collider at 30 ab −1 integrated luminosity in v 2 c tH vs v 2 c H plane, shown in figure 8. Both v 2 c tH vs v 2 c H depend on tan β. Here we choose the range of tan β such that it satisfies the constraints from flavor physics according to ref. [93]. This rules out some parameter regions as shown in figure 8 by blue regions. We note that the significance in the 2HDM is generally larger than that of the singlet, triplet and quadruplet models due to typical enhancement from the Yukawa couplings, and it is very likely to observe a significant deviation from the SM signal. We also find that, unlike the singlet and triplet, signal significances in the 2HDM are much more enhanced compared to the ones in the SM. The plots also show that the contours of significance of two types of 2HDMs are different despite the coupling to up-type quarks being the same in both Type I and Type II, the reason being that we are using the bbγγ final state and the branching ratio of h → bb are different between the two versions of the 2HDM. From figure 7 and figure 8, if we limit ourselves in these models with all the heavy particles integrated out, the di-Higgs process puts additional constraints on the scalar model  parameters. Our analysis in figure 7 and figure 8 shows that the Complex singlet and 2HDM (triplet and quadruplet) scalar models are the most (least) sensitive, among those resulting from the models under consideration, to the collider search. As a consequence, the di-Higgs process probes the allowed region of c H , and thus the Higgs scalar couplings in the UV models.

Exploring parameter region in UV models
We project the sensitivity of the Wilson coefficients into the parameter space corresponding to the models under consideration. In the real singlet model, the parameter space of the effective coefficients allowed is indicated by the light blue region in figure 7, can be probed with S/ √ B more than 25, while in the complex singlet model, the Wilson coefficients resulting from integrating out the complex singlet can be probed to S/ √ B values higher than even 40. In figure 9, we show the possible reach of the model parameters (λ HS , g HS ) in the real singlet model, and (λ HS + λ HS /2, g HS ) in the complex singlet model, given 30 ab −1 luminosity data set. One can see that, most of the region in the singlet and triplet models are within the 1σ uncertainty band for S/ √ B reach for the SM, so that they are hard to differentiate from the SM.
The 2HDM, owing to its preservation of custodial symmetry, resides on the line c HD = 4c H = 0 (up to the assumptions made in this paper, that is a tree-level dimensionsix analysis). Therefore, the Higgs coupling measurements and the electroweak precision tests do not place strong constraints on the model parameters. On the other hand, the  di-Higgs signature starts to provide a strong constraint on c H . In figure 10 we show the significance contour on the model parameter Z 6 vs tan β plane for Type-I model and Type-II model with the 30 ab −1 luminosity. Note that when Z 6 = 0, the SM limit is recovered (see table 3). We also find that in the Type-II model, for negative Z 6 and large tanβ (left top corner in the right plot in figure 10), the significance approaches to the SM value. This is because the decreasing of the Higgs to b quark coupling reduces the Higgs to b decay branching ratio, which ameliorate the increasing of the di-Higgs production rate.
In the real and complex triplet models, both c HD and c H in the EFTs obtained by integrating out real and complex triplet models are very tightly constrained as shown by the vertical dashed lines, shown in figure 7 (right panel). These vertical darker and lighter green lines represent the 3σ bounds allowed by the Higgs data global fit on the Wilson coefficient linear combination of c HD − 4c H for the real and complex triplet model respectively. The reason that these stringent bounds only exist for the triplets and not the singlets is that the coefficient c HD is connected with custodial symmetry breaking and is tightly constrained by the electroweak precision parameter T .
As table 3 denotes, the c HD and c H are tightly related for the triplet models and therefore the stringent bounds on c HD translate into stringent bounds on the c HD − 4c H as well. In the case of the singlet models, there are no couplings of the singlets to the gauge bosons resulting in c HD being identically zero as indicated in table 3, liberating them from these constraints suffered by the triplet models. As a result of these, c H is also strongly constrained from the small allowed values of c HD −4c H , as shown in figure 7 (right panel) S/ B Quadruplet Y=3/2 30 ab -1 Figure 11. The discovery potential of the model parameters (λ H3Φ,M ) in the quadruplet. The dashed black contours correspond to the S/ √ B values for an integrated luminosity of 30 ab −1 . The blue region is excluded by constraints from the electroweak precision tests. The orange and green regions are within the SM 2σ uncertainty with an assumption of the percentage theoretical uncertainty of di-Higgs production cross-section equal to 4% and 10% respectively. loosely constrained due to c H ∼ g 2 M 4 λ HΦ . Therefore, it is very hard for us to extract the Higgs scalar couplings from the c H operator, because the deviation of the Higgs coupling from the SM value is very small in the triplet case.
For the quadruplet model, at dimension-six, only the Wilson coefficient of Q H operator is non zero. However, we include the c HD generated by dimension 8 operator because it is strongly constraint by EWPD. In the left plot in figure 11, the allowed region for two types of quadruplet models are denoted by two lines with different slopes. The reason can be seen from table 3, the c H and c HD are correlated, all proportional to the coupling |λ Φ3H | 2 . So given a fixed cut off scale M , both g (1) HHH and g (2) HHH can be parameterized by a single parameter |λ Φ3H | 2 . In the right plot in figure 11, we find that the allowed parameter space from the global fit to EWPD for quadruplet models is tightly constrained, and almost becomes a point near the SM value. In figure 9, we show the significance of the model parameter λ Φ3H vs new physics scale M varies with the 30 ab −1 in contours, while the blue region is excluded by the constraint on c HD from EWPD. One could observe that, the T-parameter constraint on λ φ3H is very sensitive to the cutoff scale, the reason is that the c HD is generated by dimension-eight operator so that it is proportional to the fourth power of (v/M ).
Our collider analysis demonstrates that the potential of the 100 TeV collider in probing the Wilson coefficients resulting from the five scenarios considered here is very promising with the 2HDM. The singlet, triplet and quadruplet models on the other hand are restricted due to electroweak precision measurements and their effective coefficients will have less sensitivity. These restrictions also manifest in the constraints on the deviation of the triple

JHEP05(2018)061
Higgs couplings in such models owing to the direct relation between c H and the triple and quadruplet Higgs coupling as shown in eq. (2.28).
An interesting consequence of our analysis is that, due to the difference in the allowed region for each model under the theoretical bound and the global fit constraints, it is possible to differentiate the 2HDM model from singlet, triplet and quadruplet models with the observation of a large deviation of the signal rate from the SM expectation. If a future experiment detects a significantly larger signal rate compared with the expected SM model value, then it should favor the presence of an extended scalar sector consisting of the 2HDM the assumptions of this work. If the future experiment does not detect a significant deviation from the SM expectation, then one may have hard time to differentiate SM from all the models considered here as well as models where the wilson coefficients are induced at loop level. Both a reduction in the theoretical uncertainty estimation and higher luminosities will be needed to make a more precise measurement of the di-Higgs signal rate.

Conclusions
We began by motivating a study of the dimension-six Higgs self-interaction operator Q H = (H † H) 3 in the Standard Model effective field theory because of its importance to studying the nature of electroweak symmetry breaking and the nature of the electroweak phase transition. We noted that the largest Wilson coefficients can be obtained by considering extended scalar sectors which are the only models which admit a tree level Q H operator. After identifying all possible SU(2) L representations which allow for a tree level Q H along with the corresponding hypercharge Y we wrote the ultraviolet complete Lagrangians for each model. Finally, assuming that the new scalars are heavy we integrated them out of the theory obtaining the dimension-six effective Lagrangian at tree-level. Of the seven models which generate Q H at tree level all but two generate more than one effective operators. Those which generate only Q H at dimension-six are plagued by strong constraints coming from the dimension-eight T -parameter operator. This helps put into context that any study performed by shifting a single parameter of the model is not model independent and, in the case of shifts in the self coupling of the Higgs coming from UV physics generating dimension-six operators as leading effects, not well justified.
After identifying the full set of tree-level dimension-six operators for the extended scalar sectors we proceeded to consider the constraints on the effective theories given single Higgs measurements from run I of the LHC as well as the electroweak oblique parameters S and T . In order to fully take advantage of the single Higgs measurements we also derived the Wilson coefficient for the effective Higgs coupling to photons although it enters at loop level after integrating out the new heavy charged scalars. For the constraints from the single Higgs measurements we implemented the tool Lilith to perform a global fit to the Wilson coefficients other than the Q H . Having constrained the Wilson coefficients we then projected those constraints back into the parameters of the ultraviolet models deriving relations between various parameters of the models.
It is the multi-Higgs measurement instead of the single Higgs measurement which determines the Wilson coefficient of the Higgs self-interaction operator Q H . As such, we

JHEP05(2018)061
have investigated the dependence of the coefficient c H on the di-Higgs production cross section, and studied simulations of the di-Higgs process for the proposed future 100 TeV collider. We then obtained the sensitivity contours of the Wilson coefficients in the general effective theory framework as the luminosity varies at the future 100 TeV collider. Finally, we reduced the g (1) HHH − g (2) HHH plane, for various cases of the UV complete extended scalar sector models, to its subspace for the cut-off scale of 2 TeV and in the perturbative regime for dimensionless Higgs couplings and demonstrated that most of these regions can be probed to the statistical significance of more than 5σ using the di-Higgs signatures in a future 100 TeV collider.
We have converted the discovery reach of the Q H operator into the Higgs potential parameters in seven UV models. Among the models considered, the Higgs self coupling in the singlet and doublet models could have large deviation from the standard model prediction, while the triplet and quadruplet models can only have very small deviation, due to the strong correlation between the T parameter and the Wilson coefficient of the |H| 6 operator. We showed that for the projected data collected for an integrated luminosity of 30 ab −1 at the proposed 100 TeV collider, the trilinear Higgs coupling in the all scalar models with a single heavy scalar integrated out could be fully explored. If a significant deviation in the trilinear Higgs coupling is observed, it will rule out the possibilities of the triplet and quadruplet models. On the other hand, if there is only small or no deviation, it will strongly constrain the Higgs potential parameters in the singlet and doublet models. Therefore, combined with electroweak precision data, the di-Higgs search can effectively differentiate singlet and 2HDM models from triplet and quadruplet models, within the framework of effective field theory of new scalar models.
Overall, the di-Higgs process provides a unique opportunity to probe the c H operators which cannot be obtained by single Higgs phenomena and the electroweak precision tests. These future experimental measurements on the Higgs self-interaction operator Q H will provide important information on the shape of the scalar potential under the assumption that the new scalars are very heavy. This provides a method complementary to direct phenomenological searches to find evidence of additional scalars in these extended scalar models.

JHEP05(2018)061
where we may freely substitute V iλ i → h for the Higgs boson. Considering only amplitudes which grow with the square of the center of mass energy S in the above cited works the authors found that the operator Q H is not bounded by unitarity considerations for 2 → 2 scattering, and that the operators Q H and Q HD only result in unitarity violation for the purely longitudinal case. Note that as the 2HDM does not generate either Q H or Q HD at leading order in the Y 3 expansion it does not generate operators which violate unitarity with growing S. It was found that for one operator non-zero at a time the bounds were given by, It should be noted these constraints indicate the largest allowed values of the two operator coefficients allowing for cancellations between them and not that all values within these bounds will be simultaneously allowed. We may then consider the largest √ S at which our EFT is valid (i.e. perturbatively unitary): For di-Higgs processes we consider √ S up to 1 TeV, therefore these bounds indicate we should not consider c H or c HD larger than (8/TeV) 2 , this bound is well outside the perturbative region of the UV models considered. For example in the case of the real scalar singlet, implies, for M ∼ 1 TeV, that g 11 TeV. For the largest allowed values of g scattering in the UV theory is not perturbative unless √ S > g, but for √ S > g the EFT approach we adopt in this study breaks down.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.