Triple Higgs boson production to six b-jets at a 100 TeV proton collider

We investigate the production of three Higgs bosons at a proton-proton collider running at a centre-of-mass energy of 100 TeV, all of which decay into b-jets. This final state encapsulates by far the largest fraction of the total cross section of triple Higgs boson production, approximately 20%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20\%$$\end{document}. We examine, by constructing detailed phenomenological analyses, two scenarios: (i) one in which the triple and quartic Higgs boson self-couplings are modified independently by new phenomena with respect to their Standard Model (SM) values and (ii) an extension of the SM by a gauge-singlet scalar that could drive first-order electroweak phase transition, within the context of the so-called xSM. In the former, we find that competitive constraints of O(1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {O}}(1)$$\end{document} can be placed on the quartic coupling and in the latter we demonstrate that it will be possible to obtain important information on the structure of the extended scalar sector.


Introduction
In the past decade of operation of CERN's Large Hadron Collider (LHC), the landscape of particle physics has changed dramatically. The discovery of the Higgs boson and the lack of stark signals of new phenomena around the TeV scale are defining characteristics of this new era. In the years to come the Higgs boson is set to become itself a tool for exploration and discovery. This will be particularly true at the future circular collider (FCC), which is planned to be hosted in a 100 km tunnel, envisioning an ensemble of e + e − , e + p and pp collider programmes through towards the end of the 21st century [1][2][3][4][5]. Taken together, all of these programmes aim to a e-mail: apapaefs@cern.ch b e-mail: gtx@nikhef.nl c e-mail: m.zaro@nikhef.nl map the properties of the Higgs boson and the electroweak gauge bosons with an accuracy order(s) of magnitude better than today and to improve by almost an order of magnitude the discovery reach for new particles.
A particular "flagship" target of the FCC will be the investigation of the Higgs potential, through the measurement of the Higgs boson's (h) self-interactions that can be written, post-electroweak symmetry breaking (EWSB), as: where v 0 246 GeV is the Higgs vacuum expectation value (vev), m h 125 GeV is the Higgs boson mass and the selfcouplings take the values λ 3 = λ 4 = m 2 h /2v 2 0 ≡ λ SM within the SM. Legacy LHC measurements are expected to provide an O(1) measurement of the triple coupling, λ 3 , with respect to its SM value [6,7], and no significant direct information on the quartic self-coupling λ 4 . On the other hand, several studies have demonstrated the potential of the protonproton programme of the FCC (the FCC-hh), to constrain the triple coupling to within a few percent of the SM value, particularly through the production of Higgs boson pairs [1,[8][9][10][11][12][13][14][15]. Several studies have also hinted that constraints are possible on the quartic coupling at the FCC-hh, either indirectly in double Higgs boson production [16,17], or directly through triple Higgs boson production [18][19][20][21][22][23][24][25][26]. Up until now, in the case of the latter process, the following final states have been considered: The sum of all these channels represents less than 10% of the total branching ratio of hhh. In the present article, we investigate for the first time, to the best of our knowledge, the process that encapsulates by far the largest branching ratio: the case in which all three Higgs bosons decay into bottom quarks (bb), resulting in complex final states involving six b-jets. 1 In addition to understanding EWSB, non-standard Higgs boson self-couplings might provide the first experimental evidence of extra gauge-singlet scalars at the weak scale. These new scalar particles could "catalyse" electroweak phase transition, turning it into a violent, out-of-equilibrium event accompanied by massive entropy production (a firstorder transition), enabling electroweak baryogenesis and thus explaining the observed matter-antimatter asymmetry, see e.g. [36][37][38][39][40][41]. Evidence of such phenomena in multi-scalar production processes could materialise, for example, even in the case where the mixing of this new scalar and the "SMlike" Higgs boson is small. Indeed, current limits put an upper bound to the mixing angle that leads to cos θ 0.85 which at the end of the high-luminosity run of the LHC this is expected to be 0.95 [42]. First indications of the existence of these singlets could arise in resonant SM-like Higgs boson pair production for example, either at later stages of the LHC or during the FCC-hh lifetime. Such signals, along with the measurement of the SM-like Higgs self-coupling through non-resonant Higgs boson pair production, may not be sufficient to understand the nature of the additional singlet scalar. The production of three of these scalar particles, such as triple SM-like Higgs boson production, the main object of this article, could provide essential additional information both on the triple scalar couplings and on the quartic couplings. We demonstrate that this is possible by employing the six b-jet final state that maximises the cross section. 2 The article is organised as follows: in Sect. 2 we discuss the setup used and describe the phenomenological analysis in the context of triple SM Higgs boson production. In Sect. 3 we discuss the constraints that can be obtained in the anomalous coupling picture, where the self-couplings are rescaled with respect to the SM values, and in Sect. 4 we investigate in explicit benchmark scenarios, the potential for discovering triple Higgs boson production in the presence of a singlet scalar that can viably generate a first-order electroweak phase transition, taken from [44]. We conclude in Sect. 5. In "Appendix A" we provide investigations of relevant uncertainties entering our analysis. 1 We note that the equivalent final state in Higgs boson pair pair production, leading to 4 b-jets, has been considered extensively in both phenomenological and experimental studies, see e.g. [27][28][29][30][31][32][33][34][35]. 2 We would like to note here that the six b-jet final state might be also interesting in the context of hh + Z , see e.g. [43], or any triple neutral boson final state.

