Freeze-in through portals

The popular freeze-out paradigm for Dark Matter (DM) production, relies on DM-baryon couplings of the order of the weak interactions. However, different search strategies for DM have failed to provide a conclusive evidence of such (non-gravitational) interactions, while greatly reducing the parameter space of many representative models. This motivates the study of alternative mechanisms for DM genesis. In the freeze-in framework, the DM is slowly populated from the thermal bath while never reaching equilibrium. In this work, we analyse in detail the possibility of producing a frozen-in DM via a mediator particle which acts as a portal. We give analytical estimates of different freeze-in regimes and support them with full numerical analyses, taking into account the proper distribution functions of bath particles. Finally, we constrain the parameter space of generic models by requiring agreement with DM relic abundance observations.


I. INTRODUCTION
By now there is solid evidence for the existence of Dark Matter (DM) from a plethora of observations such as galaxy rotation curves, structure formation, the Cosmic Microwave Background spectrum or gravitational lensing [1]. These observations comprise one of our few precious evidences of physics beyond the Standard Model. Unfortunately, all the present proofs for DM stems from its gravitational effects and thus, we remain ignorant of the particle nature of DM, i.e., its mass and interactions with the rest of the known particles, crucial ingredients so as to be able to embed it in a complete theory. This has led to a proliferation of many different DM candidates with radically distinct phenomenologies and genesis mechanisms, with the most popular one being a Weakly Interacting Massive Particle (WIMP) [2]. The source of the WIMP popularity can be attributed to a combination of the Higgs hierarchy problem, whose different explored solutions typically require new weakly-interacting particles not much above the electroweak scale, and the so-called "WIMP miracle", referring to the fact that the correct DM relic thermal abundance can be obtained through these candidates.
However, other mechanisms are equally viable and should also be explored in case DM turns out not to be a WIMP, lest we miss its signals by concentrating exclusively on the WIMP paradigm. An interesting alternative, motivated partially by the failed DM searches, is the case in which the coupling of DM to the visible sector is very suppressed. In this scenario, DM does not thermalize with the visible sector and so it tends to approach its final density from below, increasing it with increasing cross section (in contrast to the situation in the WIMP framework). This scenario has been recently referred to as freeze-in [3].
Several mechanisms have been proposed in the past describing an out-of-equilibrium production of DM. In [4] a model with sterile neutrinos as DM candidates was analysed, where these particles are populated from the thermal bath through oscillations with the SM neutrinos suppressed by small mixings. In [5], a scenario with a gauge scalar singlet DM candidate is studied, where the out-of-equilibrium DM genesis is produced from the decays of the Higgs bosons present in the thermal bath. Similar alternatives have been analysed in [6] and [7]. Models where the DM candidate is produced from processes like (b → bχ), wherẽ b and b are thermal bath particles and χ the DM, have also been studied in the literature.
An example of this has been analysed in [8] in the context of Supersymmetry, and afterwards in [3] several combinations of the masses mb and m χ are analysed in detail. In [9] a case of a gravitino DM production, which is actually dominated by high temperatures has been considered [10][11][12][13][14]. An interesting model which is also sensitive to higher temperatures is described in [15], where a DM candidate whose mass is larger than the reheating temperature is studied. A scenario where the portal is massless has been analysed in [16], while an alternative where the portal is heavier than the reheating temperature has been recently proposed in [17].
In this work we will extend a possibility that was only briefly discussed in [3], namely that the weak interaction between DM and the visible sector is mediated by a portal. This is complementary to the case where the mediator is part of the thermal bath, which, as we have commented above, has been extensively analysed in the literature. The present study contains and generalises specific examples already present in the literature, such as the works presented in [16][17][18][19]. We identify the regions of parameter space where the portal (say, particle P) is not in thermal equilibrium with the thermal bath, thus the scattering process (bb → χχ) can dominate over the decay process (P → χχ) when populating the DM. We obtain analytical estimates of the predicted relic abundance coming from typical models, which we classify according to the mass of the portal.
In section II we present the formalism used to study the evolution of the DM number density in the freeze-in regime. Section III is devoted to describe in detail the approximate analytical solutions of the Boltzmann equation, where the classification of different freeze-in regimes is introduced. We cross-check and fit our analytical estimates with full numerical results in section IV, and comment on the possible phenomenology of some of these models in section V, before concluding in section VI.

