New Physics in multi-Higgs boson final states

We explore the potential for the discovery of the triple-Higgs signal in the $2b2l^{\pm}4j+E\!\!\!/$ decay channel at a $100$ TeV hadron collider. We consider both the Standard Model and generic new-physics contributions, described by an effective Lagrangian that includes higher-dimensional operators. The selected subset of operators is motivated by composite-Higgs and Higgs-inflation models. In the Standard Model, we perform both a parton-level and a detector-level analysis. Although the parton-level results are encouraging, the detector-level results demonstrate that this mode will be really challenging. However, sizable contributions from new effective operators can largely increase the cross section and/or modify the kinematics of the Higgs bosons in the final state. Taking into account the projected constraints from single and double Higgs-boson production, we propose benchmark points in the new physics models for the measurement of the triple-Higgs boson final state for future collider projects.


Introduction
After the discovery of the Higgs boson h (125 GeV) at the LHC [1,2], measurements of the Higgs self-couplings become crucial for our understanding of fundamental particle physics. In the Standard Model (SM), the Higgs boson has three types of interaction: (1) the interactions with electroweak gauge bosons (W ± and Z); (2) the Yukawa interactions with fermions; (3) the triple and quartic self-interactions. A measurement of the last type of interaction would complete the phenomenological reconstruction of the Higgs potential [3] and thus should lift our knowledge about electroweak symmetry breaking (EWSB) to a new level. Furthermore, Higgs self-interactions could be related to the problems of baryogenesis [4] and vacuum stability [5][6][7].
In the SM, the Higgs potential is written as where H = (G + , 1 √ 2 (v + h + iG 0 )) T is the Higgs doublet, and G ± , G 0 are the unphysical Goldstone bosons associated with spontaneous EWSB in a renormalizable gauge. This potential has a minimum for the Higgs-field vacuum expectation value v = 2|µ|/λ ≈ 246 GeV. After EWSB and switching to unitarity gauge, the Higgs self-interactions take the following form which corresponds to a triple-Higgs self-coupling g hhh = 3 2 λv and a quartic Higgs selfcoupling g hhhh = 3 2 λ, respectively. The parameter λ can be determined by measuring the Higgs mass m h , since λ = 2m 2 h v 2 . In the SM, the Higgs potential is thus completely fixed after the measurement of m h ≈ 125 GeV. However, the story could be different if new physics can contribute to the Higgs self-interactions. Independently measuring the triple and quartic couplings of the Higgs boson via double and triple-Higgs final states is an essential project for future collider experiments.
Deviations from the SM that manifest themselves prominently in double and triple-Higgs final-state processes are expected for various new-physics scenarios. In order to study the Higgs potential in a largely model-independent way, we will parameterize new physics beyond the SM (BSM) in terms of an effective field theory (EFT). This systematic method captures the essence of a wide class of BSM models. It is well suited to collider studies that require exclusive Monte Carlo simulations.
For concreteness, we will describe BSM Higgs physics in terms of the strongly-interacting light Higgs (SILH) version [8] of the EFT approach [9,10]. The operators in this choice of basis are designed to directly correspond to low-energy effects of specific BSM Higgs-sector realizations, including composite Higgs models [8,[11][12][13][14] and the Higgs inflation model [15]. We will consider operators up to dimension 6. Nonvanishing coefficients for some of those, such as ∂ µ (H † H)∂ µ (H † H), can substantially enhance multi-Higgs production rates and/or modify final-state kinematics.
In the SM, the leading order (LO) for the production of one or more Higgs bosons in gluon-gluon fusion involves one-loop diagrams. The calculation of higher order corrections becomes quite a challenge. Most of these calculations [16][17][18][19][20][21][22][23][24][25] are based on effectivetheory methods, working in the limit of infinite top-quark mass. Regarding effects of finite top-quark mass, only NLO QCD corrections to single-Higgs production are known analytically [26,27]. One way to estimate finite top-quark mass effects is series expansion, which can work well for single-Higgs production [23,28,29] but converges poorly for double-Higgs production [30]. Recently, NLO QCD corrections for double-Higgs production with full topquark mass dependence have been calculated numerically [31,32]. The results show large differences in kinematical distributions compared to the prediction of the infinite top-mass limit.
The triple-Higgs self-coupling g hhh can also be measured at a future lepton collider through the double Higgs-strahlung process e + e − → Zhh or the vector-boson fusion process e + e − → ννhh. It has been shown that g hhh can be measured within 27% accuracy at the luminosity-upgraded ILC [58]. At a low-energy machine, such as the 250 GeV CEPC, the triple self-coupling could be determined indirectly via the loop corrections to the ZZh vertex [59,60].
By contrast, a measurement of the quartic self-coupling g hhhh is a real challenge at the LHC, since at √ s = 14 TeV the cross section of gg → hhh is only O(0.01) fb [61,62]. Alternatively, one can consider pp → Zhhh, but that cross section is also tiny [63]. This problem cannot be solved at a lepton collider either, because the cross section for e + e − → Zhhh is only O(0.1) ab at a √ s = 1 TeV machine [64], too small for a measurement.
The proposals for future pp colliders have motivated the study of the process gg → hhh at high energy. The cross section of gg → hhh at a 100 TeV hadron collider can be estimated to be about 3 fb if NLO corrections are accounted for [65], which makes it at least possible to observe the final states of this process. The discovery potential of decay channels hhh → bbbbγγ [66,67] and hhh → bbbbτ τ [68] has been explored. It turns out that the discovery of three-Higgs final state through these channels is challenging, and an extreme high quality detector is needed.
In this paper, we investigate the sensitivity of the decay channel hhh → bbW W * W W * → 2b2l ± 4j + / E, which has not been carefully analyzed in the literature before. We also examine how new physics can contribute to triple-Higgs production. We will consider the effects of a set of dimension-6 effective operators to the cross section and kinematics of Higgs bosons in the final state. Especially, we extend the study of Ref. [15] to the triple-Higgs production case, where the effects of derivative operators on the kinematics of Higgs bosons in double-Higgs production were explored. We also study the projected bounds for all relevant couplings in the EFT at the LHC and at a future 100 TeV pp collider. This paper is organized as follows. In Sec. 2, we briefly introduce the EFT Lagrangian as appropriate for our study and relate our parameterization to particular models that are of interest in the context of new Higgs-sector BSM physics. In Sec. 3, we present a Monte Carlo (MC) analysis of hhh → bbW W * W W * → 2b2l ± 4j + / E in the SM, and investigate the discovery potential and identify challenges of this channel. In Sec. 4, we describe the calculation of triple-Higgs production in the context of the EFT with dimension-six operators in detail and present our numerical results. We conclude this paper with a discussion of our findings in Sec. 5.

Effective Lagrangian up to dimension-6 operators
It has been accepted for a long time that new-physics effects associated with a characteristic scale higher than the energy of the processes under study, can be conveniently expressed in terms of a low-energy EFT. This is a local Lagrangian which includes an infinite series of operators of dimension greater than four, constructed as monomials of fields and organized in terms of the canonical dimension. The operators may incorporate only the unbroken Lorentz, electromagnetic and colour symmetries [69]. However, our knowledge of flavor data, electroweak precision data, and Higgs properties strongly suggests to furthermore implement the power counting of EWSB and thus build operators out of classically gaugeinvariant combinations under the full electroweak symmetry. Up to dimension four, this reproduces the SM. The set of operators up to dimension six was introduced in Ref. [9] and has been reworked to a minimal basis in Ref. [10]. Adopting this as a phenomenological model implies rather generic assumptions on the flavor and gauge structure of the underlying fundamental theory.
In the present context, we are more specifically interested in the possibility that Higgs self-couplings act as primary probes to new-physics effects, while other SM fields are affected only by secondary corrections. This notion is realized by scenarios where the Higgs field acts as the only SM field with sizable couplings to a new sector. Specific models with this property have been proposed, e.g., in Refs. [11,13]. A general discussion can be found in Ref. [8] where the resulting effective low-energy Lagrangian, expanded up to dimension six, has been introduced as the SILH Lagrangian. As expected, and confirmed in Ref. [14], this Lagrangian is equivalent to the basis of Ref. [10], but the assumptions of Ref. [8] on the underlying dynamics suggest a hierarchy between induced tree-level and loop-level coefficients that allows for dropping part of the operator set and thus keeping a more economical number of phenomenological parameters. If we follow this line of reasoning, we can adopt the SILH Lagrangian as the basis of the present phenomenological study. We supply a more detailed discussion below in Sec. 2.1.
For the actual applications in later sections, we can focus on the interactions of the physical Higgs field h, after EWSB and expressed in unitarity gauge. The Lagrangian reduces to Here we confine to the CP conserving operators and omit the CP violating operators. In the SM, we have a 1 = λ 3 = λ 4 = 1 and a 2 = a 3 = κ 5 = κ 6 = c 1 = c 2 = 0. It is understood that the corresponding terms have been removed from L SM , such that they are not double-counted.
Another set of models that couple the Higgs sector to new physics is provided by certain models of inflation. As we show below in Sec. 2.2, this effectively results in the same Higgs Lagrangian, Eq.(2.1). In Sec. 2.3 we briefly review the relation to the EFT version of Refs. [9,10] as it has been applied to the Higgs sector in Ref. [70]. Finally, it can be shown that in a framework that implements a non-linear realization of electroweak symmetry, the result is again equivalent to SILH if equivalent assumptions on coefficient hierarchies are taken [71].
In summary, the phenomenological Lagrangian (2.1) provides a robust parameterization of new physics in the Higgs sector under the condition that no new on-shell states appear in the kinematically accessible range.

The SILH Lagrangian in relation to composite Higgs models
The relevant part of the SILH Lagrangian [8,14], including operators up to dimension six, has the form It includes all the CP-conserving gauge-invariant operators up to dimension six with pure Higgs interactions and Higgs-gauge boson interactions. Some operators such as H † HW µν W µν are not included here since they can be generated by integration by parts from the other operators. There are further operators with fermions coupling to the Higgs, which are omitted here. There is only one dimension-5 operator allowed by the SM gauge symmetry, up to Hermitian conjugation and flavour assignments: (H i ) T C(H j ). It gives rise to the neutrino Majorana mass and violates lepton number, so we do not include it, either.
The SM Higgs may appear as a composite pseudo Nambu-Goldstone (NG) boson associated with some enlarged symmetry beyond the SM. The Lagrangian L SILH then emerges at low energy via spontaneous breaking of that symmetry. Since any terms in the Higgs potential will violate the shift symmetry of this NG-boson Higgs, the coefficients above are all suppressed by the small breaking in relation to the compositeness scale f , i.e., carrying a ξ = v 2 f 2 factor. m ρ , g ρ stand for the characteristic mass and coupling of a strongly coupled sector, respectively, and c i ∼ 1.
We focus on the first five operators in Eq. (2.5), since they are the relevant operators for the hadron-collider processes that we want to study. The first three terms in L SILH contribute to the Higgs potential. They contain only two independent terms, as can be  [13,80]. Some notation is from Ref. [8,81]. Note that c 1 and c 2 are sensitive to the detailed construction of the models. We consider c g as roughly of order 1. For the relation between our conventions and other conventions used in the literature, cf. Appendix A.
verified by applying the equations of motion. After EWSB, the SILH potential reduces to the effective potential of Eq. (2.1). We list the relations between Eq. (2.1) and Eq. (2.5) in Table 1. Note that we have the relation κ 5 = κ 6 , since the associated terms come from the same operator c H 2f 2 ∂ µ H † H ∂ µ H † H . The rest of the operator coefficients can be measured at future electron-positron colliders, via W -pair production, Z-pair production, and Z-Higgs production [72,73].
Regarding hadron-collider measurements, the coefficient c g is accessible via the pp → h process at the LHC. Run-1 data have constrained c g /m 2 ρ ∼ 10 −6 [74]. Bounds for the coefficients c H , c y and c 6 are currently much weaker [75]. It is expected that the highluminosity LHC will yield bounds c H ξ ∈ [−0.044, 0.035] and c y ξ ∈ [−0.020, 0.008] for the top quark [76]. The coefficient c H ξ can be further constrained to O(10 −3 ) at a future e + e − collider, and the tests for c y can be extended to b, c quarks, and leptons [77]. The coefficient c 6 contributes to the triple and quartic Higgs self-couplings only, so the bounds on c 6 will stay relatively weak for both LHC and a future lepton collider.
We may also consider two more specific composite Higgs models [8,13], dubbed as MCHM4 and MCHM5, respectively. Both models result in the SILH Lagrangian as their low-energy EFT. They contain extra fermions, which are in representations 4 and 5 of an assumed global SO(5) symmetry, respectively. We adopt the notation from Ref. [78]. The SILH coefficient values are The current LHC constraints and electroweak precision data imply f ≤ 550 GeV and v 2 /f 2 ≤ 0.2 [79]. Later, we will study the projected constraints from the LHC and a 100 TeV collider.

Operators from Higgs inflation
In this section, we demonstrate how an equivalent set of dimension-six operators arises from the standard Higgs inflation paradigm [82][83][84][85][86]. We incorporate a non-minimal coupling of the Higgs field to gravity and work in unitarity gauge where H = (0, h/ √ 2). The gauge interactions are more complicated in this scenario; we ignore them for now and just focus on the Higgs potential. In the Jordan frame, the Lagrangian has the form We consider ξ in the range 1 √ ξ ≪ 10 17 , in which M M PLanck . We perform a conformal transformation from the Jordan frame to the Einstein frame, This transformation will give rise to derivative terms in Higgs potentials. We furthermore redefine Then the action in the Einstein frame is given by where the potential becomes In the standard Higgs inflation paradigm, h takes large values h M Planck / √ ξ (or χ √ 6M Planck ) during inflation and plays the role of the inflaton. We have the expressions This allows the potential to be exponentially flat at large h to produce a viable inflaton potential.
When the value of h is near the origin as today, we can approximate h χ and Ω 2 1, so the potential for the field χ generates a potential for the SM model Higgs field plus corrections at O(ξ 2 /M 2 Planck ). For the purpose of this collider study, we thus replace χ by h. Plugging Eq. (2.10) into Eq. (2.11) and omitting higher order terms, we arrive at (2.14) Note that after EWSB, replacing h → h + v yields similar extra terms as in Eq. (2.1).

Alternative Parameterization of the Higgs boson self-interaction operators
Another representation of the set of gauge-invariant dimension-6 operators which can modify the Higgs self-interactions, has been studied in Ref. [70] O The operator O 1 was considered in Ref. [75] and can safely be neglected. In the subset Therefore the kinetic term of the Higgs field is modified to This means that the Higgs field should be rescaled by h → ζh, where ζ = (1+f 2 v 2 /Λ 2 ) −1/2 . After EWSB and choosing unitarity gauge, the Lagrangian reduces to (2.1) as before, where the coefficients (a 1 , λ 3 , λ 4 , κ 5 , κ 6 ) of Eq. (2.1) can be expressed in terms of just two independent parameters:x = x 2 ζ 2 , (2.21) 3). With this definition, the rescaling factor ζ can be rewritten as ζ = (1 −x) 1/2 . The relations between our parameters and those in Ref. [15] are listed in Table 2.
At the time when the measurements that we discuss in the present work can be carried out, we should expect that data exist that set significant bounds on the parametersx andr.
1. a 1 (x) is related to the direct measurement of the top Yukawa coupling, and its value is expected to become determined within 5% precision at the high-luminosity LHC [87], via measuring the tth production rate. At a 100 TeV collider, the Yukawa coupling can be pinpointed down to a precision 1% [88] by measuring the ratio between the tth and ttZ production rates.
Our operators Operators in Ref. [15] Relations Table 2. Parameter relationship between our convention and that in Ref. [15].
2. Another bound onx is obtained from the measurement of Higgs-gauge couplings [75], since they become universally rescaled by ζ. A future e + e − Higgs factory can constrain |x| at the 1% level [89]. Since there are many other dimension-6 operators which can contribute to the gauge-boson kinetic terms, we nevertheless takex as a free parameter in our later analysis.
3. The parameterr can only be constrained by double-Higgs or triple-Higgs production.
Concerning double-Higgs boson production, the bound will be around 40 ∼ 100% at the HL-LHC at most. At a 100 TeV hadron collider, (x,r) will become more strongly constrained by double-Higgs production. As shown in Ref. [15], the bounds onx and r will be of the order 2 ∼ 5% and 4 ∼ 13%, respectively.
3 Detailed analysis of the 2b2l ± 4j + / E channel in the SM We study triple-Higgs production in high-energy proton-proton collisions, pp → hhh, where one Higgs boson decays into a bb pair while the two other Higgses decay into W W * . The semi-virtual W pairs can subsequently decay semileptonically, h → W W * → νjj.
The dominant partonic contribution to the pp → hhh signal is gluon-gluon fusion, gg → hhh. This process involves one-loop diagrams. As we did for our previous work [67], we compute the production matrix element at LO with MadLoop/aMC@NLO [90]. We take the parton distribution functions from CTEQ6l1 [91]. For phase-space evaluation and exclusive event generation, we interface the production process with VBFNLO [92][93][94].
Background event samples are generated by MadGraph 5 [95,96]. Since we require a bb pair, the dominant background is caused by top-quark pairs in association with electroweak bosons, namely pp → h(W W * )tt and pp → ttW − W + . Both classes of processes can lead to the same final state as the signal. To veto further background from Z bosons, we restrict the analysis to same-sign leptons in the final state, l + l + or l − l − .
We list the calculated cross sections of signal and backgrounds at 100 TeV in Table 3. In the absence of a complete NLO calculation for the signal, we adopt the K-factor of 2.0 that was obtained in Ref. [66] for Higgs pair production. For the H(W W * )tt background, we use K = 1.2 [97]. The K-factor for ttW − W + at 100 TeV is taken 1.3 from Ref. [56]. In the Ref. [98], a K-factor around 1.2 was obtained while the total cross section σ NLO is given as 1.3pb, which is around 1.4 times larger than our LO cross section. The derivation is mainly attributed to the choice of renormalisation and factorisation scales, i.e. our choice of K-factor equal to 1.3 is consistent with the results given in Ref. [98] after taking these uncertainties into account.
We ignore all background from h+jets, hh+jets and W ± W ± +jets, since the cross sections of those processes are negligible compared to the h(W W * )tt background. Furthermore, we observe that the total cross section of the background bbW − W + W − W + is essentially exhausted by the resonant contribution ttW − W + . Therefore, we approximate the former process by the latter with subsequent top-quark decay, which considerably simplifies the calculation.
We have three comments on the background processes hhjj in the SM and new physics models.
• In SM, the hhjj final state receives contribution for heavy-quark loop and vector boson fusion, while the former is dominant. Currently, the cross section of loop-induced processes with 2 jets can be calculated by interfacing GoSam [99] or OpenLoops [100] to Madgraph5 [101] or Herwig7 [103]. We use Madgraph5 to compute the cross section of top quark loop induced pp → hhjj at a 100 TeV collider. After imposing the MLM matching [102] and using cuts P t (j) > 20 GeV and η(j) < 5, we obtain an inclusive cross section 620 fb, which is around 128 times larger than the cross section σ(hhh) of the signal processes gg → hhh. Meanwhile, by using Madgraph5 [96], we find that the cross section of VBF with √ s = 100 TeV is 34 fb.
• It is known that when the b tagging efficiency is taken as 0.7, the rejection rate of light jets can reach 0.1% or so. Since we required one(two) tagged b jets in our preselection cuts, therefore the background gg → hh + 2jets is suppressed by a factor 10 −3 ( 10 −6 ) or so. After imposing b taggings and the decay branching fraction of h → bb, we find that the signal cross sectionbbhh is around 0.52(0.29) σ(hhh), while the cross section of background hh + 2jets is 0.13(0.13 × 10 −3 ) × σ(hhh) or so. Obviously, when n b ≥ 2 is imposed, it is safe to neglect this type of background in the SM.
• In the new physics models we will consider below, the background process of hhjj can have extra contributions from higher dimensional operators. When the cross section is 2 ∼ 5 magnitude orders smaller than the signal process, we can neglect it safely. In the cases when such a background is greatly enhanced or in the cases the signal process gg → hhh is greatly suppressed by the higher dimensional operators to such a degree that the cross sections of them are comparable, the background of hhjj should be included in the analysis. Table 3 shows a yield of 642 signal events in this final-state channel for 30 ab −1 integrated luminosity. However, without further selection there are ∼ 10 7 background events. Clearly, it is a challenge to observe triple-Higgs production through this channel. In the following subsection we discuss observables and selection methods for suppressing background and raising the signal/background ratio to an acceptable level.  Table 3. Cross sections of signal and background for the 2b2l ± 4j + / E final state in the SM. The expected number of events corresponds to 30 ab −1 integrated luminosity.

