Cosmological beam dump: constraints on dark scalars mixed with the Higgs boson

Precision cosmology provides a sensitive probe of extremely weakly coupled states due to thermal freeze-in production, with subsequent decays impacting physics during well-tested cosmological epochs. We explore the cosmological implications of the freeze-in production of a new scalar $S$ via the super-renormalizable Higgs portal. If the mass of $S$ is at or below the electroweak scale, peak freeze-in production occurs during the electroweak epoch. We improve the calculation of the freeze-in abundance by including all relevant QCD and electroweak production channels. The resulting abundance and subsequent decay of $S$ is constrained by a combination of X-ray data, cosmic microwave background anisotropies and spectral distortions, $N_{\rm eff}$, and the consistency of BBN with observations. These probes constrain technically natural couplings for such scalars from $m_S \sim$ keV all the way to $m_S \sim 100$ GeV. The ensuing constraints are similar in spirit to typical beam bump limits, but extend to much smaller couplings, down to mixing angles as small as $\theta_{Sh} \sim 10^{-16}$, and to masses all the way to the electroweak scale.


I. INTRODUCTION
A pragmatic approach to searching for new physics is to focus on generic interactions that have the potential to be detected experimentally with current or upcoming technology. Classifying the interactions of new neutral states with the Standard Model (SM) according to the dimensionality of the couplings, there are only three 'portal' operators that are unsuppressed by a new energy scale. The so-called scalar, vector and neutrino portals could provide the leading connection to a hidden or dark sector, motivated by considerations of neutrino mass and dark matter, but possibly comprising a rich structure of yet-unseen particles and forces [1].
The three portals have recently been under intense experimental scrutiny (see e.g. Refs. [2,3]), with a forthcoming program to increase sensitivity into unexplored regions of the parameter space. While collider and beamdump experiments provide sensitivity to relatively large portal couplings, astrophysical phenomena and cosmology can provide complementary reach to much weaker couplings. Constraints generically arise as follows: thermal production of new states in the very early Universe can occur with a sub-Hubble rate (a process often called 'freeze-in'), which necessarily leads to a small but nonnegligible abundance of such particles in the thermal bath. If the lifetime of these particles is large, they may survive to later epochs, and decay during or after Big Bang Nucleosynthesis (BBN) altering the light element yield. Longer lifetimes may lead to decays during or after the formation of the Cosmic Microwave Background (CMB), potentially altering the detailed features observed in precision CMB experiments.
The origin of these cosmological constraints is reminiscent of the detection strategy behind a generic particle beam dump experiment. Typically, very energetic particles in the beam initiate the production of exotic states in the target, which then propagate through a rock or dirt filter, and decay/scatter in a relatively backgroundfree environment inside the detector, thus generating a signal of anomalous energy deposition. In the cosmological setting, the analogue of the initial beam-on-target is the stage of the very early hot Universe, the analogue of propagation through a dirt filter is the long stage of subsequent expansion and cooling, as the Universe evolves into a well-understood stage associated with BBN or the CMB, which is then a direct analogue of the calorimetertype detector that measures abnormal energy deposition. Therefore, it is appropriate to name this method of studying rare long lived particles the cosmological beam dump.
Cosmological constraints of this kind were first applied to the heavy neutral lepton (HNL) portal, and were considered in a number of publications [4][5][6], resulting in stringent constraints on sterile neutrino degrees of freedom N . Cosmological constraints on the ultra-weak regime of the dark photon parameter space (a new particle A µ that has an F µν F µν coupling to photons) were explored in Ref. [7] (see also Ref. [8]). For the very small coupling constants relevant for cosmological probes, dark photons never thermalize and the abundance is determined by freeze-in production via inverse decay reactions. Subsequent energy injection from dark photon decays can alter the path of BBN and the CMB, and agreement with precision observations excludes certain disconnected regions in the parameter space. Unlike the case of HNL, thermal production of dark photons may not exhaust all channels, as bosonic states can also be copiously produced during inflation. This extra production channel has, however, a wide range of possible outcomes depending on the Hubble scale during inflation. In that sense, the limits presented in [7] are conservative.
A similar study can be carried out for the scalar portal. Unlike the cases of HNL and dark photons where the leading portal operator is of dimension 4, the neutral scalar S can have a dimension three coupling to the Higgs bilinear, SH † H. This represents the only superrenormalizable portal that exists between the SM and any potential dark sector. Moreover, the radiative corrections to the scalar mass created by such an interaction can be naturally under control [9], ensuring technical naturalness of a small m S . There has been significant attention paid to this interaction over the last few years, due for example to the idea of cosmological self-tuning of the Higgs mass through scanning via the small SH † H interaction [10].
With the SH † H interaction alone, S is guaranteed to be produced in the early universe and its subsequent decay may occur during cosmological epochs that are sensitive to energy injection, as was pointed out in Ref. [11] in an application to BBN. It is therefore interesting in our search for new physics to investigate the phenomenology of this interaction in the super-weak regime, in which the abundance is determined by the freeze-in mechanism. The electroweak era can be identified as the main contributor to the freeze-in abundance of S, at temperatures where t(t), W ± , Z, h are thermally excited, due to the preferential coupling of the SM Higgs to heavy particles. This was first recognized in [12] in the context of the quadratic S 2 H † H interaction for m S > 1 MeV, but the same conclusion holds for the SH † H interaction. This is markedly different from N and A freeze-in production [13,14], where the peak occurs at noticeably different temperatures T max (for kinematically accessible particles): The goal of the present work is to determine the cosmological constraints on S, due to its Higgs portal coupling. This requires a computation of the cosmological S abundance due to freeze-in production. Existing estimates of S freeze-in have considered QCD production via top quarks, leading to Y S ∼ 1.6 × 10 12 θ 2 [8], where θ is the mixing angle between S and the SM Higgs, and also a lower bound Y S 2.9 × 10 9 θ 2 on the abundance from Primakoff and Compton processes at low temperatures T < 20 GeV, in the context of relaxion-Higgs mixing [15]. (Further in-depth considerations of the freezein production due to quadratic and linear couplings of S were performed very recently in [16], and in a more generic setting in [17].) Our analysis of the S abundance and decays indicates that QCD and electroweak processes are both significant, and the conclusions are summarized below: • Freeze-in yield: The tree-level freeze-in production of S is computed for all electroweak and QCD channels, with a T -dependent electroweak vev v(T ) used as the first approximation of the relevant thermal effects and to provide an estimate of the pre-cision of the calculation. Solving the Boltzmann equation numerically, and incorporating a full set of QCD cross sections results in a reduction of the total QCD yield relative to the channel analyzed in [8]. We also assess the accuracy of the Maxwell-Boltzmann approximation in the production calculation, leading to the following result for the abundance from QCD and electroweak processes for the m S M W case: This estimate becomes more uncertain for masses m S ∼ 100 GeV, i.e. for masses close in value to the temperature/energy scale of the electroweak phase transition.
• Decay rate: We find that there is a sizeable uncertainty in the constraints for m S in the QCD range, due to the poorly known S decay rate to pions and kaons. We show two decay profiles in this case. Additionally, we improve the calculation of the S → γγ decay rate, which is important for m S < 1 MeV, by incorporating the light quark contribution via mesonic loops, thus decreasing (or increasing) the decay rate by a factor of 4 over the rate used when u, d, s are assumed to be massive (or are neglected).
• Early decays: Ref. [8] performed a thorough analysis of the BBN constraints, but noted that their analysis of early decays for m S < 2m π did not consider energy density considerations due to the large stored energy in the S bath. We include a treatment of early decays, transitioning from the freezein abundance to the thermalized freeze-out relic, and consider the impact on the relative neutrino and electromagnetic energy baths.
In what follows, we first review the model and describe its features in Sec. II. The freeze-in abundance calculation is described in Sec. IV, including details of several subtleties. We provide a complete scan of the S parameter space at small mixing angles in Sec. V, with details of the cosmological constraints updated in this work. The results are summarized in Fig. 1, which shows that precision cosmology provides an efficient probe of the parameter space of the model many orders of magnitude inside the region where it is "technically natural" (i.e. not plagued by the issue of fine tuning). Finally, we conclude the paper with a general discussion of the robustness of these results and final conclusions in Sec VI. Several technical results are relegated to Appendices.

II. THE SUPER-RENORMALIZABLE HIGGS PORTAL MODEL
We consider a subset of the minimal Higgs portal model, the super-renormalizable Higgs portal. The scalar 1. An overview of the excluded parameter space for the super-renormalizable Higgs portal scalar, including the updated constraints from this work due to the diffuse X-ray background (XRay), CMB anisotropies, spectral distortions, N eff and BBN. Constraints from new short-range forces (Force) [52][53][54][55] and stellar cooling (Stellar) [51] from other authors are also shown. We also display the projected SHiP sensitivity [2] and an estimate of supernova (SN) constraints [56].
part of the SM Lagrangian involving the Higgs doublet H is augmented by a mass term for the singlet S and a dimension three interaction: The A term induces a small mixing angle θ between the physical excitations S and h. At linear order in A, the mixing angle is given by and leads to Yukawa interactions between S and SM particles, equivalent to the SM Higgs boson interactions rescaled by the suppression factor θ. In the unitary gauge for the broken electroweak phase, after diagonalizing to find the physical states h and S, we have the scalar potential which exhibits the hhS and hhhS contact interactions. The S sector could include additional self interactions, e.g. λ 3 S 3 . Such self interactions would combine with the A-term and contribute to the S freeze-in production via Higgs boson decays h → SS. Large self interactions can also influence the S metastable abundance after freeze-in by maintaining thermalization of the dark sector prior to subsequent decays [18,19]. We will neglect this type of interaction and focus on the pure freeze-in regime of Lagrangian (3).

A. S decay rate
The S decay rate has well-known theoretical uncertainties associated with mesonic decay channels in the mass range 2m π < m S < 4 GeV [20]. We follow Ref. [21] and use two different decay models to demonstrate the magnitude of the theoretical uncertainty in the final S freeze-in parameter space. The baseline decay model matches low-energy theorems near the pion threshold to a ππ phase-shift analysis above 600 MeV by the CERN-Munich group [22] up to m S 1.4 GeV, and interpolates to m S 2.5 GeV where analytical results are expected to be valid [23]. For comparison, the spectator model uses perturbative results up to the c-quark threshold. In this case, the decay rate into pions is given by low-energy theorems and the kaon or η meson contributions are estimated by rescaling the muon branching ratio appropriately [2,24,25]. The S lifetime in this mass range for the two decay models is shown in Fig. 2. (See also the very recent work [26].) Below the electron threshold, a Higgs-like particle decays to 2 photons through a loop of heavy particles. The leading order decay rate is found by summing over the massive charged particles entering the loop [27], where C is a loop function, given explicitly in Appendix A. In this prescription, the light quark degrees of freedom are incorporated through their explicit breaking of chiral symmetry, and the associated mass of pions, kaons and eta mesons [28], i.e. through virtual loops of pions and kaons [29]. Adding the contributions from all SM particles, for m S 2m e we find where different scenarios of (a light) + (b heavy) quarks are shown. For the case of 2(3) light flavours the pion (pion and kaon) loops are taken into account, while for 0 light flavours they are neglected. The true physical value should be close to the 3 + 3 or 2 + 4 scenarios. Since the difference in decay rate between the two cases, O(4%), is negligible for the analysis of new physics, we simply choose C = 50/27.

III. SCALAR MIXING IN THE COSMOLOGICAL THERMAL BATH
In vacuum, the relevant SM masses are generated via the Higgs mechanism and are proportional to the electroweak vacuum expectation value (vev) v. In the cosmological thermal bath, and in particular near the electroweak symmetry restoration temperature, long-range interactions are screened by the plasma. Particles effectively develop a thermal mass as a representation of this screening. The mass of a particle at a given temperature T can generally be written as [30] where m 0 is the zero-temperature mass that depends on the vev and m T is the thermal mass. Note that the vev depends on T , so that m 0 also has a temperature dependance. A simple analytic formulation of the high-T Higgs thermal mass parameter in the effective potential is given by [31] where Inserting the additional term c h 2 T 2 h 2 into the Lagrangian (3) to generate the equivalent thermal mass, which predicts an electroweak symmetry restoration at the critical temperature T c 140 GeV. 1 The Higgs thermal mass (10) applies to the SM eigenstate prior to mass diagonalization. 2 After this diagonalization of h−S mixing, we obtain a temperature-dependent mixing angle, which incorporates the leading dependence on temperature for small m S . However, this expression also signals the presence of a thermal resonance when m S ∼ m h (T ), which can arise on scanning T for m h (T ) min m S m h,0 . Lattice results indicate the Higgs thermal mass drops to m h (T ) min ∼ 15 GeV at the electroweak crossover [33].
The apparent divergence in (13) at m S = m h (T ) is resolved by thermal broadening, which amounts to replacing the factor of 1/(m 2 h − m 2 S ) with a Breit-Wigner propagator for the intermediate metastable Higgs in the rest frame of the thermal bath. This is conveniently derived by considering the thermal rate Γ S at which S approaches equilibrium, given by Γ S = −ImΠ S /E where E is the particle energy. Γ S in turn is related to the S production rate Γ prod by a Boltzmann factor, Γ S = Γ prod − Γ dest = (e E/T − 1)Γ prod . The S self energy takes the form 1 This value is O(10%) smaller than the full SM value of T SM c 160 GeV from lattice simulations [32,33]. 2 As we are considering very small mixing, the effect of S on the SM thermal masses is negligible. Similarly, the S thermal mass will be m 2 S (T ) ∼ θ 2 T 2 and thus can be neglected for this study.
where Π h is the Higgs self energy. Computing the imaginary part in the on-shell limit, and with ReΠ h = m 2 h,T , leads to allowing us to read off the thermally broadened mixing angle, where Γ h is the Higgs width, or more generically, damping rate. The zero temperature width Γ h,0 = 4.07 MeV gives a reasonable approximation for the decay rate, since for the parameter regime of interest here, these decays occur late in the cosmological evolution when the temperature is low. Notice, however, that at temperatures around the electroweak scale, the damping rate (set by interactions with top quarks and weak gauge bosons) is expected to scale as Γ h (T ) ∝ T , and Γ h (T ∼ m W ) Γ h,0 . The effect of thermal broadening at high temperatures can be relevant for the epoch of freeze-in production. A density plot of θ eff (T ) from (16) is shown in Fig. 3. As T → T c , simulations suggest that m h (T ) drops rapidly near T c to 15-20 GeV [33], potentially allowing a resonance for any m S > 15 GeV, an effect that is not well captured by our v(T )-scaling model. The potential importance of the thermal resonance will be considered in more detail in the next section.
We conclude this section by noting an apparent discontinuity in the behaviour of the mixing angle at zero temperature and at temperatures close to the phase transition, assuming for simplicity that m S is a small parameter. From (4), it follows that θ ∝ A/(λ H v), while at finite temperature θ eff (T ), from (16), scales quite differently as θ eff (T ) ∝ Av(T )/(couplings × T 2 ). As v(T ) approaches zero, these two formulae have completely different behaviour. Taken at face value, this suggests that vertices with Feynman rules proportional to v(T ) will not contribute at all in the electroweak symmetric phase. However, this is only true at tree-level and the surviving diagrams are generated at higher order in perturbation theory and do not vanish in the limit of v(T ) → 0. These higher order corrections are discussed further in App. B, and example diagrams in the higher order expansion are shown in Fig. 14. Therefore, very near the phase transition, the thermally corrected mixing angle (16) will not provide an adequate description of thermal effects, and a complete treatment of thermal loops would be necessary.

