Freeze-in Dark Matter from Secret Neutrino Interactions

We present a simplified freeze-in dark-matter model in which the dark matter only interacts with the standard model neutrinos via a light scalar. The extremely small coupling for the freeze-in mechanism is naturally realized in several neutrino portal scenarios with the secrete neutrino interactions. We investigate possible evolution history due to the hidden sector with dark temperature: the dark sector would undergo pure freeze-in production if the interactions between the dark-sector particles are negligible, while thermal equilibrium within the dark sector could occur if the reannihilation of the dark matter and scalar mediator is rapid enough. We investigate the relic abundance in the freeze-in and dark freeze-out regimes, calculate evolution of the dark temperature, and study its phenomenological aspects on BBN and CMB constraints, indirect-detection signature, as well as the potential to solve the small scale structure problem.


Introduction
Nowadays, the existence of dark matter has been firmly supported by numerous astrophysical and cosmological observations. However, despite indisputable evidence for its existence, the nature of dark matter remains a mystery. Indeed, all the observational evidence, such as the observations on the galactic rotation curve [1], the bullet cluster [2] and the cosmic microwave background (CMB) [3], depends on the gravitational interaction of dark matter only. Due to the lack of knowledge on its non-gravitational interaction, the properties of dark matter, such as its mass, spin, its couplings to the Standard Model (SM), and whether it is composite or fundamental, remain unknown, nor do we know anything about its production mechanism. Understanding the nature of dark matter would not be so urgent if it plays an insignificant role in the universe. However, dark matter occupies roughly a quarter of the total energy density of the universe today, which is about five times the energy density of the content from the Standard Model (SM), it is the dominating species that governs the expansion of the universe in a large portion of the cosmic history, and it seeds the formation of structure which gives rise to the galaxies and the clusters of galaxies that we observe today. Therefore, resolving the mystery of dark matter is one of the biggest challenges facing particle physics, astrophysics and cosmology today.
Among the numerous attempts to tackle the problem of dark matter, the weakly interacting massive particle (WIMP) is perhaps the most studied dark-matter candidate. In the standard WIMP paradigm, the mass and the interaction strength between the WIMPs and the SM particles are often assumed to be at the weak scale. The WIMPs are assumed to be in thermal equilibrium with the radiation bath in the radiation dominated (RD) epoch, but freeze out of the thermal equilibrium when the interaction rate falls below the expansion rate. It is often referred to as the WIMP miracle that the combination of weak scale mass and coupling strength together with the thermal freeze-out mechanism successfully gives rise to the correct dark-matter relic abundance, Ω DM 0.26. However, aside from its minimality and strong motivation from many models of the physics beyond the Standard Model (BSM), such as theories of supersymmetry, decades of null result in both direct and indirect detections, as well as collider searches, have made the standard WIMP scenarios more and more constrained [4]. Therefore, the current situation has invoked a lot of interest in other dark-matter scenarios.
Indeed, by modifying some of the defining features of the thermal WIMP paradigm, other viable dark-matter scenarios could also be constructed. For example, dark matter need not be produced via the thermal freeze-out mechanism. In fact, the interaction between the dark matter and the SM particles could be very feeble such that thermal equilibrium is never established between the dark matter and the SM thermal bath. Alternative production mechanisms suitable for this kind of scenarios, such as the misalignment production and the freeze-in mechanism [5], are equally motivated. Another possibility is that dark matter could reside in a dark sector (DS) in such a way that its production and/or late-time phenomenology could depend on the internal dynamics within the entire dark sector rather than being solely determined by the interaction between the thermal bath and the dark matter itself. For instance, models from the non-minimal dark-sector framework often depend on the collective behavior of an entire ensemble of dark-sector particles [6][7][8][9]. One could also modify the standard Λ cold dark matter (ΛCDM) cosmology and consider darkmatter production in non-standard cosmologies, such as producing thermal dark matter during a matter dominated universe before the Big-Bang Nucleosynthesis (BBN), see [10,11] for reviews and [12][13][14] for recent progress. By considering different possibilities, the relevant constraints and the potential signatures could be drastically different.
In this paper, we are motivated by the possibility that the interaction between the dark matter and the SM thermal bath is very feeble such that the dark-matter particles could only be populated gradually via the annihilation of the SM particles. To be specific, we consider a hidden sector consisting of a dark scalar mediator φ and a fermionic dark matter χ, and the hidden sector is only connected to the SM via a neutrino portal. Typically, there are two kinds of neutrino-portal dark-matter models -the t-channel models with the ν − χ − φ interaction [15][16][17][18][19][20][21][22][23][24][25][26][27], and the s-channel models with the ν − ν − φ interaction [28,29], in the freeze-out scenarios. While the t-channel models have been investigated in details in recent years, in this paper, we consider the s-channel models in which the dark matter χ only interacts with the SM neutrinos through an s-channel exchange of φ in the freeze-in regime.
The neutrino portal could naturally give rise to a tiny coupling between φ and the SM neutrinos in the low-energy effective theory, which is suitable for the hidden-sector freeze-in production of dark matter. Although the dark matter can never thermalize with the SM thermal bath in freeze-in scenarios, we shall see that the dynamics within the dark sector might be able to significantly affect the dark-matter production. In particular, if the interaction rate between χ and φ is negligible, the dark-matter production depends only on the energy injection from the SM thermal bath. However, if such interaction is sufficiently rapid, the reannihilation of χ and φ could establish a period of thermal equilibrium within the dark sector. In the latter case, the dark-sector particles would eventually break away from the thermal equilibrium either by directly freezing out of the dark thermal equilibrium as their interaction rate drops below the expansion rate, or by entering a phase called the Quasi-Static equilibrium (QSE) in which the annihilation rate of χ is much larger than that of φ and is balanced with the injection from the SM sector before finally freezing out of it. Ultimately, the relic abundance of the dark matter is sensitive to the differences in the evolution history. A sizable interaction strength within the dark sector would allow rapid annihilations of χ inside dark-matter halos at the present epoch. Such late-time annihilations of χ followed by a rapid decay of φ could produce a substantial amount of neutrinos and thus would be constrained by the current neutrino observations. In addition, a sufficiently large coupling between χ and φ could be exploited to solve problems on small scale structure at late time, making our model a potentially viable candidate for the self-interacting dark matter (SIDM, see e.g. Ref. [30] for review and the references in it). This paper is organized as follows. In Sec. 2, we set up the neutrino-portal dark-matter model and discuss several examples for UV completion. In Sec. 3, we study the evolution of the dark sector in different regimes and discuss the consequence of having different evolution history. In Sec. 4, we discuss the phenomenological aspects of our model, including the constraints on the decay of the dark scalar φ, the potential for our model to solve the small scale structure problems, and the indirect-detection constraints. Finally, we conclude in Sec. 5.