Parton-level analysis
We simulate the Higgs boson decays that lead to the final state 2b2l ± 4j + / E by using the DECAY package provided by MadGraph 5. Here we do not consider any parton shower effects, which will be discussed in section 3.2. The transverse momentum (P t ) distributions of the visible particles and missing transverse energy (MET) are shown in Fig. 1. In this figure, the objects are sorted by P t . On the one hand, one can expect that the b quarks are harder than the light quarks, since they originate from a Higgs boson decay directly. On the other hand, the decay chain h → W W * → jjlν leads to soft leptons and light jets, especially when they are coming from the off-shell W bosons.
In Fig. 1(b) and Fig. 1(c), we observe that the P t distributions of the softest leptons and jets peak around 10 GeV, which might make it a challenge to successfully reconstruct these objects with the currently planned detectors. Since the signal contains only two neutrinos, MET should not be too large. As illustrated by Fig. 1(d), MET peaks around 50 GeV, somewhat below half the Higgs boson mass.
Because there are two unobserved neutrinos in the final state, their mothers being either on-shell or off-shell W bosons, it is not convenient to fully reconstruct the Higgs bosons. A partial reconstruction should nevertheless be possible. In order to extract this information, it is crucial to correctly associate the mother Higgs bosons with their decay products. Here we encounter a problem of combinatorics, which leads to a 12-fold ambiguity. To simplify the problem, we assume that both b quarks can be tagged correctly, so only the light quarks can be reassigned and the ambiguity reduces to 6-fold.
To find the correct combination of the visible particles from Higgs boson decays, we examine the following four alternative reconstruction methods at parton level: 1. The decay chain h → W W * → jjlν suggests that the lepton and the hadronically decayed W boson should have a small angular separation ∆R(l, W jj ). Since there are two Higgs bosons with this decay chain, the sum of ∆R 1 (l, W jj ) + ∆R 2 (l, W jj ) should be minimal. We choose a combination with minimal value of this observable.
2. The semileptonic Higgs invariant masses can be computed from the visible particles; we denote them as m vis h1 (l, jj) and m vis h2 (l, jj). We choose a combination which minimizes their sum.

