Forbidden frozen-in dark matter

We examine and point out the importance of a regime of dark matter production through the freeze-in mechanism that results from a large thermal correction to a decaying mediator particle mass from hot plasma in the early Universe. We show that mediator decays to dark matter that are kinematically forbidden at the usually considered ranges of low temperatures can be generically present at higher temperatures and actually dominate the overall dark matter production, thus leading to very distinct solutions from the standard case. We illustrate these features by considering a dark Higgs portal model where dark matter is produced via decays of a scalar field with a large thermal mass. We identify the resulting ranges of parameters that are consistent with the correct dark matter relic abundance and further apply current and expected future collider, cosmological, and astrophysical limits.


Introduction
Attempts to explain the presence and abundance of dark matter (DM) in the Universe often involve making various assumptions about the history of the very early Universe. The simplest and most natural one is to assume that, at high enough temperatures a DM particle is in thermal equilibrium with the plasma of Standard Model (SM) particles, which ensures that its density is given by Maxwell-Boltzman statistics. At some point in the expansion and cooling down of the Universe, DM undergoes a well-known freeze-out mechanism, which determines its subsequent population in the Universe. The freeze-out mechanism has been particularly popular because it requires a minimum amount of rather natural assumptions and, for reasonable values of parameters of specific particle candidates in the class of weakly-interacting massive particles (WIMPs), it is often able to produce the observed abundance of DM in the Universe. Furthermore, it does so in a manner that is insensitive to the condition of the Universe after inflation, thus effectively separating the high temperature regime from the one responsible for dark matter production.
However, it has long been known that, in addition to freeze-out, some other DM production mechanisms exist and can in fact play a dominant role in achieving the observed relic density. One particularly well-motivated example involves sub-eV axions that, due to their tiny interactions, are mainly produced not thermally but via the well-known misalignment mechanism; for recent reviews see, e.g., [1,2]. This mechanism was later extended to the case of ultra-light vector boson in [3,4]. Furthermore, extremely weakly interacting massive particles (usually referred to as E-WIMPs or super-WIMPs) are often predicted by many well-motivated extensions of the SM, for instance a gravitino in scenarios based on local supersymmetry (SUSY) or an axino in SUSY models of axions; see e.g., [1] for a recent review. If stable, they are potential candidates for dark matter in the Universe. However, due to their exceedingly feeble interactions, their population after inflation is negligible -assuming that their decoupling temperature is higher than the reheating temperature T R -since they never reach thermal equilibrium with the SM plasma, and the freeze-out mechanism is ineffective. Instead, they can be generated through so-called freeze-in [5] from scatterings and decays of some other particles.
A key feature of such "frozen-in" dark matter scenarios is that, while all SM particles remain in thermal equilibrium since the Universe reheats after inflation, the DM particle χ is absent in the early Universe and never reaches equilibrium with the SM plasma. Its production is mediated by some particles that typically remain in equilibrium with the plasma. Once the temperature drops below the mediator mass, DM production essentially stops and its relic density freezes-in.
In freeze-in scenarios, specific features and the final relic abundance of DM often depend on the details of a specific beyond-the-SM (BSM) model. In models with either the gravitino or axino as DM, their freeze-in production is typically dominated by non-renormalizable interactions at high temperatures in the case of scattering or at low ones in the case of decays [6,7]. On the other hand, in models where DM production involves for instance a light mediator, the low-temperature production dominates over the high-temperature one, thus separating again the physics of inflation from the one of dark matter [5,8,9].
In this article, we will consider a previously neglected case that some mediator field S -that could be a scalar, vector boson or a fermion -is not only in equilibrium with the thermal bath, but also develops a substantial thermal mass. That is, at sufficiently high temperatures the mass m S,T of the mediator deviates significantly from its "vacuum" one m S , i.e., the mass is dominated by thermal effects. Such an effect has recently been studied for instance while considering thermal photon decays [10]. The population of DM particles χ is assumed to be initially absent when the Universe reheats after inflation, but is generated by the decays of S. If at high enough temperatures the thermal mass of the mediator becomes sufficiently large, the possibility opens up that, when m S,T > 2m χ the decay S →χχ becomes allowed, while at T = 0 it was kinematically forbidden. As we will show, this opens up a new regime for DM production which we will call "forbidden frozen-in dark matter".
This kind of effect we believe was first identified for gravitino [11] and subsequently axino production [12]. The production rate of singlet fermions from the decay of scalar fields in a plasma including thermal corrections was calculated in [13] and applied to the case of right-handed neutrino DM. More recently it was also described in a more generic context in [14,15].
In this paper, we take a closer look at the "forbidden freeze-in" regime and identify its main phenomenological features. We further show that the equilibrium assumption of the mediator can be relaxed as long as S obtains a sizeable thermal correction to its mass, e.g., when it is chemically decoupled from the SM plasma, but remains in kinetic equilibrium with itself via self-scatterings. Interestingly, albeit perhaps as expected, the ensuing phenomenology is found to depend strongly on the dimension of the operator controlling the mediator decay into DM pair. For dimension-four operators, the production is dominated at low-temperature regime and peaks at m S,T ∼ 2m χ . Additionally, a striking feature is that, in this regime the relic abundance is ultimately almost insensitive of the DM mass, while the coupling responsible for DM production typically takes significantly larger values than in the standard freeze-in case. For mediator decays through higher-dimensional operators, on the other hand, the production is dominant at high temperatures and therefore depends on the reheating temperature. Furthermore, we argue that, since the forbidden freeze-in regime is a generic property, it might be worth exploring it in models of the freeze-in mechanism of DM production, e.g. [16][17][18][19][20][21][22][23][24][25][26][27].
As a specific realisation of the case presented above, we examine an explicit Higgs portal scenario, where the dark Higgs boson, kept in equilibrium with the SM fields through a quartic mixing term with the SM Higgs, can decay to a light (GeV-scale) Dirac fermion dark matter at a strongly suppressed rate. The thermal mass is predominantly generated by the dark Higgs self-coupling, enabling it to easily reach a thermal-mass dominated regime. Since Higgs portal scenarios are typically constrained by a variety of limits, we briefly review them and apply them to the considered model provide in order to identify new regions that are allowed by forbidden freeze-in. We will assume the dark Higgs boson is originally in equilibrium with the SM thermal bath. If the quartic mixing is low enough, it may also be produced through freeze-in, see e.g., [28] for more details on this setup. This scenario is similar to case III of [14], which closely resembles our model as the particle content is similar. The main difference, however, is the vacuum expectation value (VEV) structure. In our model the portal particle has a zero VEV ( S = 0) throughout the early Universe, and develops one only through its mixing with the Higgs boson after electroweak phase transition. This allows us to isolate the pure forbidden freeze-in regime without the impact of the SM Higgs VEV, thus simplifying the analysis and exploring the forbidden freeze-in independently. We additionally explore a different mass region than [14], which leads to a distinct phenomenology for the portal particle.
The paper has the following structure. In section Sec. 2 we briefly review the calculation of thermal mass of a scalar boson and proceed to describe in detail the mechanism of "forbidden freeze-in" through thermal mass effects. In Sec 3 we con-sider as an example an explicit Higgs-portal model in which the scenario can natually be realised, and briefly examine various criteria to ensure its consistent implemention. We then proceed to a full numerical study of the predicted relic density and describe various aspects of our scans and results, as well as the effect of applying relevant astrophysical and collider constraints.
2 Freeze-in with a thermally induced mass 2.1 Thermal mass in the early Universe As mentioned above, in this article we study the freeze-in production of DM via some mediator decays that are energetically allowed solely in a thermal bath. We expect this to occur in general, since frozen-in DM is usually assumed to be produced by particle species which are in thermal equilibrium with the SM plasma, and which should therefore develop a thermal mass correction [29][30][31] in the early Universe, similarly to the SM particles [32]. Moreover, it is this effective mass that allows "forbidden" decays to occur, as is the case for instance for plasmons (thermallydressed photons in a medium) that can decay to neutrinos [33].
Generally, at high temperatures applicable to the early Universe the thermal mass of a particle is proportional to the temperature. As this effect will be critical in realizing our forbidden freeze-in scenario, below we briefly review the case of a scalar mediator field S. In general, a scalar field features a self-interaction term, which implies that it does not need to interact very strongly with the rest of the plasma in order to develop a sizeable thermal mass. In the following we assume a self interaction term for S of the form The self-energy diagram, shown in Fig. 1, can then be readily evaluated at a finite temperature T , leading to the self-energy term where Π S corresponds to the corrected mass of S, i.e., m S,T 2 = m S 2 + Π S , and we have denoted β = T −1 , ω n = 2nπβ −1 , and ω 2 k = k 2 + m S 2 . The sum over n is evaluated by a standard procedure: 1 by transforming it to an integral over a complex quantity ω while introducing a function which has poles corresponding to ω n and unit residue. One obtains where we identify the first term as the T = 0 one-loop correction to m S , and the second one (denoted Π (T ) S henceforth) as the correction due to the finite temperature of the medium with f B ≡ e ω k β − 1 −1 the Bose-Einstein phase-space distribution.
The appearance of the phase-space distribution function regulates this otherwise quadratically divergent integral since it introduces a natural "cut-off" energy proportional to the temperature. The final result scales quadratically with temperature: Π (T ) S ∼ T 2 . In the high temperature limit, we can therefore neglect the m S contribution to ω k and arrive at In this limit, since the vacuum one-loop contribution is expected to be small compared to the tree-level one, we can neglect all T = 0 contributions and obtain an estimated form of the mass of S, It is well known, though, that naive perturbation theory does not work well when finite temperature effects are included (for examples see [29][30][31]). This can be seen by calculating the thermal correction using m 2 S → λ 24 T 2 , i.e., by re-summing the socalled "daisy" diagrams, where one would expect to get a correction of order at least O(λ 2 ). However, this is not the case in finite temperature calculations, since such diagrams induce correction O(λ 3/2 ), which may be important especially for larger values of the self-interaction coupling. We have explicitly checked that for λ 1 this re-summation leads to at most a 20% variation in the thermal mass. We will thus use the approximate result eq. (2.3) throughout this paper.

Freeze-in and mediator decay
We are interested in estimating the final relic density of a DM particle χ interacting extremely feebly with the Standard Model particles. The key assumption is that χ was never in thermal contact with the SM sector during the thermal history of the Universe, nor was it ever produced through some other means in the post-inflationary period, e.g., during reheating. Our assumed dominant dark matter production mechanism will be a suppressed decay of a bath particle S into a dark matter pair. More precisely, following the standard lore, we will assume the presence of a strongly suppressed decay channel with a small decay rate Γ χ (i.e., such as to unable one to overproduce or thermalise the χs). Assuming a boson mediator and neglecting Pauli blocking/Bose-Einstein enhancement factors, the Boltzmann equation governing the density of dark matter particle in an expanding universe is then (see, e.g., [34] for a complete recent treatment) given bẏ where M is the amplitude (summed over all internal degrees of freedom "idof") for the decay process (2.4), and the integration is over the standard phase space , and similarly for dΠ χ and dΠχ. Without loss of generality regarding the operator generating the decay S →χχ, we can rewrite the squared amplitude from the decay rate Γ χ as idof s where λ is the usual Källén/triangle function and J S is the spin of the mediator. It can then be shown (see, e.g., [34]) that under some general assumptions (i.e., a negligible initial number of DM particles, entropy conservation, and Maxwell-Boltzmann distributions for the plasma), the evolution of the DM yield (Y DM = nχ+nχ s ) is given by where h (g) are the relativistic degrees of freedom associated with the entropy (energy) density, 2 and K 1 (x) the modified Bessel function of the first kind. Defining x ≡ m S T and focusing for simplicity on J S = 0, the evolution of the yield becomes An important comment at this point is that, while in the standard freeze-in case Γ χ can be considered to be a number which factors out of the x dependence, this is not the case for forbidden freeze-in where the presence of a thermal mass m S (T ) needs to be accounted for. Let us first review in the rest of this section the standard freeze-in case where the thermal dependence of the mass can be neglected.
Assuming that m S > 2m χ and slowly varying relativistic degrees of freedom (which is the case for T 1 GeV), we can calculate the yield today (Y DM,0 ), by integrating from the reheating temperature (T R m S , x → 0) down until today (T 0 m S , x → ∞). 3 We then obtain the relic abundance in the form (2.12) where we evaluate g and h at the "mean" value of x during the DM production. 4 In Fig. 2a we show Y DM /Y DM,0 as a function of x for various values of m S , where we see that the production of DM essentially stops at the freeze-in temperature T FI ∼ m S 7 , as can be seen from the figure. That is, since typically S decouples at temperature T FO ≈ 20 m S (i.e., freeze-out), the calculation holds. However, if S decouples earlier than expected, the relic abundance of χ can be considerably smaller (if S decays rapidly to SM particles) or larger (if S decays predominantly to DM particles). In both cases the coupled system of Boltzmann equations describing the evolution of both S and χ has to be solved.
A different behavior is expected, however, when DM particles are produced via non-renormalizable operators, since the corresponding production rate increases with the temperature [5,36]. As an example, consider DM production via a 2 → 2 process which occurs due to a dimension-d operator. At high temperatures, all masses should be irrelevant, so the matrix element squared for the process can be written as function  of the center-of-mass energy √ŝ as which (assuming constant g and h) can be integrated from (2.13) where it is apparent that the high-temperature contributions dominate for n > 0 (i.e., d > 4). In the case of d ≤ 4, we expect DM production to be dominated at low temperatures (around m S , as denoted previously). Thus, these features should be treated in a case-by-case way, since the masses of the particles play an important role, and so the actual structure of the matrix element is needed.

Large thermal mass and forbidden freeze-in
Let us now turn to the case with a large thermal mass. In order to determine its effect on the freeze-in mechanism we shall assume for concretness that the scalar mediator mass takes the form 5 (2.14) An important consequence of eq. (2.14) is that the decay S →χχ can become kinematically allowed at large temperatures even if m S < 2m χ . This feature will determine the forbidden freeze-in regime. Before discussing this regime it is worthwhile to note that in some models dark matter production at early times is dominated not by decays but by the 2 ↔ 2 processes. 6 In such cases the thermal effects typically introduce only a correction, the significance of which is very model-dependent. 7 We will explicitly address the role of 2 ↔ 2 production processes for the Higgs portal model in Sec. 3, while for the following general discussion we restrict ourselves to the cases when they are subdominant.
As an example, let us consider the case m S = 0, i.e., when m S,T = αT . Assuming that the temperature is large enough so that at some early time m S,T > 2m χ is satisfied, our aim is to solve in this case the Boltzmann equation (2.11). Defining z as We observe two very different types of behavior depending on the dimension d of the operator that mediates the decay of S. In the case when d > 4, the right-hand side of eq. (2.16) increases with temperature and is therefore dominant at high temperatures close to the reheating temperature T R . Thermal effects in this case only provide a modification to the standard freeze-in through higher-dimensional operators, as is the case for gravitino or axino DM produced in scatterings of particles in the thermal plasma [37]. On the other hand, when d ≤ 4 most of the production takes place at temperatures around the dark matter mass. Indeed, DM production in this case increases at low temperature but stops when the decay becomes kinematically forbidden at αT = 2m χ . The production is thus dominated by temperatures close to m χ (or higher for small α).
In the following we will study in more details both cases to obtain a closed approximate form for the final dark matter aboundance when possible.