The setup
In what follows, we generate parton-level events either at leading order or next-to-leading order by using MadGraph5 _aMC@NLO [45,46] and shower/match them via the MC@NLO method [47] where appropriate, via the HERWIG (7.1.5) parton shower [48][49][50][51][52][53]. We include modeling of the hadronization and the underlying event but no detector effects beyond geometry. We use the parton density function set NNPDF23_lo _as_0130_qed [54] throughout the chain of event generation.
For the analysis, we cluster final-state particles with transverse momentum p T > 100 MeV into anti-k T jets [55] with radius parameter R = 0.4 via the FastJet package [56]. We use the HwSim package [57] for HERWIG to write out event files for each sample in a custom compressed ROOT format [58] and to perform the phenomenological analysis.

Differential distributions in pp → hhh at 100 TeV
In the present subsection we investigate the form of differential distributions for the hhh signal within the SM. Variations of the shapes of these distributions due to the effect of new phenomena are considered in the respective sections below: in Sect. 3 we show variations due to different values of the anomalous couplings and in Sect. 4 we show variations in the presence of a singlet scalar at various masses. We refer the reader to [19] for additional distributions, including a comparison to Higgs boson pair production at 100 TeV.
We show in Fig. 1, the invariant mass of (any) two or all three Higgs bosons reconstructed from Monte Carlo truth with no cuts applied, M hh and M hhh , respectively. The former Fig. 1 The invariant mass of any two (hh, dotted green) or all three (hhh, solid red) Higgs bosons in triple Higgs production at the FCC-hh at 100 TeV, reconstructed from Monte Carlo truth with no cuts applied Fig. 2 The transverse momentum of Higgs bosons ordered from hardest to softest (solid red, dotted green, dashed yellow, respectively), in triple Higgs production at the FCC-hh at 100 TeV, reconstructed from Monte Carlo truth with no cuts applied peaks at ∼ 300 GeV whereas the latter at ∼ 600 GeV. In Fig. 2 we show the Monte Carlo truth transverse momentum of the Higgs bosons ordered from hardest to softest. The transverse momentum distributions peak at ∼ 200 GeV, ∼ 150 GeV and ∼ 50 GeV from hardest to softest, respectively.

