Collider searches of scalar singlets across lifetimes

Spin-0 singlets arise in well-motivated extensions of the Standard Model. Their lifetime determines the best search strategies at hadron and lepton colliders. To cover a large range of singlet decay lengths, we investigate bounds from Higgs decays into a pair of singlets, considering signatures of invisible decays, displaced and delayed jets, and coupling fits of untagged decays. We examine the generic scalar singlet and the relaxion, and derive a matching as well as qualitative differences between them. For each model, we discuss its natural parameter space and the searches probing it.


Introduction
Light spin-zero singlets are ubiquitous in models of New Physics (NP). They can have important phenomenological roles such as serving as a portal to a Dark Sector [1] and rendering the electroweak phase transition first order to enable electroweak baryogenesis [2,3]. In many cases, the phenomenology associated with such NP can be encompassed in the minimal renormalizable extension of the Standard Model (SM) obtained by adding one spin-zero singlet φ [4]. We consider this model as a benchmark, assuming all other new degrees of freedom are sufficiently heavy or weakly coupled to the SM particles.
Despite its simple setup, the singlet extension brings about a rich phenomenology related to the Higgs, by opening the exotic decay channel h → φφ, if kinematically allowed (see e.g. Ref. [5]), and by reducing the couplings of the Higgs boson to SM particles via singlet-Higgs mixing. This applies equally to scalars and pseudoscalars, though in the latter case the φ-Higgs mixing requires breaking of CP. The phenomenological implications reach far beyond Higgs-related observables, as the singlet inherits the couplings of the Higgs to the SM particles, suppressed by the mixing angle. Therefore, the interactions of the singlet can, depending on its mass m φ , be probed across the precision, luminosity and energy frontiers. The various signatures of the singlet include its effects on atomic physics, tests of the equivalence principle, Dark Matter (DM), beam dump experiments, rare meson decays, collider signatures as well as astrophysical and cosmological observables, see e.g. Refs. [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] and references therein. In this work, we focus on the collider searches in the mass range of 3 GeV ≤ m φ ≤ m h /2.
In addition to the above generic renormalizable extension, we consider the relaxion framework [23], which offers an alternative approach to the gauge hierarchy problem, and can also provide a DM candidate [24,25], facilitate baryogenesis [26], and address other hierarchies of the SM [27]. The relaxion is a naturally light pseudo-Nambu-Goldstone field, whose variation in the early Universe induces the scanning of the Higgs mass. The key ingredient of the relaxion mechanism is the so-called backreaction potential, which stops the relaxion evolution at the observed Higgs mass value. The backreaction potential is responsible for the interactions between the Higgs and the relaxion which are relevant for the collider phenomenology, realizing the Higgs portal structure similarly to the generic singlet extension discussed above [16,17]. The mass range examined in this work, however, represents only the extremely heavy region of the full relaxion parameter space.
The two main parameters in our focus are the Higgs-relaxion mixing and the h 2 φ 2 coupling. We demonstrate that, while having many similarities with the generic portal models, the relaxion is more constrained, but at the same time allows for larger values of the mixing angle than in the generic portal scenarios. This feature occurs because the naturalness considerations, which can be used to constrain the portal parameter space, can be automatically violated by the dynamics of relaxion field [28].
While collider constraints on promptly decaying relaxions or light scalar singlets were derived in Ref. [29], here we focus on the range of couplings that makes the scalar sufficiently long-lived such that it does not decay promptly. We take indirect constraints from global Higgs coupling fits as well as direct searches for invisible Higgs decays and displaced and delayed signatures into account. Moreover, we investigate the implications of the bounds on untagged Higgs decays on singlets decaying in the detector. For each method, we evaluate the potential of various hadron and lepton colliders to probe natural parameter space.
The paper is structured as follows. In Sec. 2 we present general properties of the singlet extension. In particular, we consider the renormalizable singlet in Sec. 2.1 and the special case of the relaxion in Sec. 2.2. Sec. 3 contains the collider bounds ordered by the scalar lifetime. We compare these complementary bounds in Sec. 4 before concluding in Sec. 5.

Singlet extension of the Higgs sector
In the following we present the phenomenological features of the real scalar singlet extension of the SM Higgs sector. After discussing the general properties of the scalar, we derive the masses and the relevant couplings for the renormalizable Z 2 breaking model, as well as for the non-renormalizable relaxion framework. This direct comparison will allow us to make a mapping from one model to the other.
The most general extended scalar potential of the Higgs doublet H and a new singlet scalar Φ is given by (2.1) where V (Φ) and µ 2 (Φ) are functions of the field Φ, whose exact forms depend on the model.
We do not consider direct couplings between Φ and other SM states besides the Higgs. In general, both fields can be written in the unitary gauge as where v = 246 GeV and φ 0 are their respective vacuum expectation values (VEVs), and their dynamical degrees of freedom are denoted by h and φ.
If Φ is not protected by an unbroken Z 2 symmetry, the singlet φ mixes with the Higgs.
In this case, φ inherits the Higgs couplings to the SM particles, suppressed by the mixing angle sin θ. At the same time, the couplings of the Higgs boson to the SM particles are reduced by a global factor cos θ. The mass matrix elements are Since for small mixing angles, the mass and interaction eigenstates are approximately the same, we use the symbols h and φ for both states throughout the paper and drop the label of the physical mass eigenstates used above. We denote the physical masses by m φ and m h , respectively, and their expressions will be given in the following sections. Due to the mixing, the φ production from and decay into SM particles equal those of a SM Higgs boson with the respective φ mass, modified by the mixing angle.
In addition, if the mixing angle or µ 2 (φ 0 ) are non-zero 1 , the Higgs couples to a pair of singlets. We denote this coupling of the mass eigenstates by c hφφ , which receives the following contributions where we use the shorthand notation s θ ≡ sin θ and c θ ≡ cos θ. In the two concrete models considered below, this expression simplifies substantially, especially in the limit of small mixing. When 2m φ < m h , the Higgs can decay via the new channel h → φφ with a decay width of (2.7)