IV. COSMOLOGICAL PRODUCTION VIA FREEZE-IN
The cosmological production rate of a new species S, due to 2 → 2 interactions, is given by the Boltzmann  (16), showing the location of the thermal resonance for T < 140 GeV. The peak of the resonance defines the resonance temperature Tres as a function of mS. The behavior of m h (T ) follows (9) with the naive v(T ) model, but with an additional T -dependent contribution added to ensures that m h (T ) tracks down to the minimum value of m h (T )min ∼ 15GeV near the crossover transition, as suggested by lattice simulations [33]. A finite Higgs damping rate of 0.05T was also added for illustration.
wheres is the entropy density, Y ≡ n S /s, while Λ = represents the thermal distribution of each species and |M| 2 is the spin-summed squared amplitude. In the Maxwell-Boltzmann (MB) approximation for the freeze-in mechanism Λ → f MB T . The sum goes over various multiplicity factors, such as spin and color. For 4 different species Eq. (17) takes the form [35] where s min = Max (m 1 + m 2 ) 2 , (m 3 + m 4 ) 2 , and while σ = σ 12→3S is the standard cross section averaged over initial state degrees of freedom, while g 1(2) are the spin and color multiplicity factors for initial particles.
The total S yield is found by summing all possible 12 → 3S interactions where 1, 2 and 3 are SM particles. Production channels of the form 12 → SS are suppressed by an extra factor of θ and are neglected. Since S preferentially interacts with massive particles, we anticipate a large number of possible production channels around the electroweak scale. We classify the different channels by their asymptotic behaviour in the EW unbroken phase. According to the Goldstone boson equivalence theorem [36,37], in the v 2 /s → 0 limit the behaviour must be determined by the corresponding Goldstone bosons interactions. Expanding the Higgs doublet in the form we find that the only interactions producing S in the symmetric phase will be 2 → 2 scattering channels, Fig. 4. We can therefore categorize the S-producing interactions as follows: • QCD production, which includes all diagrams with gluons and top quarks such as tg → tS.
• Yukawa annihilation, which includes the 4 reactions contributing to t R Q L → HS.
• Compton-like scattering, which includes reactions with a quark scattering off a boson in the form of t R H → Q L S.
• Gauge boson scattering, which includes the reactions purely with electroweak bosons and Higgses.
We segment the production calculation into two regimes, first for T < T c with v/(v − v(T )) > 1, and then for T > T c , where the vacuum expectation value is negligible and the dimensionful SM couplings proportional to v vanish. T close to T c can be treated by continuity. In all instances, we compute the cross sections at tree-level, with a few phenomenological improvements justified below.
In the broken phase, we incorporate the the first thermal corrections by explicitly varying the EW vev as in (12) and treating all SM masses as vev-dependent variables and dropping the T 2 term in Eq. (9). As we will see in Sec. IV A, we expect this approximation to hold for v(T ) gT , i.e. until the temperature is high enough that the thermal masses become dominated by plasma contributions. For m h (T ), the T 2 term is retained for consistency in the definition of v(T ) and to make sure θ(T ) does not have an unphysical divergence for small m S .
In the symmetric phase, we retain the quark masses in the cross sections and promote them to thermal masses acquired from the QCD plasma [34] which affects the kinematic phase space available for interactions. The Higgs doublet components all obtain the Higgs thermal mass (10). We neglect the gauge boson transverse mass. From a finite-temperature point of view, the magnetic thermal mass of a non-abelian SU(N) gauge boson vanishes at one-loop [38]. A non-vanishing value is generated at higher order as a non-perturbative quantity m 2 T ∼ (g/3π)gT [39,40], which is sub-leading compared to the other masses.
In the intermediate regime where a full finitetemperature calculation is needed v(T ) g T ≤ T c , we extrapolate from the two limiting regimes to obtain an uncertainty band for the model. In either case, we obtain results for the relic density that are consistent to within a factor of 2, which is acceptable for the problem at hand.
Retaining the top quark Yukawa coupling y t , the electroweak couplings g, g and the Higgs self-coupling λ H as the only non-zero coupling constants, the yields from each non-vanishing production channel in the m S m h limit are compiled in Table I. In total, on including the v(T ) model, we obtain the following result for the abundance from QCD and electroweak processes: for m S well below the thermal resonance region. The quoted uncertainty band corresponds to whether or not we cut off the production at T ≤ v/g s 121 GeV or we push the extrapolation to T ≤ T c 140 GeV. Later in this section, we will also make an estimate of the precision of the Maxwell-Boltzmann approximation, which will enlarge the precision band somewhat (see Eq. (43)). Note that the total yield from QCD reactions is found to be Y (QCD) S (6.3 − 8.2) × 10 10 θ 2 , a factor of approximately 20 smaller than the value quoted by Berger et al. [8] 3 . As listed in Table I, we also find many other channels with electroweak gauge bosons that contribute at or above the level of these QCD-induced reactions. We use FeynCalc [41,42] to compute the relevant cross sections, and they are listed for completeness in App. C, primarily in the limit of large Mandelstam s. The respective emissivity as a function of temperature is shown in Figs. 5 and 6.
Thus far, we have discussed the parameter range where m S is parametrically smaller than the weak scale. In this regime, with production occurring mostly at weak scale temperatures, the prediction for Y S is approximately independent of m S . However, this will change once m S is increased to a scale comparable to the thermal masses of particles in the SM bath. In Fig. 7, we plot the resulting value of Y S (m S ) that follows from generalizing the 2 → 2 production mechanisms discussed above to finite m S . We note that once m S reaches a few tens of GeV, the 2 → 2 production channels may no longer be dominant, and other processes such as resonant oscillations (1 → 1 production), and inverse decays (2 → 1 produc- if extrapolated to the phase transition. The yield from each production category in the symmetric phase is shown as Y sym cat with the range displaying the total yield for T above Tc or extrapolated down to T gv(T ). tion), may also contribute significantly to Y S . We discuss such contributions in separate subsections below.

