keV sterile neutrino dark matter from singlet scalar decays: basic concepts and subtle features Academic research paper on "Physical sciences"

lopscience

¡opscience.iop.org

Home Search Collections Journals About Contact us My IOPscience

keV sterile neutrino dark matter from singlet scalar decays: basic concepts and subtle features

This content has been downloaded from IOPscience. Please scroll down to see the full text. View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 187.161.85.189

This content was downloaded on 26/06/2016 at 07:49

Please note that terms and conditions apply.

/ournal of Cosmology and Astroparticle Physics

An IOP and SISSA journal

keV sterile neutrino dark matter from singlet scalar decays: basic concepts

and subtle features j

Alexander Merle and Maximilian Totzauer ^

Max-Planck-Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, 80805 München, Germany o

E-mail: amerle@mpp.mpg.de, totzauer@mpp.mpg.de

Received February 23, 2015 Accepted April 20, 2015

Published June 8, 2015 O

Abstract. We perform a detailed and illustrative study of the production of keV sterile neutrino Dark Matter (DM) by decays of singlet scalars in the early Universe. In the current study we focus on providing a clear and general overview of this production mechanism. For the first time we study all regimes possible on the level of momentum distribution functions, C )

which we obtain by solving a system of Boltzmann equations. These quantities contain |_\

the full information about the production process, which allows us to not only track the evolution of the DM generation but to also take into account all bounds related to the spectrum, such as constraints from structure formation or from avoiding too much dark radiation. In particular we show that this simple production mechanism can, depending on the regime, lead to strongly non-thermal DM spectra which may even feature more than one peak in the momentum distribution. These cases could have particularly interesting consequences for cosmological structure formation, as their analysis requires more refined tools than the simplistic estimate using the free-streaming horizon. Here we present the mechanism including all concepts and subtleties involved, for now using the assumption that the effective number of relativistic degrees of freedom is constant during DM production, which is applicable in a significant fraction of the parameter space. This allows us to derive analytical results to back up our detailed numerical computations, thus leading to the most comprehensive picture of keV sterile neutrino DM production by singlet scalar decays that exists up to now.

Keywords: particle physics - cosmology connection, dark matter theory, cosmological neutrinos, cosmology of theories beyond the SM

ArXiv ePrint: 1502.01011

j Article funded by SCOAP3. Content from this work may be used

A f .v, rdVhe terms,of the Creative Commons Attribution 3.0 License. doi:10.1088/1475-7516/2015/06/011

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Contents

1 Introduction 1

2 Illustrative discussion of the setting and overview of the results 4

3 The kinetic equations in the early Universe 6

4 Analytical results for the distribution functions and relic abundance 9

4.1 The FIMP-regime 11

4.2 The WIMP-regime 12

4.3 Calculating the relic abundance 13

5 Aspects of structure formation 13

5.1 Treatment of the free-streaming horizon 14

5.2 Free-streaming horizon failing — an explicit example featuring twin peaks 15

6 Results and bounds 17

6.1 The Dark Matter abundance 17

6.2 Results for the free-streaming horizon 22

6.3 The dark radiation bound 24

6.4 Other bounds and constraints 26

6.5 Assessing the validity of approximating the relativistic d.o.f. as constant 27

7 Conclusions & outlook 28

A Details on the kinetic equations 29

A.1 Kinetic equation for the scalar 30

A.2 Kinetic equation for the sterile neutrino 32

2 O M U1

B Some background on the numerical calculation of the free-streaming horizon 33

1 Introduction

Despite great advances in our understanding, our Universe holds mysteries that are yet to be resolved. One of the biggest questions is about the identity of the so-called Dark Matter (DM) which — in terms of cosmic energy density — outweighs ordinary matter by a factor of about five [1, 2]. While historically (motivated by, e.g., supersymmetry) our best guess for DM was a Weakly Interacting Massive Particle (WIMP) with a mass of a few hundred GeV and roughly weak interaction strength, by now we have unfortunately not seen a clear signal of such a particle. Even worse, attempts for direct detection disfavour big parts of the parameter space that was deemed "natural" [3-6]. While WIMPs cannot be excluded, we are nevertheless at a point where we should seriously think of alternatives [7].

After all, there are several possibilities left for what DM could be. Alternative ideas range from very light scalars such as axions [8] over non-standard fermions in supersymmetry (such as gravitinos [9] or axinos [10]) up to very heavy exotic options like WIMPzillas [11]. In this work, we concentrate on a candidate motivated by neutrino physics, a sterile neutrino vs of (typically) a mass of a few keV. While our natural guess would be for the sterile neutrino mass to be much larger, it is in reality not bound and there exist several consistent theoretical frameworks in which sterile neutrinos can be very light [12, 13]. Sterile neutrino DM has originally been proposed by Dodelson and Widrow (DW) [14] who suggested that, although sterile neutrinos would not have interactions strong enough to keep them in thermal equilibrium in the early Universe, they could nevertheless be produced gradually by the thermal plasma due to their admixtures to active neutrinos. Although the DW mechanism was at that time consistent with all bounds [15], by now we know that it is excluded by ^^ structure formation: it produces too hot DM [16]. Taking the sterile neutrino mass to larger values is not possible due to the X-ray bound on the DM decay vs ^ v + y, where v is any active neutrino (see ref. [17] for a comprehensive discussion and ref. [18] for probably the most recent collection of limits).1

Several other production mechanisms for keV sterile neutrino DM have been proposed. If a primordial lepton asymmetry is present in the early Universe, the active-sterile transitions could be resonantly enhanced, as proposed by Shi and Fuller (SF) [28]. This mechanism ^___

produces a relatively cold DM component in addition to the thermal DW part, which leads k * to a colder spectrum in accordance with all bounds [17, 29, 30] — however, note that also

this mechanism is in some tension with data if the 3.5 keV line is taken seriously [31]. In theo

ories with an extended gauge group, the sterile neutrinos could equilibrate and be produced via generic freeze-out, if the resulting overabundance is corrected by a sufficient amount of 5 _ entropy production [32, 33] — however this scenario, even though not fully excluded, is on quite general grounds threatened by big bang nucleosynthesis [34].

An alternative is non-thermal production of keV sterile neutrinos via particle decays. C

This possibility is discussed extensively in the literature, with the decaying particle in most 1_^

cases being a scalar: variants range from inflaton decay [35, 36] over a general scalar singlet that could itself either freeze-out [37, 38] or freeze-in [39-41] to the case where the scalar is electrically charged [42]. More general possibilities exist as well, for example the production from pion decays [43], from Dirac fermions [44], or from light vector bosons [45, 46]. A benefit of this mechanism is that the velocity spectrum of the keV neutrinos produced by decays is quite generally known to have a tendency to be colder than that from other mechanisms [31, 38, 39, 47, 48].

While several cases of decay production are discussed, the treatments available in the literature involve some crude estimates and approximations. Even though a non-thermal DM spectrum could have interesting and unexpected consequences for cosmological structure formation, many references work on the level of rate equations and only estimate the spectrum, if at all. Instead, given that the distribution function is the most decisive quantity and that it can be computed from first principles with reasonable effort, its determination should be put

xIn 2014 two groups have independently reported a tentative signal of an X-ray line at 3.5 keV [19, 20] which would, if interpreted as originating from sterile neutrino decay, indicate a particle with a mass of 7.1 keV and an active-sterile mixing angle 6 of roughly sin2 (20) ~ 7 • 10-11. However, up to this date there is still a discussion on-going in the literature about whether or not this line is in fact a real signal, see refs. [21-27]. Given that the signal, if it exists, is at the moment still not at a statistically highly significant level, we do not take any side here but instead suggest to wait until more data is collected, to hopefully either strongly confirm or clearly refute a line signal.

on more solid grounds. This is partially done for the production of keV neutrinos by a general singlet scalar entering thermal equilibrium [38], however, only in the approximation where the effective number of relativistic degrees of freedom g* is constant. While this approximation is not usable in all of the parameter space, it is nevertheless a good approximation for large enough production temperatures and, after all, without this assumption it would not be possible to obtain analytical estimates. Yet, ultimately a fully numerical treatment will be necessary.

In this paper, we will start this endeavour by giving a comprehensive discussion of the production of keV sterile neutrino DM via the decays of general singlet scalars, which themselves are produced via a Higgs portal. As shown in refs. [38-40], depending on the coupling there exist different regimes where the scalar itself could e.g. be a WIMP or a feebly interacting massive particle and decay in or out of thermal equilibrium. We will generalise the previous treatments and derive the full set of approximate formulas for all regimes possible, but again under the assumption of g* = const. We furthermore perform a fully numerical study of the allowed regions in the parameter space and obtain distribution functions for

all relevant parameter points, which in principle allows us to determine all relevant DM properties for a given point in the parameter space. We determine the regions where the correct DM abundance is obtained for a given sterile neutrino mass, thereby taking into

account all relevant bounds. ^___

Structure formation properties of DM candidates can be estimated using the so-called free-streaming horizon Afs. With this quantity, we show that different estimates are possible which may lead to quite different results. We clearly illustrate that AFS, in spite of being a

standard estimator, can lead to inconsistent results depending on the approximations used

to calculate it.

The current paper serves as an illustration of the general principles and its purpose is to give a clear overview of the relatively complicated and subtle details involved. As such, some aspects lie beyond the scope of the current work. For example, we neglect active-sterile o mixing and thus completely disregard DW production, which in reality cannot be switched off for non-zero active-sterile mixing (and which on the contrary is desired for the related X-ray signal), in favour of analytical results. Our approximation is still good for very small active-sterile mixing, but it is nevertheless too simplified and will be dropped (as well as the assumption g* = const.) in a future purely numerical investigation of all possible cases [49]. Also the full set of implications for cosmological structure formation cannot be obtained using Afs, so that a numerical computation of the structure formation properties is needed. While we could in principle add this to the current work, it would lead away from our main point and furthermore prolong the paper even more, so that we have decided to decouple this study, too [50]. Our goal is to ultimately deliver a fully comprehensive study of scalar decay production of keV sterile neutrino DM, in order to set the stage to discriminate it from alternative mechanisms by future data. Current bounds indicate that this may be possible in the not-too-far future [31], so that this endeavour should be pursued in particular if the 3.5 keV line signal survives.

This paper is organised as follows. We start with an illustrative discussion in section 2, which is supposed to give an overview of our considerations and results at a relatively nontechnical level. We then introduce the main equations to be solved in section 3. In section 4, we present analytical approximations for all cases where they can be obtained. After a discussion of the aspects relevant for cosmic structure formation, section 5, we finally present the numerical results of our study in section 6, along with a discussion of all bounds and of

the validity of our considerations. We conclude in section 7. Throughout the paper we try to keep the discussion on a minimally technical level, however, all technical details relevant for an inclined reader are exposed in appendices A and B.

2 Illustrative discussion of the setting and overview of the results

This section mainly serves the purpose of describing in simple terms what we have done in our study. Before entering the technical details, we give an overview not only of the setting we are working in but also of some illustrative results. This will hopefully make the paper more accessible and prevent the reader from getting lost in the technical details.

Our basic setting consists of a real singlet scalar field S which must somehow couple to |

the right-handed neutrino fields N. The most generic coupling doing this job is a Yukawa term 2SNcN with coupling strength y which, if the scalar develops a non-zero vacuum expectation value (VEV) (S}, leads to a Majorana mass mN = y(S), where we have assumed only one generation of right-handed neutrinos for simplicity. Thus, the minimal set of ingredients we need in addition to the standard model (SM) consists of exactly these fields. Our minimal Lagrangian is

L = LSM +

iN/N + 1(0^S)(6»S) - ySNcN + h.c.

^scalar + , (2.1)

Vscaiar = -¿hHfH - 1 S2 + Ah(HfH)2 + AjSS4 + 2A(HfH)S2. (2.2)

This potential could easily lead to a VEV (S} « GeV — TeV, thus motivating a relatively small Yukawa coupling y ~ 10-9-10-5 in order for the mass of the sterile neutrino to be in the keV-range.

Not only can the above Yukawa coupling lead to a sterile neutrino mass, but it will also be responsible for sterile neutrino DM production. Several processes could be thought of: the leading contributions are the reactions SS o hh and S ^ NN, the first of which couples the scalar field S to the thermal plasma and the second of which produces sterile neutrinos

2Given that the astrophysical X-ray bounds push the active-sterile mixing down to tiny values [17, 18], this approximation is not necessarily bad. However, an additional contribution from the DW production can modify some of the results obtained here, which is why we will drop this assumption in future works and turn to a purely numerical treatment [49, 50].

3For possible issues related to the breaking of this symmetry by a non-vanishing VEV (S), see [39] and

references therein. Note also that giving up the Z4 symmetry allows for extra terms in the Lagrangian,