Event generation
The simulation of final states containing up to six coloured objects remains a challenge to this day, even at tree level. In the present study we provide initial estimates by considering the efficiency of a phenomenological analysis on (bb)(bb)(bb) final states. We stress here that we have simulated the QCD-induced (bb)(bb)(bb) exactly at tree level. 3 For backgrounds which arise from charm-jets or light jets being mis-identified as b-jets (i.e. the reducible backgrounds), we have estimated the cross sections and assumed the analysis efficiencies to be identical to the equivalent process with b-quarks, factoring out the mis-identification rates.
For the irreducible backgrounds, i.e. those that constitute the "real" (bb)(bb)(bb) final states, we have considered processes that contain three bosons (either a Higgs boson or a Z boson) that each then decay into (bb): hh Z, h Z Z, Z Z Z. We have included the loop-induced gluon-fusion component in the case of h Z Z and Z Z Z. 4 Furthermore, we have considered backgrounds with either one or two bosons plus (bb) that originate from QCD interactions: h Z + (bb), hh + (bb), 3 In particular, this was made possible thanks to the technique of Ref. [46] that performs a Monte Carlo over helicities. 4 All processes incorporate the full effect of spin correlations, implemented via MadSpin [59] that follows the method of [60], apart from gg → h Z Z and gg → Z Z Z. This is not expected to have a significant impact on the analysis efficiencies. Z Z +(bb) and Z +(bb)(bb), h +(bb)(bb). Of the aforementioned processes, Z + (bb)(bb) and h + (bb)(bb) turn out to be the largest contributors to total background cross section. For the latter process, h + (bb)(bb), we also consider the Higgs effective theory contributions (i.e. including the effective interaction ggh), which constitute approximately 3/4 of the cross section. 5 However, we have found that the largest background component by far is the pure QCD production of (bb)(bb)(bb). We did not generate large samples of events for the reducible backgrounds that would originate from mis-tagging light or charm jets to b-jets; instead, we considered estimates of their contributions to the h+jets and QCD multi-jet backgrounds, see Sect. 3.2. The details of the analysis are presented in the next subsection.
We have simulated the triple Higgs boson signal at (loopinduced) leading order and the quark-anti-quark-initiated component of the tri-boson processes at next-to-leading order. We have generated samples of 10 4 events for all signal processes, except for the SM hhh, for which we generated 10 5 events to obtain statistically reliable estimates of the significance.
All other processes have been simulated at leading order. To take into account the higher-order corrections, we multiply all leading-order cross sections by a K -factor of 2. The size of the higher-order corrections is well-motivated for the hhh signal by approximate calculations, see [62]. In Appendix A, we provide variations of the K -factor for the backgrounds to take into account this uncertainty, while given that a full NLO computation would be needed, we do not consider effects due to shapes. We have imposed generationlevel cuts on processes that involve quarks of QCD origin. We list these cuts in Table 1.
We emphasise that the simulation of processes with more than three final-state legs at next-to-leading order is an essential aspect that should be addressed in future studies at higherenergy hadron colliders, as such final states will become increasingly common.

Analysis details
We give here the details of the phenomenological hadronlevel analysis that are common between the different new physics scenarios that we consider. We ask for the events to contain exactly six identified bjets with transverse momentum p T > 45 GeV. We ask for these jets to lie within a pseudo-rapidity of |η| < 3.2 and we also ask for the distance between any two b-jets to satisfy ΔR > 0.3. The latter choice is simply to bring all processes on equal footing, given that the backgrounds that contain QCD-initiated b-quarks also obey a generation-level cut of ΔR > 0.2. We consider the potential impact of reducing the pseudo-rapidity coverage for the identified b-jets on our conclusions in "Appendix A". For each of the 15 possible arrangements I = {i j, kl, mn} of the six b-jets into pairs we construct the observable: where M qr is the invariant mass of the b-jet pairing qr in the arrangement of pairings I and m h is the Higgs boson mass. Given that it is challenging to determine experimentally the charge of the b-quarks that initiated the b-jets, we consider the minimisation of the χ 2 observable over all the possible pairings. The arrangement of pairings I that gives the minimum of χ 2 , which we call χ 2 min , defines the three "reconstructed Higgs bosons", h i r , for i = {1, 2, 3}. For this specific combination we calculate the absolute difference with the Higgs mass and order from smallest to larger: (Δm min , Δm mid , Δm max ). We impose cuts on the observables χ 2 min , Δm min , Δm mid and Δm max . Furthermore, we impose cuts on the transverse momentum of the hardest, second hardest and softest reconstructed Higgs boson, p T (h i r ) for i = {1, 2, 3}. We also impose cuts on the distances between the reconstructed Higgs bosons, ΔR(h i r , h j r ). Finally, we ask for the distances between the two b-jets that comprise the reconstructed Higgs bosons, ΔR bb (h i ), to satisfy certain upper bounds. The values of the cuts on these observables are summarised in Table 2. 6