A. Infrared divergences
When calculating various 2 → 2 production processes using the simplified v(T ) approach, we encounter additional complications due to the infrared sensitivity of the production cross sections. There are two types of infrared divergences in the interactions that require special care, both present in the channel with the largest yield contribution, bW → tS, schematically shown in Fig. 8.
At lower temperatures, where the vacuum cross sections are clearly applicable, the emission of a soft S is enhanced by the near-on-shell t and W mediators. Removing the S emission, the inverse decay process bW → t is kinematically allowed. Since we are considering an ar- bitrarily light scalar, m S m W , m t , the 2 → 2 reactions creating the S have a kinematic cutoff s ≥ (m t + m S ) 2 , which approaches the propagator singularity 1/(s − m 2 t ) as m S → 0. This type of divergence is regulated by the finite width of the propagator. For this channel, we promote the denominator of the t quark and W boson propagators to their Breit-Wigner equivalent is consistent with our v(T ) model throughout the calculation and the SM values are Γ 0 t = 1.4 GeV and Γ 0 W = 2.1 GeV [43]. The resonances are further broadened by thermal effects. Multiple schemes for calculating cross sections with unstable particles have been proposed beyond the simple substitution (24). The basic Breit-Wigner is technically incompatible with gauge invariance and Ward identities [44][45][46], a problem that can lead to dramatic inconsistencies in the small-angle scattering away from the reso-nance [44]. Explicitly comparing the cross section with and without the substitution (24) away from the resonance, we find (25) which justifies the use of Breit-Wigner propagators in this reaction near the singular point. The second class of infrared divergences appear near T c , where v → 0. The exchange of a massless spin-1 particle in the t(u)-channel generates a well-known collinear divergence in forward (backward) scattering [47]. Therefore, with the approximation m 2 W = 1 4 g 2 v 2 (T ), when the mass of m W approaches zero, all the t-channel exchange diagrams are necessarily enhanced. Once again, thermal effects will come to the rescue and stabilize these divergences.
In the context of S freeze-in production, the total cross sections with the t or u-channel gauge boson propagators do not fall off as 1/s in the high energy limit (see App. C  for expressions). For the example shown in Fig. 8, we have At low temperatures, this asymptotic behaviour does not matter as large values of s are exponentially suppressed by the energy available in the initial particle distributions, which clearly diverges as v(T ) → 0. Conceptually, this IR divergence should be regularized in the same fash-ion as the scalar QED example. In particular, g 2 v 2 (T )/4 will receive an additional m 2 W,T (T ) ∝ g 2 T 2 temperaturedependent correction. Thus the Coulomb-like enhancement near T → T c obtained from a simple extrapolation of the vacuum cross sections, with vev-dependent masses (as in Eq. 21), signifies the breakdown of our calculation as thermal effects are not incorporated. The formal strategy to deal with collinear IR divergences in thermal field theory has been laid out by Braaten and Yuan [50] in the weak coupling limit g 1. We will not perform this full calculation, but simply use the limit gT v(T ) as the boundary of validity for our tree-level cross sections.

B. Resonant S production
As seen in Fig. 3, the mixing angle has a physical resonance when m h (T ) m S . Near resonance, the h → S oscillation may become efficient, and contribute to the overall yield Y S (m S ). Below, we are going to show that the contribution of the resonance is not important for production of very light S particles, while it can contribute significantly, starting at m S in the range of few tens of GeV. In practice, there is a significant uncertainty in the behaviour of m h (T ), and the lowest value m h (T ) can acquire as a function of temperature. Recent lattice simulations suggest that near T c the thermal Higgs mass m h (T ) drops rapidly to 15-20 GeV, in a manner reminiscent of a second-order phase transition. In principle, this allow the resonance to arise for any m S above this scale. Physically, the resonance arises when the virtual Higgs that rotates into S is allowed to go on-shell, and the corresponding mixing angle develops a Breit-Wigner form (16) associated with the Higgs thermal width.
Although our primary interest is in lower values of m S , where the resonance is not present, it is interesting to consider the enhancement associated with the thermal resonance. Since Γ h m h even after accounting for thermal broadening, we can use the narrow-width approximation (NWA) to estimate the S yield from the resonance. Taking the NWA in (16), we obtain where m h (T res ) is the temperature derivative of m h (T ) evaluated at resonance. Substituting into (18), we find the simplified integral where dn h (T, E) is the Bose-Einstein distribution for the Higgs boson. Notice that the overall damping rate Γ h for the Higgs boson drops out of this formula, and its main uncertainty is encapsulated in the value for T res and m h (T res ). Evaluating the remaining integral (and using simplified Maxwell-Boltzmann statistics in the process) we arrive at an analytic estimate for the h → S oscillation-induced abundance, where all thermal quantities are to be evaluated at T = T res , m h (T res ) = m S . The simplicity of (28) is deceiving. Depending on the assumed behaviour of m h (T ), the results can vary substantially. It is possible, however, to conclude that if one takes the most extreme behaviour, m h (T ) = (v(T )/v) × m h,vac , which should conservatively over-estimate resonant contribution, the result is still quite small for small m S . In particular, we find On the other hand, for m h 100 GeV, our results indicate that Y S,res can reach ∼ 4 × 10 11 θ 2 and become comparable, or even larger than the non-resonant contributions. Interestingly, for m S as high as 100 GeV, the uncertainty in the resonant contribution becomes smaller, due to the fact that the resonance occurs at temperatures significantly lower than the cross over temperature.) Another interesting observation is that for m S m h the actual cosmological constraints are weaker than for m S = m h . In that limit, Eq. (28) is not applicable, and one has to retain the proper thermal damping rate for the production calculation. The point is that all constraints are very asensitive to lifetime of S, and the effect of close-to-resonance mixing on the decay rate is very pronounced, leading to significant shortening of lifetime, and relaxation of the bounds despite enhanced production.
To ensure that our constraints are conservative, we focus on m S below the weak scale and do not include the resonant contribution to Y S .

