The Electro-Weak Phase Transition at Colliders: Confronting Theoretical Uncertainties and Complementary Channels

We explore and contrast the capabilities of future colliders to probe the nature of the electro-weak phase transition. We focus on the real singlet scalar field extension of the Standard Model, representing the most minimal, yet most elusive, framework that can enable a strong first-order electro-weak phase transition. By taking into account the theoretical uncertainties and employing the powerful complementarity between gauge and Higgs boson pair channels in the searches for new scalar particles, we find that a 100 TeV proton collider has the potential to confirm or falsify a strong first-order transition. Our results hint towards this occurring relatively early in its lifetime. Furthermore, by extrapolating down to 27 TeV, we find that a lower-energy collider may also probe a large fraction of the parameter space, if not all. Such early discoveries would allow for precise measurements of the new phenomena to be obtained at future colliders and would pave the way to definitively verify whether this is indeed the physical remnant of a scalar field that catalyses a strong first-order transition.


Introduction
One of the biggest scientific questions that could be answered by the largest next-generation experiments is "What was the nature of the cosmological electro-weak symmetry breaking?" [1][2][3]. If the electroweak phase transition (EWPT) was strongly first order, a property that we shall discuss below in detail, signatures would potentially be detectable at space-based gravitational wave observatories such as LISA [2] or Decigo [4], whose sensitivities peak at the millihertz to the decihertz range, frequencies that would occur if a phase transition took place at the electro-weak scale [5][6][7]. 1 The Standard Model (SM) of particle physics by itself predicts the electro-weak transition to be a smooth crossover [12]. Therefore, a strong first-order electro-weak phase transition (SFO-EWPT) dictates physics beyond the SM. Although the nature of the EWPT is an interesting inquiry in itself, the requirement of new physics that enable it can also provide the right conditions to explain the baryon asymmetry of the Universe [13,14]. Next-generation collider experiments are expected to dramatically increase sensitivity to states that interact with the electro-weak sector, making such new phenomena a prime target.
To design a collider that can answer both quantitative and qualitative questions about the nature of the EWPT, it is imperative to address the fact that there exists a multitude of models that can catalyse a SFO-EWPT [10,. Keeping to models that are reasonably minimal, the real singlet scalar field extension of the SM is an ideal test case, as it has no direct gauge interactions, making detection at a collider significantly challenging. Therefore, a collider powerful enough to falsify or confirm the relevant parameter space of this model should also be powerful enough to make a qualitative statement about the cosmological EWPT [1,[52][53][54][55][56].
To investigate the nature of the EWPT, one has to confront the breakdown of conventional analysis techniques of the EWPT [57]. In particular, the scale dependence and gauge dependence of thermal parameters, including the actual strength of a phase transition, can potentially result in theoretical uncertainties large enough to qualitatively alter conclusions derived using such techniques. Here, we strive to include the theoretical uncertainties when approaching the question of what collider specifications are needed to uncover the nature of the EWPT. The dramatic nature of these uncertainties was recently discussed in detail in context of the SM effective field theory [57]. The problems discussed in ref. [57] are even further accentuated if the theory catalysing the EWPT involves large couplings and scale hierarchies.
The prospective high-energy upgrade of the Large Hadron Collider (HE-LHC) and the Future Circular Collider (FCC) at proton centre-of-mass energies of 27 TeV and 100 TeV, respectively, are the focus of our studies. To draw our conclusions, multiple detection channels of the new scalar resonance that appears in the real singlet model were considered. Specifically, parameter-space points that are difficult to detect in decays of the new scalar to two SM-like Higgs bosons, become easier to detect in its decays to vector bosons. Thus far in hadron collider analyses of the real singlet extension of the SM, the latter channels have been mostly neglected. 2 We find that including decays to vector bosons dramatically improves the reach of a collider. Finally, we note that the theoretical uncertainty appears to grow with the mass of the new scalar particle.
The paper is organised so as to summarise our methods and key results in its main portion, while deferring several technical aspects of interest to appendices. The structure is as follows: in section 2 we outline the model that forms the focus of our investigations, the real singlet scalar field extension of the SM. As a complement, appendix A provides further details on the features of the model. In section 3 we describe the method that we employ to calculate the order of the electro-weak phase transition, including the treatment of theoretical uncertainties that leads to our phase-space segmentation into categories which encapsulate varying degrees of theoretical uncertainties. Appendix B provides details on the form of the one-loop potential. In section 4 we provide a summary of the collider information that we employ in our phenomenological analysis, deferring a host of detailed information to appendices C, D, E and F, where we discuss current collider constraints, high-luminosity LHC prospects, electronpositron collider constraints and future hadron collider analyses, respectively. We present the phenomenological results, obtained for a 100 TeV and extrapolated to a 27 TeV collider in section 5. There, we also give a selection of benchmark points in the form of parameter sets with a selection of associated relevant observables. In addition, in appendix G we briefly discuss results originating from the alternative parameter-space parametrisation and in appendix H we study the degree of fine tuning present in the model. We provide our conclusions and outlook in section 6.

Standard Model Augmented by a Real Singlet Scalar Field
We focus on the SM extended by a real singlet scalar field. Then, the most general form of the scalar potential that depends on the Higgs doublet field, H, and a gauge-singlet scalar field, S, is given by (see, e.g. [21,52,[61][62][63][64][65][66][67][68]): where the interactions proportional to K 1,2 constitute the Higgs "portal" that links the SM with the singlet scalar. Note that we do not impose a Z 2 symmetry that would preclude terms of odd powers of S. Such terms are often key in catalysing a tree-level barrier between the electro-weak symmetric and broken phases, thus resulting in a stronger transition. After electro-weak symmetry breaking (EWSB) occurs, the Higgs doublet and the singlet scalar fields both attain vacuum expectation values (vevs) v 0 and x 0 , respectively. To obtain the physical states, we expand about these: H → (v 0 + h)/ √ 2, with v 0 246 GeV and S → x 0 + s. Inevitably, the two states h and s mix through both the Higgs portal parameters K 1 and K 2 as well as the singlet vev and hence they do not represent mass eigenstates. Therefore, upon diagonalising the mass matrix one obtains two eigenstates, 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.1 GeV. We will only consider m 2 > m 1 here. 3 All the couplings of h 1,2 to the rest of the SM states are simply obtained by rescaling:  Figure 1: The branching ratio BR(h 2 → h 1 h 1 ) plotted against BR(h 2 → ZZ), over the parameter-space points of the present study for the various categorisations (see section 3.3). We show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × ×). decay (and hence h 2 → W + W − ) never becomes less frequent than O(few %), it is has the potential to be overall a more powerful probe of the parameter space than h 2 → h 1 h 1 , which can attain small BR values, down to O(10 −3 ). We also note that there are exist viable SFO-EWPT parameter-space points with m 2 < 2m 1 , where h 1 h 1 production is not kinematically allowed. In appendix A we show an alternative view of the complementarity between the h 1 h 1 and ZZ final states, through the ratio of BR(h 2 → h 1 h 1 )/BR(h 2 → ZZ) versus the mass of h 2 , in the bottom-right panel of fig. 8.

Calculating the order of the phase transition
Perturbative methods of treating the effective potential at finite temperature suffer from both gauge [70] and a break down of perturbation theory [71]. The breakdown of perturbation theory arises due to Linde's infamous infrared problem -during the phase transition as the expansion parameter is not the actual coupling but the coupling multiplied by the mode occupation instead: gn B ∼ gT /m (where T is the temperature, n B is the mode occupation, m is a mass and g is some coupling), which clearly diverges for m → 0. A probe of how well perturbation theory is working is to measure the scale dependence of the theory. If perturbation theory is converging quickly, we can define our potential at one loop order and the scale dependence from running should approximately cancel the scale dependence at one loop, any remaining scale dependence is formally higher order. A large sensitivity therefore shows that the theory is particularly sensitive to higher-loop corrections which points to the failure of perturbation theory to converge quickly.
The scale dependence is the dominant theoretical uncertainty if one keeps within the perturbative region of the theory [57]. A partial solution to the scale dependence is to perform a resummation of the masses. Doing so makes the uncertainty in the critical temperature of the SM effective theory non-trivial but manageable, in contrast to the uncertainty in the gravitational wave amplitude generated from a first-order transition which is at the multiorder of magnitude level [57]. In the present paper we consider bosons that couple to the SM with potentially large couplings -in fact past papers that studied the EWPT, considered portal couplings of K 2 ∼ O(10) [52]. For that reason we expect the theoretical uncertainties even in the strength of the phase transition to be very large.
The gauge dependence can be avoided via an expansion in the Planck constant, , instead of the usual loop expansion, but unless one expands to second order in , the scale dependence is unmanageable [72]. 4 Going to two loops is cumbersome and we will not pursue it here as a solution. Another solution is to integrate out the heavy "Matsubara" modes, which results in a theory that is effectively three-dimensional [74][75][76][77][78][79][80][81][82][83][84][85][86]. The dimensionally-reduced theory is manifestly gauge-invariant and includes all-orders of resummation, resulting in a substantial reduction in the uncertainties [57]. For the SM extended by a real singlet scalar field, there can be two light dynamical modes, rendering the analysis a theoretical challenge that we leave to future work.
Here we instead advocate to analyse the electro-weak phase diagram within a given model via an approach that balances convenience, while limiting theoretical uncertainties. Given that we are interested in assessing the order of the phase transition and not gravitational wave phenomenology, the uncertainties are often large but nonetheless manageable 5 in the usual four-dimensional perturbative calculation, where the gauge and scale dependence are treated as theoretical uncertainties.
To assess the question whether a collider rules out a SFO-EWPT, we consider the maximum value of the order parameter, requiring for a given parameter point with physical quantities matched at the Z pole and with the renormalisation scale and gauge parameters varied. In the above, T C is the critical temperature at which the theory has a degenerate ground state and ∆h is the difference in the Higgs field between the two vacua at the critical temperature. By contrast, to assess the discovery prospects we instead look at the minimum value, requiring when varying over the same range of the scale and gauge parameters. It is therefore useful to classify which points might be strongly first order, as they sometimes are when vary scale and gauge parameters, and when they are robust to theoretical uncertainties. The specific implementation of this and how we categorise our points we discuss in section 3.3.

Gauge and Scale Dependence
To capture the theoretical uncertainties coming from the gauge dependence, we construct the effective potential in the general covariant gauge (also known as Fermi gauges) rather than the usual R ξ gauge [87]. The R ξ gauge uses a different gauge for every value of the scalar field(s). It is not a priori clear that this is possible [88,89]. To construct the potential we follow the techniques of ref. [87]. The potential at one loop in the covariant gauge is given by where V tr is the tree-level potential, V CW the Coleman-Weinberg (CW) term, V T the thermal term, Q is the renormalisation scale and the ξ i (i = W, B) are the gauge parameters. The remaining details of the effective potential at finite temperature in the covariant gauge are given in appendix B.
To minimise scale dependence, we use "Arnold-Espinosa" resummation [73] of the masses and allow all parameters in the effective potential to run. We use the program SARAH [90] to derive the one-loop renormalisation group equations (RGEs) and we vary the RGE scale, Q, that defines the coupling and the CW potential, by an order of magnitude: m Z /2 < Q < m Z × 5, where m Z is the Z boson mass.

Numerical calculation of the phase transition
PhaseTracer is a publicly-available C++ package [91] that traces the thermal evolution of the effective potential using the algorithm developed in ref. [92]. It was designed to be robust, 6 with an improved treatment of thermal functions [93] and discrete symmetries.
We have developed a module for PhaseTracer for the real scalar singlet extension of the SM in the covariant gauge with "Arnold-Espinosa" resummation of thermal masses. Since the choice of gauge only enters the effective potential through the field-dependent masses, the covariant gauge can be readily implemented in PhaseTracer by simply providing expressions for these. We can then treat the gauge parameters as inputs we vary.
The requirement that we have a SFO-EWPT typically requires large portal couplings, which implies that the unphysical scale dependence at one loop can become quite large, even at zero temperature. This, in addition to the unphysical gauge dependence of the minimum calculated at one loop, motivates a generous tolerance for the zero-temperature vev value of the Higgs field. In particular, if we wish to claim that a certain collider can exclude the possibility of the cosmological phase transition being strongly first order, we choose to be liberal and not too aggressive in excluding points that could be potential candidates. We discuss our analysis of the phase structure in detail in the next sub-section.

Parameter
Range To derive parameter-space points that satisfy SFO-EWPT, we perform random scans over the five free parameters in the potential. We generate points with viable zero-temperature phenomenology, leaving a large leeway for theoretical uncertainties. Table 1 presents our scan range, which was performed for O(10 6 ) points. We impose a cutoff on dimensionful parameters, at around 2 TeV, as well as on the dimensionless portal coupling K 2 , in order to lie well below the perturbativity bound. This is because the theoretical uncertainties become unmanageable for larger couplings and large scale hierarchies. Note that our scan does not represent a uniform sampling of the parameter space: this implies that our results should not be interpreted so as to indicate densities of points in parameter space, but rather to map out its boundaries. We create a set of eight points for each parameter point with {µ, ξ W,Z } ∈ {(m Z /2, m Z , 2m Z , 5m Z ), (0, 3)}. We require the Higgs boson mass to equal 125.1 GeV at µ = m Z , ξ W,Z = 0 and we impose a restriction on the one-loop value of the sine of the mixing angle: | sin θ| < 0.25. Each of the eight points is individually passed to PhaseTracer for further processing. PhaseTracer calculates the phases and the transitions between them. The transitions are classified according to the value of ρ = φ C /T C , where φ C is the Higgs field vev and T C is the critical temperature.
PhaseTracer has multiple advantages [91] over CosmoTransitions [92]. However, PhaseTracer does not calculate the nucleation rate at present. We therefore consider approximate criteria in categorizing a strong first order phase transition. If there exist 7 phase transitions with 7 This condition is to limit the cases where the tunneling rate is too slow for the phase transition to complete.
We leave an analysis of the tunneling to future work. For now the reader is referred to refs. [29,94] for further discussion on this topic. T c > 30 GeV and ρ > 1, and no other transitions with ρ ∈ [0.1, 1.0] 8 with higher T c than the ρ > 1.0 transition, then we label a point as satisfying the "SFO condition". Additionally, the maximum value of ρ, denoted as ρ max , is kept for each of the eight scale/gauge variations. Furthermore, if for an individual point there exists a Higgs field vev within v 0 = 246±30 GeV, that also corresponds to a transition to the absolute (i.e. deepest) minimum of the potential, we label it as satisfying the "first type of vev condition". The subset of points that satisfy the "SFO condition" are considered for additional post-processing. We also save the value of the Higgs field vev at the deepest minimum for each of these points. During the postprocessing of the points, we also consider an additional property, pertaining to the group of eight scale/gauge variations: if at least one point has for the deepest minimum v 0 < 246 GeV and another has for the deepest minimum v 0 > 246 GeV, we label it as satisfying the "second type of vev condition". We use these conditions to categorise the points. To best take into account the theoretical uncertainties, we consider two types of classification. The first type of classification separates points into: • Ultra-Conservative points: all eight of them satisfy the SFO and all eight also satisfy the first type of vev condition.
• Conservative points: all eight of them have to satisfy SFO and at least one point has to satisfy the first type of vev condition. Ultra-Conservative points are excluded from this category.
• Centrist points: at least one point within the eight satisfying SFO and first type of vev conditions simultaneously. (Ultra-)Conservative points are excluded from this category.
• Liberal points: at least one point within the eight satisfies SFO and any other satisfies the first type of vev conditions. Centrist and (Ultra-)Conservative points are excluded from this category.
The second classification separates points into two mutually-exclusive categories: • Tight points: all eight of them have to satisfy SFO and the group has to satisfy the second type of vev condition.
• Loose points: at least one point within the eight satisfying SFO and the group has to satisfy the second type of vev condition.
In the main part of the article we show results for the first type of classification (GWA-PUCons, GWAPCons, GWAPCentr and GWAPLib as defined above, respectively), with a selection of results for the second type (GWAPTight, GWAPLoose, respectively) shown in Appendix G. Our results show that the "Centrist" and "Loose" categories behave in a similar fashion, whereas the "Liberal" category of the first classification has no correspondence in the second and allows for larger theoretical uncertainties.
We show two-dimensional projections of the parameter-space points obtained in our scan, for the five free parameters that we scanned over in fig. 2. Notice that the theoretical uncertainties get large when there is either a large scale hierarchy (including a large ratio of dimensionful parameters) or a large dimensionless coupling. Specifically if λ s 2 or K 2 4 then the theoretical uncertainties begin to get large. This is precisely what is to be expected if the uncertainties mostly arise due to the unphysical scale dependence coming from the breakdown of perturbativity. In fig. 3 we show the resulting maximum value of ρ max over the eight scale and gauge variations for each of the parameter-space points. In both figs. 2 and 3 we show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × ×), see the next section and the relevant appendices for further clarification of these conditions. We emphasise that, in all of the phenomenological studies that follow, we employ the "central" parameter-space point in the group, i.e. with {µ, ξ W,Z } = {m Z , 0}, irrespective of whether it satisfies the above conditions.

Summary of collider constraints
There exist several categories of collider constraints that affect the parameter space of the real singlet extension of the SM. These come through direct searches for heavy scalar resonances (i.e. h 2 ), SM-like Higgs boson signal strength measurements (i.e. h 1 ) and electro-weak precision observables. To assess the potential of future colliders for exclusion or discovery of the real singlet model in the context of a SFO-EWPT, at the end of the lifetime of the high-luminosity run of the LHC (HL-LHC), we consider the most up-to-date collider constraints coming from these categories.
To impose current constraints coming from heavy Higgs boson searches and Higgs boson measurements, we employ the HiggsBounds [95][96][97][98] and HiggsSignals [99][100][101] packages. In addition, we consider constraints coming from resonant Higgs boson pair production, h 2 → h 1 h 1 and also incorporate the latest 13 TeV ATLAS and CMS SM-like Higgs boson global signal strengths, µ = σ measured /σ SM , that are currently not included in HiggsSignals. Further details on these can be found in Appendix C.
For the HL-LHC, we consider the prospects for the Higgs boson signal strength measurement, as well as various analyses assessing the heavy Higgs boson prospects of the HL-LHC in final states originating from h 2 → h 1 h 1 , h 2 → ZZ and h 2 → W + W − . We combine these with extrapolations of results from 13 TeV where appropriate. For further details, see Appendix D.
Electron-positron colliders may also provide relevant constraints on the model. In particular, electro-weak precision observable measurements (EWPO) coming from LEP and future lepton colliders, such as the International Linear Collider (ILC) can also probe the effects of new scalar particles. To consider these, we follow the treatment of refs. [29,52,53]. We found the resulting constraints to be generally weaker than either the direct heavy Higgs boson searches or the Higgs signal strength measurements and we defer the detailed description and results to Appendix E. There, we also discuss phenomenological results obtained for the Compact Linear Collider (CLIC) [58], which we briefly contrast to the hadron collider results.
To derive the prospects of future hadron colliders to either discover or rule out the real singlet extension of the SM in regards to SFO-EWPT, we have constructed detailed phenomenological analyses at the Monte Carlo level, considering the processes h 2 → h 1 h 1 , h 2 → ZZ and h 2 → W + W − at a 100 TeV proton collider in the mass range m 2 ∈ [200, 1000] GeV. These were inspired by current LHC analyses, with appropriate modifications at 100 TeV. To obtain the HE-LHC constraints we have performed an extrapolation of the 100 TeV results down to 27 TeV. These phenomenological analyses and associated results are discussed in detail in Appendix F.
To estimate the future colliders constraints on the mixing angle through signal strength measurements, we note the several signal strength and coupling measurement projections for various channels that are found in ref. [102]. These point towards a O(1%) uncertainty or better at electron-positron colliders and at a 100 TeV proton collider. In the absence of a full signal strength combination at these colliders, to remain conservative, we will assume the 95% C.L. constraint sin 2 θ < 0.01. We apply this constraint in addition to those coming from direct searches for heavy Higgs bosons at future colliders.
Throughout this article we will only show parameter-space points that pass the current constraints described in detail in Appendix C (denoted by a ), or the HL-LHC constraints described in detail in Appendix D plus future signal strength (| sin θ| < 0.1) constraints (denoted by a ). For completeness, we also show those that are currently excluded (faint × × ×) without indicating the category they fall into.

Theoretical uncertainties
To assess the severity of the theoretical uncertainties at the (zero-temperature) phenomenological level, we have calculated the variation of the one-loop mass of the SM-like Higgs boson with scale, ∆m 1 , taken in [0.5m Z , 2m Z ]. 9 The resulting values are shown in fig observe that the "Conservative" points generally possess a lower variation, ∆m 1 80 GeV, whereas the "Centrist" plus the "Liberal" points possess variations as large as 100 GeV at ∼ 600 GeV, which then reach 150-200 GeV around m 2 ∼ 900 GeV. It is also evident that the current constraints "push" the parameter-space points towards higher theoretical uncertainties, removing a large proportion of the higher-m 2 points with a relatively lower value of the variation.

Proton colliders at 100 TeV
Following the Monte Carlo-level phenomenological analyses for h 2 resonant searches we have developed, presented in Appendix F in detail, we have derived the expected statistical significance at a 100 TeV proton-proton collider, for each of the parameter-space points for a lifetime integrated luminosity of 30 ab −1 . The significances for each of the individual analy- We indicate the 5σ (discovery) boundaries by the red dashed lines. We show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × ×).
evident that the pp → h 2 → h 1 h 1 process can probe a large part of the parameter space, at an integrated luminosity of 30 ab −1 on its own, apart from a limited number of points, mostly at higher values of m 2 . The pp → h 2 → ZZ and pp → h 2 → W + W − channels yield high significances even for points with a large branching ratio h 2 → h 1 h 1 . This is due to the fact that the parameter-space branching ratios for h 2 → ZZ and h 2 → W + W − remain relatively significant, at least O(few %), even for large BR(h 2 → h 1 h 1 ), and due to the power of the gauge boson analyses themselves. In turn, the pp → h 2 → ZZ analysis performs better than the pp → h 2 → W + W − analysis that we have constructed, which fails to yield high enough significances for discovery in the intermediate mass range. The sensitivity of the ZZ analyses is driven by the ZZ → (2 )(2ν) channel, which yields the most stringent constraints on the pp → h 2 → ZZ cross section. The ZZ analysis itself represents by far the highest significance for any parameter-space point.