Renormalizable singlet
The most general renormalizable form of V (Φ) and µ 2 (Φ) is This theory, often also called a Higgs portal model, has been studied extensively in the literature, see e.g. [19][20][21][30][31][32][33] and references therein. For later convenience, we choose φ 0 = 0, which can always be obtained by a shift of the φ field. This implies t = −a hφ v 2 from the minimum condition of the scalar potential V s . In this model the mixing angle is where x state = 4va hφ /∆m 2 state , with ∆m 2 int = m 2 hh −m 2 φφ being the difference of the diagonal entries of the mass matrix before diagonalization, and ∆m 2 phys = m 2 h − m 2 φ being the difference of the physical mass eigenvalues 2 . The approximation in the last step holds for vλ h . This corresponds to a large mass splitting between 1 For φ0 = 0 this requires that µ 2 contains a φ 2 term whereas for φ0 = 0 higher powers of φ are also a valid solution. These terms can be explicit or can arise from an expansion in φ. 2 Hence the difference in the squared physical masses can be expressed as ∆m 2 phys = ∆m 2 the singlet and the Higgs and a small mixing angle. The physical masses are For |a hφ | vλ h the masses are approximated as where the approximations rely on a small mixing and a large splitting of the diagonal entries of the mass matrix, exactly as the approximation in Eq. (2.10). Using in Eq. (2.6) the explicit expressions for V (Φ) and µ 2 (Φ) given in this section, we obtain the explicit coupling c hφφ where the approximation holds for small mixing and makes use of Eq. (2.10), and we define We use this as a parameter in the phenomenological investigations.

Theoretical bounds on the parameter space
The relevant phenomenology is described by the four physical parameters m φ , s θ ,λ hφ , and a φ . The parameters s θ [30,34],λ hφ , and a φ contribute to m φ , the former two at tree-level and a φ via a φ-loop. Therefore, their viable ranges are bounded by naturalness and depend The upper naturalness bound on λ hφ is then given by As we will see in the specific case of the relaxion, such naturalness bounds may be violated by orders of magnitude as a consequence of the cosmological evolution of the fields.