C. Production via Inverse Decays at large mS
We have concentrated on 2 → 2 and resonant production modes of the S scalar, which are dominant for m S below the weak scale. However, at much larger values of m S , there is also an inverse decay, or 2 → 1 type production channel, that we comment on briefly in this subsection.
The treatment is simplest when the mass of S is asymptotically larger than the weak scale. The decay is then predominantly to longitudinal W W , ZZ and to hh pairs, or equivalently into four pairs of real scalars. The total width is Production is governed by the very same width, and appropriately modifying previous results for dark photons [7], we obtain the corresponding estimate for the yield, Parametrically, this result scales as θ 2 M Pl m S /v 2 or M Pl A 2 /m 3 S , where M Pl is the Planck mass, while numerically we find, Y S,ID (m S = 1 TeV) 2.5 × 10 12 θ 2 .
Notice the slightly larger overall numerical coefficient, which results from less phase space suppression for the inverse decay process compared to 2 → 2 processes. Jumping forward to consider potential cosmological sensitivity in this high m S regime, we note that the best chance of constraining the model is provided by BBN (as the lifetime is too short for other probes). Normalizing the decay width to the most sensitive lifetime window, Γ S = 1/(1000 s), we have In spite of the large abundance, this falls about one order of magnitude short of the current best sensitivity [84].
In practice, the sensitivity to singlet scalars is enhanced below their decay threshold to weak bosons, m S < 2m W , where the decay width is set by the Yukawa coupling of b-quarks. Correspondingly, the same lifetime is achieved through a parametrically larger value of θ, which also translates into a larger abundance, and as a consequence, tighter BBN constraints. For m S ∼ 100 GeV, the inverse decay processes is subject to significant uncertainties. In particular, it is not entirely clear what the true kinematic threshold is for W W, ZZ, hh → S production via the 2 → 1 mechanism. For example, the thermal mass of longitudinal W 's is expected to be of order m W,L ∼ (0.5 − 0.6) × T at temperatures around the electroweak scale. Thus, for m S ∼ 100 GeV, it is difficult to determine for how long the W W → S process is kinematically accessible, which renders predictions for the inverse decay processes very uncertain, and sensitive to the details of thermal physics near the electroweak cross over.
Bearing these uncertainties in mind, and given our focus on the low m S range, to be conservative we will retain only the 2 → 2 production channels in analyzing the constraints below. Note that an in-depth analysis would be required to achieve a higher precision calculation of the resonant and inverse decay production channels of O(100 GeV) singlet scalar bosons.

D. Thermalization of the S sector with the SM
The freeze-in abundance of S applies for mixing angles sufficiently small that the production rate remains below the Hubble rate, Summing all the production channels, we find that the freeze-in relic abundance obtained with Eq. (18) is valid for θ θ therm ∼ 10 −6 .
Larger mixing angles ensure complete thermalization with the SM bath before S decouples and the relic abundance is simply given by the standard freeze-out paradigm. In this case, Y S is maintained at its relativistic equilibrium value until m S becomes nonrelativistic, T m S , or the coannihilation rate becomes inefficient and S decouples with its freeze-out abundance. Since S interacts dominantly with heavy particles, S remains relativistic while the coannihilating partners become nonrelativistic and the annihilation efficiency is exponentially lowered by the phasespace suppression of the other particles. Thus, S freezes out according to Eq. (36) and the abundance only depends on the number of relativistic degrees of freedom g at the decoupling temperature. Above the QCD confinement scale T QCD ∼ 200 MeV, g varies by at most a factor of 2, in the range of g ∼ 205/4 − 427/4. A conservative estimate for the thermalized S relic abundance is therefore

E. Validity of the Maxwell-Boltzmann approximation
In the classic case of WIMP freeze-out, the decoupling temperature of the species is typically in the nonrelativistic regime T decoupl ∼ m/20. The statistical ensemble of particle energies is well described by the Maxwell-Boltzmann (MB) distribution, which allows for an analytical simplification of the phase-space integrals in the Boltzmann equation. In the freeze-in scenario considered here, this simplification is not necessarily justified and we must verify its validity. We derive in App. D the analytical 3-dimensional expression to be numerically integrated for the S abundance including the correct statistical distribution for all particles.
Instead of proceeding with the full treatment, we can verify the MB approximation with a simpler integration. Keeping the exact statistical distributions in the Boltzmann equation, we obtain [35] with where the initial energies were rewritten in terms of

or Bose-Einstein (BE) distribution of species i and the + (−) is chosen for bosons (fermions) in the last term.
The MB approximation (18) arises as an analytic solution in the MB limit f 1,2 = e −E1,2/T and (1 ± f 3 ) → 1. We should stress that Eq. (39) is not mathematically correct as E 3 should have been included in the cross section phase-space integration over the end products. This integration is in general non-trivial and includes an additional angular dependence with s (see App. D). Nonetheless, we can use Eq. (39) as an estimate to the full result.
To obtain the first correction beyond the MB approximation, we can expand It is important to notice that the first order correction for the initial particles is equal to the MB limit of the (1 ± f ) term. At first order in e −Ei/T , we have where κ i = ±1, with +(−) for bosons (fermions). As expected, the bosonic distribution enhances the overall yield, while the fermionic distribution decreases it. In principle, E 3 is a function of √ s and the angular kinematics for the final state particles. As mentioned, the 1 ± f 3 factor should be included in the annihilation cross section, modifying σ. However, we know by conservation of energy that E 1 + E 2 = E 3 + E 4 , with the following bounds on E 3 To estimate the range of possible yield values from the first correction to the full quantum distribution, we can integrate Eq. (39) with Eq. (41) for each of the E 3 extremum values. The potential spread in total emissivity in shown in Fig. 9 with an estimated range, The total error for the MB approximation is thus expected to be within a factor of 2. The first order correction band in the symmetric phase lies completely below the MB value, because the top quark thermal masses in the QCD plasma dominate and suppress the available Fermi-Dirac distribution phase space.

V. COSMOLOGICAL CONSTRAINTS
Having determined the freeze-in abundance for small mixing angles, we can now place the minimal set of bounds on the S parameter space with θ θ therm ∼ 10 −6 . We update and improve the cosmological constraints partially presented in both Ref. [8] and Ref. [15] and the final results for the low-θ parameter space are shown in Fig. 1. The cosmological constraints that depend on Y S are discussed in the following subsections. Below m S 5 keV, the strongest constraint on S comes from stellar energy loss [51] and even lighter scalars in the sub-eV range are constrained by 5 th force experiments [52][53][54][55] (as discussed in [9]). We also show the projected sensitivity from the SHiP beam dump experiment [2] and an order-of-magnitude estimate of supernova energy loss [56] (which should be modified to account for in-medium effects [51]). Above the pion threshold, we find strong sensitivity to the S decay model. The colored exclusion regions presented in Fig 1 utilize the baseline decay model, which has a decay width which is larger than or equal to that in the spectator model. The baseline model therefore provides more conservative results due to the reduced abundance for a fixed lifetime.