Proton colliders at 27 TeV
We indicate the 5σ (discovery) boundaries by the red dashed lines. We show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × ×).
Following the extrapolation described in appendix F.2, we show equivalent plots for the significances obtained in the h 1 h 1 and gauge boson channels in fig. 6, for the case of the HE-LHC 27 TeV proton-proton collider with an integrated luminosity of 15 ab −1 . It is evident that a 27 TeV collider could be powerful enough to exclude or discover a large fraction of the SFO-EWPT parameter space, with reductions in significance of O(a few). We emphasise, however, that this result should be validated by a full phenomenological analysis, taking into account the signal and background kinematics at 27 TeV.

Summary
It is evident that both a 27 TeV collider with an integrated luminosity of 15 ab −1 and a 100 TeV collider with 30 ab −1 may be able to exclude or discover the whole of the parameterspace points. We note that since the 27 TeV result originates in an extrapolation down from 100 TeV/30 ab −1 , it is essential that future studies verify this directly through detailed phenomenological analyses that take into account the changes in kinematics at 27 TeV directly. Nevertheless, it is clear that a 100 TeV machine possesses the power to potentially quite efficiently discover the real singlet scalar extension of the SM. The very high statistical significances found all over the parameter space, particularly in the ZZ final state, point to an early discovery. This would allow for further precise measurements of the properties of the heavy scalar to be obtained during the lifetime of a 100 TeV machine, so as to definitively verify whether indeed this is the physical remnant of a real singlet scalar field that catalyses SFO-EWPT.