resulting in more processes that can contribute to equilibrating the scalar. In this paper, we stick to the symmetry in order to obtain a minimal parameter space for our exploratory analysis. Considerations of taking into account terms linear and cubic in S in the potential can be found in [38].

which is very similar to what had been used previously [38-40]. Here, Lv denotes the part of the Lagragian giving mass to active neutrinos. The details of this mass generation are in fact not very important for our DM production mechanism, as is the number of fermion generations. However, in the current work, we make the additional assumption of vanishing

active-sterile mixing. In most settings, there will be couplings between active and sterile „___

neutrinos that in reality cannot be switched off. But, since in this paper we want to present analytical results wherever possible, we take this mixing to be exactly zero and in turn have

no DM production from the DW mechanism.2 The scalar potential VScalar takes the most |—1 general form compatible with an assumed global Z4-symmetry:3

from the decay of S. In principle, also processes like SS ^ NN would be possible, but the corresponding rate is proportional to y4 which is negligible for a sufficiently small Yukawa coupling y. We also neglect the inverse reaction NN ^ S which is suppressed due to the heavily suppressed phase space originating in the kinematics of any 2-to-1 process.

Let us add that, in general, there will be mixing between the scalar S and the SM Higgs field after electroweak symmetry breaking, which we have completely ignored. However, given that this mixing is proportional to the generally small Higgs portal coupling A and that it is also suppressed if the singlet scalar is considerably heavier than the Higgs, which is the case that we are investigating here (cf. section 4), taking into account the mixing would not at all change our results. Note that this would change if one used singlet scalar masses very close to or below the Higgs mass [40], which is why this simplifying assumption must C | be dropped if we are to extend our considerations to lighter singlets. On the other hand, in that limit we would need to drop future assumptions used in this work (in particular the assumption that g* is constant), so that it makes sense to postpone these considerations to further work [49].

With these assumptions, it is clear that every scalar present in the early Universe will ^^ either decay into two sterile neutrinos or undergo pairwise annihilation into pairs of Higgs bosons. Depending on the exact values of the couplings, there are different regimes possible. CTi For example, if the Higgs portal coupling A is small enough and the initial abundance4 of the scalars is zero, then the scalar will never enter thermal equilibrium but it will only be produced occasionally from the plasma. In a more modern language, this mechanism would be called freeze-in [51, 52], and the corresponding particle would be called a feebly interacting o massive particle (FIMP). In this regime, the annihilation into Higgs bosons can be neglected |—1 since its reaction rate will be suppressed by the square of the (tiny) scalar density. Hence, the frozen-in abundance of scalars will ultimately be translated into a relic abundance of sterile neutrinos with a particle number just twice the one of the scalars at freeze-in (in a co-moving volume) — irrespective of the size of the Yukawa coupling y as long as it is large enough that all scalars have decayed by now. Another example would be the case where the scalar couples to the Higgs strongly enough to be in thermal equilibrium. Then the argument of doubling the number of particles would still be valid if the decay width of the scalar is sufficiently small for the decay to become effective only after freeze-out. However, if the decay proceeds already while the scalar is in equilibrium, further contributions will be present making the whole picture more complicated. We will discuss all cases in detail in section 4, where we present analytical estimates for the momentum distribution functions f and yields Y wherever possible.

We later on turn to a numerical computation of the DM production. To arrive at the final DM relic abundance, we need to solve a system of Boltzmann equations to compute the momentum distribution function f (p, t) of the DM particle. Ultimately this distribution function contains all relevant information: one can e.g. use it to compute the final DM relic abundance, but it also encodes the information about the evolution of the momentum spectrum with time, i.e., how many particles exist per momentum interval at any given temperature of the Universe. This allows to not only track the course of DM generation in detail, but it furthermore gives information about the velocities of the DM particles at the epochs relevant for the formation of structures (i.e., galaxies and galaxy clusters) in space.

4We usually use the term abundance to denote a particle number density. Depending on the context, the term relic abundance will be used for the particle number density or the corresponding energy density after the process discussed is complete.

To give a snapshot of the results to be presented, we show in figure 1 a plot of the lines of correct DM abundance for several different sterile neutrino masses, in a plane of the parameters CHp and Cr which, as we will explain later, are nothing else than an effective Higgs portal coupling and an effective decay width in convenient units. We have augmented the plot by some example evolutions of the DM production and the underlying distribution functions for some specific parameter points marked in the central figure, in order to illustrate what is behind our computations. The purpose of figure 1 is to give a graphical illustration of what will be presented in this paper. All the plots displayed will be discussed in great detail in section 6, where also the terminology, colour-codes, and labelling will be carefully outlined so that, while advancing with the paper, the reader will ultimately be enabled to get a full understanding of figure 1.

3 The kinetic equations in the early Universe

L[f] = C[f]. (3.1)

The fundamental embodiment of every particle species in the early Universe is their distribution function f (p, t) in momentum space. This quantity does contain all the information we need to deduce the cosmological impact of the species under consideration, so that the determination of f is paramount. Due to isotropy, we will always assume that the distribution functions f only depend on the moduli of the associated 3-momenta. To compute a ^—^ distribution function, we need to solve the corresponding Boltzmann equation:

The left-hand side of this equation contains the so-called Liouville operator L:

d d -----

L = ^ - Hp—, (3.2)

where p is the modulus of the particle's 3-momentum and H is the Hubble function. The collision term C[f] on the right-hand side can be interpreted as a source term, encoding the interaction of the species of interest with itself and with the other particle species present in the plasma. This term contains all the information about the processes which contribute to the production of the species under consideration, and which accordingly shape the resulting momentum distribution function. Collision terms can be relatively lengthy, which is why we report the explicit form of C[f] only in appendix A. In this section, however, we prefer to give a more intuitive explanation of how to obtain the distribution functions of the sterile neutrinos.

Since DM production happens in the very early Universe, we only need to consider the era of radiation dominance. During this epoch, the Hubble function can be written as H(T) = T2/M0, where T is the temperature of the plasma and M0 is a function that implicitly depends on time via the evolution of the degrees of freedom g*:

Mo = (4^)1/2 = 7.35g-1/2 x 1018 GeV. (3.3)

Introducing dimensionless variables,

x = p/T and r = ms/T, (3.4)

1010-'1 0.001 0.010 0.100 1 10 100

r=ms/T

Figure 1. Lines of correct abundance with example distributions and evolutions.

the Liouville operator from eq. (3.2) can be recast into the following form:

r dr dT d ( jr. dr ^ \ a

L = dTdt^ + 1 -j.i^J fix- (35)

Throughout this work, we will stick to the assumption that the numbers of effective degrees of freedom (g*, g*S) are constant until the production of sterile neutrinos is completed, as done e.g. in [38, 39]. This assumption is absolutely necessary to obtain analytical results and, in fact, it is not a too bad approximation in a large fraction of the relevant parameter space, as we will illustrate in section 6.5. Nevertheless we will drop it in later works where we are going to present more realistic and purely numerical studies [49, 50].

If the number of effective degrees of freedom does not change during the period of interest, the advantage of the variables x and r becomes obvious. Accordingly, the dynamics o of the scalar and the sterile neutrino are given by the following set of equations:

dfs = dT dt s + Cs + Cs ) (36)

= ini^F^hh-rSS + CSS^hh + CS^NN) , (36)

Or dr dTV hh^SS ' SS * 'N N

dfN = dT dt CN (3 7)

=dT dTCS^NN . (3.7)

In eqs. (3.6) and (3.7), the upper indices on the collision terms mark the species the kinetic equation of which is governed by this term, while the subscripts describe the actual process. Thus, CS^hh describes the depletion of scalars due to annihilations into pairs of Higgses |—

while CSh^SS describes the reverse process. In turn, CS^NN describes the depletion of |_\

scalars due to decays and CSN NN encodes the creation of sterile neutrinos from the decays of scalars. Note that the collision terms contain information about the kinematics, too, so that CSN NN and CSS NN differ by more than just a sign. For a detailed derivation and explicit expressions, see appendix A.

Since we approximate g* as constant during the time of production and since the matterradiation equality only takes place at temperatures of O (1 eV), i.e., long after DM production is completed, we consistently use the time-temperature relation dT = -HT in eq. (3.7). Using this as well as the explicit form of the collision terms as derived in appendix A.1, we find the kinetic equation for the scalar:

f S (x, r) dr

yj x2 + r2

x2 + r

—Chp exp (-jx2 + r2)F(x,r) -Cprf(x,r)

00 a max

- -1 CHPfS(x, r) 2n J dx x2 J dcos 6fS(x,r)G(x, cos d, r)

= Q(x, r) - P(x, r)fS - R(x, r) Ir [fS] fS ,

where we have defined

Chp exp (-Vx2 + r2)F(x,r)

Q(x,r) =-^— '-, (3.9)

4nvx2 + r2

P (x,r) = 1 , (3.10) Vx2 + r2

R(x,r) = ^H^, (3.11)

4^x2 + r2

» «max

Ir[fs] = y dX X2 J dcos dfs(x,r)Q(X, cos 0,r). (3.12)

with Xmin = ||x — r/2/(4x)||. It is this equation that will be absolutely central for us: once we have managed to determine the distribution function fs of the scalar from first principles, we only need to plug it into eq. (3.14) to determine the distribution function fN of the sterile neutrino. As to be expected, the distribution function fs of the scalar will look differently depending on which regime we are looking at (e.g., the scalar being an early decaying WIMP compared to a late one). This will be translated by eq. (3.14) into different resulting sterile neutrino distribution functions. Note that, transforming the momentum variable into an energy, eq. (3.14) perfectly coincides with [35, eq. (8)] and [38, eq. (10)]. Note also that the master equation for the sterile neutrino distribution is decoupled from the one for the scalar only because of our assumption that the reaction NN ^ S is negligible.

4 Analytical results for the distribution functions and relic abundance

In this section, we present the limiting cases that can be treated analytically. Some of the results can be found in the literature, see in particular ref. [38], while others are completely

For the explicit forms of the kinematic functions F and G and the definition of amax, see appendix A.

In eqs. (3.8) to (3.11), we have used two important dimensionless auxiliary quantities:

| the effective decay width: Cr = M m; , (3 ^

[the effective (squared) Higgs portal: Chp = t^ . ^

In the remainder of this work it will turn out convenient to use these quantities to span the parameter space of our setting. Note that, during the DM production process, we assume Mo to be constant by virtue of eq. (3.3). Hence the interpretations of Chp as an effective ISJ Higgs portal coupling and of Cr as an effective decay width are appropriate, in the sense Q that the dependence of M0 on g* does not change their behaviour and hence they indeed play practically the same roles as the fundamental Lagrangian parameters behind them. For simplicity, we may often refer to CHP as Higgs portal and to Cr as decay width.

Turning to the sterile neutrino distribution function, which is our main quantity of interest, we can use the assumption of neglecting the back-reaction NN ^ S to find a very simple form for the corresponding kinetic equation. Using eqs. (3.6) and (3.7), see

appendix A.2 for details, one can derive a very intuitive master equation for the distribution function of the sterile neutrino in terms of that of the scalar:

rr J2 r ro X

fN(X,r)= dr/2Cr — dx -_fs(X,r/), (3.14)

Jo X iimin VX2 + r/2

new. Basically there exist two main cases, the scalar itself can either enter thermal equilibrium in the early Universe and thus act similarly to a WIMP or it can be only very feebly interacting, thus undergoing freeze-in like a FIMP:

1. The FIMP-regime. This limiting case is characterised by a small Higgs portal Chp and an (almost) arbitrary decay constant Cr. The scalars that are produced via freeze-in subsequently decay into sterile neutrinos. However, for the abundance itself the time of the decay is not very relevant as long as it happens before matter-radiation equality, which is why the exact value of Cr does not matter much in this regime.5 It will, however, play a role for the spectrum itself, cf. section 6.

2. The WIMP-regime. For large enough Higgs portals CHP the scalar will thermalise, i.e., f" | it enters thermal equilibrium and any information about its initial abundance is lost

since the abundance will always be the thermal one. The particle remains in equilibrium until the interaction rate for the process hh o SS drops below the expansion rate of the Universe and the scalar decouples from the thermal bath. During the course of these events, depending on the exact value of the decay width Cr, the scalar can decay into sterile neutrinos at various stages:

5Note that, for this regime, we will always assume a zero initial abundance for the scalars (and trivially for the sterile neutrinos). If there was a non-negligible initial abundance, this would change our results since the primordial scalars would then add to the ones produced via freeze-in. However, given that there is no reason for such an initial abundance to be present and that we do not see much value in speculating how it could possibly have been produced, we stick to the conservative viewpoint and only produce scalars from the freeze-in mechanism itself.

On the other hand, one could argue that sterile neutrinos and/or singlet scalar fields could quite generically couple to the inflaton field (see refs. [35, 36] for examples concerning the former case). In such scenarios, assuming that inflation is the correct theory in the first place, our scenario might be modified considerably. However, such couplings are only compulsory if the SM gauge group and all other low-energy symmetries are the only ones up to very high scales and even then, from a model building point of view, there exist various possibilities to strongly suppress certain couplings, e.g. by locating the various fields on different branes. While our considerations could in principle be extended to include this point, this would add further complications of which it is however unclear whether they exist, or not. We thus stick to the most minimal setting and disregard any primordial abundance of sterile neutrinos and/or the singlet scalar, as well as the related coupling to the inflaton field.

(a) In-equilibrium decay. If the decay width Cr is large enough and if the Higgs portal CHP is large, too, the sterile neutrinos are produced already very early from the decays of equilibrated scalars. The sterile neutrinos arising from the decays of the few residual scalars that are present after freeze-out can then practically be neglected since the large Higgs portal guarantees a small relic abundance of frozen-out scalars.

(b) Out-of-equilibrium decay. The scalar couples to SM particles strongly enough to

enter thermal equilibrium. However, if the decay width is sufficiently small, the —^ production of sterile neutrinos during the time of equilibrium is negligible and it is sufficient to only take into account the late decays of the scalars after their freeze-out. In this regime, the scalar itself

acts as a DM-like species, but it is unstable and decays before it can significantly contribute to the energy budget of the Universe.

(c) Intermediate regime. For intermediate values of the decay width Cr, neither of the above limiting cases is applicable and it is a priori not clear how well the combination of both of them can correctly describe the situation. On the other hand,

this case can be of particular interest since it results into a distribution function with two intrinsic momentum scales. This may open up interesting possibilities to tackle the well-known small scale problems of cosmological structure formation [53-56]. We will therefore treat an explicit example for this case numerically in section 5.2.

We will now discuss these different regimes one by one. In all cases, we will take the limit £ = mh/mS ^ 0, cf. appendix A.1. Even for £ « 0.3, the net effect on the distribution function is of the order of a few per mille. We want to keep the mass mS of the scalar far from the electroweak (EW) phase transition to ensure that scalars S are only produced by Higgs bosons. For masses mS < vEW, a new range of interaction becomes important like the production of EW gauge bosons from scalars S, SS o VV. This scenario is discussed in [40] on the level of abundances, and we will take these modifications into account on a further more technical study where the focus will be put on even more realistic numerical results [49].

4.1 The FIMP-regime

For sufficiently small Higgs portals CHP (i.e., for A ^ 10-6 [38, 39]), the scalar never enters equilibrium and hence one can — for vanishing initial abundance — neglect in eq. (3.8) the term R(r, x) Ir [fS] fS which scales as fS. The remaining kinetic equation is an ordinary differential equation which can be solved analytically. For a vanishing initial abundance of the scalar, the resulting distribution function is given by:

, „ fr , >exp ( —\/p2+x2 )

fs(r,x)= Cup dp ;

J0 VP2+x2

ePy/ P2+x2

P+\/ p2 +x2 er^r2+x2 V r+Vr2+x2

where K1 is the first modified Bessel function of second kind, cf. eq. (A.10). Plugging eq. (4.1) into eq. (3.14) in order to get an analytical result for the distribution function of the sterile neutrino is not very instructive, since there is no simple form of that expression. However, the abundance of the sterile neutrino can be computed analytically for late times, r ^ <x>. Setting Cr to zero corresponds to a stable scalar. With this choice, eq. (4.1) can be integrated rather easily and one obtains the (hypothetical) relic abundance of the stable scalar. Since all frozen-in scalars will however ultimately decay into two sterile neutrinos for a non-vanishing decay width Cr, the abundance of sterile neutrinos will then just be twice the abundance of the would-be-stable scalar [39]. The result for the yield Y = n/s, with n and s the particle number and entropy densities, respectively, is given by

2 O M U1

Yn(r ^ œ) =

135 Chp 64n2 g*(Tprod) '

Here, Tprod denotes the temperature at the time of production.6

6Of course the time of production is subject to some ambiguities in its definition. Both freeze-in/freeze-out of the scalar and its subsequent decay are continuous processes, the time scales of which are determined by Chp and Cr. It is hence convenient to define the production time as the point when the abundance of sterile neutrinos has passed some threshold fraction of the final abundance, which we take to be 95%.

4.2 The WIMP-regime

In-equilibrium-decay. If both the Higgs portal Chp and the decay width Cr are large, sterile neutrinos are efficiently produced from the decays of scalars already while being in equilibrium. This case is covered at length in [38]. However, our analytical expression differs by a constant factor, which is why we will sketch the most important steps needed to derive the result. If the scalar is in equilibrium, we know its distribution function exactly. Accordingly, the authors of [38] use a Bose-Einstein (BE) distribution to capture the quantum nature of the scalar. Since the whole set of equations governing the dynamics was however derived using the Maxwellian approximation in the Boltzmann equation, one might wonder whether it would be more consistent to again use a Maxwell-Boltzmann (MB) distribution for the scalar. We will exemplify both cases in order to illustrate that the difference between them is irrelevant.

Plugging the BE and MB distributions into (A.20) yields:

8Cr / dzx^Z-1 log (y^—r) (BE)

fN l e-*.^ erf ^/X^ ' (4-3)

8Cr / dzx^Z—Ye-xz =-2/X- " e-xz^ ^zT—ï (MB)

where we have introduced the variable zr = r2/(4x2) + 1 for convenience.7 Integrating fN (x,r) over d3x and again taking the limit r ^ to allows to calculate the yield for late times:

Yn (r ^to) = ( SZ (5) ^ (BE) . (4.4)

( ) I itgj-d) (MB) ( )

gw (MB) ■

Both results only differ by a factor of Z(5) « 1.0369, which justifies the use of either distribution. Our result in the BE case is however larger by a factor of 5/2 compared to the one C )

reported in [38]. While one may easily forget powers of two in these computations, we could |_

not trace any step where a factor of 5 could possibly be introduced, making us confident that our results are correct.

Out-of-equilibrium decay. This limiting case describes a scenario where the scalar is in equilibrium and ultimately freezes out. The decay width Cr is so small that practically no sterile neutrinos are produced before the scalar decouples from the plasma. Only after the scalar has frozen-out, the production of the sterile neutrinos sets in. As in [38], we make the approximation that the scalar has a thermal distribution until it freezes out instantaneously at r = rFo. In this case the kinetic equation, eq. (3.8), can be solved:

(,- \Crx2/2 ___

r) e-M^-o V^)/2,

rFO + y r2O +x2/

where feq(x,rFo) is the equilibrium distribution of S at freeze-out. It could again be taken to be BE or, more consistently, MB. In the case of BE, our expression coincides with [38,

7Note that, as we had already mentioned, we have neglected the mixing between the two physical scalars, which is a very good approximation in our case. However, in order to simplify the comparison of our results to the ones obtained in ref. [38], it is of course necessary to apply the same approximation to the results from that reference.

eq. (43)]. The final abundance of sterile neutrinos can in this limiting case again be calculated from doubling the abundance of scalars at freeze-out. Given as a function of rFO, the expression for the yield is

YN(r - = , 4 " [ d« ^ - rF0 ( = ^¿O^ for ,

"( ) 4n4g,(Tp,odn e< - 6 V 4n4g,(Tp,od) ), ()

with 6 = 1 (6 = 0) for the BE (MB) case. Also here the numerical difference between both versions is fairly small for realistic values of rF0.

Intermediate regime. If Cr is in an intermediate regime, neither the production before |

nor the one after freeze-out completely dominates the sterile neutrino distribution function, such that no simple analytical treatment is possible. In this instant we have to rely on a ^^ purely numerical treatment. We will present a sample case for this regime in section 6. Nevertheless this intermediate regime may be of special interest as it features two intrinsic scales.

4.3 Calculating the relic abundance

So far we have shown formulae to compute the yield at r ^ To convert this into the commonly quoted closure parameters, we first have to multiply the yield by today's entropy ISJ density of s0 = 2891.2 cm-3 [57] and then compare the DM energy density today to the critical density of the Universe. We can write

nDM = mN Yn (r:T)S0, (4.7) 5

Pcrit/h2 -----

where Pcrit/h2 = 1.054 x 10-2MeVcm-3 [57]. With this conversion formula at hand, it is o

straightforward to transform the analytical estimates for the yield Yn in sections 4.1 and 4.2 |_\

into expressions for the DM relic abundance, which can then be compared to the observed 1a range obtained by the Planck collaboration, Qdm^2 = 0.1188 ± 0.0010 [2].

5 Aspects of structure formation

Up to now, we have only been concerned with the DM relic abundance, i.e., the amount of DM in the Universe. However, for a viable DM candidate it is not only important to be sufficiently abundant in the Universe, but it must also allow for structures such as galaxies to form. Clearly, this imposes constraints on the momentum distribution function f (p, t) of the DM candidate under consideration: it must not be hot, i.e., it is not allowed to be too relativistic at the time of structure formation when galaxies form (more precisely, the allowed amount of hot DM is at most a tiny 1% of all DM in the Universe [58, 59]).

The form of the distribution function f (p, t), i.e. the spectrum of the DM candidate, can be constrained from the observed matter distribution in the cosmos: the evolution of spatial inhomogeneities depends on the spectrum of the DM particles and therefore it ultimately constrains the DM production mechanism. Since it is computationally impossible to run a simulation of structure formation for any possible distribution function, a commonly used indicator that can be calculated easily for a given f (p, t) is the so-called free-streaming horizon Afs. This quantity describes the average distance a DM particle would have travelled after

expected. Since fN(x) does not change after the production process is complete, the factor of r2mN/m| in the square root increases. Note also that, once r2mN/m| is greater than the value(s) of x around which f is concentrated, eq. (5.2) converges to the non-relativistic expression, (v) ^ T/mN (x) = (p) /mN.

For our numerics, we follow the approximations usually found in the analytical approach [16], namely we will assume the Universe to be completely radiation dominated until some temperature Teq ("matter-radiation equality"), where the Universe switches to being completely matter dominated. The last epoch of vacuum dominance is irrelevant: even though this period dominates the past of the Universe on an absolute time scale (the matter-vacuum equality being located at a redshift of Z — f a/^ m ~

2.2, corresponding to a time

of roughly 3 x 109 a [60]), the velocities of the sterile neutrinos in this epoch are so tiny that the resulting contribution to the free-streaming horizon is negligible. Commonly, the evolution of the degrees of freedom (d.o.f.) is implemented by an additional dilution factor of £1/3 = [#*s(Tprod)/g*s(To)]1/3 « (106.75/3.36)1/3 « 3.17, by which Afs is rescaled to account for entropy dilution (cf. [38, 39]). Note that, although we approximate « const. during DM production, we nevertheless have to take into account the entropy dilution until today, since there is no justification for the above assumption to be valid through the entire history of the Universe. Departing from the above approximation also during DM production would modify the dilution factor £1/3, but not too drastically because of the presence of

production without collisions and not subject to gravitational clustering. In fact this quantity can serve as a good estimator for a length below which the formation of structures is heavily suppressed. For that reason, the free-streaming horizon AFS is commonly used in the literature to discard DM models with a spectrum that is too hot to explain the filamentary structure of the large scale matter distribution in the Universe by preventing galaxy-sized objects of being formed.

5.1 Treatment of the free-streaming horizon

In this section, we will show that the free-streaming horizon itself is subject to substantial uncertainties in its definition, which will make clear why we later present some of our results in two different versions. Moreover we will demonstrate that even our simple one-component C | DM model can produce highly non-thermal momentum distribution functions, which may even feature two intrinsic momentum/velocity scales. In such a case, the average momentum (and hence the free-streaming horizon) cannot be expected to lead to sensible conclusions. The free-streaming horizon is defined by [16]:

. _ fT' (v(T)) dt

Afs = 4«, W dTdT' ("> (

where Tprod is the temperature at which production can be seen as complete (in the sense discussed before) and T0 is today's temperature. Here, (v(T)) is the average velocity of the ISJ sterile neutrinos that can be calculated from the distribution function:

°° 3

I dx — x fN (x, r)

n ■\/x2+r2(mN/ms )2 '

(v(T)) = °-Vroo d 2. ( )-. (5.2)

J0 dxx2fN(x,r) _____

From eq. (5.2) it becomes clear that (v) converges to unity, i.e., to the speed of light, as f ) long as fN is concentrated around values of x = p/T » ^Jr2mN/m| = mN/T, just as

the third root. We will later on perform an estimate of the validity of our approximation, cf. section 6.5. Note that, however, in this formalism the dependence on (Tprod) is quite mild, such that it is safe to use the SM number of d.o.f. = 106.75 even though the new particles contribute as well to some degree, depending on how strongly suppressed their true distribution functions are compared to a thermal one.

There is, however, an issue with the analytical approach. The above treatment does not capture the physical fact that, in reality, the evolution of the d.o.f. also enters in the time-temperature relation inside the integral in eq. (5.1), but this can be taken into account rather easily in a numerical evaluation of the integral. We therefore compute Afs in a second (numerical) version, in order to take the full evolution of the d.o.f. into account. Thereby, our numerics uses a set of fitted analytical formulae for the evolution of the d.o.f. [61]. For |

more technical details on this numerical integration, see appendix B.

In order to compare to the results already present in the literature, we will follow both approaches in parallel. To get an idea about whether the DM can be classified as cold, warm, or hot for a certain set of parameters, we take AFS = 0.1 Mpc to mark the boundary between hot and warm DM. This choice is relatively common in the literature (see, e.g., refs. [6265]), and it corresponds to the size of a typical dwarf satellite galaxy, which yields a sensible physical motivation. The boundary between warm and cold DM in turn is smooth but it is clear that AFS should be considerably smaller in this case, so that we choose AFS = 0.01 Mpc, i.e., one order of magnitude smaller than for the hot/warm boundary. Of course there is some arbitrariness involved in these choices, but given that the free-streaming horizon in ——. itself can only yield an indication, the actual error introduced is not as serious as it may seem at first sight. In general the free-streaming horizon can only serve as an order-of-magnitude estimate, and it clearly should not be used to prematurely discard unclear results. Ultimately,

only more advanced computations of structure formation can assess whether scenarios with „___

borderline free-streaming horizons should be discarded, or not [50]. For recent developments beyond the commonly used free-streaming approach, see e.g. [66]. ^^

5.2 Free-streaming horizon failing — an explicit example featuring twin peaks

As mentioned before, even our simple one-component DM model can feature a distribution function with two intrinsic scales, namely in the case where the relic abundance of sterile neutrinos is produced from the decay of equilibrated scalars and the decay of frozen-out scalars in comparable amounts. We present here the exemplary case of Chp = 104 and Cr = 10-3'5, which yields the correct relic abundance for a (relatively large) sterile neutrino mass of about 73 keV. Fixing the scalar mass to be 1 TeV and the number of degrees of freedom at high tempertures to g* = 106.5, this would correspond to values of the Lagrangian couplings of (y, A) = (4.7 x 10-9, 8.3 x 10-5).

Figure 2 shows the distribution function of the sterile neutrino for different values of the time parameter r. One can clearly see how early production from the plasma populates the lower comoving momenta (dubbed Thermal part — although the resulting distribution may not be perfectly thermal), while the late contributions mainly originate from the decay of the frozen-out scalars (Decay part). The inset in the plot shows the enlarged region between the two peaks.

It is obvious that in this case the average momentum does not at all capture the characteristics of the distribution function. According to our estimates using AFS and the average value (x) = (p) /T ~ 16.6, this point in the parameter space corresponds to a scenario where the sterile neutrinos are on the borderline between being hot and being warm (cf. section 6).

Figure 2. Example of the evolution of a distribution function of sterile neutrinos. One can clearly distinguish two momentum scales (global maximum at x\ ~ 1.5 and second, local maximum at x2 ~ 26). The mean comoving momentum (x) is located at (x) ~ 16.6. The standard deviation sj(x2) — (x)2 is approximately 26, which illustrates that the mean value (x) contains practically no information.

However, in fact the low momentum ("thermal") part with (x)iOw ~ 2.5 contributes practically as strongly as the high momentum ("decay") part with (x)high ~ 35.7, where in both cases we have approximated the respective peaks with the individual results from eqs. (4.3) and (4.5), respectively. Note that this splitting introduces some numerical uncertainties, since the expression in eq. (4.5) is quite sensitive to small deviations in rFo, which in turn suffers from some arbitrariness in the exact definition (due to freeze-out being a process with a small but finite temporal extent rather than an immediate effect).

Estimating the two corresponding free-streaming horizons, cf. eq. (5.1),

( a numerical VAFS,thermal,

\ estimate \

FS,thermal

numerical VAFS,decay ,

estimate FS, decay y

(0.05, 0.01) Mpc (0.7,0.1) Mpc,

the left peak tends to yield cold/warm DM while the right one would indicate hot DM in both cases. The full distribution function corresponds to (AF^Otl?1, ^ftotd) ~ (°.27, 0.05) Mpc, which is perfectly in between the two individual estimates and consistently indicates a case at the borderline of warm/hot DM. Even though the splitting into two distinct parts introduces some extra numerical uncertainty with respect to the values for the complete distribution function, these values clearly illustrate the issue with using the free-streaming horizon as an estimator.

This is precisely one of the cases where more detailed studies about structure formation have to be performed in order to obtain a final answer. With the simple estimate of the freestreaming horizon alone, such a scenario should not be prematurely discarded. In any case,

O H I-1

Figure 3. Lines of correct abundance (dark colours) and of sizeable but insufficient abundance (faint colours) in the CHP-Cr plane for different values of the sterile neutrino mass. In addition, the Tremaine-Gunn and overclosure bounds are displayed, see text for details.

it is worthwhile noting that a one-component DM model can have two similarly important momentum scales in its distribution function, which may open up new possibilities to tackle the small-scale problems from structure formation simulations.

O H1 (J1

6 Results and bounds

In this section we will present our detailed results and we will also address possible bounds on the scenario as well as the validity of our considerations.

6.1 The Dark Matter abundance

Let us first discuss the DM abundance, which is displayed in figure 3 (similarly to the one shown in section 2). In this plot, the dark coloured bands mark the regions in the CHP-Cr plane where a sterile neutrino of a given mass yields an abundance in accordance with the 3a range from the Planck 2015 data [2]. For example, the dark red lines correspond to a sterile neutrino with a mass of mN = 1 keV. Here, the lines on the left correspond to the FIMP regime while the ones on the right correspond to the WIMP regime. We also mark, by the neighbouring fainter lines, the regions where a sizeable but not sufficiently large abundance is generated.

There are two important bounds displayed in figure 3. Let us start with the so-called Tremaine-Gunn (TG) bound [67], which is based on the idea that any collection of identical fermions must have a certain minimum phase space density. Applying this bound to the observed dwarf satellite galaxies leads to a lower bound of roughly mN > 0.5 keV on the sterile neutrino mass [68]. Thus, smaller masses are ultimately forbidden by the Pauli exclusion principle. In our plots this bound is marked by the gray dashed line around the white area

Distribution Function

Chp=10_1 Cr=10_1

<x>^2.67

0.001 0.010 0.100 1 10 100

K therm

Particle number densities for:

10-4 0.001 0.010 0.100 1 10 100

r=ms/T

Figure 4. The evolutions of the distribution function (left) and of the abundances (right) for a point corresponding to the scalar being a FIMP.

which in particular cuts into the parameter space for values of CHP close to one. It basically indicates where the relic density for mN = 0.5 keV equals the upper 3a bound [2]. The second bound comes from the fact that the DM should not "overclose" the Universe, i.e., its energy density fraction QDm should be smaller than one. Since in the figures we marked the lines of correct abundance for different masses mN, the resulting forbidden regions do in fact also depend on the mass of the DM particle. However, given that there is a model-independent lower value for the mass from the TG bound, at least for this smallest mass of mN = 0.5 keV the overclosure region marks an absolute bound. Not too surprisingly, this forbidden region is smaller than that excluded by the TG bound, and in our plots it is marked by the gray patches in the upper right corners.

As already indicated in sections 2 and 4, depending on the exact values of the parameters different regimes are possible. Let us discuss a few numerical examples. Starting with the FIMP case, in figure 4 we illustrate an example point (CHP,Cr) = (10-1,10-1) which corresponds to this regime. Taking again the benchmark value mS = 1 TeV and using g* = 106.5 at high temperatures, the effective couplings translate into (y, A) = (8.4 x 10-8, 2.6 x 10-7) on the Lagrangian level. On the left panel, we depict the evolution of the sterile neutrino distribution function fN(r) with the time parameter r. As one can see, most of the abundance is produced around the time r ~ 1. Soon afterwards the production ceases such that even for very late times, the distribution hardly changes (as soon as r ~ 10, the distribution is practically identical to the final one). The distribution exhibits a clear peak whose maximal value is very close to the mean momentum over temperature, (x) ~ (p/T) ~ 2.67. This means that this sterile neutrino distribution is colder than a thermally produced one, for which this number would be equal to 3.15 [35]. However, this point nevertheless turns out to be in the hot DM region, cf. section 6.2. This is mainly due to the fact that this point in parameter space requires a mass of the sterile neutrino of mN « 2 keV in order to fulfil the relic abundance constraint. This low mass in turn leads to a long time of highly relativistic free-streaming.

INJ O M U1

The evolution of the abundances m, i.e., the integrals over the distribution function divided by T3, is depicted on the right panel of figure 4. We display both the abundances of the scalar and of the sterile neutrino. Starting with the scalar (solid black line) we can see that, as to be expected from a generic FIMP, the abundance of the scalar is gradually built up with increasing time parameter r. However, it never reaches a thermal abundance, as can be seen by comparing the black curve to the hypothetical one for a scalar in thermal equilibrium (dotted gray line). If the scalar was stable, its abundance would not anymore change after the freeze-in is completed, cf. dashed gray line, and it would in practice act as some type of DM. However, given that the scalar is unstable, once it does decay its abundance decreases and instead a sizeable abundance of sterile neutrinos is built up (red solid line). Indeed, because of each scalar decaying into two sterile neutrinos, the final abundance of sterile neutrinos is exactly twice the one of the would-be-stable scalar ^^ for late times (numerically we obtain (r = 250)/HS(Cr = 0, r = 250) ~ 2.03, in excellent agreement with the expectation). Furthermore, we can use eq. (4.2) to cross-check our numerics, and both results agree nearly perfectly with each other, within a deviation of only 2.8% in this case. Note that this value is not a measure of the quality of our numerical ^^ methods since we do not know a priori in which part of the parameter space the analytical results approximate the exact result to a desired accuracy. We expect the deviation to become CTi smaller as Cr further decreases. In fact, on the edge of our parameter space where Cr = 10-4, the analytical result is reproduced with deviations well below 1%. k »

Let us now turn to the WIMP cases. For a large enough decay width, an equilibrated scalar can decay while being in thermal equilibrium. This regime is in fact only realised in a o relatively small corner of the parameter space, but one point which is a good example for this behaviour is (Chp, Cr) = (104,10-L5) — corresponding to (y, A) = (4.72 x 10-8, 8.3 x 10-5) for our standard reference values of mS = 1TeV and g* = 106.5 — as depicted in figure 5. Glancing at the distribution function (left panel) first, the evolution appears to be relatively similar to the FIMP example just discussed — although the distribution looks slightly flatter

at early times. This distribution also seems to be slightly colder than a thermal one, as to i_^

be expected [31, 38, 39, 47, 48]. However, also this point will turn out to correspond to hot DM.

The evolution of the abundances reveals the difference to the FIMP case, cf. right panel of figure 5. Here, it is clearly visible that the scalar (solid black line), although starting with a vanishing initial abundance, equilibrates very quickly and then follows the thermal curve (dotted gray line).8 During this time, the scalar is highly relativistic and thus its decay is in fact not very efficient. However, given that its abundance is thermal and thus very large, the few occasional decays are sufficient to gradually built up a sizeable abundance of sterile neutrinos.9 The scalar remains in equilibrium for a relatively long time, until r ~ 5, and if it was stable it would just resemble the generic behaviour for frozen-out WIMPs (cf. dotted gray curve). However, given that the scalar decays relatively quickly, the frozen-out abundance does not survive and is converted into sterile neutrinos (solid red line).

We can again compare the abundances for late times, which in this case gives nN(r = 250)/HS(Cr = 0, r = 250) ~ 327.6. This is vastly different from the previous result, but this behaviour is easy to understand: as long as the scalar is in equilibrium, it may decay more or less arbitrarily fast without its abundance being affected, because it is constantly regenerated by the thermal plasma. This happens very efficiently, so that a very large number

8We have explicitly checked that this happens independently of the initial abundance, as it should.

9In fact, one could equally well interpret this case simply as the sterile neutrino itself being a FIMP.

Particle number densities for: CHP=10' Cr=10"L5

0.001 0.010 0.100 1 10 100

10-4 0.001 0.010 0.100 1 10 100

r=ms/T

Figure 5. The same as figure 4, but for the WIMP case with decay in equilibrium.

R 0.01

o.ooi o.oio o.ioo 1

/iC n s(Cr=0)

Particle number

densities for:

CHP=104 Cr=10"4 n N/ ^ therm nS

10-4 0.001 0.010 0.100 1

r=ms/T

10 100

Figure 6. The same as figure 4, but for the intermediate WIMP case.

of sterile neutrinos is produced while the scalar is still in thermal equilibrium. Of course, for the frozen-out abundance alone, the factor of two would again be present — but this part makes up only about 0.6% of the final sterile neutrino abundance. Hence, this case indeed corresponds very well to the limit of only having scalar decays in equilibrium. Using the approximation obtained in eq. (4.4), we indeed obtain a final abundance which agrees with the numerical result within 4.5%.

When we go down to smaller values of Cr, we will reach the intermediate WIMP regime where some scalars decay in equilibrium while another non-negligible fraction decays only after freeze-out, thereby contributing to the final sterile neutrino abundance in similar amounts. The example displayed in figure 6 is (CHP,Cr) = (104,10-4), which corresponds to the bot-

tom right corner of the parameter space considered, and to (y, A) = (2.7 x 10-9, 8.3 x 10-5) for the reference values for mS and g*. Having a look at the distribution function first, see left panel, one can see that a twin peak structure is visible — just as seen in the example for the free-streaming horizon failing, cf. figure 2. Looking at the evolution with the time parameter r, it is visible that for early times the left (lower momentum) peak gradually builds up, while the right (higher momentum) peaks is only generated much later. This behaviour already suggests that the left peak arises from decays in equilibrium while the right peak is generated by late decays with the scalar S already having frozen out and thus being out of equilibrium. This notion is supported by the late peak being the one corresponding to higher momenta: the freeze-out of the scalar happens at a temperature close to its mass, while the decay is a gradual process which takes some time to happen after the scalar has become non-relativistic. Thus, some energy is "stored" in the scalar mass and, once the ^^ scalars decay, the characteristic energy scale of the resulting sterile neutrinos will be of the order of the scalar mass — which may be considerably larger than the temperature of the Universe at that time.

Glancing at the right panel we can see a matching behaviour in the abundances. While f the scalar is in equilibrium, it gradually builds up a sizeable abundance of sterile neutrinos. However, this process ceases to be efficient once the scalar turns non-relativistic, thereby dropping in abundance and ultimately freezing out. Although the scalar abundance is much ^—^ smaller than it was during the equilibrium time, now that the particles are non-relativistic the decays become efficient and completely translate the abundance of frozen-out scalars into sterile neutrinos, where the particle number is again doubled (but only for the frozen-out part).

Finally, we have in section 4.2 also discussed the case where the scalar practically fully decays out of equilibrium. On the other hand, given that the case just discussed was already located at the very edge of the parameter space, this final case does not seem to make sense at all. However, this does not mean that the decay solely out of equilibrium does not ^^ exist, but only that it is not very accurate to treat it under the assumption g* ~ const., as we will demonstrate in section 6.5. We still want to briefly discuss a toy example of this case, for parameter values of CHP = 104 and Cr = 10-5, which is depicted in figure 7 and which corresponds to Lagrangian level couplings of (y, A) = (2.4 x 10-10, 8.3 x 10-5) in our benchmark scenario.

Even though the double logarithmic scale in both panels might be misleading, the main contribution now comes from the late decay of frozen-out scalars. This can be seen easily when considering the sterile neutrino abundance n« in the right panel. While the scalar is in equilibrium, a small abundance of sterile neutrinos builds up until a plateau is reached for r ~ 10. Subsequently, the decay sets in and the relic abundance of scalars is converted into sterile neutrinos at high momenta. The final abundance exceeds the value of the intermediate plateau by more than a factor of 20. Also, the average momentum (x) is strongly dominated by the high momenta, just as expected. Since the production during equilibrium is proportional to CHP, we would have to lower this parameter by several orders of magnitude to make the first peak vanish in a log-log plot. As already argued, such a case would be far beyond the validity of our assumption of a constant g* during production and will hence not be considered in this work.

Up to now, we have mainly been concerned with the DM abundance, but we have not yet shown whether the DM produced is in accordance with structure formation. A discussion of this point will be given in the next subsection.

s~ 0.01

Particle number densities for: CHP=104 Cr=10"5

0.001 0.010 0.100 1 10 100 1000

0.100 r=ms/T

Figure 7. The same as figure 4, but for the Out-of-equilibrium WIMP case.

6.2 Results for the free-streaming horizon

As discussed in section 5.1 we present the free-streaming scale 1) as calculated by our numerical approach, fully taking into account the evolution of a(t) (cf. appendix B), and 2) following the estimates put forward in [39, eq. (20)]. In figure 8, we again display the bands in the CHP-Cr plane reproducing the correct relic abundance (just as in figures 1 and 3), but this time colour-coding whether the sterile neutrino is hot (red), warm (purple), or cold (blue) according to the definitions in section 5.1. The scalar mass is chosen to be mS = 1TeV, however, the results depend only mildly on the mass of the scalar [39, 40, 47]. This is true since in our computation the only dependence of the scalar mass enters through the effective number of d.o.f. which are a function of physical temperature. The strongest physical dependence on the mass of the scalar is still implicit in the definition of the parameters Chp and Cr. If we lowered the scalar mass to some hundred GeV, the result would be altered only by a few percent.

In figure 8 it is clearly visible that our numerical results disfavour more of the parameter space than the analytical estimates do. In particular, everything but the FIMP case is under tension in that analysis. However, we want to emphasise once more that these results should be interpreted with care. First of all, the numerical approach also suffers from (mainly systematical) uncertainties: the simplified truncation of the time-temperature relation into two distinct regimes (either purely radiation dominated or purely matter dominated) will differ from the exact time-temperature relation, see appendix B for details. Second, as we will see in section 6.5, in the region where the sterile neutrinos are particularly hot (i.e., for small decay constants), the approximation of a constant number of d.o.f. during production becomes less reliable, which also affects the calculation of AFS. Third, as discussed in section 5.1, the free-streaming horizon — only taking into account average properties of the spectrum — cannot capture all features of structure formation and can hence only serve as an indication. More detailed analyses can be done using the so-called transfer function, i.e. the square root of the linear power spectrum of matter perturbations, which can in turn be constrained by data from the Lyman-a measurements. Such studies have been performed in [31] for sterile

2 O M U1

(a) Numerical results for the free-streaming horizon.

O H1 (J1

0.1 1 10 102 Higgs portal CHp

(b) Same as figure 8a but using the analytical estimate.

Figure 8. Comparison of numerical results and analytical estimates of the free-streaming horizon, i.e., with and without taking into account the evolution of a(t) in the computation of the integral for Afs • Note that the strong dependence of these results on mS is hidden in the definition of CHP and Cr. However, there is a residual weak dependence on mS since the time-temperature relation used to calculate the free-streaming horizon is sensitive to the absolute temperature scale (as opposed to CHP and Cr). Hence we explicitly state the benchmark value of mS = 1 TeV even though the result will not change dramatically for scalar masses varying by a factor of a few. For this reason and for the sake of being able to compare to the other plots more easily, we also renounce labelling the axes with Lagrangian level couplings.

neutrino dark matter produced from the DW-mechanism, from the SF-mechanism, and from a simplified version of scalar decay, using the Boltzmann solver CLASS [69] to obtain the transfer functions from an extended Press-Schechter approach. A similar study taking into account the subtleties of sterile neutrino DM from scalar decay as discussed in this paper is subject of on-going work [50]. First indications of this analysis seem to confirm the conclusions drawn from figure 8a rather than those of figure 8b, which puts the scenario of freezing-out scalars under severe tension and might also indicate a comparatively large lower mass limit of the sterile neutrino of roughly 20 keV for the case of freeze-in, which might be an interesting finding in particular in the context of the tentative 3.5-keV line [19, 20]. However, in order to give a definitive answer, we will extend the current study [49, 50], as already discussed in the introduction.

6.3 The dark radiation bound

10Note that this is slightly different to the standard case found in the literature, where the dark radiation component is typically highly relativistic, see e.g. [71, 72], whereas in our case we do not know a priori whether this is the case and have to take this subtlety into account by subtracting the rest energy which is negligibly small in the highly relativistic limit. Alternatively, one can estimate the contribution of non-thermal DM to dark radiation by using the Lorentz factor [73].

Somewhat related to the bound from structure formation is the bound on the effective number of neutrinos, Neff. In the SM, this number is equal to is 3.046, the small deviation from 3 ' arising due to the effects of the reheating at e+e- decoupling [70]. However, if the sterile neutrinos in our setting are too hot, they effectively act as radiation and could in that case ^^ also contribute to the deviation ANeff of Neff from its standard value.

We can calculate the contribution of the sterile neutrinos to ANeff by comparing the ^—^ kinetic part of their energy density to the energy density pili-m of a perfectly relativistic (i.e. massless) fermionic species at the same temperature in equilibrium:10

A"e»(T) ^ ff = 76? № L Mf( mX/T )2 - 1>N ) . (6.1) )

The factor of 2 in the denominator of the first term is due to the fact that our distribution function already contains both particle and antiparticle while Neff is constructed in a way to reproduce the number of families, i.e. 3 (up to the aforementioned small corrections) and not

a value of 6. Note also, that the factor (T/Tv)4 accounts for the fact that, once the reaction |_1

e+e- o yy freezes out, the photons get reheated while the neutrinos have already decoupled from the plasma, such that no energy is transferred to neutrinos by the annihilation of the electron-positron pairs. This is also true for sterile neutrinos, however, the temperature of their distribution (if one can define this quantity at all, given that the distribution may be highly non-thermal) is implicitly contained in the quantity /n(x,T). Yet, the relative reheating of the photons compared to sterile neutrinos is the same as that compared to active neutrinos, at least if the small corrections due to weak interactions are neglected. Thus we can simply use the factor (T/Tv)4 in eq. (6.1), where Tv is indeed the temperature of the active neutrinos. In this case, the value of (T/Tv)4 rises from unity to (11/4)4/3 « 3.85. Hence we include the factor of 3.85 for late times after electron-positron-annihilation while we drop it for temperatures above about 1 MeV.

Using eq. (6.1), we can — for every point (CHP,Cr) — simply fix the sterile neutrino mass such as to reproduce the correct relic abundance and then calculate ANeff at any given

temperature T. Again, all the information lies in the distribution function f«(x,T), where x = p/T: the higher the distribution peaks at large momenta, the bigger the contribution to the extra radiation will be. If, on the other hand, the sterile neutrino abundance was tiny, this would also be reflected in fN(x, T) and the integral in eq. (6.1) would yield a vanishingly small result.

In general, we have information on ANej from two different epochs in the history of the Universe: during big bang nucleosynthesis (BBN), at TBBN = 4MeV,11 the formation rate of nuclei depends on the overall expansion rate of the Universe which, in turn, depends on its overall radiation content [57, 77, 78]. Thus, if we do not want to spoil BBN, we have to respect an upper bound on the amount of extra radiation at BBN. While the Particle Data Group still cites a relatively old bound, ANeBBN < 1.5@95% C.L. [79, 80], newer versions exist: ANeBBN < 1@95%C.L. [81], ANbbn < 0.93@95%C.L. [82], and ANeBBN < 0.85@95% C.L. [83]. We have in our plots adopted the most stringent limit, to illustrate that not even the strongest constraint does influence our results in a significant way. A seemingly more stringent constraint can be obtained from the measurement of the cosmic microwave

nThe beginning of BBN happens at a temperature of a few MeV [74], and we take 4MeV as example which is known to "reset" conditions to how they were prior to BBN [75, 76]. Given that the sterile neutrinos keep cooling down until the end of BBN and even further, taking such an early temperature corresponds to some extend to a pretty conservative limit. However, we also have to take into account that the main temperature dependence is factored out of ANeff per definition, so that TBBN = 4 MeV is in fact not much more conservative than taking a value of 1 MeV, or similar.

background (CMB), which decouples at TCMB ~ 0.26 eV [84], since the CMB spectrum also depends on the expansion rate of the Universe and thus on the radiation content [85]. The bound for that time is ANe^MB < 0.32@95% C.L. [2]. (7^

In our analysis, we take into account both bounds and calculate ANeff (TBBN) as well as ^—^ ANeff(TCMB). Although the BBN bound on ANeff appears to be weaker than the one from the CMB measurement, one has to take into account that BBN happens much earlier than the CMB decoupling. Thus, given that the sterile neutrinos produced from scalar decays cool down as time goes by, it may very well be that their contribution to the radiation content at BBN is much more significant than later on (and, indeed, this will turn out to be the case here). Such settings with a type of dark radiation that contributes differently at BBN time than later on are known in other contexts, too (see, e.g., refs. [86-89]). Some example contour lines for ANej are displayed in the Chp-Ct plane in figure 9, both for BBN

(upper panel) and CMB (lower panel). As a guide for the eye, we have again displayed the |_\

lines of correct abundance, cf. figure 3, since we fix the mass of the sterile neutrino entering the computation of ANeff via eq. (6.1) by the constraint of reproducing the observed relic abundance. As can be seen, at BBN there could be a non-negligible contribution of the sterile neutrinos to ANeff (TBBN), and there is even a small excluded region which would violate the bound (marked by the red areas at the bottom centre of the plot). This region of a too large contribution arises only for very small decay widths Cr, i.e., the scalars must be very long-lived and inject highly energetic sterile neutrinos into the Universe at relatively late times. However, we would in any case exclude this region from any serious consideration because it would, trivially, fall far into the region where the DM is hot in any case, cf. figure 8. At CMB, on the other hand, there is not even a serious constraint left since the sterile neutrinos have cooled down by then and only have comparatively small momenta.

Thus, even though there is in principle a contribution of the sterile neutrinos to ANeff, no strong constraint arises from it and there is no threat to our production mechanism.

Higgs portal

(a) Deviation of the effective number of neutrinos from the SM value at BBN and resulting excluded region (red patches).

Higgs portal

(b) Same as figure 9a at the temperature of CMB decoupling. Figure 9. Deviation of the effective number of neutrinos from its SM value of 3.046.

O H1 in

6.4 Other bounds and constraints

In this section, we will discuss the two remaining conditions which may affect the DM production mechanism proposed. It turns out that both of them are no actual problems: the first one would only affect regions so far away from the interesting part of the parameter

space that they do not play a role in practice. The second problem is more of a "theoretical" nature, i.e., while it may be important to take into account, it can be easily circumvented in concrete settings and may only appear to be problematic if the Lagrangian presented in eq. (2.1) was viewed as "theory of everything" valid up to the Planck scale. However, for completeness we would like to at least briefly mention these points.

In general, collider bounds could also restrict the parameter space of our model, since after all the Higgs portal coupling A could be used to produce two singlet scalars from two SM-like Higgses. Using the limits on the mixing between a scalar singlet and the Higgs boson as proposed in [90], the bounds on A are far above even the largest value of CHP we show in our plots. This holds true even if the VEV of the scalar singlet is larger than its mass by two orders of magnitude. This behaviour can be intuitively understood by taking into consideration that even a value of CHP = 104 corresponds to rather small couplings A, due to the large factor M0/mS involved in its definition. Thus, in practice, we do not need to be bothered by any current bounds from colliders. However, at least in principle, there is an upper bound on CHP which may become important if our study was extended to considerably h^

larger values of the Higgs portal coupling. This conclusion still holds when confronted with updated analyses [91].

There is an orthogonal problem which is related to the symmetry breaking resulting from the scalar potential in eq. (2.2). As we had already explained, the shape of this potential is ^—^ determined by an underlying discrete Z4 symmetry. While we only use this symmetry for simplicity and we could even skip it without too drastic consequences, at least in principle it would be broken by a VEV (S} of the scalar field. Since the potential has two perfectly degenerate minima, different parts of the Universe could then have their ground states in different minima, thereby leading to so-called domain walls [92]. The existence of such objects would have considerably changed the history of the Universe and they are hence problematic. However, there are many possible solutions to this problem, since the slightest difference in energies of the two vacua could be generated by all kinds of physics, in which case the walls would decay exponentially quickly [93-97]. This can also happen in the case at hand [39]. I_^

6.5 Assessing the validity of approximating the relativistic d.o.f. as constant

As stated at the very beginning, we followed the assumption that g*, g*S are constant during the DM production process. This assumption impacts on the form of the spectrum, cf. eq. (3.5), and thereby also on the implications for structure formation. Even in the analytical estimate of the free-streaming horizon, it makes a difference if the dilution factor £1/3 is calculated from a starting point of g* = 106.75 or at some lower value. In order to assess this approximation, we have for every point in the CHP-Cr plane computed the temperature when the production of sterile neutrinos is completed (i.e., when the abundances surpasses 95% of its maximum value) and we plot the comparison to the SM number of d.o.f. of the primordial plasma at that temperature. To this end, we again assume a scalar mass of 1 TeV.

In our case, the extension of the SM will never contribute more than 1 + 2 x 7/8 = 22/8 = 2.75 (i.e., less than 3% of the SM value of 106.75) units to the count of d.o.f. since this would be the maximum possible contribution of scalars and sterile neutrinos both being present with a thermal abundance and relativistic velocities. Since we expect the SM d.o.f. to be much larger during the time of production, we can interpret the further d.o.f. as small additional perturbation in the background of the SM.

Figure 10. Deviation of the number of d.o.f. at the time of DM production from the maximum value of g* = 106.75. This deviation can be interpreted as a check of the goodness of the approximation of g* = const. during production, which was used to simplify eq. (3.5).

Hence the approximation of constant d.o.f. can be assessed by checking the evolution of the SM background. If the number of SM d.o.f. at end of production is still close to the maximum value of g* = 106.75, the approximation can be seen as adequate. Figure 10 shows the deviation of the number of d.o.f. at time of production from the maximum value of 106.75. Comparing to figures 8, in most of the interesting parameter space our approximation is not too bad. In fact, for the cooler regions in the FIMP case, the approximation is even excellent. Only for the relatively hot regions the estimate of the deviation compared to the exact treatment is larger than 10%, but this region is in any case not the favoured one. Since, after all, the Boltzmann treatment of DM production in any case cannot be expected to yield sub-percent accuracy [98], this means that the approximation is in fact not too bad.

O H1 (J1

7 Conclusions & outlook

We have presented a fully comprehensive study of the production of keV sterile neutrino Dark Matter in the early Universe by singlet scalar decays. The current paper lays the foundation for several follow-up considerations. Aiming at a clear overview first, we have applied some approximations that enabled us to present analytical results in addition to a detailed numerical computation, which we used to further back up certain limiting cases. Based on these initial considerations, we have derived the system of Boltzmann equations to be solved at the level of momentum distribution functions, and we have furthermore introduced a very efficient parametrisation to do so. After presenting our analytical results and discussing some aspects related to cosmological structure formation, we turned to present our numerical results. Our numerical solutions for the distribution functions have not only provided a comprehensive picture of where the observed Dark Matter abundance can be obtained in the large parameter space investigated, but they have also allowed us to account for all bounds

applicable. This is the first time that such detailed results have been obtained for the production mechanism at hand, and we have shown how to fully exploit the information contained in the distributions. We have in particular found situations in which highly non-trivial distribution functions featuring more than one momentum scale can result from the simple decay mechanism presented, which could be very interesting for cosmological structure formation. In such cases, the simple minded estimate of structure formation properties using only the free-streaming horizon will fail, as we have explained by an illustrative example. While preliminary results of our on-going work (beyond what is presented in this paper) indicate that the numerical estimate of the free-streaming horizon yields a rather comprehensive picture despite the difficulties involved, we nevertheless have to postpone a definitive conclusion to the pending results.

In general, we have not only found very good agreement between our analytical and ^^ numerical results, but we have also shown that the assumptions applied (in particular the effective number of degrees of freedom being constant during Dark Matter production and the singlet scalar having a mass larger than that of the Higgs boson) are good in a significant fraction of the parameter space, although of course certain regimes exist in which the error introduced gets unacceptably large. Thus, the natural next step will be to extend our (numerical) considerations to the regimes in which the approximations applied are not CD

valid [49], which will in particular allow us to treat considerably smaller singlet scalar masses. ___

We will furthermore extend our studies to investigate the detailed implications for structure k » formation [50], which we have clearly shown to require a more advanced machinery than the free-streaming horizon. Ultimately, we aim to provide a fully comprehensive study of keV o sterile neutrino Dark Matter production by scalar decays, so that this production mechanism can be put to the acid test to determine how good an alternative to resonant production it truly is.

Acknowledgments CID

We would like to thank Ninetta Saviano for useful discussions and we are particularly grateful to Aurel Schneider for giving us many insights into the details of cosmological structure formation and for kindly providing us with preliminary results of our collaborative work. AM furthermore acknowledges partial support by the European Union FP7 ITN-INVISIBLES (Marie Curie Actions, PITN-GA-2011-289442). MT acknowledges financial support from the Studienstiftung des deutschen Volkes.

A Details on the kinetic equations

The collision term for a species $ in contact with other species via some interaction of the generic type $ + a + b + ... o a + ft + ... is given by

C[/*] = f \dPadPb ... dP^dP/ ... x (2n)V4) (p$ + pa + pb + ... - pa - p/ - ...) x |M|2 L

X [fafp .../*(1 ± fa)(1 ± fb) ... - fafb ... (1 ± fa)(1 ± //J) ... (1 ± /*)[ .

Some remarks about eq. (A.1) are in order:

1. The quantity EPx denotes the energy of particle x and is hence given by Ex = a/pX+m2.

2. The internal degrees of freedom of a species are denoted by

3. We have introduced a symbol for the invariant phase-space element:

dP' = * <A'2>

4. The plus signs apply in the case of bosons and the minus signs in the case of fermions.

5. The squared matrix element is defined following the convention in [84, chapter 5], i.e. the squared matrix element contains all relevant symmetry factors and averages over

both the initial and the final state spins. Ç |

CL+ss(?) = — --Î--^ 4A2(2n)

hh^ss(q) 2Eq J J J (2n)3 2EÇ, (2n)3 2EP (2n)3 2Ep ( ) (A.3)

Cs (J)_4^2 W dV , (E + Eq)2 - qq' cos ^ - 4mH {Eq+Eq/)/T (A 5) Chh+SS(q)=8n 2Ej (2n)3 2E,V (E, + E,)2 - qq' cos d 6 q - (A.5)

In order to arrive at eq. (A.5) by integrating out the phase spaces in p and p' in eq. (A.3), we have used the standard phase space integral:

d3p d3p' (2n)3 2EP (2n)3 2Ej

-(2n)4S<4> (q + q/ - p - pO = / ^± (^)■ («)

A.1 Kinetic equation for the scalar

In this appendix, we derive the kinematic functions F and G in their exact form. We also illustrate some pedagogical steps to ease the derivation of the relevant collision terms. ' ^

Let us start by constructing the collision terms in the variables x and r, cf. eq. (3.4), f—^ under the further simplification that we approximate all factors in eq. (A.1) arising from the final states by unity, i.e. we neglect the bosonic/fermionic nature of the particles, which is a good approximation provided that their energy is much larger than the temperature of the plasma. There are several processes that can contribute to the production of the singlet scalars in the early Universe.

We start with the collision term describing the production of a pair of scalars SS from a pair of SM-like Higgses hh: I—*

___^_4\2(2n)4

m m \3,

x 5 (Eq + Eq, - Ep - Ep) ¿(3) (q + j - j- jf) fhq(q)fhq(ql).

In eq. (A.3), the momenta q and q' (p and p') belong to the Higgs bosons in the initial state (to the scalars S in the final state). In eq. (A.3), we have explicitly inserted the squared matrix element |M|2 = 4A2.

Using an argument based on detailed balance [99], we can make the following replacement in eq. (A.3):

/WW) = /Sq(p)/SV). (A.4)

Note that this can be shown quite easily in an explicit way in the case of a Maxwell-Boltzmann approximation that we are using, exploiting only energy conservation.

Accordingly, we can also use the explicit form of a Maxwellian distribution to simplify eq. (A.3) further. Integrating out the phase space in p and p', we obtain

In our case, the centre-of-mass velocity is given by

Pcm =1 (Eq + Eq )2 - qq' cos Q - 4m2 (A7)

Ecm = 2^ (Eq + Eq/ )2 - qq' cos Q . (.)

With this it is straightforward to arrive at eq. (A.5) from where the definition of F can directly be read after changing to the variables x and r:

amax / . —q-

f a — v

F(x,r, £) = 2n J dX X2 J dcos Q

VX2 + r2 \

(VX2 + r2 + Vx2 + r2) 2 — xX cos Q — 4£2r2

(VX2 + r2 + Vx2 + r2) — xxcos Q f" |

leads to

I" (Vx2 + r2 + VX2 + r2)2 - 4£2r2 1~| amax = min 1, max - 1, --:--- . (A.9)

where £ = m2/ms (cf. section 4). The maximum allowed value for the cosine of Q comes from the trivial constraint that the argument in the square root must be non-negative. This

(Vx2 + r2 + VX2 + r2)2 - 4£2r21

For scalar masses much larger than the Higgs mass, i.e. £ ^ 1, we can simplify eq. (A.8):

F(x, r, £ < 1) = 4nKi(r), (A.10)

where K1 is the first modified Bessel function of second kind. |_\

With this result at hand, the first part of the kinetic equation for the scalar, cf. eq. (3.6), reads:

(r, x) dT di s Mo 1 A2 . r~0-, . _

s dr( )" drdrc22*ss = ms6An4exp (-^+^)F(x,r). (A.11) I

For the process of a pair of scalars annihilating into a pair of Higgs bosons, the collision term reads:

,s (q) = I 1 W d3q' d3p d3p' 4A2 ' " q = 2Eq)J (2n)3 2Eq/ (2n)3 2Ep (2n)3 2Ep/ 8n

x (2n)V4) (q + q' - p - p') /s(q)/s(q')

= _/s(q) 4A2 f d3q' f (q') /(Eq + Eq/)2 - qq'cosQ - 4m2

2Eq 16W (2n)3 2Eq/ /s(q ^ (Eq + Eq/ )2 - qq' cos Q . (.)

Again, we can directly infer the definition of G as in eq. (A.14):

t «max

G (x,r, £) = 2n J dX X2 J dcos Q

VX2 + r2 \

(VX2 + r2 + \/x2 + r2) - xXcos Q - 4£2r2 (VX2 + r2 + Vx2 + r2) 2 - xX cos Q

(A.13)

Note that the only difference between eqs. (A.13) and (A.8) is the exponential factor in the integral, stemming from the equilibrium distribution in eq. (A.3). The entity amax is again defined as in eq. (A.9).

Thus, the second part of the kinetic equation of the scalar (again cf. eq. (3.6)) reads:

d/|S^hh(r,x) _ dT dt

= -; C

/s(x,r) /d3x/s(x,r)G(x,r).

dr dr dT ~ss^22 ms VX^Tr2 64n4'

(A.14)

Note that the collision term in eq. (A.14) contains the actual distribution function of the scalar on the right-hand side, yielding an integro-differential equation for the distribution function /s we are interested in.

The term describing the decay of a scalar into a pair of sterile neutrinos is constructed as:

(q) = -

2d3p 2d3p' 1

2J J (2n)32Ep (2n)32Ep, 2

T^y EpEp'

p ■ p'

EpEp' j

(2n)V4) (q-p-p') /s(q) (A.15)

= -/s (q)r mEs,

where we have explicitly inserted the decay width r = ^i^r. This collision term can be interpreted intuitively: it is clear that the rate at which scalars are depleted due to decays must be proportional to the time-dilated decay width r, which involves an additional boost factor owing to the physical momentum q of the scalar, and that it must be proportional to the amount of scalars present at that particular momentum, i.e. to /s(q).

With this, the third part of the kinetic equation of the scalar, is then given by:

(r,x) _ dT dt S — — --C

dr dr dT ~ ms VX^+T2'

A.2 Kinetic equation for the sterile neutrino

We also want to show the most import steps to derive at the kinetic equation of the sterile neutrino. The source term looks like:

j r2 /| (x,r). (A.16)

CS^NN — 2 x

1 rr 2d3p' d3ps

2E~Jj (2n)32Ep/ (2n)3 2Ep

-(2n)4^ — Ep/) ¿^(pS-p-p>) 2 |M|2/s(ps,t),

where the matrix element for the decay reads:

M|2 = 1 y2p ■ p' = 2y2EpEp/

p ■ p'

EpEp' j

(A.17)

(A.18)

According to our conventions for symmetry factors, we average over initial and final states.

Integrating out the phase space in p', we can get rid of the spatial ¿-distribution in eq. (A.17). Doing so, we can replace p ■ p' by p ■ ps - p2. The scalar product p ■ ps is restricted by the kinematics of the process. Using -1 < cos[<(p,ps)] < 1, where <(a, b) is the angle between the 3-vectors a and b, this gives the constraint

— pmin

(A.19)

which implies the lower boundary in eq. (A.20).

Hence, using eq. (3.7), the kinetic equation of the sterile neutrino is finally given by

— 2Cr-,

yj x2 + '

/s (x,r)

(A.20)

2 O M U1

which we use as our master equation, cf. eq. (3.14).

B Some background on the numerical calculation of the free-streaming horizon

This appendix briefly summarises some technical details implemented in our numerics in order to evaluate the integral occurring in the definition of the free-streaming horizon, cf. eq. (5.2). We treat the thermal history of the Universe in a way where there are two distinct periods of interest, namely a purely radiation-dominated one which is followed by an immediate turnover into complete matter domination at (Teq,teq) = (1.48eV,6.04 x 1011 s). Note that these quantities stemming from our rather crude approximation agree fairly well with the values that can be found in the literature.

During radiation dominance, we can infer the time-temperature relation from eq. (3.3) using the evolution of the d.o.f. as given in [61]. A relation between the scale factor and the temperature can be established by solving for T in the equation of conservation of comoving entropy density,

^g*sT3a3 = const. = so , (B.1)

normalising today's scale factor to unity. In our numerics, we implement the fitted evolution of the d.o.f. as presented in [61].

During matter dominance, integrating the Friedmann equation gives a relation between ^_^

the cosmological time and the scale factor which reads

t = teq + 2Mpi(Pm)-1/2(a3/2 - <2) . (B.2)

This can be converted into a time-temperature relation by virtue of eq. (B.1).

We have checked that our numerical treatment of the time-temperature relation reproduces other known benchmark points (like the time of CMB-decoupling) to a very reasonable ——^ accuracy. __^

References

[1] Planck collaboration, P.A.R. Ade et al., Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571 (2014) A16 [arXiv:1303.5076] [inSPIRE].

[2] Planck collaboration, P.A.R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, arXiv:1502.01589 [inSPIRE].

[3] XENON100 collaboration, E. Aprile et al., Dark matter results from 225 live days of XENONIOO data, Phys. Rev. Lett. 109 (2012) 181301 [arXiv:1207.5988] [inSPIRE].

[4] LUX collaboration, D.S. Akerib et al., First results from the LUX dark matter experiment at the Sanford Underground Research Facility, Phys. Rev. Lett. 112 (2014) 091303 [arXiv:1310.8214] [inSPIRE].

[5] SuperCDMS collaboration, R. Agnese et al., Search for low-mass weakly interacting massive particles with SuperCDMS, Phys. Rev. Lett. 112 (2014) 241302 [arXiv:1402.7137] [inSPIRE].

[6] CRESST-II collaboration, G. Angloher et al., Results on low mass WIMPs using an upgraded CRESST-II detector, Eur. Phys. J. C 74 (2014) 3184 [arXiv:1407.3146] [inSPIRE].

[7] H. Baer, K.-Y. Choi, J.E. Kim and L. Roszkowski, Dark matter production in the early universe: beyond the thermal WIMP paradigm, Phys. Rept. 555 (2015) 1 [arXiv:1407.0017] [inSPIRE].

[11 [12 [13

[19 [20

[21 [22 [23 [24 [25

[26 [27

L.D. Duffy and K. van Bibber, Axions as dark matter particles, New J. Phys. 11 (2009) 105008 [arXiv:0904.3346] [inSPIRE].

F.D. Steffen, Gravitino dark matter and cosmological constraints, JCAP 09 (2006) 001 [hep-ph/0605306] [inSPIRE].

K.-Y. Choi, J.E. Kim and L. Roszkowski, Review of axino dark matter, J. Korean Phys. Soc. 63 (2013) 1685 [arXiv:1307.3330] [inSPIRE].

E.W. Kolb, D.J.H. Chung and A. Riotto, WIMPzillas!, hep-ph/9810361 [inSPIRE].

K.N. Abazajian et al., Light sterile neutrinos: a 'white paper, arXiv:1204.5379 [inSPIRE].

A. Merle, keV neutrino model building, Int. J. Mod. Phys. D 22 (2013) 1330020 [arXiv:1302.2625] [inSPIRE].

S. Dodelson and L.M. Widrow, Sterile-neutrinos as dark 'matter, Phys. Rev. Lett. 72 (1994) 17 [hep-ph/9303287] [inSPIRE].

S. Colombi, S. Dodelson and L.M. Widrow, Large scale structure tests of warm dark matter,

Astrophys. J. 458 (1996) 1 [astro-ph/9505029] [inSPIRE]. h^

A. Boyarsky, J. Lesgourgues, O. Ruchayskiy and M. Viel, Lyman-a constraints on warm and ( )

on warm-plus-cold dark matter models, JCAP 05 (2009) 012 [arXiv:0812.0010] [inSPIRE].

L. Canetti, M. Drewes, T. Frossard and M. Shaposhnikov, Dark matter, baryogenesis and neutrino oscillations from right handed neutrinos, Phys. Rev. D 87 (2013) 093006 [arXiv:1208.4607] [inSPIRE].

A. Merle and V. Niro, Influence of a keV sterile neutrino on neutrinoless double beta decay: ^

how things changed in recent years, Phys. Rev. D 88 (2013) 113004 [arXiv:1302.2032]

[inSPIRE]. I—

E. Bulbul et al., Detection of an unidentified emission line in the stacked X-ray spectrum of OH galaxy clusters, Astrophys. J. 789 (2014) 13 [arXiv:1402.2301] [inSPIRE]. «___,

A. Boyarsky, O. Ruchayskiy, D. Iakubovskyi and J. Franse, Unidentified line in X-ray spectra of the Andromeda, galaxy and Perseus galaxy cluster, Phys. Rev. Lett. 113 (2014) 251301 [arXiv:1402.4119] [inSPIRE]. |—

S. Riemer-Sorensen, Questioning a 3.5keV dark matter emission line, arXiv:1405.7943 | I

[inSPIRE].

M.E. Anderson, E. Churazov and J.N. Bregman, Non-detection of X-ray emission from sterile neutrinos in stacked galaxy spectra, arXiv:1408.4115 [inSPIRE].

A. Boyarsky, J. Franse, D. Iakubovskyi and O. Ruchayskiy, Checking the dark matter origin of 3.53keV line with the Milky Way center, arXiv:1408.2503 [inSPIRE].

T.E. Jeltema and S. Profumo, Discovery of a 3.5keV line in the galactic center and a critical look at the origin of the line across astronomical targets, arXiv:1408.1699 [inSPIRE].

A. Boyarsky, J. Franse, D. Iakubovskyi and O. Ruchayskiy, Comment on the paper "Dark matter searches going bananas: the contribution of potassium (and chlorine) to the 3.5 keV line" by T. Jeltema and S. Profumo, arXiv:1408.4388 [inSPIRE].

E. Bulbul et al., Comment on "Dark matter searches going bananas: the contribution of potassium (and chlorine) to the 3.5keV line", arXiv:1409.4143 [inSPIRE].

D. Malyshev, A. Neronov and D. Eckert, Constraints on 3.55keV line emission from stacked observations of dwarf spheroidal galaxies, Phys. Rev. D 90 (2014) 103506 [arXiv:1408.3531] [inSPIRE].

X.-D. Shi and G.M. Fuller, A new dark matter candidate: nonthermal sterile neutrinos, Phys. Rev. Lett. 82 (1999) 2832 [astro-ph/9810076] [inSPIRE].

[30 [31 [32

[33 [34 [35 [36 [37 [38 [39 [40 [41 [42 [43 [44 [45 [46 [47 [48 [49

A. Adulpravitchai and M.A. Schmidt, A fresh look at keV sterile neutrino dark matter from frozen-in scalars, JHEP 01 (2015) 006 [arXiv:1409.4330] [inSPIRE].

K.N. Abazajian, Resonantly produced 7keV sterile neutrino dark matter models and the properties of Milky Way satellites, Phys. Rev. Lett. 112 (2014) 161303 [arXiv:1403.0954] [inSPIRE].

L.A. Popa, A. Caramete and D. Tonoiu, Constrains on dark matter sterile neutrino resonant production in the light of Planck, arXiv:1501.06355 [inSPIRE].

A. Merle and A. Schneider, Production of sterile neutrino dark matter and the 3.5keV line, arXiv:1409.6311 [inSPIRE].

F. Bezrukov, H. Hettmansperger and M. Lindner, keV sterile neutrino dark matter in gauge extensions of the standard 'model, Phys. Rev. D 81 (2010) 085032 [arXiv:0912.4415] [inSPIRE].

M. Nemevsek, G. Senjanovic and Y. Zhang, Warm dark matter in low scale left-right theory, JCAP 07 (2012) 006 [arXiv:1205.0844] [inSPIRE].

S.F. King and A. Merle, Warm dark matter from keVins, JCAP 08 (2012) 016

[arXiv:1205.0551] [inSPIRE]. P

M. Shaposhnikov and I. Tkachev, The vMSM, inflation and dark matter, Phys. Lett. B 639 (2006) 414 [hep-ph/0604236] [inSPIRE].

F. Bezrukov and D. Gorbunov, Light inflaton hunter's guide, JHEP 05 (2010) 010 CD

[arXiv:0912.0390] [inSPIRE].

A. Kusenko, Sterile neutrinos, dark matter and the pulsar velocities in models with a Higgs svnglet, Phys. Rev. Lett. 97 (2006) 241301 [hep-ph/0609081] [inSPIRE].

K. Petraki and A. Kusenko, Dark-matter sterile neutrinos in models with a gauge singlet in the Higgs sector, Phys. Rev. D 77 (2008) 065014 [arXiv:0711.4646] [inSPIRE]. I—

A. Merle, V. Niro and D. Schmidt, New production mechanism for keV sterile neutrino dark matter by decays of frozen-in scalars, JCAP 03 (2014) 028 [arXiv:1306.3996] [inSPIRE].

Z. Kang, FImP miracle of sterile neutrino dark matter by scale invariance, arXiv:1411.2773 [inSPIRE]. |—

M. Frigerio and C.E. Yaguna, Sterile neutrino dark matter and low scale leptogenesis from a charged scalar, Eur. Phys. J. C 75 (2015) 31 [arXiv:1409.0659] [inSPIRE].

L. Lello and D. Boyanovsky, Cosmological implications of light sterile neutrinos produced after the QCD phase transition, Phys. Rev. D 91 (2015) 063502 [arXiv:1411.2690] [inSPIRE].

A. Abada, G. Arcadi and M. Lucente, Dark matter in the minimal inverse seesaw mechanism, JCAP 10 (2014) 001 [arXiv:1406.6556] [inSPIRE].

D. Boyanovsky, Clustering properties of a sterile neutrino dark matter candidate, Phys. Rev. D 78 (2008) 103505 [arXiv:0807.0646] [inSPIRE].

B. Shuve and I. Yavin, Dark matter progenitor: light vector boson decay into sterile neutrinos, Phys. Rev. D 89 (2014) 113004 [arXiv:1403.2727] [inSPIRE].

K. Petraki, Small-scale structure formation properties of chilled sterile neutrinos as dark matter, Phys. Rev. D 77 (2008) 105004 [arXiv:0801.3470] [inSPIRE].

F. Bezrukov and D. Gorbunov, On the applicability of approximations used in calculation of spectrum of dark matter particles produced in particle decays, arXiv:1412.1341 [inSPIRE].

A. Merle, A. Schneider and M. Totzauer, A fully numerical investigation of keV sterile neutrino dark matter produced from singlet scalar decays, work in progress.

[50 [51 [52 [53 [54

[55 [56

[57 [58

[59 [60 [61 [62 [63

[65 [66 [67 [68 [69

A. Merle, A. Schneider and M. Totzauer, Structure formation properties of keV sterile neutrino dark matter produced from singlet scalar decays, work in progress.

J. McDonald, Thermally generated gauge singlet scalars as selfinteracting dark matter, Phys. Rev. Lett. 88 (2002) 091304 [hep-ph/0106249] [inSPIRE].

L.J. Hall, K. Jedamzik, J. March-Russell and S.M. West, Freeze-in production of FIMP dark matter, JHEP 03 (2010) 080 [arXiv:0911.1120] [inSPIRE].

A.A. Klypin, A.V. Kravtsov, O. Valenzuela and F. Prada, Where are the missing galactic satellites?, Astrophys. J. 522 (1999) 82 [astro-ph/9901240] [inSPIRE].

M. Boylan-Kolchin, J.S. Bullock and M. Kaplinghat, The Milky Way's bright satellites as an apparent failure of LCDM, Mon. Not. Roy. Astron. Soc. 422 (2012) 1203 [arXiv:1111.2048] [inSPIRE].

J.F. Navarro, C.S. Frenk and S.D.M. White, A universal density profile from hierarchical clustering, Astrophys. J. 490 (1997) 493 [astro-ph/9611107] [inSPIRE].

I. Ferrero, M.G. Abadi, J.F. Navarro, L.V. Sales and S. Gurovich, The dark matter halos of dwarf galaxies: a challenge for the LCDM paradigm?, Mon. Not. Roy. Astron. Soc. 425 (2012) 2817 [arXiv:1111.6609] [inSPIRE]. CD

Particle Data Group collaboration, K.A. Olive et al., Review of particle physics, O^

Chin. Phys. C 38 (2014) 090001 [inSPIRE].

K. Abazajian, E.R. Switzer, S. Dodelson, K. Heitmann and S. Habib, The nonlinear cosmological matter power spectrum with massive neutrinos: the halo model, Phys. Rev. D 71 (2005) 043507 [astro-ph/0411552] [inSPIRE]. O

R. de Putter et al., New neutrino mass bounds from Sloan Digital Sky Survey III data release 8 | I photometric luminous galaxies, Astrophys. J. 761 (2012) 12 [arXiv:1201.1909] [inSPIRE].

M. Carmeli, J.G. Hartnett and F.J. Oliveira, The cosmic time in terms of the redshift, Found. Phys. Lett. 19 (2006) 277 [gr-qc/0506079] [inSPIRE].

O. Wantz and E.P.S. Shellard, Axion cosmology revisited, Phys. Rev. D 82 (2010) 123508

[arXiv:0910.1066] [inSPIRE]. |_\

P. Colin, V. Avila-Reese and O. Valenzuela, Substructure and halo density profiles in a warm dark matter cosmology, Astrophys. J. 542 (2000) 622 [astro-ph/0004115] [inSPIRE].

W.B. Lin, D.H. Huang, X. Zhang and R.H. Brandenberger, Nonthermal production of WIMPs and the subgalactic structure of the universe, Phys. Rev. Lett. 86 (2001) 954 [astro-ph/0009003] [inSPIRE].

M. Viel, J. Lesgourgues, M.G. Haehnelt, S. Matarrese and A. Riotto, Constraining warm dark matter candidates including sterile neutrinos and light gravitinos with WMAP and the Lyman-a forest, Phys. Rev. D 71 (2005) 063534 [astro-ph/0501562] [inSPIRE].

S. Das and K. Sigurdson, Cosmological limits on hidden sector dark matter, Phys. Rev. D 85 (2012) 063510 [arXiv:1012.4458] [inSPIRE].

A. Schneider, Structure formation with suppressed small-scale perturbations, arXiv:1412.2133 [inSPIRE].

S. Tremaine and J.E. Gunn, Dynamical role of light neutral leptons in cosmology, Phys. Rev. Lett. 42 (1979) 407 [inSPIRE].

A. Boyarsky, O. Ruchayskiy and D. Iakubovskyi, A lower bound on the mass of dark matter particles, JCAP 03 (2009) 005 [arXiv:0808.3902] [inSPIRE].

D. Blas, J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS). Part II: Approximation schemes, JCAP 07 (2011) 034 [arXiv:1104.2933] [inSPIRE].

[70 [71

[74 [75 [76 [77 [78 [79 [80 [81 [82

[84 [85

A.D. Dolgov, Neutrinos in cosmology, Phys. Rept. 370 (2002) 333 [hep-ph/0202122] [inSPIRE].

H. Vogel and J. Redondo, Dark radiation constraints on minicharged particles in models with a hidden photon, JCAP 02 (2014) 029 [arXiv:1311.2600] [inSPIRE].

J. Hasenkamp, Dark radiation from the axino solution of the gravitino problem, Phys. Lett. B 707 (2012) 121 [arXiv:1107.4319] [inSPIRE].

D. Hooper, F.S. Queiroz and N.Y. Gnedin, Non-thermal dark matter mimicking an additional neutrino species in the early universe, Phys. Rev. D 85 (2012) 063513 [arXiv:1111.6599] [inSPIRE].

J.P. Kneller and G. Steigman, BBN for pedestrians, New J. Phys. 6 (2004) 117 [astro-ph/0406320] [inSPIRE].

M. Kawasaki, K. Kohri and N. Sugiyama, MeV scale reheating temperature and thermalization of neutrino background, Phys. Rev. D 62 (2000) 023506 [astro-ph/0002127] [inSPIRE].

S. Hannestad, What is the lowest possible reheating temperature?,

Phys. Rev. D 70 (2004) 043506 [astro-ph/0403291] [inSPIRE]. p

S. Sarkar, Big bang nucleosynthesis and physics beyond the standard model, Rept. Prog. Phys. 59 (1996) 1493 [hep-ph/9602260] [inSPIRE].

B.D. Fields, P. Molaro and S. Sarkar, Big-bang nucleosynthesis, Chin. Phys. C 38 (2014) 339 [arXiv:1412.1408] [inSPIRE].

R.H. Cyburt, B.D. Fields, K.A. Olive and E. Skillman, New BBN limits on physics beyond the standard model from 4He, Astropart. Phys. 23 (2005) 313 [astro-ph/0408033] [inSPIRE].

E. Lisi, S. Sarkar and F.L. Villante, Big bang nucleosynthesis limit on Nv, Phys. Rev. D 59 (1999) 123520 [hep-ph/9901404] [inSPIRE]. |—

G. Mangano and P.D. Serpico, A robust upper limit on Neff from BBN, circa 2011, Phys. Lett. B 701 (2011) 296 [arXiv:1103.1261] [inSPIRE].

Y.I. Izotov, T.X. Thuan and N.G. Guseva, A new determination of the primordial He abundance using the He

I X10830Â

emission line: cosmological implications, Mon. Not. Roy. Astron. Soc. 445 (2014) 778 [arXiv:1408.6953] [inSPIRE]. |_\

R. Cooke, M. Pettini, R.A. Jorgenson, M.T. Murphy and C.C. Steidel, Precision measures of the primordial abundance of deuterium, Astrophys. J. 781 (2014) 31 [arXiv:1308.3240] [inSPIRE].

E.W. Kolb and M.S. Turner, The early universe, Front. Phys. 69 (1990) 1 [inSPIRE].

M. Archidiacono, E. Giusarma, S. Hannestad and O. Mena, Cosmic dark radiation and neutrinos, Adv. High Energy Phys. 2013 (2013) 191047 [arXiv:1307.0637] [inSPIRE].

W. Fischler and J. Meyers, Dark radiation emerging after big bang nucleosynthesis?, Phys. Rev. D 83 (2011) 063520 [arXiv:1011.3501] [inSPIRE].

R. Foot, Mirror dark matter cosmology — predictions for Neg [CMB] and Neg [BBN], Phys. Lett. B 711 (2012) 238 [arXiv:1111.6366] [inSPIRE].

J.L. Menestrina and R.J. Scherrer, Dark radiation from particle decays during big bang nucleosynthesis, Phys. Rev. D 85 (2012) 047301 [arXiv:1111.0605] [inSPIRE].

P. Di Bari, S.F. King and A. Merle, Dark radiation or warm dark matter from long lived particle decays in the light of Planck, Phys. Lett. B 724 (2013) 77 [arXiv:1303.6267] [inSPIRE].

D. Lopez-Val and T. Robens, Ar and the W-boson mass in the singlet extension of the standard 'model, Phys. Rev. D 90 (2014) 114018 [arXiv:1406.1043] [inSPIRE].

T. Robens and T. Stefaniak, Status of the Higgs singlet extension of the standard model after LHC run 1, Eur. Phys. J. C 75 (2015) 104 [arXiv:1501.02234] [inSPIRE].

[92 [93 [94 [95 [96 [97 [98

Y.B. Zeldovich, I.Y. Kobzarev and L.B. Okun, Cosmological consequences of the spontaneous breakdown of discrete symmetry, Zh. Eksp. Teor. Fiz. 67 (1974) 3 [inSPIRE].

J. Preskill, S.P. Trivedi, F. Wilczek and M.B. Wise, Cosmology and broken discrete symmetry, Nucl. Phys. B 363 (1991) 207 [inSPIRE].

F. Riva, Low-scale leptogenesis and the domain wall problem in models with discrete flavor symmetries, Phys. Lett. B 690 (2010) 443 [arXiv:1004.1177] [inSPIRE].

G.R. Dvali and G. Senjanovic, Is there a domain wall problem?, Phys. Rev. Lett. 74 (1995) 5178 [hep-ph/9501387] [inSPIRE].

G.R. Dvali, A. Melfo and G. Senjanovic, Nonrestoration of spontaneously broken P and CP at high temperature, Phys. Rev. D 54 (1996) 7857 [hep-ph/9601376] [inSPIRE].

S.E. Larsson, S. Sarkar and P.L. White, Evading the cosmological domain wall problem, Phys. Rev. D 55 (1997) 5129 [hep-ph/9608319] [inSPIRE].

K. Hamaguchi, T. Moroi and K. Mukaida, Boltzmann equation for non-equilibrium particles and its application to non-thermal dark matter production, JHEP 01 (2012) 083 [arXiv:1111.4594] [inSPIRE].

P. Gondolo and G. Gelmini, Cosmic abundances of stable particles: improved analysis, Nucl. Phys. B 360 (1991) 145 [inSPIRE].

U1 —