The Neutrino Portal Model
In the s-channel neutrino model, we consider a dark sector which consists of a Dirac fermionic χ and a light scalar mediator φ. While the dark-matter candidate χ is directly coupled to φ only, the mediator φ is in addition coupled to the SM neutrinos. To be specific, the effective Lagrangian for a dark sector described above can be in general written as where g χ andḡ χ are the couplings between the mediator φ and the dark matter χ, while g i ν andḡ i ν are the effective couplings of φ to the SM neutrinos with the superscript i being the flavor index. In our model, we do not assume any non-vanishing vacuum expectation value (VEV) for the φ field. Consequently, 1) φ does not mix with the SM Higgs boson and could thus evade all the collider constraints; 2) the only way for the dark sector to interact with the SM sector is through the light mediator φ. However, it is possible that the Higgs VEV might contribute to the mass of φ, which we implicitly include in the m φ term.
To produce dark matter through freeze-in in the early universe, we generically require the couplings g i ν andḡ i ν to be very tiny in order to avoid establishing thermal equilibrium between the dark and the SM sectors. Such small couplings could be naturally generated through appropriate UV models. In what follows, we shall discuss the UV completion for the general Lagrangian above in two limiting scenarios -the one with a pure scalar φ and only pure-scalar interactions; and the one in which φ is a pseudoscalar, and only pseudoscalar interactions exist. We shall refer to these two scenarios as Scenario I and Scenario II. The interacting Lagrangian in Eq. (2.3) can be rewritten by settingḡ χ(ν) = 0 or g χ(ν) = 0 for the two scenarios, respectively: For Scenario I, we are motivated by the type-I seesaw and introduce the right-handed neutrinos N in addition to the scalar φ. We assume the masses M N M EW with M EW being the electroweak scale. At high energy scale, the Lagrangian can be written as Around the electroweak scale, one can integrate out the heavy fields N i and obtain the dimension-6 operators: 1 After electroweak spontaneous symmetry breaking, tree level interactions between the scalar φ and the SM neutrinos can be generated as follows: where v is the SM Higgs VEV. Note that since we assume M N M EW , g ν is naturally suppressed.
For Scenario II, the φ − ν − ν interaction for a pseudo-scalar φ could be generated in the so-called minimal Majoron model [31][32][33][34]. In this type of model, we still introduce three right-handed neutrinos N i but in addition a SM-singlet complex scalar field Φ above the electroweak scale. At high energy scale, the interacting Lagrangian can be written as We then assume that Φ develops a VEV and spontaneously breaks a global U (1) symmetry at a scale higher than M EW . Therefore, we can write Φ = ϕ+v Φ +iφ, where v Φ is the VEV, the pseudoscalar φ is the Nambu-Goldstone boson, and we assume v Φ v. We further assume that the U (1) symmetry is not exact, but explicitly broken by a soft-breaking mass term in the potential of Φ. As a result, the pseudoscalar φ acquires a small mass and becomes a pseudo Nambu-Goldstone boson. After the symmetry breaking, first term in Eq. (2.8) generates the Majorana mass term and the interacting Lagrangian becomes 9) with M N = g Φ v Φ , where we have transform N to the four component notation for a Majorana field such that N c = N . We can then integrate out the heavy filed N , and the effective Lagrangian can be expressed as Finally, after the electroweak spontaneous symmetry breaking, (2.11) Again,ḡ ν is naturally suppressed. We point out that it has also been studied when the pNGB φ becomes a dark-matter candidate in Refs. [35][36][37][38][39][40].
Having demonstrated how Scenario I and Scenario II could be generated from UV theories, in what follows, we shall simply use the effective Lagrangian in Eq. (2.4) when discussing production of dark matter. Note that, in the above discussion, we choose not to specify the coupling strength between χ and φ. As a consequence, while a highly suppressed g ν orḡ ν could ensure a sufficiently feeble interaction between the SM sector and the dark sector during the dark-matter production, by tuning the couplings g χ orḡ χ , we are able to study interesting dark-sector dynamics and phenomenologies in several different regimes. To conclude this section, we list all the Feynman diagrams relevant for our scenarios in FIG. 1.