Selected benchmark points
We present a selection of benchmark points within the "(Ultra-)Conservative" and "Centrist" categories that we have defined in section 3.3. 10 These represent the (very) "low" and "medium" theoretical uncertainty categories of points. Tables 2 and 3 show these points in increasing mass of the h 2 from left to right. The top panel shows the real singlet-extended potential parameters and the bottom panel shows some useful derived quantities: the maximum value of ρ max over the variations of scale and gauge parameters, sin θ, the width of h 2 , the significant branching ratios of h 2 and the gluon-fusion production cross sections of h 2 at 13, 27 and 100 TeV. The cross sections are given at NNLO+NNLL according to the CERN "Yellow Report 4" [104] at 13 TeV and calculated at N 3 LO+NNLL as described in Appendix F.1 via the ihixs program (v2.0) [105] at 27 TeV and 100 TeV. We also show the branching ratio to top quarks, which was not part of our analysis. It is nevertheless of interest, as it could provide the only direct way to measure the coupling of the h 2 to fermions.

Conclusions
We have performed a multi-channel analysis of the electro-weak phase transition at future colliders utilising the real scalar extension of the SM as the test model. Considering decays into vector bosons dramatically improves the reach of colliders to the point that a substantially "weaker" collider, particularly in terms of the integrated luminosity, may be needed to probe the nature of the electro-weak transition than previously thought. This is true even  when theoretical uncertainties are taken into account, dramatically expanding the potentially relevant parameter space. However, it is not clear that the entire relevant parameter space can be probed, as theoretical uncertainties become unmanageable for large scale hierarchies and/or large portal couplings. Improved techniques, such as a higher-loop expansion, an improved resummation or dimensional reduction, will be needed to verify if the real singlet extension of the SM can indeed catalyse the eletro-weak transition within its viable parameter space. Furthermore, an improved analysis could probe larger portal couplings than we were able to consider here. Nevertheless, we emphasise that the present treatment of theoretical uncertainties is generous, and we expect our conclusions to be robust.
Finally, if a new heavy scalar particle is discovered early on at any future collider experi-   ment, as our studies indicate, the challenge would then be to comprehend whether it is indeed the remnant of a new scalar field related to the electro-weak phase transition. This endeavour, part of the so-called "inverse" problem, should be pursued in future investigations.

