Fusing Vectors into Scalars at High Energy Lepton Colliders

We study vector boson fusion production of new scalar singlets at high energy lepton colliders. We find that CLIC has the potential to test single production cross-sections of a few tens of attobarns in di-Higgs and di-boson final states. In models with a sizeable singlet-Higgs mixing, these values correspond to a precision in Higgs couplings of order 0.1% or better. We compare our sensitivities with those of the LHC and interpret our results in well-motivated models like the Twin Higgs, the NMSSM and axion-like particles. Looking forward to even higher energy machines, we show that the reach of muon colliders like LEMMA or MAP overcomes the one of future hadron machines like FCC-hh. We finally study the pair production of the new scalar singlets via an off-shell Higgs. This process does not vanish for small mixings and will constitute a crucial probe of models generating a first order electro-weak phase transition.


Introduction
The legacy of the LHC on the Standard Model (SM) Higgs boson will not be enough to fully assess the structure of the Higgs sector. This provides a strong motivation for future lepton colliders, which is somehow independent of other possible LHC findings. An additional important motivation for these machines comes from their potential to directly produce and test new physics states. It is the purpose of this paper to study the interplay of these two motivations in a concrete yet general model, focussing on high-energy lepton colliders (HELCs).
The different proposals of lepton colliders can be classified in low-energy, such as FCC-ee [1,2], CEPC [3,4], and ILC in its current design [5,6], and high-energy, such as CLIC [7] in its stages at a center-of-mass energy of 1.5 TeV (Stage II) and 3 TeV (Stage III). Other futuristic examples of HELCs are high-energy circular muon colliders to be possibly built at CERN [8] and/or at the muon accelerator facility at Fermilab [9]. Muons could be produced from p scattering on a target (MAP [9]) or e + scattering on a target (LEMMA [10][11][12]), see Ref. [12][13][14] for the attainable luminosities with these technologies. Even more futuristic HELCs are linear electron colliders like AWAKE [15] and ALEGRO [16], where the electrons are accelerated through proton-driven plasma wakefield acceleration PWFA (see for example Ref. [17], and Ref. [13] for the related luminosity). We summarise in Figure 1  Low energy e + e − colliders are known to be wonderful machines to achieve very precise indirect measurements of the properties of the Higgs sector, well beyond the reach of HL-LHC (see Ref. [18] for a recent discussion). As a concrete example FCC-ee could probe Higgs couplings to SM gauge bosons at the 0.1% level, by collecting 10 ab −1 at the peak of Higgs-strahlung rate [19]. In HELCs like CLIC the production cross-section for the SM Higgs will instead be dominated by the W W -fusion process, with a rate given in the high-energy limit by (see e.g. Ref. [20]) where s is the center of mass energy, v = 246 GeV and m h 125.1 GeV is the Higgs boson mass. Since the W W -fusion rate depends only logarithmically on s, while the Higgs-strahlung rate is suppressed by 1/s, the cross-section in Eq. (1) gets larger than σ(eē → Zh) for √ s 400 GeV. Therefore, by means of the W W -fusion process, HELCs can partially overcome the smaller luminosity compared to the low energy colliders and obtain a similar precision in the Higgs couplings to gauge bosons, as discussed in Ref. [21] (see also Ref. [19]).
The purpose of this work is to assess the HELCs capabilities to discover new physics directly, i.e. to produce and detect new particles, in particular if weakly coupled to the Standard Model and therefore at the edge of, or beyond, the reach of the HL-LHC. In order to discuss the interplay between precision tests of the Higgs sector and the direct exploration of new physics, we focus on models where new physics is coupled to the SM Higgs sector. To keep the discussion minimal, but retaining all the features discussed so far, we consider a model where only a new real scalar singlet is added to the SM.
The new state mixes with the Higgs such that i) it induces a tree-level shift in the Higgs couplings, ii) it can be singly produced analogously to the Higgs, and iii) it can be pair-produced via the quartic Higgs portal coupling. Having set the framework, the main question we want to address in this paper is the following: To what extent will CLIC, and HELCs in general, directly test new physics inducing deviations in the Higgs couplings to gauge bosons of a few percent or smaller?
A first answer to this question comes from the study of single production. From a full collider study, focusing on the dominant decays of the scalar singlet into di-Higgs and di-boson final states, we find that searches for resonances in four b-jets are going to be the dominant discovery/exclusion channel at CLIC. 1 Interestingly, CLIC can test new resonances well beyond the capabilities of HL-LHC and down to couplings correlated to a deviation in the Higgs couplings smaller than 0.1%. We then carry out a similar analysis at muon colliders with center-of-mass energies of 6 or 14 TeV, like LEMMA or MAP, and find that direct searches for extra scalars would be more effective than the corresponding ones at a 100 TeV pp collider. We then recast our results in two explicit models featuring an extra scalar singlet, the NMSSM [23] and Twin Higgs [24], that are well-motivated as they address the hierarchy problem of the Fermi scale. We also comment on to what extent di-boson searches at HELCs can probe the heavy mass regime of an axion-like particle that couples only to electroweak (EW) gauge bosons, a scenario which is notoriously challenging from the phenomenological point of view [25][26][27].
A second answer comes from the study of pair production. This becomes very important when the coupling controlling the single production is parametrically suppressed, as it happens when the new scalar is odd under an approximate Z 2 symmetry. In this limit, pair production of the new scalars constitutes a direct test of the Higgs portal coupling. Models of this type have been considered as interesting benchmarks for electroweak baryogenesis (see Ref. [28] for a review), because the coupling of the singlet with the Higgs can induce a first order electroweak phase transition (FOEWPT), see e.g. Ref. [29][30][31]. Depending on the decay length, the singlet pair production can lead to final states with multiple gauge or Higgs bosons, with displaced vertices in the tracker/muon chamber, or with missing energy. We find that, in all of the cases above, CLIC (and a fortiori muon colliders) will sensibly ameliorate the HL-LHC reach and, most importantly, has the potential to entirely probe the region where a FOEWPT is possible. This paper is organised as follows. In Section 2 we introduce the model and the relevant formulae. Section 3 is dedicated to the discussion of the collider study for single production, and Section 4 to its application to concrete models. In Section 5 we present our findings for double production. We conclude in Section 6.