Higher-dimensional case
Let us first assume that d > 4 in which case at high temperatures the thermal mass of S dominates and we can write its decay rate in the form where again n = d − 4, and γ Sχ a dimensionless factor that depends on the nature of this operator. In the high-temperature regime where the approximation (2.17) is justified, the abundance equation becomes Since the production is dominated by the high temperature contribution, it is straightforward to integrate this equation, between z = 1 (the decays are kinematically not It is clear that, for d > 4 the dominant contribution comes from the regime of high temperatures (z R → 0). An important consequence of the thermal effects included here is the fact that two-body decays can significantly alter the predictions of the scenario mentioned before (which was akin to the so-called ultraviolet freeze-in scenario advocated, e.g., in [36]). That is, even if the decays S →χχ are allowed in the vacuum, the appearance of the thermal mass of S still plays a dominant role at high enough reheating temperature since in this case DM production is most efficient at high temperatures. Furthermore, comparing eq. (2.13) with eq. (2.19), we can see that since α < 1, the later tends to be generally less efficient. 8 Therefore, we conclude that the DM production via the forbidden freeze-in, in general, requires larger couplings in order to reproduce the observed relic abundance.

Four-or three-dimensional case
In the four (or three) dimensional case, most of the production is expected to take place at low temperatures, as can be seen from eq. (2.19) where the contribution from z = z R drops out (unless α is so small that the production happens close to the reheating temperature). More precisely, it takes place at around the time when the decay S →χχ stops. Thus, we expect the production to be dominated at time scale corresponding to the temperature at which m S,T ∼ 2m χ . This actually implies that up to an order one function, the decay rate satisfies Γ χ ∝ m χ . While it is therefore not possible to fully simplify the decay rate without specifying the details of the interaction, we can straightforwardly observe from eq. (2.18) that the abundance will be proportional to 1/m χ , thus implying that, up to order one corrections, the final relic density will be independent of the dark matter mass, as mentioned before.
As an example and in order to obtain a closed form for the final relic density, let us assume that: S is a scalar field, the dark matter candidate χ is a Dirac fermion and the Lagrangian contains an Yukawa interaction between S and χ, (2.20) The bath particle S decay width to dark matter is then given by The evolution of the yield is then given by In Fig. 2b we show the evolution Y DM /Y DM,0 of the number of DM particles as a function of z for the thermal mass case. It is similar to the standard case, apart from the point when the production stops, i.e., at m S,T = 2m χ . Assuming that the relativistic degrees of freedom do not vary rapidly during the production of the χs, we can integrate eq. (2.22) to obtain where again g and h are evaluated at z . 9 As we pointed out earlier, the relic abundance of χ becomes (mostly) independent of its mass, with any m χ dependence coming from z . This, and the suppression due to the α 4 , will result in relaxed constraints for the Yukawa coupling, with respect to the standard freeze-in, where 9 In this case, z is defined as  Ωh 2 scales predominantly linearly with the DM mass. Notice furthermore that in the case where the temperature correction never dominates (i.e. α T < m S ), the relic abundance is given by eq. (2.12) with the decay width (2.21), which is the standard freeze-in case, as expected.
Finally, let us conclude this section by presenting some numerical results in the case where both m S and m S,T play an important role as the temperature varies. In this case one has to calculate Y DM,0 by including both mass terms. That is, the evolution of Y DM as in eq. (2.7) needs to be solved, with m S,T given by eq. (2.14), numerically.
An example of typical dependence of Ωh 2 on m χ for the production of DM due to the decay of S, is shown in Fig. 3a. The two extreme cases of α = 0 (standard freezein) and m S = 0 (dominance of the thermal corrections to the mass) are shown by dashed blue and orange lines, respectively, while the exact numerical result is shown in solid grey. Notice that the transition between the two limits happens suddenly at m χ ≈ m S /2 which is where the blue line terminates since S →χχ becomes forbidden in the vacuum.
In Fig. 3b we present the Yukawa coupling y χ as a function of α that give the observed Ωh 2 for the scanned range of masses 10 MeV ≤ m S , m χ ≤ 1 TeV, hence overlapping regions between the two regimes may correspond to completely different values of the masses. We observe two distinct regimes: the region of standard freeze-in where m S > 2m χ is marked in blue, while the forbidden freeze-in region of m S < 2m χ is marked in orange. The shape of the forbidden freeze-in band in Fig. 3b is a simple consequence of the α 2 y χ dependence of Y DM,0 in eq. (2.23). As already noted in the d > 4 case, in the forbidden freeze-in regime one requires either larger self-interaction α of the mediator to generate a larger thermal mass, or a stronger interaction coupling between DM and the mediator, since the DM production is not as efficient as the standard case (as also shown in Fig. 3a). An important comment is that the transition between the two regimes, which happens for m S ∼ 2m χ occurs typically in a mass range of order (2m χ − m S ) ∼ αm S , which become very narrow for small α. 10 3 Forbidden freeze-in and the Higgs portal In this section we explore an explicit realisation of the general mechanism described above. We focus on a Higgs portal model, which is an archetype for a wide class of DM models where the dark sector is connected to the visible sector by a scalar mediator mixing with the SM Higgs boson.

The model
We introduce a real scalar "dark Higgs" boson field S, which is not protected by a Z 2 symmetry and hence can decay into Standard Model fields through its mixing with the SM Higgs boson. A dark matter candidate is taken to be a Dirac fermion that couples to the dark Higgs boson through a small Yukawa coupling y χ . The corresponding part of the Lagrangian thus reads with the dark Higgs boson potential term defined as 11 where H denotes the Standard Model Higgs boson doublet. The total scalar potential 10 This corresponds to the case where the mass difference preventing the decay of S into two DM particles is of the same order as the thermal contribution to m S at the typical scale T ∼ m S . In particular, in Figure 3b the forbidden region shown in orange do not probe this tuned transition regime in details for small α. 11 Notice that several other operators can be written within our symmetries, including a trilinear coupling S 3 and Yukawa couplings to left and right components of the dark matter fermion. We will neglect the trilinear in the following and enforce an exact χ-number global symmetry to fix the latter to zero.

At low temperatures (T 160 GeV), both the Higgs and dark Higgs fields
develop a non-zero vacuum expectation value (VEV), so that H = 1 √ 2 0 h + v and S → v S + S. 12 In the limit where A v the calculation simplifies significantly and the minimization conditions for the scalar potential in term of λ H and v S can be easily obtained as Furthermore, we can rotate the scalars to their eigenvalue basis, i.e. h, S → R h, S , where R is a rotation matrix parametrised by the small angle θ given by where we have used the masses of h and S (at T = 0) defined by The branching ratio of the Higgs decay to invisible particles is constrained to be [38] smaller than 0.19, which translates to λ HS 10 −2 . Furthermore, note that while we have supposed that the trilinear term λ 3 S 3 was negligible in our original Lagrangian, 13 the shift by v S re-introduces such a term as λ S 3! v S S 3 . For consistency, we will therefore further require that this contribution is negligible with respect to µ S , leading to the condition (3.6) Notice that this also automatically ensures that the shift in the SM Higgs boson quartic coupling λ H is negligible in eq. (3.3). An interesting feature is that the dark Higgs boson is extremely long-lived at low mass. When only its decays into a lepton pair are kinematically allowed, and assuming µ S ∼ m S , we obtain (3.7) 12 Since S plays a crucial role in the production of DM before and after EW phase transition, we just denote the VEV-shifted dark Higgs boson as S in order to avoid changing the notation when dealing with different temperature regimes.  The relative importance of these processes is to large extent determined by the hierarchy of the highlighted couplings y and λ HS as well as by the mixing angle θ.
As we will see in the next section, such long lifetime are severely constrained by astrophysical limits and beam dump limits. For simplicity, we will therefore typically restrict ourselves to m S > 100 MeV in the following. 14 The relevant processes determining the evolution of number densities of S and χ in this model are: i) the direct mediator decay S →χχ, ii) the mediator decay to SM particles due to its mixing with the SM Higgs boson, and iii) the annihilation of S to SM particles, as well as all the inverse reactions. The Feynman diagrams for these processes are given in Fig. 4. 15 The direct S →χχ decay width is given by eq. (2.21) and is suppressed by the very small Yukawa coupling y χ . The decay of S to SM particles is given by Γ(S → SM) = θ 2 Γ h→SM (m S,T ), where the Γ h→SM (m S,T ) is the total width of the SM-like Higgs boson with mass m S,T . We implement using the results taken from [39][40][41] and a direct evaluation for leptonic decay at low masses.
The S annihilation cross section as a function of the Mandelstam variable s reads It can be of the order of the standard WIMP annihilation cross-section, or smaller. This is because S is unstable and therefore its number density right after freeze-out can be much larger than for standard WIMP. We will assume that either λ HS or the mixing angle θ are large enough to ensure that S was in equilibrium at very early times (see discussion in the next section).
Apart from the processes shown in Fig. 4, additional 2 ↔ 2 processes can in principle play a role in the production of χs and/or their early-time thermalization with the SM plasma. These are: SS ↔χχ, hh ↔χχ and the co-annihilation process Sh ↔χχ. The first one has s−, t− and u−channel contributions which are proportional to θ 2 λ 2 HS y 2 χ and y 4 χ , respectively. The second and third have only s−channel diagrams proportional to A 2 y 2 χ and λ 2 HS y 2 χ , respectively. It is clear that all of theses 2 ↔ 2 processes are strongly suppressed with respect to direct S decays due to phase space suppression exhibited by the 2−body phase space of the former channels. However, in the deeply forbidden regime (i.e., for very small λ S ), when the decay is kinematically allowed only at very high temperatures, all the aforementioned channels could in principle play some role in the evolution of χ. In light of this, we have implemented all of the above processes in the numerical approach presented in the next section and checked explicitly that for the parameter ranges covered by our scan these processes indeed can be safely neglected in solving the evolution equations of S and χ number densities.