A Heavy Higgs boson decay modes
Since we concentrate on the scenario m 2 > m 1 , no new decay modes appear for the h 1 . For m 2 ≥ 2m 1 , the h 2 → h 1 h 1 decay mode opens up, with width: where λ 112 is the h 1 − h 1 − h 2 coupling contained in the scalar potential after electroweak symmetry breaking, V (h 1 , h 2 ) ⊃ λ 112 h 1 h 1 h 2 , and m 1 , m 2 are the masses of the h 1 and h 2 particles, respectively. For completeness, we give the full list of the tree-level triple couplings between the scalars h 1 and h 2 , representing terms of the form [106]): where we have defined c θ ≡ cos θ and s θ ≡ sin θ.
The total width of the h 2 scalar is given by: where where Γ SM xx (m 2 ) corresponds to the SM-like width of a scalar boson of mass m 2 into the final state xx, i.e. in the limit λ 112 → 0. The values Γ SM xx (m 2 ) were likewise obtained from the CERN YR4 and interpolated. We note that owing to the current constraints on sin θ, the width of h 2 for viable parameter-space points is small compared to its mass. Therefore, throughout this paper, we have assumed that Γ h 2 m 2 . In fig. 7 we show the one-loop value of sin θ, the most crucial ingredient when calculating the collider production cross sections of h 2 . We note that this can also attain negative values, which dominate the lower-m 2 region of parameter space. Although the sign of sin θ does not affect single scalar production, it has an impact on multi-scalar production processes. The portion of parameter space above | sin θ| ∼ 0.16 is excluded by current experimental results.
In fig. 8, we show the branching ratios of the h 2 for the parameter-space points that we have generated. The "(Ultra-)Conservative" points exhibit a relatively high branching ratio to h 2 → h 1 h 1 , between ∼ 0.2 − 0.75, whereas the "Centrist" plus "Liberal" points span values from O(10 −2 ) up to ∼ 0.9. The lower-right plot of fig. 8 shows the ratio of the rate to h 2 → h 1 h 1 over that to h 2 → ZZ. 11 It is clear that the gauge boson final states can play a crucial in discovering or excluding the existence of h 2 , particularly in the higher-mass regions (m 2 > 600 GeV), where the decay to h 1 h 1 is sub-dominant. We show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × ×).

B One-loop corrections to the effective potential in the covariant gauge
In this appendix we give details about how to derive the effective potential at zero and finite temperature in the covariant gauge (also known as Fermi gauges). We follow the prescription given in ref. [87]. We begin by writing the Higgs doublet components as It then becomes useful to write a vector of the five real scalars, including the singlet component which we denote as s, Similarly, we write a vector of gauge bosons 3) The one-loop Lagrangian can be written as: Taking the determinant of the matrix Σ, it is straightforward to derive the mass eigenvalues that enter the one-loop corrections to the tree-level potential: In the above, the field-dependent top and gauge boson masses take their usual form. The four Goldstone-like scalar masses 12 and two physical modes are, respectively, Finally, the multiplicities for the scalar masses are n 1,+ = n 1,− = 2 and n 2,± = n h,± = 1.
The "Arnold-Espinosa" method involves resumming the bosonic masses which results in the addition of Daisy terms to the effective potential where we have suppressed gauge, field and scale dependence in the arguments. The scalar Debye masses are given by The scalar masses that go into the Daisy resummation term of the potential are simplym 2 = m 2 i,± + Π h for the four Goldstone-like mass terms. For the physical states, the tilded masses are the eigenvalues of the matrix (B.28) 12 They reduce to the Goldstone modes that one derives in the Landau gauge when working in the R ξ gauges.
To remarkable accuracy, one can write the two gauge boson Debye masses as (B.30) C Current constraints

C.1 Current constraints through HiggsBounds and HiggsSignals
To incorporate an array of current collider constraints into our analysis, we employ the Effectively, HiggsBounds and HiggsSignals impose constraints through direct SM-like Higgs boson signals (h 1 ) and direct searches for heavy scalars (h 2 ). The SM-like Higgs boson signals impose constraints through the reduction of the rate via cos θ and the h 2 searches impose constraints on sin θ and the branching ratio to a particular final state, according to eq. 2.3. In HiggsBounds, the most constraining analyses have been found to be the CMS combination of SM-like Higgs boson searches and measurements of its properties using 7 and 8 TeV proton-proton centre-of-mass data corresponding to an integrated luminosity of 5.1 fb −1 and 12. [121], both constructed in the cases of narrow resonances. To construct a constraint for a given parameter-space point of the real-singlet extended SM, we compare to the rescaled 13 TeV cross section calculated at next-to-next-to-to-leading order in QCD with next-to-nextto-leading logarithmic resummation (NNLO+NNLL), which appear in the CERN "Yellow Report 4" [104].

