Forbidden frozen-in dark matter

TL;DR

This paper explores a novel dark matter production mechanism via freeze-in, emphasizing the role of large thermal corrections to mediator masses at high temperatures, which can dominate dark matter generation.

Contribution

It introduces the significance of thermally corrected mediator decays in dark matter production, expanding the standard freeze-in paradigm with new parameter space considerations.

Findings

Thermal corrections can enable forbidden mediator decays to dominate dark matter production.

The model identifies parameter ranges consistent with observed dark matter abundance.

Constraints from collider, cosmological, and astrophysical data are applied to the new mechanism.

Abstract

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.

Tables

Table 1. Table 1 : Ranges of the parameters of the model analysed in this scan. Dimensionful quantities are given in GeV.

Parameter	Description	Range	Prior
${\mu }_{\chi }$ (GeV)	Dark matter Lagrangian mass	$\mathrm{\hspace{0.17em}0.005},\mathrm{\hspace{0.17em}50}$	Log
${\mu }_S$ (GeV$)$	Dark Higgs boson Lagrangian mass	$\mathrm{\hspace{0.17em}0.100},\mathrm{\hspace{0.17em}50}$	Log
$A$ (GeV)	Trilinear mixing	${\mathrm{\hspace{0.17em}10}}^{-8},{\mathrm{\hspace{0.17em}10}}^{-2}$	Log
${\lambda }_{HS}$	Quartic mixing	${\mathrm{\hspace{0.17em}10}}^{-8},{\mathrm{\hspace{0.17em}10}}^{-2}$	Log
$y_{\chi }$	Dark matter Yukawa	${\mathrm{\hspace{0.17em}10}}^{-14},{\mathrm{\hspace{0.17em}10}}^{-8}$	Log
${\lambda }_S$	Dark Higgs self-coupling	${\mathrm{\hspace{0.17em}10}}^{-4},\mathrm{\hspace{0.17em}1}$	Log

Equations

$${\mathcal{L}}_S=-\frac{{\lambda }_S}{4!}{S}^4.$$

$${\Pi }_S=\frac{{\lambda }_S}{2\beta }\sum _{n=-\infty }^{\infty }\int \frac{d^3\vec{k}}{(2\pi {)}^3}\frac{1}{{\omega }_n^2-{\omega }_k^2},$$

$${\Pi }_S=i\frac{\lambda }{2}\int \frac{d^{4}k}{(2\pi {)}^4}\frac{1}{k^2-m_S{}^2+iϵ}+\frac{\lambda }{2}\int \frac{d^3\vec{k}}{(2\pi {)}^3}\frac{f_B({\omega }_k)}{{\omega }_k},$$

$${\Pi }_S^{(T)}=\frac{\lambda }{24}{T}^2.$$

$$m_{S,T}{}^2\approx {\Pi }_S^{(T)}=\frac{{\lambda }_S}{24}{T}^2.$$

$$S\to \bar{\chi }\chi ,$$

$${\dot{n}}_{\chi }+3H{n}_{\chi }=\int d{\Pi }_{S}d{\Pi }_{\chi }d{\Pi }_{\bar{\chi }}\times \frac{1}{e^{E_S/T}-1}\times (2\pi {)}^4{\delta }^4(P_S-P_{\chi }-P_{\bar{\chi }})\sum _{\mathrm{ido}{\mathrm{f}}^{\prime }\mathrm{s}}|\mathcal{M}{|}^2,$$

$$\sum _{\mathrm{ido}{\mathrm{f}}^{\prime }\mathrm{s}}|\mathcal{M}{|}^2=(2J_S+1)\frac{8\pi m_S^3}{\sqrt{\lambda (m_S^2,m_{\chi }^2,m_{\chi }^2)}}{\Gamma }_{\chi },$$

$$-HsT{\delta }_h^{-1}\frac{d{Y}_{\mathrm{DM}}}{dT}=\frac{(2J_S+1){\Gamma }_{\chi }}{{\pi }^2}{K}_1(m_S/T)m_S^{2}T.$$

$$s\equiv$$

$${\delta }_h\equiv$$

$$H=$$

$$\frac{d{Y}_{\mathrm{DM}}}{dx}=\left(\frac{{\Gamma }_{\chi }(m_S,m_{\chi })}{5.93\times 10^{-19}\mathrm{GeV}}\right){\left(\frac{1\mathrm{GeV}}{m_S}\right)}^2\frac{K_1(x)x^3}{\sqrt{g}h}{\delta }_h.$$