Methods
The percentage of correctness Table 4. Methods for determining the correct combinations of (l, j, j) and their percentages of correctness.
set an upper bound on the Higgs mass, so we choose a combination which minimizes mT 2.
4. Since mT 2 should have a value close to the Higgs mass m h = 126 GeV, we choose a combination which minimizes |mT 2 − m h |.
These methods and their associated percentages of correct assignment in a simulated event sample are listed in Table 4. The effect of realistic b-tagging efficiency will be discussed in the next subsection. We observe that in a parton-level analysis, the method that relies on the quantity |mT 2 − m h | has the best performance, approaching 100% probability for correct particle assignment in the reconstruction.

Detector-level analysis
To obtain a hadronic event sample, we use the parton-shower and hadronization modules of Pythia 6.4 [109]. For jet clustering, we use the package FASTJET [110] with the anti-k t algorithm [111] and cone parameter R = 0.5. To veto the large number of soft jets from initial-state radiation, only jets with P t > 20 GeV are accepted.
The multiplicity distribution of jets is plotted in Fig. 2(a). Both signal and background in the MC sample provide six jets at parton level, which explains the peak of n j around 6 in Fig. 2(a). In Fig. 2(b), we show the P t distributions of the six leading jets in the signal event sample. The 1st to 4th jet exhibit similar distributions as at parton level, but the 5th and 6th jet P t distributions have different shapes with respect to their parton-level counterparts.
There are two simple reasons for this result: (1) the softest quark in Fig. 1(c) typically has P t only around 10 GeV while most of the low-P t jets are vetoed by our P t > 20 GeV cut; (2) jets from initial-state radiation can easily be as hard as 20 GeV at a 100 TeV collider. So the 5th and 6th jet are more likely produced by initial-state radiation than by Higgs boson decays. Fig. 2 illustrates the challenge of reconstructing the soft jets generated by the multi-Higgs signal. Another important problem is the reconstruction of leptons. We assume that the future detector can reach a better efficiency in reconstructing leptons than possible today (95% for P t > 5 GeV), so it becomes feasible to find the soft lepton as shown in Fig. 1(b). But in order to reject huge QCD background, we need isolated leptons. To find a suitable isolation condition, we investigate the angular separations between two leptons (∆R(l, l)) and between leptons and jets (∆R(l, j)), respectively. The minimum-value distributions of these two observables at hadron level are displayed in Fig. 3. On the one hand, min ∆R(l, l) tends to have a large value, and only 10% of the events have min ∆R(l, l) < 0.5. On the other hand, almost 50% of the events have min ∆R(l, j) < 0.2. This makes it difficult to isolate the leptons from the jets.
To study the detector effects, we use DELPHES [112,113] to perform a detector simulation for the generated event samples. The setup of DELPHES is similar as in Ref. [67], with the following modifications: 1. The b-tagging efficiency is assumed to be a constant b = 0.7, and mistagging rates are 0.1 and 0.001 for charm and light jets, respectively. The pseudorapidity for b (c, jet) is required to be η < 5.0, respectively.
2. As described above, the jets are clustered by FASTJET with a cut P t (j) > 20 GeV.