Relic density and numerical study
In light of the above discussion the coupled computation of the freeze-out of S and the freeze-in of χ is performed under the assumptions that: i) χ had negligible abundance after reheating and had not reached chemical equilibrium, ii) S was in chemical equilibrium at early times and remained in kinetic equilibrium for all the temperatures relevant for the production of χs. 16 In practice, the assumption made in the numerical code is that the above conditions are satisfied up to x = 0.1, where we define x ≡ m S /T . For x < 0.1 it is assumed that S traces its equilibrium value while the evolution of χ is given by eq. (2.11), starting from the reheating temperature T R assumed to be given by x R = 10 −9 . We checked explicitly that assuming different T R does not change the result. For x > 0.1 the coupled system of the Boltzmann equations for the number densities of S and χ is numerically solved, including all the relevant processes discussed above. 17 16 Kinetic equilibrium is an extremely good assumption in the parameter space studied in this work since away from the Higgs boson resonance elastic scatterings of S off particles of the SM plasma are much more frequent than annihilations of S. In a different model where this assumption would be violated one would be required to solve also for the temperature of S or even its full phase space density, see [42]. This would also bring additional complication to the forbidden freeze-in case as the thermal mass of S would need to be computed out of equilibrium. In fact, even if S is still in kinetic equilibrium (with the SM plasma or with itself), but already chemically frozen-out, the thermal mass would not be given by eq. (2.3). However, this caveat has no implications for our results since in the studied model the forbidden freeze-in happens at large enough temperatures where S is still in equilibrium. 17 This is done to ensure that the χ production from S decay takes into account possible deviations from chemical equilibrium of S. As stated before, this does not affect the forbidden freeze-in regime in our model, but it does some part of the parameter space of the standard freeze-in. For discussion and explicit forms of suitable Boltzmann equations see e.g. [34].
Within this setup there are several possible regimes leading to the correct DM abundance. In the following we first show some representative examples of the evolution of the yields of S and χ for different regimes and then present and discuss the results of our scan of the parameter space of the model.