C.3 Current Higgs boson signal strength constraints
The on-peak SM-like Higgs boson measurements directly constrain the mixing angle. For the current measurements, these constraints are implemented in our analysis through the HiggsSignals package. The latest 13 TeV ATLAS and CMS SM-like Higgs boson global signal strengths, µ = σ measured /σ SM have been found to be µ ATLAS = 1.11 +0.09 −0.08 [122] and µ CMS = 1.02 +0.07 −0.06 [123], with 80 fb −1 and 137 fb −1 of data, respectively. These measurements are currently not included in the HiggsSignals constraints. Since µ = µ SM cos 2 θ in the SM extended by a real scalar, the model cannot accommodate measurements of µ > 1. Therefore, the fact that current experimental central values both lie in µ > 1 results in more stringent constraints than expected. In particular, at 95% C.L. (2σ), the ATLAS measurement would imply that sin 2 θ = −0.11 +0. 18 −0.16 and the CMS measurement, sin 2 θ = −0.02 +0.14 −0.12 . Since sin 2 θ > 0, taking only the positive values results in sin 2 θ < 0.07 and sin 2 θ < 0.12 for the ATLAS and CMS measurements respectively, at 95% C.L.. Combining the two measurements by taking their mean and their errors in quadrature results in sin 2 θ < 0.05 for the combined current ATLAS and CMS constraint at 95% C.L.. We impose this constraint in addition to the HiggsSignals constraints.

D.1 HL-LHC Higgs boson signal strength
Both ATLAS and CMS have presented projections for the high-luminosity phase of the LHC, foreseen to collect 3000 fb −1 of data at 14 TeV [124,125]. In ref. [124], the ATLAS collaboration provides a projection on a global coupling modifier, which they call κ, of δκ = 2.2% for the "halved systematics" scenario (table 3 of ref. [124]). Since the signal strength µ is proportional to κ 2 , the expected uncertainty on µ would be δµ = 4.4%. The CMS analysis of ref. [125] does not provide a global signal strength or coupling modifier estimate, therefore we take their "S2" projection for the gluon-fusion Higgs boson production mode, with an expected uncertainty on µ of δµ = 3% ( fig. 3 of [125]). Combining the two in quadrature would result in a total ultimate uncertainty estimate at the HL-LHC of δµ ∼ 2.7%. This would imply a 95% C.L. constraint on the mixing angle of sin 2 θ < 0.054. This is in fact almost identical to the current constraint. This is due to the fact that the current measured central value of the signal strength is µ > 1, providing a more stringent constraint on the mixing angle than expected. Thus, in the absence of new phenomena, as the central value of the signal strength measurement µ → 1, it is likely that the current constraint will become weaker at the HL-LHC.

D.2 HL-LHC Heavy Higgs boson searches
Ref. [126] presented an extrapolation of the results of ref. [127] at 13 TeV, of searches of a heavy resonance decaying to a pair of Z bosons. This was done in the narrow-resonance limit, matching the analysis of the present article. The extrapolation focuses on a single final state ZZ → (2 )(2q) (for representing electrons or muons and q quarks). This final state is more sensitive in the high-mass region at LHC energies and therefore the extrapolation only provides constraints for m 2 ≥ 550 GeV. If m 2 < 550 GeV, we will assume that no information is given by this analysis on h 2 → ZZ. A parameter f VBF is defined as the fraction of the electroweak production cross section with respect to the total cross section. The results are given in two scenarios: f VBF floated, and f VBF = 1. In the expected result, the two scenarios correspond to the gluon fusion (ggF) and vector-boson fusion production (VBF) modes, respectively. Here we consider the floating VBF scenario, corresponding to the ggF mode. We consider the "YR18" systematics scenario presented in ref. [126], where the systematic uncertainties are halved with respect to the values given in ref. [127]. This corresponds to the lower-left plot of fig. 173 in ref. [126]. Here and in the rest of this article, we have "digitised" the constraints of the figure using the package EasyNData [128]. As in section C.2, we compare to the rescaled 13 TeV cross section at N 3 LO.
To estimate the expected HL-LHC limit for m 2 < 550 GeV through h 2 → ZZ, we naively extrapolate the results of the 13 TeV analysis of ref. [129] that employed the ZZ → 4 , → (2 )(2q) and → (2 )(2ν) channels for 35.9 fb −1 of data. This is done simply by rescaling the expected 95%.C.L. cross section by the square root of the ratio of luminosities, L curr. /L HL−LHC , where L curr. is the luminosity used to derive the constraints of the ATLAS and CMS analyses and L HL−LHC = 3000 fb −1 . Further details of this analysis are outlined in section F.1.2, where we re-purpose it for our future proton collider studies.

D.2.2 HL-LHC
The ATLAS note [130] has presented the HL-LHC prospects for narrow-width diboson resonance searches in the W W → ( ν)(qq) final state. The analysis considered a 14 TeV centreof-mass energy and therefore we use the NNLO+NNLL ggF Higgs boson cross sections of ref.
[104] to deduce whether a parameter-space point is excluded or not. For this purpose, we have digitised the lower-left plot of fig. 5 in ref. [130]. As before, the h 2 → W + W − analysis considered relatively high-mass resonances, m 2 ≥ 500 GeV, therefore we will assume that it provides no information on parameter-space points with m 2 < 500 GeV. To accommodate points with m 2 < 500, we extrapolate the ATLAS analysis of ref. [131], where the final state W W → (eν)(µν) was investigated at 13 TeV with 36.1 fb −1 of data. Further details of this analysis are provided in F.1.3.

D.2.3 HL-LHC resonant Higgs boson pair production, h 2 → h 1 h 1
To derive the approximate expected constraints at the HL-LHC originating from searches of resonant SM-like Higgs boson pair production, we perform a naive extrapolation of the ATLAS [120] and CMS [121] analyses considered in section C.2.