Relaxion
Unlike the generic Higgs portal model considered above, the relaxion scenario is designed to solve the SM hierarchy problem and is therefore much more constrained and predictive.
First, we briefly summarize the cosmological relaxation mechanism [23], considering the relaxion potential of the form Here Λ is a UV cutoff,M is the height of the backreaction potential (see below) and f is the relaxion oscillation scale. During its evolution, the relaxion scans the Higgs mass parameter µ 2 (Φ) from a large and positive value ∼ Λ 2 v 2 down to negative values. This scanning is a result of the slow-roll potential V (Φ), which is controlled by the small dimensionless coupling g, and r > 1/16π 2 which is bounded from below by the requirement of technical naturalness [35]. Once µ 2 (Φ) becomes negative, the Higgs gets a VEV and thereby activates a backreaction potential ∝ cos(Φ/f ), which eventually stops the rolling of the relaxion at a value φ 0 , where v(φ 0 ) = 246 GeV (see [36] for a recent discussion of the stopping mechanisms).
Such a theory naturally generates a large hierarchy between the electroweak scale and Λ, solving the SM naturalness problem 4 . In the following, we require The backreaction mechanism is model-dependent, and its most general potential is where we chose j = 2 and assume the minimal scenario of [23] (for alternative scenarios see 4 The relaxion does not solve the gauge hierarchy problem up to the Planck scale, and thus requires a UV completion to provide the needed Λ M Pl [37][38][39], and also to produce a large relaxion field excursion [40,41]. e.g. Refs. [23,35,[42][43][44][45] Making these substitutions in Eqs. (2.10) and (2.12), and omitting the term suppressed by the small coupling g, we obtain where we neglect small corrections to the Higgs mass m h . We notice that all the other couplings can be expressed as functions of m φ and s θ aŝ This means that this relaxion model has only two free parameters relevant for collider phenomenology, i.e. two less than the generic singlet case. The triple scalar coupling c hφφ can then be written as Hence, in contrast to the renormalizable singlet extension that hasλ hφ and a hφ as additional parameters, in the relaxion model this coupling is fully determined by m φ and s θ . Thus, the viable phenomenological parameter space is more limited and the model is more predictive.

Theoretical bounds on the parameter space
Naively, the general naturalness bound on s θ obtained in Eq. (2.19) applies also to the relaxion model. However, following Refs. [25,28], the dynamical evolution of the relaxion can fix the value of φ 0 at such a position that the two contributions to the relaxion mass in Eq. (2.32) cancel each other to a high precision, leading to a larger allowed value for sin θ for a given mass. In the following, we denote the number of a minimum by n.
First minimum The degree of such a cancellation is maximal in the first local minimum of the relaxion potential. There, in the limit ofM √ λ h v , the relaxion mass and mixing angle are given by (see Appendix A) 5 The mixing is maximized for maximalM and minimal f , namely f = Λ. (2.39) Generic minimum As mentioned above, the degree of tuning decreases if the relaxion stops in a later minimum. This may happen either through quantum fluctuations or by classical rolling [46,47]. In the limit of small tuning in a far minimum, n 1, sin(φ 0 /f ) ∼ cos(φ 0 /f ) ∼ O(1) and naturalness arguments lead to an estimate of the minimal value of 5 Obtaining m φ in the GeV range necessitates a large value ofM , and therefore the limitM √ λ h v is justified, see Eq. (2.32).
the mixing angle. In this limit, the mass can be approximated as while the mixing angle reads ExpressingM through the relaxion mass, and using the lower bound on f leads to a lower bound on the mixing, For the relaxion in such a minimum, and also for generic untuned Higgs portal models, the maximal mixing is given by Eq. (2.19). All the sin θ bounds derived in this section are valid up to order-one factors and thus should not be taken as exact.
Combined constraints: the relaxion band As follows from the above discussion, for each mass m φ there is a relaxion-specific lower and upper bound on sin θ.
where the mixing angle of the renormalizable singlet is natural (unnatural). The plotted λ hφ (m φ , sin θ) of the relaxion is defined in Eq. (2.33), within the natural sin 2 θ range from Sec. 2.2.2.

Collider bounds on (long-lived) scalar singlets
We present bounds on scalar singlets for a broad range of their lifetime. This necessitates a combination of various search strategies. Central to them is the lifetime which is shown in Fig. 2 for the relevant masses and mixing angles. For • short lifetimes, untagged Higgs decays into a pair of singlets lead to strong indirect bounds; • intermediate lifetimes, displaced vertex (DV) searches and strategies based on timing information probe a broad range of the parameter space; • long lifetimes, the singlet escapes the detector and can account for invisible signatures.
We compare these bounds to the ones previously studied from direct searches in Z decays and from associated Zφ production. The presented bounds are based on singlet pair production via Higgs decays (h → φφ). The production via singlet-Higgs mixing is negligible for the parameter region considered here, for details see Appendix B. Our bounds apply to the general singlet extension of Sec. 2.1. We will point out which regions of the displayed parameter space can also be realized by the relaxion.

Fits of untagged and invisible Higgs decays
Beyond the Standard Model (BSM) physics can modify the tagged Higgs branching ratios both by modifying the Higgs couplings to SM particles by κ x = c x /c SM x , and by introducing new decay channels for the Higgs, depleting the relative SM contribution to the total decay width [29,48] The BSM particles produced in these Higgs decays can either decay visibly, or remain invisible. While searches for Higgs decays with missing energy directly constrain the invisible branching BR inv , these search results can also be used as a tagged category in a fit. In contrast, the final states of the visible BSM Higgs decays (e.g. light jets) are generally not included in the list of tagged visible decays (such as h → τ τ, bb, V V, ..., explicitly displayed e.g. in Tab. 1 of Ref. [49]). Hence they remain untagged 6 , and the corresponding Higgs branching BR unt is not determined by any specific search, but by the uncertainties of the tagged channels. Therefore, global fits of the Higgs coupling modifiers κ x to measured signal strengths times tagged branching ratios normalized to the SM prediction), together with searches for invisible Higgs decays, allow to constrain the Higgs decay width into BSM particles 7 , 6 For the implications of the direct searches for h → φφ further decaying promptly into four visible particles, e.g. h → 4b, h → bbτ τ , see Refs. [29,[50][51][52]. 7 The SM contributions to Γinv and Γunt are subtracted.
The global Higgs fits performed in the scope of the European Strategy Update [53] present results for the future hadron colliders HL-LHC, LHeC, HE-LHC and FCChh, as well as for the lepton colliders ILC, CLIC, CEPC and FCCee running at different energy stages. Here we apply the results from the so-called kappa-2 scenario that treats BR inv and BR unt as free parameters for each collider individually. In addition, it has several independent κ x whereas in the general singlet and the relaxion models there is only one overall κ ≡ cos θ, see also Ref. [29]. Furthermore, in the region of intermediate and high Fig. 2), all φs decay inside the detector, hence BR inv = 0, and fitting only two parameters, κ and BR unt , would be sufficient. In the opposite case of very small sin θ, the Higgs couplings to SM particles become SM-like (κ 1), and fitting only BR inv would be enough. Hence, the multiparameter fits used in Tab. 1 and in Fig. 3 give rise to conservative bounds on this actually more predictive model, defined by less parameters. To evaluate the gain in sensitivity by fitting only the needed parameters, we also include the dedicated fit results performed for the CLIC stages [50], see the lower part of Tab. 1.
A combination of the ATLAS and CMS data collected in Run-1 results in a limit on BR BSM < 20% [48] (which can applied be as a conservative bound on BR unt ), comparable to that of ATLAS alone in Run-2 of BR unt < 21% [54]. The strong result of the Run-1 combination, despite the smaller summed luminosity, is due to the fit of only a global κ and BR BSM . A Run-2 combination or a dedicated 2-parameter fit will be able to exclude further parameter space based on the already existing data.
In Fig. 3 we show the constraints on the m φ -sin 2 θ parameter plane of the relaxion. In addition, we show in gray the natural relaxion band, whose upper and lower sin θ limits are discussed in Sec. 2.2. The experimental limits and projections result from requiring where the partial width Γ h→φφ ∝ c 2 hφφ is given in Eq. (2.35), and the total Higgs width in the SM is Γ SM tot = 4.1 MeV [55]. The contours form horizontal and vertical asymptotes determined by the sin 2 θ and m φ contributions to c hφφ , respectively. When neglecting the kinematical mass dependence of Γ h→φφ (for m φ m h /2) and the BSM contribution to the total width, the location of the asymptotes for the relaxion can be approximated as The shaded blue area is already ruled out by Run-1 of the LHC, excluding natural mixing angles of heavy relaxions above m φ 18 GeV. As indicated in Tab. 1, this Run-1 bound is in fact on BR BSM , whereas for such large values of sin 2 θ all relaxions decay inside the detector. Hence, a specific bound on BR unt will exclude also lighter relaxions.
The strongest bound will be reached by the FCChh, excluding m φ 10 GeV and sin 2 θ 3 · 10 −3 . As indicated by the dash-dotted yellow lines, the fit of only BR BSM , assuming all κ x = 1, for CLIC leads to a significant improvement of the bound compared to the multi-κ fit at CLIC 8 . Dedicated fits for the FCCee and FCChh could have the potential to close the high-mass relaxion window above few GeV.
The situation is different for the general singlet model where λ hφ is a free parameter and BR h→φφ varies with the choice of λ hφ . For a larger value of λ hφ than the one predicted within the relaxion framework, the bounds from untagged Higgs decays can become even stronger, whereas they get reduced to the sin 2 θ dependence in Eq. In general, while the bounds on BR BSM hold for arbitrary values of sin θ, the more specific bounds on BR unt are valid as long as the decay length is significantly smaller than the detector size. Conversely, the bounds on BR inv apply to decay lengths clearly exceeding the detector size.

Displaced jets
The singlet can be detected in searches for Higgs decays into displaced jets if it is sufficiently long-lived, but still decays in the detector. ATLAS searches [56][57][58] and FCC-ee projections [59] provide upper bounds on the branching ratio BR h→φφ as a function of the proper decay length cτ φ for a few singlet masses 9 . We transform them into upper limits on λ hφ as a function of sin 2 θ, for the corresponding mass points given in the analyses, shown in Fig. 4. The dashed lines show the upper limit on λ hφ from naturalness, see Eq. (2.22).
While for m φ = 5 GeV the ATLAS searches do not constrain any natural parameters of the singlet model, for higher masses the searches already probe parts of the natural parameter region. In contrast, FCC-ee will access natural parameter space for all masses.
The displayed FCC-ee bounds show the combination of the two analysis strategies from Ref. [59], and therefore span a larger range of sin 2 θ. The CLIC sensitivity to a long-lived scalar singlet via displaced vertex searches was studied in Ref. [61] and is included in our overview plot in Fig. 9. The comparison shows that CLIC and FCCee provide a comparable sensitivity. 8 Strictly, this fit is applicable only for vanishing sin 2 θ, but in any case the exclusion contour of CLIC380 (CLIC3000) reaches only sin 2 θ 3 · 10 −3 (1.5 · 10 −3 ) corresponding δκ ≡ 1 − cos θ 1.5 · 10 −3 (7.6 · 10 −4 ), i.e. just below the resolution of κ, see Tab. 1. Hence we use this fit as an illustration of the gain of sensitivity in a suitable fit of such a predictive model 9 For Higgs decays into complex singlets at the LHeC, see the recent Ref. [60].

Delayed jets
A powerful strategy to search for long-lived particles was recently presented in Ref. [62], allowing to detect displaced vertices in the CMS tracker 10 . This proposal utilizes the timing detector layer, to be installed at the high luminosity (HL)-Large Hadron Collider (LHC) [64], to identify secondary vertices by the delayed arrival, ∆t, of the light decay products, compared to the arrival time expected for a directly travelling SM particle.
An initial state radiation (ISR) jet is used to time-stamp the collision. Ref. [62] provides the bounds for the benchmark scalar masses of m φ = 10 GeV and 50 GeV at the HL-LHC. In order to determine the mass dependence of the experimental reach, we simulate Higgs events at the LHC and FCC-ee, using MadGraph5 [65] at leading order (LO), where 10 The new proposal in Ref. [63] evaluates the sensitivity of the High Granularity Calorimeter of the CMS detector upgrade to the same type of h → φφ decays. While the conservative estimate yields bounds comparable to the timing bounds of Ref. [62], only the analysis assuming a displaced track trigger could improve them. the Higgs decays by h → φφ, and each scalar decays through φ → jj. Subsequently, we implement the search strategy presented in [62], reproduce its results, and apply it to the additional mass points. For the FCC-ee, we assume a (hypothetical) timing detector comparable to the one planned for the HL-LHC. The detection efficiency is mostly affected by demanding a long time delay of the jet produced in the singlet decay, related to the singlet's path through the detector, along with requiring the singlet to decay between the inner tracker and the timing layer. Hence, the selection criteria for this search are mainly geometrical. Therefore, for each event kinematics and for each jet j in the event, we find the range of lab frame singlet decay lengths l φ for which an event will be accepted. Since the detection of a single delayed jet is sufficient, each event is then weighed by the event      to decay within the allowed region, which is calculated from an exponential distribution More details on the calculation, as well as on the resulting efficiencies and the expected upper limits on BR h→φφ as a function of cτ φ can be found in Appendix C.
The interpretation of these bounds in terms of the singlet parameters λ hφ and sin 2 θ is shown in Fig. 5. While the HL-LHC probes natural values of λ hφ for m φ > 5 GeV, at the FCC-ee this is the case only for slightly higher masses. As this analysis has almost zero background in the signal region of ∆t > 1 ns (for details see Ref. [62]), its sensitivity is determined by the number of Higgses. Therefore, the HL-LHC appears to perform better than the FCCee. Since it is the hadronic environment at the HL-LHC that necessitates this restrictive cut on ∆t, the FCCee can allow for a looser cut, and the limit presented here based on the HL-LHC cut is conservative.

Searches for invisible Higgs decays
If the proper decay length of the scalar is larger than, or comparable to, the size of the detector, the scalar may give rise to missing energy. Global Higgs coupling fits set strong bounds on BR h→inv [53]. These can be interpreted as bounds on λ hφ in the limit of vanishing where the sum runs over all h → φφ events passing the selection criteria when an infinite decay length is assumed, p is the momentum of each scalar, L is the distance the scalar travels inside the detector, and the indices {1, 2} mark the two scalars produced in the Higgs decay. A conservative estimate of the rescaled bound can be given by minimizing r for each search. For LHC searches, which require a large missing p T , this can be approximated by where L T is the transverse detector size and p miss T is the minimally required missing transverse momentum. For lepton colliders, such as the FCC-ee with a lower √ s = 240 GeV, a better approximation is given by setting the energy of each scalar to m h /2, as the Higgs is produced at low momentum, i.e. r consv.
For a more precise estimate of the bounds, we determine r for each search in Tab. 2.
We use MadGraph5 [65] to simulate the leading signal process in each search at LO. We then apply their selection cuts, and obtain the distribution for each signal mass, and subsequently obtain r following Eq. (3.6). The signal processes and selection cuts applied are summarized in Tab. 3. The resulting r for the HL-LHC and FCC-ee is shown in Fig. 6 as a function of cτ φ . For a given m φ and cτ φ , r is larger for the HL-LHC because the L/p distributions peak at lower values than for the FCC-ee.
The CMS bounds as well as the HL-LHC and FCCee projections on λ hφ and sin 2 θ are  shown for different values of m φ in Fig. 7. In general, each contour has a horizontal and a vertical asymptote, driven by the limit on BR h→φφ and by the lifetime, respectively. As a consequence, the horizontal asymptotes are hardly mass dependent (apart from m φ = 50 GeV which is near the decay threshold), whereas the reach in sin 2 θ is larger for low m φ -owing to the longer lifetime. While for m φ = 5 GeV no natural parameter space is probed, for m φ = 10 (15) GeV only FCCee (and HL-LHC) access the natural parameter space, and for higher masses this is also achieved in the present CMS analysis.
For the FCChh, the vast amount of produced Higgses can result in a very strong upper limit on the invisible branching ratio. Ref. [66] reports for a luminosity of 30 ab −1 an expected sensitivity of a direct search to BR h→inv 3 · 10 −4 , i.e. similar to the result from a global fit of BR h→inv ≤ 2.4 · 10 −4 [53]. The asymptotic limit on λ hφ for vanishing sin 2 θ can be approximated as This translates into the asymptotic bound on λ hφ for m φ = 5 GeV (50 GeV) of λ hφ ≤ 2.3 · 10 −4 (2.9 · 10 −4 ) using the fit result, hence stronger than the limit of the direct searches for h → inv at the FCCee, and probing natural values of λ hφ throughout this mass range.
The approximate FCChh bounds are included in Figs. 8 and 9 up to values of sin 2 θ for which all singlets can be safely assumed to decay outside the detector.    Table 3. Signal processes and selection cuts applied in the calculation of the fraction r of invisible signal events. The p j T cuts refer to the leading (subleading) jet.

Overview
Having presented details about each search strategy in the previous section, here we compile them for comparison, to highlight their complementarity and to evaluate the probed parameter regions, both for the general singlet and the relaxion.
In Fig. 8  . In contrast, for the relaxion, the accessible λ hφ within the natural band of sin 2 θ is confined to the dark blue line that extends to larger sin 2 θ than in the renormalizable singlet model, see Fig. 1.
For both models, λ hφ only impacts the decay of the Higgs into a pair of singlets, i.e. the number of produced φs, whereas sin θ mainly determines their lifetime τ φ , and only contributes to BR h→φφ for high sin θ. Because the specific decay of φ does not play a role, the shape is entirely determined by the λ hφ and sin θ contributions to the coupling c hφφ in BR h→φφ . The green vertical lines represent the LEP1 bound [71] for the rare Z → φ decay, and the GigaZ and TeraZ projections we obtained by rescaling with the ratio of produced Z bosons, or the bound on e + e − → Zφ at LEP2 [72] and ILC [61] which are stronger than the respective Z-decay constraint for m φ = 50 GeV [29].
The natural parameter space of the general singlet with m φ = 5 GeV has not been probed yet. Only small fractions of it can be probed by timing and displaced searches, as well as by fitting the untagged and BSM Higgs width and by searches of rare Z-decays.
For the higher masses considered, all investigated bounds contribute to probing the natural parameter space, mainly because the upper naturalness bounds increase with the mass.
Considering the relaxion at m φ = 5 GeV, so far only the Z-decays at LEP1 marginally constrain the upper end of the natural region, which can be further probed by the same process at GigaZ, and excluded by TeraZ. Furthermore, untagged Higgs decays at future colliders are sensitive to the natural relaxion parameters. The heavier relaxion examples are already excluded by the BSM Higgs decays at the LHC1.
In Figs. 9 and 10 we show the bounds in the m φ -sin 2 θ plane for the singlet scalar and for the relaxion, respectively. For the singlet scalar, we set the coupling λ hφ = m 2 φ /v 2 =λ max hφ , hence λ hφ could be even larger. For the relaxion, the value of λ hφ is given by Eq. (2.33).
In addition to the bounds discussed above, we also show the direct bound for m φ < 5 GeV from B → Kµµ at the LHCb [16,73,74]. Furthermore, we translate the uncertainties δκ of the Higgs coupling modifier in global fits 11 into model-independent bounds on sin 2 θ that are independent of m φ and λ hφ . The strongest bound stems from δκ Z at the FCChh (see Tab. 1), and is shown in Fig. 10, but omitted in Figs. 8 and 9. From Fig. 10 we see that relaxions heavier than ∼ 18 GeV are already excluded by the current LHC bounds on BSM Higgs decays. Rare Z-decays from LEP1 probe parts of the natural parameter space of the relaxion for m φ 5 GeV, but the bound from the BSM Higgs branching at the LHC Run-1 is stronger than this LEP1 bound for m φ 15 GeV. The best bounds from untagged Higgs decays will come from the FCChh, and can exclude relaxions above m φ 8 GeV. On top of that, TeraZ can exlude relaxions of m φ 3 GeV. 11 We obtain the approximate 95% CL bound on sin 2 θ from the provided 68% CL bound on δκ with κ = 1 + δκ by sin 2 θ (95) 1 − (1 + r δκ (68) ) 2 where r = q

Conclusions
In this work, we exploit the sensitivity of the exotic Higgs decay channel h → φφ to parameters of the relaxion and singlet models, taking into account existing searches and global fits, as well as projections for future colliders.
We discuss the renormalizable, non-Z 2 -symmetric singlet extension of the SM, focusing on the exotic Higgs decay h → φφ via the triple scalar coupling c hφφ . The collider phenomenology is determined by the four parameters m φ , sin θ,λ hφ , and a φ . Beyond the usual naturalness bound on the mixing angle, we present naturalness bounds onλ hφ and a φ and investigate their implication on the physical parameter space. Moreover, we provide a matching between the singlet parameters and those of the relaxion. Here, the absence of a Z 2 symmetry is pivotal to accommodate the linear slow-roll relaxion potential. The h 2 φ 2 term in the singlet model maps onto the first term of the expansion of the backreaction potential. We extend the naturalness relaxion band to higher masses relevant at colliders, where it is described by only two parameters, m φ and sin θ, which determine λ hφ . Consequently, the relaxion model is both more constrained and predictive than the renormalizable singlet extension.
The lifetime of φ, given by sin θ and m φ , is the crucial handle in determining the kind of search strategy that sets the strongest bound. We study various lifetime dependent strategies. In particular, we evaluate the limits from global coupling fits on the new Higgs branching ratio into BSM, split into untagged and invisible final states; interpret the searches for the Higgs decaying into displaced jets in terms of the singlet model; exploit the time delay of jets originating from the φ decay to derive bounds in the region of intermediate lifetime; constrain the region of low sin 2 θ by searches for invisible Higgs decays and pay attention to the range where the decay lengths are of the detector size such that only a fraction of the particles actually gives rise to the invisible signature.
Our main phenomenological findings are: • For m φ = 5 GeV, only a small fraction of the natural singlet parameter space can be probed. For higher masses, larger coupling values become natural and the LHC has already excluded parts of it.
• The FCC can probe almost the complete considered parameter region by combining TeraZ, FCCee and FCChh, unless λ hφ is much smaller than used here.
• The natural range for relaxions heavier than 18 GeV is already excluded by searches for untagged Higgs decays at the LHC. The FCCee has the potential to exclude relaxions down to 8 GeV using the same strategy. Only the search for rare Z decays at TeraZ will be able to exclude the full mass range for heavy relaxions with m φ > 3 GeV.

Acknowledgments
We thank Gilad Perez for the initial collaboration, many valuable discussions, and helpful comments on the draft. Moreover, we thank Zhen Liu for useful information on timing In the following we set r = 1 for simplicity, given that the exact expression for the relaxion mass only mildly depends on r. As the relaxion rolls down its potential before stopping, during each relaxion period ∆φ = 2πf the maximal (in absolute value) slope of the oscillatory potential V br changes by where θ = φ /f denotes the relaxion angle at which the V br slope is maximized within the given 2πf period, i.e. the inflection point of the periodic potential. ∆V br isM 2 /Λ 2 suppressed with respect to the V br overall size at the stopping point (A.1). Close to the final minimum, θ can be found using Eq. (2.32) for m 2 φ V (φ ), and solving by a Taylor expansion of µ 2 (φ) neglecting the term suppressed by g. After the first minimum is formed, the slope of the periodic potential can overcompensate V sr only by ∆V br . After the n-th minimum it can do so by n ∆V br . Correspondingly, the slope of the overall potential is given by the same value We therefore know the position of the inflection point φ and its slope V (φ ). They can be used to find the properties of the closest minimum φ 0 located before φ . The value of V (φ 0 ) can be expressed as a Taylor series around φ Note that V is obtained from the effective relaxion potential V eff after integrating out the Higgs boson, h 2 → −µ 2 (φ)/λ h , which is given by and consequently all the related parameters of the theory. In particular, the relaxion mass can be approximated as As we see, the mass is proportional to V (φ ), which itself carries a factorM /Λ. This is precisely the reason why the relaxion mass is suppressed with respect to the naive estimate.
In this paper, we are interested in the corner of the parameter space where the relaxion reaches its maximal possible masses, which requires takingM v. In the limitM √ λ h v, applicable within the relaxion mass range considered in this work, the relevant expressions simplify to Inserting the relaxion angle θ 0 into the general expression for the relaxion mass in Eq. (2.32), we see that the small relaxion mass appears as a result of a fine cancellation between two contributions. Note that this also means that the loop corrections, otherwise subleading, may contribute sizeably to the relaxion mass. This, however, should not change qualitatively the results that we have derived, as the presence of the relaxion mass suppression is directly linked to the slow growth of the periodic barriers amplitude-∆V br /V br 1-the feature which is not expected to be altered by the loop effects.
For completeness we also write down corresponding expressions in the opposite limit, M v, relevant for lighter relaxion, which were derived in Ref. [28] 12

B Estimating singlet production via Higgs mixing
For small values of the coupling λ hφ s 2 θ m 2 h /(2v 2 ), the branching ratio BR h→φφ is proportional to sin 4 θ, cf. Eq. (2.17). If in addition sin θ is small, the Higgs almost never decays into a pair of scalars. On the other hand, the production of scalars via their mixing with the Higgs only scales as sin 2 θ and becomes the dominant production mechanism if λ hφ is small. However, if a sufficiently long lifetime is required in order to have a handle for the considered analyses, we estimate in the following that production via mixing yields only few events making a dedicated search difficult.
The number of scalars produced via mixing is given by n mix = Lσ φ s 2 θ , where L is the luminosity and σ φ is the production cross section of a Higgs boson with mass equal to m φ .
Since detecting a dijet resonance at low mass is extremely challenging, we will consider only the searches for displaced jets or missing energy. To obtain a displaced or invisible signature, we need cτ 1 cm ( 1 µm) for the HL-LHC (FCCee) which translates into sin 2 θ 10 −9 ( 10 −5 ) using Fig. 2 for m φ = 5 GeV. A higher value for m φ would be helpful in an analysis, but at the same time require even smaller mixing angles and also imply a smaller production cross section for kinematical reasons.
The HL-LHC will collect a luminosity of L = 3·10 6 pb −1 . The production cross sections for a light Higgs at the LHC are below 100 pb for all modes except for gluon fusion without p T requirement [29]. A leading order parton-level estimate with MadGraph5 for φ + j production at 14 TeV with a very mild p T > 20 GeV requirement for the scalar yields σ φ ≈ 120 pb. Therefore the HL-LHC can only produce n HL-LHC mix 0.4 scalars. Consequently, even before significant selection cuts no events will be available for an analysis.
The FCCee on the other hand will collect L = 5·10 6 pb −1 and the dominant production mode for a light Higgs at FCCee is Higgs-strahlung with a cross section of about 0.6 pb [29] with p φ T > 10 GeV. Therefore n FCCee mix 30. Considering more selective cuts on top of the minimal example cut applied here as well as the detector acceptance and e.g. leptonic Z decay branching ratios, it will be impossible to have a sufficient number of scalars left for an analysis.
Here we argued why we focus only on φ production via Higgs decays. However, progress in detecting promptly decaying low-mass resonances may provide a new channel for singlet and relaxion searches [29]. Especially bb, τ τ or µµ decays from production via mixing may allow to constrain the parameter regions where λ hφ is negligible.

C Timing of delayed jets
The crucial requirement of the analysis proposed in Ref. [62] is that a jet leaving no track in the inner tracker hits the proposed timing layer with a delay ∆t > 1 ns with respect to a (hypothetical) SM jet, going directly from the interaction point to the same location on the timing layer. This signature can be achieved by a particle that is invisible to the inner detector and decays into SM hadrons between the inner tracker and the timing layer. The delay then is a result both of the lower velocity of the heavier decaying particle, and of the displacement of the secondary decay in which the hadron is produced. For this reason, the acceptance probability of a given event is dominated by the geometrical trajectory of the decaying scalar and its decay product, once the kinematics is determined. Namely, once the four-momenta of the scalar and the jet are set, the lab-frame decay length of the scalar determines the secondary vertex position, the position in which the final jet hits the timing layer, and the overall time delay.
Since the analysis only requires at least one delayed jet, we can consider the four final state jets from the decay chain h → φφ → 4j independently. Then, for a jet in a given event, we can find the range of allowed lab frame decay lengths of the scalar l φ,j , for which the jet will be accepted as signal. If this range is non-empty, we can assign a weight w j , calculated as the probability to obtain l φ,j within the allowed range, given that the proper decay length is cτ φ , as in Eq. (3.5). The probability for the whole event to be accepted is then given by event = 1 − (1 − w 1 )(1 − w 2 )(1 − w 3 )(1 − w 4 ).
In the following we will explain the computation of the allowed range of l φ,j . As described above, the scalar needs to decay between the outer radius of the inner tracker L 1 and the outer radius of the timing layer L 2 . For CMS L 1 = 0.2 m and L 2 = 1.17 m [62], and for the FCCee we assume L 1 = 0.127 m and L 2 = 2.1 m [75]. Thus, the distance the scalar may travel before decaying is constrained by l L 1 ≤ l φ,j ≤ l L 2 , given by where θ φ is the polar angle between the beam axis and the three-momentum of the considered scalar. In addition, we demand that the displaced jet does not cross the inner radius L 1 towards the beam axis, as it will leave a signature in the tracker. We thus require l min φ,j = L 1 sin θ φ max 1, − sign(cos ϕ) | sin ϕ| , where ϕ ≡ ϕ φ − ϕ j , and ϕ φ and ϕ j refer to the azimuthal angles of the scalar and the jet, respectively.
The main selection criterion of the search is the time delay of the decay product, which is a result of the displaced vertex. The delay is defined as where l φ is the distance traveled by the scalar before it decays, l j is the distance traveled by the decay product (a jet, in our case) to the timing layer, and l SM is the distance a hypothetical SM particle would travel directly from the interaction point to the timing layer. The velocities of the particles are denoted by β φ , β j and β SM in units of the speed of light c. Because the SM hadrons are light, β SM = 1 to a good approximation. By demanding that the delayed jet hits the timing layer at radius L 2 , and by setting l SM = | l φ + l j |, the time delay can be expressed solely as a function of the event kinematics and l φ . As the time delay has at most one maximum as a function of l φ , the allowed decay should lie between l ∆t max , l ∆t max , given by solving Eq. (C.3) for l φ with the required minimal time delay. Note that Eq. (C.3) can be brought to a 4th-degree polynomial form in l φ , and thus its roots can be found analytically. The temporal resolution of the timing layer is simulated by assigning normally distributed time stamps to the displaced jet hit δ j and to the SM-ISR hit δt ISR , smeared by σ = 30 ps [64], and requiring ∆t th ≤ ∆t + δt j − δt ISR , where ∆t th = 1 ns is the minimal time delay set by the analysis.
Lastly, the decay product of the scalar should hit the timing layer at L 2 within the length of the detector (in theẑ direction), where we set |z max | = 2.6 m at CMS and |z max | = 2.3 m at the FCCee. If the scalar is produced at z 0 , then the z position of the hit of its decay product is Z ≡ l φ cos θ φ + l j cos θ j − z 0 , (C. 4) which is yet again completely determined by the event kinematics and l φ (we set z 0 = 0 for simplicity, as the variations in the exact primary vertex position are negligible compared to the detector length). Therefore, imposing −|z max | ≤ Z ≤ |z max | and solving for l φ yields another set of constraints. Note that since Z can have at most one extremum as a function of l φ , there may be at most two disconnected allowed ranges of l φ satisfying the requirement above.
The final range of allowed decay lengths is then set by the union of the constraints given by the conditions above. For each allowed continuous range of l φ , w is calculated by If the union has two or more disconnected regions, their contribution to w should be summed. The resulting bounds on the Higgs branching to a pair of scalars and the efficiency of the search, both as a function of the lifetime, are presented in Fig 11.     . BR(h → φφ) and efficiency as a function of cτ φ for a search for delayed jets.