II. DM YIELD FROM FREEZE-IN
The evolution of the DM (χ-particle) number density, in the case where finite temperature effects are neglected 1 , is described by the usual Boltzmann equation which can be expressed as: when considering a generic 2-to-2 process bb ↔ χχ where the DM (with number density n χ ) is produced from (and annihilates to) bath particle pairs bb. The factor a in the LHS is the scale factor; the index i runs over the four particles, and f i are the thermal distribution functions, which for particles in thermal equilibrium will be given by Upper (lower) signs correspond to bosons (fermions). Furthermore, P in(out) are the incoming (outgoing) 4-momenta of the process and |M| 2 is the amplitude of the process, summed over all spins.
The freeze-in scenario assumes that the initial abundance of DM (at reheating epoch) is negligible. The thermal bath then starts populating the DM through interactions sufficiently suppressed for the DM not to thermalize with the bath. If this is the case, then the backreaction annihilation term in Eq. (1), proportional to f χ fχ(1 ± f b )(1 ± fb), can be safely neglected, simplifying the expression to where we have made the approximation of m b s and we have taken Maxwell-Boltzmann distribution functions (f i ≈ e −E i /T 1). In section IV we evaluate how good this approximation is when we cross-check our analytical estimates with the full numerical results which take into account the appropriate distribution functions instead. Following [21], the integrals over initial 3-momenta p b and pb have been re-expressed to the variables E + ≡ E b + Eb, where σ ≡ σ bb→χχ is the unpolarised cross-section and g b are the degrees of freedom of the b-particles. K 1 is the order-1 modified Bessel function of the second kind.
Taking into account that a −3 d(n χ a 3 )/dt = −sHT dY χ /dT , where Y χ is the comoving number density, or yield, (Y ≡ n/s), s the entropy density, and H the Hubble parameter, we can finally express the DM relic density as: or where the 0-subindex refers to the values today, and we have considered a symmetric scenario where nχ = n χ . Here T R is the reheating temperature, which acts as the initial condition in scenarios where the reheating epoch is assumed to be instantaneous. M P l is the Planck mass and g * (g s * ) are the energy (entropy) density effective degrees of freedom. The dependence of σ(s) on s is the essential ingredient to know the behaviour of the DM yield Y χ . Here we concentrate on theories for which the DM is generated through "portal" interactions, where these portals are particles directly interacting with both the bath and the DM, whereas the DM only interacts directly with the portal. The dominant processes populating the DM sector are shown in Fig. 1. In what follows, we will refer to the couplings λ BB for bath-to-bath, λ χP for DM-to-portal interactions, and λ BP for the the bath-to-portal interactions, where either λ BP or λ χP should be small for the DM to be out of equilibrium. The decay of the mediator depicted in case b) corresponds to the case that was discussed in more detail in Ref. [3]. This process will tend to dominate when the portal thermalizes with the thermal bath and its mass M is such that 2m χ < M < T R . However, as we will show, there are regions of the parameter space where the portal particle does not thermalize either and production processes via freeze-in of the type depicted in a) can dominate and lead to the correct relic abundance for DM. Thus, as a complementary view to Ref. [3], in this work we will concentrate on freeze-in via a portal particle, of the type depicted in a), discussing the allowed parameter space to obtain the correct relic abundance as a function of the mediator mass, which will characterize different regimes with distinct dependence on the parameters. For each regime, we will derive approximate analytical expressions and compare our estimates with full numerical evaluation of the integral in Eq. (1), taking into account the appropriate distribution functions for the corresponding particles.