Evolution of number densities
In Figures 5-7 we present the yields of S and χ for some characteristic cases. In all following figures the green dashed lines correspond to Y S while the solid lines to Y DM with the blue color indicating standard (non-forbidden) regimes and the beige one forbidden regimes. For completeness, the light gray area highlights the evolution of the yields during the time before the electroweak phase transition (EWPT). In all the plots the different shadings of the lines correspond to the variation of the most relevant parameter for a given regime, as indicated in the figures.
The simplest case is the usual freeze-in, where m S > 2m χ and Y DM gradually grows, with most of the production happening around T ∼ m S . This is shown in Fig. 5a. In this case the final relic abundance of χ is insensitive to any variations in the self-coupling λ S due to the fact that the thermal effects are important only for T m S , which is a very short (in real time) period. Thus, the thermal mass of S has a very small impact on the result in the standard freeze-in regime, as expected. Additionally, note that the equilibrium number density of S is also affected only at early times due to thermal corrections, as they shift the value of m S,T .
In Fig. 5b we show a typical case of forbidden freeze-in, where an opposite behaviour can be seen. The production is active only at small x and is both stronger and terminates later for larger values of λ S . In this forbidden regime the final DM abundance is therefore very sensitive not only to value of y χ but also the self-coupling of the mediator. Another point worth stressing is that one does not need large values of λ S to get a sizable effect, so the opening of the forbidden decay due to thermal effects is in fact a generic feature of the freeze-in mechanism. Figure 6a shows a case of a transition between the standard and the forbidden regimes. For fixed m S = 100 GeV we vary m χ and see that, as expected, around the transition the result is very sensitive to precise value of the DM mass. In the forbidden regime increasing m χ further leads to only very mild change in the relic abundance, i.e., the yield Y DM is inversely proportional to m χ , in agreement with eq. (2.23). This approximate DM mass independence of the relic density is an distinct feature of the forbidden freeze-in scenario.
In Fig. 6b a slightly different mechanism is shown. It occurs when nominally this would be a standard freeze-in case with m S > 2m χ but, due to the EWPT and its effect on the mass of S (which arises when the SM Higgs gets its VEV due to the presence of the mixing quartic coupling λ HS ), there appears a temporary regime where S →χχ is not allowed and the χ production is blocked for a while. However, if the self-coupling λ S is large enough the thermal mass overcomes the suppression due to the EWPT and re-opens the decay. This is an example of a situation when the thermal mass has a large impact on the relic abundance even in the standard freeze-in regime of m S > 2m χ . A scenario like this is close to what was studied, in a more general context, in ref. [14].
Finally, in Fig. 7 we show for completeness examples of cases where the χ production is dominated by the late-time decay of S. These cases are not directly related to the main focus of this work but are present in some regions of the parameter when we scan of the full model and therefore important in their own right. In these cases the complete evolution of both S and χ is crucial. In Fig. 6a the final DM abundance is determined by the branching fraction of the S decays to χ and to SM particles which in the plot is parametrised by the value of the trilinear coupling A. For smaller values (corresponding to a weaker mixing with the SM Higgs boson), DM particles constitute a larger fraction of S-decay products. Figure 7b shows a situation where the details of the freeze-out of S strongly affect its abundance that is then transferred to the χs via (rare) decays. This also shows the potential impact that the choice of λ HS can have on the final relic abundance of DM. Note that in this plot different lines correspond to different relation between x and T due to electroweak symmetry breaking contribution to m S which depends on λ HS . Around the EWPT the T -dependence of the VEV causes a temporary regime where in the standard case the S →χχ is forbidden and χ production is blocked. However, if the self-coupling λ S is large enough the thermal mass overcomes the suppression of m S,T due to the EWPT and re-opens the decay.

