Sterile Neutrino Dark Matter Production in the Neutrino-phillic Two Higgs Doublet Model

Sterile Neutrinos with a mass in the keV range form a good candidate for dark matter. They are naturally produced from neutrino oscillations via their mixing with the active neutrinos. However the production via non-resonant neutrino oscillations has recently been ruled out. The alternative production via Higgs decay is negligibly small compared to neutrino oscillations. We show that in the neutrino-phillic two Higgs doublet model, the contribution from Higgs decay can dominate over the contribution from neutrino oscillations and evade all constraints. We also study the free-streaming horizon and find that a sterile neutrino mass in the range of 4 to 53 keV leads to warm dark matter.


Introduction
The Standard Model (SM) of particle physics is very successful, but it fails to explain neutrino mass and dark matter (DM). Dark matter accounts for about one quarter of the energy density of the Universe, five times more than ordinary matter, but its origin is unknown. A good candidate for the dark matter are sterile neutrinos with a keV-scale mass and tiny mixing with the active neutrinos, which is a simple extension of the SM [1,2]. In contrast to standard cold dark matter, they are generally warmer with a larger freestreaming horizon. Thus they are candidates for warm dark matter and suppress structure at small scales addressing the missing satellite problem [3][4][5] and possibly explaining the velocities of pulsars [6,7].
There are many ways of producing sterile neutrinos: (i) they can be produced through neutrinos oscillations in the early Universe via a small mixing with the active neutrinos [8][9][10]. This mechanism is already excluded by observation [11], but the bounds can be avoided, if there is a large enough primordial lepton asymmetry and sterile neutrinos are produced via resonant oscillations [12]. (ii) Another well-studied alternative is non-thermal production via decay of a scalar field in thermal equilibrium [6,[13][14][15][16][17], or a scalar produced via the freeze-in mechanism [18], which subsequently decays to sterile neutrinos [17,[19][20][21].
Recently several alternative production mechanisms from decay have been suggested such as the production from the decay of pions [22], Dirac fermions [23], light vector bosons [24], or a condensate formed during inflation [25,26]. (iii) Finally the keV sterile neutrinos could have been in thermal equilibrium and their abundance diluted by production of entropy [27][28][29][30][31].
In any model with mixing between active neutrinos and the keV sterile neutrinos a fraction of the sterile neutrino abundance will be generated via neutrino oscillations. The mixing is generally induced after an electroweak doublet scalar obtains a vacuum expectation value (vev). The Yukawa interaction however induces a second contribution: The electroweak doublet scalar can decay into a SM lepton and a sterile neutrino. In a model with one Higgs doublet, this contribution is always negligible compared to the contribution from neutrino oscillations [32], because the vev, v = 174 GeV, and thus the mixing is sizeable. However this does not hold anymore in models with multiple Higgs doublets. The vev of one of the Higgs doublets might be tiny, smaller than O(MeV), and consequently the production via Higgs decay might dominate. Sterile neutrino dark matter with a keV-scale mass has been previously considered in a two Higgs doublet model in Ref. [33] and the production of the required number density via the decay of an electroweak doublet has been studied in Ref. [34] in the scotogenic model of neutrino mass [35].
We consider a two Higgs doublet model, where one of the electroweak doublet scalars exclusively couples to the sterile neutrino and study the production of keV sterile neutrino DM via the decay of this electroweak doublet in detail. The main result is the momentum distribution function for the sterile neutrino DM and the free-streaming horizon, which we use to determine the relevant parameter space where the keV sterile neutrino constitutes warm dark matter. This mechanism can be easily embedded in a seesaw [33,36] or radiative [34,35] neutrino mass model. The paper is organised as follows: In section 2, we introduce the two-Higgs doublet model and discuss the mass spectrum. The produced sterile neutrino DM abundance is discussed in section 3. Section 4 is dedicated to the free-streaming horizon of the sterile neutrinos and we briefly comment on the effective number of relativistic degrees of freedom in section 5. In section 6, we discuss the possibility that the scalar doublet obtains a tiny vev, thus the sterile neutrinos can decay and explain the observed an X-ray line at 3.55 keV [37,38]. We conclude in section 7. Technical details are collected in the appendices.

Neutrino-phillic Two Higgs Doublet Model
We consider a two Higgs doublet model with a second scalar doublet H ν which exclusively couples to the sterile neutrino N and the left-handed lepton doublet L. This is guaranteed by a Z 2 symmetry under which the new fields, N and H ν , are odd, but all SM particles are even. A coupling to other fermions is strongly constrained by flavour-violating processes. The most general Yukawa interactions in the lepton sector are given by and the most general scalar potential is defined in the usual way where we used a field redefinition of H ν to absorb the complex phase of λ 5 . The Yukawa couplings y LN,α are generally complex. After the SM Higgs doublet obtains a vev, v 2 = m 2 1 /λ 1 , electroweak symmetry is broken and we decompose the fields in terms of their components The scalar masses at leading order are given by where h describes the observed Higgs boson at m h = 125 GeV [39,40]. As long as H ν does not obtain a vev, there is no mixing between the different states. We will comment on this possibility in section 6. The active neutrinos obtain mass in the usual way. Both Dirac and Majorana mass terms are possible. For example, one can introduce the Weinberg operator [41], LLHH, to generate the neutrino mass. We will not discuss it further, because it does not affect the production of the keV sterile neutrino N .

Dark Matter Production
Freeze-in production of sterile neutrino DM with a second Higgs doublet has been studied in the scotogenic model [34]. Here we focus on keV sterile neutrinos and do not specify the mechanism of neutrino mass generation explicitly. The keV sterile neutrino can be produced via the decay of the scalar fields (k, K 0 , K ± ) while they are in thermal equilibrium. We will present a crude simple calculation using the Maxwell-Boltzmann approximation and neglect Pauli-blocking. The result will result in a good order of magnitude estimate. In our discussion, however, we will use the more accurate result using the distribution function, which can be found in App. C. Assuming that inverse decays can be neglected, the yield Y (T ) = n(T )/s(T ) of the sterile neutrino can be obtained from the Boltzmann equation, where s is the entropy density of the Universe, H(T ) is the Hubble parameter at a given temperature and γ N 1 (T ) is the thermally averaged sterile neutrino production rate, where X = k, K 0 , K ± and l is a SM lepton. Following the freeze-in calculation in Ref. [18], we can integrate the equation and obtain the final yield of sterile neutrinos after freeze-in where we defined the typical mass scale This allows us to rewrite the Friedmann equation during the radiation dominated epoch as using the Planck mass M P l and the usual definition of the effective entropy degrees of freedom g s * and the relativistic degrees of freedom g ρ * We can obtain an approximate analytic solution to this equation by extending the integration over all positive values of x, i.e. x min → 0 and .

(3.8)
In order to obtain simple analytic results, we made several (crude) approximations: (i) We used the Maxwell-Boltzmann approximation and thus also neglected Pauli-blocking of the neutrinos. (ii) We extended the integration boundaries in (3.3) to obtain the leading order result in Eq. (3.7). This is justified by noting that freeze-in is typically dominated by processes around T ∼ m X [18]. (iii) We assumed that the effective number of relativistic degrees of freedom for entropy, g s * , and energy, g ρ * , do not change during the production of dark matter, which is reasonably well satisfied for scalar masses m X 100 GeV. (iv) We neglected 2 ↔ 2 scattering processes like X 0 + ± → N + W ± and X ± + ν → N + W ± . These processes are subdominant compared to two body decays of X due to phase space suppression. (v) We assumed that the particles X are in thermal equilibrium until freeze-in occurs at x f i ∼ 2 − 5 [18] and thus T f i 20 GeV for m X 100 GeV. For sufficiently large Higgs portal couplings 1 λ 3,5 , the scattering of the scalars X with b quarks, X + b → X + b, will keep the scalars X in thermal equilibrium similar to the SM Higgs down to T ∼ 5 GeV. (vi) Finally we use the usual vacuum decay rate and do not take into account finite temperature effects, which has been studied e.g. in Ref. [42]. We expect these corrections  3 × 10 The contour labels are the effective coupling α |y LN,α | 2 . Figure 1: Contour plots with fixed DM relic abundance Ω N 1 h 2 = 0.1199 [43].
to be small, because the Yukawa coupling is small and the sterile neutrino a gauge singlet. In the following we will, however, use the more accurate result in Eq. (C.1). .
It has been obtained by solving the Boltzmann equation for the distribution function without using the Maxwell-Boltzmann approximation and including Pauli-blocking. The detailed calculation is outlined in App. B and C. The number of degree of freedom of the scalar fields are given by g k = 1, g K 0 = 1, g K + = 2 and the decay widths Γ X are given in App. A. We finally obtain the relic abundance of the sterile neutrino dark matter using Ω N 1 = m N s 0 Y ∞ N 1 /ρ cr with the critical energy density ρ cr = 3H 2 M 2 P l /8π and find . (3.10) Taking the limit of equal scalar masses, m k m K 0 m K + ≡ m kk and T d,X ≡ T d , the relic abundance only depends on three parameters: the two masses m kk , m N and the effective coupling α |y LN,α | 2 . Fixing the dark matter relic abundance to the observed best fit value for dark matter, Ω DM h 2 = 0.1199, by Planck [43], we show in Fig. 1 two contour plots illustrating the dependence on the three parameters. We find that the effective coupling has to be of order 10 −9 for keV sterile neutrino masses in the range 2 − 100 keV and scalar doublets with electroweak-scale masses. Although there is no explanation for the smallness of the Yukawa coupling, it is technically natural. Larger scalar masses generally require smaller couplings to compensate for the suppression by the scalar mass m kk . Larger sterile neutrino masses generally require larger effective couplings or larger scalar masses, because the DM abundance is proportional to the ratio m N /m kk .

Free Streaming Horizon
The free streaming horizon characterises the scale below which perturbations are suppressed in the power spectrum [44]. It is defined by the average distance a particle travels without any collisions where v is the average velocity, t in denotes the time when the sterile neutrino is produced and t 0 the time today. In order to evaluate it, we have to find the average velocity of the sterile neutrino where t nr (T N,nr ) denotes the time (temperature of the sterile neutrinos) when neutrinos become non-relativistic, i.e.
We can define the average momentum in terms of the momentum distribution function 2 f (p, t) of the sterile neutrinos, which can be obtained from the Boltzmann equation In analogy to the treatment in Refs. [15,17], we introduce the dimensionless quantity which allows to write the distribution function as (4.8) The scalar particle mass (decay width) is denoted m X (Γ X ) and we assumed that the effective number of relativistic degrees of freedom, g ρ * , remains constant during production, T d,X is the decay temperature of particle X and we neglected the back reaction from inverse decays. See appendix App. B for a derivation of the distribution function. As the integrand is exponentially suppressed for y 1, we can take the upper limit of the integral to infinity for temperatures T N m X / √ 8x N . This is generally justified at late times, (4.9) In this approximation the function x 2 f N (x) has only one maximum atx 1.54 and it falls off very quickly away from the maximum. A typical value for the momentum is thus of order xT N . This justifies the approximation for small temperatures. The average momentum is given by and thus lower compared to a sterile neutrino in thermal equilibrium with p T = 3.15T N . We find that the sterile neutrinos become non-relativistic at a temperature The corresponding temperature of the SM thermal bath T nr is obtained using the usual entropy dilution given in Eq. (C.15) with g s * (T d ) = g ρ * (T d ) = 110.75 and g s * (T nr ) = 3.94 and thus the time t nr , when the sterile neutrinos become non-relativistic, is determined by (4.12) Thus we can finally evaluate the integral for the free-streaming horizon by evaluating it piece-wise in the different regions set by t nr and t eq , respectively r FS = √ t eq t nr a eq 5 + ln t eq t nr 0.047 Mpc 10 keV m N . (4.13) where t eq = 1.9 × 10 11 s and a eq = 8.3 × 10 −5 . A detailed derivation can be found in Refs. [19,20,46]. Note that we do not have the additional entropy dilution factor, because it is already included in t nr in Eq. (4.12). Note, that the free-streaming horizon does not (strongly) depend on the mass of the heavy scalar or the effective coupling, which determine the sterile neutrino abundance. There is only an implicit dependence via the effective degrees of freedom at decay. A different scalar mass will lead to a different decay temperature, but for scalar masses above 100 GeV, the number of effective degrees of freedom stays almost constant. We show the free-streaming horizon in Fig. 2. The red region indicates when the free-streaming horizon becomes larger than 0.1 Mpc and the keV sterile neutrinos can be considered as hot dark matter. The blue region indicates the region when keV sterile neutrinos are in the cold dark matter regime. We take as benchmark value r F S 0.01 Mpc following Refs. [19,20]. We find free-streaming horizons within the desirable range for warm dark matter for a sizeable fraction of the parameter space. More precisely, the free-streaming horizon is in the warm dark matter range for keV sterile neutrino masses of 4 to about 53 keV. For example for m N = 7.1 keV, the keV sterile neutrino mass fitting the recently claimed X-ray line observation at 3.55 keV [37,38], we obtain r F S 0.06 Mpc, which is well within the range of warm dark matter.

Effective Degrees of Freedom
With the introduction of an additional sterile neutrino, one might wonder about its contribution to the effective relativistic degrees of freedom. Its contribution can be quantified by [17] where we subtracted the non-relativistic energy density contained in the mass of the sterile neutrino. The energy density ρ N is given in Eq. (C.8). Thus we find that there are no additional relativistic degrees of freedom in the non-relativistic limit r N ≡ m N /T N x, but in the ultra-relativistic limit (r N x) the additional effective relativistic degrees of freedom are Hence we find for a temperature T ν = T BBN 4 MeV shortly before the onset of big bang nucleosynthesis ∆N eff (T BBN ) ∼ 10 −3 with g s * (T BBN ) = 10.75, g s * (T d ) = 110.75, y LN,α ∼ 10 −8 , the scalar mass m kk ∼ 500 GeV, m N ∼ 10 keV and hence a negligibly small contribution to the effective number of neutrinos. At recombination the temperature is well below the mass of the keV sterile neutrinos, m N , i.e. r N x, and thus ∆N eff (T rec ) 0.

Sterile Neutrino Mixing and the X-ray Line
The scalar doublet scalar H ν may obtain a vev, H ν = 0 v ν T , similar to the SM Higgs. A vev will induce mixing between the sterile neutrino and the active neutrinos and lead to a new contribution to the active neutrino mass via the seesaw mechanism ∆m ν,αβ = m N sin θ α sin θ β . (6. 2) The mixing allows the production of the sterile neutrinos via neutrino oscillations [8][9][10], which can be estimated using the approximate formula [47] Ω N,osc h 2 0.2 × α sin 2 θ α 3 × 10 −9 m N 3keV In addition the sterile neutrino can decay into a photon and an active neutrino generating an X-ray line. This already constrains the mixing angle. See Ref. [11] for different constraints on a keV sterile neutrino which is produced via oscillations. Last year two independent groups observed an X-ray line at 3.55 keV [37,38], which can be explained by the decay of sterile neutrino DM to a photon and a neutrino with an active-sterile mixing, α sin 2 (2θ α ) 7 × 10 −11 . However, the observation is still debated [48][49][50][51][52].
The vev v ν can be naturally small and satisfy the X-ray bound, if it is induced via a small, possibly complex, Z 2 soft-breaking term [33,36], after the electroweak symmetry is broken. We obtain for the vev v ν .

(6.5)
There are no charge-breaking minima at leading order in µ 2 12 . The scalar masses will receive small corrections proportional to µ 2 12 appearing at the second order in the scalar mass in Eq. (2.4).
As the production via non-resonant neutrino oscillations can only give a subdominant contribution to the abundance of keV sterile neutrinos [11] 3 , the production is dominated by scalar decay. Thus we can translate the limit on the mixing angle from X-ray observations m kk [GeV] into a limit on the vev, v ν . We find v ν 1−10 MeV depending on the sterile neutrino mass m N 2 − 100 keV and the scalar mass m kk 100 − 1000 GeV. Hence the contribution to the active neutrino masses in Eq. (6.2) is much smaller than the solar mass scale. The small splitting might explain the existence of pseudo-Dirac neutrinos, if the active neutrino mass originates from a Dirac mass term. Taking the claimed hint for an X-ray line at 3.55 keV [37,38], we show in Fig. 3 the dark matter abundance as a function of the scalar mass m kk and the effective Yukawa coupling α |y LN,α | 2 fixing m N = 7.1 keV and the active-sterile mixing α sin 2 (2θ α ) 7 × 10 −11 . The red band indicates the 2σ-allowed region around the best-fit value measured by Planck [43]. The required effective Yukawa couplings are of order (5 − 22) × 10 −9 for scalar masses m kk 100 − 1000 GeV.

Conclusions
In the Standard Model of particle physics extended by one singlet Majorana fermion, keV sterile neutrinos are usually produced via neutrino oscillations and the decay of the electroweak Higgs doublet scalar into the keV sterile neutrino is a subdominant contribution. However, we showed in this paper, that the decay of an electroweak doublet scalar can become the dominant production mechanism in a two Higgs doublet model, if the vev of the electroweak scalar doublet is small or even vanishes. Neutrino oscillations will only account for a minor additional contribution to the keV sterile neutrino abundance. Thus we do not consider a vev of the electroweak scalar doublet for the first part of the paper and only indicate in Sec. 6 the possible changes, when the electroweak doublet scalar obtains a vev. As long as the vev is small enough, the produced abundance can be approximated by the result for a vanishing vev. The vev will lead to a mixing of the sterile neutrino with the active neutrinos. This renders it unstable via a two body decay allowing to search for the keV sterile neutrinos using X-ray line searches.
We explicitly derived an analytic expression for the momentum distribution of the keV sterile neutrino for late times, studied its free-streaming horizon and briefly commented on its contribution to N eff , the effective number of neutrinos in the early Universe, which is neglibly small. This production mechanism leads to a cooler spectrum of the sterile neutrino. The range for the sterile neutrino to be the warm dark matter is in between 4 and 53 keV.
The mechanism to produce sterile neutrinos via the decay of a Higgs doublet can be easily embedded in models of neutrino mass, which naturally explain the smallness of neutrino mass. This requires the addition of two or more sterile neutrinos. If the second Higgs doublet obtains a tiny vev, neutrino mass is naturally suppressed, while the Yukawa couplings of the other neutrinos can be relatively sizeable [33,36]. For a vanishing vev of the second Higgs doublet, neutrino masses can be induced via the radiative seesaw in the scotogenic model [35] and dark matter is produced via freeze-in [34].
Note added: A recent publication [55] pointed out additional thermal corrections to the production rate, which become relevant at high temperatures. We do not expect any significant corrections, because keV sterile neutrinos are dominantly produced at low temperatures, T m X . The inclusion of these effects goes beyond the scope of this work, which focused on a simple analytic discussion of this novel keV sterile neutrino production mechanism. A future numerical discussion should include these effects and also include the other neglected subdominant corrections discussed in Sec. 3.

A Decay Widths
The decay width of the scalar fields are given as Hence the decay widths in the limit m N m X can be approximated by

B Dark Matter Momentum Distribution Function
In case of the production of sterile neutrinos from decay of the real scalar k, the collision term is given by with the Lorentz-invariant phase space element which is denoted dΠ ν (dΠ k ) for neutrinos (the real scalar). The number of degrees of freedom are g k = 1 for the scalar and g ν = 2 for the neutrino. The distribution functions of the scalar (neutrino) are f k (f ν ) and the spin-averaged matrix element is given by Using the decay width Γ k |y LN | 2 m k /32π in the limit of negligible final state masses, m N = m ν = 0, we can rewrite the matrix element and generalise it to any of the three considered scalars using Eq. (A.4) Neglecting all terms proportional to f N and taking the ultra-relativistic approximation E N p N , we obtain Note that we do not neglect Pauli-blocking, i.e. do not approximate 1 − f ν 1. The integration over the neutrino momenta p ν is easily evaluated using the δ-function of the momenta. As the integrand only depends on E ν and thus the scalar product p X · p N , which can be evaluated via the delta function of the energies the collision term can we written as Using the dimensionless variables we rewrite the Liouville operator in the radiation dominated epoch assuming g ρ * to be constant L[f N ] = Hr ∂f N ∂r (B.10) and the collision term Hence the distribution function can be obtained by a simple integration If the scalar and the neutrino are in thermal equilibrium, they are described by a Bose-Einstein distribution function f −1 BE (y) = e y − 1 and Fermi-Dirac distribution f −1 F D (y) = e y + 1, respectively. We can perform the integral over y X , Changing the variables of integration to y = r 2 /8x N we obtain 2 e y dy (B.16) We will employ the integral representation of the incomplete polylogarithm 4 The integral representation of the incomplete polylogarithm is also known as the incomplete Bose-Einstein and Fermi-Dirac integrals, more precisely 17) which are well-defined for Re(s) > 0, to rewrite the integral in terms of known functions. The usual polylogarithm is obtained for b = 0, i.e. Li s (z) ≡ Li s (0, z). The first summand of Eq. (B.15) can be integrated directly using the integral representation of the polylogarithm in Eq. (B.18), while the second term requires integration by parts to obtain fractions, which can then be integrated using Eq. (B.18). Finally the distribution function can be expressed at late times with r 2 8x N by A comparison to the result without Pauli-blocking, shows that the Pauli-blocking leads to the correction As expected the correction is negative and becomes larger for smaller x N . This agrees with the distribution function given in Ref. [56]. In the limit x N 0 we obtain the Maxwell-Boltzmann result which agrees with the result obtained from taking the Maxwell-Boltzmann limit in Eq. (B.13).
where we used the incomplete Γ function The final expression for the distribution function at a finite temperature is given by We are mainly interested in the distribution function after the sterile neutrino DM production has finished and thus we will use the zero-temperature limit f 0 N (x N ) in the following.
The generalisation to multiple scalar fields is simply the sum of the different contributions The function x 2 f 0 N (x) has a maximum atx 1.54 and falls off very quickly away from the maximum. A typical value for the momentum is thus of orderx. After the sterile neutrinos are produced via the decay, they are completely decoupled from the thermal bath of the SM particles and generally have a different temperature T N compared to the thermal bath of SM particles. Obviously the distribution function of the sterile neutrinos does not resemble a thermal Maxwell-Boltzmann, Fermi-Dirac or Bose-Einstein statistic, but represents a non-equilibrium distribution function.

C Temperature of the Sterile Neutrino Dark Matter
In the following we will calculate all relevant quantities for the sterile neutrino dark matter sector for sufficiently late time and low temperature T N m X /2x N using the distribution function in Eq. (B.27). In order to simplify the notation we will also drop the subscript N from p N and x N = p N /T N . The number density and the average momentum are given by where we used the integral formulas ∞ 0 x n Li s (e −x/a )dx = a n+1 Γ(n + 1)ζ(n + s + 1) to evaluate the integrals. The result for the number density in Eq. (C.1) is consistent with the calculation of the yield Y = n/s in Eq. (3.7) using the integrated Boltzmann equation in the Maxwell-Boltzmann approximation. Using the expression for the number density and the distribution function we can determine how the temperature of the sterile neutrinos depends on the scale factor [57]. The sterile neutrinos are not interacting and thus the only effect is due to the cosmic expansion. After the abundance of the sterile neutrino is frozen-in, n(x)a 3 does not change and we can relate the number density n at time t to the number density n at a later time t , n(x)a 3 = n (x )(a ) 3 . Then we obtain from the number density of the sterile neutrinos with momenta between p and p + dp the following relation between the distribution functions at different times where we used the dependence of the momentum on the scale factor, p ∝ a −1 , to relate x and Thus the distribution function f N at a later time t has the same form with a temperature (C.7) The expression for the energy density is similarly obtained from with r N ≡ m N /T N , where we took the ultra-relativistic or non-relativistic limit, respectively, in the last step. The pressure is given by 61ζ(6)−64ζ (7) 29ζ (5) showing that the sterile neutrinos behave like matter for temperatures much below m N , i.e. r N 1. Similarly the entropy density is given by 29ζ(5)−32ζ (6) 26ζ(4)−32ζ(5) n 3.28n r N x (1 + w) r N n r N x . (C.11) Thus we find that the entropy density scales like T 3 N in the ultra-relativistic limit for r N x and observe that the entropy of the sterile neutrino dark matter sector is separately conserved using Eq. (C.7).
We can use the conservation of entropy to obtain a relation between the temperature of the SM thermal bath T and the temperature T N of the sterile neutrino dark matter sector. Consequently we parameterize the entropy density similarly to the SM thermal bath s = 2π 2 45 g s * T 3 and s N = 2π 2 45 g s * ,N T 3 N , (C. 12) we find for the entropy degrees of freedom in the sterile neutrino sector, g s * ,N , g s * ,N = X 5π 3/2 56 5 g ρ * (T d,X ) g X Γ X M P l m 2 X . (C.13) Using the conservation of entropy g s * (T d )T 3 d a 3 d = g s * (T )T 3 + g s * ,N T 3 N a 3 (C. 14) and assuming that all scalars decay at the same time T d , we find for the temperature T N , In the non-relativistic limit, we can employ Eq. (C.7) to relate the temperature T N to the temperature T N,nr , T N = T N,nr a nr a (C. 16) and ultimately the temperature of the SM thermal bath T .