The Freeze-In Production of Dark Matter
As we have discussed in the previous sections, we are interested in the regime where the interaction between the dark-matter particles and the SM neutrinos is too feeble to ever establish thermal equilibrium between them, and the dark sector is mainly populated by a gradual injection of energy from the SM thermal bath, i.e. via the freeze-in mechanism. Similar to the freeze-out scenario, it is well know that the final yield of dark-matter particles is insensitive to the initial condition in freeze-in scenarios. Therefore, without loss of generality, we shall also assume that the initial energy density of the dark-sector particles is negligible. To be concrete, the evolution of dark-matter number density is determined by the following Boltzmann equation where T refers to the temperature of the SM thermal bath, and T χ,φ are the temperatures that the dark-sector particles χ and φ might have during their evolution which are not necessarily identical. In the equation above, we have implicitly assumed that the phase-space distribution of χ and φ takes the form f χ,φ ∼ e −(E−µ χ,φ )/T χ,φ with their own temperature and chemical potential. Besides, the condition of detailed balance, f eq ν (p ν 1 )f eq ν (p ν 2 ) = f eq χ (p χ 1 )f eq χ (p χ 2 ), is applied to convert n eq ν (T ) 2 σ νν→χχ Tthe term that describes the production of dark matter from neutrino annihilation, into n eq χ (T ) 2 σ χχ→νν T which involves only the number density of χ at the temperature of the SM thermal bath. In general, if the scattering rates are too slow, φ or χ might not be able to reach kinetic equilibrium -they might not have a thermal distribution with a well defined temperature. Nevertheless, we shall keep using the thermally averaged cross-section assuming that it is a good approximation when the typical kinetic energy of the dark-sector particles is T χ,φ .
As it is already discussed in Ref [41][42][43], depending on g χ and/orḡ χ and the darksector number densities, the production of dark-sector particles can be categorized into two different regimes: • The pure freeze-in regime, where the couplings between χ and φ are so small that interactions between them are negligible. As a result, χ and φ can never reach thermal equilibrium.
• The reannihilation regime, where the couplings g χ and/orḡ χ are large enough that the annihilations process χχ ↔ φφ starts to play a significant role in the evolution The classification of production mechanism for our freeze-in dark matter model. g χ andḡ ν are the pseudo-scalar couplings of the mediator φ to the dark matter particle χ and the SM neutrino respectively.
of n χ and n φ . This could also allow χ and φ to reach thermal equilibrium during a certain period of time in the evolution of the universe with a common temperature denoted as T D which is smaller than the temperature of the SM thermal bath T .
In what follows, we describe these two regimes in detail while focusing on Scenario II in which φ is a pseudoscalar case, as the discussion can be easily carried over to Scenario I in which φ is a pure scalar. We also associate the FIG. 2 for a summary for the different regime with the corresponding scale of the relevant couplings.

Pure Freeze-in Regime
As mentioned above, in the pure freeze-in regime, not only the dark matter is decoupled from the SM thermal bath, the thermal equilibrium between dark-sector constituents χ and φ is also never reached. Therefore, during the production of dark-matter particles χ, one can make the following approximations: 1) the number density n χ is so small that the annihilation process χχ → νν can be safely neglected; 2) the particle exchange between dark-sector species φφ ↔ χχ is also negligible. With these approximations, the Boltzmann equation Eq. (3.1) can be approximated as where T represents the temperature of the SM thermal bath, and we have used the condition n eq χ (T ) n χ (T ). To connect with the observed dark-matter relic density today, it is often convenient to use entropy conservationsa 3 = const. wheres is the entropy density and a is the scale factor, and define the comoving number density or the yield of dark matter Y χ ≡ n χ /s which eventually becomes a constant after the freeze-in process concludes. Therefore, to reproduce the expected present-day dark-matter relic density, the final yield The above expression can thus be used to constrain our model parameters. In what follows, we shall estimate the yield from the pure freeze-in process. We first note that where Γ φ is the decay width of φ, s is the Mandelstam variable, and E χ 1 , E χ 2 represent the energy of the dark-matter particles in the annihilation process. In the high-temperature regime (T m χ ), most of the annihilations take place with s 4m 2 χ , whereas in the low-temperature regime (T m χ ), typically s ∼ 4m 2 χ . Taking these approximations, the thermally averaged cross-section can be expressed as follows assuming m χ m φ : In the pure freeze-in regime, since the νν → χχ process is the only production channel for the dark-matter particles , the condition that the dark matter never reaches thermal equilibrium with the SM generally requires that where H is the Hubble expansion rate. In the radiation dominated epoch, we have where g is the effective numbers of relativistic degrees of freedom for the energy density, and M P is the reduced Planck mass. Since the number density in Eq. (3.5) is always exponentially suppressed in the non-relativistic regime, we simply require Eq. (3.6) to be satisfied in the relativistic regime, which implies To proceed further with the estimate, we rewrite the Eq. (3.2) as which we can directly integrate from a large initial temperature T i m χ to the presentday temperature T now : Using our previous discussion, the integrand Obviously, the Boltzmann suppression factor which appears in the non-relativistic regime suggests that the contribution to the entire integral from this part, i.e. m χ > T ≥ T now , is always subleading. At the same time, the 1/T 2 behavior in the relativistic regime also indicates that the production of dark matter is very small when T m χ . In fact, most of the contribution to the integral in Eq. (3.10) actually comes from the part where T ∼ m χ . Therefore, we find as expected that the dark-matter production in the pure freeze-in regime is largely insensitive to the initial conditions.
Based on our discussions, we can then make following approximations and roughly estimate the final yield: in which g ,s is the effective numbers of relativistic degrees of freedom for entropy density. x ≡ m χ /T . In the left panels, we show a comparison between the expansion rate H and the dark-matter production rate n eq χ (T ) σ χχ→νν v T . As expected, the dark-matter production rate is always smaller than H. However, despite the fact that it is always decreasing as the SM temperature drops, it is closest to H at T ∼ m χ , which supports the approximation we use in Eq. (3.12). After the temperature drops below m χ , a quick down turn appears as the neutrinos no longer have enough energy to create the dark-matter particles. As a result, we see in the right panel that the yield grows monotonically from a vanishingly small initial value, and then asymptotes to its maximum near T ∼ m χ . Note that the sudden jump for both H and Y χ around x ∼ 5 − 10 is not a numerical artifact but merely due to the sudden change in g and g ,s when the temperature drops to T ∼ O(100) MeV, i.e., when the QCD phase transition occurs.
Comparing the pseudoscalar with the pure scalar cases, it is obvious that for couplings of the same size, the final yield in the latter is slightly smaller. This is simply because the factor s 2 in the numerator of the thermally averaged cross-section in Eq. (3.4) becomes s(s − 4m 2 χ ) in Scenario I (see Eq. (A.1) for more details), which results in a suppression in the production rate as T approaches m χ . [GeV] Examples in the pure freeze-in regime. The upper panels correspond to Scenario II where φ is a pseudoscalar, and the lower panels correspond to Scenario I where φ is a pure scalar. In the left panels, the blue curve shows the evolution of the dark-matter production rate n eq χ (T ) σ χχ→νν v T , whereas the dash-dotted purple line is the decay rate of φ. The dashed red curve is the evolution of the Hubble expansion rate H. In the right panels, we show the evolution of Y χ as a function of x ≡ m χ /T . The dotted horizontal line indicates the observed relic density of the dark matter today.

Reannihilation Regime
In the pure freeze-in regime, the interactions between the dark-sector particles are so weak such that both χ and φ barely annihilate. As a consequence, the χχ ↔ φφ process has little impact on the evolution of dark-matter number density. However, it is also possible that while the entire dark sector is decoupled from the SM sector, the annihilations of χ and φ proceed with a non-negligible rate. Indeed, since the interaction strength between χ and φ is governed by the couplingḡ χ in Scenario II (or g χ in Scenario I), by increasinḡ g χ while adjustingḡ ν accordingly in order to prevent thermalizing with the SM sector, as required in Eq. (3.8), the dark sector could establish its own thermal equilibrium for a period of time in the history of the universe as it is populated by the SM thermal bath. This thermal equilibrium is characterized by a common temperature T D for both χ and φ. In this situation, Eq. (3.1) can be rewritten as (3.13) in which the the thermally averaged cross-section for χχ ↔ φφ is evaluated at T D . The thermal equilibrium in the dark sector requires the interaction rate between the dark-sector particles n eq χ (T D ) σ χχ→φφ T D > H(T ). For a sizable cross-section σ χχ→φφ , this condition could be immediately realized if the initial particle number density n χ or n φ is large enough. By contrast, if the initial number density is small, the dark sector could start without thermal equilibrium but slowly build it up as the SM sector keeps populating the dark-sector particles. Once the thermal equilibrium is established within the dark sector, one would expect the number density n χ to track n eq χ (T D ), since the last term in the Boltzmann equation above always tends to bring back any deviation from the thermal equilibrium. This is indeed true if the second term on the right-hand side of Eq. (3.13) is not large enough such that the χχ → φφ process can always be balanced by the φφ → χχ process. However, if this SM source term is sufficiently large, n χ might be able to grow significant larger than n eq χ (T D ). Therefore, the annihilation rate of χ could overwhelm that of φ, and this could occur even when the latter is still larger than the Hubble rate. In this situation, Eq. (3.13) can be approximated as (3.14) The equation above simply means that while the SM thermal bath is always sourcing the dark matter χ, the annihilation of χ always tends to deplete this particle injection. As a result, the QSE phase often occurs in which the last two terms in Eq. (3.13) are balanced against each other. To be precise, the QSE phase is defined by The dark thermal equilibrium or the QSE persists until n χ (T ) σ χχ→φφ T D H(T ), and then the dark matter experiences a dark freeze-out. However, unlike the standard thermal freeze-out scenario, the comoving number density of the dark matter might not be fixed at its value around the dark freeze-out, because the production of χ from the SM thermal bath could still make a significant contribution to the final yield even after the dark freeze-out. In any case, differences caused by having different initial number densities for χ or φ would eventually be washed out by the thermalization within the dark sector -the freeze-in mechanism is also independent of initial conditions in the reannihilation regime.
As mentioned above, Eq. (3.13) depends on the dark-sector temperature T D , which, unlike the SM temperature, is not known a priori. In order to find T D , we numerically solve the Boltzmann equation for the total energy density in the dark sector: where the energy density ρ D ≡ ρ χ + ρ φ , the pressure P D ≡ P χ + P φ , and the thermally averaged energy transfer rate P a···→i... gives the energy transferred to or from the dark sector through the process a + · · · → i + . . . . We shall show its explicit form in Appendix B. Note that the factor of 2 in the second and third term on the right-hand side takes into account the fact that two dark-sector particles are produced via each of the corresponding processes. Again, we have made the assumption that the annihilation rate of the dark-sector particles into the SM neutrinos is negligible. Such assumption is valid as long as T D T , which implies the number density of either χ or φ is much smaller than that of the SM neutrinos, and is also consistent with the second law of thermodynamics -the dark sector can never be hotter than the SM thermal bath which populates it in the first place. We shall see later that the decay and inverse decay could potentially thermalize φ with the SM neutrinos. However, this has little impact for the yield of χ, if any.
While in thermal equilibrium, the dark-sector energy density is only a function of the dark-sector temperature T D once the masses of χ and φ are fixed. Therefore, at every step of the numerical simulation, T D can be obtained by solving in which ξ χ = 2 and ξ φ = 1 are the number of degrees of freedom for χ and φ. After obtaining the time dependence of the dark sector temperature T D (t), one can then insert it back to Eq. (3.13) and solve for n χ . In FIG. 4, we present our benchmark results in the reannihilation regime including the evolution of the energy-transfer rates, the particle exchange rates, as well as the dark-matter yield. The upper two panels show the cases in which the dark sector undergoes a period of dark thermal equilibrium and then freezes out. By contrast, the cases with the additional QSE phase before the dark freeze-out are shown in the lower panels. Comparing with the pure freeze-in cases in FIG. 3, we have enhanced the coupling of φ to the dark matter, and at the same time lowered its coupling to the neutrinos in order to avoid thermalization with the SM thermal bath. The couplings are tuned such that the cases in Scenario II could have a final yield of χ similar to that of the dark matter today. At the same time, for cases with (or without) the QSE phase, we choose to use couplings of the same size just to demonstrate that the cases in Scenario I tend to have a smaller final yield due to the same reason that we have mentioned in the pure freeze-in cases.
In the left panels, we see for all the cases, most of the energy injected into the dark sector is from the νν → χχ channel at early times (solid blue curve). However, the contribution from this channel decreases drastically after the SM temperature drops below the mass of the dark-matter particles, as we can tell from the quick down turn of the solid blue curve for x 1, and then the inverse decay νν → φ (solid green curve) starts to dominate. On the other hand, the energy transfer rate of the annihilation process νν → φφ (solid orange curve), is always negligibly small. These observations are consistent with the expectations from our scenario in whichḡ ν (or g ν ) is very tiny, since the amplitude of the νν → φφ process is proportional toḡ 2 ν (or g 2 ν ), whereas the amplitudes of the νν → χχ and νν → φ processes are proportional toḡ νḡφ andḡ ν (or g ν g φ and g ν ), respectively.
Comparing the cases with and without the QSE phase, we notice that because both the coupling of φ to the dark matter and to the neutrinos are larger in the cases with the QSE phases, more energy is injected from the SM in these cases. As a result, in the showing the cases without (with) the QSE phase. The left columns shows the rate of the energy transfer from the SM sector as well as the rates of energy-density dilution due to the expansion of the universe. In the middle column, we compare the expansion rate with the interaction rates relevant for the dark-matter number density. In the right column, we show the evolution of Y χ and Y eq χ (T D ) -the yield of χ if it is in thermal equilibrium with φ, which we have colored to indicate the corresponding T D . The evolution of Y QSE χ (T )) is also plotted in the lower two panels for the cases with the QSE phase. middle and the right panels, while the upper two cases only reach thermal equilibrium around x ∼ 10 −3 , the lower two cases have already established thermal equilibrium long before that, and have acquired much more dark-matter particles when the dark sector is in thermal equilibrium. As a consequence, the evolution after reaching thermal equilibrium is also drastically different.
In the first two cases, Y χ tracks Y eq χ (T D ) until n eq χ (T D ) σ χχ→φφ T D H, which occurs around x ∼ 0.02 − 0.05. At this time, the temperature in the dark sector T D is already smaller than m χ as one can tell from the color variation in Y eq χ (T D ). As we have mentioned above, after the dark freeze-out, the yield Y χ is not fixed yet. Similar to the pure freeze-in scenario, the solid blue curve in the middle panel shows that the production of χ from the SM neutrinos is still significant until the SM temperature is much smaller than m χ . Therefore, with the injection from the SM thermal bath continues to contribute, the yield only approaches its asymptote after x 1. In the bottom two cases, Y χ tracks Y eq χ (T D ) at the beginning, and can even grow larger than its final asymptotic value because of a largerḡ ν (or g ν ). These "excess" dark-matter particles are rapidly depleted by annihilating into φ. As a result, as T D approaches and goes below m χ , Y χ grows increasingly slower and eventually starts to decrease (see the green part of Y eq χ (T D )). Note that, due to a larger coupling between χ and φ, Y χ can still track Y eq χ (T D ) even when χ is already non-relativistic. As T D keeps dropping, the inverse process φφ → χχ becomes less and less important due to kinematics. However, sinceḡ ν (or g ν ) is larger in these cases, the dark-matter annihilation rate n χ σ χχ→φφ v T D can still be larger than H. This is clearly shown in the middle panels as the solid brown curve peels off from the dash-doted black curve. Consequently, the dark matter breaks away from the thermal equilibrium and enters the QSE phase in which Y χ tracks Y QSE χ (T ). Finally, when n χ (T ) σ χχ→φφ v T D H(T ), the QSE can no longer be maintained. Then, χ freezes out, and the yield approaches its final value.