$$\Omega h^{2}\approx 2.8\times 10^{8}\,\dfrac{m_{\chi}}{{\rm GeV}}Y_{\rm DM,0}{}\approx\left(\dfrac{\Gamma_{\chi}(m_{S},m_{\chi})}{4.5\times 10^{-28}\,{\rm GeV}}\right)\left(\dfrac{1\ {\rm GeV}}{m_{S}}\right)^{2}m_{\chi}\left(\dfrac{1}{\sqrt{g}\;h}\right)\!\Big{|}_{x=\langle x\rangle}\,,$$

$$⟨x⟩\equiv \frac{{\int }_0^{\infty }dx{x}^3{K}_1(x)\times x}{{\int }_0^{\infty }dx{x}^3{K}_1(x)}\approx 3.4.$$

$$|\mathcal{M}{|}^2\approx {\gamma }_d{\left(\frac{\sqrt{\hat{s}}}{\Lambda }\right)}^{2n},$$

$$\frac{d{Y}_{\mathrm{DM}}}{dx}\approx \frac{1}{512{\pi }^5}\frac{{\delta }_h}{Hsx}\frac{m_S}{x}{\int }_0^{\infty }d\hat{s}\frac{{\hat{s}}^{n+1/2}}{{\Lambda }^{2n}}{K}_1\left(\frac{\sqrt{\hat{s}}}{T}x\right),$$

$$Y_{\rm DM,0}{}\approx\dfrac{x_{\rm R}^{1-2n}-x_{0}^{1-2n}}{2n-1}\left(\dfrac{4^{n}\,n!\,(n+1)!\,\gamma_{d}}{2.34\times 10^{-15}}\right)\left(\dfrac{m_{S}}{\Lambda}\right)^{2n}\left(\dfrac{1\ {\rm GeV}}{m_{S}}\right)\left(\dfrac{1}{\sqrt{g}\;h}\right)\!\Big{|}_{x\sim x_{\rm R}}\,,$$

$$m_{S,T}{}^2\approx m_S{}^2+{\alpha }^2{T}^2.$$

$$z\equiv \frac{2m_{\chi }}{\alpha T},$$

$$\frac{d{Y}_{\mathrm{DM}}}{dz}=\left(\frac{{\Gamma }_{\chi }(m_S,m_{\chi })}{5.93\times 10^{-19}\mathrm{GeV}}\right){\left(\frac{1\mathrm{GeV}}{2m_{\chi }}\right)}^2\frac{{\alpha }^4{K}_1(\alpha )}{\sqrt{g}h}{\delta }_{h}z.$$

$${\Gamma }_{\chi }\sim \frac{{\gamma }_{S\chi }}{16\pi }{m}_S{\left(\frac{m_S}{\Lambda }\right)}^{2n}=\frac{{\gamma }_{S\chi }}{16\pi }{\alpha }^{2n+1}{\left(\frac{T}{\Lambda }\right)}^{2n}T,$$

$$\frac{d{Y}_{\mathrm{DM}}}{dz}=\left(\frac{{\gamma }_{S\chi }}{2.96\times 10^{-17}}\right){\left(\frac{2m_{\chi }}{\Lambda }\right)}^{2n}\left(\frac{1\mathrm{GeV}}{2m_{\chi }}\right)\frac{{\alpha }^4{K}_1(\alpha )}{\sqrt{g}h}{\delta }_h{z}^{-2n}.$$

$$Y_{\rm DM,0}{}=\frac{z_{\rm R}^{1-2n}-1}{2n-1}\left(\dfrac{\alpha^{4}\;K_{1}(\alpha)\;\gamma_{S\chi}}{2.96\times 10^{-17}}\right)\left(\dfrac{2m_{\chi}}{\Lambda}\right)^{2n}\left(\dfrac{1\ {\rm GeV}}{2m_{\chi}}\right)\left(\dfrac{1}{\sqrt{g}\;h}\right)\!\Big{|}_{z\sim z_{\rm R}}\,.$$

$${\mathcal{L}}_{\mathrm{int}}=-y_{\chi }\bar{\chi }\chi S.$$

$${\Gamma }_{S\to \bar{\chi }\chi }=\frac{y_{\chi }^2}{8\pi }\frac{{\left(m_S^2-4m_{\chi }^2\right)}^{3/2}}{m_S^2}.$$

$$\frac{d{Y}_{\mathrm{DM}}}{dz}={\left(\frac{{\alpha }^2{y}_{\chi }}{3.86\times 10^{-9}}\right)}^2\left(\frac{1\mathrm{GeV}}{2m_{\chi }}\right)K_1(\alpha )\frac{{\left(1-z^2\right)}^{3/2}}{\sqrt{g}h}{\delta }_h.$$