A. Diffuse X-ray background
Many present-day satellites observe the galactic and extra-galactic photon spectrum in various wavelength bands and provide upper bounds on the luminosity of X-ray or Gamma-ray emission. These bounds apply for example to photons from decaying or annihilating dark matter. Below m S = 2m e and at small mixing angle, the lifetime of S is longer than the age of the Universe, and therefore the model is constrained by these observations. In particular, Ref. [57] derived the lifetime constraint on scalar dark matter particles decaying into 2 photons in the 4 keV < m S < 10 GeV mass range assuming τ S τ universe . We directly rescale their constraint from the HEAO-1 [58] and INTEGRAL [59] satellites to obtain an exclusion band for 4 keV < m S < 1 MeV with 10 16 sec τ S 10 22 sec displayed as X-Ray in Fig. 1.

B. CMB Anisotropies
Precision measurements of the temperature and polarization anisotropies in the CMB by the WMAP [60] and Planck [61] satellites provide strong constraints on energy injection that can ionize cosmic neutral hydrogen after recombination [62][63][64][65][66][67]. The raised ionization fraction at lower redshifts allows for delayed photon interactions, which modifies the visibility function that weighs the probability of last scattering for a given CMB photon at a specific time. This effectively damps the high-l tail of the TT power spectrum and increases the low-l E-mode polarization [62,63].
At redshift z dep , the efficiency of energy deposition in the cosmic plasma by an energetic electron-positron pair or photons injected at an earlier redshift z inj > z dep has been determined in Ref. [68] and updated in Ref. [69]. This update provides the energy fractions that go into ionization, excitations, heating and emission of low-energy photons. Given the process-dependent and z-dependent ionization efficiency, comparing the modified power spectra to the CMB data is computationally intensive. In practice, principal component analysis of modified recombination histories shows that a decaying particle is well described by a constant deposition efficiency taken at z dep = 300 [69,70]. We can then simply utilize the derived constraints for decaying particles in Ref. [7] and translate to the current model with where f eff = f (z = 300) is the ratio of energy absorbed leading to ionization over energy emitted at z dep = 300. In the mass range where S → µ + µ − is the main decay channel, we solve for f eff by integrating over the electron from muon decay, which decreases the ionization efficiency by a factor of 3 due to neutrinos radiating away energy. We repeat the procedure for the decay chains S → π + π − → µ + µ − νν, S → π 0 π 0 → γγγγ and find that it is well approximated by evaluating the decay products at their average energy from the decay. We evaluate the efficiency of kaons by weighting the main branching ratios and the decay products by their average energy, percolating down to their final e ± − γ − ν spectra. Above the di-charm threshold, the light quarks, charm quark and gluon all have similar deposition efficiencies that lie somewhere between those of electrons and muons [67].
In general, for m S 10 GeV, the efficiency tends to approach the muonic case [71]. We adopt the same ionization efficiency as muons for conservative results. The overall f eff for S with τ S = 10 14 sec is shown in Fig. 10 for both the baseline and spectator decay models. We will not extrapolate the CMB constraints down to lifetimes τ S < 10 13 s because f eff is not numerically stable   Fig. 13). The solid lines and shaded areas represent the parameters excluded in the baseline decay model. Dashed lines refer to the projected sensitivity to spectral distortions of a PIXIE-like detector [72] and the thin gray line exhibits the would-be PIXIE sensitivity in the spectator S decay model.
for decays before recombination [69] and the on-the-spot approximation at z dep = 300 fails to represent the correct physics for short lifetimes [70]. The excluded band is shown in orange within Fig. 11 for the baseline model, with the would-be exclusion region for the spectator decay model delimited by a thin gray line.

C. CMB Spectral distortions
While energy injection after recombination may be observed in the CMB as variations in the anisotropies, ear-lier energy injection can induce spectral distortion of the blackbody distribution (see Ref. [73] for a recent review), and can be used as probes of decaying particles [74,75].
Cosmological thermalization is very efficient at arbitrarily early times and the photon plasma becomes susceptible to incomplete re-equilibration of its spectrum for energy injected at z z µ 2 × 10 6 . The CMB photons are still efficient at redistributing their energy across the energy spectrum, but double-Compton scattering and Bremsstrahlung interactions that adjust the number of photons become inefficient. The bath thus develops a non-zero chemical potential in its high-energy tail resulting in the µ-distortion. At lower redshifts, z z µy 5×10 4 , Compton scattering between electrons and photons fails to maintain both species at a common temperature. The photon bath inherits a reduced temperature at low energies while high frequencies receive a relative gain in temperature, a phenomenon called the Compton y-distortion [73].
Distortion due to arbitrary energy injection can be approximated by [76] y = 1 4 where z rec = 1000 and the normalized injected electromagnetic energy is In this expression, dE dtdV is the total energy injected and Br em is the branching ratio to electromagnetic end products. In a radiation-dominated Universe, the y-distortion can be evaluated analytically, with the current entropy density s 0 = 2891 cm −3 , the current photon energy density ρ γ0 = 0.26 eV cm −3 , a time normalization of t 0 = 2.4 × 10 19 sec and where the integral I(Γ) is defined as The measured bounds from COBE/FIRAS [77] and the projected sensitivity from a PIXIE-like detector [72] are COBE/FIRAS: |y| ≤ 1.5 × 10 −5 |µ| ≤ 9 × 10 −5 , PIXIE: |y| ≤ 2 × 10 −9 |µ| ≤ 1 × 10 −8 . We approximate the electromagnetic branching ratio by weighting the average energy carried by end products from initial decays at rest with their respective branching ratios from S. Since the averaged energy carried away by the electron in a muon decay is E e /m µ = 0.35 [11], we have Br S→µ + µ − em = 0.35. The electromagnetic fractions for heavier decay products can be found from the corresponding fractions of their lighter decay products. We find the following electromagnetic energy injection ratios where we neglected the kaon decay channels that contribute less than 10% of the kaon decay width. The total electromagnetic branching ratios for the S decay channels are where 2/3 (1/3) of pions are charged (neutral), 1/2 (1/4+1/4) of kaons are charged (neutral short + long) and we assume the electromagnetic yield of high energy quarks and gluons of 0.45 is maintained for c-quarks. The total Br em as a function of m S is shown for the baseline and the spectator decay models in Fig. 12. The excluded regions from COBE/FIRAS are shown in Fig. 11, with a robust conservative overlap between all probes in the 1 MeV < m S < 2m µ mass range. A PIXIE-like detector would not change the constraints in the m S < 2m µ mass range, but has the potential to increase the sensitivity range to m S 8 GeV, with a sensitivity band that somewhat depends on the S decay model.