III. ANALYTICAL RESULTS
For the analytical estimates of Eq. (4) it is useful to note the different limits of the Bessel In Fig. 2, we depict the region of integration where the integrand is not exponentially suppressed. Essentially for √ s T the DM production is negligible due to the huge Boltzmann suppression, which means that the integral over s can be estimated introducing a cut-off close to s T 2 . A naive estimate for the cut-off is s max 9T 2 , since 100. Thus, beyond √ s/T = 3 the contribution to the integral is expected to be negligible.
On the other hand, the suppression from K 1 may be balanced by an enhancement from the cross section σ (see Eq. (4)). However, we will see below that this rough estimation is actually rather accurate. But we will keep the cut-off parameter B such that s max = (BT ) 2 free in order to compare with the exact numerical results and choose the value of B that best reproduces them, so that the analytical approximation can become an accurate proxy of the full numerical simulation.
In this section we will adopt a generic cross section for an s-channel process given by: in order to discuss the different regimes for the mediator mass M .

B
A C Before reheating

FIG. 2: Relevant parameter space of DM production by freeze-in.
A. Heavy mediator: M > T R We will refer to a heavy mediator of mass M when M > T R , where T R is the reheating temperature, which acts as a cut-off for the integral over the temperature (see Eq. 4). In this case, the cross-section has the following dependence with s: where we have assumed all the particles other than the portal to have masses m 2 i s.
As before, λ BP is the coupling of the visible (SM) sector to the portal, whereas λ χP is the coupling between the portal and the DM. From inspection of Eq. (8) it can be noted that in this case the DM production is dominated by the largest temperatures, given that the cross-section grows with s. We thus expect a direct dependence of the relic abundance on the reheating temperature. Indeed, for the relic abundance we obtain: This actually constitutes a special case of a more general result, where the DM production happens through an effective, non-renormalizable operator of dimension N and thus suppressed by a scale Λ 4−N . By dimensional analysis, the relic abundance must behave as: Recently, a model with this characteristic has been presented in [17], which has been dubbed NETDM. The idea is that a GUT framework, e.g., SO (10), can naturally provide a very heavy portal through which the Standard Model populates the DM. For several cases of SO(10) breaking patterns, mediator masses can be larger than 10 10 GeV, which require, according to Eq. (9), large reheating temperatures in order to obtain the correct relic abundance.
B. Light mediator: M < 2m χ We define as the light regime a portal whose mass M is such that M < 2m χ . In this case the cross-section can be approximated by the following expression: Contrary to the heavy-mediator case, now the DM production is dominated by the lowest temperatures, given the energy dependence shown in Eq. (11). Thus, since the DM production always stops at the "freeze" time (T 2m χ ), a priori one expects that the final yield depends directly on m χ . Indeed, the relic abundance is for this case: Note that the m χ -dependence of the yield (proportional to the term in squared brackets above) cancels when computing the relic abundance Ωh 2 . It implies that for the light mediator the relic abundance is much less sensitive to the DM mass (as compared to the heavy and intermediate mediator cases). Indeed, the dependence on m χ in this regime only stems from the effective degrees of freedom g * (m χ ).
A particular model with these characteristics has been analysed in [16], where a "dark photon" which acts as a massless mediator of a feeble interaction between the SM and DM is considered.
There is a range of mediator masses in between the regimes described above, namely 2m χ < M < T R . In this region, the pole of the mediator propagator will be accessible during the DM production and thus we can split the integration region in three zones with qualitatively different behaviours as depicted in Fig. 2. Within region A, the mediator mass can be considered heavy and the integrand is given by the same expression as in the heavy mediator case, while in region B the mediator can be considered light and the integrand is given by the same expression as for the light regime. Thus, the integrals over A, and B can be approximated much in the same way as seen above. These results are typically smaller than the contribution from the peak (region C) around √ s M assuming that Γ M , the Breit-Wigner peak can be approximated using the relation: leading to The integral over the small width of the region C is then simply taken care of by the delta function and the remaining temperature integral for the yield is Since in this regime the mediator is heavier than the DM particles, the mediator can decay directly into DM. If the mediator thermalizes with the bath, its decay would instead dominate the DM production through freeze-out of the mediator and its subsequent decay to DM as described in Ref. [3] 2 . However, we will show that for large regions of the parameter space in which the intermediate mediator scenario provides the correct relic density, the mediator does not thermalize and its decay is not such an effective way of increasing the DM population.
As an example we will consider the relevant case of a vectorial mediator (e.g. a massive dark photon) coupling to the SM bath via kinetic mixing with the photon. In this scenario, the mediator will couple to SM electrons 3 with λ BP = λ 0 e, where e is the electron charge and λ 0 the mixing between the mediator and the photon. Three processes could in principle lead to the thermalization of the mediator: its direct production through coalescence in a collision of electron and positron (notice that the cross section of this process is proportional to δ( √ s − M )); its production in e − e + annihilation in association with a photon; or via inverse Compton scattering with a dark photon instead of a photon in the final state. The rates for these processes can be found in Ref. [22].
In Fig. 3 we compare these three production rates with the Hubble rate. For λ BB = e, λ BP = 10 −11 e and M = 10 TeV. We can see that all production processes are at least 10 orders of magnitude smaller than the Hubble rate for any temperature, thus preventing thermalization. As we will show with our numerical results in Sect. IV, these choices of the parameters can lead to the correct DM abundance. For M = 10 TeV the coupling would need to be increased by about 6 orders of magnitude in order for the production rates to increase above the Hubble rate and reach thermalization. At temperatures somewhat higher than the mediator mass a spike in the production rates of the mediator appears. This spike corresponds to the temperature at which the mediator mass is equal to the thermal mass of the photon, leading to a resonantly enhanced mixing between the two [22].