T D Evolution
We comment here on the evolution of dark-sector temperature T D in the the reannihilation regime. The curves for T D which are obtained by solving Eq. (3.17) are shown explicitly in FIG. 5. Note that, we only include two curves for T D because there is no visible distinction between the scalar and pseudoscalar cases. It is worth emphasizing that, though the solution to Eq. (3.17) is unique and exists in the entire regime we are considering, T D is only physically meaningful when the dark-sector particles χ and φ are in thermal equilibrium. We therefore shade the region in which the dark thermal equilibrium is reached. Within this region, we can see that slopes of both T D curves are less steep than the that of T . This is because while the energy density of the SM thermal bath is only subject to the expansion of the universe and some minor impact from the change of g , the dark-sector energy density is in addition always sourced by the SM. Therefore, its temperature T D decreases slower than T . After χ and φ decouple, the momentum distributions of χ and φ evolve independently. As we can see from the left panels of FIG. 4, the inverse decay of neutrinos into φ starts to surpass the annihilation into a pair of χ around x ∼ 0.1, which is shortly before or after φ and χ decouple. After that, most of the energy injection from the SM thermal bath goes The yellow curve is obtained from the cases with the QSE phase, whereas the blue curve is from cases without it, and there is no visible distinction no matter φ is a scalar or pseudoscalar. Since the temperature in the dark sector is only physical when χ and φ are in thermal equilibrium, we roughly shaded the region in which the dark thermal equilibrium is established. As a reference, we use the red curve to show the temperature of the SM thermal bath.
into φ. However, the energy that goes into φ and cannot be rapidly converted into χ since χ is already non-relativistic when it decouples. The energy density in the dark sector is thus dominated by φ after decoupling. Therefore, assuming that the kinetic equilibrium is still maintained for φ, and that chemical potential is negligible, the T D curves to the right of the shaded region largely reflect the evolution of the energy density of φ, i.e. T D ≈ T φ since ρ D ≈ ρ eq φ (T D ). However, T D always tends to be an overestimate for T φ in this region as ρ χ > ρ eq χ (T D ) after χ decouple. Eventually, when H(T ) ∼ Γ φ , φ decays rapidly into neutrinos. For the choice of m φ andḡ ν (or g ν ) in the examples, the decay occurs around x ∼ 10 − 20, when φ starts to transition from relativistic to non-relativistic. Therefore, though both T D curves eventually merge with T , which indicates that φ could thermalizes with the SM sector through decay and inverse decay, its energy density is already suppressed at that time -it is only the neutrinos living on the high momentum tail of the Fermi-Dirac distribution that are able to create the φ particles, which, in turn, quickly decay and deposit its energy back to the SM thermal bath.