D. Relativistic degrees of freedom (Neff )
The total relativistic energy density in the Universe at decoupling is well constrained by the CMB. The temperature of the photon bath determines its contribution to the total radiation energy density. Any additional component is parametrized in N eff , the effective number of neutrinos with a temperature T ν = (4/11) 1/3 T γ . The Planck collaboration measurement of N eff = 3.04±0.33 [61] at 2σ is in agreement with the SM predicted value of 3.046 [78].
The constraints on early injection τ S < 1 sec of S → {e + e − , γγ} or S → µ + µ − have been derived in Ref. [21] including neutrino decoupling effects, which we apply using the freeze-out abundance from Eq. (37) in the mass range 10 MeV m S < 2m µ . For lower masses, the S lifetime at the border of the exclusion band is longer than τ S = 100 s, the maximal range derived in Ref. [21]. For longer lifetimes, we simply compare the energy density of the S sector with the SM energy densities in the neutrino and EM baths prior to the energy release (at t = xτ S ) and assume an instantaneous decay. N eff can then be estimated by comparing the relative energy density in the neutrino and EM baths where ρ S γ = Br em ρ S (t = xτ S ), ρ S ν = (1 − Br em )ρ S (t = xτ S ) and the electromagnetic/neutrino energy partition of the end products Br em is taken from Fig. 12. We choose x = 1/10 to match the constraints for τ S < 100 s. This choice of x = 1/10 is conservative. The nonrelativistic S energy density decreases less with time than the relativistic SM energy density and larger values of x would yield stronger constraints. Using the 2σ range for N eff = 3.04 ± 0.33, we find the constraints labelled N eff in Figs. 1 and 11.

E. Big Bang Nucleosynthesis
The synthesis of light nuclei during BBN is well understood, with the final abundance of stable light elements in good agreement with predicted values (see for example [79][80][81] for recent overviews and discussions of the discrepancy with 7 Li). The concordance between predicted and observed abundances of 4 He, 3 He and D can be used to constrain electromagnetic and hadronic energy injection (see Refs. [82][83][84] for reviews).
The impact of decaying particles in the BBN era depends on the ability of the decay products to efficiently interact with light nuclei, which varies as the BBN reaction network evolves and the universe cools down. Early mesons decays were thoroughly discussed in Ref. [21], effectively increasing the n/p ratio before they freeze out, thus raising the 4 He yield above its observational limit. After most neutrons have converted to 4 He, the negatively charged mesons, π − and K − can dissociate the copious 4 He, producing lighter 3 He, T, D, n and p that are fed back into the reaction network [11]. This mechanism was suggested to decrease the 7 Li prediction by reducing the amount of 7 Be and resolving the discrepancy with observations [11], but incidentally also raising the D/H ratio above 3×10 −5 from inefficient D burning, which is now excluded by observations [85,86]. Beyond τ S 10 4 s, the mesonic interaction rate with ambient nuclei is suppressed below the decay rate by the dilution due to expansion. The mesons instead have enough survival time to decay away in a shower of electromagnetic energy. Photodissociation of nuclei becomes efficient when photons have cascaded down below the e + e − thermal pair creation energy E th = m 2 e /22T . These γ-rays can photodissociate D with a binding energy of E bind D = 2.2 MeV at t ∼ 5 × 10 4 s and similarly for 4 He with E bind 4 He = 19.8 MeV at t ∼ 4 × 10 6 s [87]. We implement mesonic decays of S → ππ and S → KK (charged and neutral) by weighting the freeze-in abundance by their respective ratios. We can then apply the constraint on early decays from Ref. [21], constraining the S parameter space by an overproduction of 4 He shown in blue in Fig. 13. Moreover, we use the D/H < 3 × 10 −5 limit from Ref. [11], utilizing their stopped pion and kaon analysis for conservative results, displayed as the orange region on the S parameter space. Longer lifetimes with electromagnetic shower constraints are shown in green. We weight the S abundance by the electromagnetic branching ratio from Fig. 12 and compare the value with the electromagnetic injection limit from Ref. [84]. The upper protruding band is from a decreased D/H ratio while the lower green region is from an increased 3 He/D ratio. We have assumed 100% decays to kaons in the uncertain region, from m S ∼ 1.4 GeV to the di-charm threshold. Note that the BBN constraint on S from mesonic decays has large uncertainties due to the unknown decay spectrum. This is illustrated by the rather different exclusion region if the spectator decay model is assumed instead, as shown by the thin gray line in Fig. 13.
Above the di-charm threshold, the S decay products are well understood and the model allows for more accurate predictions. BBN is sensitive to the total energy injected via quark pairs, with a negligible dependence on the quark flavor [84]. The sensitivity is mostly dominated by the number of hadrons produced in the hadronic shower, with a m 0.3 S dependence on the total energy input (assuming a nonrelativistic S). Constraints from quark injection therefore increase with a smaller m S [84]. We use the bb limits from a 30 GeV initial particle from Ref. [84] for cc and bb injection, both rescaled by their respective branching ratio. The resulting exclusion is shown in red in Fig. 13, with the upper θ range constrained by 4 He overabundance and the lower θ values excluded by a D overproduction.
Electromagnetic injection below the pion threshold will also affect the BBN network, through decays to muons and electrons with a nontrivial dependence on the S mass [8]. However, the S decay rate is much weaker and requires θ values close to S thermalization in order to decay during the active BBN epoch. In similar regions of parameter space, there are also significant constraints from the limits on N eff , derived in section V D, along with the spectral distortion results from section V C. These already exclude the parameter space with τ S ∼ 0.1 − 10 13 s for 2m e < m S < 2m µ . We therefore do not need to perform a detailed analysis of the BBN constraints below the di-pion threshold and, by the same token, the potential solution to the 7 Li problem tagged in Ref [8] is also ruled out.