IV. NUMERICAL RESULTS
In this section we cross-check the analytical estimates made above with a full numerical analysis for the three regimes defined. For the sake of illustration, we take two type of amplitudes of DM production from the thermal bath (assumed for concreteness here to be the SM fermions), corresponding to a vector interaction: and a scalar interaction: where Γ is the total decay width of the mediator of mass M , taken here to be λ 2 χP M/4π and λ 2 χP M/8π for vector and scalar interaction, respectively. We have also assumed for simplicity only an interaction with electrons (i.e., m b = m e ), but the conclusions would be similar for more complex models.
After solving (5)  We compute the resulting DM abundance as a function of the DM mass in three different ways: analytically, numerically using the MB approximation, and numerically using the correct distribution functions for the bath particles. This is done for the three different regimes: heavy, light and intermediate mediator. The analytical estimation has been performed through Eqs. (9), (12) and (15), corrected for with appropriate factor to take into account the assumed vector (scalar) type of interaction: 1/3π (1/4π), 3/8π (3/16π) and 1/3π (1/4π), respectively. We show the results in Fig. 5. For the heavy regime (upper left panel) we used M = 10 14 GeV for the mediator mass, λ BP = 1 for the bath-to-portal coupling, λ χP = 1 for the portal-to-DM coupling, and T R = 10 10 GeV for the reheating temperature. The first thing to note is a very reasonable agreement between the numerical results using both the appropriate distribution functions (solid lines) or the MB approximations (dashed lines), for the two models considered in Eqs. (16) and (17). The discrepancy turns out to be around 15 %.
(1) should not appear, since here we work under the assumption that DM particles are way out-of-equilibrium, which translates to f χ 0. Note also that in the more complete way of solving Eq. (1) would be to solve the complete integro-differential equation. However given the order of effective bath-to-DM couplings we need to obtain good relic abundances, the assumption of neglecting f χ works extremely well.  Regarding the analytical estimates (dotted lines), we provide in Table I   We see that the fitting values are close to the naive estimation B 3 made above, and they are useful to provide a very accurate analytical approximation to the exact numerical results and reproducing the correct parametric behaviour, as can be seen in Fig. 5. others. Notice that, since in the light regime Ω DM h 2 is independent of M and very weakly dependent on m χ , a coupling λ BP ∼ 10 −11 is required to obtain the correct relic abundance.