Model Phenomenologies
In this section, we discuss general considerations on the phenomenological aspects of our model. We notice that, at early times, the existence of the dark scalar φ might contribute to a non-negligible fraction of energy density in the universe, and its decay would inject energy into the SM thermal bath. Therefore, our model could be constrained by measurements on BBN and CMB. At late times, the interaction within the dark sector could allow a large annihilation rate inside dark-matter halos which would produce a lot of neutrinos. Thus our model also receives constraints from the current observations on astrophysical neutrinos. Moreover, the self-interaction between dark-matter particles could also suppress the formation of small structures. As a result, our model could potentially be exploited to solve the small scale structure problems. In what follows, we discuss these aspects in turn based on the general effective Lagrangian in Eq. (2.3) instead of referring to any specific UV completion.

Dark Scalar Decay
The existence of φ and its decay might lead to measurable effects by changing the expansion rate of the universe and injecting SM neutrinos to the visible sector.
On one hand, φ could decay out of equilibrium, i.e. Γ φ /H(m φ ) < 1 (or equivalently T dec < m φ ), where T dec is the SM temperature when H = Γ φ . If this occurs before BBN, it simply reheats the universe by producing entropy into the visible sector since neutrinos are still tightly coupled to the SM plasma. Therefore, other than slightly changing the temperature evolution of the SM thermal bath and diluting the relic density of χ, the decay of φ would not leave any measurable effect on BBN, even if φ has a significant amount of energy density before decaying. However, if such out-of-equilibrium decay occurs during BBN but before the neutrino decoupling, it could potentially disturb the formation of light elements by changing the way the SM temperature evolves. In addition, if φ decays after the neutrino decoupling, it could contribute to the effective number of neutrino species N eff which can be constrained by CMB measurements as ∆N eff ≡ N eff − 3.045 [44]. Therefore, to avoid upsetting BBN and CMB, we basically need φ and its decay products to have a negligible contribution to the total energy density of the universe when T 10MeV. On the other hand, if Γ φ /H(m φ ) > 1 (or T dec > m φ ), the decay and inverse decay would bring φ into thermal equilibrium with the SM neutrinos. If this occurs before neutrino decoupling, φ then becomes an additional component of the visible plasma and could contribute as much as one relativistic degree of freedom to the total energy density when it is relativistic. This could lead to significant increase of the expansion rate during BBN, which would affect the prediction for the abundances of light elements. If φ decays in equilibrium after neutrino decoupling, the decay products could then lead to a measurable ∆N eff in CMB. It could also affect BBN if the production of light elements is not concluded yet. Therefore, observational constraints generally require φ to be non-relativistic when BBN starts, i.e. m φ 10 MeV. In our model, using Eq. (A.7), the SM temperature at which the decay occurs is roughly estimated to be MeV . MeV allows φ to be non-relativistic when BBN starts, therefore its energy density is already negligible. For the latter cases, the fact that T dec > 10 MeV ensures that no BBN or CMB bounds would be violated.

Indirect Detection
The decay of φ could also be exploited in the indirect detection of dark matter. In the late-time universe, though the interaction rate between the dark matter χ the SM particle is extremely small, χχ → φφ followed by rapid decay of φ → νν could potentially produce a signal that can be probed by astrophysical neutrino observation experiments. These experiments usually give bounds on the dark matter annihilation cross-section σv . In the non-relativistic regime, σv can be expanded into powers of the relative velocity v: where the angle brackets represent the average over the local dark-matter distribution at the signal source. In our model a and b can be expressed as: According to the result in Ref. [45], we find that the most stringent bound on the annihilation cross-section σ χχ→νν v for m χ ∼ 1 GeV comes from Super-Kamiokande (SK) experiment [46], which translates into a 10 −24 cm 3 · s −1 and b 10 −16 cm 3 · s −1 . We plot this constraint againstḡ χ and g χ in FIG. 6 where the gray region is excluded by the SK experiment. Our benchmark points in FIG. 4 all seat infinitely far away along the four colored lines which cannot be shown on the log-scale plot as they all have one of the neutrino couplings g χ orḡ χ being zero. Nevertheless, our benchmarks all survive from the indirect detection constraint.
FIG. 6. We plot the indirect detection constraint from the Super-Kamiokande. The excluded region is colored in gray. Different colored lines are related to the four benchmark points in FIG. 4. As the benchmark all have one of the neutrino couplings vanished, they seat at infinitely negative g χ orḡ χ along the lines in this log-scale plot.