VI. DISCUSSION AND CONCLUSIONS
Interacting light scalar bosons are usually associated with a familar technical naturalness problem, m scalar ∼ (couplings) × Λ UV , resulting in the dilemma of justifying fine-tuning. The super-renormalizable portal offers an exception: the mass of the scalar can be natural both at tree and loop level, as long as m S is larger than the trilinear coupling A, or in terms of the mixing angle, m S > θ × (100 GeV). Only a very few particle physics experiments (notably flavour changing decays b → sS, s → dS [89]) and the 5th force experiments [9] can probe the parameter space with natural values of couplings. In contrast, cosmology within the ΛCDM framework provides a sensitive probe of this model, with couplings reaching down to θ ∼ 10 −16 , which in terms of the coupling to electron translates to g Se = θ × (m e /v) ∝ 10 −21 . (Incidentally, gravitons have approximately the same coupling strength to nonrelativistic electrons.) In this paper, we have attempted to take full advantage of the avail-able cosmological sensitivity, by improving the calculations of metastable S abundance produced through the unique super-renormalizable portal the SM has to offer: We have revisited the calculation of freeze-in production of the Higgs portal scalar and shown that electroweak gauge bosons and the t − b quarks have nonnegligible contributions to the total yield. We find the largest yield from top quark coalescence with a soft S emission, regulated by the finite top quark width. Even for m S m h , the favoured interactions with heavy particles push the bulk of the production just below the electroweak phase transition, in the 10 GeV ≤ T ≤ T c temperature range. Improvements in the precision of the calculation require the full finite-temperature quantum field theory machinery that models the phase transition and includes plasma screening effects. This incremental effort is substantial and not necessarily justified for new physics searches, especially because we estimate the current level of precision to be a factor of ∼ 2.
Control of the abundance calculation has allowed us to provide a full parameter space scan of the low mixing angle constraints, from sub-eV masses to 100 GeV. The cosmological constraints that depend on the primordial S abundance appear for m S > 10 keV. The Higgs portal scalar has well-defined decay products below the pion threshold. It allows for rigorous constraints in the 10 keV < m S < 2m µ range that fully cover 6 to 10 orders of magnitude in mixing angle. The low θ boundary is probed by the late decays seen in the diffuse X-Ray background and the CMB. In both cases, the lifetimes exceed the sensitive injection time window. This means that the number of decaying particles scales as dN/dt ∝ Y S Γ S ∝ θ 4 . We note at this point that the dependence on (coupling) 4 is a familiar one: it appears in any beam dump experiment searching for a decay within a detector in the limit where the decay length is much larger than the size of the experiment. It also appears in the "cosmological beam dump" constraints considered here, for the case of extremely long-lived particles. The change in S abundance by a factor of 2 would change the θ sensitivity by a factor of 2 1/4 , which would not be visible on the scale of the final log-log plot. For this reason, improvements in the calculation of Y S would provide a minimal gain in precision and our approximate framework with near-vacuum cross sections is justified, at least for m S < 2m µ . Similarly, the upper boundary for large θ is set by the freeze-out relic abundance, rather than the freeze-in abundance. Early cosmological energy injection has a logarithmic dependence on the S primordial metastable abundance [21]. This log dependence renders the exclusion region robust to variations of a few in Y f−o S , confirming the accuracy of our conservative estimate, Eq. (37), without the need for a relativistic freeze-out computation.
Above the pion threshold, in the 2m π < m S 2m c mass range, we have derived the CMB anisotropy, spectral distortion (plus the forecast for a PIXIE-like detector), and BBN constraints for two different decay models, with significant differences in the exclusion bands. Especially near the kaon resonance, m S ∼ 1 GeV, the sensitivity of each probe can vary by an order of magnitude in θ. Given this poorly defined feature, we stress that improvements in the determination of the S mesonic decay width are needed before finite-temperature abundance calculations are warranted to improve the accuracy of cosmological probes. For higher masses, where S predominantly decays to quarks, the large θ sensitivity boundary diverges sharply for short lifetimes in the m S Y S vs τ S plane [84]. Thus the θ constraint is insensitive to minor changes in Y S . For τ S > 10 4 s, however, the required energy stored in the dark sector, m S Y S , is almost flat as a function of increasing lifetimes. Changes in the freeze-in yield by a factor of 2 can thus be compensated by a factor of 2 1/4 in the mixing angle. We therefore conclude that cosmological constraints on S are robust to variations of Y S by a factor of a few, except where S decays primarily to mesons, where the determination of the decay width results in a larger uncertainty.
The derived constraints are relatively conservative with respect to the structure of the dark sector provided that a quartic interaction λ S S 2 H † H does not thermalize the dark sector. The only additional requirement for the existence of these constraints is that S decays visibly and not into some additional stable dark state. In this case, the freeze-in production mechanism provides the minimal metastable abundance S can have for a given lifetime. Any additional interactions with other states will increase its population, until it thermalizes with the SM. Large self-interactions can potentially dilute Y S before it decays, which would reduce the S abundance by a factor of ln(a τ S /a f−i ) [18]. Given that the freeze-in relic abundance is set at T f−i ∼ 10 GeV, this is a negligible factor of a few for early decays probed by BBN or N eff , but can potentially be more than an order of magnitude in Y S for lifetimes relevant to the CMB and X-ray measurements. Large self-interactions could therefore reduce the low θ sensitivity by a factor of a few. where β = 1/T and ω is the energy of the boson in the loop. For a massless boson, this integral is infrared divergent ∝ dω ω 2 , a well-known feature of finite-temperature corrections [34]. However, accounting for screening with the effective thermal mass m 2 h,T = c h T 2 , the regulated corrections scale as follows Γ µν ZZSm ∝ (Ag 2 / √ c h ) g µν . Notice that despite the associated loop suppression, this thermal loop correction is finite in the limit v(T ) → 0. Therefore, near the phase transition, a more complete treatment of thermal loops would be necessary. and the bosonic scattering cross sections are angle definitions (p i · p j ≡ p ij ) p 12 = E 1 E 2 − p 1 p 2 cos α, (D8) p 23 = E 2 E 3 − p 2 p 3 (cos α cos θ + sin α sin θ cos β), p 13 = E 1 E 3 − p 1 p 3 cos θ, (D9) The argument of the last delta function can be expressed as a function of β f (β) = p 2 4 − m 2 4 (D13) = ω + 2 (p 2 p 3 cos α cos θ + p 2 p 3 sin α sin θ cos β − p 1 p 2 cos α) , with Q=m 2 1 +m 2 2 +m 2 3 −m 2 4 . The β integral can evaluated using f (β)= − 2p 2 p 3 sin α sin θ sin β, forcing β → β 0 , where cos β 0 = − ω + 2 (p 2 p 3 cos α cos θ − p 1 p 2 cos α) 2p 2 p 3 sin α sin θ (D14) is found by solving f (β 0 ) = 0. There are actually two β 0 solutions given by sin β 0 = ± 1 − cos 2 β 0 . Since everything is symmetric in β (all dot products are cos β-dependent and the ∂f /∂β factor that appears in the dominator is an absolute value), we can simply use the positive root and multiply by 2. Hence, we get We can now focus on the angular integrations.
Since |M| 2 only consists of simple functions of cos α, the cos α integration can be performed straightforwardly for each process. Since a ≤ 0, the integration bounds are set by the real-valued criterion, between which the quadratic function is positive. Given the roots cos α ± = (−b ∓ √ b 2 − 4ac)/(2a), notice that we always have −1 ≤ cos α − and cos α + ≤ 1. The θ integral can be performed in a similar fashion. Requiring cos α ± to be real implies the condition b 2 − 4ac ≥ 0, and cos θ ± = − Q + 2p 2 2 + 2γ ∓ 2p 2 Q + p 2 1 + p 2 2 + p 2 3 + 2γ 2p 1 p 3 , where we use the shorthand γ = E 1 E 2 − E 1 E 3 − E 2 E 3 and we have These 2 integrals can be carried out analytically for each reaction. The final integral to be performed numerically for the emissivity is 3-dimensional, and a step function guarantees that the phase-space is physical and cos α ± is real-valued,