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-6 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=3Y_H$ and $Y=Y_H$, 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 dihiggs 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 $3$ ab$^{-1}$ luminosity data we will be able to explore the Higgs potential parameters in all seven models.


I. 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. 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 tri-Higgs 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 the higher dimensional 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 . These scalars should carry some SU (2) L charge which allows for a tree-level Q H along with a corresponding hypercharge Y . 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 dihiggs production process in a model-independent and predictive way.
Many new physics models with SM-compatible single Higgs phenomena could exhibit dihiggs 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 dihiggs process at the LHC and future colliders is the only direct way to measure the Wilson coefficient of the effective Q H operator. The dihiggs 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 dihiggs 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 dihiggs cross section may be increased considerably making measurement a possibility at the proposed 100 TeV collider [21,22]. In this paper we investigate the dihiggs 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.
We proceed by, in Section II, 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 III 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 IV we study these EFTs' impacts on dihiggs production at the proposed 100 TeV collider. Finally our conclusions are found in Section V.

II. THE EFFECTIVE LAGRANGIAN
We consider all ultraviolet (UV) complete models which include one additional heavy scalar charged under SU (2) L which generate, after integrating out the new scalar, dimension-six operators affecting the tri-Higgs 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 hypercharge Y = 0 and Y = −1 respectively, and finally complex scalar quadruplets 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. [23,24]. 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: 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 [25]. The operators which are relevant to our analyses are: 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 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.

A. 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 [23,26,27]. The Lagrangian, neglecting SM terms, is given by: L U D  Type I: Φ2 Φ2 Φ2  Type II: Φ1 Φ2 Φ1  Lepton-Specific: Φ1 Φ2 Φ2  Flipped: Φ2 Φ2 Φ1   TABLE I. List of Fermion couplings used for various Types of 2HDM. 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 dimension-six operators Q H and Q H which affect the tri-Higgs 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 II and III we neglect such corrections.

B. Complex Scalar Singlet
For the complex scalar singlet [6] we consider the case of Y = 0, 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. Integrating out Φ and Φ † gives the effective Lagrangian: Again we induce corrections to the |H| 4 vertex as well as the effective operators Q H and Q H .

C. Two Higgs Doublet Model
Of the many extended scalar sectors studied in the literature the two Higgs doublet model is the most well studied, reviews on the status of the model from the UV perspective have a long history (some extensive reviews include [15][16][17]), the two Higgs doublet model has also recently been studied in the EFT framework in great detail [27][28][29] including comparisons between the phenomenological aspects of both the UV complete and EFT frameworks at tree and one-loop levels [30,31].
We begin in the "Higgs basis", where the doublets have already been rotated to a basis where the physical CP even state is the observed 125 GeV Higgs. This rotation is performed by rotation of H 1 and H 2 by the angle β. We follow the notation of [28]. Note the Yukawa couplings are entered generically and later will be recast in terms of each of the four "types" usually considered to evade flavor changing neutral currents when we write the EFT. These various types considered are outlined in Table I.
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 realization 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.

D. Real Scalar Triplet
The real scalar triplet model [18,32,33] has been studied in the literature with ambitions of making the electroweak phase transition first order, e.g. in [34], with the possibility of the neutral component being a dark matter candidate [9], as well as from an EFT point of view in [23,35].
The relevant Lagrangian is given by, Integrating out the heavy triplet then gives the effective Lagrangian: 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 [36], 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 [37][38][39][40]. 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. 1, yielding the new form of Eq. 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 III.

E. 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, which after applying the equation of motion from Eq. 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.
For the two quadruplet models we follow the notation of [41], the UV Lagrangian is then given by: Integrating out the quadruplet leads to the simple EFT, Noting that for a quadruplet we expect a contribution to the T -parameter. This operator does not occur at dimensionsix, but does at dimension-eight. Deriving only the dimension-eight operator contributing to the T -parameter yields: Here we have confirmed the sign of [42]. We will see in the case of Y = Y H we obtain a different sign from this work.
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 [42] 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. 22 and 24.