Isolated leptons are defined by Ref. [113]
where l is a lepton. The sum in the numerator runs over particles with transverse momenta above P min t = 0.1 GeV within a cone with radius R = 0.5, except for l. A lepton is classified as isolated if I(l) < 0.1. Fig. 4 shows the number of b jets and isolated leptons after detector simulation. Since both signal and backgrounds include two b jets, it is easy to understand the similarity of the shapes in Fig. 4(a). However, only 10% of the signal events are found to include two leptons ( Fig. 4(b)), which makes it difficult to separate signal from background. The small value of min ∆R(l, j) in typical events explains this result.
To further enhance the signal over background ratio, we apply three preselection cuts: 1. The number of b jets is required to be n b ≥ 1. One might worry about the background hh+2 jets. It is found that it can only contribute around 5 events, which can further be reduced to 3 by the cut |m bb − m h | < 58 GeV and is much smaller than the other two types of background given in Table (5), it is safe to omit it here.
2. To veto background from a Z boson, we require two same-sign leptons, as discussed above. Note that this also removes triple-Higgs signal events which decay to oppositesign leptons.
3. The number of light jets is required to be n j ≥ 4.
We are interested in three observables: (1)  Nevertheless, we can try to suppress background by applying cuts on the above observables. The efficiencies of each cut are listed in Table 5. The significance of the signal in the cut-based method finally amounts to just 0.02, which is clearly much worse than could be expected from the parton-level calculation. We conclude that in the SM, a discovery of triple-Higgs production through this channel will be extremely challenging.   Table 5. Efficiencies of cuts as described in the text, for a total integrated luminosity of 30 ab −1 ..