Anomalous self-couplings
We first consider a scenario in which the triple and quartic couplings are modified independently of each other. This "agnostic" anomalous coupling approach does not necessarily represent a physically viable theory, but allows for an investigation of the possible constraints that can be obtained for SM-like triple Higgs boson production. We thus consider interactions of the form: where the coefficients c 3 and d 4 represent the modifications of the triple and quartic Higgs boson self-interactions respectively. Assuming that the Yukawa couplings to the top and bottom quarks remain unchanged, these interactions will induce changes to the main production channel for triple Higgs boson production, that proceeds through gluon fusion, mediated by heavy quark loops. Example Feynman diagrams are shown in Fig. 3, together with their scaling with the coefficients c 3 and d 4 .
In Fig. 4 we show a variation of the cross section at a 100 TeV proton collider, normalised to the SM value. Evidently, variations of the triple self-coupling via c 3 produce larger changes than equivalent variations with d 4 . A fit of the cross section on this plane yields a polynomial in c 3 and d 4 which is quartic in c 3 and quadratic in d 4 . This is because there exist diagrams with two insertions of the triple self-coupling c 3 in triple Higgs boson production (Fig. 3d), whereas there are only diagrams with at most a single insertion of d 4 (Fig. 3c) at this order. The dependence of the cross section on c 3 and d 4 , normalised to the SM cross section, was fitted as: Example Feynman diagrams contributing to Higgs boson triple production via gluon fusion in the Standard Model, taken from [19]. The vertices highlighted with blobs indicate either triple (blue) or quartic (red) self-coupling contributions. In the (c 3 , d 4 ) model, a produces matrix elements unmodified by the anomalous couplings (at this order), b produce terms The interference of such diagrams produces the terms that appear in Eq. 4

Fig. 4
The cross section for triple Higgs production with modified self-couplings (λ 4 = λ SM (1 + d 4 ) and λ 3 = λ SM (1 + c 3 )) at the FCChh at 100 TeV, normalised to the SM value. The black star indicates the SM point tion. For example, in the context of the SM effective field theory (EFT) at D = 6, setting the relation d 4 = 6c 3 , one reproduces to a good approximation, the fit of [19]. 7 We show in Fig. 5 the normalised invariant mass distribution for the triple Higgs boson system for a few extreme values of the modifications of the quartic or triple coupling. It is clear that the anomalous couplings can modify substantially not only the cross section but also the distributions. This implies that for the same analysis cuts, the efficiency values will vary across the (c 3 , d 4 )-plane. The efficiency ranges from ∼ 0.5% to ∼ 2%. A fit of the analysis efficiency (≥ 0), for the cuts given in Table 2, over the (c 3 , d 4 )-plane yields: where some coefficients are zero up to the uncertainty obtained by the Monte Carlo sample employed. The cross section after cuts over the (c 3 , d 4 )-plane can be obtained by convolving Eqs. 4 and 5.

Background processes
Some of the background processes will be affected at the order we are considering by the rescaling of the selfcouplings of Eq. 3, an effect that should be taken into account in a future analysis. However, we found that the processes that are affected at leading order by the anomalous couplings, i.e. those of the form hh + X , where X = Z or bb, constitute sub-permille contributions to the sum of all backgrounds after our analysis cuts are applied (see results of Table 3). Therefore we do not consider these variations in our anal- Table 3 The processes considered in the six b-jet analysis, for the Standard Model. The second column shows the generation-level cross sections with the cuts (if any) as given in the main text. The Z bosons were decayed at generation level and hence the cross section is given with the Z branching ratios applied. The third column shows the starting cross section for the analysis, including the branching ratio to (bb)(bb)(bb), with a flat K -factor of K = 2.0 applied to all tree-level processes as an estimate of the expected increase in cross section from leading order to next-to-leading order. The fourth column gives the analysis efficiency and the final column gives the expected number of events at 20 ab −1 of integrated luminosity at 100 TeV. The results are given for perfect b-jet tagging efficiency. The label "ggF" implies that it is gluon-fusion initiated ysis, instead only considering their SM counterparts as an order-of-magnitude estimate.
It is also evident that in Table 3 we have only included irreducible processes, those that are identical at parton level in flavour content to the signal: (bb)(bb)(bb). As discussed previously, the degree of the contamination from reducible backgrounds, those that come from the mis-identification of light jets or charm-jets to b-jets, can be estimated by assuming that the efficiency of the analysis is identical to that of the equivalent irreducible ones. Explicitly, we will assume e.g. that the probability of a (bb)(bb)(cc) event passing the analysis cuts is identical to (bb)(bb)(bb), multiplied by the probability that two charm jets are mis-identified as b-jets. We will assume that the probability of a charm-jet being mis-identified as b-jet is P c→b = 0.1 and that of light jets is P j→b = 0.01, and that these values are independent of the b-tagging efficiency which we will take to range from perfect (100%) to the "worst-case scenario" of 80%, see "Appendix A". 8 Table 4 shows the starting cross sections of the main reducible processes and the estimated contribution to the total cross section of the equivalent irreducible process, QCD six b-jet production by taking into account appropriate rescaling with powers of P c→b and P j→b . Given our results, the reducible six-jet QCD backgrounds are expected to contribute O(10%) to O(30%) of the total tagged six b-jet back- Table 4 The reducible background processes considered in the six b-jet analysis. The second column shows the generation-level cross sections with the cuts identical to the ones applied to the irreducible processes ( Table 2). The third column shows the cross section after the mis-tagging rates have been applied. We only consider processes equivalent to QCD 6 b-jet production. We do not consider process that contain mis-tagged light and charm jets at the same time ground, for perfect b-tagging to P b→b = 0.8, respectively. Therefore it is clear that the contributions are sub-dominant with respect to the irreducible process and from here on we absorb them in the overall uncertainty of the cross section estimates, the effect of which is also examined in "Appendix A".