H. Summary of EFTs
Finally after deriving the corresponding EFTs for each model we may construct a table with the Wilson coefficients for each operator for each model considered. We summarize the renormalization of the (H † H) 2 term in Table II and the Wilson coefficients of the dimension-six operators in Table III. While it appears that of all the theories the 2HDM is the only which does not generate a correction to the renormalizable (H † H) 2 , this is a reflection of neglecting terms suppressed by Y 3 /M 2 , these corrections are generated first at O(Y 3 /M 2 ). 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 W -and 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 16 and 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 was particularly interesting as studies such as [43] only allow the SM coupling λ to vary. Our work indicates that such a study corresponds to a very specific 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 Type II: Lepton-Specific: Flipped: TABLE III. 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 QHD 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.
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. It is useful to project these effective Lagrangians into Lorentz forms relevant to the dihiggs 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 III. 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 II in terms of λ of Eq. 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, [44][45][46].
Below EWSB expanding out the Lorentz forms we find (employing the unitary gauge): 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. 27 are given by: Note in Eq. 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 g (2) HHH operator by a field redefinition of h, however as pointed out in [47] 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 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 [43,[48][49][50][51][52].

III. HIGGS COUPLING MEASUREMENTS AT THE LHC
In this section we consider important constraints on our EFTs in Section II. 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.

A. 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. [46,53,54] we have the elements of the anomalous dimension matrix for each of these operators: 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. 30 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. [55] and [44]): 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 [56] we have, 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. [57]), however for our study we need only consider c H and c HD as discussed above. Therefore performing a simplified chi-square 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.

B. Higgs Diphoton Rate
In Section II, 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. [58]. For the UV complete models considered in Section II we find the wilson coefficients in Table IV. 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, 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. 27. Here the loop functions A 1/2 (τ ) and A 1 (τ ) are defined in Ref. [58].

C. Higgs Global Fits
The Run-I Higgs measurements [59][60][61][62] provide constraints on some Wilson coefficients in the effective Lagrangian. For convenience we reproduce our effective Lagrangian below EWSB here: The corresponding Wilson coefficient dependence can be found in Eq. 28 while the Wilson coefficients for each model can be found in Tables II, III, and IV. 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 µ i extracted from the Higgs coupling data. Since the Higgs discovery global fits to the effective operators relevant to Higgs physics have become an important area of research [44,63,64] and recently they have gone beyond simple inclusion of signal strengths to inclusion of kinematic variables and off-shell measurements [65,66]. They have also been considered in scenarios where EWSB is not linearly realized [67][68][69]. 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 [70].
In Lilith, all the Run I LHC Higgs measurements [59][60][61][62] are taken into account, and a likelihood statistical procedure is performed to obtain the constraints on the signal strengths. 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 These Wilson coefficients are typically small due to suppression by v 2 M 2 . However from Subsection III A 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: We also obtain that v 2 * (c HD − 4c H ) = −0.09256 ± 0.8731, which by Eq. 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 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 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.  In the global fitting procedure, all the Wilson coefficients are assumed to be independent. We know from Section II 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 III. 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 Fig. 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, with the benchmark value M = 1 TeV assumed. At the same time, we also show the central values and errors for the model parameters in Tab. V. These plots exhibit similar features. First, the Higgs-Higgs-scalar coupling g/v or Z 6 is constrained to be O(1 − 10) 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/v. 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 Fig. 2 and Tab. V are very loose. Even the perturbativity constraint, shown as the blue dashed lines in Fig. 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 dihiggs 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 The central values and 1σ errors of the model parameters for each UV complete model. In all models, we use a benchmark value M = 1 TeV. We also limit the range of the dimensionless Higgs couplings to be less than ±4π. 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 dihiggs 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.

