Multi-Higgs Boson Production with Anomalous Interactions at Current and Future Proton Colliders

We investigate multi-Higgs boson production at proton colliders, in a framework involving anomalous interactions, focusing on triple Higgs boson production. We consider modifications to the Higgs boson self-couplings, to the Yukawa interactions, as well as new contact interactions of Higgs bosons with either quarks or gluons. To this end, we have developed a MadGraph5_aMC@NLO loop model, publicly available at https://gitlab.com/apapaefs/multihiggs_loop_sm, designed to incorporate the relevant operators in the production of multiple Higgs bosons (and beyond). We have performed cross section fits at various energies over the anomalous interactions, and have derived constraints on the most relevant anomalous coefficients, through detailed phenomenological analyses at proton-proton collision energies of 13.6 TeV and 100 TeV, through the 6 $b$-jet final state.


Introduction
Measuring the interactions of the Higgs boson [2][3][4][5][6] with the rest of the Standard Model (SM), with itself (i.e. its self-interactions), and with potential new phenomena, is a crucial task for future collider experiments.There are multiple theoretical justifications for the importance of the Higgs field and how it interacts; for example, the Higgs doublet bilinear, H : H, is the only dimension-2 (D " 2) singlet operator within the SM.This kind of operator readily allows for D " 3 or D " 4 interactions with one, or two, new scalar singlet fields, that may constitute part of a "hidden sector" of hitherto undetected particles.Therefore, the Higgs boson could be our primary window to this, potentially "dark", sector, that could comprise part, or all, of dark matter.Another reason of the importance of the Higgs boson's interactions is that they may illuminate the mechanism of baryogenesis -in particular, if this process occurred during the electroweak phase transition, the period during which the Higgs field acquires a vacuum expectation value (VEV), then the Higgs field may be one of the protagonists of electroweak baryogenesis (see, e.g.[7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]).A central rôle in this mechanism is played by the Higgs potential, which, within the SM, is introduced in a rather ad-hoc manner.The form of this potential can, at least in principle, be directly probed by measuring the Higgs boson's self-interactions through multi-Higgs boson production processes: the triple self-coupling can be probed through Higgs boson pair production, and the quartic through triple production.
Novel Higgs boson interactions, to new or to SM particles, may arise at the electroweak scale, but they may also appear at higher scales, Opfew TeVq.If this is the case, we can parametrize our ignorance using a higher-dimensional effective field theory (EFT), see, e.g.[22][23][24].Neglecting lepton-number violating operators, the lowest-dimensionality EFT that can be written down consists of D " 6 operators.Upon electroweak symmetry breaking, the Higgs boson would acquire a VEV and these operators would result in several new interactions, as well modifications of the SM interactions of the Higgs boson.Several of the operators in the physical basis of the Higgs boson scalar (h) would then have coefficients that are correlated, according to D " 6 EFT.These correlations, however, may be broken by even higher-dimensional operators (e.g.D " 8), particularly if the new phenomena are closer to the electroweak scale.Therefore, it may be beneficial to lean towards a more agnostic, and hence more phenomenological, approach and, while still remaining inspired by D " 6 EFT, consider fully uncorrelated, "anomalous" interactions of the Higgs boson with the SM.These kind of interactions can also be achieved within the framework of the electroweak chiral Lagrangian including a light Higgs boson, where free coefficients for all Higgs couplings are obtained [25][26][27][28][29][30][31][32][32][33][34][35].This framework provides a consistent EFT justification of the usual "κ" formalism.The description of the Higgs boson's interactions in terms of uncorrelated coefficients is the approach that we pursue here.
Within the realm of self-coupling measurements, and beyond, Higgs boson pair production has received considerable attention from both the experimental (e.g. ) and theoretical (e.g. ) point of view, particularly following the discovery of the Higgs boson.Owing to its much smaller SM cross section at hadron colliders [90], triple Higgs boson has understandably received much less attention.To the extent of our knowledge, it was first investigated in [125], where it was demonstrated that the prospects for the measurement of the quartic self-coupling will be very challenging, both at the LHC and a future 'VLHC' at a center-of-mass energy of 200 TeV.Subsequent studies of triple Higgs boson production at future colliders [126][127][128][129][130][131][132][133][134][135], with the knowledge of the value of the Higgs boson mass, and the prospects for Higgs boson pair production measurements at hand, further quantified the difficulty of this measurement, nevertheless demonstrating that some useful information may be obtained via the process.The situation at LHC energies is much more dire, as expected, requiring either new phenomena for sufficient enhancement [136], or very large new contributions to the Higgs boson's self-interactions [137].Up to now, all of the phenomenological studies addressing triple Higgs boson production, apart from discussions of the cross section modifications in an EFT in [138], have considered only the effects of modifications of the self-couplings.However, due to their loop-induced nature, proceeding via heavy quark loops, multi-Higgs boson production processes are highly sensitive to modifications of the heavy quark Yukawa couplings, as well as to additional quarkor gluon-Higgs boson contact interactions that generically arise in EFTs.Such effects have been studied within the context of Higgs boson pair production, see, e.g.[139,140].In the present article we study, for the first time to the best of our knowledge, the combined effect of modifications to the self-couplings, Yukawa interactions, and new heavy-quark-or gluon-Higgs boson contact interactions in triple Higgs boson production, within the aforementioned EFT-inspired anomalous coupling picture.To achieve this, we have constructed a specialized anomalous-coupling model within the MadGraph5 aMC@NLO event generator that incorporates the full interference effects between contributing diagrams.
The paper is organized as follows: in section 2 we introduce the phenomenological Lagrangian that we employ in our study, inspired by an EFT extension of the SM.In section 3 we discuss the details of the Monte Carlo implementation of the MadGraph5 aMC@NLO model, briefly addressing its usage for Higgs boson pair and triple production.In section 4 we discuss our phenomenological analysis of triple Higgs boson production, focusing on cross section fits, existing constraints on the anomalous couplings, and the extraction of limits on the most relevant coefficients through the 6 b-jet final state.We present our conclusions and outlook in section 5.In addition, in appendix A we describe the D " 6 EFT Lagrangian relevant to Higgs boson production processes, i.e. the the linear EFT approach, and in appendix B we discuss the nonlinear approach through the electroweak chiral Lagrangian.In appendix C we discuss aspects of the validation of the Monte Carlo model.In appendix D we discuss relevant Feynman diagrams that appear within our framework, up to two operator insertions.Finally, in appendix E we discuss limits obtained with up to three powers of anomalous operators (at the matrix element-squared level) instead of the two presented in the main text, and in appendix F we present supplementary cross section fits at various proton collision energies.

Phenomenological Lagrangian
To allow for additional freedom, we have modified the D " 6 EFT Lagrangian relevant to Higgs boson processes, derived in detail in appendix A (eq. A.11), breaking the correlations between several of the interactions that originate in the D " 6 EFT.This yields more flexibility when exploring multi-Higgs boson production experimentally.Nevertheless, D " 6 EFT can be trivially restored by reinstating the appropriate relations between the coefficients (see below).We also provide additional EFT justification for this approach, inspired by the electroweak chiral Lagrangian including a light Higgs boson in appendix B.
To achieve this, we have defined new anomalous couplings d 3 and d 4 , as modifications of the triple and quartic self-couplings, c g1 and c g2 as interactions of two gluons and one or two Higgs bosons, respectively, and c f 1 , c f 2 and c f 3 as interactions between fermion-antifermion pairs, and one, two or three Higgs bosons, respectively, leading to the following phenomenological Lagrangian: To recover the D " 6 EFT Lagrangian relevant to the Higgs sector (i.e.eq.A.11), one would need to simply set: The implementation of this study further modifies the above phenomenological Lagrangian to match more closely the LHC experimental collaboration definitions: where we have taken λ SM " m 2 h {2v 2 .The CMS parametrization is then obtained by setting: κ λ " p1`d 3 q, k t " c t1 , c 2 " c t2 , c g " c g1 , c gg " c 2g and the ATLAS parametrization by c hhh " p1 `d3 q, c ggh " 2c g1 {3, c gghh " ´cg2 {3.The Lagrangian of eq.2.2 encapsulates the form of the interactions that we employ for the rest of our phenomenological analysis.

Loop ˆTree-Level Interference
The implementation of the Lagrangian of eq.2.2 in MadGraph5 aMC@NLO (MG5 aMC) [141,142] follows closely the instructions for proposed code modifications found in [143]. 1 These modifications essentially introduce tree-level diagrams in the form of "UV counter-terms", that are generated along with any loop-level diagrams, allowing the calculation of interference terms between them.The model of the present article, created through this procedure, has been fully validated by direct comparison to an implementation of Higgs boson pair production in D " 6 EFT in the HERWIG 7 Monte Carlo, and by taking the limit of a heavy scalar boson for those vertices that do not appear in that process.See appendix C for further details of the latter effort.The necessary modifications to the MG5 aMC codebase, 2 as well as the model can be found in the public gitlab repository at [1]. 3

Higgs Boson Pair Production
Higgs boson pair production contains a subset of all diagrams relevant to triple Higgs boson production, minus one Higgs boson insertion, and therefore we begin by investigating this process, so as to provide a simple example.
MG5 aMC allows for calculation of the interference terms via one insertion of the operator at the squared matrix-element level (|M| 2 ): i.e. schematically: |M| 2 " ‚1 `‚c i .The classes of diagrams that are included in addition to the SM ones (fig.1), are those that appear in fig.8, found in appendix D.
Alternatively, one can obtain: the interference terms of one or two insertions with the SM, and the squares of one-insertion terms (equivalent to two powers of the anomalous couplings) via: where schematically this would correspond to |M| 2 " ‚1 `‚c i `‹c i c j .The classes of diagrams of fig. 9 (appendix D) appear in addition.

Triple Higgs Boson Production
For the investigation of triple Higgs boson production, which constitutes the focus of this article, we included the interference terms of one or two insertions with the SM and the squares of one-insertion terms through:

generate g g > h h h [noborn=MHEFT QCD] MHEFT^2<=2
The results with up to three powers of the anomalous operators, more appropriate within the framework of the electroweak chiral Lagrangian (appendix B), are presented in appendix E.
As in the case of Higgs boson pair production, our current treatment would schematically correspond to |M| 2 " ‚1 `‚c i `‹c i c j .The classes of the various diagrams that appear are shown in fig. 2 for the SM, and in figs.10 for one insertion, and 11 for two insertions, in appendix D.

Cross Section Fits
By employing our Monte Carlo implementation of eq.2.2, We have performed cross section fits at a proton-proton collider at 13, 13.6, 14, 27, and 100 TeV, over all coefficients relevant to either Higgs boson pair or triple production.Two-dimensional projections of the fits for E CM " 13.6 TeV appear in fig. 3 for triple production.In this figure, for each corresponding plot, all other coefficients have been set to zero, for the sake of simplicity.
The numerical values of the polynomial coefficients, yielding the cross section in the form σ{σ SM ´1 "  (gg hhh)@13.6TeV, normalized to SM value Figure 3: Fit of the cross section for triple Higgs boson production at 13.6 TeV, normalized to the SM value.For each combination of couplings, the other couplings have been set to zero for simplicity.Each change of color in the contours represents a shift of a factor of 0.5ˆthe SM value. is ´0.0703, i.e. read off the second row, third column.Using these coefficients, one can construct the cross section for any given value of the anomalous couplings.
All the fits for the signal processes, and subsequent simulations, have been performed using the MSHT20nlo_as118 PDF set [144] and the default dynamical scale choice (option 3) in MG5 aMC, which corresponds to the sum of the transverse mass divided by 2.
We note here that the expected contribution of bottom-quark loops in the SM, both at LHC energies and at 100 TeV, is expected to be Op0.1%q.Therefore, anomalous couplings of the Higgs boson to the bottom quark are included merely for completeness, and furthermore, this ensures that we can safely neglect charm quark contributions in our analysis.

Other Constraints on Anomalous Couplings
The majority of the anomalous couplings that appear in the phenomenological Lagrangian of eq.2.2 are already tightly constrained by other processes that involve the interactions of gluons, top and bottom quarks with the Higgs boson.The two exceptions that are not presently constrained are the anomalous interactions of three Higgs bosons and two top or bottom quarks, with relevant coefficients c t3 and c b3 , as well as the anomalous modification to the Higgs boson's quartic interaction, related to the d 4 coefficient.While it is beyond the scope of the present study to perform a full fit, involving several processes and constraints, with their associated correlations, it is important to provide order-of-magnitude estimates for the two scenarios that we examine: the 13.6 TeV high-luminosity LHC, and the 100 TeV FCC-hh at the end of its lifetime.
We consider the constraints on c t1 , c b1 and c g1 that would arise within the "kappa-0" scenario, as they are defined in [145] (table 3).For the HL-LHC we consider those labeled "HL-LHC", and for the 100 TeV FCC-hh, we consider those projected after the combined results of "FCC-ee240+FCC-ee365", "FCC-eh", and "FCC-hh" for the FCChh, collectively labeled "FCC-ee/eh/hh".For the modifications to the Higgs boson triple self-coupling, we use the values quoted for the "di-Higgs exclusive" results for the "HL-LHC", shown in table 12 of [145].For c g2 and c t2 , we assume approximate constraints derived from the results of the analysis of [140].Due to the lack of concrete studies in the literature for c b2 , to the best of our knowledge, we have assumed that the constraint is identical to that of c t2 , obtained from [140].For all the coefficients, we assume that the central value of these future measurements will be at zero.The assumed constraints are summarized as percentage uncertainties in table 2, along with the corresponding references, for convenience.

Statistical Analysis
The goal of the present phenomenological analysis is two-fold: on the one hand, we wish to determine if triple Higgs boson production itself can be observed at the end of the HL-LHC or FCC-hh lifetime, and on the other, we wish to calculate the constraints that would be obtained on the anomalous couplings, given the assumption the SM is the true underlying theory.The first goal corresponds to the null hypothesis being the absence of triple Higgs boson production, whereas the latter corresponds to the null hypothesis being SM triple Higgs boson production, i.e. triple production with all the anomalous coupling coefficients set to zero.
For the first of these goals, i.e. to determine the significance, Σ, of observing triple Higgs boson production for a certain anomalous coupling parameter-space point, we optimize the signal versus background discrimination according to our phenomenological analysis (described in the following section, 4.4), to obtain the maximum significance over the set of cuts for SM triple Higgs boson production.Since we expect the number of signal events to be much lower than that for the background, we employ the following formula for the significance: where S and B are the expected number of signal and background events, respectively, at a given integrated luminosity L, and α is a factor that we employ to model the systematic uncertainty present in the background estimation.During our search for the maximum SM triple Higgs boson production significance, we used the value α " 0, i.e. we optimized the quantity S{ ?B. Subsequently, all results are presented with this optimal set of cuts, obtained separately for each collider, while setting the systematic uncertainty factor α " 0.05.
We have derived the expected "evidence" (3σ) and "discovery" (5σ) limits on the observation of triple Higgs boson production on the pc t3 , d 4 q-plane, where the novel information coming from triple Higgs boson production is expected to arise.Furthermore, we have derived the 68% confidence level (C.L.) (1σ) and 95% C.L. (2σ) expected limits on the same plane, given that the SM hypothesis is true (i.e.SM triple Higgs boson production).While doing so, we allow the coefficients d 3 and c t2 to vary.These will be constrained primarily via Higgs boson pair production, and since this is a challenging process in itself, will not be substantially constrained compared to the gluon-Higgs anomalous interactions, or the top-anti-top-Higgs anomalous coupling, c t1 , see table 2 for a quantification of this statement.In addition, we do not expect triple Higgs boson production to be competitive to Higgs boson pair production in terms of constraints on d 3 and c t2 , so we "marginalize" over these two.Note that, due to the fact that we are assuming minimal flavour violation, as discussed in appendix A, as well as for the sake of simplicity, we do not expect the effect of the c b3 coefficient to be as significant as that of c t3 , and therefore we have set c b3 to zero in the phenomenological analysis.
To derive the evidence and discovery limits for triple Higgs boson production, we assume that the number of signal and background events follow gaussian distributions, and define the uncertainty on the background number of events as: Then, the probability value (p-value) to obtain a given number of signal events, Sptc i uq, corresponding to the set of anomalous coupling coefficients tc i u, is given by: To take into account the approximate constraints of table 2, that would be present, ignoring possible correlations that will exist between constraints, we multiply the p-values by a gaussian prior factor: at a value c i of each of the anomalous coefficients d 3 and c t2 .We then marginalize by integrating over the d 3 and c t2 directions.Since this is done over a grid, the integral is represented by a sum over the grid points: where ∆d 3 and ∆c t2 are the bin widths in the d 3 and c t2 directions, respectively.The probability is then converted to χ 2 via the Python χ 2 inverse survival function for two degrees of freedom (scipy.stats.chi2.isf),and the 3σ and 5σ limits are found by drawing the iso-contours corresponding to χ 2 ´χ2 min " p11.8, 28.7q, respectively.The results are shown in Fig. 4, in section 4.5.
For the determination of the constraints obtained for the null hypothesis being SM triple Higgs boson production, we folow a similar procedure.In this scenario, we assume that the number of events for both signal and background follow a gaussian distribution, as we did before, and hence the uncertainty on the number of total expected SM events is given by: where S SM represents the expected number of events for SM triple Higgs boson production.Then the p-value to obtain a given number of events corresponding to an anomalous coupling parameter-space point, Sptc i uq, given that the SM is the truth, is given by: This p-value is then multiplied by the prior factors of eq.4.4, and summed as before: PM pc t3 , d 4 q " ÿ pd 3 ,c t2 q ∆d 3 ∆c t2 P C pd 3 qP C pc t2 q P ptc i uq , (4.8) followed by a conversion to the χ 2 value as before.The 68% (1σ) and 95% (2σ) confidence level (C.L.) intervals are found by drawing the iso-contours corresponding to χ 2 ´χ2 min " p2.28, 5.99q, respectively.The results of this analysis are shown in Fig. 5, in section 4.5.
To obtain the one-dimensional limits on either of the d 4 or c t3 couplings alone, we take the procedure one step further: we sum over the additional marginalized coupling as in eq.4.8, and convert the p-value into a χ 2 value.For the triple Higgs boson searches, we calculate the 3σ and 5σ evidence and discovery limits by finding the points that correspond to χ 2 ´χ2 min " p9.00, 25.0q, respectively.For the expected constraints on the d 4 and c t3 coefficients, given that the SM is true, we calculate the 1σ and 2σ intervals by finding the points that correspond to χ 2 ´χ2 min " p0.99, 3.84q, respectively.The results for the evidence and discovery points, and for the expected limits given that the SM is true, are shown in tables 5 and 6, respectively, in section 4. 5.
In regard to the use of the Gaussian approximation in our analysis, we note here that, while the expected number of events for the case of the SM signal is low and, strictly speaking, Poisson statistics should be employed in that case, we would like to emphasize that this is not the case at the 1, 2, 3 and 5σ exclusion/discovery boundaries that we have deduced.To be concrete, for the HL-LHC case, to obtain a 1σ significance, we would need, given our estimates for the expected number of background events (table 4, top) B " 419.6, the following number of signal events: For this expected number of signal events, the Gaussian approximation would yield a 68% confidence-level interval of 29.32 ˘5.41, whereas the Poisson approximation would yield 29.32 `4.48  ´6.48 .Given the other uncertainties that enter our phenomenological analysis, it is therefore reasonable to assume that events are Gaussian-distributed at the boundaries of the exclusion/discovery regions.

Phenomenological Analysis
To obtain constraints on anomalous triple Higgs boson production at proton colliders, following the statistical procedure outlined above, we have performed a hadron-level phenomenological analysis of the 6 b-jet final-state originating from the decays of all three Higgs bosons to b b quark pairs.We closely follow the analysis of Refs.[132,136].Partonlevel events have been generated using the MG5 aMC anomalous couplings implementation constructed in the present article, with showering, hadronization, and simulation of the underlying event, performed via the general-purpose HERWIG 7 Monte Carlo event generator [146][147][148][149][150][151][152][153].The event analysis was performed via the HwSim framework addon to HERWIG 7 [154].No smearing due to the detector resolution or identification efficiencies have been applied to the final objects used in the analysis, apart from a b-jet identification efficiency, discussed below.We expect such effects to produce differences in binned jet observables of Op10%q at the LHC, see for example [155,156], and anticipate substantial improvement to this at the FCC-hh.
The branching ratio of h Ñ b b will be modified primarily due to c t1 , c g1 , indirectly through modifications to the h Ñ gg and h Ñ γγ branching ratios, and directly through c b1 .
To take this effect into account, we employed the eHDECAY code [157].The program eHDECAY includes QCD radiative corrections, and next-to-leading order electroweak corrections are only applied to the SM contributions.For further details, see Ref. [157].We have performed a fit of the eHDECAY branching ratio h Ñ b b, and we have subsequently normalized this to the latest branching ratio provided by the Higgs Cross Section Working Group's Yellow Report [158,159], BRph Ñ b bq " 0.5824.The fit is then used to rescale the final cross section of pp Ñ hhh Ñ pb bqpb bqpb bq.The background processes containing Higgs bosons turned out to be subdominant with respect to the dominant QCD 6 b-jet and Z+jets backgrounds, and therefore we did not modify these when deriving the final cross sections.
For the generation of the backgrounds involving b-quarks not originating from either a Z or Higgs boson, we imposed the following generation-level cuts for the 100 TeV proton collider: p T,b ą 30 GeV, |η j | ă 5.0, and ∆R b,b ą 0.2.The transverse momentum cut was lowered to p T,b ą 20 GeV for 13.6 TeV, except for the QCD 6 b-jet background, for which we produced the events inclusively, without any generation cuts. 5The selection analysis was optimized considering as a main backgrounds the QCD-induced process pp Ñ pb bqpb bqpb bq, and the Z+jets process (represented by Z `pb bqpb bq), which we found to be significant at LHC energies.
The event selection procedure for our analyses proceeds as follows: as in [132], an event is considered if it contains at least six b-tagged jets, of which only the six ones with the highest p T are taken into account.A universal minimal threshold for the transverse momentum, p T,b , of any of the selected b-tagged jets is imposed.In addition a universal cut on their maximum pseudo-rapidity, |η b |, is also applied.We subsequently make use of the observable: r8, 8, 8s GeV ∆R bb ph i q ă r3.5, 3.5, 3.5s r3.5, 3.5, 3.5s ∆Rph i , h j q ă r3.5, 3.5, 3.5s r3.5, 3.5, 3.5s p T ph i q ą r0.0, 0.0, 0.0s GeV r200.0,190.0, 20.0s GeV Table 3: Optimized cuts determined for the phenomenological analysis.The indices i, j can take the values i, j " 1, 2, 3.For the cut ∆Rph i , h j q the three pairings correspond to p1, 2q, p1, 3q, p2, 3q.The indexed elements should be read from left to right in increasing order.The last two rows refer to cuts over light jets.
We further refine the discrimination power of the χ 2,p6q min variable by using the individual mass differences ∆m " |m qr ´mh | in eq.(4.10), sorting them out according to ∆m min ă ∆m med ă ∆m max , and imposing independent cuts on each of them.We also consider the transverse momentum p T ph i q of each reconstructed Higgs boson candidate.These reconstructed particles are also sorted based on the value of p T ph i q, on which we then impose a cut.Besides the universal minimal threshold on p T,b , introduced at the beginning of this section, we impose cuts on the three b-jets with the highest transverse momentum p T,b i , for i " 1, 2, 3.The set of cuts p T,b 3 ă p T,b 2 ă p T,b 1 is the second most powerful discriminating observable in our list.Finally, we also considered two additional geometrical observables.The first of them is the distance between b-jets in each reconstructed Higgs boson ∆R bb ph i q.The second one is the distance between the reconstructed Higgs bosons ∆Rph i , h j q themselves.
Our optimization process then proceeds as in [132,136]: we sequentially try different combination of cuts over the observables introduced above on our signal and background samples until we achieve a significance above 2 or when our number of Monte Carlo events is reduced so drastically that no meaningful statistical conclusions can be derived if this number becomes smaller (this happens for instance when for a given combination of cuts, we are left with less than 10 Monte Carlo events of signal or background).
For 100 TeV, the observables are optimized in the following order, (i) p T,b , (ii) |η b |, (iii) χ 2,p6q , (iv) ∆m min , ∆m med , ∆m max , (v) p T,b 1 and (vi) p T,b 2 .By testing different cut combinations over the variables above we reach a SM triple Higgs boson significance of S{ ?B » 2.33σ without b-jet tagging or S{ ?B » 1.43σ including the b-tagging factor (0.85).We note that this is an improvement over the previous result of [132], where the corresponding significance was around 1.7σ without b-tagging.For 13.6 TeV the optimization of the cuts follows a similar path, however after applying the cut on χ

5.07
Table 4: The lists of processes considered during our phenomenological analysis, along with their respective cross sections to the 6 b-jet final state.The efficiencies ε analysis and number of events N cuts L , correspond to those obtained after applying the set of cuts given in table 3. A b-jet identification efficiency of 0.85 (for each b-jet) has also been applied to obtain the number of events.For the HL-LHC we considered an integrated luminosity of L " 3000 fb ´1, and for for the FCC-hh a luminosity of L " 20 ab ´1.To approximate higher-order corrections, a K-factor K " 2 has been included for all processes, with respect to the leading-order cross section.The background processes marked with "LI" represent loop-induced contributions that have been generated separately.
with Op10q Monte Carlo events for background, and therefore the full procedure stops at a significance of only 0.23σ without b-jet tagging or 0.14σ with b-jet tagging (0.85).The concrete values of the cuts applied on the different kinematic variables depend on the center of mass energy of the collider under study, and are summarized on table 3. We note that the optimal set of cuts for the HL-LHC requires more central b-jets than that at the FCC-hh.The complete list of backgrounds, their initial cross sections, efficiency of cuts and expected number of events at each of the full collider luminosities, are given in table 4. It is interesting to note that the Z `pb bqpb bq is in fact the dominant one at the HL-LHC, whereas the reverse is true for the FCC-hh.It will be important to validate this result in a future detailed experimental study that consider the full effects of the detectors.
To employ the results of the SM analysis over the whole of the parameter space we are considering, we have performed a polynomial fit of the efficiency of the analysis on the signal, ε analysis phhhq, at various, randomly-chosen, combinations of anomalous coefficient values.In combination with the fits of the cross section, and the fit of the branching ratio of the Higgs boson to pb bq, we can estimate the number of events at a given luminosity, for a given collider for any parameter-space point within the anomalous coupling picture, which we dubbed Sptc i uq in our discussion of the statistical approach we take, in section 4.3.The main results of our two-dimensional analysis over the pc t3 , d 4 q-plane are shown in figs.4 and 5.In particular, fig. 4 shows the potential "evidence" and "discovery" regions (3σ and 5σ, respectively) for triple Higgs boson production at the high-luminosity LHC on the left (13.6 TeV with 3000 fb ´1 of integrated luminosity), and at a FCC-hh (100 TeV, with 20 ab ´1) on the right.Evidently, very large modifications to the quartic self-coupling are necessary for discovery of triple Higgs boson production at the HL-LHC, ranging from gg hhh@100 TeV, L=20000 fb 1 , syst.= 5.0%The one-dimensional analysis' results, presented in tables 5 and 6, for the "evidence" and "discovery" potential, and exclusion limits, respectively, reflect the above conclusions.For instance, it is clear by examining table 5, that the HL-LHC will only see evidence of triple Higgs boson production in the 6 b-jet final state only if d 4 has modifications of |d 4 | " Opfew 10q, and will only discover it if |d 4 | " Op100q.On the other hand, there could be evidence or discovery of Higgs boson triple production if |c t3 | " Op1 ´10q.The 1σ and 2σ exclusion regions are much tighter, as expected, with |d 4 | " Op10q at 1σ or 2σ at the HL-LHC, improving somewhat at the FCC-hh, and |c t3 | " Op0.1 ´1q, both at the HL-LHC and FCC-hh.

Conclusions
We have conducted a phenomenological analysis of the triple Higgs boson production at the high-luminosity LHC (13.6 TeV with 3000 fb ´1) and a proposed future circular collider colliding protons at 100 TeV (20 ab ´1).For this analysis, we have constructed a model to be used with the MadGraph5 aMC@NLO event generator, which is publicly available at [1].By examining the 6 b-jet final state, within a signal versus background analysis, we have concluded that interesting constraints can be obtained on the most relevant coefficients, d 4 and c t3 , representing modifications to the Higgs quartic self-coupling, and additional top-Higgs interactions, respectively.These are presented, for the HL-LHC and the FCC-hh, over the pc t3 , d 4 q-plane in figs.4 and 5, showing the evidence/discovery potential of the triple Higgs boson production process itself, marginalized over the d 3 and c t2 coefficients, and the 1/2σ exclusion regions arising from the process on that plane.The one-dimensional evidence/discovery regions over either d 4 , or c t3 at both colliders are given in table 5, and the possible constraints extracted in table 6.
The results of our study demonstrate the importance of including additional contributions, beyond the modifications to the self-couplings, when examining multi-Higgs boson production processes, and in particular triple Higgs boson production.We are looking forward to a more detailed study for the HL-LHC, conducted by the ATLAS and CMS collaborations, including detector simulation effects, and the full correlation between other channels.From the phenomenological point of view, improvements will arise by including additional final states, e.g.targetting the process pp Ñ hhh Ñ pb bqpb bqpτ `τ ´q, or by performing an analysis that leverages machine learning techniques to maximize significance. 6n summary, we believe that the triple Higgs boson production process should constitute part of a full multi-dimensional fit, within the anomalous couplings picture.

A D " 6 Effective Field Theory
The phenomenological Lagrangian employed in this study can be motivated by considering the D " 6 effective field theory (EFT) Lagrangian of ref. [139].The discussion that follows has been adapted for the purposes of the present article.
New Physics associated to a new scale Λ " v can be described in a model-independent way by augmenting the Lagrangian of the SM with all possible gauge-invariant operators of mass dimension D ą 4, where the leading effects arise from D " 6 operators (neglecting lepton-number violating operators, irrelevant to our study).Working at this level, the extension of the SM that is relevant for multi-Higgs boson production reads: where α s is the strong coupling constant and α 1 " g 1 2 {4π.The full set of D " 6 operators that can be formed out of the SM field content was first obtained in [22] and reduced to a non-redundant minimal set in [23].Here, we employed equations of motion to move to the basis used in [24,160] and then imposed constraints from precision tests to neglect a class of operators whose effect is already constrained to be at most 1% with respect to the SM, following [161][162][163][164][165]. Including these operators would have a negligible numerical impact on the analysis, given the experimental and theoretical errors.
Precision measurements also lead to the approximate restrictions [24] c HB 16π 2 " ´cHW 16π 2 " ´cB " c W , (A.2) which we will employ in the following.Thus our setup corresponds to a restricted stronglyinteracting light Higgs (SILH) Lagrangian [166] where c T has been set to zero and the relations (A.2) are used.
We note here that we do not include the effects of the chromomagnetic dipole operator, which has been shown be significant in processes such as pp Ñ t t, pp Ñ h, pp Ñ t th, [167][168][169] as well as Higgs boson pair production [170].Nevertheless, the chromomagnetic operator is subleading from the point of view of weakly interacting UV dynamics [171], and we therefore leave it to future work.
The terms contributing to multi-Higgs production via gluon fusion at the LHC at D " 6 EFT are: where the first line includes the relevant Standard Model terms that will receive corrections from D " 6 operators.For the purposes of the current study, we set C H " 0 in what follows.We also assume minimal flavour violation [172], which leads to the coefficients of the Yukawa-like terms in the last row of eq.A.3 being proportional to the (SM-like) Yukawa couplings.
Examining the relevant terms in the scalar potential and expanding about the electroweak minimum, we get: pv `hq 6  8 Omitting terms with h n , n ą 4, and constant terms we arrive at Finally, we focus on the fermion-Higgs boson interactions that receive contributions from where f " t, b, ..., with f L,R the left-and right-handed fields, and the first term comes from the SM whereas the second term is a dimension-6 contribution.Expanding, we obtain The first line gives the expression for the modified fermion mass, and we can re-express eq.(A.7) in terms of this: The final term that we need to consider is where the omitted constant term can be absorbed into an unobservable re-definition of the gluon wave function.
We thus obtain the following interactions in terms of the Higgs boson scalar h, relevant to Higgs boson multi-production: where we have explicitly written down the contributing components of the Q L doublets.

B The Electroweak Chiral Lagrangian Including a Light Higgs Boson
The D " 6 Lagrangian of eq.A.3 represents a "linear" realization of an effective field theory, including the leading terms in terms of "canonical dimensions" D, and assuming baryon-and lepton-number conservation.This is also sometimes referred to as "SMEFT".An alternative way to organize the EFT is by chiral dimensions, which yields a "nonlinear" realization.
In this generalization of the SM, the gauge interactions are unchanged (at leading order), but general anomalous couplings are introduced for the physical Higgs boson.To do this in a consistent, gauge-invariant way, the scalar fields of the theory have to be decomposed into the three Goldstone fields φ a , described by: where and where T a are the generators of SUp2q.The polar decomposition of the Higgs doublet H is then given by: where h is the physical Higgs boson.Under electroweak gauge transformations SUp2q L ˆUp1q Y : such that the h is invariant, and its couplings can be consistently modified.The transofmrations g L , and the Up1q Y subgroup of g R are gauged, so that the covariant derivatives are given by:7 New dynamics may appear in the Higgs sector at a new physics scale f " OpTeVq.This could be, e.g.composite, pseudo-Goldstone Higgs particles [173][174][175][176][177], but also by other models with a modified Higgs sector with either weak or strong couplings.We can then parametrize deviations of the SM couplings by a quantity ξ " v 2 {f 2 , where v " 246 GeV is the electroweak scale.Anomalous contributions, with respect to the SM, of order ξ in the Higgs couplings will generically lead to a cut-off Λ " 4πf in the effective description of the new Higgs dynamics.This picture might be supplemented by TeV-scale (order f ) new degrees of freedom (non-standard fermions, extra pseudo-Goldstone bosons), integrated out in the EFT at the electroweak scale v.The EFT can then be organized in full generality [178] as a double expansion in ξ " v 2 {f 2 and f 2 {Λ 2 " 1{16π 2 , which are the two dimensionless parameters that can be formed out of the three relevant scales v, f and Λ.They are both small under the condition v !f !Λ.The expansion in ξ amounts to an expansion of the Lagrangian in operators of increasing canonical dimension.The expansion in f 2 {Λ 2 , corresponds to a loop expansion or, equivalently, to an expansion in terms of increasing chiral dimension (χ).To leading order in chiral dimensions, and to first order in ξ, the SM effective Lagrangian can be written in nonlinear notation as For generality, we have included generic flavor matrices Ȳf arising at NLO.In scenarios with minimal flavor violation [172], we can assume Ȳf 9Y f , as we did for the D " 6 EFT (i.e.linear) scenario.
When applying the chiral Lagrangian formalism to processes that arise only at oneloop level in the SM, such as h Ñ gg, h Ñ γγ or h Ñ Zγ, terms at NLO (chiral dimension 4, or loop order 1) become relevant, in addition to loops with modified couplings from eq.B.6.For the complete list of NLO operators we refer the reader to [34].The relevant terms here are " e 2 F µν F µν h, " eg 1 F µν Z µν h and " g 2 s ⟨G µν G µν ⟩ h.In summary, the terms from the electroweak chiral Lagrangian relevant to multi-Higgs boson production, up to gg Ñ hhh, are given by [159,179]: ı ´" m t v 3 pc ttt q tL t R h 3 `mb where we have again taken λ SM " m 2 h {2v 2 .The identifications that yield the phenomenological Lagrangian employed in this article, eq.2.2, are: Almost all vertices relevant in triple Higgs boson production are present in pair production.Therefore, these were validated using the HiggsPair model in HERWIG 7 Monte Carlo event generator, with minor modifications, in the gg Ñ hh process.Validation of the new t thhh contact interaction was performed by considering the processes t t Ñ hhh at 1 TeV in the t t center-of-mass frame without a PDF involved.The processes were calculated both in the HEFT and in a model with a heavy scalar (H) that couples to t t and hhh only.This implies taking the limit of the effective field theory directly and checking whether the effective vertex functions as expected.The matching of the coefficient of eq.2.2 with the singlet model, e.g. of [18], implies that c t3 " 2v 2 {M 2 H , when the the quartic coupling between the heavy scalar and the three Higgs bosons is set to λ 1112 " 1 and the mixing angle θ " π{2 such that the SM Higgs boson is decoupled.Figure 7 shows the ratio of the anomalous t thhh interaction cross section over the heavy scalar cross section for various masses of the heavy scalar, chosen to be much higher than the center-of-mass energy.

Cross Section Ratios
Validation of the tthhh contact interaction HEFT/SING@tt hhh@1 TeV

Figure 7:
The ratio of cross sections between for t t Ñ hhh between the anomalous interaction (HEFT) and the heavy scalar (H) descriptions.See main text for further details.We present here for completeness, the evidence/discovery boundaries and the 1σ and 2σ limits for the HL-LHC and FCC-hh, obtained with up to three powers of the anoma-HL-LHC 3σ HL-LHC 5σ

F Cross Section Fits for Triple and Double Higgs Boson Production at Various Energies
This appendix provides the coefficients of the polynomial function for the cross section σ{σ SM ´1 "

Figure 1 :
Figure 1: Example Feynman diagrams for leading-order gluon-fusion Higgs boson pair production in the Standard Model.

Figure 2 :
Figure 2: Example Feynman diagrams for leading-order gluon-fusion Higgs boson triple production in the Standard Model.

4 SM 4 SM 5 Figure 4 :
Figure4: The 3σ evidence (black solid) and 5σ discovery (red dashed) curves on the pc t3 , d 4 q-plane for triple Higgs boson production at 13 TeV/3000 fb ´1 (left), and 100 TeV/20 ab ´1 (right), marginalized over the c t2 and d 3 anomalous couplings.Note the differences in the axes ranges at each collider.

d 4 "
125 for c t3 " ´8, to d 4 " ˘40 for c t3 " 0 and then down to d 4 " ´200 for c t3 " 12.The situation is greatly improved, as expected, at the FCC-hh, where the range of d 4 is reduced to d 4 " 40 for c t3 " ´1.5, and to d 4 " ´20 for c t3 " 1.0.It is interesting to note that the whole of the parameter space with c t3 Á 1.0, or with c t3 À ´1.5 is discoverable, at the FCC-hh at 5σ.For the potential 68% (1σ) and 95% C.L. (2σ) constraints of fig.5, the situation is slightly more encouraging for the HL-LHC, with the whole region of d 4 Á 40 or d 4 À ´60 excluded at 95% C.L..The corresponding region at 68% C.L. is d 4 Á 20 and d 4 À ´30.For c t3 , it is evident that all the region c t3 À ´2 and c t3 Á 5 will be excluded at 95% C.L. and c t3 À ´1, c t3 Á 4 at 68% C.L.. On the other hand, the FCC-hh will almost be able to exclude the whole positive region of d 4 for any value of c t3 at 68% C.L..This will potentially be achievable if combined with other Higgs boson triple production final states.For the c t3 coupling, both the constraints reach the Opfew 10%q level for any value of d 4 .

Figure 6 :
Figure 6: The t t Ñ hhh process used to validate the implementation of the t thhh vertex.The process is shown in the model with a new heavy scalar (H) and in the limit of M H " ?ŝ.

Figure 8 :
Figure 8: Example Feynman diagrams with one EFT operator insertion contributing to Higgs boson pair production.

Figures 8 and 9 Figure 9 :Figure 10 :Figure 11 : 4 SM 5 Figure 12 :
Figures 8 and 9 represent the Feynman diagrams for either one or two insertions of the operators used in the present article in Higgs boson pair production.Figures 10 and 11 represent the Feynman diagrams for either one or two insertions in the context of Higgs boson triple production.

4 SM 2 Figure 13 :
Figure 13: The 68% C.L. (1σ, black solid) and 95% C.L (2σ, red dashed) limit on the pc t3 , d 4 q-plane for triple Higgs boson production with up to triple powers of anomalous operators at the matrix element-squared level, at 13 TeV/3000 fb ´1 (left), and 100 TeV/20 ab ´1 (right), marginalized over the c t2 and d 3 anomalous couplings.Note the differences in the axes ranges at each collider.
where c i P td 3 , d 4 , c g1 , c g2 , c t1 , c t2 , c t3 , c b1 , c b2 , c b3 u, are shown for E CM " 13.6 TeV in table 1 for triple production.For Higgs boson pair production at E CM " 13.6 TeV, see table9, in appendix F. The coefficients for both processes, at E CM " 13, 14, 27 and 100 TeV can be found in appendix F. The table should be read as follows: each value shown corresponds to the coefficient multiplying the product of the couplings of the corresponding first column and last row, respectively.As an example, with reference to table 1, the coefficient of the term proportional to d 4 d 3

Table 1 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson triple production, in the form σ{σ SM ´1 "ř i A i c i `ři,j B ij c i c j, where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 13.6 TeV.

Table 5 :
Figure 5: The 68% C.L. (1σ, black solid) and 95% C.L (2σ, red dashed) limit on the pc t3 , d 4 q-plane for triple Higgs boson production at 13 TeV/3000 fb ´1 (left), and 100 TeV/20 ab ´1 (right), marginalized over the c t2 and d 3 anomalous couplings.Note the differences in the axes ranges at each collider.The 3σ evidence and 5σ discovery limits on for triple Higgs boson production, for the c t3 and d 4 coefficients at 13 TeV/3000 fb ´1, and 100 TeV/20 ab ´1, marginalized over c t2 , d 3 and either d 4 , or c t3 .

Table 7 :
The 3σ evidence and 5σ discovery limits on for triple Higgs boson production with up to triple powers of anomalous operators at the matrix element-squared level, for the c t3 and d 4 coefficients at 13 TeV/3000 fb ´1, and 100 TeV/20 ab ´1, marginalized over c t2 , d 3 and either d 4 , or c t3 .

Table 8 :
The 68% C.L. (1σ) and 95% C.L (2σ) limits on c t3 and d 4 for triple Higgs boson production with up to triple powers of anomalous operators at the matrix element-squared level, at 13 TeV/3000 fb ´1, and 100 TeV/20 ab ´1, marginalized over c t2 , d 3 and either d 4 , or c t3 .lousoperatorsat the matrix element-squared level.Since we are only considering the set td 4 , c t3 , d 3 , c t2 u in our phenomenological analysis, there can only be diagrams with up to two insertions, and the third power only comes from interference terms of these diagrams with single-insertion operators.This can be generated within MG5 aMC, by: Figures12 and 13are to be compared to figs. 4 and 5, respectively, and tables 7 and 8 to tables 5 and 6, respectively.As can be clearly observed by examining the plots, as well as by the one-dimensional quantitative limits in the tables, this treatment does not substantially alter either the expected evidence/discovery boundaries, or the limits.

Table 9 :
where c i P td 3 , d 4 , c g1 , c g2 , c t1 , c t2 , c t3 , c b1 , c b2 , c b3 u, for Higgs boson pair and production triple production.The coefficients for Higgs boson pair production for E CM " 13.6, 13, 14, 27, and 100 TeV are shown in tables 9, 10, 12, 14, 16, respectively, and for triple production in tables 1, 11, 13, 15, 17.Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson pair production, in the form σ{σ SM ´1 "ř i A i c i `ři,j B ij c i c j, where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 13.6 TeV.

Table 10 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson pair production, in the form σ{σ SM ´1 " ř i A i c i `ři,j B ij c i c j , where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 13 TeV.

Table 11 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson triple production, in the form σ{σ SM ´1 " ř i A i c i `ři,j B ij c i c j , where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 13 TeV.

Table 12 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson pair production, in the form σ{σ SM ´1 "ř i A i c i `ři,j B ij c i c j, where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 14 TeV.

Table 13 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson triple production, in the form σ{σ SM ´1 " ř i A i c i `ři,j B ij c i c j , where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 14 TeV.

Table 14 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson pair production, in the form σ{σ SM ´1 " ř i A i c i `ři,j B ij c i c j , where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 27 TeV.

Table 16 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson pair production, in the form σ{σ SM ´1 " ř i A i c i `ři,j B ij c i c j , where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 100 TeV.

Table 17 :
Polynomial coefficients, A i (second column only) and B ij , relevant for the determination of the cross section for leading-order Higgs boson triple production, in the form σ{σ SM ´1 " ř i A i c i `ři,j B ij c i c j , where c i P td 3 , c g1 , c g2 , c t1 , c t2 , c b1 , c b2 , u, at E CM " 100 TeV.