Possible Solution to the Small Scale Structure Problem
Finally, we investigate the potential for our model to solve problems on small scale structures. It is well known that though the standard ΛCDM cosmology is very successful in predicting the large scale structure of the universe, it fails to predict correct structures that are compatible with the observational data at the small scales. The most prominent small-scale problems are the follows: 1) the core-cusp problems -the dark-matter density profile in the center of a galaxy predicted from simulations with collisionless cold dark matter (CDM) typically scales as r −1 ("cuspy"), while the observations show that both cored (dark matter density ∼ r 0 ) and cuspy centers exist [47]; 2) the missing satellites problem -the observed number of satellite galaxies in the local group is one to two orders of magnitude fewer than the prediction from CDM simulations; the too-big-to-fail problemwhen matching the most luminous observed satellites of the Milky Way with the most massive subhalos from CDM simulations, the subhalos appear to be too massive to host those observed galaxies 2 . These problems motivate the consideration on the self-interaction of  dark-matter particles which could potentially solve the problem by lowering the density in the central region of a galaxy as compared to the CDM simulations. An excellent review on this topic is provided by the authors of Ref. [30], where a rough bound on the dark-matter self-scattering cross-section over its mass σ/m χ can be derived by fitting the astrophysical data. As pointed out in Ref. [30], σ/m χ ∼ O(1) cm 2 /g can solve the small-scale problems, and a transition in the scaling of the velocity dependence of the self-interacting cross-section around O(10 3 ) km/s is crucial for reconciling the collisionless nature of the dark matter on large scales with the collisional behavior on small scales. On the other hand, Ref. [49] points out that the viscosity cross section is a better proxy than the center-of-mass cross-section σ in characterizing the self-interacting strength of the dark matter when dealing with small scale structure problems as it emphasizes the momentum transfer in the perpendicular direction and better describes the speed of energy equalization. Therefore, we calculate σ V /m χ in our model and find that, at zero velocity, to leading order in the expansion with respect to m φ /m χ Notice that it depends on g χ only as the terms containingḡ χ are suppressed by at least For each m φ , we further select a value for g χ to ensure that σ V /m χ ∼ O(1) cm 2 /g at low velocity in order to solve the small scale structure problem, and two extreme values forḡ χ to demonstrate its weak effect on σ V . Obviously, for the benchmark masses we select in FIG. 4, O(1) value is needed for g χ in order to solve the small scale problems while roughly satisfying the constraints on large scales (e.g. the Bullet Cluster constraint). However, a g χ of this size would lead to a long period of dark thermal equilibrium or QSE, making the yield of dark matter too low at the dark freeze-out. In addition, since g ν and/orḡ ν are bounded from above by the requirement on freeze-in production, increasing g ν and/orḡ ν could not solve the problem with our choice of mass parameters. On the other hand, decreasing m φ to ∼ 10 MeV could allow a smaller g χ and is thus able to solve the small scale problems and obtain the correct relic abundance simultaneously. However, for suitable values of g ν andḡ ν (roughly the same order as the values shown in FIG. 3 and FIG. 4 since it is not sensitive to m φ ) that yields the correct relic abundance, the decay of φ would violate the BBN constraint by presenting a sizable ∆N eff . Thus, we find the only way to make our model a potentially viable SIDM candidate while satisfying the constraints on the relic abundance and BBN is to increase m χ while adjusting g χ accordingly to maintain an appropriate σ V /m χ , since a larger m χ would enable χ to freeze out earlier before getting severely Boltzmann suppressed.
In general, we see that the conditions for solving the small scale structure problems (requiring g χ ∼ O(1) and preferring lighter m φ ) tend to act against the requirement to produce the correct dark-matter relic abundance in the freeze-in regime (favoring smaller g χ to prevent a late dark freeze-out) and the bounds on the decay of φ (favoring m φ 10 MeV).

Conclusion
In this paper, we study a simplified model describing the neutrino-portal dark matter χ.
In this model, the dark matter χ does not couple directly to the SM particles. Instead it interacts with the SM neutrinos via an s-channel exchange of the light scalar mediator φ. Such neutrino-portal interaction could naturally give rise to a tiny coupling between the dark and the SM sector by relating the smallness of the coupling to the smallness of the neutrino mass, and therefore is appropriate for the freeze-in production of dark matter which requires an extremely small interaction rate between the two sectors. We then study the possible UV completions of our model and the production of dark matter in the early universe focusing on interesting dynamics that could arise in the dark sector. We also investigate several potential phenomenological constraints associated with our model. Our main results are summarized in the following: • We point out two UV origins of the simplified model, namely, the type-I seesaw and the Majoron model, that can naturally induce a small pure scalar and pure pseudoscalar coupling between φ and neutrinos, respectively.
• We study the freeze-in production of χ in the pure freeze-in regime as well as in the reannihilation regime. For the former, we analytically estimate the relic abundance and present two benchmark results in FIG. 3 by solving the Boltzmann equation for the number density n χ numerically. For the latter, we study the dark-sector dynamics in detail and further identify two different scenarios -the one with the dark thermal equilibrium only and the one with the additional QSE phase. Typically, the one with the QSE phase requires a larger rate of injection from the SM sector (larger g ν and/or g ν ). We find that since the Boltzmann equation for the n χ depends on the dark-sector temperature T D , we need to solve the Boltzmann equation for the dark-sector energy density ρ D in order to obtain the evolution of T D , and use T D to solve the Boltzmann equation for n χ . We present the other four benchmark results which include both scenarios for both the scalar and the pseudoscalar case. We also present and analyze the evolution of T D in these cases.
• We find that the decay of φ could potentially be constrained by BBN and CMB. Nevertheless, we point out that, for our benchmark results, a relatively heavy mediator with m φ ∼ 50 MeV could avoid the BBN and CMB constraints while simultaneously give the correct relic abundance.
• The Super-Kamiokande experiment is possible to constrain our model by searching for dark matter annihilating to mediator φ followed by φ decaying to SM neutrinos. We show in FIG. 6 that our benchmarks are all allowed by the Super-Kamiokande experiment.
• Since our model could allow sufficiently large self-interaction between dark-matter particles, we analyze the potential for our model to solve the small scale structure problems. We find since there is a tension between the required properties for a successful SIDM candidate and the constraints on both BBN and the relic abundance, the choice of mass parameters that we choose in the benchmarks cannot solve the problems on small scale structures. Nevertheless, it is possible to reconcile the problem by increasing the mass of χ.
Finally, we would like to point out several caveats. First of all, in our examples for UV completion, we have not considered the possibility that the heavy right-handed neutrino N could also be a dark-matter candidate. Such possibility in general depends on the evolution history of the universe at high scales. We argue that, in scenarios in which the temperature of the universe has never been higher than M N , or the mass of the inflaton is smaller than M N , one would not expect a non-negligible population of N . Our analysis thus implicitly assumed an appropriate cosmological history before the scales relevant for the production of dark matter. Scenarios in which both N and χ could play the role of dark matter are nevertheless worth studying.
Second, although we have confined ourselves within the freeze-in production of dark matter through the neutrino portal, if the couplings g ν and/orḡ ν are large enough, the dark-matter particles could thermalize with the SM neutrinos and undergo thermal freezeout from the SM thermal bath. However, we do not consider this scenario in our paper as we are motivated by the possibility that the interactions between the dark matter and the SM particles are extremely feeble.
Third, as we have pointed out before, in the reannihilation regime, the interaction rate between the dark-sector particles depends on the dark-sector temperature T D . However, T D is only physical when the χ and φ particles are in thermal equilibrium. In our analysis, we obtain T D from the energy density of the dark sector and use it to evaluate the dark-sector interaction rate throughout. Therefore, our approach relies on the assumption that the thermally averaged cross-sections evaluated at T D is a good approximation for the actual interaction strength near but beyond the dark freeze-out. Far away from the dark freezeout, the interaction within the dark sector becomes negligible, and the evolution of the dark matter number density is dominated by the injection from the SM sector and the Hubble expansion.
Furthermore, when calculating the thermally averaged cross-sections and the energy transfer rates, we have assumed that the dark-sector particles χ and φ all follow a thermal distribution. While this is true when the dark sector is in thermal equilibrium, it might not be true if the dark sector is not, as the self-scattering rate might be too small to reach kinetic equilibrium. Therefore, it is possible for the phase-space distributions of χ and φ to be highly nonthermal as the SM thermal bath keeps injecting particles into the dark sector while the existing dark-sector particles are constantly redshifting. This would in principle affect the calculation of the collision terms. The nonthermal features in the phase-space distribution could also leave observable imprints on the matter power spectrum [9,50,51] if some of the χ particles could have non-negligible momenta during the structure formation. A thorough study on these effects would require solving the Boltzmann equations at the level of the phase-space distribution for both χ and φ which we leave for future work.