Triple-Higgs production with dimension-6 operators
Given the dim prospects for observing triple Higgs production in the pure SM, we may ask the question about SM extensions that enhance the production rate such that the process becomes observable at a 100 TeV collider. In that case, such an observation would not just indicate a significant deviation from the SM, but at the same time provide a measurement of new BSM parameters. We work in the context of the genuine Higgs-sector BSM models that we have introduced above, conveniently parameterized by the SILH Lagrangian with dimension-six operators, or, alternatively, by the effective Higgs Lagrangian in unitarity gauge, Eq. (2.1). In the SM, the production process gg → hhh proceeds via top-quark loop diagrams coupled to Higgs bosons and involves triple and quartic Higgs couplings, where we are obviously most interested in the quartic-coupling contribution. The various anomalous couplings generated by Eq. (2.1) modify all contributing loop Feynman graphs, and furthermore direct Higgs-gluon couplings can appear which are induced either from the underlying theory directly or emerge from loop-diagram renormalization via operator mixing. Therefore, we redo the calculation of the production process at one-loop order with all new parameters included.
While the observation of the triple-Higgs process ultimately would determine a particular combination of the new parameters, sizable values for those will definitely also affect other, more easily accessible processes such as Higgs production in association with a top quark and double-Higgs production. Any actual measurement of the EFT parameters will involve a fitting procedure that takes all available information into account. Nevertheless, the fact that the quartic Higgs coupling appears only in triple-Higgs production indicates that the current process will contribute independent, and potentially essential information.
If the triple-Higgs final state is to become observable, we have to allow for EFT parameter values that distort the amplitudes rather drastically, at least in the high-energy or high-P t regions of phase space. We should worry about unitarity of the amplitudes and consistency of the EFT. Physically, we expect a dampening effect from strong rescattering of intermediate real top quarks into multiple Higgs bosons. While the effects of strong rescattering have extensively been studied in the linear EFT context for vectorboson scattering [114,115], no results are available for processes involving top quarks and Higgs bosons. Power counting suggests that the dimension-six operators that we consider in this work are affected to a lesser extent than the dimension-eight operators considered in Ref. [114]. We also note that in the SM as a weakly interacting theory, the Higgs mechanism tends to suppress electroweak production cross sections by orders of magnitude in relation to the bounds enforced by unitarity, so there is a significant margin for enhancing event yields in BSM models. For the current study, we take the EFT unmodified over the complete parameter space and defer a study of constraints and relations imposed by unitarity to future work.

Calculation
In this section, we describe our calculation of the one-loop induced production amplitude in the presence of the new parameters of the unitarity-gauge effective Lagrangian, Eq. (2.1). We use the package Madgraph5/aMC@NLO [95,96] for calculating the loop diagrams and, subsequently, evaluating phase space and generating event samples. To this end, we have implemented the model described by Eq. (2.1) as a UFO model file. Note that our choice of unitarity gauge for the electroweak symmetry does not affect the QCD loop calculation, since the color symmetry is implemented in a renormalizable gauge as usual. The program reduces the one-loop Feynman integrals to scalar integrals in four dimensions, employing techniques such as the OPP method [116]. The difference between the D-dimensional and 4dimensional expressions that arises in the calculation yields additional rational terms [117]. These are identified as R1 terms associated with D-dimensional denominators, and R2 terms associated with D-dimensional numerators. All R1 terms are automatically generated as a byproduct of the reduction method, while the R2 terms must be calculated manually [118]. We have performed this calculation and supplied the results as effective tree-level vertices in the UFO model file [119].
In particular, we obtain the R2 terms that amount to contact interactions of a pair of gluons with one to three Higgs bosons: The coefficients depend on the EFT parameters a 1 , a 2 , and a 3 . Since these terms are required to restore the exact QCD symmetries in the calculated amplitude, by themselves they manifestly violate gauge invariance. We have verified that the complete renormalized one-loop result does respect gauge invariance, a convenient cross-check of the calculation.
Besides these loop-induced contributions, gluon fusion into Higgs bosons also receives contribution from contact interactions between gluons and Higgs bosons that do not exist in the SM. As mentioned above, the inclusion of such contact interactions is required by loop-induced operator mixing in the EFT, but could also originate from independent BSM contributions. Technically, we implement them as independent R2 terms, so that Madgraph5/aMC@NLO will sum them together with loop-induced contribution.

Cross sections of gg → hhh and Kinematics
The amplitudes of the process pp → hhh are constructed from the Feynman diagrams in Fig. 4.4. For illustrating the method, we take the terms that depend on a 1 , κ 5 , κ 6 , λ 3 and λ 4 . Complete results are given in the Appendix B.