Results for anomalous triple Higgs boson production
As a result of the analysis described in Sect. 2.4, we show in Fig. 6 the expected significance that would be obtained on the (d 4 , c 3 )-plane for an integrated luminosity of 20 ab −1 and assuming perfect b-tagging. This result demonstrates that the Given that the constraints on the triple self-coupling at the FCC-hh will reach the percent level, we also consider a scenario in which c 3 = 0, allowing for variations of the quartic self-coupling through d 4 . The resulting one-dimensional significance is shown in Fig. 7 for the case of perfect b-tagging. The constraint in this scenario would then be, at 95% confidence level (i.e. 2σ ), d 4 ∈ [−1.7, 13.3] as indicated by the red dashed lines in the figure. We defer the equivalent plots with reduced b-tagging efficiencies and the range of the pseudorapidity of b-tagging to "Appendix A". We note that the significance for the SM triple Higgs boson production (d 4 = 0, c 3 = 0) is ∼ 1.7σ , up to the Monte Carlo uncertainties.

Triple Higgs boson production in the presence of a singlet scalar
The discussion of the so-called xSM and the study of this section follows from [44]. A more detailed discussion of the model and its relation to strong first-order electroweak phase transition is discussed therein. See also [66] for a similar study. The most general form of the xSM that depends on the Higgs doublet, H , and a gauge-singlet scalar, S, is given by (see, e.g. [38,40,44,67,68]): where the interactions proportional to a 1,2 constitute the Higgs "portal" that links the SM with the singlet scalar. We follow the study of [44] in retaining all of the parameters, i.e. we do not impose a Z 2 symmetry that would preclude terms of odd powers of S. After EWSB, the Higgs doublet and the singlet scalar both attain vevs: with v 0 246 GeV and S → x 0 + s. Inevitably, the two states h and s mix through both the Higgs portal parameters a 1 and a 2 as well as the singlet vev. Diagonalising the mass matrix, one obtains two eigenstates, denoted by h 1 and h 2 , where: where θ is a mixing angle that can be expressed in terms of the parameters of the model. For θ ∼ 0, h 1 ∼ h and h 2 ∼ s. We will identify the eigenstate h 1 with the state observed at the LHC, and hence set m 1 = 125 GeV.
All the couplings of h 1,2 to the rest of the SM states are simply obtained by rescaling by: with X X any SM final state. This allows for constraints to be imposed on cos θ through the measurements of Higgs signal strengths. We concentrate on the scenario m 2 ≥ 2m 1 , allowing for resonant h 2 → h 1 h 1 , with no new decay modes appearing for the h 1 . 9 The triple couplings between the scalars h 1 and h 2 , representing terms of the form , are given by: where we have defined c θ ≡ cos θ and s θ ≡ sin θ . The quartic couplings, representing terms of the form V (h 1 , h 2 ) ⊃ λ i jkl h i h j h k h l , i, j, k, l = {1, 2}, are given by: The above couplings will lead to processes with multiple h 1 and h 2 in the final state.
In [44], the authors studied parameter-space points, satisfying conditions on the scalar sector of the xSM that lead to strong first-order electroweak phase transition (SFOEWPT). They then derived benchmark points taken from this allowed set that leads to enhanced resonant Higgs boson pair production, i.e. h 2 → h 1 h 1 , considering the phenomenological consequences, i.e. whether enhanced h 1 h 1 would be observed at future colliders, including a 100 TeV proton collider. One of the main conclusions was that such a collider could probe nearly all of the viable SFOWEPT-viable parameter space through this process, leading to a potential discovery of the xSM.
Here we consider the benchmark points of [44] in the context of (SM-like) triple Higgs boson production, pp → h 1 h 1 h 1 , which can potentially lead to a measurement of both the triple and quartic couplings in the xSM, in the event of discovery. Furthermore, there could be fine-tuned points in the xSM that lead to some of the scalar couplings being small. In that scenario, triple Higgs boson production could conceivably provide an alternative route for discovery of the xSM. We show in Tables 5 and 6 in the next section the parameters for the benchmark points, which are labelled in [44] as "B1max" to "B11max" and "B1min" to "B11min".