IV. DIHIGGS PRODUCTION AT THE 100 TEV COLLIDER
The measurement of the triple Higgs coupling using non-resonant dihiggs production at both the LHC and future 100 TeV collider has been studied in great detail in the literature which was recently reviewed in [22]. Among all the channels for the Higgs decay final state, the bbγγ channel [71][72][73][74][75][76][77][78] 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 Fig. 3. Due to the destructive interference between the triangle and box diagram for the dihiggs 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 can be measured to only 30% uncertainty given that it is consistent with the SM expectation [71]. However, at a 100 TeV collider with 3 ab −1 luminosity, the SM value of the triple Higgs coupling can be discovered with a significance of S/ √ B = 8 [75] with a very well constructed di-photon invariant mass. In our EFT framework, the effective Lagrangian relevant to the dihiggs 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 tri-Higgs couplings, we also have derivative triple-Higgs couplings, which may contribute differently to the distribution compared with solely nonderivative couplings. According to Fig. 3, the parton amplitude of the dihiggs production g(p 1 )g(p 2 ) → h(p 3 )h(p 4 ) via the gluon fusion process is 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. [79]. Correspondingly, the differential cross-section for dihiggs 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 dihiggs cross section about the SM value of g 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. Therefore, 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 . For c tH = 0, the anomalous Higgs fermion coupling g HHt in Eq. 49 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. 50, 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. 50 and Eq. 48: 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 c tH = −0.4, the cross section will reach some minimum value between g  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.
In order to perform our simulations we begin by using FeynRules [80] to generate an UFO model file adding the effects of the dimension-six operators in Eq. 45. We then modify the model file to include the full triangle and box form factors as computed in [81]. Then we implement MadGraph 5.2.4.3 [82] to generate events. We use Pythia 6 [83] for the parton shower and the FCC card in Delphes 3.4 [84] 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. HHH , g HHH ) plane, the plots from left to right correspond to three different value of ctH = 0, 0.4, −0.4. We adopt the NNLL matched NNLO SM dihiggs cross section: 1.75 pb [22].
We refer to the cuts applied while generating the events in MadGraph/Delphes as preselection cuts in the Table VI. They are as follows 3 : 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 [85]. 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 HHH v) plane are as follows: p T and di-photon p T in Fig. 5 indicate that the signal is favored for p T values larger than 150 GeV and 140 GeV respectively. We expect such cuts to help enhance signal to background ratio due to the fact that the b-jets and photons belonging to the signal come from a decay of heavy particles like the Higgs. 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 VI. Using the above methodology, we scan the g (1) HHH space and the corresponding space in terms of the Wilson 0.0 0.5 1.0 1.5 2.0 5.
FIG. 6. The significance of the dihiggs process as a function of the tri-linear Higgs coupling λHHH = −g HHH /v assuming that the derivative Higgs coupling g (2) HHH is zero. The vertical axis on the left hand side is the result of 3 ab −1 integrated luminosity, while the vertical axis on the right hand side corresponds to integrated luminosity of 0.3 ab −1 . Statistical significance will be lowered due to systematics.
coefficients c H and c HD − 4c H . The range we scan is obtained by a naive estimation using the following method. First, we use the Eq. 28 and Table III to express the Higgs and scalar quartic coupling in terms of g (1) HHH and g (2) HHH and the remaining parameters in the UV model. Then for each model the constraint is obtained by demanding all the dimensionless parameters of the UV model to be within the range ±4π and couplings with mass dimension to be in the range ±500 GeV, given the cutoff scale M = 1 TeV. The derived range of the g (1) HHH has been presented in Fig. 4. In this region, we believe, our effective theory is well-defined. We will use this parameter region to explore the reach of a 100 TeV collider at 3 ab −1 for probing the five different UV models.
We first investigate the sensitivity of the tri-linear 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 value λ SM at the future collider. Compared with the work in [75], we obtain comparable significance of about 8.25σ for the SM dihiggs production. We exhibit the significance for 3 ab −1 and 0.3 ab −1 with a zero derivative interaction in Fig. 6. The different scales on the left and right vertical axes indicate the statistical significances at the integrated luminosities of 3 ab −1 and 0.3 ab −1 respectively. Systematic uncertainties will lower these values and are not considered as they are not known for the proposed future collider. Based on the plot, with 0.3 ab −1 data, one can exclude values of the tri-linear Higgs coupling λ HHH less than 0.7λ SM . While with 3 ab −1 data, we expect to probe all values of the λ HHH of interest.
Turning on the derivative Higgs coupling g HHH will change the significance of the dihiggs signatures. 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 (2) HHH exist. In Fig. 7 we present the reach of the 100 TeV collider with integrated luminosity of 3 ab −1 in the space of g (1) 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 Fig. 7 are not independent. Their values are connected by Eq. 28, where the contours in the right panel are essentially rotated around the SM values as governed by the Eq. 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 Fig. 8. c tH = 0 corresponds to tan β → ∞, which is outside the experimental bounds on tan β in 2HDMs. Fig. 7 also shows the allowed parameter regions in singlet, triplet and quadruplet models, which overlap within the significance contours. In these models, according to Table III, 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 HHH /v vs g (2) HHH plane, the shaded region is constraint by dimensionless couplings in the Lagrangin within the range ±4π, couplings with mass dimension within the range ±500 GeV and cutoff scale M = 1 TeV. The darker green and light green dotted lines on the left panel correspond to the Wilson coefficient constraints from the Higgs coupling measurements and the T parameter in the real and complex triplet models. The Red line and magenta line corresponds to quadruplet model with Y = 1/2 and 1/3 respectively. The plot on the right hand side are in the v 2 cH vs v 2 (cHD − 4c H ), the shade region are allowed by all the global fit constraints. The SM limit is located at (0, 0) with significance 8.25σ in the right plot. Both plots are with ctH = 0.
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.
For the 2HDM, c tH is non-zero. We demonstrate the significance for 100 TeV collider at 3 ab −1 integrated luminosity in v 2 c tH vs v 2 c H plane, shown in Fig. 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. [86]. This rules out some parameter regions shown in Fig. 8. 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. We also find that, unlike the singlet and triplet, signal significances in the 2HDM are always enhanced compared to 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 Fig. 7 and Fig. 8, if we limit ourselves in these models with all the heavy particles integrated out, the dihiggs process puts additional constraints on the scalar model parameters. Our analysis in Fig. 7 and Fig 8 shows that the real 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 dihiggs process probes the allowed region of c H , and thus the Higgs scalar couplings in the UV models.
Finally, 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 Fig. 7, can be probed with more than 5σ statistical significance. In the complex singlet model, the Wilson coefficients resulting from integrating out the complex singlet can be probed to higher than 7σ statistical significance. In Fig. 9 be explored respectively given 3 ab −1 luminosity data set.
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 dimension-six 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 dihiggs signature starts to provide a strong constraint on c H . Similar to the singlet models, the 2HDM parameters can be probed with a statistical significance of 10σ. In Fig. 10 we show the discovery potential of the model parameter Z 6 vs tan β varies for Type-I model and Type-II model with the 3 ab −1 luminosity. Note that when Z 6 = 0, the SM limit is recovered (see Table. III). We find that with 3 ab −1 luminosity, all the parameter space of interest can be explored for Type-I model and Type-II model. In Type-II model, the top left corner corresponds to the vanishing b quark to SM-like Higgs coupling due to the correction from effective field theory, so there are some regions the significance is smaller than the 8.25 σ limit found for the SM.
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 Fig. 7 (left 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 III 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 III, 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 Fig. 7 (right panel). However, the dimensionless Higgs potential parameters, such as λ HΦ and λ, are still very 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 Fig. 7, the allowed region for two types of quadruplet models are denoted by two lines with different slopes. The 10.

22.
22.  reason can be seen from Table. III, 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 HHH and g (2) HHH can be parameterized by a single parameter |λ Φ3H | 2 . In the right plot in Fig. 7, one can find that the allowed parameter space by the global fit to EWPD for quadruplet models is tightly constrained, and almost becomes a dot near the SM value. In Fig. 9, we show the discovery potential of the model parameter λ Φ3H vs new physics scale M varies with the 3 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 real and complex scalar singlet models. Triplet models 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 Higgs couplings in such models owing to the direct relation between c H and the triple and quadruplet Higgs coupling as shown in Eq. 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 also possible to differentiate the 2HDM and singlet models from 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 compare with the expected SM model value, then it should favor the presence of an extended scalar sector consisting of the 2HDM or singlet under the assumptions of this work. On the other hand, if the future experiment detect a significantly smaller signal rate than the expected SM model value, then it suggests a UV completion of the SM EFT corresponding to the singlet models. We note, however, that both of these statements depend on a lack of new light propagating states which affect the dihiggs signal, an assumption implicit in our use of the EFT framework. 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. Furthermore, depending on the level of deviation one may be able to differentiate between complex singlet model and real singlet model.

V. 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 have investigated the dependence of the coefficient c H on the dihiggs production cross section, and studied simulations of the dihiggs 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 1 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 dihiggs 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 3 ab −1 at the proposed 100 TeV collider, the tri-linear Higgs coupling in the all scalar models with a single heavy scalar integrated out could be fully explored. If a significant deviation in the tri-linear 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 dihiggs 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 dihiggs 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.
Note added: As this work was being prepared a similar publication was submitted to the arXiv [42]. This work contributed to the discussions contained in [42] by expanding on discussions of the 2HDM, including the derivative triple-Higgs coupling corresponding to the coefficient g (2) HHH in Eq. 27 (an approach more suitable for comparison with one-loop analyses and linearly realized gauge symmetry), and a detailed discussion of constraints from EWPD and single Higgs measurements, as well as fully simulated dihiggs production at the proposed future 100 TeV collider.