$$Y_{\mathrm{DM},0}={\left(\frac{{\alpha }^2{y}_{\chi }}{5\times 10^{-9}}\right)}^2\left(\frac{1\mathrm{GeV}}{2m_{\chi }}\right)K_1(\alpha ){\left(\frac{1}{\sqrt{g}h}\right)}_{z=⟨z⟩},$$

$$⟨z⟩\equiv \frac{{\int }_0^{1}dz{\left(1-z^2\right)}^{3/2}\times z}{{\int }_0^{1}dz{\left(1-z^2\right)}^{3/2}}\approx 0.34.$$

$${\mathcal{L}}^{\mathrm{DM}}$$

$$V_{HS}=\frac{{\mu }_S^2}{2}{S}^2+\frac{{\lambda }_S}{4!}{S}^4+AS{H}^{†}H+{\lambda }_{HS}{S}^2{H}^{†}H,$$

Full text

1 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., Baer:2014eja ; Marsh:2015xka . This mechanism was later extended to the case of ultra-light vector boson in Nelson:2011sf ; Arias:2012az .

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., Baer:2014eja 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_{\rm 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 Hall:2009bx 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 $\chi$ 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 Ellis:1984eq ; Covi:2001nw . 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 McDonald:2001vt ; Hall:2009bx ; Blennow:2013jba .

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 Dvorkin:2019zdi . The population of DM particles $\chi$ 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_{\chi}$ the decay $S\to\bar{\chi}\chi$ 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 Rychkov:2007uq and subsequently axino production Strumia:2010aa . The production rate of singlet fermions from the decay of scalar fields in a plasma including thermal corrections was calculated in Drewes:2015eoa and applied to the case of right-handed neutrino DM. More recently it was also described in a more generic context in Baker:2017zwx ; Bian:2018mkl .

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}{}\sim 2m_{\chi}$. 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. Aoki:2015nza ; Ayazi:2015jij ; Shakya:2015xnx ; Tsao:2017vtn ; Dedes:2017shn ; Bae:2017dpt ; Duch:2017khv ; Biswas:2018aib ; Bhattacharyya:2018evo ; Goudelis:2018xqi ; Belanger:2018sti ; Abdallah:2019svm .

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., Heeba:2018wtf for more details on this setup. This scenario is similar to case III of Baker:2017zwx , 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 ($\langle S\rangle=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 Baker:2017zwx , 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 consider 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 Das:1997gg ; Kapusta:2006pm ; Bellac:2011kqa in the early Universe, similarly to the SM particles Fodor:1994bs . Moreover, it is this effective mass that allows “forbidden” decays to occur, as is the case for instance for plasmons (thermally-dressed photons in a medium) that can decay to neutrinos Braaten:1993jw .

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

[TABLE]

The self-energy diagram, shown in Fig. 1, can then be readily evaluated at a finite temperature $T$, leading to the self-energy term

[TABLE]

where $\Pi_{S}$ corresponds to the corrected mass of $S$, i.e., $m_{S,T}{}^{2}=m_{S}{}^{2}+\Pi_{S}$, and we have denoted $\beta=T^{-1}$, $\omega_{n}=2n\pi\beta^{-1}$, and $\omega_{k}^{2}=\vec{k}^{2}+m_{S}{}^{2}$.

The sum over $n$ is evaluated by a standard procedure:111Details can be found in the literature, e.g., Das:1997gg ; Kapusta:2006pm ; Bellac:2011kqa . by transforming it to an integral over a complex quantity $\omega$ while introducing a function which has poles corresponding to $\omega_{n}$ and unit residue. One obtains

[TABLE]

where we identify the first term as the $T=0$ one-loop correction to $m_{S}$, and the second one (denoted $\Pi_{S}^{(T)}$ henceforth) as the correction due to the finite temperature of the medium with $f_{B}~{}\equiv~{}\left(e^{\omega_{k}\beta}-1\right)^{-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: $\Pi_{S}^{(T)}\sim T^{2}$. In the high temperature limit, we can therefore neglect the $m_{S}{}$ contribution to $\omega_{k}$ and arrive at

[TABLE]

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$,

[TABLE]

It is well known, though, that naive perturbation theory does not work well when finite temperature effects are included (for examples see Das:1997gg ; Kapusta:2006pm ; Bellac:2011kqa ). This can be seen by calculating the thermal correction using $m_{S}^{2}\to\dfrac{\lambda}{24}T^{2}$, i.e., by re-summing the so-called “daisy” diagrams, where one would expect to get a correction of order at least $\mathcal{O}(\lambda^{2})$. However, this is not the case in finite temperature calculations, since such diagrams induce correction $\mathcal{O}(\lambda^{3/2})$, which may be important especially for larger values of the self-interaction coupling. We have explicitly checked that for $\lambda\lesssim 1$ this re-summation leads to at most a $20\%$ variation in the thermal mass. We will thus use the approximate result eqs. (2.3),  throughout this paper.

2.2 Freeze-in and mediator decay

We are interested in estimating the final relic density of a DM particle $\chi$ interacting extremely feebly with the Standard Model particles. The key assumption is that $\chi$ 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

[TABLE]

with a small decay rate $\Gamma_{\chi}$ (i.e., such as to unable one to overproduce or thermalise the $\chi$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., Belanger:2018mqt for a complete recent treatment) given by

[TABLE]

where $\mathcal{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 factors $d\Pi_{S}\equiv\displaystyle\frac{d^{3}p_{S}}{(2\pi)^{3}2E_{S}}$, and similarly for $d\Pi_{\chi}$ and $d\Pi_{\bar{\chi}}$. Without loss of generality regarding the operator generating the decay $S\rightarrow\bar{\chi}\chi$, we can rewrite the squared amplitude from the decay rate $\Gamma_{\chi}$ as

[TABLE]

where $\lambda$ 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., Belanger:2018mqt ) 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_{\rm DM}{}=\frac{n_{\chi}+n_{\bar{\chi}}}{s}$) is given by

[TABLE]

with

[TABLE]

where $h$ ($g$) are the relativistic degrees of freedom associated with the entropy (energy) density,222To obtain our numerical results we use the ones provided in Drees:2015exa . and $K_{1}(x)$ the modified Bessel function of the first kind. Defining $x\equiv\dfrac{m_{S}}{T}$ and focusing for simplicity on $J_{S}=0$, the evolution of the yield becomes

[TABLE]

An important comment at this point is that, while in the standard freeze-in case $\Gamma_{\chi}$ 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_{\chi}$ and slowly varying relativistic degrees of freedom (which is the case for $T\gtrsim 1~{}{\rm GeV}$), we can calculate the yield today ($Y_{\rm DM,0}{}$), by integrating from the reheating temperature ($T_{\rm R}\gg m_{S}{}$, $x\to 0$) down until today ($T_{0}\ll m_{S}{}$, $x\to\infty$).333If DM is produced mainly at temperatures at which the assumptions are violated, $Y_{\rm DM,0}{}$ can obtained numerically from eq. (2.11). We then obtain the relic abundance in the form

[TABLE]

where we evaluate $g$ and $h$ at the “mean” value of $x$ during the DM production.444This value is defined as

$\langle x\rangle\equiv\dfrac{\int_{0}^{\infty}dx\,x^{3}K_{1}(x)\times x}{\int_{0}^{\infty}dx\,x^{3}K_{1}(x)}\approx 3.4\,.$

In Fig. 2(a) we show $Y_{\rm DM}{}/Y_{\rm 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_{\rm FI}\sim\frac{m_{S}}{7}$, as can be seen from the figure. That is, since typically $S$ decouples at temperature $T_{\rm FO}\approx 20\,m_{S}$ (i.e., freeze-out), the calculation holds. However, if $S$ decouples earlier than expected, the relic abundance of $\chi$ 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 $\chi$ 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 Hall:2009bx ; Elahi:2014fsa . As an example, consider DM production via a $2\to 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 $\sqrt{\hat{s}}$ as

[TABLE]

with $n=d-4$. The corresponding Boltzmann equation is

[TABLE]

which (assuming constant $g$ and $h$) can be integrated from $x_{\rm R}=m_{S}/T_{\rm R}$ to today ($x_{0}$). The result is

[TABLE]

where it is apparent that the high-temperature contributions dominate for $n>0$ (i.e., $d>4$). In the case of $d\leq 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.

2.3 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 form555In the case where $S$ gets its thermal mass due to the self interactions (2.1), $\alpha^{2}=\dfrac{\lambda_{S}}{24}$ (for $\lambda_{S}<1$).

[TABLE]

An important consequence of eqs. (2.14),  is that the decay $S\rightarrow\bar{\chi}\chi$ can become kinematically allowed at large temperatures even if $m_{S}{}<2m_{\chi}$. 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\leftrightarrow 2$ processes.666This is the case, e.g., for gravitino or axino dark matter at high $T_{R}$; for a review see, e.g. Baer:2014eja and references therein. In such cases the thermal effects typically introduce only a correction, the significance of which is very model-dependent.777Indeed, the dominant $2\leftrightarrow 2$ process should involve the coupling between the mediator and the thermal bath, since it is assumed that it is in equilibrium contrary to the dark matter. A proper estimation of this effect can thus be done only on a model-dependent basis, as studied later in Sec. 3 for the Higgs portal case. We will explicitly address the role of $2\leftrightarrow 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}{}=\alpha T$. Assuming that the temperature is large enough so that at some early time $m_{S,T}{}>2m_{\chi}$ is satisfied, our aim is to solve in this case the Boltzmann equation (2.11). Defining $z$ as

[TABLE]

we obtain

[TABLE]

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 eqs. (2.16),  increases with temperature and is therefore dominant at high temperatures close to the reheating temperature $T_{\rm 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 Roszkowski:2017nbc . On the other hand, when $d\leq 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 $\alpha T=2m_{\chi}$. The production is thus dominated by temperatures close to $m_{\chi}$ (or higher for small $\alpha$).

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

[TABLE]

where again $n=d-4$, and $\gamma_{S\chi}$ 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

[TABLE]

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 allowed for $z\geq 1$) and $z=z_{\rm R}\equiv\dfrac{2m_{\chi}}{\alpha T_{\rm R}}$ to obtain

[TABLE]

It is clear that, for $d>4$ the dominant contribution comes from the regime of high temperatures ($z_{\rm R}\to 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 Elahi:2014fsa ). That is, even if the decays $S\to\bar{\chi}\chi$ 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 eqs. (2.13),  with eqs. (2.19), , we can see that since $\alpha<1$, the later tends to be generally less efficient.888 Usually $2\to 2$ processes involve higher powers of couplings, and they are often subdominant to decays. For higher dimensional operators, however, both $1\to 2$ and $2\to 2$ processes can involve similar powers of the couplings. This is the case where this argument is applicable. 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 eqs. (2.19),  where the contribution from $z=z_{\rm R}$ drops out (unless $\alpha$ 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\to\bar{\chi}\chi$ stops. Thus, we expect the production to be dominated at time scale corresponding to the temperature at which $m_{S,T}{}\sim 2m_{\chi}$. This actually implies that up to an order one function, the decay rate satisfies $\Gamma_{\chi}\propto m_{\chi}$. While it is therefore not possible to fully simplify the decay rate without specifying the details of the interaction, we can straightforwardly observe from eqs. (2.18),  that the abundance will be proportional to $1/m_{\chi}$, 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 $\chi$ is a Dirac fermion and the Lagrangian contains an Yukawa interaction between $S$ and $\chi$,

[TABLE]

The bath particle $S$ decay width to dark matter is then given by

[TABLE]

The evolution of the yield is then given by

[TABLE]

In Fig. 2(b) we show the evolution $Y_{\rm DM}{}/Y_{\rm 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_{\chi}$.

Assuming that the relativistic degrees of freedom do not vary rapidly during the production of the $\chi$s, we can integrate eqs. (2.22),  to obtain

[TABLE]

where again $g$ and $h$ are evaluated at $\langle z\rangle$.999In this case, $\langle z\rangle$ is defined as

$\langle z\rangle\equiv\dfrac{\int_{0}^{1}dz\;\left(1-z^{2}\right)^{3/2}\times z}{\int_{0}^{1}dz\,\left(1-z^{2}\right)^{3/2}}\approx 0.34\,.$

As we pointed out earlier, the relic abundance of $\chi$ becomes (mostly) independent of its mass, with any $m_{\chi}$ dependence coming from $\langle z\rangle$. This, and the suppression due to the $\alpha^{4}$, will result in relaxed constraints for the Yukawa coupling, with respect to the standard freeze-in, where $\Omega h^{2}$ scales predominantly linearly with the DM mass. Notice furthermore that in the case where the temperature correction never dominates (i.e.  $\alpha\,T<m_{S}{}$), the relic abundance is given by eqs. (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_{\rm DM,0}{}$ by including both mass terms. That is, the evolution of $Y_{\rm DM}{}$ as in eqs. (2.7),  needs to be solved, with $m_{S,T}{}$ given by eqs. (2.14), , numerically.

An example of typical dependence of $\Omega h^{2}$ on $m_{\chi}$ for the production of DM due to the decay of $S$, is shown in Fig. 3(a). The two extreme cases of $\alpha=0$ (standard freeze-in) 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_{\chi}\approx m_{S}{}/2$ which is where the blue line terminates since $S\to\bar{\chi}\chi$ becomes forbidden in the vacuum.

In Fig. 3(b) we present the Yukawa coupling $y_{\chi}$ as a function of $\alpha$ that give the observed $\Omega h^{2}$ for the scanned range of masses $10\;{\rm MeV}\leq m_{S}{},m_{\chi}\leq 1\;{\rm 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_{\chi}$ is marked in blue, while the forbidden freeze-in region of $m_{S}{}<2m_{\chi}$ is marked in orange. The shape of the forbidden freeze-in band in Fig. 3(b) is a simple consequence of the $\alpha^{2}y_{\chi}$ dependence of $Y_{\rm DM,0}{}$ in eqs. (2.23), . As already noted in the $d>4$ case, in the forbidden freeze-in regime one requires either larger self-interaction $\alpha$ 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. 3(a)). An important comment is that the transition between the two regimes, which happens for $m_{S}{}\sim 2m_{\chi}$ occurs typically in a mass range of order $(2m_{\chi}-m_{S}{})\sim\alpha m_{S}{}$, which become very narrow for small $\alpha$.101010This 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\sim m_{S}{}$. In particular, in Figure 3(b) the forbidden region shown in orange do not probe this tuned transition regime in details for small $\alpha$.

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.

3.1 The model

We introduce a real scalar “dark Higgs” boson field $S$, which is not protected by a $\mathcal{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_{\chi}$. The corresponding part of the Lagrangian thus reads

[TABLE]

with the dark Higgs boson potential term defined as 111111Notice 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 $\chi$-number global symmetry to fix the latter to zero.

[TABLE]

where $H$ denotes the Standard Model Higgs boson doublet. The total scalar potential is $V=V_{HS}-\mu\,H^{\dagger}H+\dfrac{\lambda_{H}}{2}\,\left(H^{\dagger}H\right)^{2}$.

At low temperatures ($T\lesssim 160\;{\rm GeV}$), both the Higgs and dark Higgs fields develop a non-zero vacuum expectation value (VEV), so that $H=\dfrac{1}{\sqrt{2}}\left(\begin{matrix}0\\ h+v\end{matrix}\right)$ and $S\to v_{S}+S$. 121212Since $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. In the limit where $A\ll v$ the calculation simplifies significantly and the minimization conditions for the scalar potential in term of $\lambda_{H}$ and $v_{S}$ can be easily obtained as

[TABLE]

Furthermore, we can rotate the scalars to their eigenvalue basis, i.e.$\left(\begin{matrix}h,S\end{matrix}\right)\to R\left(\begin{matrix}h,S\end{matrix}\right)$, where $R$ is a rotation matrix parametrised by the small angle $\theta$ given by

[TABLE]

where we have used the masses of $h$ and $S$ (at $T=0$) defined by

[TABLE]

The branching ratio of the Higgs decay to invisible particles is constrained to be Belanger:2013xza smaller than $0.19$, which translates to $\lambda_{HS}\lesssim 10^{-2}$. Furthermore, note that while we have supposed that the trilinear term $\lambda_{3}S^{3}$ was negligible in our original Lagrangian,131313For example, this term can shift the thermal mass of $S$ by a factor of $\mathcal{O}\left[\lambda_{S}\left(\dfrac{v_{S}}{\mu_{S}}\right)^{2}\right]$. the shift by $v_{S}$ re-introduces such a term as $\dfrac{\lambda_{S}}{3!}v_{S}S^{3}$. For consistency, we will therefore further require that this contribution is negligible with respect to $\mu_{S}$, leading to the condition

[TABLE]

Notice that this also automatically ensures that the shift in the SM Higgs boson quartic coupling $\lambda_{H}$ is negligible in eqs. (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 $\ell$ pair are kinematically allowed, and assuming $\mu_{S}\sim m_{S}{}$, we obtain

[TABLE]

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.141414Note, though, that strictly speaking one could still satisfy the above bounds while keeping the $\tau_{S}\sim 0.1$ s, for very low values of $\lambda_{S}$. The parameter space is however extremely restricted experimentally, as we will see in Sec. 3.3.

The relevant processes determining the evolution of number densities of $S$ and $\chi$ in this model are: i) the direct mediator decay $S\rightarrow\bar{\chi}\chi$, 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.151515Note that for heavy mediators a decay/annihilation channels to $h$ could be open resulting in additional two processes $S\rightarrow hh$ and $SS\rightarrow hh$, governed by $A$ and $\lambda_{SH}$, respectively. These are not relevant in our analysis, which is focused on the regime $m_{S}{}<m_{h}$ for the temperatures around the freeze-out temperature of $S$, and therefore are not included. The direct $S\to\bar{\chi}\chi$ decay width is given by eqs. (2.21),  and is suppressed by the very small Yukawa coupling $y_{\chi}$. The decay of $S$ to SM particles is given by $\Gamma(S\rightarrow{\rm SM})=\theta^{2}\;\Gamma_{h\to{\rm{SM}}}(m_{S,T}{})$, where the $\Gamma_{h\to\rm{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 Dittmaier:2011ti ; Fradette:2017sdd ; Winkler:2018qyg and a direct evaluation for leptonic decay at low masses.

The $S$ annihilation cross section as a function of the Mandelstam variable $s$ reads

[TABLE]

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 $\lambda_{HS}$ or the mixing angle $\theta$ 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\leftrightarrow 2$ processes can in principle play a role in the production of $\chi$s and/or their early-time thermalization with the SM plasma. These are: $SS\leftrightarrow\bar{\chi}\chi$, $hh\leftrightarrow\bar{\chi}\chi$ and the co-annihilation process $Sh\leftrightarrow\bar{\chi}\chi$. The first one has $s-$, $t-$ and $u-$channel contributions which are proportional to $\theta^{2}\lambda_{HS}^{2}y_{\chi}^{2}$ and $y_{\chi}^{4}$, respectively. The second and third have only $s-$channel diagrams proportional to $A^{2}y_{\chi}^{2}$ and $\lambda_{HS}^{2}y_{\chi}^{2}$, respectively. It is clear that all of theses $2\leftrightarrow 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 $\lambda_{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 $\chi$. 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 $\chi$ number densities.

3.2 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 $\chi$ is performed under the assumptions that: i) $\chi$ 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 $\chi$s.161616Kinetic 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 Binder:2017rgn . 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 eqs. (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. 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\equiv m_{S}{}/T$. For $x<0.1$ it is assumed that $S$ traces its equilibrium value while the evolution of $\chi$ is given by eqs. (2.11), , starting from the reheating temperature $T_{\rm R}$ assumed to be given by $x_{\rm R}=10^{-9}$. We checked explicitly that assuming different $T_{\rm R}$ does not change the result. For $x>0.1$ the coupled system of the Boltzmann equations for the number densities of $S$ and $\chi$ is numerically solved, including all the relevant processes discussed above.171717This is done to ensure that the $\chi$ 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. Belanger:2018mqt .

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 $\chi$ for different regimes and then present and discuss the results of our scan of the parameter space of the model.

3.2.1 Evolution of number densities

In Figures 5–7 we present the yields of $S$ and $\chi$ for some characteristic cases. In all following figures the green dashed lines correspond to $Y_{S}$ while the solid lines to $Y_{\rm 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_{\chi}$ and $Y_{\rm DM}$ gradually grows, with most of the production happening around $T\sim m_{S}{}$. This is shown in Fig. 5(a). In this case the final relic abundance of $\chi$ is insensitive to any variations in the self-coupling $\lambda_{S}$ due to the fact that the thermal effects are important only for $T\gg 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. 5(b) 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 $\lambda_{S}$. In this forbidden regime the final DM abundance is therefore very sensitive not only to value of $y_{\chi}$ but also the self-coupling of the mediator. Another point worth stressing is that one does not need large values of $\lambda_{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 6(a) shows a case of a transition between the standard and the forbidden regimes. For fixed $m_{S}{}=100$ GeV we vary $m_{\chi}$ 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_{\chi}$ further leads to only very mild change in the relic abundance, i.e., the yield $Y_{\rm DM}$ is inversely proportional to $m_{\chi}$, in agreement with eqs. (2.23), . This approximate DM mass independence of the relic density is an distinct feature of the forbidden freeze-in scenario.

In Fig. 6(b) a slightly different mechanism is shown. It occurs when nominally this would be a standard freeze-in case with $m_{S}{}>2m_{\chi}$ 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 $\lambda_{HS}$), there appears a temporary regime where $S\to\bar{\chi}\chi$ is not allowed and the $\chi$ production is blocked for a while. However, if the self-coupling $\lambda_{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_{\chi}$. A scenario like this is close to what was studied, in a more general context, in ref. Baker:2017zwx .

Finally, in Fig. 7 we show for completeness examples of cases where the $\chi$ 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 $\chi$ is crucial. In Fig. 6(a) the final DM abundance is determined by the branching fraction of the $S$ decays to $\chi$ 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 7(b) shows a situation where the details of the freeze-out of $S$ strongly affect its abundance that is then transferred to the $\chi$s via (rare) decays. This also shows the potential impact that the choice of $\lambda_{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 $\lambda_{HS}$.

3.2.2 Scan setup and results

A numerical scan of the model parameter space has been conducted using MultiNest Feroz:2008xx to direct the scan towards values of the relic density within $2\sigma$ of the standard result from the Planck Collaboration Aghanim:2018eyx $\Omega h^{2}=0.1198\pm 0.0012$ that we set as an allowed range.181818We used an additional $10\%$ theoretical uncertainty on our numerical results. The private code BayesFITS, automatically created using routines from SARAH staub_sarah_2008 ; Staub:2012pb ; Staub:2013tta 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 $\lambda_{S}$. It also requires very low values of the Yukawa coupling; otherwise DM is overproduced. In contrast, the forbidden regime is highly sensitive to $\lambda_{S}$, as expected. Indeed, the smaller the self-coupling, and therefore the thermal mass, the earlier the production stops and therefore the larger $y_{\chi}$ 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 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 $\lambda_{SH}$. This is a direct consequence of eqs. (3.5), .

3.3 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 $\theta$, hence mediating the decay of $S$ to SM particles (if kinematically allowed). Since $\theta$ is suppressed by powers of $v_{S}/v$, the dark Higgs boson $S$ is typically long-lived, as shown in eqs. (3.7),  – particularly for low masses. In this case bounds from both colliders and fixed target experiments Beacham:2019nyx , and for longer life-time, from astrophysics Fradette:2017sdd 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}\sim 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 Krnjaic:2015mbs are derived from the requirement that the cooling rate from dark sector particles do not exceed the neutrinos one Burrows:1986me ; Burrows:1987zz ; Raffelt:1996wa .191919While we use here the results from Krnjaic:2015mbs , this should be considered only an order of magnitude calculation. Note, however, that for the parameter space presented in Sec. 3.1, this bound is not directly relevant as can be seen in Figure 9.

• •

Bounds from enforcing a successful big bang nucleosynthesis. We use the recent bounds from Fradette:2017sdd which are derived from the same Lagrangian as in Sec. 3.1. In the lower mass range (below the $\pi$-meson mass threshold) the dominant bounds are derived by constraining the entropy injections from the $e^{+}e^{-}$/ $\mu^{+}\mu^{-}$ 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\sim 1/7$ due to the $p\leftrightarrow 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 eqs. (3.6), , dark Higgs bosons abundance typically freezes-out earlier than in Fradette:2017sdd which implies that the relativistic abundance is maintained for larger masses. Altogether, modifying $\lambda_{HS}$ only changes the limits by an $\mathcal{O}(1)$ factor, as can be seen in Fradette:2017sdd . 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}/\tau_{S})$ where $t_{p/n}\sim 2.6s$ is the freeze-out time of the proton/neutron ratio.202020 This behaviour is clearly illustrated in Figure 2, from ref. Fradette:2017sdd

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 Bergsma:1985qz , 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 Winkler:2018qyg . Note that we have included the projected limits from SHiP at $2\times 10^{20}$ proton-on-target Alekhin:2015byh 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 Feng:2017vli . 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^{+}\rightarrow\pi^{+}\nu\nu$ studied in the E949 experiment Artamonov:2008qb . 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 Aaij:2016qsm . In both cases, we use the recasted bound from Winkler:2018qyg .

Note that in the long term, several planned experiment have the potential to greatly improve the limits in this mass range Beacham:2019nyx . LHC-based experiments, such as FASER, MATHUSLA Curtin:2018mvb or CODEX-b Gligorov:2017nwh are particularly interesting in that the decay of Higgs boson mediated through the quartic mixing $\lambda_{HS}$ can significantly enhanced the detection prospects as they are not tied to the mixing angle per-se but only to $\lambda_{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 Beacham:2019nyx .

4 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\leq 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.

Acknowledgments

LD, DK and LR are supported in part by the National Science Centre, Poland, research grant No. 2015/18/A/ST2/00748. AH is supported in part by the National Science Centre, Poland, research grant No. 2018/31/D/ST2/00813. LR is also supported by the project “AstroCeNT: Particle Astrophysics Science and Technology Centre” carried out within the International Research Agendas programme of the Foundation for Polish Science financed by the European Union under the European Regional Development Fund.