Neutrino Portal to FIMP Dark Matter with an Early Matter Era

We study the freeze-in production of Feebly Interacting Massive Particle (FIMP) dark matter candidates through a neutrino portal. We consider a hidden sector comprised of a fermion and a complex scalar, with the lightest one regarded as a FIMP candidate. We implement the Type-I Seesaw mechanism for generating the masses of the Standard Model (SM) neutrinos and consider three heavy neutrinos, responsible for mediating the interactions between the hidden and the SM sectors. We assume that an early matter-dominated era (EMDE) took place for some period between inflation and Big Bang Nucleosynthesis, making the Universe to expand faster than in the standard radiation-dominated era. In this case, the hidden and SM sectors are easily decoupled and larger couplings between FIMPs and SM particles are needed from the relic density constraints. In this context, we discuss the dynamics of dark matter throughout the modified cosmic history, evaluate the relevant constraints of the model and discuss the consequences of the duration of the EMDE for the dark matter production. Finally, we show that if the heavy neutrinos are not part of the thermal bath, this scenario becomes testable through indirect detection searches.


Introduction
The origin of dark matter (DM) is still one of the most important open problems in cosmology and particle physics. Dark matter is required to explain the galaxy rotation curves, the anisotropies of the Cosmic Microwave Background (CMB) and the observed structure of the Universe on large scales. However, despite the large number of viable DM candidates (long-lived, cold, and sufficiently abundant), the nature of DM remains unknown (see Ref. [1] for a review).
Weakly Interacting Massive Particles (WIMPs) are the most popular DM candidates. These particles attained thermal equilibrium with the cosmic plasma in the early Universe and when their interactions with the Standard Model (SM) particles could no longer keep up against the expansion of the Universe, they decoupled from the thermal bath, yielding a frozen-out abundance. However, the lack of evidence of such DM candidates in detection experiments, the strong constraints on this framework [2] and the absence of new particles at the LHC motivate exploring alternative scenarios.
An interesting alternative to the WIMP paradigm is provided by the freeze-in mechanism [3][4][5]. According to this scenario, the DM abundance results from decays and annihilations of particles of the SM thermal bath and, due to small couplings between the dark and the visible sectors, the DM never reached chemical equilibrium with the cosmic bath. Since the interactions between DM and the SM particles are so feeble, this kind of DM candidates is called Feeble Interacting Massive Particles (FIMPs). The small couplings needed for the freeze-in allow to easily evade the stringent observational constraints, which makes the model appealing, although more difficult to test.
Another way to evade experimental constraints on DM is to assume a non-standard cosmological history of the Universe. Although we know that the Universe was radiationdominated at the time of the Big Bang Nucleosynthesis (BBN), nothing prevents us to assume that another component dominated the Universe at early times, and it is interesting to study the consequences for the DM production. An early matter-dominated era (EMDE) would take place if some pressureless fluid -such as inflaton candidates [6], meta-stable particles [7][8][9] and moduli fields [10][11][12] -dominated the energy density of the universe prior BBN. In fact, when the EMDE is over, the dominant matter content might decay into SM degrees of freedom [13], reheating the visible sector and diluting the DM abundance. If DM was initially coupled to the thermal bath, the freeze-out needs to happen earlier than in the usual radiation-dominated era to overcome such dilution and agree with the relic density constraints. So, the interaction strengths of WIMPs to the thermal bath need to be weaker in the context of an EMDE. On the other hand, in order to overcome the dilution, FIMPs would only need to interact stronger to SM particles. This topic of DM production in models with a non-standard cosmology has gained increasing interest lately [8,[14][15][16][17][18][19][20][21][22][23][24][25][26][27][28].
In addition to DM, the neutrino sector constitutes another intriguing piece of our Universe. The discovery of neutrinos' mixing and masses, which are not included in the Standard Model (SM) of Particle Physics, is another reason to think of theories beyond the Standard Model (BSM). Thus, finding a scenario where these two phenomena are connected is an interesting possibility. In the literature, one can find several works that explore the neutrino portal dark matter, in the context of freeze-out [29][30][31][32][33] and freeze-in mechanisms [34][35][36][37][38][39][40]. However, the study of a neutrino portal DM in a non-standard cosmology was only considered in the context of WIMPs [9,41].
In this sense, we will consider in this paper, for the first time, the freeze-in production of dark matter through the neutrino portal, when the Universe was dominated by a matter content at early times. The hidden sector of our model contains a fermion (χ) and a complex scalar (S). Depending on the masses, either χ or S can be a DM candidate with a Z 2 symmetry ensuring the stabilization of the DM candidate. We implement the Type-I Seesaw mechanism for generating the masses of the SM neutrinos, considering three heavy neutrinos, N , which also mediate the visible and the hidden sectors. Therefore, the same couplings that provide neutrino's masses are also responsible for the DM phenomenology. The heavy neutrinos may or may not be coupled to the thermal bath, and we consider both possibilities separately in order to understand the impact of this hypothesis on our parameter space.
In terms of the cosmological evolution of the Universe, we assume that, after inflation, the inflaton decays but there is a matter component whose energy density is larger than the radiation energy density, taking over the evolution of the Universe. In this way, the Universe undergoes an early phase of radiation, followed by a matter-era domination until, at some temperature above the Big Bang Nucleosynthesis (BBN), this matter component decays and reheats the visible sector. We study the evolution of the Hubble parameter in each phase and the impact of this non-standard cosmology on the freeze-in production of the DM candidate.
In this paper, we present a discussion of the dynamics and phenomenology of this dark matter model, exploring the possibility of probing it through both direct and indirect detection. In regards to direct detection, we compare the spin-independent dark matternucleus scattering cross-section with bounds coming from Xenon-1T [42] and projections of Xenon-nT [43]. Indirect signals of WIMPs in the context of neutrino portal involve the internal bremsstrahlung [44][45][46][47] and cascade decays of mediators (see, for instance, [47]). Although challenging, indirect detection of FIMPs is a possibility that has been recently explored [48][49][50]. As an interesting consequence of an EMDE, we explore the sensitivity of the indirect detection of FIMP annihilations to heavy neutrinos, which subsequently decay into SM states. Our work is therefore a contribution to the recent and promising literature regarding the phenomenology of FIMP candidates, which involve constraints from direct [51][52][53] and indirect [48,50] detection searches for dark matter, collider [14,[54][55][56][57] and accelerator [53] experiments and dark matter self-interaction [22,52]. Frozen-in species are also subject to current cosmological bounds [58]. This paper is organized as follows: in Section 2 we introduce the particle content of the model, while in Section 3 we parametrize the EMDE and discuss its impact on the Hubble rate. Then, in Section 4, we discuss the freeze-in production of DM, as well as the theoretical constraints restricting our parameter space. We present the viable parameter space in Section 5 and show the prospects for probing this scenario through direct and indirect detection in Section 6. Finally, in Section 7, we present our conclusions.

The Model
In this work, we consider the Standard Model (SM) content and introduce an extra fermion, χ, and a scalar, S, both comprising the hidden/dark sector. The dark sector particles are odd under a dark Z 2 symmetry, ensuring the stabilization of the fields, while the SM fields are even, and the dark matter candidate will be the lightest particle of the dark sector. In this model, S and χ are not in thermal equilibrium with each other and, for simplicity, we assume that the dark scalar S does not acquire a vacuum expectation value (vev). In addition to this, the model contains three heavy neutrinos, N , responsible for mediating the interactions between the visible and the hidden sectors and generating the neutrino mass through Type-I Seesaw mechanism [59][60][61][62][63]. The Lagrangian of the model is as following: where the L SM corresponds to the SM Lagrangian, and the other terms are given by: corresponding to the Lagrangian of the hidden sector, containing the kinetic and mass terms of χ and S (Eq. (2.2)), the term responsible for the generation of neutrinos masses (Eq. (2.3)) and the term leading to the interactions between the hidden and the visible sectors (Eq. (2.4)). In this model, Y ij ν stands for the Yukawa coupling matrix (the structure will be discussed later), L i are the left-handed leptons, with i = 1, 2, 3 being the generation index, the subscript denotes the interaction basis, m χ , m S and m N are the χ, S and N masses, respectively, andH = iσ 2 H * , with H being the SM Higgs doublet: where v is the Higgs vev, h is the physical Higgs, G ± and G 0 are the Goldstone bosons from symmetry breaking and will become the longitudinal modes of W ± and Z respectively, and H ± and H 0 are the charged and neutral components of the Higgs doublet before EWSB (H 0 is complex). In general, we can include |S| 2 |H| 2 terms in the scalar potential. However, since we want to focus on the neutrino portal, such scalar portal interactions are temporarily ignored, and the exact form of the scalar potential V (S) is, then, irrelevant.

Type-I Seesaw
After electroweak symmetry breaking (EWSB), the Lagrangian in Eq. (2.3) provides masses and mixings for neutrinos through Type-I seesaw mechanism. The information from low energy neutrino measurements (neutrino mixings and masses/mass differences) is embedded in the interaction matrix Y ij ν , which is parameterized in the Casas-Ibarra scheme [64]: where U PMNS is the PMNS matrix containing three mixing angles (θ 12 , θ 23 , θ 13 ) and three phases (δ CP , α 1 , α 2 ) and is parametrized as where c ij ≡ cos θ ij and s ij ≡ sin θ ij , and P = diag(e iα 1 , e iα 2 , 1). The value of these angles and phases are taken from the recent global fitting results [65] 1 . m 1/2 ν/N represent the diagonal matrices with square root of the eigen-masses ( m i ν/N ) in the diagonal entries and R is an extra complex orthogonal matrix (R T R = I) parameterized by three complex angles. In fact, Eq. (2.6) is the generalization of the well-known seesaw formula y ∼ √ mν m N v . In the three generations case, the complex orthogonal matrix R is also important to determine the interaction patterns. If we consider large phases in the complex angles, the interaction can be highly enhanced while still keep neutrino masses light. However, in order to focus on the impact of a non-standard cosmological history, we will only consider the most trivial case where R = I.
Before EWSB, there is no neutrino mixing, and the calculation is done directly with the interaction basis, whereas after EWSB the interaction in Eq. (2.3) will induce the mixing among neutrinos: where ν i /N i are in the interaction basis, ν i m /N i m are in the mass basis and ψ ≡ γ 0 Cψ * is the Lorentz Covariant Conjugate (LCC) in the convention of [66]. The 6 × 6 mixing matrix N is parameterized by with, up to O (m ν /m N ): (2.13) Figure 1. The history of the Universe with an early matter era.

The Early Matter Era
An early matter era is a generic modification to the standard history of the Universe coming from many extensions of the standard model (SM) of particle physics. With such an era, the cosmic history can be simply described as in Fig. 1. After inflation (which ends at the post-inflationary reheat temperature T RH ), the ultra-relativistic SM species produced from inflaton decay dominate the total energy density and we have therefore a radiationdominated (RD) era. Let us assume that a BSM field M was once part of the thermal bath and decoupled while ultra-relativistic, so that its number density remains n M = γ M ζ(3)/π 2 T 3 afterwards. In this case, the ratio between its energy density, ρ M , and the energy density of the bath species, ρ R (T ) = π 2 /30g e (T )T 4 ∝ a −4 , dominated by radiation, would grow with the expansion of the Universe once M becomes non-relativistic (ρ M ∝ a −3 ): where γ M and m M are the internal degrees of freedom and the mass of the BSM field M , respectively, g e is the energetic degrees of freedom parameter, a is the scale factor of the Universe and T is the temperature of the SM thermal bath. When ρ M (T ) > ρ R (T ), at some temperature T i , M starts to dominate the cosmic expansion until some later time (at a temperature T e ), leading the so-called early matterdominated era (EMDE). If M is not to be a cosmic relic, it has to decay completely, and since it was once thermalized with the SM species, it should decay at least partially into the thermal bath. After T e , the decay of M into radiation becomes efficient and makes ρ R to increase, producing entropy until M completely decays, at the so-defined reheat temperature T r . This period is termed as entropy production (EP). After that, the Universe starts to be radiation-dominated again and evolves as in the usual case. Given our current understanding of the cosmic thermal history, the synthesis of light elements, the Big Bang nucleosynthesis (BBN), took place around the MeV scale in a radiation-dominated Universe, and therefore we must ensure that T r 4MeV > T BBN [67][68][69][70].
Since the expansion rate of the Universe depends on the total energy density it contains, the presence of M could in principle affect the evolution of any species through the Hubble rate: with M P 2.4 × 10 18 GeV the reduced Planck mass.
For the purpose of studying the DM phenomenology, it is convenient to know how H(a) depends on temperature. The temperature-scale factor relations during each period are different, as we review in what follows.
During a RD era, the Hubble parameter is the usually considered one: During an isentropic EMDE, the energy density of the non-relativistic matter component M is: where s is the entropy density, g s is the entropic degrees of freedom parameter, Y M ≡ n M /s is the yield of the matter component, the variables with superscript i are evaluated at the beginning of the EMDE, at T i , and Y i M is the initial abundance of the matter component, which is a free parameter in our scenario. With ρ M dominating the energy density, we have: From this expression, it is clear that for temperatures above T r , H EMDE > H RD and, consequently, the Universe expands faster. The EP period is more involved, since in this case both ρ R and ρ M can change significantly. As long as M is much more abundant than the decay products, it is well-known that such an out-of-equilibrium decay can raise the total entropy of the bath species in a comoving volume S = sa 3 [13]: with B R the branching ratio of the decay of M into radiation. From the First Law of Thermodynamics in an expanding Universe, it follows: where w i is the ratio between the energy density and the pressure of a fluid i (w M = 0 for matter and w R = 1/3 for radiation). In Fig. 2, we show the main physics at play in the non-standard thermal history we consider in this work, leaving a more detailed exposition to Appendix B, where the reader can also find the initial conditions we have used. In the left panel, we present the solutions of Eq. (3.6) and Eq. (3.7) for the evolution of SM radiation (magenta) and matter (blue) contents, and of entropy (yellow), as functions of the scale factor normalized to the reheat temperature. When ρ M > ρ R , we have the early matter era, which might finish with a brief period of entropy production. We differentiate between the early radiation-dominated (ERD) era, before the EMDE, and the usual RD one since any relic produced by then might be diluted after the EP period. For simplicity, we assume that the usual radiation era follows just after the EP period, without any other dilution events. In the right panel, we present the evolution of the temperature of the SM bath T (a) (red) and the Hubble rate H(a) (grey). The dotted magenta and blue curves show the Hubble rate containing only radiation and matter, respectively. Details of this figure will be clear in what follows. Although solving the coupled differential equations Eq. (3.6) and Eq. (3.7) is needed for a precise description of the EP period, useful results can still be obtained with some reasonable approximations.
Assuming that the decrease of ρ M is negligible even during the EP ( d(ρ M a 3 ) da ≈ 0) and that M still dominates the energy density, we have: Then, the energy density of the radiation during the EP is solved to be: On the other hand, temperature is defined through the energy density of radiation, ρ(T ) ∝ T 4 . Thus, instead of the usual relation T ∼ a −1 , valid for an isentropic expansion, we have T ∼ a −3/8 while entropy is being produced (see the right panel of Fig. 2). Finally, the Hubble parameter during the EP is: We can further simplify this expression by identifying T r as the temperature at which the decay width Γ M of M is comparable to the Hubble parameter: Γ M = κH RD (T r ). In terms of the dimensionless variable κ, we have: (3.11) As a consequence of the entropy production, the number densities will be diluted. Hence, we would like to identify a variable that can quantify the dilution. From Eq. (3.6), we have: Integrating Eq. (3.12) for a sufficiently long EP period and assuming From the equation above, we can see that the matter component needs to be sufficiently long lived for our purposes. By further using Γ M = κH RD (T r ), we have: (3.14) Thus, we define a dilution factor ∆ as: Notice that the continuity of the Hubble rate among different periods also provides useful qualitative relations. Requiring: Eq. (3.3), Eq. (3.5) and Eq. (3.11) lead to: With these relations, we can fully determine the cosmic history with an EMDE by only specifying T r and T i , with the properties of the matter component M embedded in the latter. In summary, the Hubble rates that we will use in the DM relic density calculation for different periods are:

Freeze-in Production in an Early Matter Era
In the freeze-in mechanism [3,4], the DM interaction strength with the thermal bath fields is so weak that DM cannot reach thermal equilibrium with them. Therefore, if DM particles were not initially present in the SM thermal bath, which might also contain BSM fields, the bath species could be able to produce DM without back reaction. The freeze-in of such FIMP DM candidates would happen whenever kinematically possible: while thermal bath particles are abundant enough (for temperatures above their Boltzmann suppression) and have enough energy to produce FIMPs (for temperatures above DM Boltzmann suppression). As opposed to the freeze-out production of WIMPs, the FIMP relic abundance is proportional to the annihilation cross section or the partial decay width into DM. If we consider the simplest case where a dark matter particle χ is produced through a heavy resonance N decay channel, N → χS, with N coupled to the SM thermal bath, the DM relic abundance from this decay is generically approximated as: where m χ and m N are the masses of χ and N and Γ N →χS is the partial decay width of N into χ and S. Particles which are not coupled to the thermal bath might be much less abundant than the SM species. In this work, we consider both χ and S as FIMP candidates. The mediator N between the FIMPs and the SM fields (Higgs bosons and leptons) is considered separately as thermalized and non-thermalized with the SM bath. As it will become clear in Section 6, this is a crucial point for the phenomenology of our scenario.
In what follows, we present the processes contributing to the production of our FIMP candidates and we give an approximate expression of the relic density of a generic FIMP candidate taking into account an early matter-dominated era. We then discuss the conditions for the heavy neutrinos to be thermalized and the consistency conditions for the freeze-in mechanism to hold. . Processes contributing to the freeze-in production of our dark matter candidate, considered as either χ or S.

Reaction Rate Densities
The evolution of dark matter is governed by the Boltzmann fluid equation for its total number N DM = n DM a 3 : where R DM (t) is the reaction rate density, accounting for all the production and loss of dark matter particles in a comoving volume a 3 . Since we are concentrated here in the freeze-in mechanism, we hereafter regard R DM (T ) as just production rate densities.
In this work, we are going to explore both dark scalar (S) and dark fermion (χ) as dark matter candidates. Whenever kinematically allowed, the channels contributing for their production before EWSB, depicted in Fig. 3, are: After EWSB, we have: Notice that we are neglecting the contributions of the processes S → N χ and χ → N S for the productions of χ and S respectively since they are non-thermal and their initial densities are negligible (n χ , n S n N , n ν i , n H ). For a process 1 → 23, in which species 1 is thermalized with other species than 2 and 3, the rate density for the production of species 2 is approximately given by where Γ 1→23 is the decay width of species 1 into species 2, 3. It is a good approximation to consider the number density of the decaying field n 1 as a Maxwell-Boltzmann distribution: where γ 1 is the internal degrees of freedom of species 1 and K i is the modified Bessel function of second kind and order i.
In the case of the decay N →χS, the rate density is thus given by: where, for convenience, we define the dimensionless parameters For a process 12 → 34, in which species 1 and 2 are thermalized between themselves, the rate density for the production of species 3 is given by: where the approximation holds for initial states with Maxwell-Boltzmann statistics at zero chemical potential. The symmetrization factor S A is the Källen function, |M| 2 is the squared amplitude of the process, E i is the energy of the particle i, f eq i is the phase space equilibrium distribution function of the particle i and s is the center-of-mass energy squared.
The integrated squared amplitudes of the processes of the kind l i H → χS, with an s-channel exchange of N , are all given by: (4.9) When the initial states from the thermal bath have enough energy to produce N onshell, we can use the narrow width approximation (NWA) to have an approximate expression for the rate in the resonant region: (4.10) Far from the resonance, the s-channel processes contribute with the following production rate: where the only approximation is m l m H (T ). Here, we take into account the thermal corrections for the Higgs mass according to the high-temperature expansion approximation of the finite temperature scalar potential [71]. The mass parameter is 16 (3g 2 + g 2 ) + 1 4 y 2 t + 1 2 λ, g and g are the electroweak SU (2) L and hypercharge U (1) Y gauge couplings, y t is the top Yukawa coupling and λ is the Higgs quartic self-coupling. Before EWSB, H 0 and H ± have the same mass which, at high temperature, is approximately √ c h T , whereas after EWSB, the Higgs develops a vev which varies with temperature and G ± and G 0 become the longitudinal modes of W ± and Z, respectively. Then, h, W ± and Z masses evolve with the vev accordingly. We have defined the following integral: accounting for a numerical factor of order one for the most of our parameter space.
In the case of the t-channel, in the limit s , the contribution for the production rate of χ or S reads: where (4.14) and (4.15) Since neutrinos are the only ultra-relativistic species in the freeze-in processes, and their Fermi-Dirac statistics do not provide a significant suppression in the rates, we safely use the Maxwell-Boltzmann approximation in our numerical code.

Contributions to FIMP Relic Density
In terms of the dark matter yield Y DM ≡ N DM /S = n DM /s, the evolution of the dark matter abundance given in Eq. (4.2) reads, in general: The second term in equation above accounts for the dilution in the dark matter abundance due to the entropy production which might take place at the end of an early matter era. In this case, this equation becomes coupled to the set of Eq. (3.6) and Eq. (3.7). As soon as dark matter is produced, its number of particles in a comoving volume, N DM , becomes constant. Assuming no other entropy production period for temperatures below T r , the yield of DM, and therefore its relic density, becomes also constant.
By definition, the relic density of dark matter today reads: and T 0 the present CMB temperature. In the case of the FIMP dark matter candidate, the reaction rate density R DM contains only its production term and does not depend on the dark matter yield itself, as there is no back reaction into the thermal bath. A fair approximate expression for the relic density today, taking into account the early matter era, can be therefore found from the results of Section 3, as we describe in what follows.
In the presence of an EMDE, it is crucial to determine the DM yield accumulated before reheating: is the yield of DM at T r . Under entropy production (from T e to T r ), we can define a different comoving yield given byỸ DM ≡ N DM Φ = n DM Φa −3 , with the nearly constant dimensionless quantity Φ = a 3 ρ M /T r [72]. From Eq. (3.9), we can see that: and therefore its relations with Y χ at T r and T e read: During an EP period, the evolution ofỸ DM is given by [73]: The yield of dark matter today has therefore the following contributions: Notice that all information regarding the specific FIMP model is encoded in the production rate density R DM (T ) 2 . From the expressions above we can extract important information regarding the nature of the freeze-in. If the dominant contribution to the production rate depends on temperature through some power law T n , we can predict whether the freeze-in would happen at the highest or at the lowest scale available of a given cosmological era [75,76]. The kind of freeze-in is therefore referred to as infrared (IR) or ultraviolet (UV) with respect to that era. In the case of a 1 → 2 or a resonant 2 → 2 process, the freeze-in happens at temperatures close to the decaying field or mediator mass. On the other hand, any production channel is only possible for temperatures above the Boltzmann suppression of the heaviest particle involved. In summary, the freeze-in temperature T F I can be determined as follows:  In our analysis, we neglect the subdominant contribution of production during the post-inflationary reheating period (above T RH ), as it would be only important for models featuring R χ ∝ T n with n > 12 or with mediator masses above T RH .
In Fig. 4 we show the solution of the coupled system of Boltzmann fluid equations for the evolution of χ (Eq. (4.16)), the generic matter content M and radiation R, as discussed in Section 3 and detailed in Appendix B. We choose m N > m χ in the left panel and m N < m χ in right panel, with detailed values shown in the plots. We report our results for two possibilities regarding the dilution factor ∆. The case of a short EMDE era with ∆ = ∆ 1 = 2×10 2 is indicated by the dashed grey arrows (for temperatures from T i = 1 GeV to T e = 1.2 × 10 −2 GeV) and corresponds to the dashed black curves, whereas the case of a long EMDE is indicated by the solid grey arrows (for temperatures from T i = 10 14 GeV to T e = 7.6 GeV), with ∆ = ∆ 2 = 2 × 10 16 , corresponding to the solid black curves. In the latter case, we show explicitly the contributions of the decay, s-channel and t-channel processes for the relic density of χ (red, blue and green dotted curves respectively). In both cases, the reheat temperature of the EMDE is set as T r = 4×10 −3 GeV, so that the dilution of the relic density can be observed from the corresponding values of T e until T r . The decay process becomes inefficient at temperatures just below the decaying field mass, so that its contribution for the relic density levels off around T ∼ m N . As we can see, the on-shell production of N from lH annihilations around T ∼ m N lead essentially to the decay contribution. Finally, we can also see that the t-channel becomes ineffective around the Boltzmann suppression of N .

Thermalization of N during the DM Freeze-in Production
Since the freeze-in happens when thermal bath species produce dark matter in out-ofequilibrium processes, it is important to know under which conditions the heavy neutrinos thermalize with the SM particles during the DM freeze-in production.
The chemical equilibrium of heavy neutrinos with the SM bath would have been mainly driven by decays and inverse decays involving Higgs and leptons (N ↔ Hl or H ↔ N l) [77,78], whichever kinematically available, since they are proportional to |Y ij ν | 2 . In our analysis we are not going to consider the subdominant t-channel annihilations N N → ll and N N → HH, both proportional to |Y ij ν | 4 . To determine the thermalization condition, we compare the heavy neutrino or the Higgs decay width with the Hubble parameter in the corresponding phase of the Universe's evolution (Eq. (3.18)-Eq. (3.21)). Thus, when Γ > H(T ) at temperatures relevant to the freeze-in production, the heavy neutrinos are thermalized. As a consequence of the heavy neutrinos being thermalized, all the processes considered in the previous section contribute  for the production of χ or S. Otherwise, the heavy neutrinos are not abundant enough as to effectively decay and annihilate via t-channel into FIMPs and only the s-channel annihilation of Higgs or gauge bosons and leptons can be considered 3 .
In Fig. 5, we show representative cases of the ratios of decay widths Γ N →Hl (T ) (magenta curves) and Γ H→N l (T ) (cyan curves) to the Hubble rate H(T ), as well as the temperaturedependent Higgs mass m H (T ) (yellow curves) as function of the inverse of temperature. Three cosmological scenarios are considered: • without an EMDE (solid curves); 3 In this case, N will be produced via the freeze-in mechanism and then it will contribute to the DM relic abundance through its decay. We temporarily ignore this contribution. Solving the coupled Boltzmann equations is needed for a concrete treatment, which is beyond the scope of this work.
We can recognize in Fig. 5 different slopes in the Γ/H ratios, which are mainly due to the different temperature dependencies of the Hubble rate. Along the solid lines, the Universe is always RD. The abrupt changes in the slope observed in the dashed and dotdashed curves happen once the Universe goes from an EMDE to an EP period, towards the slope expected for a RD era, below the respective T r values. On the other hand, the height of each curve for the same cosmic period depends on the Yukawa coupling Y ij ν , which is proportional to √ m N .
We check whether Γ/H > 1 over temperatures within the vertical gray lines, from max(m N , m S , m χ )/10 to 100×max(m N , m S , m χ ), which is an interval relevant for the DM freeze-in production. Roughly below T ∼ m N , the abundance of N becomes Boltzmannsuppressed. Hence, we set a conservative lower bound on T at m N /10, as indicated by the vertical dotted red line, if the decay of N (and back reaction) is relevant for its thermalization. Once the abundance of the Higgs boson becomes Boltzmann-suppressed, though, scattering with leptons could keep N in the thermal bath. Hence, whether or not N is thermalized during DM freeze-in production highly depends on the relation among m N (which sets the Yukawa coupling), max(m N , m S , m χ ) (which determines the freeze-in production temperature) and T i /T r (which determines the cosmic history). Panels (a) and (b) of Fig. 5 display two limiting cases regarding the mass of the heavy neutrinos: m N m S , m χ and m N m S , m χ , respectively. When m N m S , m χ (panel (a)), the DM freeze-in production happens mainly at temperatures close to m N . A relatively large mass of N will result in a large Yukawa coupling which, in turn, yields a corresponding large decay width. Hence, it is easier for N to thermalize with the cosmic bath and a long enough EMDE is required to avoid the thermalization (dashed and dotdashed curves). On the other hand, when m N m S , m χ , the DM freeze-in production happens at temperatures much higher than m N . In this case, the decay width of the heavy neutrino/Higgs is suppressed by the Yukawa coupling, which hampers the thermalization of N , regardless the duration of the EMDE (panel (b)). However, for light N , relatively small χ and S masses can help achieving the heavy neutrino thermalization during the DM freeze-in production, by lowering the freeze-in temperature (solid curve), as shown in the panel (c) of Fig. 5. In general, a long EMDE hampers the thermalization process, since the Universe expands faster during this period, as we have pointed out before. Nevertheless, choosing a different T r (recall that this is the temperature at the end of the EMDE), the thermalization of N can still be achieved over DM freeze-in production temperatures, as shown by the dot-dashed curve in the panel (c) of Fig. 5.
In summary, for heavy N , the thermalization is easily attained and, therefore, a long EMDE is required to avoid it. In the scenario where N is light, there are two possibilities: if χ and S are heavy, it is hard to achieve thermalization, no matter how long the EMDE lasts. On the other hand, if χ and S are also light, the thermalization is reached even for a long EMDE, as long as it ends around the temperature at which the freeze-in production becomes efficient.  Figure 6. Ratios between the DM production and the Hubble rates for the decay (red), s-channels (blue) and t-channels (green), in the case of small and large entropy production (dashed and solid lines, respectively). In the left (right) panel we have m N > m χ , m S (m N < m χ , m S ). Notice that, the longer the EMDE, the easier to satisfy the out-of-equilibrium conditions, allowing for larger couplings in the freeze-in regime.

Freeze-in Conditions
Let us now consider in more detail the conditions for the freeze-in regime to hold. Besides the assumption of negligible initial abundance, the production of bath species from FIMPs needs to be avoided. To ensure that FIMPs are never as abundant as bath species, our parameter space is such that the reaction rates of FIMP production from the thermal bath fields i (R(T )/n i (T )) were always slower than the cosmic expansion rate (H(T )).
If the heavy neutrinos were always coupled to the thermal bath at the relevant temperatures for the freeze-in processes, all the following conditions need to be satisfied: (4.23) Notice that we will have different conditions for the different eras in the cosmic history, so that we need to consider the corresponding expression for H(T ) at each period (Eq. (3.18)-Eq. (3.21)).
In Fig. 6, we show the ratios between the production rates and the appropriate Hubble expansion rate for the scenario in which m N = 1000 GeV > m χ , m S (left panel) and m N = 1 GeV < m χ , m S (right panel), with m S = 900 GeV and m χ = 50 GeV. In this plot, we have set the reheating temperature to be T r = 4 MeV and considered two possible values of T i : T i = 1 GeV (dashed curves), corresponding to an entropy production of ∆ 2 × 10 2 , and T i = 10 14 GeV (solid curves), corresponding to an entropy production of ∆ 2 × 10 16 . Therefore, in the dashed curves the Hubble rate is dominated by radiation, while in the solid curves it is dominated by matter.
In the case of the decay (red curves), the production rate is given only by Γ N →χS , which does not depend on the temperature. The slightly different slopes we see are only due to the different Hubble rates. The blue curves correspond to the s-channel contributions, so that R χ/S /(n H H) = n H σv /H. Around T ∼ m N , as we have said, N is produced on-shell and n H σv /H ∼ Γ N →χS /H and, for lower temperatures, n H σv /H becomes ineffective due to the Boltzmann suppression on n H . The t-channel contributions shown in green, with R χ/S /(n N H) = n N σv /H, become ineffective due to the Boltzmann suppression of the initial or final states. While in this figure we have integrated Eq. (4.8) without any approximation, Eq. (4.10), Eq. (4.11) and Eq. (4.13) were found to be fair approximations.
In the case where the heavy neutrinos were never coupled to the thermal bath, they were always much less abundant than SM species. Therefore, their decays and t-channel annihilations were always negligible with respect to the SM s-channel annihilations and only the condition [II] of Eq. (4.23) is to be ensured.
In our numerical scans of Section 5 and Section 6, the conditions Eq. (4.23) are always satisfied. As an example of how they constrain our parameter space, let us consider the decay channel [I]. At temperatures T m N /10, the decay channel becomes negligible due to the Boltzmann suppression on the abundance of N . For this reason, the out-of-equilibrium condition for the decay N →χS can be set at the minimum temperature T ≈ m N /10. Considering the parameters of Fig. 6, the out-of-equilibrium condition [I] reads: Thermal equilibrium between χ and S could be achieved via χS → Sχ and χχ → SS, with t-channel exchanges of heavy neutrinos. Since these processes are proportional to λ 4 χ , which we are already constraining with the freeze-in conditions, we can safely assume that our FIMP candidates χ and S do not constitute a decoupled thermal bath.
We can therefore conclude that the possibility of a long EMDE allows for out-ofequilibrium processes with larger couplings. Interestingly, as we are going to see in the next section, larger couplings are needed for keeping the correct value for the FIMP relic density in the case of longer EMDE. Having seen how the early matter era affects the relic density of our dark matter candidate, in the next section we show the parameter space of our model providing the right amount of relic density either for χ or S.

Study of the Parameter Space
In this section, we study the regions of our free parameter space -comprised of m χ , m N , m S , λ χ , T r and T i -providing the correct present relic density of Ω 0 h 2 0.12 [79] for the fermionic dark matter candidate χ. Results for the scalar dark matter do not significantly differ from those shown here. As we have discussed in Section 4, non-thermalized heavy neutrinos are always much less abundant than the SM species and are not able to efficiently produce dark matter. In this work, we roughly take this into account by considering the contributions of decays and t-channels for the evolution of the DM relic density only when the heavy neutrinos thermalize reasonably before the freeze-in time. Otherwise, we consider only the s-channel contributions since SM species are the dominant ones. Carefully taking this issue into account, by tracking the coupled evolution of χ and N , is beyond the scope of this work.
In order to understand the impact of this assumption, we show in Fig. 7 and Fig. 8 the contours of λ χ providing the correct relic density of χ in the plane m χ -m N , considering respectively only s-channels and all the channels. In these plots, we do not check the thermalization of heavy neutrinos nor the freeze-in conditions. We investigate the effect of the duration of the EMDE by considering three different cases: • without an EMDE (gray curves); • EMDE from T i = 10 3 GeV to T e 110 MeV, with T r = 10 MeV (blue curves); • EMDE from T i = 10 14 GeV to T e 24 GeV, with T r = 10 MeV (green curves).
The hierarchy between m N , m S and T i is studied by setting m S = m χ in the left panels and m S = 100m χ in the right panels. The region for m χ > m N is of particular interest since it is when the DM t-channel annihilations relevant for indirect detection searches become possible, as we explore in the next section.
Let us first focus on Fig. 7, where we only consider the contribution from s-channels in the freeze-in process. In order to understand the features in Fig. 7, the approximations for the s-channel reaction rate density in different regions can be helpful. In the resonant region, where the mediator N can be produced on-shell (m N > m S + m χ and m N > m H (T ) + m m h = 125 GeV), from Eq. (4.10), we have: where we have neglected some irrelevant factors related to the phase space and distributions, and Γ N is the total width of N , which depends on both |λ χ | 2 and |Y ij ν | 2 (see Appendix A for a detailed expression of Γ N ). On the other hand, in the non-resonant region, where the mediator N can only be produced off-shell (m N < m S + m χ or m N < m h = 125 GeV), from Eq. (4.11), we have: In the case of a RD Universe, along the gray and the blue lines 4 , the DM yield Y 0 χ in Eq. (4.21) can be estimated by only y RD , integrating from m N (resonant region) or m χ (non-resonant region) to T RH . Therefore, it follows, approximately: where we have used |Y ij ν | ∼ √ m ν m N /v and, in the resonant region, for small enough |λ χ |, From the two approximations in Eq. (5.3), we can see that, for grey and blue lines, in the resonant region, the slope for each equal-λ χ contour will be 1, while in the non-resonant region, it will be 0 (almost no dependence on m χ ). It is also interesting to notice that, in the resonant region, a smaller coupling is required to achieve the correct relic abundance, when compared to the non-resonant region, due to the enhancement of the cross-section.
For much larger T i , as along the green lines, the contribution to the DM yield from the EMDE is dominant. Thus, we have: where, for the resonant region, we have used an unknown function to represent the dependence on λ χ , since, for large λ χ , the NWA is no longer a good approximation. Besides this caveat, from the two approximations in Eq. (5.4), we can see that, for the green lines, in the resonant region, the slope for each equal-λ χ contour will be -2, whereas in the non-resonant region it will be -1/2. Also note that, along the green lines, the continuity of the couplings between resonant region and non-resonant region is also an indication of the failure of NWA. Let us now turn our attention to the case where also decays and t-channels are considered (Fig. 8). While analytic approximations would be much more complicated as the t-channel contribution is more complex, we can take important information from this figure.
We can clearly see that increasing values of T i make the t-channel to play an increasing role in the freeze-in, as λ χ needs to be larger. When m N is heavy enough so that Y ij ν is sizable, the s-channels (∝ |λ χ | 2 |Y ij ν | 2 ) dominate over decays (∝ |λ χ | 2 ) and t-channels (∝ |λ χ | 4 ). We can observe this along the gray and on-shell region of blue contours. For T i = 1 TeV (blue), we see that the t-channel starts to dominate only in the (Y ij ν -suppressed) non-resonant region (its slope is different from the s-channel only cases). In the case of T i = 10 14 GeV (green), the need for much larger couplings makes the t-channel to always dominate in our parameter space. In this case, though, λ χ does not need to be extremely large as in the case where only s-channels contribute to the freeze-in.

Phenomenology
Having discussed the dynamics of the model throughout the modified cosmic history and determined how the dark particle can account for the present dark matter abundance, in this section, we explore different methods to probe the model, which includes searching for signatures in direct and indirect detection experiments. For this purpose, we scan the parameter space for two cases: Case A, where χ is the dark matter; and Case B, where S is the dark matter. The scanned ranges for all input parameters are listed in Tab. 1. Note that the three heavy neutrinos are chosen to be degenerate and the Yukawa interaction matrix, Y ij ν , is fully determined by the heavy neutrino mass and the complex orthogonal matrix R, which is chosen to be the identity I. Additionally, the freeze-in conditions are satisfied. The coupling λ χ is chosen to provide the observed dark matter relic density for each parameter point. As pointed out in the previous section, although a huge DM dilution requires couplings λ χ > O (1), in our plots, we consider only the parameter space where  Table 1. The scan ranges for each input parameter in all cases. Note that Y ij ν is fully determined by m N and R = I, and λ χ is chosen to give the observed dark matter relic density and is required to be less than 4π. λ χ < 4π.

Prospects for Direct Detection
The direct detection experiments aim at identifying the deposited energies when the incident dark matter particle bombards the target nucleus of the detector and, then, is scattered off. In this model, the direct detection relevant vertices (Z − χ − χ and Z − S − S) are induced from the loops shown in Fig. 9. Since the vertex Z − S − S is suppressed by momentum transfer when considered in the context of direct detection, we will not discuss the direct detection approach for Case B. For Case A, the effective Z − χ − χ coupling is expressed as: L ⊃ g Zχχ Z µ χγ µ P L χ, where: where g is the SU (2) L coupling, c W ≡ cos θ W (s W ≡ sin θ W ) is the cosine (sine) of the Weinberg angle, C 0 and C 00 are the loop functions in the convention of LoopTools [81] and 0 − corresponds to the momentum transfer square in the scattering which is negative and close to zero. Note that the finiteness of the above expression is guaranteed by the unitarity of the mixing matrix N. Focusing on the spin-independent (SI) DM-nucleus scattering cross section, we have: where σ SI 0 is the SI DM-nucleus scattering cross section, σ SI χN is the SI DM-nucleon scattering cross section, which is usually used to compare experimental results with different target isotopes, µ χN is the reduced mass between DM and the nucleus (mass number A, charge number Z), µ χp is the reduced mass between DM and the nucleon, i.e. µ χN = mχm N mχ+m N and µ χp = mχmp mχ+mp and g V q are the corresponding vector current couplings of the quarks with Z-boson: The SI cross section σ SI χN is shown in Fig. 10, where each point provides the observed DM relic density in a scenario with (colored points) or without (grey points) EMDE. From Fig. 10, it is clear that the SI cross section is enhanced when the Universe undergoes an early matter era. It is also possible to see that the case where the heavy neutrinos are not thermalized with the cosmic plasma is more promising in terms of the possibility of being detected by direct searches in the future. This is due to the fact that, when the heavy neutrinos are not part of the thermal bath, the s-channel exchange of N is the only mode that produces DM efficiently. Given that this channel is suppressed by the Yukawa couplings, the dark matter coupling λ χ must be larger to achieve the observed DM relic abundance without compromising the freeze-in conditions. However, both cases presented in Fig. 10 cannot be constrained by the most stringent bounds to date from XENON-1T [42] which, together with the projections from XENON-nT [43], are also shown in Fig. 10 for comparison. Even considering the ultimate liquid xenon direct detection experiment DARWIN [82], or the liquid argon direct detection experiment ARGO [83], which are the most promising experiments with a sensitivity for SI cross-sections σ SI χN down to ∼ 10 −49 cm 2 (being able to reach the neutrino floor and to detect or exclude WIMPs with masses above 5 GeV), it would be hard to probe this model in the near future. However, note that, in our analysis, we use R = I as a conservative choice to focus on the effect only coming from cosmology. Other form of R giving different neutrino interaction pattern might improve the sensitivity.

Indirect Detection
As we have seen in Section 6.1, although the perspectives of probing our model through direct detection are slightly better than in a scenario where the freeze-in production occurs in the usual radiation-dominated Universe, it is still hard to find it in the near future (see Fig. 10). The main reason for this is that the DM interaction strength with the SM (via the heavy neutrino) is still very small, suppressed by m ν /m N as well as by loop, which hampers the DM detection in direct experiments.
Nevertheless, space and ground-based telescopes can place stringent bounds on the dark matter annihilation cross-section and lifetime (see [84] for a recent review). In our case, if dark matter is heavier than the heavy neutrino (m DM > m N ), it may annihilate into heavy neutrinos that further decay into SM particles (νh, νZ, l ± W ∓ , l ± W ∓ γ,ννν, νl + l − ) [85][86][87][88], generating neutral and charged cosmic rays after hadronization and parton showers. In particular, over the past years, excesses of gamma-rays has been observed by several experiments, including INTEGRAL/SPI [89], Fermi-LAT [90] and H.E.S.S. [91]. This excess of high-energy photons can be sourced by DM annihilation, and, therefore, it can be a possible DM signature. An interesting feature of the neutrino portal DM is that, although a gamma-ray line signature is not expected, a distinct continuum gamma-ray signal might be present, as studied in [47,86] for the thermal DM case. In this section, we present the annihilation cross-section of DM into heavy neutrinos and show the constraints placed on it by Fermi-LAT and H.E.S.S..
In the non-relativistic regime (s ≈ 4m 2 χ/S (1+v 2 /4), with v the relative velocity between DM particles), the leading order of the annihilation cross-sections for the processesχχ → N N and S * S → N N are given by In Ref. [87], the authors constrain a generic dark matter annihilation cross-section into a pair of right-handed neutrinos in the context of the type-I seesaw, based on the amount of gamma-rays that can be produced via their two and three body decays. The constraints come from the Fermi-LAT and H.E.S.S. observations of Milky Way's dwarf spheroidal Figure 11. The annihilation cross section of the dark matter into heavy neutrino with (color points) and without (grey points) EMD for both Case A (upper panels) and Case B (lower panels). The heavy neutrino is considered to be thermalized (left panels) and non-thermalized (right panels) respectively. The blue and red lines indicate the constraints from Fermi-LAT and H.E.S.S. respectively [87]. galaxies and center, respectively. They have considered one-flavor right-handed neutrinos separately with masses between 10 GeV and 1 TeV. For m χ ∼ 10 GeV, they have found σv 4 × 10 −27 cm 3 /s (Fermi-LAT), while the most stringent constraint for heavier dark matter comes from H.E.S.S., with σv 4 × 10 −26 cm 3 /s for m χ ∼ 2 TeV.
In Fig. 11, we show the annihilation cross section into heavy neutrinos pairs versus dark matter masses for both Case A and Case B with (colored points) and without (grey points) EMDE with both N thermalized (left panels) and non-thermalized (right panels). Similar to the direct detection, with EMDE, larger λ χ is needed to compensate the faster expansion and then to achieve the observed dark matter relic abundance. The largest annihilation cross section that can be achieved with EMDE is roughly 10 −19 cm 3 /s for both Case A and Case B, with N non-thermalized. As pointed out already in Section 6.1, the case where heavy neutrinos are non-thermalized with the cosmic bath provides enhanced relevant crosssections. This happens because the most efficient DM production mode is the exchange of N via the s-channel, which is suppressed by Yukwawa couplings and, therefore, the coupling between N and the DM candidate must be larger to attain the observed DM abundance. The limits from Fermi-LAT (blue line) and H.E.S.S. (red line) [87] for m N = 10 GeV are also shown in Fig. 11 as a guidance. Notice, though, that since all points in this plot are for m N < m χ , we expect comparable bounds on our parameter space for different values of m N . Also, as we discussed in Section 4.3, when N is much lighter than S/χ, the thermalization of N is more difficult. Thus, there are fewer points in the left panels of Fig. 11 than in the right panels. We come therefore to our most important result: reasonably sizable FIMP couplings, allowed by long enough EMDE scenarios which, in turn, avoid the thermalization of heavy neutrinos, can already be tested by the current indirect detection experiments.
Future searches for gamma-ray and cosmic-ray signals with existing (H.E.S.S.-II [92]) and planned (CTA [93], GAMMA-400 [94]) experiments might become able to probe DM annihilations cross-sections into SM particles as low as order 10 −28 cm 3 /s in the light dark matter mass regime [44][45][46][95][96][97][98][99]. Indeed, potential gamma-ray lines from DM annihilation might be finally confirmed or ruled-out [97]. Also, the constraints on DM annihilation strongly depend on the halo profile which is subject to uncertainties. In particular, the existing limit on DM annihilation in the Galactic Center can be as low as 10 −28 cm 3 /s for contracted DM profiles [100].
Interestingly, in the context of an early matter-dominated era, matter perturbations start growing linearly with the scale factor, which in principle have testable consequences [101][102][103][104]. Even a short period of linear growth could allow for perturbation modes to enter the nonlinear regime during the radiation era, leading to the early formation of very small and dense microhalos of DM particles 5 [105][106][107][108]. A long enough EMDE could also lead to the generation of detectable gravitational waves [109,110].
The presence of the EMDE-induced microhalos nowadays would boost the DM annihilation rate relevant for indirect detection searches by many orders of magnitude, depending on the duration of the EMDE and DM free-streaming length. The gamma-ray signal that could come from DM annihilation inside the microhalos is similar to the signal from DM decay, as it depends on the dark matter density instead of its square [105]. As a consequence, the current lower bounds on dark matter lifetime might be translated into upper bounds on dark matter annihilation cross-sections. Considering the Fermi-LAT observations of the isotropic gamma-ray background, the authors of Ref. [111] set upper bounds on dark matter annihilation into bb inside unresolved microhalos of the order σv 10 −32 cm 3 /s, for m DM ∼ 10 GeV. They have also shown that the Fermi-LAT constraints on DM annihilations within the Draco dwarf galaxy are complementary to their conservative bounds, σv 4 × 10 −31 cm 3 /s, and could be used to distinguish them from dark matter decays, in the case of a future discovery.

Conclusions
In this work, we have studied the freeze-in production of FIMP dark matter through the neutrino portal with a modified cosmological history. In terms of the particle physics content of the model, we considered a scenario where, in addition to the SM sector, there is a hidden sector comprised of a fermion, χ, and a complex scalar, S, with the lightest one being the DM candidate. The mass of the SM neutrinos is generated by the Type-I seesaw mechanism and we consider three degenerated heavy neutrinos, N , mediating the interactions between the SM and the hidden sectors. We explored the possibility of the heavy neutrinos being in thermal equilibrium with the cosmic bath, as well as the case where they are non-thermalized, to assess the impact of these hypothesis on the model. In regards to cosmology, we assumed that an EMDE took place and dominated the evolution of the Universe for some period between inflation and BBN, making it to expand faster. Then, this matter component decayed into SM degrees of freedom, reheating the visible sector, and the usual radiation-dominated phase took over. After characterizing the evolution of the Universe with an EMDE for a generic FIMP model, we computed the present DM abundance through the neutrino portal, via the freeze-in mechanism, evaluated the conditions for the freeze-in regime to hold and checked the thermalization condition of the heavy neutrinos.
We found that an EMDE requires larger FIMP couplings to achieve the observed DM relic abundance. The reason is twofold: 1) to overcome a faster expansion, as during an EMDE, as well as an entropy production period, particles need to interact faster; 2) to overcome the dilution due to an EMDE with an initial DM overproduction.
The most important consequence of having larger couplings in the context of a FIMP model is the possibility of making it testable. In this work, we investigated the detectability of our scenario in direct and indirect detection experiments. We have found that an EMDE is able to significantly enhance the cross-sections relevant for both direct and indirect searches. Moreover, when N is not part of the thermal bath, the only efficient production mode of FIMPs is via the s-channel exchange of N in Higgs-lepton annihilations. Since it is suppressed by the Yukawa couplings, the couplings between N and the dark matter candidates need to be larger in order to agree with the relic density constraints while still respecting the freeze-in regime. Therefore, non-thermalized heavy neutrinos manage to enhance even more the relevant cross-sections.
Regarding direct detection, the loop-induced nuclear recoil due to the fermionic FIMP is always much below the current bounds imposed by XENON-1T, as well as the prospected ones from XENON-nT.
As expected for a neutrino portal model, the indirect signals from DM annihilation in dense regions is a better way of probe. In this case, with the larger couplings allowed by an EMDE, our FIMP candidates could efficiently annihilate into heavy neutrinos which then decay and produce cosmic rays. We have shown that a long enough EMDE makes possible to bring the FIMP annihilation cross-section into heavy neutrinos to the current sensitivity of Fermi-LAT and H.E.S.S. experiments. Hence, in the context of an EMDE, motivated by many extensions of the standard model of particle physics, indirect detection experiments can already test FIMPs and, in the case of a detection, it could also provide some hint about the evolution of the Universe.
In summary, the proposed model offers a viable DM candidate whose experimental signatures can be already tested through indirect detection experiments, showing that freezein DM can be probed at present. We should stress that both the EMDE and the nonthermalization of the heavy neutrinos play a key role on our model: in the case where the heavy neutrinos are not in thermal equilibrium with the cosmic plasma and the Universe undergoes an early-matter dominated phase, larger DM couplings are required to attain the observed DM abundance, which translates into a richer phenomenology.
For the s-channel process, the following squared amplitudes contribute: and The factors of 2 account for the contribution of the anti-particles.
Regarding the t-channels contributing for the production of χ and S, we have: . (A.6)

