Triple Higgs Boson Production at the Large Hadron Collider with Two Real Singlet Scalars

We investigate the production of three Higgs bosons in the Two Real Singlet extension of the Standard Model, where the scalar sector is augmented by two additional real scalar fields which are singlets under the Standard Model gauge group. The model contains three neutral CP-even scalars, allowing for resonant production and asymmetric decay chains. We focus on the signature $p p\,\rightarrow\,h_3\,\rightarrow\,h_1\,h_2\,\rightarrow\,h_1\,h_1\,h_1$, where we identify $h_3$ as the heaviest scalar state, $h_2$ as the second heaviest and the lightest, $h_1$, as the Standard Model-like Higgs boson discovered by the Large Hadron Collider experiments. The dominant final state occurs when all three Higgs bosons decay to bottom-anti-bottom quark pairs, $h_1\,\rightarrow\,b\,\bar{b}$, leading to 6 $b$-jets. Taking into account all current theoretical and experimental constraints, we determine the discovery prospects for this channel in future runs of the Large Hadron Collider, as well as in the high-luminosity phase.


Introduction
With the discovery of a particle which complies with the expected properties of the Higgs boson of the Standard Model (SM) by the CERN Large Hadron Collider (LHC) experiments in 2012 [1,2], particle physics has entered an exciting new era. Although current experimental results agree rather well with the predictions of the SM, both experimental and theoretical uncertainties allow for new phenomena that may be observable either at current or future colliders. In the present article we focus on the particularly interesting possibility of models that extend the scalar sector of the SM by additional scalar fields that transform as singlets under the SM gauge group. Such models may provide solutions to a multitude of fundamental open questions: they could contain viable candidates of dark matter or enable mechanisms that could explain the observed cosmic matter-anti-matter asymmetry (see, e.g. [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]). They are also very rich in terms of their collider phenomenology, introducing new physical scalar states that can participate in cascade decays.
In this work, we concentrate on the triple production of scalar final states resulting from the asymmetric decay chain: where h 1,2,3 are the physical scalar states of a model with an extended scalar sector. We require that one of these scalars, specifically the h 1 boson, is identified with the 125 GeV SM Higgs particle, including agreement with all current measurements. The other scalars, however, can lie in any mass range, as long as all theoretical and experimental constraints are satisfied. As we are interested in the discovery potential of colliders that probe the TeV scale, we choose to consider scenarios with masses 1 TeV. In order to allow for the decay chain (1.1), and assuming CP conservation, the new physics model under consideration needs to contain at least three CP-even scalar states. One of the simplest ways to realise this is through models that extend the SM scalar sector by two additional singlet fields. The two real 1 singlet extension that contains three unstable physical scalars has been widely investigated in the literature, see, e.g. [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37] for recent discussions.
The LHC experimental collaborations have already largely scrutinised models which allow for several scalar particles in the final state, including searches for processes with symmetric di-scalar production via resonances, p p → h 2 → h 1 h 1 , where either h 1 or h 2 take the role of the 125 GeV SM-like scalar . Furthermore, in [47] the ATLAS collaboration also interpreted their results for the above production and decay chain for pure beyond-the-SM (BSM) scalars, i.e. neither h 1 nor h 2 assume the role of the SM Higgs boson. For models with extended scalar sectors, however, triple couplings between different mass states, λ h i h j h k , can best be probed at leading order in resonant production modes such as the decay chain (1.1). Such states have e.g. been discussed in [29,35,37,[62][63][64][65], but currently no experimental results for such searches are available.
While the investigation of the process pp → h 1 h 2 with decays into SM-like final states is an important quest as such, 2 here we plan to focus on the specific case where h 2 → h 1 h 1 , 1 Models with two real singlets or one complex singlet field are equivalent, given that potential additional symmetries are correctly translated, see, e.g. [23,24]. 2 For representative benchmark points for such scenarios, see e.g. [24].
leading to triple scalar final states as indicated above. In the SM, the production cross section for the triple Higgs boson final state is the lowest-order process to include the quartic Higgs self-coupling. At the LHC's nominal centre-of-mass energy, 14 TeV, the corresponding cross section in the SM is diminishingly small, ∼ 0.1 fb [66,67], rising up to a cross section of ∼ 5.6 fb at a 100 TeV proton-proton collider [68]. While the quartic self-coupling in the SM can also be indirectly constrained [69][70][71][72], direct determination seems to call for future high-energy proton-proton colliders [73][74][75][76][77][78] or a possible muon collider [79,80]. As discussed above, the simplest realisation that achieves (1.1) are models that extend the SM by two additional real scalar fields, which are singlets under the SM gauge group. We consider here a specific version, the "Two Real Singlet Model" (TRSM) [24], where in addition two Z 2 symmetries are imposed, leading to a reduction of the available number of degrees of freedom. In the TRSM, the gluon-fusion pp → h 1 h 1 h 1 cross section is enhanced via the resonant production of h 3 and can reach up to 140 fb at the LHC. 3 While direct searches for an SM-only triple Higgs boson production are not very promising at current centre-of-mass energies, we will show that several benchmark points of the TRSM are within a 2 − 4 σ significance range with an integrated luminosity of 300 fb −1 , reaching up to ∼ 5 σ for selected points, and can reach up to ∼ 16 σ for the full high-luminosity LHC (HL-LHC) nominal dataset of 3000 fb −1 .
This article is organised as follows: in section 2, we briefly review the model under consideration as well as the specific benchmark plane that our study focusses on. In section 3, we discuss current theoretical and experimental constraints. The event generation, cross sections and selection analysis are discussed in section 4. We present the results of our analysis in section 5. There we also present projections for the sensitivity of the full HL-LHC run for searches of heavy scalars within the TRSM into di-boson final states. Our summary and conclusions can be found in section 6.

Extending the Standard Model by Real Singlet Scalar Fields
The scalar potential of the SM can be extended by an additional sector of scalar fields that transform as singlets under the SM gauge group, leading to with the most general renormalizable expression for V singlets (Φ, φ i ) given by In this work we focus on the TRSM [24], which introduces two extra real scalar fields S and X. The number of free parameters is constrained by imposing the following discrete Z 2 symmetries: and where all SM particles transform evenly under both symmetries. The application of the discrete symmetries of eq. (2.3) reduces the scalar potential for two real singlet fields to: which is characterised by nine real couplings µ Φ , λ Φ , µ S , λ S , µ X , λ X , λ ΦS , λ ΦX , λ XS . All fields are assumed to acquire a vacuum expectation value (vev). The physical gaugeeigenstates φ h,S,X then follow from expanding around these according to: In this study we consider the broken phase in which v S , v X = 0 and v = v SM 246 GeV. Then, the discrete symmetries Z S 2 and Z X 2 are spontaneously broken and the scalars φ h , φ S , φ X mix into the physical states h 1 , h 2 and h 3 according to with the rotation matrix R given by (2.7) To simplify our discussion we have used the following notation when writing R in eq. (2.7): (2.9) Using the same notation as in [24], the entries of the first row in the matrix R are denoted as κ i ≡ R i1 for i = 1, 2, 3.
In principle, any of the three scalars can take the role of the SM-like Higgs boson resonance discovered by the LHC experiments, as long as the other parameters are set such that all experimental constraints are fulfilled. Here, however, we will focus on the scenario where the state h 1 is identified with the SM-like Higgs boson, and h 2 and h 3 are two new physical heavier scalars obeying the mass hierarchy (2.10) As previously described, there are 9 real parameters characterising the TRSM. However, the identification of h 1 as the SM Higgs boson fixes This leaves us with 7 independent parameters, which we chose as In this model all couplings for the mass eigenstates h i to SM particles are inherited from the SM-like Higgs doublet through the rotation from the gauge to the mass eigenstates, such that g i ≡ κ i g SM . For example, in a factorised approach, this leads to predictions for production cross sections of the form where σ SM (M i ) denotes the production cross section of an SM-like Higgs boson of mass M i . Furthermore, the total width of the h i scalars (i = 1, 2, 3) is given by: where Γ SM (M i ) corresponds to the width of a scalar boson of mass M i possessing the same decay modes as a SM Higgs of mass M i . The branching ratios corresponding to h i → xx, for x = h j (j = i) are then given by: where Γ SM xx (M i ) corresponds to the SM-like partial decay width of a scalar boson of mass M i for the final state xx. The scalar-to-scalar branching ratios are equivalently obtained via

Benchmark Scenario
As discussed in [24], depending on the values that the free parameters of eq. (2.12) assume, different realisations of the TRSM are possible, yielding a rich phenomenology at colliders. Here we concentrate on the "Benchmark Plane 3" (BP3) addressed in [24], which was carefully tailored to allow for a large region in the (M 2 , M 3 ) plane which obeys all current theoretical and experimental constraints, while at the same time allowing for a large h 1 h 1 h 1 decay rate. 4 BP3 is characterised by the numerical values of the parameters shown in table 1. 5

Constraints and Allowed Regions
Constraints on the TRSM have been discussed in detail in [24], and we essentially follow that description in this work. In particular, we include constraints from perturbative unitarity, the requirement that the potential is bounded from below and agreement with electroweak precision observables. Results from null searches at colliders for the additional resonances as well as agreement with the current signal strength measurements, have been tested using the HiggsBounds [82][83][84][85][86][87] and HiggsSignals [88][89][90][91] packages. We additionally made use of the ScannerS [26,29,92] code to cross-check several of the constraints discussed in this section.
In the rest of this section, we describe the constraints in further detail. 4 Note that, in addition, this rate depends on the mixing angles and additional vevs, which are fixed in BP3. 5 Note that we actually set M1 = 125 GeV in the analysis performed throughout this work.

Parameter
Value

Theory Constraints
We can derive constraints on the values that the masses M 2 and M 3 can assume by considering the perturbative unitarity of the 2 → 2 scalar scattering matrix in the TRSM. Moreover, we impose an upper limit |M i | ≤ 8 π on the eigenvalues M i of the scattering matrix M. These limits can be written in terms of the coupling constants as 6 where a 1,2,3 correspond to the roots of the following polynomial: The potential of eq. (2.4) additionally needs to be bounded from below. This requirement was implemented in the scan discussed in [24] using the conditions derived in [94,95], which we list here for completeness

Electroweak Precision Constraints
In the benchmark plane discussed here, constraints from electroweak precision observables have been imposed using the ScannerS interface, which calculates the oblique parameters S, T, U [96][97][98][99] from expressions in [100,101] and compares them to the most recent fit results of the GFitter collaboration [102], including all correlations.

Collider Constraints
To apply current constraints we employ the HiggsBounds (v5.9.0) and HiggsSignals (v2.5.1) packages. HiggsBounds takes a selection of Higgs sector predictions for any model as input and then uses the experimental topological cross-section limits from Higgs boson searches at LEP, the Tevatron and the LHC to determine if this parameter point has been excluded at 95% C.L.. HiggsSignals performs a statistical test of the Higgs sector predictions of arbitrary models with the measurements of Higgs boson signal rates and masses from the Tevatron and the LHC. HiggsBounds returns a boolean corresponding to whether the Higgs sector passes the constraints at 95% C.L. (true) or not (false). HiggsSignals returns a probability value (p-value) corresponding to the goodness-of-fit of the Higgs sector over several SM-like "peak" observables. The code contains searches up to the full LHC Run II luminosity, and we refer the reader to the documentation of the code for details [103].
For BP3, we found that searches for h 2,3 → V V from 2016 LHC Run II data [104][105][106] constrain some parts of the parameter space, in agreement with the results presented in [24].

Monte Carlo Event Generation
All the parton-level events used in the phenomenological analysis of the present article have been generated via the Monte Carlo (MC) event generator MadGraph5 aMC@NLO (v2.7.3) [107,108]. The TRSM signal MC samples were produced via a custom modification of the loop sm model to incorporate the additional scalar particles and their interactions with the SM particles. This yields a leading-order description of the signal, including the full top and bottom quark mass dependence and all interference effects between the contributing Feynman diagrams. The production of the samples for the background process, i.e. the final state that originates from the QCD production of (bb)(bb)(bb) constitutes the most challenging aspect of the event generation. Note that within the SM this entails the evaluation of 6762 Feynman diagrams. To address this challenge we heavily parallelised the event generation via the "gridpack" option provided by MadGraph5 aMC@NLO.
QCD parton showering, hadronization and underlying event simulation were all performed within the general-purpose MC event generator HERWIG (v7.2.1) [109][110][111][112][113][114][115]. Events were subsequently analysed via the HwSim module [116] for HERWIG which saves events in a ROOT compressed file format [117], with jets clustered using FastJet (v3.3.2) [118]. The anti-k T algorithm [119] with a radius parameter R = 0.4 was used to cluster jets. A detailed study of pile-up effects arising from secondary proton-proton interactions is beyond the scope of the present phenomenological study and will need to be addressed in a full experimental study that will include in conjunction a detailed description of detector effects. For a recent a discussion on the issue of pile-up mitigation and corrections, we would like to point out the reader to the detailed studies of Ref. [120], which demonstrate the degree of the effects on jet resolution and suggest approaches in the form of advanced techniques to improve on this.
To capture the detector effects, we only consider particles with transverse momentum p T > 100 MeV as being detectable. We do not consider any smearing of momenta coming from detector mis-measurements. Similarly, we do not take possible mis-identification of light or charm jets as b-jets into account. These assumptions are not expected to have a dramatic impact on the conclusions of the present study and we anticipate that a full experimental analysis will assess their effects in detail. Throughout this work, we assume a b-jet tagging efficiency of 0.7, which lies on the conservative side of 13 TeV ATLAS and CMS performance [121][122][123] and was also adopted in the studies presented in [124]. We have elected to consider a constant b-tagging efficiency with transverse momentum and pseudo-rapidity of the jets. This is justified since, for example, by examining Fig. 6 of Ref. [123], where the b-tagging efficiency appears to be relatively flat for both observables, and in particular above p T ∼ 30 GeV, with O(10%) uncertainty, which is precisely where we impose a cut on the b-jets in our analysis.

Cross Sections
We present the production cross sections for the pp → h 1 h 1 h 1 final state over BP3 in fig. 1, where in addition bounds from perturbative unitarity and the requirement for the potential to be bounded from below are shown. The cross sections displayed in this plot have been obtained following the leading-order MC description of section 4.1, which includes all gluon-fusion-initiated contributions as well as interference effects (e.g. the box diagrams . Note that, for points where indeed the h 3 resonant production contributes dominantly, one could additionally apply a K-factor to account for missing higher-order contributions, e.g. with respect to the NNLO+NNLL corrected predictions for production cross sections of an SM-like scalar with mass M 3 [81]. For our selected benchmark points within BP3, specified below, we found that these K-factors for gluon-gluon induced h 3 production are ∼ 2.5. 7 Furthermore, for all of our benchmark points we found that ∼ 93 − 99% of the cross section stems from the decay chain specified in eq. (1.1). For our analysis, we have selected specific benchmark points within BP3. The corresponding cross-section predictions for pp → h 1 h 1 h 1 as well as 6 b-quark final states are given 7 For parameter points where h3 production dominates, the total cross section is in addition sensitive to the total width of h3 and follows the scaling predicted by the narrow width approximation, i.e. σ pp→ h 1 h 1 h 1 ∼ Γ −1 3 . Therefore, percent-level differences in the width can induce similar changes in the final result.   in table 2. 8 Here we have taken the branching ratio of the h 1 to bb to be BR h 1 → bb = 0.5824 [81]. The SM background amounts to a cross section of 6.38 pb for the 6 b-quark final state from QCD-induced diagrams, including a K-factor of 2, typical for gluon-fusion processes. Additional backgrounds from electroweak processes, e.g. Z bb bb production with Z → bb, as discussed in [78], were found to be at least two orders of magnitude lower and have not been considered in our study. We expect that these will form a sub-dominant contribution with respect to the QCD background after the analysis cuts are imposed.

Selection Analysis
Our analysis has been adapted from that of ref. [78]. An event is analysed if it contains at least 6 b-tagged jets 9 with a transverse momentum of at least p T min,b = 25 GeV and a pseudo-rapidity no greater than |η b,max | = 2.5. These initial cuts are further optimised for 8 The widths for the three scalars have been calculated according to eq. (2.14), with SM-like widths from [125]. We list the corresponding values in Appendix B, together with the corresponding new physics branching ratios. 9 Since the Higgs bosons are produced with transverse momenta up to O(100) GeV, i.e. comparable to their mass, we do not expect the b-jets to frequently merge into a singlet jet and therefore we focus only on the "resolved" 6 b-jet scenario.
each of our signal samples, which are characterised by different combinations of M 2 and M 3 . We then select the 6 b-tagged jets with the highest transverse momentum and form pairs in different combinations, with the aim of first reconstructing individual SM-like Higgs bosons, h 1 , and subsequently the two scalars h 2 and h 3 . To this end, we introduce two observables:  The above procedure still allows for different approaches in the combination strategy, on which we briefly comment in Appendix C. We then "identify" candidates for the scalars h 2 and h 3 with the pairing configurations I min and J min which minimise χ 2,(4) and χ 2,(6) respectively, as described above. Note that this procedure does not guarantee that I min indeed reconstructs to h 2 ; in fact, we found this to be the case in about 40% on average for all benchmark samples, being slightly higher than a "blind guess" that would lead to a probability of 1/3. Based on the invariant mass of the b-jet combinations entering in I min and J min , we define two additional observables m inv 4b and m inv 6b . We wish to stress that we do not make explicit use of the values of M 2 and M 3 for the individual samples. The fact that the masses are different is however taken into account implicitly considering that we find different selection cuts depending on the concrete signal sample during the analysis. Our approach is already able to deliver a good selection performance and using additional information on the assumed values for M 2 and M 3 can only improve the selection results.
Since each pairing inside J min "defines" a Higgs boson candidate h i 1 , we determine the absolute differences between the invariant mass of each pairing and M 1 , i.e. the mass of the SM Higgs boson. Each one of these differences is sorted from minimum to maximum, (∆m min , ∆m med , ∆m max ). The size of these deviations is an indicator of how accurately the individual SM Higgs bosons are reconstructed. Since ∆m min , ∆m med < ∆m max , the maximum deviation from M 1 is precisely ∆m max . In practice we find that our selection criteria give a distribution for ∆m max which peaks at about 10 GeV in all the signal samples studied.
We also obtain the transverse momentum p T (h i 1 ) of the h i 1 candidate, constructed from the pairings inside J min . These transverse momenta are then ordered from hardest to softest and used as variables for signal and background discrimination. Similarly, we make use of the angular distance ∆R(h i 1 , h j 1 ) between the h 1 candidates h i 1 and h j 1 and additional angular cuts ∆R bb (h i 1 ) are enforced between the b-jet pairs that define each of the h i 1 candidates. The optimisation of the analysis is based on the sequential application of cuts on the different observables described previously, until the significance is numerically above the minimum threshold of 2. More concretely, we obtain the "best" selection cuts for each observable using the following order: (i) p T min,b and |η b |, (ii) χ 2,(6) and χ 2,(4) , (iii) m inv 6b , (iv) m inv 4b . We finally establish the values for the selection cuts affecting the pairings of b-jets which define J min and I min as follows: . The optimisation takes place by constructing a grid over the selection observables and exploring sequentially combinations of cuts which deliver the maximum rejection of the background while maintaining the highest acceptance for the signal. The grid is established by studying the observable distributions to deduce its limits appropriately. Specifically, we look for the maximum and minimum values that capture all the signal events. In the particular case of the invariant masses, bounds from perturbative unitarity pose an additional constraint, which allows us to define the corresponding grid. As an explicit example, the values for p T min,b and the maximum |η b | are obtained by calculating all the possible combinations inside the intervals [25,40] GeV and [1.0, 2.5] over a 20 × 10 grid, respectively. Each possible cut combination is then tested over signal and background and the significance is calculated. At this stage we keep those cut combinations which deliver a significance above 1.5. We then optimise on χ 2,(6) and χ 2,(4) in an analogous fashion, taking as starting values for p T min,b and |η b | from the best pairings obtained in the first stage. At each layer of the optimisation procedure we increase the minimum threshold for the significance. In table 3 we summarise the combination of cuts which give the best performance in our selection procedure.

Results for Triple Higgs Boson Production
In   and 2.5 for the signal, as well as the corresponding selection efficiencies ε, giving the fraction of MC events that pass the cuts. We also show the predicted statistical significances at integrated luminosities of 300 fb −1 and 3000 fb −1 .
Since the number of signal events S and the number of background events B are of the same order, S ∼ B, we employ the following definition of the statistical significance [126] sig (S, B) = 2 [(S + B) ln (1 + S/B) − S]. (5.1) To incorporate the effects of systematic uncertainties, the significance can be estimated according to [126][127][128] where σ B is an estimate of the systematic uncertainty on the total background contributing to this process. We will assume this to have the form σ B = αB, where we will set α = 0.1 to represent a 10% systematic uncertainty on the total background rates. 10 We see that already at an integrated luminosity of 300 fb −1 and in the absence of systematics, significances of up to ∼ 5σ can be achieved for some of the chosen benchmark points. Furthermore, with the full HL-LHC integrated luminosity, all points are within discovery reach, and we obtain significances up to ∼ 16 σ for selected benchmark points. Once systematic errors are taken into account the values for the significance are affected when the background is relatively large. However even for these cases, the significances for 3000 fb −1 are nearly always above 4σ.
In general, the significance that can be achieved is correlated with the h 1 h 1 h 1 production cross sections given in table 2, such that points with higher cross sections have a tendency to lead to higher significances. As production cross sections are directly correlated to the mass M 3 , in general lower masses result in higher significances. For similar masses M 3 , the mass region M 2 ∼ 280 − 300 GeV seems to yield the best results. For parameter points with similar masses for h 2 , on the other hand, significances can largely vary with the production cross section for h 1 h 1 h 1 and/or M 3 , see e.g. points B and F or G and J for comparison, where in each case a smaller mass M 3 /larger production cross section are correlated with higher significance. Note that the semi-automatised cut selection we apply, described in section 4.3, optimises each event sample separately and therefore comparisons in the multivariate parameter space are not straightforward. In a more detailed investigation of points I and F we found, e.g., that a ∼ 6% difference in a cut selection efficiency can increase the difference in significance by a factor 2. A similar behaviour can also be observed in the comparison of points I and E.
In summary, we find that in the region we consider in BP3, significances over 5 σ can already be achieved with an integrated luminosity of 300 fb −1 and that at the HL-LHC all points should be within discovery range. We would like to highlight that our full optimisation strategy and our final results for the significance can be improved by using more sophisticated analysis techniques such as machine-learning multi-variable classifiers. However, in this work we chose not to no pursue such a strategy, since we have demonstrated that it is possible to reach a meaningful threshold for the significance by solely employing an iterative selection procedure.

Other Channels at the HL-LHC
The decay modes of the h 2 and h 3 scalars directly into gauge or Higgs boson pairs can also provide signatures for exclusion or discovery in the BP3 at the HL-LHC. To investigate these, we have extrapolated various analyses assessing the heavy Higgs boson prospects of the HL-LHC in final states originating from h i → h 1 h 1 [54, 57], h i → ZZ [105,124] and h i → W + W − [130,131], for i = 2, 3. We have combined these with extrapolations of results from 13 TeV where appropriate. For further information, see the detailed analysis presented in Appendix D of ref. [22]. The expected exclusion regions for each final state, for an integrated luminosity of 3000 fb −1 are displayed in fig. 2. One can observe that the ZZ final states are by far the most powerful, being capable of excluding almost all of BP3 at the HL-LHC. In addition, the h 1 h 1 final states will achieve an exclusion of a large fraction of BP3. On the contrary, the W + W − final states are foreseen to be rather weak, excluding only a small region of BP3. The significance of the processes in providing exclusion may change in the future if additional decay channels of the gauge or Higgs bosons are considered for each of the processes. Furthermore, detailed experimental studies will be necessary to verify, and potentially improve, our extrapolated observations.
We note that at the HL-LHC, the effects of the TRSM may also be observed through the reduction of the Higgs boson signal strengths. In [124] a lower limit of κ 2 1 min = 0.933 was projected for the so-called S1 scenario [129], where LHC Run 2 systematic uncertainties were assumed. From table 1, we see that BP3 fulfils this requirement and therefore will not be affected by these measurements.

Conclusions
We have examined the triple production of SM-like Higgs bosons, resulting from the asymmetric decay chain pp → h 3 → h 2 h 1 → h 1 h 1 h 1 , within an extension of the SM by two real singlet scalar fields, the TRSM. Our study focused on a specific scenario, "Benchmark Plane 3" (BP3) of [24], where current experimental and theoretical constraints are satisfied on a large portion of the plane of masses of the h 2 and h 3 scalars, (M 2 , M 3 ). We have constructed a Monte Carlo-level phenomenological analysis at the LHC, targeting the 6 b-jet final state originating from the decays of the h 1 scalars. Our analysis demonstrates that at an integrated luminosity of 300 fb −1 , significances of up to ∼ 5 σ can be achieved for some of the chosen benchmark points on BP3. Furthermore, with the full HL-LHC integrated luminosity of 3000 fb −1 , all points that we have considered are within discovery reach, with significances reaching up to ∼ 16 σ. We have also shown that gauge or Higgs boson pair final states of the heavy scalars h 2 and h 3 could probe most of the BP3.
Our results demonstrate that a combination of all of the examined processes of the present article will be essential to discover and gain more insight into the origin of scenarios in which the new physics manifests in a similar manner to BP3. In particular, measurements of the masses of the scalars, the scalar couplings as well as the mixing angles through either single scalar production (pp → h i ), or multi-scalar production such as the pp → h 1 h 1 h 1 process of the present article, will allow measurement of the model parameters and reconstruction of the Lagrangian. This will enable model discrimination and a deeper understanding of the rôle that such new scalars play in Nature, in case they are discovered. Finally, we emphasise the fact that our analysis indicates that the triple Higgs boson final state, thought to be completely hopeless in the past, should be actively pursued at the LHC through concrete experimental analyses by the ATLAS and CMS collaborations.
A Scalar Quartic Self-couplings We define the quartic scalar self-couplings via with i, j, k, l = 1, 2, 3. We then have for a = b = c.

B Total Widths and Branching Ratios
In table 5, we list the total widths as well as decay branching ratios between the physical  scalars of the TRSM, for the benchmark points listed in table 2. The total widths have been calculated according to eq. (2.14), with SM-like widths taken from [125]. Note that the effective branching ratios might vary slightly, as they correspond to BR eff = Γ MG5 x → y z /Γ x , where Γ MG5 x → y z is the respective partial decay width as calculated by MadGraph5 aMC@NLO, while Γ x corresponds to the total decay width, which we here treat as an input parameter. For the benchmark points considered here, we however found that deviations are on the sub-percent level.

C Combinatorics for Scalar Reconstruction
Here we briefly elaborate further on the scalar reconstruction based on the different arrangements of the 6 b-jets with the highest transverse momentum in each event. As discussed in section 4.3, the aim is to determine the combination of two and three pairs of b-jets which minimise the sum χ 2,(6) + χ 2, (4) .
One important aspect of the minimisation is that the set I that defines χ 2,(4) should be a subset of the arrangement J which allows to determine χ 2, (6) .
Here we achieve our target by using the following procedure (C.2) • Subsequently, we obtain all the possible pairings for the full set of 6 b-jets and for each one of them we calculate the corresponding χ 2, (6) . Out of all the possible configurations we select the combination J B min with the smallest value for χ 2, (6) . We label this as χ Note that this procedure assumes that the h 1 bosons are produced on-shell. As discussed in section 4.3, if M 2 and M 3 are such that h 3 can be produced on-shell through the process h 3 → h 1 h 2 and subsequent h 2 → h 1 h 1 , then the configurations I min and J min will ideally correspond to h 2 and h 3 , respectively. Finally, we would like to stress that our optimisation does not assume any a priori values for the masses of the h 2 and h 3 scalars, i.e. M 2 and M 3 .