Scan setup and results
A numerical scan of the model parameter space has been conducted using MultiNest [43] to direct the scan towards values of the relic density within 2σ of the standard result from the Planck Collaboration [44] Ωh 2 = 0.1198 ± 0.0012 that we set as an allowed range. 18 The private code BayesFITS, automatically created using routines from SARAH [45][46][47] is used to interface it with the Mathematica implementing the approach discussed above which we use to evaluate the relic density. The details of the parameter ranges are given in Table 1.
In Fig. 8 we show the points in the scan that satisfy the DM relic density constraint. As before, blue colour indicates the standard freeze-in regime and the beige one the forbidden regime. It is apparent that these two regimes exhibit very distinct patterns. In particular, as discussed in a previous section, the standard freeze-in is in most cases not sensitive to the value of self-coupling λ S . It also requires very low values of the Yukawa coupling; otherwise DM is overproduced. In contrast, the forbidden regime is highly sensitive to λ S , as expected. Indeed, the smaller the self-coupling, and therefore the thermal mass, the earlier the production stops and therefore the larger y χ is needed to obtain the correct relic abundance of DM. Nevertheless, it is a new, interesting regime that is generically present in our scans and  Figure 7. Examples of yields evolution when the χ production is dominated by the late time decay of S. (a) Dependence on the trilinear coupling A, which (for fixed λ HS ) governs the branching ratio of S decay to χ and to SM particles. Here the freeze-out of S proceeds as for usual WIMP, with decoupling at x ∼ 20. (b) Dependence on the portal coupling λ HS (for fixed A). Lowering λ HS leads to smaller mass, due to the EWSB contribution, and also earlier freeze-out with larger Y S which then translates to larger χ population. Note that in this plot the relation between x and time/temperature is different for different lines.

Parameter
Description additionally leads to a freeze-in DM interacting more strongly than in the usually studied scenarios. An important comment is that, while we explicitly enforce the consistency condition (3.6) for all the scan-based plots, our choice of parameters implies that most of our points with low dark Higgs boson mass exhibit also a small quartic mixing λ SH . This is a direct consequence of eq. (3.5).   Figure 9. Experimental limits for our model, for points satisfying the observed relic density at 95%CL in the plane m S − τ S for m S < 2m χ (orange) and m S > 2m χ (blue).

Experimental limits
In dark Higgs models dark matter particles are largely out-of-reach of current experiments due to their extremely small interactions with the visible sector. The mixing of the scalars h and S induces, however, interactions of S with the SM particles which are proportional to θ, hence mediating the decay of S to SM particles (if kinematically allowed). Since θ is suppressed by powers of v S /v, the dark Higgs boson S is typically long-lived, as shown in eq. (3.7) -particularly for low masses. In this case bounds from both colliders and fixed target experiments [48], and for longer life-time, from astrophysics [40] apply. Such limits have traditionally been very well-studied. We summarise them below and in Figure 9 which indicates the most relevant ones for our setups. First of all, and apart from enforcing the proper dark matter relic density, astrophysical bounds can be divided in two main categories, and typically set an upper bound on the dark Higgs boson lifetime, or equivalently a lower limit on its mixing angle with the SM Higgs boson.
• Cooling rate of the supernovae SN1987. This limit uses the fact that the core of the nova is a thermal environment with temperature T SN ∼ 30 MeV where dark Higgs bosons can be produced and -if sufficiently feebly coupled -escape the core and lead to a faster cooling of the supernova. Standard bounds for dark Higgs boson [49] are derived from the requirement that the cooling rate from dark sector particles do not exceed the neutrinos one [50][51][52]. 19 • Bounds from enforcing a successful big bang nucleosynthesis. We use the recent bounds from [40] which are derived from the same Lagrangian as in Sec. 3.1.
In the lower mass range (below the π-meson mass threshold) the dominant bounds are derived by constraining the entropy injections from the e + e − / µ + µ − decays of the dark Higgs boson. Once dark Higgs boson annihilation/decay into hadrons becomes accessible, more stringent bounds arise from preventing neutron-proton ratio to differ significantly from 1/6 ∼ 1/7 due to the p ↔ n meson-mediated interaction. Finally, for heavy enough dark Higgs boson, direct baryon/anti-baryon production become the dominant decay channel of S. The subsequent anti-baryon annihilation with the ambient proton and neutron population further modifies the proton-neutron ratio. This limit dominates above the di b-quark threshold. An important comment is that this limits depends on the dark Higgs bosons abundance Y S , however given our restriction eq. (3.6), dark Higgs bosons abundance typically freezes-out earlier than in [40] which implies that the relativistic abundance is maintained for larger masses. Altogether, modifying λ HS only changes the limits by an O(1) factor, as can be seen in [40]. This is a simple consequence of the fact that, in order to avoid a significant modification of the p/n ratio, one relies on ensuring that the dark Higgs boson decay before BBN. The limit then roughly depends on the exponentially suppressed initial abundance Y S exp(−t p/n /τ S ) where t p/n ∼ 2.6s is the freeze-out time of the proton/neutron ratio. 20 The second class of constraints arises from colliders and beam-dump experiments, and typically sets a lower bound on the dark Higgs boson life-time.
• Limits from dark Higgs boson production and decay. Based on the original ALP searches in CHARM [53], these limits have been recently updated with a better modelling of the dark Higgs boson lifetime in the challenging region of m S around 1 GeV in [41]. Note that we have included the projected limits from SHiP at 2 × 10 20 proton-on-target [54] as a long-term prospect. Similarly, and as an example of limits from LHC-based experiments, we have included a projection for FASER phase 2 at the HL-LHC from [55]. Notice that these next generation experiments have the potential to start probing the relevant parameter space.
• Precision physics in meson decays. In the lower mass range, the dominant limits arise from the meson decay K + → π + νν studied in the E949 experiment [56]. Finally, for the heavier mass range -corresponding to intermediate masses around 1 GeV -the main constraints come from searches for visible decay of B-meson by the LHCb collaboration [57]. In both cases, we use the recasted bound from [41].
Note that in the long term, several planned experiment have the potential to greatly improve the limits in this mass range [48]. LHC-based experiments, such as FASER, MATHUSLA [58] or CODEX-b [59] are particularly interesting in that the decay of Higgs boson mediated through the quartic mixing λ HS can significantly enhanced the detection prospects as they are not tied to the mixing angle per-se but only to λ SH (hence the invisible branching ratio of the SM Higgs). Saturating the limits from invisible Higgs decay then leads to orders of magnitude improvements, particularly in the case of MATHUSLA or CODEX-b [48].