Triple Higgs boson production in the xSM
The process by which three h 1 scalars are produced via gluon fusion consists of diagrams identical to those that appear in Fig. 3, with the addition that there exist diagrams with SMlike Higgs propagators (i.e. h 1 in this case) substituted by h 2 . The strength of the interactions that appear in these diagrams is governed by the triple and quartic couplings of Eqs. 9 and 10. Note that the triple h 1 coupling, λ 111 , will also be modified in the xSM. In general there will be an intricate interference pattern between all the contributing non-resonant and resonant diagrams. Our aim here is not to provide a detailed study of these effects; instead we investigate the observability of triple Higgs boson production, pp → h 1 h 1 h 1 , in the context of the six b-jet final state, focussing on the SFOEWPT benchmark points provided in [44], which appear in Tables 5  and 6. For each point we also give the total triple h 1 production cross section as a ratio to the SM hhh, including the full (top or bottom quark) loop structure and interference effects. For comparison we have also calculated the total h 1 pair production cross section as a ratio to the SM hh. One can observe that the enhancement in h 1 h 1 h 1 production can be larger than the enhancement in h 1 h 1 .
We show in Figs. 8 and 9, the invariant mass of the three Higgs boson system and the transverse momentum of the hardest Higgs boson in triple h 1 production within the xSM for three benchmark points as well as the SM expectation for comparison. The double-peak structure that is present in the distributions is physical and is due to the possibility of either an on-shell decay h 2 → h 1 h 1 h 1 , leading to a peak in M hhh at ∼ m 2 , or an on-shell decay h 2 → h 1 h 1 with either h 2 or h 1 being off-shell in a preceding s-channel propagator, leading to the peak in M hhh at ∼ m 2 + m 1 . We note that a similar effect was pointed out in [70] in pp → hS → hγ γ , in the context of a Z 2 -symmetric singlet scalar model. Table 5 Values of the various xSM independent and dependent parameters for each of the benchmark values chosen to maximize the σ · B R(h 2 → h 1 h 1 ) at a 100 TeV proton collider, taken from [44]. The ratio of cross sections of h 1 h 1 to SM hh and of h 1 h 1 h 1 to hhh production are given in the last two columns   Table 7 shows the significance of the analysis applied to the 22 benchmark points B1max-B11max and B1min-B11min. The analysis has not been optimised for the specific features of these points, but the cuts are instead applied verbatim following those described previously (Table 2). It is quite likely, as was shown in [44], that the h 2 mass will be known during the lifetime of the FCC-hh through the observation of resonant production pp → h 2 → h 1 h 1 . 10 This information could then be employed in the analysis to enhance the significance of the h 1 h 1 h 1 further, particularly taking into account the "double-peak" structure that we have pointed out in Sect. 4.2. Furthermore, cuts affected by the changes in 10 See also [71] for an investigation of resonant Higgs boson pair production in the SM with an additional gauge-singlet scalar.
the transverse momentum distributions as well as the angular distances can be subject to further optimisation.
Given the values of the significance that we find here, it is conceivable that the pp → h 1 h 1 h 1 channel will play a crucial role in understanding the extended scalar sector in many viable scenarios of scalar gauge-singlet models that satisfy the constraints provided by requiring a SFOEWPT. 11