The Singlet Model
To set our notation, we consider here an extra real CP-even scalar degree of freedom coupled to the standard model only via renormalisable interactions 1 See also Ref. [22] for a detailed study of the search in the 4b final state at ILC, and an estimate for the 1.5 TeV CLIC.
We define the mixing angle γ as the rotation needed to go from the basis of Eq. (2), where only the Higgs couples directly to the SM fields, to the mass basis where h is the SM-like Higgs with mass m h = 125 GeV, and φ is the singlet-like state with mass m φ . We now highlight the main phenomenological consequences of the Lagrangian in Eq. (2), discussing both deviations in the SM-like Higgs couplings, and single and double production of the new scalar φ. For definiteness, in this paper we only consider the mass ordering m φ > m h .
SM-like Higgs boson. The main deviation in the Higgs couplings to vectors and fermions is generated at tree-level by the mixing γ: the Higgs signal strengths µ h are universally rescaled as When a HS and λ HS are both non-vanishing, the above deviations are uncorrelated from those in the trilinear Higgs coupling, that can in principle be larger. Under favourable circumstances, the HL-LHC could even observe deviations in double Higgs production without observing any in the Higgs couplings to SM fields, see e.g. Ref. [32]. An accurate description of the Higgs sector in our setup can also be achieved by integrating out the singlet field and computing the Wilson coefficients of the dimension-6 operators where λ H is defined as the coefficient of the |H| 4 operator. These operators predict the following relative shifts in the Higgs couplings where v = 246 GeV and g SM hhh = 3m 2 h /2v (see also Ref. [33] for a related discussion at lepton colliders). In the singlet model, the tree level and one-loop contributions to these operators read Notice that we assume negligible cubic terms for the singlet to compute the first contribution to c 6 (see Ref. [34,35]). This shows the importance of sin 2 γ in the low-energy phenomenology of the Higgs, and the possible interplay between a small mixing and a large portal coupling λ HS to get visible effects in the triple Higgs coupling.
Singlet-like φ boson. The phenomenology of φ depends largely on the presence of a non-vanishing mixing angle (i.e. of a breaking of the Z 2 symmetry), given that in this case it can be singly produced, and it can decay to SM particles. Single production cross-sections and SM decay rates for the singlet-like state φ are proportional to sin 2 γ, and read where σ h (m φ ) is the production cross-section for a SM Higgs boson of mass m φ and in the second equation we assumed the absence of non-SM decay modes. From Eq. (9) we see that the heavy singlet dominantly decays into W , Z, and Higgs bosons. The branching ratio into hh is in principle a free quantity that depends on the parameters of the scalar potential, but for heavy scalars m φ m W the potential exhibits an approximate SO(4) symmetry which implies BR φ→hh BR φ→ZZ BR φ→W W /2.
Double production of the singlet can instead proceed at tree-level even for γ = 0, through the portal coupling λ HS . We show in Section 2.2 the explicit dependence of the production rates on the masses and couplings of the model.

Scaling of physical quantities
The actual size of the parameters appearing in Eq. (2) depends on the concrete UV model under consideration, in such a way that measuring single and/or double production could provide a hint about the underlying dynamics. In order to maintain a simple description, before specialising to concrete cases in Section 4, we consider a new sector with two intrinsic parameters: a mass scale M * and a coupling g * (see also Ref. [36]). In this case we can write where the proportionality constants are numerical coefficients of O(1), unless constrained by additional symmetries and/or scales. The squared mass matrix of the two states (h 0 , S) reads where s is the vacuum expectation value (VEV) of the singlet. Under the assumption of a small mixing angle as suggested by data, γ M 12 /M 22 . The size of γ depends on the origin of the singlet VEV s. We can identify three different scenarios: 1. If the singlet develops a VEV because of its potential alone, then s ≈ M * /g * , so that the mixing angle scales as times a combination of O(1) numbers (that can be tuned to be small).
2. If the singlet gets its VEV due to the interaction with the Higgs dynamics, then s ≈ a HS v 2 /g 2 * M 2 * . In this case, the mixing angle is controlled by the size of the explicit breaking of the S → −S symmetry, where the dimensionless parameter δ = a HS /g * M * can be made arbitrarily small.
3. If, in addition to the previous case, the dynamics responsible for the Z 2 breaking terms is related to a mass scale independent of M * , then the mixing decouples as M −2 * and not as M −1 * . These situations are realised in different concrete examples discussed in the literature: Models with a moderate coupling g * with a HS = 0 and λ HS ∼ g 2 * . In this case we expect a deviation γ ≈ g * v/M * , so that the only way to comply with the bounds is by M * /g * ≡ f v. This type of decoupling is a well known feature of Twin Higgs models (see Section 4.2). Models with a weak coupling g * , such that γ ≈ g * v/M * 1 even for light states. This is the case in the NMSSM, once we identify M * with the SUSY-breaking mass of the singlet (M * =m) and g * with the coupling in superpotential (g * = λ). The only way to additionally suppress the mixing angle is to invoke a tiny s, which is achievable by neglecting the A κ S 3 + h.c. soft term and by allowing a HS ≡ A λ sin(2β) λm (see Section 4.1).
Notice also that the bounds obtained at the kinematic edge of the lepton collider, where γ quickly approaches O(1) for large masses, could be interpreted in terms of strongly coupled new physics. This region however is (and will be) strongly constrained by single Higgs production.