E.1 Electroweak precision observables
To incorporate the constraints of electroweak precision observables (EWPO), we follow the procedure outlined in refs. [29,52,53]. The existence of a real singlet scalar field will modify the Higgs field's contributions to the diagonal weak gauge boson vacuum polarisation diagrams and will induce additional contributions. The effect can be characterised via the S, T and U parameters [132] The above expression implies weaker constraints for m 2 ∼ m 1 and small mixing angles. The parameter U has a negligible dependence on the mass of the scalar particle involved, and hence we only consider modifications of S and T , i.e. ∆S and ∆T , as in ref. [53]. The relevant contributions come from the S B and T B functions found in the appendices C.1 and C.2 of ref. [132].
To extract "current" constraints, i.e. those coming from LEP, we performed a correlated χ 2 fit of the Gfitter results [133]: with correlation coefficient ρ ST = 0.91. The χ 2 is then given by: i is the measured difference of the corresponding observable from the SM expectation and the matrix (σ 2 ) −1 is the inverse of the covariance matrix: where σ S and σ T are the uncertainties on the measured differences coming from the global electroweak fit.
In addition, future lepton collider prospects of EWPO measurements are discussed in ref. [133]. Assuming no deviation from the SM expectations, the International Linear Collider (ILC) "GigaZ" option at 90 to 200 GeV e + e − centre-of-mass energy, is expected to yield: ∆S = 0.000 ± 0.017(exp.) ± 0.006(th.), ∆T = 0.000 ± 0.022(exp.) ± 0.005(th.), where the first uncertainty is experimental and the second theoretical. We combine these in quadrature in what follows.  Figure 9: The exclusion regions for a real singlet scalar particle with mixing angle θ and mass m 2 through EWPO precision observables obtained through LEP measurements (left) and future ILC/GigaZ projections (right). The curves correspond to constant χ 2 = 2.30, 6.18, 11.83, giving the 1, 2, 3σ exclusion regions, respectively. The grey dashed lines show the | sin θ| = 0.25 and | sin θ| = 0.1 boundaries, the former corresponding to the upper limit in our scans and the latter to the future collider limit from Higgs boson signal strength measurements. We note that points with sin θ 0.16 are excluded by current signal strength constraints. We show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × × × × × × × ×).
Using the above considerations we have calculated the χ 2 on the (| sin θ|, m 2 )-plane, shown in fig. 9, both using the current LEP results and the ILC/GigaZ projections. The curves correspond to constant χ 2 = 2.30, 6.18, 11.83, giving the 1, 2, 3σ exclusion regions, respectively. It is clear that by comparing the 2σ region excluded by current EWPO constraints to the 2σ region coming from future ones, that the current region appears to be more restricted than the future one. In fact, all points on the parameter space are currently incompatible at more than 1σ with LEP results through EWPO. That is due to the fact that the central values of the S and T parameters extracted LEP measurements are not centred about the SM expectations, thereby increasing the value of χ 2 . It is also clear that the constraints coming from EWPO observables will not be as stringent as the future signal strength measurements. Nevertheless, they could provide information on the overall picture, in case a signal is observed.

E.2 Direct searches for Heavy Higgs bosons at future lepton colliders
A comprehensive study of prospects of constraining the real singlet extension of the SM at the Compact Linear Collider (CLIC) was performed in ref. [58]. The article examined the decays h 2 → V V for V = W and V = Z gauge bosons, as well as h 2 → h 1 h 1 at e + e − centre-of-mass energies of 1.4 TeV and 3 TeV and concluded that CLIC would allow an up to two orders of magnitude improvement with respect to the sensitivity achievable by HL-LHC in m 2 ∈ [250, 1000] GeV. To incorporate constraints in our study, we have extracted the combined h 2 → V V constraints of fig. 7   In table 4 we show the fraction of currently-viable points in m 2 ∈ [200, 1000] GeV that is excluded by each of the CLIC analyses performed in ref. [58], for each of the categories defined in section 3.3. Evidently, at CLIC, the h 1 h 1 analyses are much more powerful in excluding points than the V V analysis. At a 3 TeV CLIC, the h 1 h 1 analysis on its own could exclude the whole of the viable parameter space. Therefore, CLIC should be able to provide strong constraints on the parameter space of the real singlet scalar model, but will most likely be unable to exclude the totality of the viable parameter space on its own, except at 3 TeV centre-of-mass energies. To the best of our knowledge, the prospects of detecting heavy scalar resonances at a 100 TeV proton collider have not been previously examined in detail. To this end, we have performed detailed phenomenological analyses at the Monte Carlo level, of the main decay channels of a heavy scalar resonance with SM-like couplings in the mass range m h 2 ∈ [200, 1000] GeV, assuming a narrow width, i.e. Γ h 2 m h 2 . 13 We have considered both gluon-fusion (GF) production and vector boson-fusion (VBF) production of the h 2 .

F Future proton colliders
All the parton-level events have been generated via the Monte Carlo (MC) event generator MadGraph5_aMC@NLO (v2.7.2) [134,135], with a custom modification of the loop_sm model to incorporate an additional scalar particle and its interactions with the SM particles. Both GF and VBF production channels have been simulated at leading order (LO). In the case of GF, we have calculated approximate next-to-next-to-next-to-leading (N 3 LO) corrections via the ihixs program (v2.0) [105], with renormalisation/factorisation scales set to µ R = µ F = m 2 /2 and the PDF set NNPDF23_nnlo_as_0119 [136]. The corrections were calculated in the Higgs Effective Field Theory, including full NLO quark mass dependence and threshold resummation up to next-to-next-to-leading logarithmic accuracy. 14 In the case of VBF we calculated the next-to-leading (NLO) corrections via fixed-order runs in MadGraph5_aMC@NLO. Next-toleading order matching and multi-jet merging (where appropriate), QCD parton showering, hadronization and underlying event simulation were all performed within the general-purpose MC event generator HERWIG (v7.2.1) [137][138][139][140][141][142][143]. We consider NLO matching through the MC@NLO method [144] and for a subset of processes we employ FxFx merging [145,146]. The Monte Carlo simulations were performed with the PDF NNPDF23_nlo_as_0119 for LO and NNPDF23_nlo_as_0118_qed for NLO. Events were analysed via the HwSim module [147] for HERWIG which saves events in a ROOT compressed file format [148], with jets clustered using FastJet (v3.3.2) [149]. The anti-k T algorithm [150] with a radius parameter R = 0.4 was chosen to be the default jet clustering algorithm, unless otherwise stated.
To capture the detector effects, we restrict the pseudo-rapidity coverage of all objects to |η| < 4 and only consider particles with transverse momentum with p T > 400 MeV as being detectable. We smear all the final state reconstructed object momenta according to ∆p T = 1.0 × √ p T for jets, ∆p T = 0.2 × √ p T + 0.017 × p T for photons [52], with p T in GeV. Muons and electron momenta are smeared according to [151]. The lepton identification efficiency is assumed to be (p T ) = 0.85 − 0.191 × exp(1 − p T /(20 GeV)) for electrons and = 0.54 for muons in |η| < 0.1 and = 0.97 otherwise. The efficiency for measuring a jet was taken to be (p T ) = 0.75 + 0.2 × p T /(30 GeV) [151,152]. The tagging of heavy-flavour jets is discussed, where applicable, in the individual analyses outlined below.
Following the construction of a set of observables for each analysis, we employ ROOT's "Toolkit for Multivariate Data Analysis" (TMVA) to obtain the optimal signal versus background discrimination. We use the "Boosted Decision Tree" (BDT) classifier method with AdaBoost. 15 To further reduce dependence on the size of our Monte Carlo samples, we have run the BDT classifier a number of times where appropriate, and obtained the mean signal cross section that yields the desired significance. We have requested that at least 100 signal events remain in all analyses.
The statistical significance, S, of the analysis is calculated as: where σ S is the signal cross section, σ B is the sum of the background cross sections, S = σ S L and B = σ B L are total signal and background events at an integrated luminosity L for a specific analysis channel, after the BDT classification has been applied and α is a fractional systematic uncertainty on the background processes, aiming to approximate the individual process systematic uncertainties. We do not include any systematic uncertainty on the signals, as it is very likely to be sub-dominant at the time of analysis of any future collider results. We consider α = 0.05 as the systematic uncertainty estimate on all the backgrounds throughout our analysis.