B Evolution of a Matter-radiation System
Here we give more details on the solution of the coupled set of Boltzmann fluid equations Eq. (3.7) and Eq. (3.6), shown in Fig. 2. From aH(a) = da/dt, valid throughout all the expansion, we can translate the evolution over time to the evolution over the scale factor. We re-scale the energy density of the matter component by Φ ≡ ρ M /T r a 3 . We also define the comoving amount of radiation by N R = ρ R a 4 , and use the dimensionless variable A ≡ aT r as evolution parameter. Notice that N M = AΦ. We therefore need to solve: where we define the recurrent constant c 1 = √ 3κ π √ ge(Ar) 3 √ 10 and assume that the production of dark matter from radiation does not change ρ R significantly.
The evolution of the entropy reads: In the case of a Universe filled with radiation and matter, the Hubble rate is explicitly: so that we have H(a) ∝ a −2 when a radiation content dominates and H(a) ∝ a −3/2 when a matter content dominates, as we can see in Fig. 2.
In this work, we assume that the initial conditions for the EMDE, S(T i ) and N R (T i ) are set by a given inflationary model. These values are found after solving the above set of equations for the inflaton-radiation system. From Eq. (B.3), we see that the initial condition for Φ is given by Φ I /A 3 I = 3M 2 P H 2 I /T 4 r . The current observational constraints on the postinflationary reheat temperature poses T RH 7 × 10 15 GeV [112]. In order to explore the scenario in which there was a large entropy production prior BBN, we are interested in the case where T RH > T i T r . As a benchmark value, we consider T RH = 7 × 10 15 GeV and the Hubble rate after inflation H I = 3.3 × 10 14 GeV [113], with A I = 1. We find T MAX 7.9 × 10 15 GeV, S i ≡ S(T i ) = 897 and N i R ≡ N R (T i ) = 1790. It is interesting to notice that the initial conditions fix the amount of entropy which is going to be produced from T e to T r : where we have used Φ = N M m M /T r , in terms of the total number N M = n M a 3 . Choosing ∆ as a free parameter, we can therefore find the initial condition Φ i ≡ Φ(T i ), which sets the moment of ERD-EMDE equality (N i R = N i M ), A i = N i R /Φ i . Let us now find out when the entropy would start to be produced, from the initial conditions. Before the decay of M , while Φ ≈ Φ i , we can solve Eq. (B.1) for N R from A e to a given A: By equaling both terms of the equation above we can find when the entropy will start to be produced: Generalizing the approach used in Ref. [72], the temperature of the thermal bath during an EP period goes as A −3/8 (see Fig. 2) and is approximately given by: As expected from Eq. (3.17), the hierarchy between T e and T r also depends on the initial condition for the matter content: