Sterile neutrino Dark Matter production from scalar decay in a thermal bath

We calculate the production rate of singlet fermions from the decay of neutral or charged scalar fields in a hot plasma. We find that there are considerable thermal corrections when the temperature of the plasma exceeds the mass of the decaying scalar. We give analytic expressions for the temperature-corrected production rates in the regime where the decay products are relativistic. We also study the regime of non-relativistic decay products numerically. Our results can be used to determine the abundance and momentum distribution of Dark Matter particles produced in scalar decays. The inclusion of thermal corrections helps to improve predictions for the free streaming of the Dark Matter particles, which is crucial to test the compatibility of a given model with cosmic structure formation. With some modifications, our results may be generalised to the production of other Dark Matter candidates in scalar decays.


Motivation
The motion of celestial bodies over a vast range of scales, ranging from sub-galactic to supercluster scales, deviates from the prediction of general relativity and the observed distribution of baryonic matter. In addition, it seems impossible to explain the formation of the observed structures in the universe from the gravitational collapse of the primordial density perturbations observed in the CMB. All attempts to consistently explain these phenomena on all scales by a modification of the laws of gravity have failed so far. All known data are, however, in excellent agreement with the hypothesis that significant amounts of non-luminous Dark Matter (DM) are present in the observable universe. This hypothesis is not only consistent with the observed movement of stars, galaxies and clusters as well as structure formation in the early universe (with the known laws of gravity), it also gives an explanation for the apparent gravitational lensing of light from distant galaxies and quasars. If one takes this viewpoint, then these observations can be used to map the distribution of DM in the universe. The simplicity of this picture, the excellent agreement with data from very different physical environments and the lack of any known alternative explanation make the DM hypothesis very compelling. It is, however, not clear at this stage what the DM is made of.
The particles that compose the DM must be massive, electrically neutral, long lived and almost collisionless. The only candidates that fulfil these requirements within the SM are neutrinos. The known neutrinos, however, have masses below 1 eV and were relativistic at the time when structures in the early universe formed. The free streaming of such Hot Dark Matter (HDM) would suppress the formation of small scale structures in the universe, which is in contradiction with observations, see e.g. [1] for a summary. Hence, DM must be made of one or more new particles. One obvious candidate would be heavier neutrino mass eigenstates, which could be consistent with structure formation constraints. A comprehensive overview of particle physics scenarios that realise this idea and constraints on them from experiments, astrophysics and cosmology are given in [2]. In the following we recap only the basic facts that are relevant for the present work.

Heavy sterile neutrinos
Speculations about the existence of heavy neutrinos N have been motivated for various other reasons, see [3] for a review. Adding new leptons with gauge charges to the SM requires to add a whole fourth generation of particles to ensure cancellations of all anomalies. This scenario is strongly constrained by experimental data, hence the heavy neutrinos should be "sterile", i.e., singlets under the SM gauge groups. In the strict sense, we use the term sterile neutrinos for singlet fermions that mix with the known SM neutrinos ν L , which are usually identified with right handed neutrinos ν R . However, many authors use the term in a more loose sense for any fermionic gauge singlet (or heavy neutral lepton), regardless of whether or not it mixes with ν L . Since the scalar decay production mechanism we discuss here does not rely on a mixing of N with ordinary neutrinos, our results apply in both cases.
Seesaw type heavy neutrinos -The probably strongest motivation for sterile neutrinos in the strict sense, i.e. right handed neutrinos ν R that do mix with ν L , are the observed neutrino oscillations, which have been awarded with the 2015 Nobel Prize in physics. In the SM neutrinos only exist with left handed (LH) chirality, while all other fermions come with LH and right handed (RH) chirality. An explicit mass term for LH neutrinos ν L is forbidden by gauge symmetry, and the Higgs mechanism cannot generate a Dirac mass termν L m D ν R for them without existence of a RH counterpart ν R . If ν R exist, they are allowed to have an explicit Majorana mass termν R M M ν c R because they are gauge singlets. Here ν c R = Cν R T . If ν R exists, the mixing between the known "active" (SU(2) charged) neutrinos ν Lα (with α = e, µ, τ ) and the "sterile" singlets ν RI (where I = 1 . . . n labels the families of RH neutrinos) generates a Majorana mass term for the ν Lα and explains the neutrino oscillations via the seesaw mechanism [4][5][6][7][8][9].
The seesaw mechanism predicts the existence of n new neutrino mass states N I that are almost identical with the singlets ν RI , but contain a small admixture θ αI ν c Lα of the doublet fields. One can describe them by Majorana spinors Here θ = m D M −1 M , V N is the equivalent of the Pontecorvo-Maki-Nakagawa-Sakata matrix in the sterile sector and U N is its unitary part. We can in good approximation take V N = U N = 1. 1 The active-sterile mixing θ leads to a θ-suppressed weak interaction of the N I , and this is the only way how they couple to the SM at low energies. There have been different suggestions what the role of the N I in cosmology and particle physics could be. Their CP-violating interactions could, for instance, generate the observed baryon asymmetry of the universe [11] via leptogenesis during their decay [12] or production [13,14] in the early universe. For very low eigenvalues M I of M M , they could be responsible for the LSND [15], Gallium [16] and reactor [17] neutrino oscillation anomalies and/or act as extra relativistic degrees of freedom in the early universe ("dark radiation"). Last but not least, the N I are obvious DM candidates: They inherit all advantages of the known neutrinos, but their masses are given by the eigenvalues of M M and can be much heavier. Which of these roles the N I are able to take strongly depends on the magnitude of the M I , see e.g. [3,18]. For masses below the electroweak scale, various experimental constraints exist [10,19,20], and the N I can be searched for in near future experiments [20][21][22][23][24][25][26][27][28][29][30].
Observational tests -The idea that heavy neutrinos compose the DM can be tested in various different ways, leading to significant constraints on their properties. Let us in the following focus on one singlet fermion N with mass m N that should compose the DM. The most model independent constraint on its properties can be found by applying phase space considerations [31] to DM dense regions, which imposes a lower bound on m N . The precise value of this bound depends on the N phase space distribution function f N , which depends on the way how N were produced in the early universe. It usually lies in the range of a few keV [32], 1 − 2 keV can be used as a conservative estimate.
If N is a sterile neutrino in the strict sense, we can identify N ≡ N 1 in (1). Then N at least interacts via its total mixing angle U 2 ≡ α |θ α1 | 2 . This makes decays N → ννν into neutrinos possible, and the N -lifetime ∝ U −2 m −5 N must be at least comparable to the age of the universe. Moreover, the radiative decay N → γν [33,34] predicts an emission line of energy m N /2 from DM dense regions in space [35]. Non-observation of this emission imposes an upper bound on U 2 for given m N . In 2014, two independent publications reported the detection of an unidentified emission signal at 3.5 keV that could be interpreted in this way, though this claim is disputed [36][37][38][39][40][41][42][43]. For relatively small masses, the mixing U 2 can also lead to a signal in an upgraded version of the KATRIN experiment [29,44]. A direct detection [45,46], on the other hand, would be very challenging [47].
The free streaming of N and its effect on structure formation provide another way to constrain the parameter space. These depend strongly on the N phase space distribution function f N , which is determined by the way the heavy neutrinos are produced in the early universe. The precise shape of f N can be found by solving a set of momentum dependent kinetic equations for the coupled system composed of the scalar, N and the thermal plasma, see e.g. [48][49][50][51][52][53]. A detailed discussion can be found in section 5 of [2]. Doing this requires knowledge of the N -production rate. In this work we calculate thermal corrections to the production rate due to interactions with the primordial plasma or a pre-existing N -population.
Sterile neutrino distribution function -The framework of nonequilibrium quantum field theory allows to systematically compute thermal corrections to the production rate of heavy neutrinos. Some elements of nonequilibrium quantum field theory are briefly summarised in the appendix A. The field theoretical equivalent to the classical phase space distribution f N q is given by the function f N (q 0 ) in the Kadanoff-Baym ansatz (111) for the heavy neutrino propagator when evaluated at the quasiparticle pole Ω N q , Due to the feeble interaction of the singlets N , we can for all practical purposes approximate Ω N q ω N q = (q 2 + m 2 N ) 1/2 . Our main results for the rate of N -production, which can be expressed in terms of one-dimensional integrals, are valid for arbitrary f N . That is important because f N is a dynamical quantity that changes throughout the production, and the production rate in each moment in time is affected by f N in that moment.
In addition to the general expressions, we provide illustrative analytic results by employing a simple parametrisation, 2 The ansatz (3) has a simple physical interpretation. b = 1 corresponds to "kinetic equilibrium", i.e., the distribution function has the same shape as a Fermi-Dirac distribution, but a different normalisation. For b = 1 the average momentum is shifted with respect to the background plasma temperature T : larger b correspond to a smaller effective Ntemperature. This could e.g. be realised if a N -population that was produced at earlier times has been diluted by entropy injection into the plasma after it decoupled. One can constrain a and b by the requirement that N with distribution f N at a given temperature T make up for a fraction r of the observed DM density Ω DM , where s = 2π 2 g s T 3 /45 is the radiation entropy density with the effective number g s of degrees of freedom, m P is the Planck mass and H the Hubble parameter. All quantities with a subscript 0 refer to present day values. Plugging in Ω DM = 0.268 [54] for b = r = 1 yields a 5 × 10 −8 GeV/m M , which implies a < 0.1 for any m N that is consistent with phase space analysis bounds.

Singlet fermion production
Thermal N production via active-sterile mixing -If N is a sterile neutrino in the strict sense (1), then a minimal amount of N is produced thermally via the mixing θ with active neutrinos [55] (Dodelson-Widrow mechanism [56]). The expected free streaming classifies this contribution as Warm Dark Matter (WDM). The efficiency of the thermal production can be resonantly enhanced by the Mikheyev-Smirnov-Wolfenstein (MSW) effect if there are considerable lepton asymmetries µ α in the plasma [57]. The MSW resonance does not only affect the total number of produced N -particles, but also their momentum distribution f N , which determines their free streaming length and thereby affects the formation of structures. The requirement to produce the right amount of DM from mixing alone imposes an upper and a lower bound on U 2 for given m N . The upper bound, which requires the N -density to remain far below its equilibrium value, turns out to be weaker than the bounds from emission line searches and structure formation. The lower bound is crucial, but due to the MSW effect it depends on the value of the µ α during N -production.
The requirement to explain the existence of small scale structures observed in Lyα absorption in quasar spectra puts an upper bound on the free streaming of DM. This already appears to exclude the scenario where all DM is composed of thermally produced N for µ α = 0 (i.e. in absence of a MSW resonance) [52,58]. However, for µ α = 0 the momentum distribution tends to be "colder" and can be consistent with structure formation [59][60][61][62] if m N > 3.3 keV [63]. It is difficult to work out all the details of the MSW-enhanced production because the N production happens to peak at temperatures at which quarks start to form hadrons. Moreover, the MSW resonance produces nonthermal spectra, for which it is difficult to simulate structure formation [61]. In spite of significant progress [64][65][66][67] considerable uncertainties remain. At present it seems likely that thermally produced sterile neutrinos can explain all the observed DM consistent with the formation of small scale structure only if their production is resonantly enhanced by lepton asymmetries µ α . The allowed range for m N reaches from ∼ 3.3 keV to a few tens of keV, depending on µ α . The mixing must be smaller than U 2 < 10 −8 − 10 −12 , depending on m N . 3 Note that N with such small mass and mixing can avoid constraints from the Cosmic Microwave Background or Big Bang Nucleosynthesis [19]. However, they give no significant contribution to the generation of active neutrino masses in the seesaw mechanism. Hence, one needs at least two heavier siblings N 2 and N 3 to explain the two observed active mass splittings in the seesaw model. Interestingly, this minimal scenario is able to simultaneously explain neutrino oscillations, the observed DM density and the baryon asymmetry of the universe [14,69,70]. Alternatively, there can of course be an additional source of active neutrino masses.
Production by gauge interactions -At low energies, N in (1) only interacts with other particles via its mixing U 2 with active neutrinos. This need not be true at high energies in the early universe, as N are charged under some spontaneously broken gauge group in many extensions of the SM. In left-right symmetric theories, for instance, the ν R (and hence N ) are charged under a RH SU(2) gauge group. These interactions bring them into thermal equilibrium in the early universe, which would lead to a too large DM density.
There are different ways to avoid this problem. The equilibration is avoided if the wouldbe freezeout temperature of the RH gauge bosons exceeds the (p)reheating temperature of the universe. This temperature is unknown; it may be constrained by CMB observations [71,72]. If the N equilibrate, their abundance can be made consistent with observation if the number g * of degrees of freedom in the primordial plasma changes significantly after N -freezeout (e.g. because of a phase transition). Another (related) possibility is that the N are diluted by a release of entropy into the primordial plasma after their freezeout [73][74][75][76].
N production in scalar decays -N -particles can also be produced non-thermally in the decay of a heavier particle. This is the scenario under consideration here. The heavy particle may be the inflaton [77,78], another scalar singlet [49,50,[79][80][81][82], or a charged scalar [83]. Note that the minimal coupling to the SM already includes some production from the decays of scalars, namely pions [84][85][86] and SM Higgs particles [87,88], which is, however, sub-dominant compared to the thermal production from mixing. Higgs decay may give a significant contribution in models with an additional leptophilic (or, more generally, DM-philic) Higgs doublet, which can decay into N at much higher rate [89][90][91]. The N production may happen while the scalar field is in equilibrium (which is usually the case if it is charged) or while it "freezes in", leading to somewhat different phenomenology [53]. Since the decay mechanism does not rely on active-sterile mixing to produce the N , it does not impose a lower bound on U 2 and can avoid any constraints from searches for emission lines. Moreover, it seems to lead to relatively cold DM spectra that are consistent with structure formation [52]. In the present article, we compute thermal corrections to the abundance and momentum distribution of N -particles.

Goal and outline of this work
In this paper we calculate thermal corrections to the production rate of sterile neutrinos N in the decay of massive scalar fields φ (gauge singlet) and Φ (charged). For these considerations, it is irrelevant whether or not N mixes with ordinary neutrinos; our results apply to any massive gauge singlet fermion that is produced in the decay of scalars. If N is a sterile neutrino in the sense of (1), then there is an additional population of N that is produced thermally via the Dodelson-Widrow mechanism. However, at the temperatures where scalar decays usually happen, this thermal production by mixing is negligible 4 and we can set U 2 = 0. Then (1) simplifies to N ν R + ν c R , and N approximately does not interact with the SM.
In the literature it is usually assumed that the production rate of N particles with a given momentum p can be obtained from vacuum decay rate of scalar particles into heavy neutrinos. The momenta of the heavy neutrinos are uniquely fixed in this 1 → 2 decay. However, in reality the decay happens in the hot and dense primordial plasma that filled the early universe, and it is well-known that the damping rate of scalars receives thermal corrections in this regime [92]. There are several reasons for this. First, quantum statistics can enhance or suppress the decay rate, depending on whether the final state particles are bosons or fermions. Second, the dispersion relations of quasiparticles in a dense plasma usually differ from those of particles in vacuum. This can change the phase space. Third, at high temperature the decay scalar → N + something is not the only process that can produce N ; inelastic scatterings and Landau damping can also contribute to the rate.
In this work we calculate thermal corrections to the DM production rate from first principles of nonequilibrium quantum field theory. This in principle is a very complicated calculation because the production of N involves the scalar field, which is potentially far from thermal equilibrium during production (if it happens during freeze-in). To avoid calculations with nonequilibrium propagators of the scalar in the loop, we use a trick and calculate the thermal damping rate for the scalar-quasiparticles with momentum q due to the interactions with N . This can be done without specifying the thermodynamic state of the scalar field. The total number of N -particles and their momentum distribution can be found by plugging the production rate into a set of Boltzmann equations for N and the scalars. Here we assume that the N -particles are produced in the decay of the scalarparticles. If they are produced in the decay of the condensate or "classical field" ϕ = φ (or Φ † Φ 1/2 ), the effective masses and couplings are generally ϕ-dependent. Thermal corrections to the dissipation rate in this case have e.g. been studied in [93] and references therein.
Compared to previous calculations, our method allows to include all quantum statistical factors and a proper treatment of the dispersion relations in the plasma. It also systematically includes processes other than the decay to leading log accuracy in the gauge coupling. For instance, if Φ is electrically charged and decays into N and a charged fermion Ψ as Φ → N Ψ, then there are also scattering processes Ψγ → ΦN (s-channel) and γΦ → N Ψ (t-channel) as well as their inverse that change the number of N and Φ particles, or the final state may radiate off a photon (Φ → N Ψγ). Although these are of higher order in some coupling constant, they can become important at high density because of the Bose enhancement that photons are subject to. The same happens if Φ is charged under a non-Abelian group and e.g. has electroweak interactions. This has not been taken into account in past calculations. In this work we calculate the leading order corrections to the production rate due to these effects.
The paper is organised as follows. In section 2 we calculate the thermally corrected rate of N -particle production from neutral scalar decays. In section 3 we calculate the thermally corrected rate of N -particle production from U(1)-charged scalar decays. The equivalent rate for SU(N) charged scalars can be obtained from that by a few simple replacements. We also compute the leading log contribution to the N -production from scatterings in the plasma. The rates can in general only be determined numerically. We give analytic approximations that hold for the regime where the decay products are relativistic. In section 4 we discuss our results and conclude. In the appendix A we briefly recall the methods from nonequilibrium quantum field theory required for this calculation, for more details see [94] (in general) and [95,96] (specifically for the approach used here).

Scalar singlet decay
We first treat the case of production from the decay of a real singlet scalar field φ. That is, we consider the Lagrangian Here h is the Higgs doublet andh = h * with being the SU(2) antisymmetric tensor. L is the SM lepton doublet, F is a matrix of Yukawa couplings, the potential V (φ) includes a mass term m 2 φ φ 2 /2 and y is the Yukawa coupling that mediates the φ decay into N . Recall that N = N c is a Majorana spinor. Since we are not interested in production via active sterile mixing, we immediately set F = 0. L φint contains interactions of φ with other fields than N . Such interactions must exist in order to produce φ in the early universe. L φint may e.g. include a Higgs portal termλφ 2 h † h, Yukawa interactions or other scalar interactions.

Heavy neutrino production rate
In nonequilibrium and thermal field theory, gain-and loss rates for the occupation numbers of weakly coupled fields are related to the discontinuities of self-energies evaluated at the quasiparticle pole [92], see appendix A. Using the finite temperature optical theorem and cutting rules [97], one can interpret these rates in terms of decays and scatterings involving the external particles as well as the "cut" propagators in the loop [98,99]. In principle, the most straightforward way to calculate the number and momentum distribution of Nparticles of a given momentum p that are produced in a given time interval is to compute the N -self energies Σ ≷ N in the closed time path formalism and obtain rates of the form (134). This, however, requires evaluating Feynman diagrams with nonequilibrium propagators for both, the decaying scalar φ and N , in the loop. The interactions of the singlet fields N are typically feeble, which justifies to use the Kadanoff-Baym ansatz (111) and apply the narrow width approximation (Γ N q = 0) for propagators inside the loops. For the scalar field it is not obvious that these approximations are always justified because it usually has other interactions in addition to the coupling to N . In the present work we are interested in scenarios where N -production is dominated by decays. In this regime we can use a trick to avoid the scalar propagator in the loops: Instead of directly calculating the production rate of N , we compute the damping rate of φ due to interactions with N and use it to calculate the number of produced N -particles via (14).
Kinetic equations for the occupation numbers -The occupation numbers for φmodes can be characterised by a function f φq (t), which follows the effective kinetic equation (137), Here Γ ≷ φq are gain-and loss rates and Γ φq = Γ > φq − Γ < φq can be identified with the total thermal damping rate. The precise definitions of these rates are given in the appendix by (126), (127) and (125). In the following we will use the symbolΓ φq for the contribution to the thermal φ-width Γ φq that comes from diagrams with N -propagators in the loop. We can formally defineΓ φq by the equation where Γ φq is the total thermal φ-width and Γ φq is generated from the interactions in V (φ) and L φint . Analogously, one can define self energiesΠ ≷ φ and ratesΓ ≷ φq from (126) and (127). For convenience, we can rewrite (6) as Note thatf φq implicitly depends on y and f N because it has to be evaluated with the full Γ ≷ φq . If the interactions in L φint are numerous or much stronger than the Yukawa coupling (e.g.λ y), then one may expand inΓ φq /Γ where f (0) φq ≡ f φq | y=0 is independent of y and f N . The total number of N -particles is given by the rate equation The factor 2 in front of the integral is due to the fact that two N -particles are produced in each φ-decay. If the main contribution toΓ φq comes from 1 → 2 decays and their inverse, one can expressΓ and Here Ω φq is the mass shell of a φ-quasiparticle with spatial momentum q, see (118), and Ω N p the energy for a N -quasiparticle with spatial momentum p. In the second equality we usedγ φ (p, q) =γ φ (q − p, q) for Ω φq = Ω N p + Ω N p−q . We should recall that the quantities Ω φq , f φq ,f φq ,γ φ (p, q) and f N p all depend on time. In the following, we focus on the computation of the quantitiesΓ > φq andΓ < φq andγ φq , i.e., the gain-and loss rates for N -production. From (9)- (14) it is clear that the full Γ ≷ φ is needed to determine the time evolution of f φq and calculate f N p in general.
In the special case that φ predominantly couples to fields that are in equilibrium, Γ ≷ φq fulfil a relation of the type (132) and one can further approximate where f B (Ω) = (e Ω/T − 1) −1 is the Bose-Einstein distribution and we have suppressed the time dependence of all quantities. No knowledge of the full Γ φq is needed to compute f B (Ω φq ). However, Ω φq depends on time, and the parameters in L φint still affect f N p because the time evolution of f φq is governed by (9).
The damping rate -As summarised in the appendix, Γ φq can be determined from (123) by calculating the imaginary part of the retarded self-energy for φ at the quasiparticle pole, Γ φq is obtained from the y-dependent part of the self-energy ImΠ R φ . The ratesΓ ≷ φq are obtained equivalently fromΠ ≷ φ . At one-loop level, the self energiesΠ ≷ φ can be calculated from the diagram by applying finite temperature Feynman rules [100]. Here the solid lines represent the external φ and the dotted lines N -propagators. Using the expressions given in appendix A, the self-energies read In combination, this gives If N had gauge interactions while being produced in the primordial plasma, then this would not only add another efficient N -production mechanism (thermal production via gauge interactions), but would also modify the production from scalar decays [92,101]. In this case the results from example IV in section 7 of [102] may be applied. We will not repeat this calculation here. Instead we assume that N is indeed a singlet not only under the SM gauge groups, but with respect to all interactions that are relevant in the early universe.

Small Yukawa couplings y: Production in decays
It is often assumed that the coupling y is very small, y < 10 −5 − 10 −9 [53]. Then the term yφN N can be used to generate a keV scale Majorana mass term from a GeV-TeV vev of φ, without having to add a small Majorana mass N m N N/2 in (5) by hand. In this case the effects of forward scatterings (which modify the N dispersion relation in the plasma) and inelastic scatterings (which increase the thermal width of N and thereby reduce its effective lifetime) are negligible even if there is already some N -population in the plasma. Then we can approximate ρ N (p) with the free spectral density, This essentially means that we use free thermal propagators for N (instead of resummed propagators), which greatly simplifies the evaluation of the integrals in (17)- (19). It also implies that the energy for N -quasiparticles in the plasma is in good approximation Corrections are suppressed by both, the smallness of y 1 and f N 1 due to (4).
Calculation of the self-energies -We first focus on the case q = 0. The results arẽ Let us now turn to the case q = 0. In a thermal bath, Γ φq in general cannot be obtained by multiplyingΓ φq=0 with the time dilatation factor m φ /(m 2 φ + q 2 ) 1/2 . This can be understood by recalling that 1/Γ φq is the lifetime of a quasiparticle and determines its mean free path in the bath. The difference in the lifetime of a q = 0 and q = 0 quasiparticle in a bath is not only due to time dilatation, but also due to the fact that the cross section for scatterings with bath particles (which reduces the lifetime) generally depends on the energy of the colliding particles. Moreover, the quantum statistical factor (1 − 2f N ) breaks Lorentz invariance: A moving φ-particle decays into N with different momenta than one at rest, and the occupation numbers in those modes generally differ. In the approximation (20) for very small y, the effect of scatterings is negligible, but quantum statistics may still play a role. To evaluate the integrals in (17)-(19) at q = 0, we neglect the mass m N for simplicity, which seems justified for m N ∼ keV. Using δ( where p = |p|, q = |q| and x = pq/pq is the cosine of the angle between the spatial vectors p and q. The remaining δ-function can be used to evaluate the x-integral, which fixes the limits for the p-integration: It is worth to note that no assumptions about f N were required to obtain these expressions. If f N can be approximated by (3), the final integrals can be solved analytically, e.g.
If there is no pre-existing N -population present at the onset of the φ-decay, then we can set f N = 0 and all temperature dependent effects in (31) seem to disappear. This is what one would intuitively expect, and this is why the vacuum decay rate has been used in the past literature.
Production rate -From (125) it is clear thatΠ − φ should be evaluated at the φ-quasiparticle pole Ω φq . The φ-quasiparticle dispersion relation generally differs from that is, Ω φq = ω φq . This is true even if φ is a singlet because φ can have Yukawa couplings or couplings to other scalars, and V (φ) in general contains self-interactions. In order to be produced in the early universe, φ necessarily must have some interactions, and these will necessarily shift the quasiparticle mass pole. If, for example then with the thermal mass Hence, even for f N = 0 there is a thermal correction toΓ φq . To estimate this effect, we compare (24) to the vacuum decay rate (m φ m N ) and express 5Γ This can be compared to the production rate of a Dirac fermion at temperatures below its mass in section 7.1 of [102]. The only difference is an overall factor 2 for the Dirac fermion. This difference arises from the fact that, for a Dirac fermion, there are two diagrams of the type considered here, with a fermion flow in the loop going in opposite directions. More physically, a Dirac fermion has twice more internal degrees of freedom into which φ can decay. If one would calculate corrections to (24) by using dressed spectral densities instead of (20), then one would generally see further differences between the production of Dirac and Majorana fermions because the interactions of Majorana fermions are generally different. However, for small y such corrections are negligible, and the free spectral density for Dirac and Majorana fermions is the same, hence the factor 2 is the only difference. It is worth to emphasise that the expression (39) holds for any sterile neutrino distribution function f N , i.e., it does not rely on the ansatz (3).
Based on (28)-(30), one can calculate rates for arbitrary q = 0 as Here we have neglected the temperature dependence of Z. 6 With the approximation (34), one findsΓ One can combine them to expressΓ φq in terms of the vacuum decay rate for particles at restΓ 0 ,Γ The physical interpretation of the different factors in this expression is very simple.Γ 0 is simply the decay rate for φ-particles at rest in vacuum. The factor M φ /m φ takes account of the fact that the effective mass of φ-quasiparticles in the plasma is larger than in vacuum, leading to a reduced lifetime. The factor M φ /Ω φq takes account of the fact that the lifetime of a particle that moves with momentum q is larger than that of a particle at rest due to time dilatation. In vacuum, the dilatation factor would be given by m φ /ω φq ; this cannot be applied in the thermal plasma (c.f. figure 5), but the obvious generalisation M φ /Ω φq works. What is left is a quantum statistical factor that breaks Lorentz invariance (because the thermal bath singles out a reference frame with respect to which q is measured). For Inserting (46) into (14), we obtain These expressions are our main results for the case of singlet decay. In the derivation of above results, we only made use of the relation (113), which holds for any Majorana fermion. Therefore they hold for any phase space distribution function of N -particles. In particular, the results do not rely on the validity of the ansatz (3). If the ansatz (3) provides a valid approximation to the phase space distribution of N -particles, then the final integral can be solved analytically, For a = b = 1, this can be compared to the results found in [96,[102][103][104]. Our result shows that, as long as the N -occupation numbers are far below their equilibrium value (a 1), thermal effects actually enhance the production rate at high temperature because M φ m φ for T m φ / √ λ. Once some amount of DM has been produced, a = 0 does not exactly hold. One may wonder if the statistical factors f N can have an effect towards the end of the DM production. Moreover, it is interesting to see what effect a pre-existing Nabundance may have. The considerations following (4) show that the effect of f N is always negligible for a standard thermal history and if the average N -momentum is similar to T (b 1). For b 1, there is room for a sizable a, and Pauli blocking may be significant for some modes. This either requires a production mechanism that lead to a rather "cold" pre-existing N -population or a significant injection of entropy into the primordial plasma after N -decoupling. Some results are plotted in Figs. 1-5. In addition to (19) there are other contributions toΓ φq at order y 2 , but these always contain additional powers of the coupling constants in V (φ) and L φint as well as loop factors. For instance, the potential (33) gives rise to diagrams of the form , cuts through which e.g. include processes φφ → φN N . The size of these can be estimated by considering the equilibrium situation, in which case they are y 2 -suppressed corrections to a "setting sun" contribution ∼ 10 −3 λ 2 T 2 /Ω q calculated in appendix B of [93]. Using (48) and (35), one can estimate that they are always suppressed by a factor λ relative to (48) unless there is a very dense pre-existing population of N .  Figure 3: The rateΓ φq in units of the time-dilated rateΓ 0 m φ /ω φq for q = m φ , with M φ given by (35) and λ = 1, as a function of T . The blue curve is for a = 0, the red curve for a = 1. The enhancement of the rate due to the increasing thermal φ-mass kicks in when M φ /m φ 1 above T m φ 24/λ 5m φ . For a = 1, it competes with a Pauli suppression, which reduces the rate for T > m φ . Comparison with Fig. 3 shows the effect of changing λ.

Sizable Yukawa coupling y and scatterings
The use of the free spectral density (20) is based on the assumption that the N -particles essentially do not feel the primordial plasma any more after they were produced. Since the yφN N -vertex necessarily involves another N , this assumption is not only justified by the small Yukawa coupling (y 1), but also by the low density of N -particles (f N 1). The use of (20) is only questionable if both of these suppressions are avoided, i.e., if y is relatively large and there is a pre-existing N -population. We do not treat this special case here, but it is clear that one expects more significant thermal corrections toΓ φq . We expect these to be very similar to the thermal corrections found in section 7 of [102]. In that work, the produced fermion is a Dirac particle with gauge interactions. However, the thermal corrections to the spectral density of a Majorana fermion from Yukawa interactions [105] have the same structure as the gauge corrections for a Dirac fermion [106]. Finally, we would like to recall that (14) should only be used to determine f N p if N -particles are predominantly produced in 1 → 2 decays and do not re-scatter within a Hubble time. For sizable couplings y this may not be the case. If scatterings are important, the momentum dependent N -production rate should directly be computed from an appropriate nonequilibrium generalisation of (134).  Figure 5: The rateΓ φq in units of the time-dilated rateΓ 0 m φ /ω φq with M φ given by (35) and λ = 1 as a function of q. The blue curve is for a = 0, the red curve for a = 1. In the upper plot T = m φ , in the lower plot T = 100m φ . At high temperatures, the increased thermal mass M φ leads to a larger rate than the vacuum estimateΓ 0 m φ /ω φq suggests. The Pauli blocking for a = 0 is inefficient for momenta q T because the produced particles' momenta are outside the Fermi sphere.

Charged scalar decay
Let us now consider the decay of a charged scalar Φ. To be specific, we use a model described by the Lagrangian 7 Dark Matter is produced by the Yukawa interactioñ The case where Φ is part of an additional Higgs doublet h ν and carries SU(2) charge [91] can be treated in the same way. In this case Ψ should be identified with the lepton doublet L , and there are several massive scalars that decay. The decay of each of them can be treated equivalently to the present calculation. The only differences are that L is a chiral field (leading to a factor 1/2 in the production rate) and the fact that the leptons in the final state carry SU(2) charge. The way how SU(2) interactions modify the spectral density ρ(p) is the same as for U(1), except for a numerical factor in the thermal fermion mass M f defined below [106,107].

Heavy neutrino production
The kinetic equations for the occupation numbers during charged scalar decays can be obtained from (6)-(15) by the replacements φ → Φ, y →ỹ and Ω N p−q → Ω Ψp−q , where Ω Ψp−q is the energy of a Ψ-quasiparticle with spatial momentum p − q. The overall factor 2 in front of the integral in the kinetic equations (12), (14) and (15) should be dropped because only one N -particle is produced in decays Φ → N Ψ. This yields whereγ Φ (p, q) is given from The ratesΓ Φq andΓ ≷ Φq for Φ are defined analogously to the singlet case. As in the case of the singlet scalar we will approximate Ω N p ω N p = (p 2 + m 2 N ) 1/2 assuming that the coupling of N to the plasma is very feeble.
The fact that there is a charged particle in the final state makes the decay Φ → N Ψ different from φ → N N considered in the previous section. First, the Ψ-quasiparticles in the final state tend to be thermally populated in the early universe (with a Fermi-Dirac distribution f F for the occupation numbers) due to their gauge interactions, and one expects significant Pauli blocking even if the N population is low (f N 1). Second, a charged particle feels the presence of the primordial plasma. This leads to "screening", and the properties of Ψ-quasiparticles in the plasma are not the same as those of particles in vacuum. Finally, the gauge interactions of Ψ open up new production channels for N . For example, there are s-channel and t-channel scatterings γΦ → N Ψ with intermediate Φ or Ψ and "photons" in the initial state as well as their inverse processes. When N -particles are produced in decays Φ → N Ψ, their momentum distribution can still be calculated using Boltzmann equation (52). If scatterings play a significant role, (52) cannot be used, and it is more convenient to directly calculate the N -self-energies Σ ≷ N (p) to obtain f N p . However, our present approach can still be used to calculate the total DM density from a rate equation (51).
The thermal corrections can be incorporated systematically by calculatingΠ − Φ (q) to a given order in perturbation theory. In the following we calculateΠ − Φ (q) only to leading order in the small couplingỹ, but we use resummed hard thermal loop (HTL) spectral densities for Ψ. This is necessary because the "naive" loop expansion is not a consistent expansion in α at high temperatures [108,109]. While the use of resummed propagators for the interacting Ψ is mandatory, we can (assuming that f N andỹ are not too large) continue to use the free spectral density for the singlet ρ N (p) given by (20). The results obtained in this way are accurate to leading order inỹ and leading log in the gauge coupling α.
Herepγ γ γ = p i γ i /|p|. The two functions are the sum of singular contributions ρ pole ± and a continuous part ρ cont ± . The poles given by the singular parts define the energies (or dispersion relations) Ω ± of quasiparticles with momentum p, The continuous part is given by Here x = p 0 /|p| and y = 1 2 M f /|p|. The thermal fermion mass reads [106] Here C is the quadratic Casimir of the representation of the gauge group. For an U(1) interaction as considered here C = 1. If Ψ carries some kind of hypercharge Y = 1, one has to replace α → Y α. For SU(N) interactions the shape of ρ Ψ is be exactly the same as for U(1), except for the numerical factor C in the definition of M f that depends on the representation under which Ψ transforms. The mass M f is sometimes referred to as the asymptotic mass, it differs from the plasma frequency ω f at |p| = 0 by a factor √ 2. The residues Z ± are The dispersion relations Ω + and −Ω − have to be found as the solutions to There are two solutions, corresponding to two types of quasiparticles, with dispersion relations Ω + and Ω − . The former can be interpreted as a screened one-particle state, the latter are collective excitations [110]. They are often referred to as "holes" or "plasminos". There is an analytic expression for the dispersion relations Ω ± in terms of the Lambert W -function [105] with s = −e −(y −2 +1) .
Self-energies -The rateΓ Φq at leading order inỹ is given by the diagram , where the solid lines represent the external Φ, the dotted line is a free N -propagator and the dashed line a resummed Ψ-propagator. The arrow indicates charge flow. The gain and loss terms read They combine intõ The evaluation of (62)-(64) is complicated by the fact that the angle between the spatial vectors p and q appears in the complicated function ρ Ψ (p − q). This can be avoided by making a shift in the spatial part p of the integration variable p by q and introducing k = p − q. Afterwards we perform the integration over the angular variable x = qk/(qk) between q and k with the help of δ(p 2 − m 2 N ) in (20). The requirement that the zero of the δ-function must be in the interval x ∈ [−1, 1] fixes the limits of the p 0 integral to where we have again used the notation q = |q| and k = |k|. The other angular integration is trivial because of rotational invariance. Then it is straightforward to obtaiñ and thereforeΠ For the q = 0 mode, the x-integration in (62)-(64) is trivial, and the δ-function can be used to perform the p 0 -integral. This yields whereω k = (m 2 N + k 2 ) 1/2 . Note that the above can be also obtained from more general expression (68) by taking limit q → 0, which gives ω + − ω − ≈ 2qk/ω k for the integral interval, and by setting p 0 =ω k in the integrand. In (71) the δ-functions in ρ pole ± in principle allow to perform one more integration in the pole part analytically, but due to the appearance of Ω ± it is in general not possible to find the zero analytically.

Production in decays: Analytic approximations
For M Φ M f + m N the decay products have momenta ∼ M Φ /2 that are large compared to the typical energy ∼ T of particles in the plasma. In this regime, one can approximate Physically this approximation of Z ± means that the decay into holes can be neglected, which is intuitive because collective excitations are an infrared phenomenon that is only relevant for soft momenta. The screened particles with hard momenta have a dispersion relation like free particles, but with the vacuum mass replaced by the momentum independent thermal mass M f . This is also intuitive because the dispersion relation must be simple if the momentum exceeds all other scales in the problem, hence one expects that it is given by a constant mass term due to forward scattering. Neglecting the continuous part of ρ Ψ means that we only consider the contribution toΓ Φq from 1 → 2 decays and their inverse, which is reasonable at low T .
Self-energies -With the approximations (72) the integrals in (69)-(71) can be solved (75) shows that Pauli blocking becomes important as soon as T exceeds the vacuum Φ-mass m Φ even if there are no N -particles (a = 0). The reason is of course the Pauli blocking of the Ψ-particles in the final state. If we set q 0 = Ω Φq = M Φ , we can also again observe the enhancement of the rate due to the thermal mass M Φ . Since Φ is charged, it does not only receive a mass correction from its self-interaction, but also from the gauge interaction, and we have to take The term α 2 T 2 /4 arises because Φ carries a U(1) charge. If we were dealing with SM leptons and a Higgs doublet, as considered in [91], the equivalent to M f would be and the Higgs would have the thermal mass where the different constant factors arise from the hypercharges [111], cf. discussion after (58). If T is near or below the temperature where spontaneous symmetry breaking occurs, there are additional contributions that are proportional to the temperature dependent expectation value of the Higgs condensate ∝ h †h = h † h (and ∝ h † ν h ν in the leptophilic two Higgs doublet model).
The expressions (73)-(75) can be significantly simplified if one neglects m N and M f , Note that the assumption of massless final states is only justified if M Φ M f + m N at all relevant temperatures, which in the large T limit is guaranteed by (77) and (58) only if λ is sizable and Ψ does not have additional interactions. As in the singlet case, it is not possible to find an analytic expression forΓ Φq with q = 0 without additional assumptions about the phase space distribution function of N -particles. One can, however, again expressΓ Φq in terms of a one-dimensional integral. In the approximation m N = M f = 0 of massless final states and (72), the only contributing kinematic δ-function is the one corresponding to the decay Φ → N Ψ and we find where the integration limits have again been fixed by the usual kinematic considerations.  Figure 6: The ratioΓ Φq /Γ 0 in the approximation (88) for q = 0, with b = 1 and M Φ given by (77) and λ = 0, as a function of T . The blue curve is for a = 0, the red curve for a = 1.
In the upper plot the gauge coupling is chosen as α = 10 −2 , in the lower plot α = 1/2. The enhancement of the rate due to the increasing thermal Φ-mass is similar to the scalar case and would be even more prominent for λ = 0. For a = 0 the Pauli blocking can at most reduceΓ Φq by a factor 1/2, for a > 0 a stronger suppression is possible.
Production rate -From these results, we finally obtaiñ As in the previous section, we have expressed the thermal damping rate in terms of the vacuum damping rateΓ 0 . Note that, compared to the singlet case, there is a symmetryfactor 2 difference in the latter, i.e., in this section we usẽ instead of (36). For q = 0 we find With the approximation Ω 2 Φq M 2 Φ + q 2 , this reads The physical interpretation of this expression is very simple and the same as discussed following (45). From (53) we havẽ This allows to derive an expression similar to (47) for the charged scalar case, In the limit M f = m N = 0 of massless final states, (95) simplifies to which can be compared to (47). Note that in contrast to the singlet case there is no factor 2 in front of the integral because only one N -particle is produced in each decay. These expressions hold for an arbitrary sterile neutrino phase space distribution. If one assumes that f N can be approximated by the ansatz (3), then even the final integral in (92) can be solved analytically for arbitrary q, i.e.
The different factors in the results (93) and (97) have a simple physical interpretation, which is exactly analogue to the discussion following (45). At low temperatures, they reproduce the T = 0 decay rate. Once the temperature exceeds m Φ , Pauli-blocking of the Ψ-particle in the final state suppressesΓ Φq . This effectively reducesΓ Φq (for q = 0, T M Φ and f N 1 roughly by a factor 1/2). If there is a significant population of N -particles (a 1), then Pauli blocking also applies to the N in the final state as soon as the temperature approaches ∼ bm Φ . At even higher temperatures, when M Φ /m Φ 1, the thermal Φ-mass increases the rate again due to the reduced lifetime of a heavier particle. For a potential of the form (33) this roughly happens at T > m Φ (24/λ) 1/2 . This behaviour is illustrated in Figs. 6-8. Up to now we have completely neglected the contributions from the holes Ω − as well as the term ρ cont ± . At T m Φ /α this is justified because M Φ M f , so the decay products' momenta p ∼ M Φ /2 are hard with respect to the plasma. The holes Ω − are relevant only in a small temperature interval: The suppression of Z − is only lifted when the temperature is high enough that M f ∼ M Φ , but low enough to keep M Φ > m N + ω f because otherwise the decay into holes is kinematically forbidden. Outside this small interval, we do not expect the holes Ω − to play an important role.
At this point one should compare the analytic expressions obtained in the singlet and charged case, i.e., (30) to (85), (39) to (88), (45) to (92) and (48) to (97). This comparison  (97) for q = m Φ , with M Φ given by (77) and λ = 0, as a function of T . The blue curve is for a = 0, the red curve for a = 1. In the upper plot the gauge coupling is chosen as α = 10 −2 , in the lower plot α = 1/2. As in Fig. 6, the enhancement of the rate due to the increasing thermal Φ-mass is similar to the scalar case and would be even more prominent for λ = 0. For a = 0 the Pauli blocking can at most reduceΓ Φq by a factor 1/2, for a > 0 a stronger suppression is possible.  (97) for T = 100m Φ , with M Φ given by (77) and λ = 0 as a function of q. The blue curve is for a = 0, the red curve for a = 1. In the upper plot the gauge coupling is α = 10 −2 , in the lower plot α = 1/2. As usual, at high temperatures, the increased thermal mass M Φ leads to a larger rate than the vacuum estimateΓ 0 m Φ /ω Φq suggests. The effect is larger for stronger gauge coupling and would be more prominent for λ = 0.
shows that, under the present assumptions, the production rates of sterile neutrinos in singlet and charged scalar decays look almost the same and are governed by the same physical effects, namely Pauli blocking and "thermal masses" in the plasma. The only differences are the occupation numbers (the charged particles in the final state of Φ-decays are in thermal equilibrium) and the actual values of the thermal masses (which depend on the interactions that are responsible for forward scatterings and screening). There is, however, another difference. In the singlet-model (5) with tiny y, N -particles are in very good approximation only produced in φ-decays. The contribution toΓ φq from scatterings is suppressed by additional powers of the coupling y and f N . 8 Physically this means that N -particles effectively do not interact with the plasma any more after they have been produced. Mathematically it implies that the free spectral density (20) is a good approximation for ρ N . In the model (50), on the other hand, the Φ-particles and Ψparticles are charged. At sufficiently large T , their number density is large enough that they scatter frequently, and one expects that N -particles can be produced in scatterings (rather than decays). Mathematically this is reflected by the fact that the expression (54) for ρ Ψ does not only consist of the pole-contribution (56), which is equivalent to (20), but also the continuous part (57), which we have ignored so far.

Full damping rate and scatterings
As pointed out above, the decay Φ → N Ψ is Pauli suppressed for T > m Φ , leading to a factor 1/2 suppression ofΓ Φq for f N 1 and a bigger suppression for f N ∼ 1. 9 On the other hand, the continuum contribution from ρ cont ± grows with temperature. It is strongly dominated by the ρ cont − (ω k − q 0 )-term. This can be understood physically because this contribution comes from Landau damping due to scatterings, and the number of possible scattering partners in the bath increases. To incorporate all these effects, we have to solve the integral in (71) numerically. The result is shown in Figs. 9 and 10.
Regimes and scenarios -One can clearly distinguish different temperature regimes. For T m Φ , the vacuum decay rate is an excellent approximation. For m Φ T < m Φ /α the main effect of the plasma is the suppression ofΓ Φq due to Pauli blocking, and (88) is a good approximation. For M Φ m Φ the thermal correction (77) to the Φ-mass kicks in and enhancesΓ Φq again. (88) remains to be a good approximation. What happens at larger temperatures depends on the value of λ.
Let us first set λ α 2 . This situation is shown in Fig. 9. In this case the hierarchy M Φ M f holds at all temperatures because the thermal Φ-mass is larger than the asymptotic Ψ-mass. Then (72) can be applied even at large temperatures, (88) remains a valid  approximation at T > m Φ /α, and the discussion following (97) applies. That is, (92) and (95) can be used to compute f N p for all temperatures if λ α 2 . 10 The situation λ α 2 is shown in Fig. 10. In this case, the ratio M f /M Φ approaches unity for large temperatures. The approximations for Ω + and Z + in (72) start to break down at T ∼ m Φ /α, which can be seen from the fact that the blue line in Fig. 10 deviates from the gray dashed line. The reason is that (72) is based on the assumption that the produced Ψ-quasiparticles have large momenta p M f = αT /2 with respect to the plasma, which is only true for M Φ M f , m N . If this is not fulfilled, the pole contribution becomes sensitive to the details of the fermion dispersion relations Ω ± in the infrared. At the same time, the contribution from the continuum part ρ cont ± , which is negligible for all lower temperatures, starts to dominate. This contribution can be interpreted as Φannihilation in scatterings. The resulting total rateΓ Φq grows linear with T (as expected due to the increasing number density of scattering partners), but with a coefficient that is larger than predicted by (88). We have verified this for different choices of the parameters, but it is not clear at this stage if a simple analytic expression for the coefficient can be found. Note that Fig. 10 shows the result for the zero-mode q = 0. For q = 0 the validity of (72) extends to larger temperatures because the decay product's momenta in the plasma rest frame are larger. For small values of λ, our simple results therefore hold only up to temperatures T ∼ m Φ /α; at higher temperaturesΓ Φq is dominated by scatterings, and the integrals (69)- (71) have to be evaluated numerically. Since the N -momentum cannot be reconstructed from q in a scattering, one cannot use (95) to calculate f N p . One can, however, still use a rate equation (51) to compute the total DM density. If one is interested in the DM momentum distribution in the regime where N -particles are predominantly produced in scatterings, it would be more convenient to directly evaluate the fermionic self-energies Σ ≷ N with nonequilibrium propagators in the loop to determine the production rate from (134). This goes beyond the scope of the present work, which is focused on DM scalar decays.
A comment on scatterings -We have used resummed thermal Ψ-propagators to account for thermal corrections toΓ Φq . This does not only allow to take quantum statistical factors and the correct quasiparticle dispersion relations in the decays Φ → N Ψ into account, but also includes some of the contributions toΓ Φq due to scatterings. These are encoded in the continuous part of ρ Ψ . In the HTL-approximation (57), this continuous part is non-zero only for x 2 < 1. The resulting contribution toΓ Φq can be interpreted as Landau damping due to logarithmically enhanced t-channel scatterings [112]. The full ρ Ψ 10 One may wonder why the decay approximation holds for arbitrarily high temperatures. Since the number density of scattering partners grows with T , one would expect that Φ-annihilation in scatterings should dominate over Φ-decay at sufficiently high temperature. One physical explanation is that the Φquasiparticles simply are too short lived due to their large thermal mass to find a scattering partner before they decay. For λ = 0, the thermal mass and scattering rate are determined by the same gauge coupling constant α, and scatterings indeed become more efficient at some temperature, c.f. dotted line in figure  10. For λ α the large thermal mass (77) ensures that the decay rate always remains larger than the scattering rate.
(beyond HTL) would also give a contribution for x 2 > 1, which corresponds to s-channel scatterings. These are, however, not logarithmically enhanced. Therefore the use of the HTL approximation corresponds to a systematic expansion in the gauge coupling to leading log in α [112].
Additional contributions from scatterings have been studied by different authors [112][113][114][115][116][117][118] for the case that Φ is identified with the SM-Higgs and m N m Φ (while we assumed m N m Φ ). In the high temperature regime T M Φ m N , m Φ , the difference in the hierarchy of vacuum masses m N and m Φ should be negligible because the thermal mass corrections dominate the kinematics. A direct comparison is nevertheless difficult because the authors of the above works numerically calculated the inclusive rate Γ N under the assumption that Φ is in equilibrium. This situation is e.g. realised in several leptogenesis scenarios. For the case of sterile neutrino DM production considered here, the phase space distribution f Φ for Φ is a dynamic quantity that in general deviates from equilibrium (in particular in the "freeze-in" scenario). Even if we fix f Φ to a Bose-Einstein distribution and obtain a total production rate from (51), the result still depends on the choice of gauge group under which Φ is charged. While it is straightforward to generalise the analytic result (97) we found from U(1) interactions to an arbitrary SU(N) gauge group by making the replacements discussed after (58), there appears to be no simple way to obtain the contribution from multiple scatterings for general SU(N) interactions and non-equilibrium Φ from the numerical results in the literature. The discussion of the approximations made in [112,114], however, suggests that the contributions from vertex-and ladder-diagrams are not logarithmically enhanced, and that our treatment is accurate to leading log order in the gauge coupling α. In conclusion, the above considerations suggest that calculation ofΓ Φq is accurate to orderỹ 2 in the Yukawa coupling and order α 2 log α −2 in the gauge coupling.

Discussion and Conclusions
We have studied thermal corrections to the production rate of singlet fermions N in the decay of neutral scalars φ and charged scalars Φ in a hot plasma. This rate determines the total abundance and momentum distribution of Dark Matter particles in scenarios where sterile neutrino Dark Matter is produced in the decay of heavier particles. With some modifications that account for the different spin, our results may be generalised to the production of other DM candidates in decays, including gravitinos [119][120][121][122][123] or WIMPs [125][126][127].
We found that the sterile neutrino production rate receives considerable thermal corrections if the plasma temperature T exceeds the mass of the decaying scalar. Our main results include expressions (47) and (95), which can be used to determine the abundance and momentum distribution of heavy neutrinos produced in decays. They are based on the gain-and loss rates (37)-(46) for singlet decays and (86)-(94) for charged scalar decays. Our methods is not suitable to compute the heavy neutrino momentum distribution when the heavy neutrinos are produced in scatterings. It is, however, still possible to determine their total number density in such scenarios from (12), (51) and by the methods described in Sec. 3.3. If the scalar is a gauge singlet and decays as φ → N N , then the N -particles in the final state are not significantly affected by the plasma as long as their occupation numbers remain well below the equilibrium value, which is the case in many scenarios studied in the literature. In this case the main thermal correction comes from the fact that φquasiparticles in a thermal bath are screened and pick up a thermal mass. The thermal mass correction appears even if φ is a gauge singlet because it must have some interactions to be produced in the early universe. The effective mass M φ in the plasma grows with T and the coupling constant λ; for T m φ / √ λ it is much larger than the vacuum mass m φ . Since the decay rate is larger for heavier quasiparticles (Γ φq ∝ M φ ), this enhances the DM production rate significantly. Equation (45) gives the thermally corrected production ratẽ Γ φq in this case. The effective mass M φ therein can be determined once the φ-interactions are specified. The momentum distribution of the DM can be calculated from (47).
If the scalar Φ is charged under some gauge group, then at least one of the final state particles in the decay (which we call Ψ) must also carry a charge. This particle is typically in equilibrium in the early universe. This leads to Pauli blocking as soon as the temperature exceeds the scalar mass (T > m Φ ). The behaviour at higher temperatures depends on the effective masses M Φ and M f of Φ and Ψ. If M Φ M f at all temperatures, then the dominant contribution toΓ Φq at all temperatures comes from the decay Φ → ΨN , and the decay products with energies Ω Ψ Ω + and Ω N ω N are always relativistic. In this case the rateΓ Φq is given by (92), as illustrated in Fig. 9. The behaviour can easily be understood as the interplay between Pauli blocking, the enhancement of the decay rate due to the effective Φ-mass M Φ and the increased lifetime of a particle that moves in the plasma due to time dilatation (to be calculated with thermal masses):Γ Φq is suppressed for T ∼ m Φ due to Pauli blocking, but grows linearly with T for M Φ m Φ . This situation is usually realised if Φ has some other interactions (in addition to the gauge coupling that increases the mass of both, Φ and Ψ). We performed the calculation explicitly for an U(1) gauge interaction of Φ. At the order in perturbation theory considered here, it is straightforward to obtain the results for each component of a general SU(N) multiplet from (88), (92) and (97) by making the replacements described after (58). Hence, analytic results allow for a simple inclusion of thermal corrections in the regime where N is predominantly produced in quasiparticle decays of U(1) or SU(N) charged scalars. The DM momentum distribution can be obtained from (95).
If at some temperature M f becomes comparable to M Φ , then the Ψ-quasiparticles that are produced in the Φ → ΨN decay are not relativistic. In this case the decay rate becomes sensitive to the complicated infrared behaviour of fermion dispersion relations in a plasma. For M f > M Φ the decay may even become kinematically forbidden. At the same time other processes, such as inelastic scatterings, contribute to heavy neutrino production. In this scenario it seems difficult to obtain an analytic expression for the scalar damping rateΓ Φq . Qualitatively the behaviour is similar to (97), as illustrated in Fig. 10:Γ Φq is suppressed for T ∼ m Φ due to Pauli blocking, but grows linearly with T for M Φ m Φ . However, we cannot easily determine the coefficient in the relationΓ Φq ∝ T in the high temperature regime. The precise value of this factor depends on the gauge group, the coupling constant and the thermodynamic state of Φ, and it can only be determined numerically.

A Nonequilibrium quantum field theory
Correlation functions in a medium -The usual methods to calculate S-matrix elements are not suitable to describe nonequilibrium systems at large density because there is no well-defined notion of asymptotic states, and the properties of quasiparticles in a medium may significantly differ from those of particles in vacuum. However, observables can always be expressed in terms of correlation functions of the quantum fields, without reference to asymptotic states or free particles. There are two independent two point functions for each field. For a real scalar field φ these are often chosen to be the connected Wightman functions where the . . . is to be understood in the sense of the usual quantum statistical average . . . = Tr( . . .) of a system characterised by a density operator . For a fermion Ψ the corresponding definitions are Here α and β are spinor indices, which we suppress in the following. We will mostly be concerned with fermions, and the Wightman functions S ≷ will be needed to calculate Feynman diagrams. Their linear combinations have an intuitive physical interpretations. The spectral function S − roughly speaking characterises the spectrum of quasiparticles in the plasma. If Ψ is weakly coupled to a thermal bath, then S − (x 1 , x 2 ) = S − (x 1 − x 2 ) depends only on the relative coordinate even if Ψ itself is out of equilibrium and S + is time dependent [95]. Then we can express S − as the Fourier transform of the spectral density ρ(q), which can be viewed as a definition for ρ(q). In a weakly coupled theory, the pole structure of ρ(q) determines the dispersion relations in the plasma: The real parts of poles correspond to quasiparticle energies, the imaginary parts are the thermal widths of the quasiparticles, see (122) and (123). The statistical propagator S + provides a measure for the occupation numbers. The correlators fulfil the symmetry relations which can be seen from the definitions. If Ψ is a Majorana fermion, then there is an additional symmetry where C is the charge conjugation matrix. S − has the boundary condition which follows from the canonical anticommutation relations. The boundary conditions for the statistical propagator are determined by the physical initial conditions in which the system is prepared.
Kadanoff-Baym ansatz -It is reasonable to assume that all fields with gauge interactions are in thermal equilibrium at the time when N gets produced. We will adopt that assumption in this work. The scalar φ may or may not be in equilibrium, depending on whether N production during freeze-in is relevant. Nonequilibrium propagators for scalars have e.g. been studied in [95,[128][129][130]. We will, however, not need them here because φ only appears as external particle. In thermal equilibrium at temperature T , 11 the correlation functions only depend on the relative coordinate x 1 − x 2 , 12 and they are related by the Kubo-Martin-Schwinger (KMS) relations, which are most conveniently written in terms of their Fourier transforms. In absence of chemical potentials they read These allow to express the momentum space two-point functions in terms of the spectral densities, Here f F (q 0 ) = (e q 0 /T + 1) −1 is Fermi-Dirac distribution, which naturally arises as a consequences of the boundary conditions for the correlation functions. The KMS-relations (110) in principle do not hold for N because the sterile neutrinos are not in thermal equilibrium. Explicit expressions for the propagators of Majorana fermions out of thermal equilibrium have been found in [132][133][134][135][136][137], but lead to rather complicated calculations when being used in loop integrals. Since N couples weakly to φ and the thermal bath, the time scale 1/Γ N p on which the occupation numbers evolve is slow compared to the time scale of microscopic processes in the plasma. This justifies to assume that a relation of the form holds locally, which is known as Kadanoff-Baym ansatz. Here f N is a function that can be interpreted as N -distribution function. It changes on macroscopic time scales 1/Γ φq . It is known that this ansatz does not exactly hold [128,129], but it should be sufficient for the present purpose. The Kadanoff-Baym ansatz allows to postulate relations similar to (110) also for N locally in time, in spite of the fact that N is not in equilibrium, By using a single scalar function f N , we already assumed that the different helicity states of N all have the same occupation number, which seems reasonable. The function f N must fulfil the relation f N (−q 0 ) = 1 − f N (q 0 ), (113) which follows directly from the Majorana condition (107) and does not rely on any assumption about the phase space distribution function of N -particles. This can alternatively be found by demanding that there is no charge associated with N present in the plasma, which corresponds to the requirement that dq 0 2π tr[γ 0 S + (q)] = 0.
If evaluated at on-shell values q 0 = Ω N q > 0, where Ω N q is the mass shell of N -(quasi)particle, the function f N (Ω N q ) can be interpreted as the physical phase space distribution function f N q of these quasiparticles in the plasma, Brownian motion -If φ is weakly coupled to a thermal bath, then dissipative effects can be described by a damping rate Γ φq , which can be extracted from the pole structure of the scalar spectral density For a scalar φ that couples to a plasma in equilibrium, ρ φ (q) at leading order can be expressed as even if φ itself is not in equilibrium [95]. Here the retarded self-energy Π R φ (q) depends on T . LetΩ q be a pole of ρ φ (q), with Ω φq ≡ ReΩ q and Γ φq ≡ 2ImΩ q .
In weakly coupled theories one observes the hierarchy Γ φq Ω φq (119) and can make the Breit-Wigner approximation with the residue Then the dispersion relations can be obtained by solving the equation and the φ-damping rate Γ φq that we aim to calculate is given by the imaginary part of the retarded self energy at the quasiparticle pole q 0 = Ω φq , Γ φq and Ω φq depend on T due to the temperature dependence of Π R φ (q). The appearance of the imaginary part can be interpreted as finite temperature version of the optical theorem. The rate (123) is the thermal damping rate for φ-quasiparticles. The contribution to it from Feynman diagrams involving N can be used to determine production rate for Nparticles. We have assumed that the screened single particle states are the only relevant quasiparticles, i.e., collective scalar excitations [138] play no role. If the N -particles are produced from the decay of a coherent condensate ϕ ≡ φ , the rate (123) is still applicable if this production happens during oscillations near the minimum of the effective potential V(ϕ), see e.g. [93]. Analogous to ∆ ≷ , one can introduce self-energies Π ≷ φ 13 and define This allows to rewrite where are the gain and loss rate for φ-particles, which in total give the thermal damping rate Γ φq . The thermal width Γ φq does not only determine the damping of the spectral function ∆ − via (117), it also governs the relaxation of the statistical propagator to its equilibrium form, which in the approximation (119) where ∆ + q;in ,∆ + q;in and∆ + q;in are boundary conditions for the statistical propagator and its derivatives at t 1 = t 2 = 0. The choice ∆ + q;in = 1 Ω φq [ 1 2 + f φq (0)],∆ + q;in = 0,∆ + q;in = Ω φq [ 1 2 + f φq (0)]), which can be interpreted as a quantum state with well-defined quasiparticle occupation number f φq , simplifies (128) to where t = (t 1 + t 2 )/2 and f φq (t) = f B (Ω φq ) + [f φq (0) − f B (Ω φq )]e −Γ φq t is the solution to the Boltzmann equation 13 In the Schwinger-Keldysh-formalism, Π < φ (x 1 , x 2 ) corresponds to a self-energy where the first time argument lies on the "forward" part of the closed time path and the second argument lies on the "backwards" part. For Π > φ (x 1 , x 2 ) it is the other way around.
Here f B (Ω) = (e Ω/T − 1) −1 is the Bose-Einstein distribution. Finally, the modes of the mean field ϕ = φ obey an equation of motion of the form where V(ϕ) is the finite temperature effective potential. Near its minimum Γ ϕq = Γ φq [93]. Hence, Γ φq in this situation is the sole damping scale in the system [129]. With the KMS-relation Π < φ (q) = e −q 0 /T Π > φ (q) it is easy to see that the detailed balanced relation holds and the total relaxation rate Γ φq can be related to the φ-production rate Γ < φq via [92] Γ < φq = f B (Ω φq )Γ φq .
These relations hold to all orders in the couplings amongst the thermal bath's constituents [139]. Analogous to (126) and (127) one can find a damping rate Γ q for fermionic correlation functions S ± , where Ω Ψq is a fermionic quasiparticle pole and Σ ≷ are fermionic self-energies define analogous to Π ≷ φ .
Damping rate far from equilibrium -The expressions following (116) apply if all propagators Feynman diagrams that contribute to the self energies Π ≷ φ fulfil KMS-relations like (109). This is in good approximation the case if φ couples weakly to a larger thermal bath. This is not the case in the present calculation because neither the scalar field not N are in equilibrium in general. Without this assumption, the nonequilibrium correlation functions can be expressed as [128] ∆ − q (t 1 , t 2 ) = sin t 1 where t = min(t 1 , t 2 ). Note that (135) and (136) can be derived from first principles, without the Kadanoff-Baym ansatz (111). Here f φq can be interpreted as the quantum mechanical generalisation of a phase space distribution function. On macroscopic time scales, and if the damping is primarily driven by decays and inverse decays, its evolution is governed by the Markovian equation of motion ∂ t f φq (t) = (1 + f φq (t)) Γ < φq (t) − f φq (t)Γ > φq (t).
Formally this equation can be casted in the same form as (130) by definingf In the special case that φ couples to a large thermal bath, one can easily recover (130) from (138) by using (132). In this casef φq → f B is time independent, and its values for all modes q are fixed by a single parameter T . In general this is, however, not the case.