(4.6)
To determine the integrated form factors t 1 . . . t 30 , we calculate the total cross section at 480 selected points in the space of parameters (a 1 , λ 3 , λ 4 , κ 5 , κ 6 ), then obtain the numerical values of these coefficients t 1 . . . t 30 via linear regression. The resulting values are shown in Table 6. The complete set of results that accounts for all effective operators is provided in the Appendix B.  Table 7.

Numerical values of integrated form factors
To study the parameter dependence of the cross section of gg → hhh in one of the more specific models introduced in Section 2.3, we simply replace (a 1 , λ 3 , λ 4 , κ 5 , κ 6 ) by (r,x) according to Table 2, so we obtain Eq. (4.7) below. The numerical results fort 1 . . .t 14 are listed in Table 7, and the total cross section has the value σ hhh SM = 5.84 fb. (This includes a K factor of 2.0, following Ref. [65]). 3 (1 +t 1x +t 2r +t 3x 2 +t 4xr +t 5r 2 +t 6x 3 +t 7x 2r +t 8xr 2 +t 9r 3 +t 10x 4 +t 11x 3r +t 12x 2r2 +t 13xr 3 +t 14r 4 ) A visual representation of the cross-section dependence on the parameters (x,r) is shown in Fig. 6. We read off that the cross section can exceed the SM value by two orders of magnitude for reasonable variations of (x,r). In particular, ifr is fixed ( Fig. 6(a)), the cross section increases in thex < 0 region. In this region, all of the dependent parameters λ 3 , λ 4 , κ 5 , and κ 6 have the same sign, and the derivative couplings can greatly enhance the cross section. By contrast, in thex > 0 region, the contributions of λ 3 and λ 4 cancel against the terms with κ 5 and κ 6 .
The complementary plot Fig. 6(b) shows the dependence onr for fixedx. The cross section changes only mildly withr as long asx is small or positive, and forr > 0 it actually undershoots the SM value. We recall that the dominant contribution to triple-Higgs production originates from the diagram with a pentagon top-quark loop [67]. This part does not depend on the Higgs self-couplings which enter the parameterr. Only if the latter contributions become sizable and the interference is constructive, we expect a large enhancement of the cross section.  Beyond the effect on the total cross section, we may look at distortions of kinematical distributions. The P t distributions of the three Higgs bosons are shown in Fig. 7(a), Fig. 7(b), and Fig. 7(c), for three different values ofx: −0.5, 0 and +0.5, respectively. We observe that the distributions change significantly with respect to the SM reference value ifx = +0.5, especially in the large P t region. The distortion happens in the parameter region where the total cross section is not enhanced by a large factor, and it is helpful for the 2b2l ± 4j + / E channel since it should improve the reconstruction of the softest jet. For x = −0.5 the distributions do not change that much, but the analysis would benefit from the remarkable cross-section enhancement in that region. We also show the invariant mass distribution of the three Higgs bosons (Fig. 7(d)); this is also modified by the derivative operator.
In Fig. 8, we show the same observables as in Fig. 7; this timer is varied andx is fixed to zero. The distributions do not actually depend onr, since the parameter affects only λ 3 and λ 4 which are not associated with derivative couplings.

Correlations between gg → hhh and single and double-Higgs production
A measurement of triple-Higgs production would not occur in a vacuum. Apart from the genuine quartic Higgs couplings, all parameters of the EFT also enter other Higgs processes, and the discussion in Sec. 2 suggests that in typical strongly-interacting models, all parameters would receive BSM contributions. We can expect that a measurement or exclusion limit on the triple-Higgs process would add information on top of the amount of Higgs-physics data gained up to that point, and all results should be combined in the interpretation towards BSM physics. Therefore, in this section we study correlations between gg → hhh and gg → h, gg → hh in particular. We phrase the problem in terms of the following questions: • To what extent can a 1 and c 1 be determined from gg → h at the LHC (14 TeV) and at a 100 TeV collider?
• To what extent can a 2 , c 2 , λ 3 , κ 5 be determined from gg → hh at the LHC and at a   Table 9. Parameters that contribute to the particular Higgs-production processes.