F.1.2 pp@100 TeV h 2 → ZZ
The final states arising through the decay of a heavy scalar to Z bosons, h 2 → ZZ, provide a relatively clean avenue to discovery at hadron colliders, usually ranking as the second or third branching ratio. We have adapted the CMS analyses of [127] to 100 TeV, where the final states ZZ → (2 )(2ν), ZZ → 4 and ZZ → (2 )(2q) were examined.
The event selection for this final state consists of combining di-lepton Z boson candidates with a relatively large missing transverse momentum ( / p T ). We require two oppositely-charged leptons of the same flavour, each with p T ( ) > 50 GeV. We further require their combined invariant mass within 30 GeV of the Z boson mass and di-lepton transverse momentum, p T ( ) > 55 GeV. In addition require p miss T > 125 GeV. We veto events if ∆φ( / p T , any jet with p T > 30 GeV) < 0.5, where ∆φ is the difference in angle between the / p T and any jet on the plane perpendicular to the beam axis. We also require the Z boson candidate to satisfy ∆φ(Z, / p T ) > 0.5. We construct the transverse mass as: where m( ) is the invariant mass of the di-lepton system. The final set of observables that are used in the discrimination of signal versus background consists of: • the transverse momenta of the leptons that form the Z boson candidate, p T ( 1 ), p T ( 2 ), • the corresponding, di-lepton invariant mass, m( ), and transverse momentum, p T ( ), • their pseudo-rapidity distance ∆η = |η( 1 ) − η( 2 )| and their distance ∆R = ∆η 2 + ∆φ 2 , • the transverse mass m T as defined above and • the magnitude of the missing transverse momentum, / p T . As backgrounds, we consider those that can yield the 2 final state with an associated missing transverse momentum, originating from the on-shell production of ZZ, W Z, ZV V where V = W, Z, tt and W W production, all matched via the MC@NLO method to the parton shower. We do not consider the mis-identification of jets or photons as leptons, and we do not include τ leptons in either signal or backgrounds, here and whenever we consider leptons in the following analyses. These would contribute additional sensitivity when considered in future analyses.
h  → ZZ → ( )(q): The leptons that form one of the Z boson candidates are required to satisfy the same constraints as in the (2 )(2ν) analysis described above, with the only difference being p T ( ) > 100 GeV instead. In addition to the jets with R = 0.4, we cluster "fat" jets with the anti-k T algorithm with R = 0.8. The hadronic Z candidate can thus be formed either from the "resolved" R = 0.4 jets or through a fat R = 0.8 jet. The fat jets undergo "pruning" [154], with parameters β = 1.0, z cut = 0.1 and r cut = 0.5 and the ratio of the observables "2subjettiness" to "1-subjettiness", τ 21 is calculated for each jet [155]. For the fat jets that will form Z boson candidates, Z fat had we require that τ 21 < 0.6, m(Z fat had ) ∈ (70, 105) GeV, p T (Z fat had ) > 170 GeV, and the distance from any lepton to satisfy ∆R( , Z fat had ) > 0.8. For the candidates formed through the resolved R = 0.4 jets, Z res had , we require m(Z res had ) ∈ (40, 180) GeV and p T (Z res had ) > 100 GeV. We then choose the "best" candidate through the following procedure: • If there exists a fat jet candidate with p T > 300 GeV and p T ( ) > 200 GeV then it takes precedence.
• Otherwise if there's no merged candidate with precedence, move to the resolved jets and pick the highest-p T one.
• If there are no resolved candidates, but there are fat ones, pick the highest-p T one.
Contrary to the analysis of [127], we do not impose any cut on the invariant mass of the ZZ system. The final set of observables used in the BDT consists of: • the transverse momenta of the leptons that form the leptonic Z boson candidate, p T ( 1 ), p T ( 2 ), • their invariant mass m( ) and transverse momentum, p T ( ), • their pseudo-rapidity distance ∆η = |η( 1 )−η( 2 )| and • their distance ∆R = ∆η 2 + ∆φ 2 , • the invariant mass of the hadronic Z boson candidate, m(Z had ), • its transverse momentum p T (Z had ), • the transverse momentum of the combined ZZ leptonic and hadronic candidates, p T (Z lep Z had ), • the invariant mass of their combination m(Z lep Z had ), and • the distance between them, ∆R(Z lep , Z had ). As backgrounds we consider those originating from the Z+jets events, ZZ, W Z, W W and tt. The Z+jets backgrounds have been matched/merged via the FxFx method, whereas the rest have been matched via MC@NLO.
h  → ZZ → ( ): Events are considered if they contain four leptons with transverse momenta satisfying, from hardest to softest, at least: p T ( 1,2,3,4 ) > 50, 50,30,20 GeV. Further, events are only accepted if they contain two pairs of oppositely-charged same-flavour leptons and these are combined to form the Z boson candidates, with the constraint m( ) ∈ [12,120] GeV. If an event does not contain at least two Z boson candidates, it is rejected. In the case of four same-flavour leptons, if there exist two viable lepton combinations, the combination ( i j )( k l ) with the lowest value of χ 2 = (m( i j )−m Z ) 2 +(m( k l )−m Z ) 2 is chosen, forming the candidates Z 1 and Z 2 . We require that the combined invariant mass of the four leptons satisfies m( i j k l ) > 180 GeV.
The final set observables consists of: • the lepton transverse momenta, p T ( 1,2,3,4 ), • the combined lepton invariant mass m( i j k l ), • the transverse momentum of the two Z boson candidates, p T (Z 1 ), p T (Z 2 ), • their invariant masses, m(Z 1 ), m(Z 2 ), • their distance ∆R(Z 1 , Z 2 ), • the distance between the leptons that form the two candidates, ∆R( i , j ) and ∆R( k l ), • the invariant mass of the combined Z boson candidates m(Z 1 Z 2 ) and • their combined transverse momentum, p T (Z 1 Z 2 ). We consider only the dominant backgrounds, originating from non-resonant and resonant SM four lepton production, matched at NLO via the MC@NLO method. In addition, we consider the LO gluon-fusion component of four lepton production that originates from the resonant loop-induced production of two Z bosons, i.e. gg → ZZ, deemed to be important at higher proton-proton centre-of-mass energies [156]. To obtain the best limit on σ(h 2 ) × BR(h 2 → ZZ), we pick the most constraining analysis at each value of m 2 . This is in fact represented by the h 2 → ZZ → (2 )(2ν) analysis in the mass range we consider, as can be observed in fig. 10. A combination of the channels should yield a more stringent limit, but we omit it in the absence of a full study of the correlations between channels. Furthermore, the (2 )(2q) and 4 final states will allow for reconstruction of the h 2 and therefore a precise determination of its mass in case of discovery.
For the analysis of the W + W − final state we follow the ATLAS analysis of [131] focusing on the W W → (eν)(µν) final state. We require two oppositely-charged different-flavour leptons with transverse momenta p T ( 1,2 ) > 80, 60 GeV with combined invariant mass m( 1 2 ) > 60 GeV and maximum pseudo-rapidity difference |η( 1 )−η( 2 )|. We veto events that contain additional leptons with p T ≥ 15 GeV. To reject the backgrounds originating from tt, we veto events containing b-jets with p T > 30 GeV. The b-jet tagging probability was chosen to be 0.75. We also allow for the mis-identification of jets as leptons, with flat probability 5 × 10 −3 . For the events that pass the above constraints, we construct the transverse masses of the two W i 's (i = 1, 2) as: where p T ( i ) are the lepton transverse momenta, / p T is the missing transverse momentum, φ i and / φ are the azimuthal angles of the lepton i and the missing transverse momentum vector, respectively. We require both of these observables i = 1, 2 to satisfy m W i T > 60 GeV. We also construct an event transverse mass via: where the combined transverse energy of the leptons was defined as E T ( ) = | p T ( )| 2 + m( ) 2 . The final set of observables consists of: • the lepton transverse momenta, p T ( 1,2 ), • the di-lepton transverse momentum, p T ( ), • the di-lepton invariant mass m( ), • the distance between the leptons ∆R( 1 , 2 ), • their pseudo-rapidity distance ∆η = |η( 1 )−η( 2 )|, • the W transverse masses m W i T , • the event transverse mass, m T , and • the magnitude of the missing transverse momentum, / p T . As background processes we have considered the SM W W , tt, ZZ, W Z processes, simulated at NLO via MC@NLO and the W +jets process simulated via FxFx merging at NLO.
We show the resulting expected (5σ, red dashes) and exclusion (2σ, blue dashes) cross sections times branching ratio, σ(h 2 )×BR(h 2 → W + W − ), in fig. 11. The constraints obtained through the h 2 → W + W − → eνµν are comparable to those coming from the h 2 → ZZ analyses. These could be improved in future studies, e.g. by including hadronic decay modes of the W bosons. The non-resonant h 1 h 1 final state has been investigated at a 100 TeV collider in [157]. Here we investigate resonant pp → h 2 → h 1 h 1 production at 100 TeV. We focus on by far the most sensitive final state, h 1 h 1 → (bb)(γγ). Future studies may include other final states to improve on this, see e.g. [158][159][160]. We require all jets (including b-tagged) to have transverse momentum p T > 30 GeV and to lie within |η| < 3.0. The b-jet tagging probability was set to 0.75, uniform over the transverse momentum, as for the previous analysis. The jet to photon mis-identification probability was set to 0.01 × exp (p T /30 GeV), where p T is the jet transverse momentum [161]. We require that the invariant mass of the two b-jets lies in m bb ∈ [100,150] GeV and that the invariant mass of the di-photon system within m γγ ∈ [115, 135] GeV.
The final set of observables constructed for the BDT consists of: • the invariant mass of the two b-jets, m bb , • the invariant mass of the di-photon system, m γγ , • the invariant mass of the combined system of the two b-jets and the photons, m bbγγ , • the distance between the b-jets, ∆R(b, b) • the distance between the photons, ∆R(γγ), • the distance between the two b-jet system and the di-photon system, ∆R(bb, γγ), • the transverse momentum of each b-jet, p T (b 1 , 2), • the transverse momentum of each photon p T (γ 1 ), p T (γ 2 ), • the transverse momentum of the two b-jet system, p T (bb), • the transverse momentum of the di-photon system p T (γγ), the transverse momentum of the combined b-jet and photon systems, p T (bbγγ) and • the distances between any photon and any b-jet, ∆R(b i , γ j ) with i, j = 1, 2. As backgrounds we consider γγ+jets, γ+jets, by producing, respectively, γγj and γjj via MC@NLO, ttγγ via MC@NLO, bbγγ and bjγγ at LO. We also consider backgrounds originating from single Higgs boson production: bbh 1 , Zh 1 , tth 1 , where we assume that the branching ratios possess their SM values. We also consider the non-resonant part of h 1 h 1 as a background, assuming that the self-coupling maintains a value close to the SM value. In the case of discovery of a new scalar particle, the full interference pattern should be considered in a more detailed analysis. We leave this endeavour to future work. We show the expected discovery (5σ, red dashes) and exclusion (2σ, blue dashes) cross sections multiplied by the branching ratio of h 2 → h 1 h 1 in fig. 12. The magnitude of the constraints is comparable to that of the gauge boson analyses. We note here that the same process was also considered in [52] in the context of the real singlet extension of the SM. Our results appear to be comparable but somewhat different. This can be attributed to differences in the analysis approach: we have calculated the K-factors for both signal and backgrounds, we have considered the vector-boson fusion process as part of the signal process as well as a different set of observables for which we ran the BDT. We also use the HERWIG Monte Carlo event generator, hence a completely different "tune" for the parton showering and the description of non-perturbative effects. We expect the aforementioned facts, in conjunction with different parton density functions for the incoming proton beams, that contain rather large uncertainties when extrapolated to higher energies, will introduce large differences in the significances, at the level observed between our results and the results of [52]. We expect all of these issues to be evaluated in detail at the onset of operation of any future collider through more precise phenomenological analyses.