A Cross Sections
In this section, we present relevant cross-sections and decay widths that we have used in this paper based on the general effective Lagrangian in Eq. (2.3). For all the practical purposes, we shall treat the SM neutrinos as massless.
• χχ → νν where we have defined in order to simplify the expression.

B Energy Transfer
In this appendix, we present the explicit forms for the energy transfer terms in Eq. (3.16). Let us begin with the Boltzmann equation for the phase-space distribution of a particle species ψ: in which C[f ] is the collision operator which involves all possible processes for ψ, e.g. decay, inverse decay, elastic scattering, annihilation, etc, and thus depends on the phase-space distribution of all relevant species. For a specific process ψ + a + b + · · · ↔ i + j + . . . , the collision term takes the general form where dπ i = d 3 p i / (2π) 3 2E i , and "±" represents the Bose-enhancement/Pauli-blocking effects -one needs to choose "+" for bosons whereas "−" for fermions. In principle, the phase-space distributions f could take any form. For species in thermal equilibrium, f is either Bose-Einstein or Fermi-Dirac. However, for most of the practical purposes, the Maxwell-Boltzmann distribution can be used as a good approximation, and the condition f 1 is often satisfied. Therefore, one can use the approximations f ≈ e −E/T and 1±f ≈ 1 which greatly simplify the calculation.
To get the energy-transfer terms in the Boltzmann equation for the energy density, we just need to multiply Eq. (B.1) with E ψ and integrate over the phase-space element d 3 p ψ /(2π) 3 . The result is For the process ψ + a + b + · · · ↔ i + j + . . . , one can separate the contribution from the forward and backward directions, i.e. the contribution from the two terms in the square brackets in Eq. (B.2). Thus, the last term can be written formally as In the following subsections, we shall calculate the energy transfer terms for several particular processes.
• Decay: φ → νν The contribution from the decay process φ → νν is where Γ φ→νν is the decay width in the center of mass frame of φ.
• Inverse Decay: νν → φ The energy-transfer term of this process is n eq ν (T ) 2 P νν→φ = Here we use the superscript "eq" to emphasize that neutrinos are in thermal equilibrium with the rest of the SM thermal bath. Following the convention in in Gondolo's paper: Notice that this would be exactly the same with n φ (T φ )P φ→νν in Eq. (B.5) if φ is in thermal equilibrium with the SM thermal bath.
• 2-to-2 process Let us consider the 2-to-2 process 3 + 4 ↔ 1 + 2. The energy-transfer collision term associated with this process for particle 1 is where the "bar" over the amplitude |M| 2 indicates necessary average over different spin states. For our purpose, we assume that particles 3 and 4 are the states in the SM thermal bath, whereas particles 1 and 2 are dark-sector states. Since for freeze-in process, the energy transferred from the dark sector to the visible sector is negligible, we just need to calculate the energy transferred to the dark sector through the above process. Therefore, we just need to calculate the forward process -the part associated with the first term in the brackets: .

(B.16)
Using detailed balance f eq 3 (p 3 )f eq 4 (p 4 ) = f eq 1 (p 1 )f eq 2 (p 2 ) and unitarity, and inserting the flux factor F ≡ [(p 1 · p 2 ) 2 − m 2 1 m 2 2 ] 1/2 [52], we can rewrite the above integral as in which we have used the fact that v Møl = F/E 1 E 2 , and the second line is simply σ 12→34 . Comparing with the second term in the square brackets of Eq. (B.15), it is easy to notice that the calculation above is exactly equivalent to calculating the energy transfer rate of the inverse process 1 + 2 → 3 + 4 if particle 1 and 2 are in thermal equilibrium with the visible sector. This observation proves that n eq 3 (T )n eq 4 (T )P 34→12 = n eq 1 (T )n eq 2 (T )P 12→34 , (B.18) which allows us to trade the energy-transfer term of one process for that of its inverse. Following the procedure in Ref. [42], the Eq. where s 0 is the minimum of s.