Vector boson fusion
As discussed in the introduction, the advantage of HELCs is mainly due to the effectiveness of vector boson fusion as a production mode for scalar particles. Both single and double productions can be written in terms of the cross-section of the subprocess V V → φ and V V → φφ properly convoluted with the splitting functions for → V . Any differential distribution for the process eē → ννX can be written as a distribution in the invariant mass squared of the subprocesses as where we defined the effective parton luminosities C V i V j in terms of the splitting functions f V i (x). These can be computed analytically in the regime M 2 V /ŝ 1 [37,38]. Here we focus on the longitudinal polarisations, which are the only ones coupled to the extra singlet through the mixing with the SM Higgs: By inspecting the behaviour at high s, we see that the total rate of W W -fusion does not fall with energy neither for single nor for double singlet production. The total rates can be computed to be where Eq. (17) holds in the limit sin γ = 0. The formulas in Eq. (16)- (17) are extremely good approximations as long as the dominant contribution to the rates comes from kinematic configurations where M 2 V /ŝ 1. We checked that they reproduce with excellent accuracy the full result, which we compute with MadGraph5 [39,40]. This is shown in Figure 2 and we use it in all our numerical calculations.
Here and in what follows, we assume unpolarised electron beams.
The above expressions for the production rate show explicitly what is well known: W W -fusion is a powerful production channel for HELCs. At increased center-of-mass energy, other production mechanisms such as φ-strahlung and double φ-strahlung are subdominant, because they are suppressed at large s (see also Ref. [20] for a comparison). Based on these considerations we motivate our approach of just considering V V -fusion processes for the production of the scalar singlet. In our study we do not include next-to-leading orders in EW radiation, thus an uncertainty of the order of α 2 4π log(s/m 2 should be understood in all our sensitivities. 2 This is safely below the 3% level even at the 14 TeV stage of future µ-colliders.

Single production
In this section we assess the capabilities of HELCs to test the existence of new scalar particles by means of their single production in W-fusion. The total production rate as a function of the mass of the scalar has been computed in the previous section, and is displayed in the left panel of Figure 2. The dominant decay channels of φ are into pairs of vector bosons and Higgs bosons, as given in Eq. (9). We are going to study resonant production modes, in narrow-width approximation and with only visible final states, and thus we perform our analyses in the "cut-and-count" scheme. The significance of a given number of signal events N sig around the resonance peak, against a background N bkg , is defined as where α sys are the systematic and theoretical uncertainties on the SM rates. For definiteness, in what follows we always set α sys = 2%. As we will show, all our results are dominated by statistics up to systematic errors of 10% or larger. We refer to Appendix B for a precise assessment of the impact of different choices for α sys .
Before entering into the details of the analysis, to set a reference for the sensitivities, we compute the best possible reach that one would achieve in the case of negligible background. We define it as the signal cross section that results in 3 signal events  where L is the integrated luminosity. Using Eq. (16), this limit translates into an approximate sensitivity on the mixing angle Notice the logarithmic dependence on the particle mass for m 2 W m 2 φ s, explaining why our sensitivities are almost flat when compared with those obtained at hadron colliders. The aim of the following two sections is to determine how much a realistic analysis can approach the sensitivity in Eq. (20).
We now discuss the reach at different center-of-mass energies in the dominant decay channels hh, ZZ, and W W . As we show below, the sensitivities from the hh(4b) decay mode turn out to be very strong at lepton colliders. For this reason we start performing a detailed simulation of this channel, while we simply work at parton level (before showering) for the leptonic and semi-leptonic V V decays.

Decay channel φ → hh
In the model under consideration the largest individual branching fraction of the singlet is φ → hh(4b).
We look for this signal as a narrow resonant contribution over the SM background in the 4b invariant mass distribution. The same signature has been studied in [22], where the authors discuss the reach of ILC and CLIC 1.5 TeV. With respect to that work we include a full CLIC detector simulation.
Requiring W -fusion production, the principal background is the irreducible SM contribution to e + e − → 2ν4b, with a dominant component due to hh(4b) and Zh(4b). The total cross-section for this process is computed with MadGraph to be 1.8 fb (0.6 fb) at the center-of-mass energy of 3 TeV (1.5 TeV). A potentially large reducible contribution from γγ → 4b is avoided imposing cuts on the transverse momentum of the b quarks (p T > 20 GeV) and on the missing energy (E miss > 30 GeV), and turns out to be completely negligible.
We also compute the cross-sections for the signal e + e − → φ(4b)νν with MadGraph, after implementing the Lagrangian in Eq. (2) in FeynRules 2.0 [42], always working in the narrow-width Table 1. Efficiencies for signal and background in e + e − → 4b 2ν, for each individual cut applied in the analysis. The two cases m φ = 1 TeV and m φ = 2 TeV are shown, respectively, for CLIC Stage II and Stage III.
approximation for the singlet, and retaining the subdominant contribution from φ → ZZ. We use Pythia8 [43] for showering and Delphes3 [44] for detector simulation, using the configuration of the CLIC cards of Ref. [45]. We apply the VLC exclusive jet reconstruction algorithm [46] with working point R = 0.7 and N = 4 (see also Ref. [47]): this allows us to reconstruct b-jets with ∆R as small as about 0.1, well below the standard isolation cut, compatibly with the detector resolution expected at CLIC (see Appendix A for more details).
In order to select the events we proceed with the following steps: 1. We impose a cut on the transverse momentum of the jets p T > 20 GeV and on the missing energy E miss > 30 GeV in order to select events coming from W -fusion.
2. b-tagging: we require the presence of four jets tagged as b, using the loose selection criterion as implemented in Ref. [45] in order not to excessively reduce the signal efficiency.
3. h reconstruction: we identify the candidate Higgs bosons by choosing the pairing of the four bjets that gives reconstructed invariant masses of the two Higgses closest to 125 GeV, i.e. the one that minimises the quantity ( We then retain the events having two distinct b-pairs with m bb in a window of about [90,130] GeV. The exact boundaries of the invariant-mass window are chosen differently for each m φ hypothesis, in order to maximise the significance of the signal. 4. We apply a cut on the polar angle | cos θ| 0.9 of the two Higgs bosons, in order to reduce the contribution from the forward region, where the background is enhanced. The precise value of the cut is chosen for each value of the mass in order to maximise the significance.
5. φ reconstruction: we select the events with a total invariant mass of the 4b system in a window of about 0.75 m φ m 4b 1.05 m φ around the resonance peak, again optimising the cut for each signal hypothesis.  Table 1 for two benchmark cases. We verified that these numbers do not vary substantially changing the R parameter of the jet reconstruction algorithm, and changing the exact values of the kinematical cuts. For the signal, the most important effects come from b-tagging and from the Higgs mass reconstruction, which both have efficiencies in the 40%-60% range, depending on the resonance mass and on the collider energy, 3 for a total signal reconstruction efficiency of about 25%. The SM background, on the other hand, is reduced by a factor of at least a few 10 −3 , reaching up to 10 −4 for masses above a TeV at CLIC Stage III. In Figure 3 (right) we show the 4b invariant-mass distribution for the background at CLIC Stage III, with two signal examples superimposed. Notice that for singlet masses above a TeV the search becomes essentially background-free, and, after taking into account the efficiency of the signal, the limit roughly corresponds to the estimate in Eq. (19). The exclusion limits are computed from Eq. (18) requiring a significance of 1.64, which corresponds to 95% C.L. (one sided). In Figure 4 we show the results for the CLIC sensitivity in σ(e + e − → φνν) × BR(φ → hh) as a function of the mass of the singlet, and compare it with the reach in the various φ → V V channels described in Section 3.2. One can see that, despite the rather strong assumptions made to derive the other limits (see next section), the 4b channel turns out to be the best probe of scalar singlet production, if one assumes similar branching ratios into hh and ZZ. It is also evident that, due to the low number of background events over a large range of masses, the reaches depend only mildly on the collider energy and resonance mass as expected.

Decay channels φ → V V
Given the SM-like properties of the high-mass state φ, a large fraction of singlets is expected to decay into pairs of vector bosons. We consider four kinds of resonant signals: all of which allow a rather precise reconstruction of the resonance mass. Because of the excellent resolution expected at future lepton colliders, we do not perform a complete detector simulation for these channels, but we rather give an estimate of the reach obtained from the backgrounds calculated at parton level. The limits derived here thus do not stand on the same grounds as the ones of the previous section, and should be understood as optimistic estimates. Considering again W -fusion production mode, the dominant backgrounds to the e + e − → ννφ(V V ) signal originate from electroweak di-boson production. The cross-section rates for these processes, computed with madgraph, are while the corresponding values at √ s = 1.5 TeV are 18 fb and 52 fb. In the ZZ(4 ) channel, the other main source of background is e + e − → e − νW + Z (and its conjugate), in particular from tri-boson production, which however is reduced to 0.3 fb requiring the mass of the leptons to be within 5% of m Z . The same process, with a hadronically decaying Z, contributes also to ZZ(2 2ν). A few more words deserve to be spent for the hadronic modes. Given the broad shape of the hadronically decaying vectors, a non-negligible misidentification probability for Z and W should be taken into account in a realistic analysis. As a consequence, a contamination between the various V V backgrounds could become important. The dominant effect comes from the process e − γ → νW − Z and its conjugate, which is a background for both the 2 2j and the 4j channels and has a large cross-section of 163 fb at 3 TeV. This contribution could however be suppressed by requiring sufficiently hard momenta for the vectors in the final state (see e.g. Ref. [33]) and, for the ZZ signal, by a sufficiently hard cut on the dijet invariant mass. We leave a detailed study of this background to a future work and, in the spirit of giving an aggressive estimate of the sensitivities, we ignore it in what follows. The contamination between the W W and ZZ channels, on the other hand, requires a double misidentification, and is expected to be small. Therefore, in the following we assume that all W (jj) and Z(jj) will be told apart thanks to the excellent jet mass resolution at CLIC. Finally, we also assume that all backgrounds without neutrinos in the final state will become negligible after a suitable missing energy cut is imposed, and consider only the irreducible backgrounds of Eq. (21) for our analysis.
For each signal mass m φ and final state f , we select the simulated background events with an invariant mass of the two vectors that falls within a window of width ∆ f around the resonance peak, We take ∆ 4 = 5%, ∆ 2 2j = 10%, and ∆ 4j = 15% for the three channels under consideration. The first two numbers are in rough agreement with the resolutions that are achieved at the LHC for resonances reconstructed in ZZ [49], while the third number correctly reproduces our results for hh(4b) of the previous section. The expected sensitivities are then computed from the resulting number of background events, solving Eq. (18) to find the excluded number of signal events at 95% C.L. (one sided), rescaled to take into account the branching ratios of W and Z into the relevant final states. We show the sensitivities obtained this way in Figure 4, for the exclusive channels ZZ(4 ) and ZZ(2 2j), and for the combination of the channels ZZ(4 ), ZZ(2 2j), ZZ(4j), and W W (4j).

Discussion and comparison with hadron colliders
We now translate the projected sensitivities on the cross-section into a limit on the mixing angle sin 2 γ, which is trivially obtained rescaling by the cross-section of a SM Higgs with the same mass times the singlet branching ratio into hh or V V . For CLIC we compute the Higgs production cross-section using Madgraph 5 at LO, see Figure 2 left.
Dependence on BR φ→hh of the combined direct reach in hh or V V at CLIC, compared with the sensitivity to sin γ expected from Higgs couplings measurements taken from Ref. [48]. The green and blue lines are for m φ = 500 GeV at CLIC Stage II and III, respectively, the gray dashed line is for m φ = 2 TeV (at 3 TeV).
The pink lines show the precision in Higgs couplings at the three CLIC stages.
In Figure 5 we show the combined reach from the hh(4b) and V V searches described in the previous section, as a function of BR φ→hh and compare it with the expected reach in Higgs signal strengths at the various stages of CLIC. Notice that, especially for lower singlet masses, direct searches are more powerful than indirect ones at each stage. Varying the value of the BR φ→hh , the searches in the two channels hh and V V are clearly complementary.
We then compare our CLIC results with Higgs couplings measurements and direct searches at various stages of the LHC: the current run at 13 TeV, the end of the 14 TeV run with 300 fb −1 , and the high-luminosity phase with 3 ab −1 . For what concerns the Higgs couplings, we display as excluded the combined ATLAS and CMS 8 TeV constraint [50], in order to be conservative. Indeed, the 13 TeV best-fit Higgs signal strength from CMS [51] is larger than the SM one by almost two sigma. For direct searches we show the present LHC exclusions [52,53] as well as the projected sensitivities at 300 fb −1 and 3 ab −1 .
To determine future sensitivities at pp colliders, we rescale the expected sensitivity of the existing 13 TeV search [52] at higher energies and luminosities using quark parton luminosities, with a procedure analogous to the one presented in Ref. [32]. We refer the reader to Appendix B for more details on this aspect, as well as for the expected sensitivities also at the HE-LHC with center-of-mass energy of 27 TeV, and FCC-hh at 100 TeV, that will be used in the next section when comparing to muon colliders.
To translate the future sensitivities on cross-section to a reach in mixing angle, we use the SM Higgs production cross-section at hadron colliders. We include gluon fusion plus VBF plus VH production, and compute the first one with ggHiggs v3.5 [54][55][56], and the second two with Madgraph 5 at LO.
It can be seen that CLIC at 1.5 TeV and 3 TeV will be able to probe singlet masses up to about 1.2 TeV and 2.6 TeV, respectively. CLIC at 1.5 TeV stage will be more sensitive than the high-luminosity LHC up to masses of about 1 TeV, while the 3 TeV stage will be significantly more sensitive over the full mass range.  Figure 6. Exclusions at 95% C.L. in the plane (m φ , sin 2 γ). The shaded regions are the present constraints from LHC direct searches for φ → ZZ (red) and Higgs couplings measurements (pink). The reach at CLIC Stage II (green) and Stage III (blue) in φ → hh(4b) is compared with the projections for LHC in φ → ZZ with a luminosity of 300 fb −1 (solid red) and 3 ab −1 (dashed red). We have fixed BR φ→hh = BR φ→ZZ = 25%. The dashed grey lines show two different scalings of s γ with m φ , as described in Section 2.1 (g * = 1 in both cases).

Muon colliders
Among the proposals for muon colliders, the two main categories differ in the way muon bunches are produced: either from protons on target (MAP [9]), or from positrons on target (LEMMA [10][11][12]). The former case guarantees a larger instantaneous luminosity, but is complicated by the technological challenge of muon cooling. The latter case overcomes this difficulty at the price of a slightly reduced luminosity. See Ref.
[57] for an updated discussion about the state-of-art of muon colliders.
Here and in the following we simply focus on an idealised muon collider at 6 TeV and 14 TeV, respectively with 6 ab −1 and 14 ab −1 of integrated luminosity, as benchmarks attainable both by MAP and LEMMA. We determine their sensitivity to resonances decaying to hh(4b) only at the Madgraph level, because of the present lack of knowledge of the detectors that will be used at those machines.
We simulate the background processes µ + µ − → ννhh and impose the cut |η h | < 2 on the pseudorapidity of each Higgs boson, which roughly corresponds to | cos θ h | < 0.95. For every signal mass m φ that we want to test, we then take the fraction of background events that satisfies m hh = m φ ± 15%. Finally we assume an additional efficiency of 30% (as a rough estimate of b identification and other effects), and we determine the 95%C.L. sensitivities according to Eq. (18), using as usual systematics of 2%, and of course taking into account the branching ratio of h → bb. For the signals, we compute the total + − → ννφ cross-sections at LO with Madgraph, at all the machines of our interest. We then just impose the same efficiency of 30%, assuming it captures the effects of the various cuts (this assumption is to some extent supported by the study in Section 3.1).
We find that this procedure reproduces extremely well, at both 1.5 and 3 TeV, the results of the more careful detector study of Section 3.1, at least for m hh 700 GeV. We report in Appendix B more details on the validation above, as well as the sensitivities both on the mixing angle and on the production cross-section of a generic resonance decaying to hh, at lepton machines from 1.5 TeV to 14 TeV of center-of-mass energy. Since these searches are essentially background-free for large masses, they are dominated by statistical errors. We discuss the impact of systematic errors in more detail in Appendix B, also in relation with possible target luminosities for muon colliders.
Here, we show in Figure 7 the 95% C.L. sensitivities in the plane (m φ , sin 2 γ) at √ s = 6 TeV and 14 TeV, for total integrated luminosities of 6 ab −1 and 14 ab −1 , respectively. We also compare the reach of muon colliders to the one of high-energy hadron collider proposals such as HE-LHC and FCC-hh. The take-home message of this comparison is that HELCs in the very high energy regime could become very powerful discovery machines, even stronger than future hadronic colliders, at least for New Physics mostly coupled to the Higgs sector.

Single Production & Beyond the Standard Model Scenarios
In this section we discuss the implication of the CLIC reach on singlet resonances in well motivated Beyond the Standard Model (BSM) scenarios.

NMSSM
In the NMSSM, the particle content of the MSSM is extended with a singlet of the SM gauge group S, so that the superpotential reads W = W MSSM + λ SH u H d + f (S), with f a polynomial up to degree 3. The SM-like Higgs boson mass receives an extra tree-level contribution, which lifts its upper limit to where t β = tan β is the ratio between the up and down Higgs VEVs, and ∆ hh encodes the usual SUSY radiative contributions. In light of the null LHC searches for coloured superpartners, the NMSSM with a large coupling λ is a particularly attractive SUSY model from the point of view of naturalness of the EW scale, see e.g. Ref. [58][59][60][61]. Indeed, the fine-tuning needed to reproduce the EW scale is parametrically alleviated, for a fixed value of the stop and gluino masses and with respect to the MSSM, by a factor ∼ λ/g. Having λ 1.5 would however overshoot the Higgs mass, and thus introduces a new tuning problem to bring m h down to its measured value. In addition, λ 2 becomes strongly coupled at scales of O(10) TeV. 4 The NMSSM constitutes an ideal application of our sensitivity study for extra singlet scalars. To describe the scalar sector of the NMSSM we employ the parametrisation put forward in Ref. [63,64]. Then, to obtain a simple description of the scalar singlet phenomenology, we assume the extra Higgs doublet to be in the decoupled or alignment limit. This is also somehow favoured by the present LHC constraints, because the mixing of the SM Higgs with a doublet is more constrained than the mixing with a singlet, see e.g. Ref. [63]. The phenomenology of the SM-like Higgs plus the extra singlet can then be described by 4 free parameters To connect ∆ hh with parameters of more immediate physical interpretation, we employ the concise analytical expression (see e.g. Ref. [65,66]) where mt = √ mt 1 mt 2 , mt 1,2 are the physical stop masses, X t = A t − µ/t β , M t = 173 GeV is the top pole mass, m t (Q) is the running top mass, and Q t = M t mt. Such expression is accurate to the level of a few GeV in ∆ hh , which is more than enough for our purposes. The phenomenology of the Higgs plus singlet system is displayed, for the NMSSM, in Figure 8, where we fix λ = 1 as a benchmark motivated by naturalness. In the left-hand panel we let m φ and t β vary and we fix ∆ hh = 80 GeV, a value obtainable e.g. for stop masses in the range of 1-2 TeV. The precise value of ∆ hh does not affect the Higgs sector phenomenology as long as it is within O(10%) of 80 GeV. One sees that direct searches for the extra singlet, at both CLIC stages II and III, would probe a parameter space that is completely unexplored by the HL-LHC. At CLIC, direct φ searches and Higgs coupling measurements would constitute a complementary probe of the parameter space, as noticed already in Section 3.
To connect with the phenomenology of the SUSY coloured sector, in the right-hand panel we let m φ and mt vary, for t β = 5 and X t = mt. We see that the precise value of the stop masses does not have a major impact on the phenomenology of the Higgs-singlet sector. This consideration holds independently of the values t β and X t , with the only exception of large t β and very small X t (in which one recovers the MSSM problems in reproducing the correct Higgs mass). We can therefore conclude that searches for the singlet scalar will not provide significant information about the stop masses. They will, on the other hand, give a measure of the tuning in the scalar sector, which is a dominant source of tuning in the NMSSM. The contribution due to the singlet can be roughly quantified as λ 2 s 2 /m 2 Z ∝ m 2 φ /m 2 Z (where s is the VEV of the singlet). 5 Here we have not discussed deviations in the trilinear Higgs coupling in the NMSSM. They depend on more parameters than those in Eq. (23), and can reach O(50%) or more if λ 1, see Ref. [67] for a precise quantification.

Twin Higgs
In Twin Higgs models the SM Higgs sector is extended by adding the twin Higgs H B , which is a singlet under the SM and a doublet under a mirror EW gauge group SU (2) B . The twin Higgs is coupled to the SM Higgs H A via a portal coupling λ * , that realises a global SO(8) symmetry at tree-level. This is spontaneously broken to SO (7) at a scale f . The radiative stability of the construction is ensured by an approximate Z 2 symmetry between the SM and the mirror sector, so that the full content of the SM, or part of it, is doubled [68]. An explicit breaking of the Z 2 symmetry is then introduced to allow for f > v and for a viable phenomenology [24]. All in all, the scalar potential reads where κ is the SO(8)-breaking quartic and we parametrise the Z 2 -breaking contributions by σ soft and ρ hard . These correspond to a soft and a hard breaking of the discrete symmetry respectively. Notice that κ receives irreducible IR contributions from (mirror) top-loops, and that if Z 2 is broken only softly at some scale, then a small quartic ρ hard is generated at lower energies by Z 2 -breaking 1-loop thresholds.
After the spontaneous breaking of the SO(8) symmetry, we are left with two real scalars in the spectrum: the SM-like Higgs h and the radial mode φ. Their physical masses, in the limit λ * κ , σ soft , ρ hard , read  This explicitly shows that the Higgs mass is of the correct size for typical values of the parameters κ and ρ hard . In TH models the fine-tuning is parametrically reduced with respect to the ones of regular SUSY or Composite Higgs scenarios by λ H /λ * , where λ H 0.13, see e.g. Ref. [69,70]. In models where the Z 2 -breaking is mostly achieved by the quartic ρ hard , one obtains an additional gain in fine-tuning by λ H /|λ H − ρ hard |, which is maximised for ρ hard as close as possible to the SM quartic λ H . This parametric gain is however limited by the irreducible IR contributions to the SO(8)-breaking quartic κ, as discussed in Ref. [69].
The requirement to reproduce the EW scale v and the Higgs mass m h fixes 2 out of the 5 free parameters in Eq. (25). We choose the three remaining free parameters as the spontaneous breaking scale f , the physical singlet mass m φ , and the Z 2 -breaking quartic ρ hard . We then find the following exact analytical expression for the mixing angle which to our knowledge was never presented before in the literature. Other useful relations, to track the impact of ρ hard , are where g SM hhh = 3m 2 h /2v, and where we ignored all higher orders in v 2 /f 2 and λ H /λ * ≈ ρ hard /λ * . Eqs (27)- (28) show nicely that the effect of the new quartic decouples with the mass of the singlet state (or equivalently with λ * ), and therefore it could affect the phenomenology of the scalar sector only at small to intermediate m φ .
The parameter space of TH models is displayed in Figure 9 in the m φ − f plane, for two benchmark values of the hard breaking quartic ρ hard . 6 As anticipated by the analytical understanding above, the region where the impact of ρ hard = 0 is most visible is the one where m φ is small. In particular we see that a non-zero ρ hard allows the Higgs mass constraint to be satisfied at large f and small m φ . In this region the Higgs mass is mostly achieved via ρ hard . However, in the same region the fine tuning gain of the TH is limited because λ * 0.1 [69]. Figure 9 also displays the phenomenological results of Section 3, where we have extended the framework to include the invisible decays of the radial mode into W W , Z Z (all with masses m W ×f /v, because the U (1) could well be not gauged [72,73]) and t t (with mass m t ×f /v). The SO(8) symmetry implies that the invisible branching ratio asymptotises to 3/7 for m φ m t . One learns from Figure 9 that the phenomenology of the twin Higgs φ is independent on how the Z 2 -breaking is achieved, at least in the region of parameter space where the fine-tuning is ameliorated. HELCs like CLIC are expected to probe the most natural regions of TH models mainly via their precision in Higgs coupling measurements. While direct searches for the radial mode would constitute a weaker probe of the interesting region of the parameter space, they could provide precious complementary information. A similar conclusion was drawn also in [22], where the hh(4b) signature was studied.

Comments on heavy electroweak ALPs
Our results can be applied generically also to scalar resonances that are produced singly from the fusion of transverse W bosons. Resonances of this type are the so-called axion-like particles (ALPs), a quite generic category of pseudo-scalar particles coupled via ABJ anomalies to the SM gauge bosons. These arises in many theoretical models related to Dark Matter production [74], Naturalness [75][76][77] and vector-like confinement [78].
In this context we consider a somehow heavy ALP a with only electroweak anomalies and mass m a > 2m W . The effective Lagrangian for an ALP of this type reads whereF µν = (1/2) µνρσ F ρσ for any field strength, α 1 = 5α Y /3 is the GUT-normalised U (1) Y coupling constant. The scale f a is associated to the spontaneous breaking of the U (1) under which the ALP shifts, while we do not specify the origin of the explicit breaking terms m a . The above theory has a physical cut-off at a scale g * f a , where both the radial mode and new particles charged under the electroweak group appear, so that for consistency m a < g * f a . The anomaly coefficients c 1,2 depend linearly on the number of states in the UV physics which carry EW quantum numbers and are charged under the global U (1). After electroweak symmetry breaking (EWSB), one can write Independently of the underlying UV theory, the anomaly coefficients are correlated with g * . For example, in a QCD-like sector with N colours g * ∼ 4π/ √ N , and for a 'composite' ALP from that sector one roughly expects c 1,2 ∼ N from the degrees of freedom above the confinement scale. Similarly, in weakly coupled models with N f charged fermions of masses ∼ g * f a one has c 1,2 ∼ N f (see e.g. Ref. [77]), and g * 4π/ √ N from perturbativity. The arguments above indicate that a large ALP coupling to gauge bosons requires a somehow small g * , and therefore it strengthens the bound m a < g * f a .
In order to focus on a concrete model, we choose as a benchmark the photophobic ALP discussed in Ref. [27], defined by c γγ = 0. Notice that assuming a zero photon coupling is a good approximation even after radiative corrections are included. Analogous results can be obtained in more general models as long as |c 2 /c 1 | is large enough to not make the branching ratios in vector bosons subdominant compared to the ones in final states containing photons. In the photophobic ALP the WW branching ratio dominates over the ZZ with BR W W /BR ZZ 4 at high masses.
For this reason in the right panel of Figure 10 we show the reach in the plane (m a , f a ) of the W W channel at the stage II and III of CLIC as determined by the analysis in Section 3.2. The only difference compared to the cases discussed in previous sections is that the ALP couples only to the transverse polarisations of the gauge bosons. The cross sections for single production at 1.5 TeV and 3 TeV are given in the left panel of Fig. 10. These are computed for the signal e + e − → a(2V )νν with MadGraph, after implementing the Lagrangian in Eq. (2) in FeynRules 2.0.We refer to [38] for a derivation of analogous formulas to the ones presented in Eq. (16) for the transverse polarisations.
The final reach depends largely on c 2 = 16π 2 /g 2 * , which affects the rate at quadratic level, and it has been set to reproduce a benchmark with g * = 4. A relatively small g * is required to have a sizeable production rate. However, in the same regime, the EW charged UV states responsible for generating the anomaly are pushed down to be within the reach of high-energy colliders. Aware of this issue, we show the line corresponding to masses of new electroweak states of 0.5 and 1 TeV, which roughly indicates the reach of CLIC on the pair production of these particles at stage II and III respectively (see e.g. Ref. [79]). Notice that this is somehow less transparent in the parametrisations of other studies (e.g. Ref. [25][26][27]), where the loop factor in the anomaly coefficients in Eq. (29) is reabsorbed in the definition of the cut-off scale.

Double production & first order Electroweak Phase Transition
When sin 2 γ goes below 10 −3 single production becomes ineffective at CLIC. In this case, which is possibly related to an underlying Z 2 symmetry acting on the singlet, the only sizeable process can be the production of singlet pairs. This channel can give sizeable signals when the portal coupling λ HS in 0.12 0.03 0.018 7 × 10 −3 2.1 × 10 −3 Table 2. A selection of branching ratios leading to interesting final states double singlet production, assuming BR(φ → W W ) = 2 BR(φ → ZZ) = 2 BR(φ → hh) = 50%.
Eq. (2) is non negligible, as discussed in Section 2. In this section we consider the following production processes e + e − → φφνν , which are the dominant ones at HELCs and can be relevant for singlet masses below 1 TeV at the Stage III of CLIC. The numerical value of the double production cross-section for sin γ = 0 is shown in the right panel of Figure 2. In order to understand the possible sensitivities of the channels in Eq. (31) we first try to compare with the situation at the LHC, where related searches have been suggested. At LHC the VBF production can also benefit from sizeable couplings, however the gluon-fusion process is available and it dominates the total rate for light masses. This production channel has been exploited in a similar context in Ref. [31,80]. 7 The comparison between the rate at the 14 TeV LHC and 3 TeV CLIC is plotted in the left panel of Figure 11, where the rates are normalised to λ HS = 1. The typically smaller backgrounds at HELCs are likely to overcome the smallness of the rates for light masses. This suggests that lepton colliders can be an extremely useful tool to study pair production of scalar singlets that modify the Higgs potential.
Depending on the possible final states, different experimental strategies can be undertaken. A table with possible final states for a promptly decaying singlet and the relative branching fractions is given in Table 2. Given the relatively small backgrounds for visible final states with multi jets and/or leptons, we could expect a large improvement when a combination of several channels is considered. For this reason in the left panel of Figure 11 we show a line of 100 signal events at 3 ab −1 for λ HS = 1 in order to show the possible reach in mass if the singlets could be fully reconstructed. This corresponds to roughly 10 events in four vectors decaying hadronically.
Broadly speaking, we can classify the final states in three main categories depending on the values of sin γ, which controls the decay length of the scalar singlet. This gives rise to: i) prompt decays; ii) displaced decays; iii) collider-stable singlets. The right plot in Figure 11 shows the various categories as a function of sin γ. We now discuss the opportunities at HELC for these final states and the comparison with the analogous search strategies at the LHC.
i) Prompt decays are obtained in a large region of the singlet parameter space for which single production is subdominant. This corresponds to mixing angles between 10 −2 and 10 −6 -10 −8 . The latter numbers correspond to require a decay length of less than 1 mm for a singlet mass of 100 GeV and 1 TeV, respectively. In this portion of the parameter space multi-boson and multi-Higgs final states can give spectacular signals, with very low or even negligible backgrounds. 7 The total production cross-section can be computed by integrating over the gluon parton luminosities dLgg/dτ as [81] σ(pp → φφ) where the loop function is ii) Displaced decays with displacement between 0.1 cm and 100 m are obtained for mixings down to sin γ ≈ 10 −8 − 10 −11 . This would lead to spectacular signals with multiple displaced tracks in the tracker (between 0.1 and 100 cm) or in the muon chamber (between 1 and few tens of meters). This type of events has been shown to be basically background free at the current stage of the LHC (see e.g. Ref. [82,83] and the sensitivity derived for displaced decays of scalar pairs in Ref. [84]). At HELC like CLIC these channels are going to be cleaner than at the LHC because the reduction in multi-jets backgrounds is going to diminish the expected rate of fake displaced vertices.
iii) Invisible decays correspond to mixings smaller than 10 −11 . This is the so-called "nightmare scenario" where the scalar singlet can only be pair produced and decay invisibly [85]. A comprehensive study on the possible reach at the LHC has been performed in Ref. [31]. This shows that searches at the HL-LHC with √ s = 14 TeV and 3 ab −1 can probe λ HS ∼ 6 for m φ ∼ 200 GeV by combining double production channels for the singlet in VBF, gluon fusion and associated production with tt. At CLIC the dominant background for an invisibly decaying singlet pair produced in Z-boson fusion is coming mostly from σ(e + e − → e + ν e W − ) = 0.48 pb with the W decaying leptonically. By using our analytical estimate of the signal in Eq. (17), we can extrapolate at CLIC the results obtained in Ref. [86] for ILC at 1 TeV and 1000 fb −1 . For a 100 GeV singlet we get S √ B 2.3 λ 2 HS at the stage III of CLIC, by rescaling the sensitivity in Ref. [86] with the square root of the luminosity times signal cross-section, thus assuming that the background cross-sections would scale in energy roughly as the signal ones. This gives hope to cover this type of scenarios at CLIC. A more careful study for final states with missing transverse energy is left for future work.
In the next section we discuss how double production at HELCs can possibly probe a interesting region of the parameter space where a first order electroweak phase transition is induced by the singlet and the deviations in the Higgs couplings are below the best sensitivity of future lepton colliders. Figure 12. Test of the double singlet production at CLIC. The region with a first order electroweak phase transition is enclosed in the shaded regions, as discussed in the text. Isolines of 100 events refer to e + e − → φφνν at CLIC 1.5 TeV (red) and CLIC 3 TeV (blue). Contour lines of deviations in the trilinear Higgs coupling (dashed, see Section 2) and in Zh production (dotted, see [87,88]) correspond to possible CLIC sensitivities [19].