V. PHENOMENOLOGY
In the scenario explored here the bath-portal coupling is necessarily very small, which makes DM probes through direct, indirect and colliders searches rather challenging. However, the DM-portal coupling could be sizable. Indeed, as we have shown, the product of the two couplings should not be too small in order to obtain the correct relic abundance. A sizable DM-portal coupling would then change the paradigm of collisionless DM for structure formation, leading to DM self-interactions which could even have relatively long distance forces depending on the mass of the portal. The structure formation phenomenology of these models is thus altered, particularly at small scales, and can be probed through observations.
On the one hand, DM self interactions can be very directly bounded by the X-ray and lensing observations of the Bullet cluster to σ/m χ < 1 cm 2 /g [23], which clearly shows a separation of the luminous and dark matter components through a weaker scattering of the latter. DM self interactions can also affect the ellipticity of clusters. Indeed, self interactions would tend to thermalize the DM velocity spectrum and lead to more spherical shapes. These observations actually lead to the strongest constraints on DM self interactions σ/m χ < 0.02 cm 2 /g [24]. However, these bounds have been relaxed more recently through more detailed numerical simulations which show that cross sections as large as σ/m χ = 0.1 cm 2 /g [25] are agreement with all observations.
On the other hand, DM self interactions could even solve some of the experimental ogy is complementary to that of the scenario with a more strongly coupled portal in the thermal bath and very feebly interacting dark matter, which would not lead to these modifications of structure formation but would typically present more prominent phenomenology at colliders [3].
Nevertheless, the light regime can also be probed in direct detection experiments; specifically for mediators lighter than the recoil energy E r , since the scattering cross-section has an infrared divergence as E −2 r , even if the coupling is very tiny. An interesting alternative analysing electron recoils has been presented in [28], but more common experimental studies based on nuclear recoils have promising prospects, being able to test sufficiently feeble couplings in a few years from now (see e.g. [16]).
Finally, for the heavy regime the phenomenology is much more challenging given that the main phenomenology of these models takes place at the reheating period, for which there are at present no direct probes. Apart from the allowed range 1MeV T R 10 14 GeV, the lower bound coming from BBN and the upper bound from a typical prediction of chaotic inflationary models [29], there is no prospect for constraining T R better than this in the near future. Assuming that the observed DM content comes solely from a heavy-mediated candidate, the above range can be translated into

VI. CONCLUSIONS
In this work we have concentrated in the so-called freeze-in mechanism for dark matter (DM) production. In this framework, the genesis of DM happens out of thermal equilibrium, since its connection to the thermal bath is assumed to be very suppressed. Thus, if there is a portal mediating this interaction, the product of the couplings of the portal with the bath and DM must consequently be small. Here we have focused on scenarios in which the portal is also out of thermal equilibrium because the bath-to-portal coupling is suppressed, while the DM-to-portal coupling could be sizable. This scenario is complementary to the more discussed freeze-in case where DM genesis occurs from the out-of-equilibrium decays of a particle which is part of the thermal bath, characterized by a sizable bath-to-portal coupling, whereas DM is only feebly interacting.
We have performed analytical estimates of the DM relic abundance for the different regimes that can be identified according to the mass of the portal. These analytical results are based on the assumption that the distribution function of all particles follow a Maxwell-Boltzmann law, which a priori is not justified in processes for which the temperature T is much greater than the masses of all relevant particles during production. Furthermore, the resulting Bessel function has been approximated by a simplified expression in the region of interest. We have studied the size of the corrections driven by these simplifications by cross-checking our analytical estimates with a complete numerical analysis. We found that the Maxwell-Boltzmann approximation is reasonable with a discrepancy of around 15%. On the other hand, the analytical approximation, while maintaining the correct parametric behaviour, strongly depends on how the production region is approximated. By comparing with the exact numerical results, we have obtained the size of the integration region which allows to reproduce with very good accuracy these results with the simple analytical expressions derived. Thus the analytical expressions can be safely used instead of the exact numerical results with the corresponding correcting factor. Finally, we have used the exact numeric results to set constraints on the parameter space of generic models (i.e. masses and couplings of DM and the portal) so as to obtain a correct DM abundance as measured by WMAP9 [30] or more recently, Planck [31].