TeV Collider?
• To what extent can a 3 , λ 4 , κ 6 be determined from gg → hhh at a 100 TeV collider, including channels not considered in this paper?
In Table 8 we quote estimates for the theoretical and projected experimental uncertainties for the processes gg → h, gg → hh and gg → hhh. For convenience, we list the parameters that enter these processes in Table 9. The theoretical uncertainties are obtained by summing the squared uncertainties in PDF, renormalisation scales, and α s , based on current knowledge. For the process gg → h, the experimental uncertainties are mainly statistical ones which pertain to the Higgs decay h → γγ. The projected experimental bound for gg → hh at the LHC is taken from the studies of the bbγγ final state [15] and 3 2j+MET [55]. We also quote the expected experimental bound for gg → hhh which is expected from the analysis of 4b2γ final states [67] at 100 TeV. We emphasize that these estimates are derived from phenomenological studies; full simulation and experience gained in the analysis of actual data may change the conclusions significantly, such as in the expectations for the observation of Higgs-pair production at the LHC [74].
We first consider Higgs-pair production in gluon fusion, the dominant contribution to the process pp → hh. The dependence of the total cross section on the EFT parameters can be written as in Eq. (4.8). Explicitly, we have The numerical values for the coefficients f 1 − f 10 at a 100 TeV hadron collider are provided in Table 10.   Table 11. Numerical results forf 1 −f 5 at a 100 TeV hadron collider.
Turning to single-Higgs production, in Fig. 9 we show projections on the bounds of a 1 and c 1 from the process gg → h at the LHC (14 TeV) and at a 100 TeV collider, respectively. Assuming a measurement result equal to the SM prediction, the allowed ranges for a 1 and c 1 are highly correlated and are confined to be two narrow bands, one of which containing the SM reference point. The other band is centered on a mirror solution a 1 = −1, c 1 = 0.
At the 14 TeV LHC, both theoretical and experimental (statistical) uncertainties are relevant and have to be taken into account. By contrast, at a 100 TeV collider the main uncertainties will come from theory, while statistical uncertainties will become less than 0.1%. Compared with the ultimate LHC bounds, the projected accuracy on the determination of a 1 and c 1 from a 100 TeV collider will improve by a factor 2. If theoretical and parametric uncertainties can also be improved in the future, we can expect the bounds on a 1 and c 1 to tighten even further.
The correlation of the parameters in Fig. 9a and Fig. 9b indicates an obvious shortcoming of a simple cross-section analysis. The degeneracy of solutions can be lifted by examining the kinematics of the Higgs boson in the final state, and by adding the measurement of gg → hh, as we will explain later.
A future e + e − collider can improve the measurement of a 1 and c 1 beyond the reach of the LHC. This is mainly due to QCD entering the processes only beyond leading order, so statistical uncertainties will likely dominate. For instance, at the CEPC a combination of the parameters c 1 and a 1 can be determined within 1% for a collision energy √ s = 240−250 GeV and an integrated luminosity 5 ab −1 [77]. At the ILC, the parameter a 1 can be determined within 10% for √ s = 500 GeV and an integrated luminosity of 1 ab −1 , via the process e + e − → tth [121]. We add the result that at a 100 TeV hadron collider with integrated luminosity 30 ab −1 , a 1 can be further constrained by the measurement of tth to within 1% [88].
In Fig. 10, we display the expected bounds in the a 1 -c 1 plane for both the LHC 14 TeV and a 100 TeV collider. Adding in gg → hh, the linear degeneracy in the a 1 -c 1 plane that (a)   In this situation, we may include the cross section of the process gg → hhh and distinguish the second point from the rest. The second point has a production rate for the process gg → hhh that is large enough to be observed. (Here we would like to remind the reader that when we examine the correlations of a 1 and c 1 , we set the other free parameters to their SM values.) We now extend the study to the other parameters that enter gg → hh. In Fig. 11 we display the expected LHC bounds for those in four planes, namely a 2 -λ 3 , c 1 -λ 3 , c 2 -λ 3 , and κ 5 -λ 3 . All bounds are derived by requiring the cross section of gg → hh to be smaller than 140 fb, so the points inside the exclusion bounds will remain allowed by the LHC 14 TeV data, if no deviation from the SM is detected. The SM reference points are also shown. When we examine the correlations between two parameters, here and in all later discussions we set the remaining parameters equal to their SM values.  Figure 11. Projected exclusion bounds in two-parameter planes between (a 2 , c 1 , c 2 , κ 5 ) and λ 3 , extracted from the process gg → hh at the LHC 14 TeV. If the coefficient values are equal to the SM prediction, parameter values inside the contours are still allowed by the measurement. The exclusion bounds correpond to a limit of 140 fb for the cross section.
In Fig. 12, we draw the analogous bounds that we can expect from analyzing gg → hh in 100 TeV pp data. We assume that the cross section of gg → hh can be measured to a precision of 8%. This 8% uncertainty results from combining theoretical and experimental uncertainties. The allowed parameter regions shrink considerably and become pinched between two contours, in each plot. Comparing Fig. 11 and Fig. 12, we conclude that a 100 TeV collider can significantly improve the precision on a 2 , c 2 , κ 5 , and λ 3 .
We can also individually project out single-parameter bounds for each of these four parameters, shown in Fig. 13. In this approach, the parameters a 2 , κ 5 , and λ 3 can be determined with a precision close to 10%. The two-fold ambiguities in the solutions could possibly be removed by using the kinematics of final states, as demonstrated in Ref. [55]. The parameter c 2 can be determined within the range [−0.1, 0.4].
Finally, we add the exclusion bounds that we can expect from the process gg → hhh. There are three parameters, namely a 3 , κ 6 and λ 4 , which only contribute to the process These results have to be combined with the parameter exclusion regions derived from gg → hh. In Fig. 15, we show the correlations between a 3 and a 2 , c 2 , κ 5 , and λ 3 . In these plots, we overlay the limits that follow from Higgs-pair production, presented in Fig. 13, to the exclusion contours that follow from triple-Higgs production. The SM prediction is displayed for reference. Clearly, the triple-Higgs results yield weaker constraints, but they are nevertheless sensitive to a different combination of parameters and thus cut off part of the two-parameter exclusion regions. As an example, we note that adding in gg → hhh can help in resolving a two-fold ambiguity in the κ 5 -a 3 plane, cf. Fig. 15(c).

Analysis for models
Finally, we can adapt the above results to more specific scenarios, as we have introduced above in Sec. 2. In any given model, the parameters of the unitarity-gauge Lagrangian (2.1) can be related to the original model parameters. In particular, for a model with a small set of independent parameters, we can recast the analysis to a concrete prediction for the expected sensitivity to this model, in a straightforward way.  Figure 14. Projected two-parameter (left) and one-parameter (right) exclusion bounds, extracted from the process gg → hhh. The left column shows the two-parameter planes a 3 -λ 4 , κ 6 -λ 4 , and a 3 -κ 6 , while the right column displays a 3 , κ 6 , and λ 4 .

Strongly-interacting Higgs models
In the generic dimension-six SILH Lagrangian, there are four free parameters relevant to this stody, denoted by C y = c y ξ, C H = c H ξ, C 6 = c 6 ξ, and c 1 . For both MCHM4 and MCHM5 as models which reduce to SILH at low energy, there are only two independent parameters, ξ and c 1 . We apply the above analysis to those and conclude that at a 100 TeV collider, data from gg → h and gg → hh significantly constrain the allowed parameter space. This is demonstrated by Figs. 18(a)-18(b). Data from gg → h will result in an exclusion region bounded by two lines in the ξ − c 1 plane, while data from gg → hh will further reduce the allowed parameter space to two small spots in the plane.
To illustrate the added value from triple-Higgs production, we consider two benchmark points for MCHM4 and MCHM5 in Table 13. For both benchmark points we obtain a large cross section for the gg → hhh process, actually 80 and 55 times larger than that of the SM, respectively. These are examples of benchmark models that can not just be detected, but also be distinguished from each other at a 100 TeV collider, given the result   Table 13. Three representative points for the models MCHM4, MCHM5, and GHM, respectively, and the corresponding cross sections for Higgs production. that the threshold cross section value of gg → hhh for being sensitive to new physics is 30 fb or smaller [67]. (We add the caveat that the benchmark point of MCHM4 could be independently excluded by incorporating electroweak precision data due to a large ξ value.)

The Gravity-Higgs Model
The Gravity-Higgs model has only two free parameters,x andr. The analysis is straightforward, since in this model, single-Higgs production gg → h depends only on a single BSM parameterx, while double-Higgs production constrains the second parameterr. The cross section of gg → hhh is completely determined oncex andr are constrained.
The expected LHC exclusion contours in thex −r plane are depicted in Fig. 19(a). From the process gg → h we obtain a narrow band which is cut off by adding in the result from measuring gg → hh. The latter constrains the parameterr down to the range from -1.8 to 5.0 if we assume that parameter space with a cross section of σ(gg → hh) larger than 120 fb can be excluded. Fig. 19(b) shows the analogous results for a 100 TeV collider. The measurements of gg → h and gg → hh let the allowed parameter space shrink substantially compared to the LHC results. From these measurements alone, the value ofr is extracted with a two-fold ambiguity. However, for the solution with the larger value ofr (r ≈ 3) the cross section of gg → hhh is 5 times larger than nearr = 0 due to the λ 4 3 dependence in the cross section, cf. the benchmark points in Table 13. A measurement of triple-Higgs production could therefore eliminate one of the solutions.