Electroweak phase transition
It is well known that the presence of a singlet with a sizeable coupling to the Higgs field can induce a first order electroweak phase transition (FOEWPT) (see Ref. [89] for a review and references). It is therefore natural to ask whether our analysis can constrain such scenarios thanks to the pair production of singlets. Many detailed studies of the FOEWPT have been performed in the literature [85,[90][91][92][93], and numerical codes are available [94]. Here, we are willing to sacrifice numerical accuracy for analytic simplicity and we will adopt a as much as possible simplified description in order to characterise the main regions where a FOEWPT can occur. When the SM plasma has a sufficiently high temperature T , the electroweak symmetry is restored due to thermal correction to the Higgs mass. Once the temperature drops as an effect of the expansion of the Universe, tunnelling to the true vacuum might happen. In what follows we describe the main features of the two regions where a FOEWPT can occur in our model, displayed in Figure 12. To simplify our discussion, in this section we assume that the Z 2 symmetry forbidding the a S and a HS terms in the Lagrangian Eq. (2) is a good approximate symmetry, up to a small perturbation parametrised by the the mixing sin γ, following Ref. [30].
• Two-step Phase Transition. The so-called 2-step phase transition occurs when the singlet develops a VEV at finite temperature before the Higgs settles to its minimum. Later on (at smaller temperatures), this vacuum tunnels via a strong first order PT to the vacuum with the present EW VEV, which is the global minimum for T = 0. The situation can be arranged with and λ S large enough, to guarantee that the global EWSB minimum has S = 0 while another local minimum appears at h = 0 and S = 0. This condition gives a lower bound on λ S : When the above inequality is saturated the two minima are actually degenerate and the critical temperature approaches zero. All in all the region in the plane (m φ , λ HS ) where the two-step phase transition can occur is roughly given by The requirement that the singlet develops a VEV before tunnelling at finite temperature to the vacuum with S = 0 would actually impose a more stringent condition, the determination of which goes beyond the purpose of this paper. We choose to use the region defined by m 2 S < 0, as it conservatively contains the region satisfying the more precise requirement. The two different green regions in Figure 12 correspond to different perturbativity bounds on the couplings λ S and λ HS : in the light green region we require β λ HS /λ HS < 1, while in the darker green region we also require β λ S /λ S < 1. β λ HS and β λ S are the beta functions of the corresponding couplings 8 (see Ref. [95] for a discussion of the connection of this bound with the more standard requirement of not having Landau poles at the TeV scale).
• One-step Phase Transition. As already discussed in Section 2, in the Z 2 -symmetric limit we can integrate out the heavy singlet generating the dim 6 effective operators in Eq. (7). In this scenario the quartic of the singlet does not play any role and can be set to zero for simplicity. A first order phase transition can occur in the Higgs EFT if the effective Higgs quartic is negative and the potential is stabilised by the dimension 6 operator |H| 6 . In the Z 2 -symmetric case, requiring the loop corrections to the Higgs quartic to be big enough to make it negative gives roughly the region displayed in Figure 12. The precise shape of the region requires to include the full Coleman-Weinberg 1-loop contributions at zero temperature and the finite temperature corrections. We take this result from Ref. [30] which agrees sufficiently well with the computation that includes resummation of thermal loops [96]. The perturbativity constraint requires β λ HS /λ HS < 1 as above. The singlet self-coupling is irrelevant in this scenario and can be set to zero. Notice that going beyond the Z 2 symmetric case, the first term of c 6 in Eq. (7) can be made large with a sizeable mixing angle. This scenario has been shown to induce a FOEWPT in Ref. [97].
The main result of this section is Figure 12, where we show the region where a FOEWPT occurs based on the discussion above. Our approach has two main limitations. First, we did not compute the temperature where bubble nucleation occurs, and this could further shrink the 2-step region as emphasised in Ref. [98]. Second, we did not impose the condition for a fast decoupling of the sphaleron transitions inside the bubbles, v c /T c 1, which is a necessary requirement for EW baryogenesis. Such region, however, has been found to almost coincide with the blue shaded area in our Figure 12, see e.g. Ref. [30]. Therefore, we expect the region where EW baryogenesis may take place in this model to be fully contained in our shaded areas, at least for not too large values of the mixing angle sin γ. Figure 12 shows the relevance of pair production as a test of models with a FOEWPT. For this purpose we plot isolines of 100 number of events at 1.5 TeV and 3 TeV CLIC. They may be enough to test this model under the reasonable assumption that many of the multi-Higgs and multi-bosons final states would face very small to zero backgrounds (e.g. 8b), and that one could combine different channels as suggested by Table 2. As shown in Figure 11, lowering further the mixing angle would improve the reach of double production in the displaced region until the singlet would become long lived on collider scales. A dedicated analysis would be required to precisely assess the reach of CLIC for invisible final state. We also show the projected sensitivity of CLIC at 1.5 TeV and 3 TeV on triple Higgs coupling deviations (taken from Ref. [19]), and two benchmarks for deviations in Zh associated production. These deviations are generically predicted in this setup as shown from Eq. (6).
In conclusion, pair singlet production at HELCs has the potential to test the entire parameter space allowing for FOEWPT and potentially EW baryogenesis. It constitutes a complementary probe to deviations in triple Higgs couplings and in Zh associated production (see e.g. Ref. [19,99]).

Outlook and conclusions
A clean background environment and a high energy in the center of mass constitute a dream for any particle physicist. Machines satisfying both properties are High Energy Lepton Colliders (HELCs), like CLIC [7] and more futuristic proposals of muon colliders [8,9,11,100]. In this paper we made progresses in building their physics case.
HELCs allow for powerful precision tests of the Higgs and electroweak sectors, as well as for direct production of new resonances beyond the reach of current experiments. As discussed in the Introduction, the SM Higgs boson is dominantly produced in W W -fusion, which may allow a precision on the Higgs coupling at the per-mille level. Building upon this observation, we studied direct W Wproduction of new scalar resonances which are singlet under the SM gauge group and couple to the SM through the Higgs sector. These particles are very well motivated from the theoretical view point and represent an important target for future collider machines.
The first message of our study is summarised in Figure 6. HELCs like CLIC extend the reach on the resonant production of a singlet decaying into di-bosons and/or di-Higgs well beyond the HL-LHC reach, i.e. down to couplings which correspond to Higgs coupling deviations at the per-mill level and up to 1-2 TeV masses. A similar conclusion holds for more ambitious proposals like future µ colliders (like LEMMA or MAP), which stand as fantastic discovery machines to probe new physics coupled to the Higgs sector (see Figure 7 for a comparison with future hadron colliders). The consequences of the improved reach in single production are deep on the theoretical side, as we explored in Section 4. Singlet searches and Higgs coupling deviations might represent the only window to explore models of Neutral Naturalness like the Twin Higgs, where the coloured states can be beyond the LHC reach. In conventional models of Naturalness a new singlet might have other reason to be within the reach of CLIC, for example the SM Higgs mass in the NMSSM. The same searches can test heavy axion-like particles coupled to EW gauge bosons. These arise from the spontaneous breaking of approximate global symmetries and could be e.g. related with Dark Matter production and/or vector-like confinement.
The second message of our study is shown in Figure 12, where the reach in double production of singlet scalar particle coupled through the Higgs portal is compared against the allowed parameter space of electroweak baryogenesis. Depending on the decay length of the singlet different kinds of final states can be probed at CLIC. Without attempting to perform a detailed collider study for all of them, we point out that CLIC at √ s = 3 TeV has the potential to explore models exhibiting a first order electro-weak phase transition, well beyond the reach of indirect constraints from Higgs coupling deviations. The same region can give correlated signals in future interferometers for gravitational waves such as LISA [101]. We believe that our results motivate a more detailed collider study of the various channels for singlet pair production, in order to assess robustly the reach of such mode.
We hope that this study represents another little piece of motivation to push forward the quest to explore the high energy frontier.   Figure 14. Sensitivities on the signal cross-sections of a new resonance decaying into ZZ, at the high-luminosity (dashed) and high-energy (dot-dashed) LHC, and at the FCC-hh (dotted). They are derived from the present LHC sensitivity at 13 TeV with 36 fb −1 of luminosity [52] (continuous red), see text for more details. For comparison, we also show the 13 TeV limit corresponding to the sensitivity of Ref. [52] (continuos black), and the sensitivity that one would obtain at 13 TeV with 36 fb −1 , rescaling the 8 TeV searches [103]. The latter serves as a rough validation of our procedure.
rescaling the expected sensitivity of the existing 13 TeV search [52] to higher energies and luminosities using quark parton luminosities, with a procedure analogous to the one presented in Ref. [32]. At the time Ref. [32] was published only 8 TeV LHC results were available, so the sensitivities presented there were obtained rescaling the 8 TeV searches. The results of Figure 14 are thus to be considered as an update of Ref. [32]. As a rough validation of the method, we find that the 13 TeV sensitivity determined rescaling the 8 TeV results [103] with the same procedure is in reasonable agreement with the one of the actual experimental search [52] (see the green dashed and continuous red lines in Figure 14).
Future lepton colliders. In Figure 15 we plot the sensitivities on the signal cross-sections of a singlet that decays to pairs of Higgs bosons, at CLIC Stage II (1.5 TeV, 1.5 ab −1 ) and Stage III (3 TeV, 3 ab −1 ), and at a µ-collider 6 (6 TeV, 6 ab −1 ) and µ-collider 14 (14 TeV, 14 ab −1 ). They are determined at the Madgraph level, from the simulation of the background + − → ννhh at parton level, as explained in Section 3.4. The comparison of the sensitivities determined in this way (dashed blue lines), with those determined from a proper study including detector simulation in Section 3.1 (continuous gray lines), gives an explicit visualisation of the agreement between the two. In order to better assess the robustness of our analysis, we also compare the relative importance of statistical and systematic uncertainties in setting our limits and sensitivities. On the left-hand side of Figure 16 we show the statistical error of the sensitivities at CLIC and at muon colliders, for the specific values of the luminosities that we have used. On the right-hand side of the same figure, we show the value of the luminosity needed in order for the systematic uncertainties to become important. For definiteness, we use a conservative benchmark value of α sys = 10%, recalling that the resulting luminosity scales as α −2 sys . These results show that resonant hh searches at CLIC will always be statistically dominated, and thus our results are independent of the precise value of α sys used in the analysis. Similarly, for the chosen benchmark luminosities, our estimated reaches at muon colliders are largely dominated by statistical errors for masse above a TeV, while they are expected to become sensitive to systematic errors of ∼ 15% (6 TeV) and ∼ 7% (14 TeV) for low masses.  Figure 15. Comparison between the determination of sensitivities using Madgraph only (both at CLIC and future muon colliders) and the full results for CLIC, as a validation of our simplified analysis. Left: Higgs-singlet mixing angle, assuming BR(φ → hh) = 25%. Right: signal cross-section of a generic resonance produced in e + e − and decaying to hh.  Figure 16. Left: Statistical uncertainty on the background expected in resonant hh searches at CLIC and muon colliders. Right: luminosity at which the impact of statistical uncertainties on our sensitivities becomes equal to the one due to systematics (above the continuous lines a systematic uncertainty of 10% is more important than the statistical one).