Conclusion
In this article we studied the forbidden freeze-in regime. Building on a standard decay-mediated freeze-in scenario, we focused on the case where the decaying mediator field couples strongly enough to the SM thermal bath to develop a significant thermal mass at high temperature. This strongly modifies existing predictions, and in particular leads to a particularly interesting regime of forbidden freeze-in, where the decay into DM particles is kinematically forbidden in the vacuum but is allowed to proceed in the thermal bath. In Sec. 2, we described in some detail the effect of including a sizeable thermal mass of the mediator. Assuming that the main production channel of DM is the decays of a bath particle into a pair of DM particles, we showed that freeze-in can be dominant at both high and low temperatures, depending on the dimension of the operators that couple the DM to the bath particle. Although the d > 4 operators show high-temperature dominance of DM production, this is different from the standard freeze-in case at high temperatures since the dominance does not happen due to the kinematics of the production process, but due to the thermal mass of the bath particle. Comparing the forbidden with the standard case of high-temperature freeze-in, we showed that the forbidden freeze-in is generally less efficient, leading to a stronger coupling between the DM particle and the mediator. For the case of operators with d ≤ 4 we showed that the production is dominant at lower temperatures close to the DM mass. In this case the scale of DM production is insensitive to the scale of inflation and reheating, similarly to the case of standard "freeze-out". Furthermore, the relic abundance is ultimately almost insensitive to the DM mass and the coupling responsible for the DM production can take significantly larger values than in the standard freeze-in scenario.
As a concrete example we studied a scalar portal model where the DM (assumed to be a Dirac fermion) is coupled only to a scalar which in turn is coupled to the SM Higgs boson field. In Sec. 3 we showed the effect that the scalar thermal mass has on the production of DM. We studied in detail the solution of the coupled Boltzmann equations for the DM particle and the mediator and discussed various possible types of the evolution of DM relic density. We also performed a scan of the parameter space of the model at hand and presented the region where the observed relic abundance can be obtained. Focusing on the same model, in Sec. 3 we discussed its experimental search prospects. Since the DM coupling to the SM particles is expected to be extremely suppressed (due to the small Yukawa coupling and the small mixing angle between the portal and Higgs boson fields) this model can be mostly probed by searching for a long-lived scalar mediator. We showed the impact of all the relevant bounds on the parameter space, including BBN, LHCb, CHARM, as well as astrophysical bounds for the presence of a light scalar field coupled to the Higgs boson. Also, we discussed the reach of upcoming fixed-target experiments (SHiP and FASER) and showed what part of the parameter space they will be able to probe.
As we have already pointed-out, the forbidden freeze-in regime is a general feature of the freeze-in mechanism. It greatly expands the parameter space in models where otherwise the DM cannot be produced by the decays of bath particle. Therefore, the analysis performed in this work not only provides new interesting viable regions of the Higgs portal model but may also bring some insight into how the forbidden freeze-in works in general. Our results also strongly suggest that it would be interesting to re-examine the dark matter abundance in other types of freeze-in models in order to uncover their respective forbidden freeze-in regimes.