Conclusions
We have investigated triple Higgs boson production at a future proton collider with centre-of-mass energy 100 TeV, Fig. 8 The (normalised) invariant mass distribution of the three Higgs boson system in triple h 1 production within the xSM at the FCC-hh at 100 TeV, reconstructed from Monte Carlo truth with no cuts applied. We show three benchmark points from [44], as well as the SM expectation for comparison Fig. 9 The (normalised) transverse momentum distribution of the hardest Higgs boson in triple h 1 production within the xSM at the FCC-hh at 100 TeV, reconstructed from Monte Carlo truth with no cuts applied. We show three benchmark points from [44], as well as the SM expectation for comparison in the case when all three Higgs bosons decay to bottomanti-bottom quark pairs, producing six b-jets. We have constructed a detailed phenomenological hadron-level analysis including the effects of detector geometry. This analysis was applied to two scenarios: in the first one, SM-like triple Higgs boson production, we allowed for "anomalous" modifications of the triple and quartic self-couplings independently. For the SM point, (d 4 , c 3 ) = (0, 0), we demonstrated that significances of ≈ 2σ can be obtained from the six b-jet final state alone. Furthermore, we have shown that a constraint of d 4 −2 could be obtained in the case that the triple coupling is measured to be close to the SM value, c 3 ∼ 0. These results are competitive with previously studied final states, rendering the six b-jet process an important contribution to the study of the self-couplings of the SM Higgs boson. In the second scenario, we considered an extension of the SM by a gauge-singlet scalar that could drive strong first-order electroweak phase transition. We investigated the triple production of the resulting SM-like scalar in the particular six b-jet final state, for several well-motivated benchmark points compatible with strong first-order electroweak phase transition, and we concluded that large significances can be obtained for many of these. This motivates further study of the triple Higgs boson process in the context of future collider studies of scalar singlet models. Finally, we emphasise that our conclusions are affected by uncertainties due to the absence of higher-order calculations for several of the background processes and details of the performance of the detector, particularly with respect to the tagging efficiencies, acceptance rates, resolution and triggers. Once these uncertainties have been better understood, a more detailed analysis, e.g. considering the differences between the radiation pattern of the colour singlet Higgs boson and QCD, or employing more advanced multivariate techniques, could lead to higher significances. Nevertheless, we have demonstrated here by varying several parameters, that the six b-jet process will certainly constitute an important component of the study of triple Higgs boson production at a future 100 TeV hadron collider. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and therefore no relevant data has been produced.] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Appendix A: Variations and uncertainties
In Figs. 10 and 11 we show variations of the significance on the (d 4 , c 3 )-plane with 20 ab −1 , when the b-tagging efficiency is reduced from 100% (perfect), to 90% and 80%, respectively. We also show the c 3 = 0 significance over the values of d 4 in Figs. 12 and 13. The equivalent constraints at 95% C.L. on d 4 would then be, respectively, The FCC-hh detector coverage over which b-jets will be tagged might also be tighter. Maintaining perfect b-tagging within the restricted region, with the given set of cuts, we apply |η b | < 2.5 instead of |η b | < 3.2. The 95% C.L. constraint on d 4 in this case would be slightly shifted, yielding a range: d 4 ∈ [−3.1, 14.2] at 95% C.L.. The significance for the SM point (d 4 , c 3 ) = 0 is also reduced to 1.3σ .
Finally, by far the largest theoretical uncertainty is in the K -factors for the tree-level background processes. In the main part of this article, we have applied K = 2 to all of these. If this is increased to K = 3 for all tree-level back-    (1 + d 4 )) and no modification to the triple self-coupling (c 3 = 0) at the FCC-hh at 100 TeV. We have assumed an integrated luminosity of 20 ab −1 and b-tagging of 90%. The red dashed lines indicate the 2σ points for the constraints on d 4 ground processes, we would obtain ∼ 4.3 × 10 4 background events, yielding a significance for the SM point of ∼ 1.3σ . On the other hand, if this is reduced K = 1.5, the number of background events would decrease to ∼ 2.1 × 10 4 with a significance of ∼ 1.9σ . We reckon that shape uncertainties may also be important, in particular those effects due to extra radiation generated before the production of b quarks. However given the complexity related to assessing this kind Fig. 13 The significance for our analysis of triple Higgs boson production in the (bb)(bb)(bb) final state with modified quartic self-coupling (λ 4 = λ SM (1 + d 4 )) and no modification to the triple self-coupling (c 3 = 0) at the FCC-hh at 100 TeV. We have assumed an integrated luminosity of 20 ab −1 and b-tagging of 80%. The red dashed lines indicate the 2σ points for the constraints on d 4 of uncertainties, we do not include them in this explorative study.