F.2 Extrapolation of Heavy Higgs boson searches to 27 TeV
To obtain an estimate of the performance of a potential upgrade of the LHC to 27 TeV centreof-mass energy, we consider an extrapolation of the detailed analyses obtained in the previous sections. We write the collider energy (E) dependence of the cross sections in the significance, given in Eq. F.1, as To extrapolate the obtained limits through our detailed analyses to a different energy E and luminosity L , we solve Eq. F.1 for σ S (E ): Under the approximation that the impact of the analysis cuts remains the same between 100 TeV and 27 TeV, to obtain a limit on the signal cross section at the new energy E we only need to determine the change in the total background cross section σ B (E) → σ B (E ) with energy. We achieve this by assuming that the background cross sections scale in the same way as either the gg → h 2 process or a (hypothetical) qq → h 2 process. Since the analysis cuts will emphasise regions of phase space in which the backgrounds behave similarly to the signal, this should be a reasonable approximation. We assume that the gg and qq scalings provide the error estimate on this part of our extrapolation. In the derivation of the actual limits, we use the one that provides the weakest signal constraints, so as to provide a conservative performance of a 27 TeV machine under this procedure. Another noteworthy assumption is that there will remain a sufficient number of signal events at 27 TeV after the analysis cuts are applied. This assumption ignores the fact that the extrapolation is performed from higher centre-of-mass energy and luminosity. This leads to potentially optimistic results and these should be assessed in a fully-fledged phenomenological analysis at 27 TeV. The results for the 27 TeV extrapolation of the individual channels are shown in fig. 6. The ZZ final state represents the highest significance all over parameter space as for the case of the 100 TeV collider.
G Results for the alternative categorisation of parameter-space points Figures 13 and 14 present a selection of results for the alternative categorisation of "Tight" and "Loose" points described in section 3.3, found during the parameter-space scans via PhaseTracer. The results show that the "Centrist" and "Loose" categories behave in a similar fashion, whereas the "Liberal" category of the first classification has no correspondence in the second and allows for larger theoretical uncertainties.

H Fine-tuning studies
Any parameter-space point has to contain a light Higgs boson, h 1 , close enough to 125.1 GeV, if it is to be phenomenologically relevant. It is therefore reasonable to investigate whether the introduction of additional interactions requires fine tuning of the parameters to achieve the observed Higgs boson mass. To do this, we employ the Barbieri-Giudice measure [162] (see also [163]): for a set of observables O i and a set of model parameters p j . One can roughly interpret the logarithm of ∆ as the number of significant digits that need to be tuned in at least one of the parameters of the theory. We are interested in fine tuning in the Higgs boson mass and hence we consider derivatives of the one-loop expression for O = m 1 with respect to the relevant free parameters in this case.
We consider this quantity for each of the parameter-space points and plot against the oneloop value of m 2 at renormalisation scale µ = m Z . The derivatives have been calculated numerically.
The results are shown in fig. 15. It is evident that the parameter-space points that we consider do not possess substantial fine tuning. This is in line with the tree-level studies of We indicate the 5σ (discovery) boundaries by the red dashed lines. We show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × ×).
ref. [163], where it was pointed out that it is possible to have a large hierarchy of scales in theories with extended Higgs sectors, without having fine tuning.  Figure 15: The value of the one-loop fine-tuning measure plotted against the one-loop mass of h 2 , m 2 . We show points that pass the current constraints ( ) or the HL-LHC constraints plus future signal strength (| sin θ| < 0.1) constraints ( ), as well as those that are currently excluded (faint × × ×). [109] ATLAS collaboration, Analysis of ttH and ttW production in multilepton final states with the ATLAS detector, .
[110] CMS collaboration, Measurement of the associated production of a Higgs boson and a pair of top-antitop quarks with the Higgs boson decaying to two photons in proton-proton collisions at √ s = 13 TeV, .
[ [114] CMS collaboration, Measurement of ttH production in the H → bb decay channel in 41.5 fb −1 of proton-proton collision data at √ s = 13 TeV, .
[115] CMS collaboration, Evidence for associated production of a Higgs boson with a top quark pair in final states with electrons, muons, and hadronically decaying τ leptons at