Conclusion and discussion
In this paper, we have explored the discovery potential for triple-Higgs production via the 2b2l ± 4j + / E decay channel at a 100 TeV hadron collider. Despite the extremely small cross section of the signal process in the SM, the parton-level study demonstrates that the signal can be detected in principle. We find that the mT 2 variable is useful to find the correct combinations of the visible objects that originate from Higgs boson decay, and to suppress background efficiently. However, once hadronic events and detector effects are properly accounted for, extraction of the SM signal becomes a real challenge. The two main problems are that (1) the transverse momentum of the softest jet from Higgs boson decay assumes values of about 10 GeV, which makes it difficult to reconstruct; (2) since there are six jets in the final state, the lepton isolation criteria reject most of the signal events.
New-physics effects may enhance the cross section significantly and distort kinematical distributions, so if the SM is not the true theory, observation of the triple-Higgs production process becomes more likely. A measurement would then amount to a determination of BSM parameters. To discuss this possibility in a suitably generic framework, we have employed an EFT approach and added a set of dimension-6 operators that is taylored to represent new physics in the Higgs sector. In particular, we find that a sizable coefficient for a derivative operator can modify the kinematical distributions of the visible objects such that a reconstruction of the triple-Higgs signal becomes feasible. Using this information, we have investigated the potential of such a measurement to improve on knowledge which can already be gathered from single and double-Higgs production data. It turns out that while those processes are generally more powerful in constraining BSM parameters, the triple-Higgs signal nevertheless reduces the allowed parameter space and in some cases can eliminate ambiguities in the parameter determination.
We would like to point out that our EFT approach does not incorporate any new BSM particles which may be discovered in the future. Our study also treats the SM particles, especially the Higgs, as elementary degrees of freedom at the energy scale relevant for 100 TeV collider. We thus assume that, if the SM Higgs is a composite particle, the compositeness scale is higher than a few TeV, at least. If this assumption turns out to be invalid, i.e., qualitatively new phenomena become observable at lower energy, our conservative approach would have to be revised to include explicit model-dependent BSM effects in the calculation.

A Derivation of parameters relations
Expanding the SILH Lagragian in unitarity gauge and introducing the physical Higgs scalar h, the derivative operator induces the following term In effect, the kinetic term of the Higgs field is modified to where ξ ≡ v 2 /f 2 . This means that the Higgs field should be rescaled by h → ζh, where ζ = (1 + c H ξ) −1/2 . Eq. (A.1) induces two further derivative operators which translate into the relations for κ 5 and κ 6 in Table 1.
To find the relations for λ 3 and λ 4 , we have to consider the Higgs potential amended by the c 6 term: In this case, the VEV is given by −µ 2 + 2λv 2 + 3 4 c 6 ξλv 2 = 0, (A. 6) and the corresponding Higgs mass is defined by After combining Eq. (A.6) and Eq. (A.7) and rescaling the Higgs field, we obtain a modified Higgs mass With these definitions of Higgs field and Higgs mass, we can write down the h 3 and h 4 terms: which yield the relations for λ 3 and λ 4 in Table 1.
Finally, we consider the operator This term modifies the fermion mass by After this redefinition, we obtain the following Higgs-fermion interaction operators which yield the relations for a 1 , a 2 and a 3 in Table 1. An alternative way of deriving such relations is performing a non-linear transformation We compare the results of both approaches, up to O(c H ξ) terms, in Table 14 (the coefficients of all other effective operators are set to zero). Despite these differences, both transformations necessarily yield the same cross section up to O(c H ξ). However, in the non-linear transformation approach, higher-order terms such as O(ξ 2 ) are much more complex than in the rescaling approach, and we have to deal with vertices such as tthh and tthhh even if c y is zero. Therefore, we adopt rescaling rather than the non-linear transformation method for defining our phenomenological parameters.
Rescaling Non-linear B Numerical cross sections of gg → h, gg → hh, and gg → hhh The cross section for gg → h can be put as where the integrated form factors F h i and the coefficients C i,h are given in Table 15, and K h denotes the K-factor. The unit of F h i is pb.
It is found that values of F h i given in Table (15) do produce a positive definite cross section of gg → h.
The cross section for gg → hh at 14 TeV LHC and a 100 TeV collider can be put as where the integrated form factors F 2h i and the coefficients C i,2h are given in Table 16 and Table 17, and K 2h denotes the K-factor, which is equal to 2.20 for the LHC 14 TeV and 2.17 for the 100 TeV collision, respectively. The unit of F 2h i in these two tables is fb.   In order to guarantee the positive and definite results of the cross section of all points in the parameter space, the contribution of b quark should be removed from the diagrams. Otherwise, a more general parameterisation of the cross section should be introduced. Furthermore, we have used more than 5,000 points in the parameter space of a 2 , c 2 , κ 5 , and λ 3 to determine these F 2h i after taking into account the constraints on parameters a 1 and c 1 from the projected precision in the measurement of σ(gg → h). The positivity and definiteness of the cross sections are examined to be hold in a random scan in the parameter space of a 2 , c 2 , κ 5 , and λ 3 with a total number of points 10 7 . If a 1 and c 1 can significantly deviate from the values of the SM, these results might not be valid anymore.
The cross section for gg → hhh at a 100 TeV collider can be put as  where the integrated form factors F 3h i and the coefficients C i,3h are given in Table 18 and  Table 19, and K denotes the K-factor which is taken as 2.1. The unit of F 3h i is fb. We have used more than 12,000 points to determine these F 3h i . The largest absolute coefficient is F 3h 30 . In contrast, the smallest absolute coefficients are F 3h 81 and F 3h 83 . After taking into account the constraints on parameters a 1 and c 1 from the projected precision data of σ(gg → h) and the constraints on parameters a 2 , c 2 , λ 3 and κ 5 from the projected precision data of σ(gg → hh), the positivity and definiteness of the cross sections are examined to be hold in a random scan in the parameter space of a 3 , λ 4 and κ 6 with a total number of points 10 7 .