CC BY
0lopscience
¡opscience.iop.org
Home Search Collections Journals About Contact us My lOPscience
keV sterile neutrino dark matter from singlet scalar decays: the most general case
This content has been downloaded from lOPscience. 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: 80.82.78.170
This content was downloaded on 03/01/2017 at 10:01 Please note that terms and conditions apply.
You may also be interested in:
keV sterile neutrino dark matter from singlet scalar decays: basic concepts and subtle features Alexander Merle and Maximilian Totzauer
Dodelson-Widrow production of sterile neutrino Dark Matter with non-trivial initial abundance Alexander Merle, Aurel Schneider and Maximilian Totzauer
New production mechanism for keV sterile neutrino Dark Matter by decays of frozen-in scalars Alexander Merle, Viviana Niro and Daniel Schmidt
Dark matter in the minimal inverse seesaw mechanism Asmaa Abada, Giorgio Arcadi and Michele Lucente
Lepton number violation as a key to low-scale leptogenesis Asmaa Abada, Giorgio Arcadi, Valerie Domcke et al.
Thermal and non-thermal production of dark matter via Z-portal(s) Xiaoyong Chu, Yann Mambrini, Jérémie Quevillon et al.
Sterile neutrinos with secret interactions—lasting friendship with cosmology Xiaoyong Chu, Basudeb Dasgupta and Joachim Kopp
Gamma ray tests of Minimal Dark Matter Marco Cirelli, Thomas Hambye, Paolo Panci et al.
Production regimes for Self-Interacting Dark Matter Nicolás Bernal, Xiaoyong Chu, Camilo Garcia-Cely et al.
/ournal of Cosmology and Astroparticle Physics
An IOP and SISSA journal
keV sterile neutrino dark matter
from singlet scalar decays: the most j
general case A
Johannes König, Alexander Merle and Maximilian Totzauer
Max-Planck-Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, 80805 München, Germany
E-mail: jkoenig@mpp.mpg.de, amerle@mpp.mpg.de, totzauer@mpp.mpg.de
Received September 14, 2016 Accepted October 31, 2016
Published November 21, 2016 ____,
Abstract. We investigate the early Universe production of sterile neutrino Dark Matter by the decays of singlet scalars. All previous studies applied simplifying assumptions and/or
studied the process only on the level of number densities, which makes it impossible to QQ give statements about cosmic structure formation. We overcome these issues by dropping all simplifying assumptions (except for one we showed earlier to work perfectly) and by computing the full course of Dark Matter production on the level of non-thermal momentum distribution functions. We are thus in the position to study a broad range of aspects of the resulting settings and apply a broad set of bounds in a reliable manner. We have a particular focus on how to incorporate bounds from structure formation on the level of the linear power spectrum, since the simplistic estimate using the free-streaming horizon clearly fails for highly non-thermal distributions. Our work comprises the most detailed and comprehensive study of sterile neutrino Dark Matter production by scalar decays presented so far.
Keywords: dark matter theory, Lyman alpha forest, neutrino theory, power spectrum ArXiv ePrint: 1609.01289
j Article funded by SCOAP3. Content from this work may be used
A f .v, rdV. term" of the Creative Commons Attribution 3.0 License. doi:10.1088/1475-7516/2016/11/038
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
2 The basic idea: qualitative discussion
3 Technicalities
3.1 The Boltzmann equation and its solution
3.1.1 How to crack a coupled system of non-linear partial integro-differential equations in two variables
3.2 Relevant bounds
4 Results
4.1 Very light scalars: mS < mh/2
4.2 Light scalars: mh/2 < mS < mh
4.3 Heavy scalars: mh < mS
5 Conclusions & outlook
A Details on the Boltzmann equation
A.1 Collision terms
A.2 Transformation of variables
B The freeze-out of the Higgs boson
C Failure of the free-streaming horizon as measure for non-thermal Dark Matter
D Robustness of the halfmode analysis
29 ISJ
35 CTi
•—"
1 Introduction
Invoking something invisible that even our most sensitive detectors cannot detect does not sound like science. Yet, in contemporary cosmology, we have so much indirect evidence for Dark Matter (DM) that hardly any scientist doubts its existence. DM is a non-luminous form of matter, i.e., it does not interact with light — unlike everyday objects around us. Nevertheless, having entered an era of precision cosmology, we have been able to determine that DM outweighs ordinary matter by a factor of about five in the energy balance of the Universe [1]. Further observational evidence such as galaxy [2] or cluster dynamics [3] and the Bullet Cluster [4] allow us to constrain the properties of DM: we are sure that it is a form of matter (and not, e.g., modified gravity [5]), and our best guess for its identity is a new elementary particle that is electrically neutral and massive [6]. Importantly, it must have been produced in the early Universe in the right amounts and with a momentum spectrum suitable not to spoil the formation of cosmic structures.
It is generally accepted that DM was responsible for the emergence of structures in the Universe [7], as seen in elaborate N-body simulations on high performance computers [8]. Historically, the most generic DM candidates were WIMPs (Weakly Interacting Massive Particles), which appear in popular scenarios like supersymmetry. Such particles typically form cold DM (CDM), i.e., particles produced by thermal freeze-out [9, 10] which have non-relativistic velocities. The velocity of the particles strongly impacts structure formation. For example, hot DM (HDM), which is highly relativistic, is excluded as it would have wiped out all small structures in the Universe [11, 12]. But CDM may have problems, too: it possibly forms "too many" small halos, which might even have the "wrong" structure (known small scale issues include the missing satellite problem [13, 14], the abundance of isolated small halos [15], the too-big-to-fail problem [16, 17], and the cusp-core problem [18, 19]). While |
these may be cured once baryons are correctly included in the simulations (see, e.g., [20-22]), an alternative attempt is to consider non-cold DM settings. In the literature, these ideas have triggered the slightly unfortunate terminology of warm DM (WDM), i.e., DM particles with a thermal spectrum (i.e., Bose-Einstein or Fermi-Dirac) of a temperature close to their mass. This is also assumed in many astrophysical studies, from Lyman-a bounds [23-26] over galaxy formation [27-29] to N-body simulations [30-32]. However, looking at realistic settings, many DM-models feature a non-thermal spectrum, which one cannot associate any temperature with.
Among the most popular non-cold DM candidates is a sterile neutrino with a mass of a few keV, see ref. [33] for a recent collection of information on this topic. Sterile neutrinos are motivated from a particle theory point of view, as they are related to light (active) neutrinos
and possibly even involved in their mass generation, see e.g. refs. [34-51] for concrete models. Furthermore, the topic has gained some attention in recent years, because a sterile neutrino N would — with a very small rate — decay like N ^ vy, thereby producing a nearly monoenergetic X-ray photon. A detection of the corresponding line signal has been claimed by two groups in 2014 [52, 53], but it was very actively disputed — see ref. [33] for a detailed discussion and ref. [54] for the (still ambiguous) data that the Hitomi satellite could take before enduring its unfortunate fate. It may also be worth mentioning that there could be a "dip" in the cluster data [55], tending to make a non-observation of the line in stacked data more consistent.
Sterile neutrinos with a sufficiently large lifetime are excellent DM candidates provided that 1) an efficient production mechanism exists in the early Universe which 2) produces a spectrum in accordance with cosmic structure formation. The most generic idea is to produce sterile neutrinos by their small admixtures to active neutrinos, by non-resonant active-sterile transitions, a mechanism first proposed by Langacker [56] and related to DM by Dodelson and Widrow (DW) [57]. In modern terms, one would refer to these sterile neutrinos as FIMPs (Feebly Interacting Massive Particles [58]), which do not thermalise but are instead produced gradually in the early Universe via freeze-in by their feeble couplings to Standard Model (SM) particles. While this mechanism is rather simple, it is excluded by data, due to the particles being too close to the HDM limit [26, 59]; still, it is an unavoidable addendum to the spectrum produced by any other mechanism, which can modify the DM distribution function at least for sterile neutrino masses below 3keV [60].
A popular alternative is based on a resonant enhancement of the active-sterile neutrino transitions by the presence of a primordial lepton number asymmetry. First proposed by Enqvist and collaborators [61] and related to DM by Shi and Fuller (SF) [62], this mechanism has been studied actively [63-69], and it does indeed yield a spectrum colder than that
produced by DW. However, while there is no perfect agreement between the results of different groups, at least for the results that are publicly available [63, 67-69] it seems unclear whether they are in full agreement with cosmic structure formation [70-72].
On the other hand, one could produce sterile neutrinos thermally via freeze-out, if they had non-trivial charges beyond the SM, as long as the resulting overabundance is diluted by a sufficiently efficient production of additional entropy [73-75]; this, however, is constrained by big bang nucleosynthesis [76].
A completely different direction is to produce sterile neutrinos by the decays of other particles. A generic possibility is a decaying singlet scalar S, which can easily couple to sterile neutrinos after having been produced itself in the first place. Apart from this scalar being the inflaton [77-79], in which case it had been present all along in the early Universe, one can tune the Higgs portal coupling of S in such a way that it is either produced like a WIMP via freeze-out [80-82] or like a FIMP via freeze-in [83, 84]. Variants have been presented in [85-95], some also featuring other parent particles such as vectors [96-98], Dirac fermions [99], or pions [100, 101]; related aspects such as influences from inflation [102, 103]
or thermal corrections [104] have been discussed, too.
Our main goal is to close a gap in the treatment of sterile neutrino production from scalar decays. In earlier works [80, 81, 83, 84], simplifying assumptions have been applied to at all arrive at a result, such as neglecting active-sterile mixing, taking the number of relativistic degrees of freedom to be constant, and assuming a heavy scalar. This is also true for a previous paper by two of us (AM & MT) [84], which was the first to show how to numerically compute momentum distribution functions of sterile neutrinos from scalar decay. While it has been proven that neglecting active-sterile mixing is in fact a very good assumption [60], |—* going to the low-mass region of the singlet scalar — where the number of degrees of freedom is not constant — poses new technical challenges. The only treatment available for this regime was put forward in [88]. However, the authors only used rate equations, such that an actual computation of the momentum distribution function and the confrontation with structure formation data cannot be reliably performed. We will in this work present the full numerical computation on the level of momentum distribution functions instead, including a detailed treatment of all subtleties related to the many new technical aspects that arise in this regime. We will give an a-posteriori justification of some of the assumptions made in [88] (such as the Higgs staying equilibrated during S-production), while we will also reveal that others (like the Higgs degrees of freedom in the unbroken phase) are maybe less justified. We will furthermore show in great detail how to constrain the resulting spectra by cosmic structure formation.
Our study presents the first complete treatment of scalar decay production in the whole parameter space. It involves a fully general solution of the production equations on the level of distribution functions. Apart from our techniques being transferable to virtually any type of decay production, our study will guide the particle physics community towards using limits from cosmic structure formation without having to leave their comfort zone completely. We thus contribute to closing the gap between particle physics models, early DM-production, and their phenomenological consequences.
This text is structured as follows. We first present a qualitative discussion of decay production of sterile neutrinos in section 2, before explaining how to solve the evolution equations and to apply the relevant bounds in section 3. Our main results are presented in section 4, which features a thorough discussion of all relevant aspects from the production process over the distribution functions to structure formation. We conclude in section 5.
Technical aspects on the Boltzmann equation, on the evolution of the Higgs distribution, on the failure of the free-streaming horizon as a reliable estimator, and on the robustness of our half-mode analysis are discussed in appendices A, B, C, and D, respectively.
2 The basic idea: qualitative discussion
This section is intended to introduce the particle physics model used in our analysis and to give a qualitative understanding of its different regimes. Furthermore, we will justify the assumptions that go into the numerical computations.
Our setting introduces one real scalar singlet S (with mass mS) and one right-handed neutrino N (with mass mN) beyond the particle content of the SM.1 The right-handed r-neutrino is coupled to the singlet scalar via a Yukawa-type interaction with coupling y,
CD- 2 SN^N + h.c., (2.1)
while the new scalar singlet is coupled to the Higgs doublet $ via the most generic potential symmetric under a global z4-symmetry:2'3
Vscaiar = 2 mS S2 + ^ S4 + 2A ($t$)s2. (2.2)
Accordingly, the complete Lagrangian of the model reads
where Cv is the part of the Lagrangian that can give mass to the active neutrinos. We do not |—^ assume any vacuum expectation value (VEV) for S in our analysis, however, we will include the constraints which would arise of the VEV of S gave the sterile neutrinos their mass, so that the case of (S} = 0 can be qualitatively recovered from our analysis.
Note that we assume active-sterile mixing to be so small that the contribution from the DW mechanism to the production of sterile neutrinos is negligible compared to the part produced by scalar decay, based on taking into account X-ray limits on the decay of sterile neutrino DM. In ref. [60] it has been shown that the DW modification to any previously produced population of sterile neutrinos (from whichever main production mechanism) is of a few percent at most, and completely negligible for masses mN larger than about 4keV if one demands that the mixing angle does not exceed the limits set by X-ray analyses such as [106].
Our setup mainly allows us to produce scalars from their coupling to the SM with interactions of the type XSMXSM ^ SS. Subsequently, the scalars decay into sterile neutrinos via S ^ NN. We have to distinguish three different regimes:
C = CSM +
2N/N + S)(rrS) - 2SNcN + h.c.
^scalar + Cv , (2.3)
xIn general, any number of right-handed neutrinos can be assumed in this model. In the case of more than one generation, the scalar will then decay into all kinematically accessible right-handed states Ni with branching fractions given by y2/ ^k yk. If the mixing between the different right-handed states is large enough, all right-handed neutrinos will decay into the lightest state (Ni by convention) quickly and the results of this paper stay unaltered under the substitution y2 ^ ^k y^. If, however, the mixing inside the sterile sector is small, there might be additional complications due to late injection of highly energetic DM particles by the decay of the heavier states into the lighter ones.
2 Suitable charge assignments are given by S ^ —S and N ^ ±iN. For further comments on this assumption, see [83, 84] and references therein. The (rather mild) consequences of giving up this simplification are discussed in [81].
3In our analysis, we neglect the effect of strong self interactions, which would arise for rather large values
of AS. For a detailed discussion of the related phenomenology, see ref. [105].
^ production channels
\ / w+ ,S Z S ,S h
/\ )>---< '
II h "s W- Z f S h
\ y w+ s z
>-K x >■<
>-h-< h—<
III h S W- S Z S f
Table 1. Relevant production channels in regimes I—III.
I Production before the electroweak phase transition (EWPT): all four degrees of freedom of the SU(2)¿-doublet Higgs $ contribute equally to the production/depletion of scalars from/into the thermal bath.
II Production after EWPT with mS > mh/2, where mh is the mass of the Higgs after electroweak symmetry breaking. Now the Higgs and the massive gauge bosons interact with the scalar in different ways.
III Production after EWPT with mS < mh/2. This is similar to case II, the difference being that the Higgses present in the thermal plasma are now kinematically allowed to decay into pairs of scalars.
The channels corresponding to regimes I-III are listed in table 1. Note that we have omitted the diagrams with a scalar S in the t- or the u-channels, since they will contribute to the production only at order O (A3), i.e., suppressed by another small factor of A compared to the leading order. In regimes II and III, the fermionic initial state f only receives a contribution from the top quark in practice, since all other channels are suppressed by their small Yukawa couplings. Couplings to lighter fermions may appear relevant at first sight, since in that case the Higgs could potentially be on-shell, which may enhance the cross section by orders of magnitude. However, this case of an on-shell Higgs needs to be subtracted adequately to avoid double-counting of decay events of thermal Higgses [107], and they ultimately do not contribute much. Note that, at precision level, even thermal corrections to the parameters can play a role [104].
In general, production will always start in regime I and then, unless finished in that regime, proceed either in regime II or in regime III after EWPT, depending on the value of mS. In the case of small Higgs portal couplings A, the scalar will freeze in. This production mechanism is most efficient when T ~ ms. In the case of larger A, the scalar first equilibrates and then freezes out. In this scenario, the freeze-out time will also be linked to ms.
^ Iewpt
freeze-in (FIMP) freeze-out (WIMP) I
III r II
mh!2 io2 ms [GeV]
Figure 1. Temperature ranges of the plasma temperature T = TY at which the scalar singlet S is efficiently produced in the freeze-in scenario (red) and in the freeze-out scenario (light blue). We present four different scalar masses (30, 60, 65, and 500 GeV) that span the variety of possible production times. TEwpt indicates the temperature when electroweak phase transition happens, while TPT is supposed to give a rough idea of the temperature when the Higgs and gauge bosons become strongly Boltzmann-suppressed, even if still equilibrated with the remaining SM degrees of freedom.
Accordingly, there is a finite span in cosmic time (or in the temperature T, equivalently) in which the production is effectively taking place.
This time span is illustrated in figure 1 for four different masses mS. Red arrows correspond to small Higgs portal couplings A, implying that the scalar never equilibrates, but freezes in instead. In this case, we have defined the time span of production starting when 10% of the final yield Y of a would-be stable scalar is produced and ending when 90% are produced. Blue arrows indicate the time spans for scalars freezing out before decaying into sterile neutrinos. In this case, there is no physically preferred initial time. Note that, even when assuming a vanishing initial abundance, thermalisation is fast compared to the time scale of the decay. The late-end boundary (i.e., at low temperatures) of the time span is, however, defined by Y/Yeq = 10, where Yeq denotes the equilibrium yield.
One can clearly see that, in cases where mS < mh/2, the freeze-in of scalars occurs only after EWPT. This can be explained by the fact that the decay of thermal Higgses after EWPT completely dominates the production as soon as this channel is kinematically accessible [88]. For mS = 65 GeV, production sets in before EWPT and stops at temperatures of some tens of GeV. In the case of a much heavier scalar (say, mS = 500 GeV), production even ends before EWPT, since the kinetic energies of the Higgses after EWPT are insufficient to produce the heavy scalar.
Note that we need the particle distribution function of all sourcing species in the plasma, like gauge bosons or Higgses, to account for the production of scalars induced by the channels shown in table 1. In general, the dynamics of each of these species is determined by their own Boltzmann equation. Fortunately, it is not necessary to incorporate more than two equations
into the system, since we can safely assume these particles to be in thermal equilibrium, an assumption that has been readily made in [88] even for small T, but not proven to work. We have explicitly checked that the Higgs indeed stays in equilibrium at least until T « 1 GeV. Note that the usual estimate for WIMP-like particles of mass m is that freeze-out occurs at Tfreeze-out ~ m/20. The reason that the Higgs stays in equilibrium longer are inverse decay processes, and similar arguments apply to the gauge bosons and to the top-quark, cf. appendix B. Still, the number densities of all SM particles become strongly Boltzmann-suppressed for temperatures of, say, 10-20 GeV. Hence, production always ceases at these temperatures, simply because there are hardly any particle left in the bath, even if they are still in equilibrium.
3 Technicalities
3.1 The Boltzmann equation and its solution ^^
In order to compute the distribution functions of the scalar and the sterile neutrino, we employ the Boltzmann equation, L[f] = C[f], where f is the distribution function we want
to determine, L = — Hp dp (with the Hubble function H) is the Liouville operator in a i_^
homogeneous and isotropic Friedman-Robertson-Walker universe, and C contains the collision
terms describing interactions of the particles.4 In our case, we actually need to solve two —^ coupled Boltzmann equations, for the two distribution functions fN and fS of the sterile
neutrino N and the scalar S, respectively: ^—^
L fs = CS (3.1)
and LfN = CN. (3.2)
The collision terms are themselves sums of several individual terms which encode the details of the respective processes as well as information on whether they contribute to production
or depletion of the species under consideration. Depending on the regime, various diagrams (jiJ listed in table 1 can contribute to the production of scalars and hence, ultimately sterile neutrinos. Explicitly, the scalar collision term Cs consists of the following:
Regime I (T > Tewpt):
— C<MnSS + CJUNN . (3-3)
• Regime II (T < TEWPT and mS > mh/2):
CIS — + + CW+W--n-SS + + CS^NN . (3-4)
• Regime III (T < TEWPT and mS < mh/2):
CIII — Cfi + Ch-nSS . (3-5)
The labels should be quite intuitive. For example, in regime I, the term encodes both
the annihilation of two Higgs doublets $ into two scalars S and vice versa, with opposite signs
4This is already an approximation in itself because, by using the Boltzmann equation — which is a classical
equation — we do neglect quantum corrections. To include those, Kadanoff-Baym equations must be used, which are way more difficult to solve. The error caused by the approximation we make here is of 0(10%) in the abundance [108].
for the two cases. On the other hand, CS^NN will always come with a negative sign, since it only describes the decay of the scalar into sterile neutrinos. Similarly, the only difference between regimes II and III is the appearance of the Higgs decays into two scalars, C^SS, for sufficiently small scalar masses.
The sterile neutrino collision term is comparatively simple, and it is given by:
CN = cN^NN , (3.6)
but in this case with a positive sign as it creates sterile neutrinos. Here we neglected the inverse decay, i.e. the process NN ^ S, because the sterile neutrino is a FIMP, such that its abundance is always far below its would-be equilibrium abundance for large enough temper-
atures T. Hence the process NN ^ S is suppressed by the small number of sterile neutrinos present in the Universe.
Note that the collision terms will in general depend on the distribution functions f of all particle species i. However, as argued, we only need to compute the distributions fs and fN, since all SM species can be assumed to be in thermal equilibrium such that their distributions are well approximated by simple Boltzmann factors.5 To anticipate one particular example collision term, we quote from eq. (A.23) the one describing the production of sterile neutrinos of momentum p at temperature T by the decaying scalars,
rN \f l(p T) y2mS / dp/ p fS(p, T) (37)
Cn-nnfs](p,T) = J dp , (3.7)
where the lower boundary of the integral over scalar momenta p/ is given by p'min N = |p — mf |. This term is indeed positive, since it produces sterile neutrinos. As expected, this equation contains the distribution function fs of the scalar at temperature T, making it explicit that eqs. (3.1) and (3.2) are coupled and have to be solved as a system of equations. For a detailed discussion of all collision terms we refer the reader to appendix A.1.
3.1.1 How to crack a coupled system of non-linear partial integro-differential equations in two variables
Glancing at the explicit forms of the collision terms in appendix A.1, we can see that the task of computing sterile neutrino production from scalar decay is a highly non-trivial one: it necessitates the solution of a coupled system of partial integro-differential equations in two variables. In an abstract, yet intuitive manner, the system of equations we have to solve
5 As already anticipated, this is a non-trivial assumption for the SM Higgs boson. We have shown explicitly that the Higgs only freezes out at relatively small temperatures, T ~ 3 GeV, which is due to its decay into SM particles keeping it in equilibrium for a longer time than naively expected for a WIMP-like particle. This also provides an a posteriori justification of the assumption made by Adulpravitchai and Schmidt in ref. [88]. See appendix B for details.
turns out to have the following form:6
Lfs(p, T) = D(p, T)fs(p, T) + J dp'KS(p,p', T) [fs(p, T(p', T) - fSq(p, Tfq(p', T)] ,
L/n(p,T) = Jdp'KN(p,p',T)fs(p',T), (3.8)
where and D are known functions directly related to the collision terms C, and feq
denotes the (possibly hypothetical, in case the interactions are too feeble) equilibrium distribution of particle i.
The question is now how to tackle such a problem. Before starting the computation we note that, while the sterile neutrino distribution function fN does depend on the distribution
6Note that, at this point, we could also have taken an alternative route. The part of the integral that only contains Boltzmann distributions does not depend on fs and could hence, in principle, be integrated numerically. This would result in a numerical function of the form F(p, T) on the right-hand side of the first eq. (3.8). However, this strategy would have at least two drawbacks: 1) first, the structure of the equation would be much less clear, since several crucial dependencies would just be "hidden" inside a numerical function F(p,T). 2) Furthermore, this procedure would ultimately force us to subtract two numerically computed integrals from each other, which is numerically less accurate than first computing the difference of the integrand functions before evaluating the integral over their difference. We thus stick to the strategy outlined in the main text.
7We exploit the one-to-one correspondence between cosmic time t and plasma temperature T.
function fs of the scalar, the converse is not true. The reason is that, as we had argued, we can neglect any inverse processes involving sterile neutrinos. This already yields to a considerable simplification, allowing us to solve the first eq. (3.8) independently and to then insert the result obtained into the second equation. However, even then, we are still left with |— a non-linear partial integro-differential equation in fs which is highly non-trivial to solve numerically (and impossible to solve analytically). We will thus need to play some tricks, in order to tame this equation.
The first step is to perform a transformation of variables, (t,p) ^ (r,£), such that C )
the left-hand side of the first eq. (3.8) only involves a derivative with respect to a single |_i
variable. As shown in appendix A.2, this is possible for r = f (t) being an arbitrary function f independent of the scalar momentum p, and £(p, t) = g(p) simultaneously being an arbitrary function g of a specific combination of time t and momentum p, which involves the ^ " scale factor a(t) and an arbitrary reference time to. Thus, in fact, there exists a whole family of functions f and g which would simplify our complicated equation, however, choosing them in a smart manner may even lead to further simplifications. In our view, the following choice of variables is particularly convenient:7 00
r = T ,
£ = 1 a(t) p = (gsM Y/3 p (3 9)
{ To fl(t(To))p V g-(T)J t , (3.9)
where gs(T) is the number of effective entropy degrees of freedom (for which we have used the numerical fit developed in ref. [109]). The arbitrary reference mass m0 and temperature T0 have been chosen by us to both equal the Higgs mass, m0 = T0 = mh, which will later prevent
us from working with extremely small or large numbers in our numerical computation. This choice of variables indeed implies certain simplifications, e.g., the Liouville operator now reads:
L= H = -<")(! + lf|, (3.10)
where / denotes a derivative with respect to the temperature T.
A valid interpretation of the new variables is to view r as "time" variable and £ as "momentum" variable. Indeed, r increases as the temperature decreases, just as the cosmic time t, so it can really be thought of as a rescaled or stretched time. In turn, £ can be thought of as something like a comoving rescaled momentum. Alternatively, one could think of it as a rescaled momentum that is red- or blue-shifted with respect to the reference temperature To. C | Indeed, the big advantage of the variable £ is that the distributions remain constant in time r as soon as the production of particles is finished. For example, if the scalar S was stable, its distribution for late time would be given by fs(£,r) = fs(£,rprod), Vr > rprod, where rProd marks the time when no S-particles are produced anymore (e.g. the freeze-out time in h^ case S was a stable WIMP). This effect translates into the sterile neutrino distribution, i.e., fN will also be constant in r once the scalar decays cease to be efficient. All redshifts of the momenta are automatically included in the definition of £.
The resulting form of the Boltzmann equations used in our implementation is
f (£-r) = rm(l—31 >n\g.(r>])ci\f<.fe.]. p-11»
where i = S, N. Indeed, we have by our substitution transformed the partial differential equation in two variables into an equation containing a differential with respect to one variable only. Of course, the collision terms in their general form of eq. (3.8) have to be transformed according to this change of variables. Note that this general form can be written down for —^ both the scalar and the sterile neutrino, but it is only necessary for the former, as the sterile neutrino evolution equation is comparatively simple once fS is known. We will therefore concentrate on the case i = S in the following.
Still, the main issue remains: the collision term is nonlinear in fS, and it furthermore includes an integral over £, which renders this equation to be of a partial integro-differential type, which are very difficult to solve. While several strategies for certain types of integro-differential equations (such as the Volterra-type) are available in the literature, the equation under consideration turns out not to be of any such type, so that we have to find a dedicated workaround. We apply a trick to get at least an (arbitrarily good) approximation by discretising the equations and substituting the integral over £/ by a finite sum over M discrete momentum values £i, i £ {1,..., M}. This transforms the original non-linear integro-differential equation (3.11) into a system of coupled non-linear ordinary differential equations for the different modes fi(r):
df (r) = DDi(T)fi (r)+ £ it(r)f (r)fj (r) — fi'eq(r)fj'eq(r)] , V i £{1,...,M} .
j=l,...,M
(3.12)
Note that the integral results in a "non-local" coupling, since it couples any mode f£ not only to some "neighbouring" modes, but to all modes instead.
Finally, we have to numerically solve the resulting system of differential equations. This step involves one more difficulty, since the most simple strategies (such as the Runge-Kutta method) fail for large parts of the parameter space due to the stiffness of the system.
Another more technical issue is that packages such as Mathematica natively employ list-based algorithms, while matrix-based ones are much more appropriate for the equations under consideration. We have thus found it more convenient to use the built-in Matlab solver ode15s [110, 111], which is particularly efficient when dealing with stiff problems. Nonetheless, one should try to tweak any solver algorithm to take advantage of our a-priori knowledge that a distribution function must never attain negative values. This way, we have been able to solve the equations for the scalar distribution functions fS efficiently. The resulting numerical functions can then be inserted into eq. (3.11) for i = N or, alternatively, one can use the "master formula" as given in eq. (16) of ref. [84], once the variables are transformed according to eqs. (3.9).
Obtaining the functions fN (r, £) was the goal we wanted to reach. These distribution C | functions form the basis to compute all properties of the DM species for a given point in the parameter space, see ref. [84] for details. Not only can we compute the DM abundance QNh2, we can also use fN (r, £) to derive predictions for cosmic structure formation, which can then be matched to observations. We will explain what to do with fN in the next subsection. For
instance, the particle number density of sterile neutrinos is computed as expected,
oo 3 oo
"» (r) = S2 / di |p2(i)/N (i.r) = S2 ^ / d«2 fN (i.r), (3.13)
where gN = 2 counts the spin degrees of freedom of the sterile neutrino N. From the particle
number density, the Dark Matter abundance follows as
„ .2 .0 ........ „„, )
QDMh = SiT^T • 0 -t/h2 . (3.14)
.('prod) pcrit/h
where n(rprod) and s(rprod) are the number and entropy desnsities at r = rprod, respectively, s0 = 2891.2 cm-3 [112] is today's entropy density, and pcrit/h2 = 1.054 • 10-2 MeVcm-3 [112] is the critical density in units of the squared reduced Hubble constant h.
Note that there are publicly available codes to compute the momentum distribution functions for resonantly produced sterile neutrinos in a similar fashion, see refs. [68, 69].
3.2 Relevant bounds
By solving the set of Boltzmann equations for the scalar S and for the sterile neutrino N, we obtain the momentum distribution function of the sterile neutrino, which contains all relevant information and must therefore be in accordance with a number of (partly model-dependent) observational or experimental constraints. In the following, we will present how this is tested and which assumptions go into these tests.
Mass of the sterile neutrino. Since we neglected active-sterile mixing (as argued in section 2 and in [60]), the momentum distribution function does not yet contain any information on m^, but only on the number density of steriles present in the Universe. We can, however, fix the mass of the sterile neutrinos by demanding that they make up all (or a certain part) of the cosmic DM energy density, cf. eq. (3.14). This will later on result into the allowed regions marked in our plots.
Tremaine-Gunn bound. The mass of the sterile neutrinos is also restricted by an absolute lower limit, the so-called Tremaine-Gunn (TG) bound [113], which arises from fermions having a maximal phase space density. When applied to astrophysical objects such as the central regions of galaxies, the resulting mass turns out to be mN > 0.5keV [114]. This bound will result into a relatively large excluded region in our plots, and it will in practice separate the regions of the scalars freezing out and freezing in.
Overclosure. In principle, even for the lowest possible mass of mN = 0.5 keV, a too large sterile neutrino number density will overclose the Universe (i.e., lead to Qtot > 1). Since this is not in agreement with observations, a certain patch of the parameter space is excluded by this bound. However, the overclosure bound turns out to be irrelevant in practice, since it is inferior compared to the TG bound.
Structure formation. A simple estimator for the predictions of structure formation is the so-called free-streaming (FS) horizon Afs [25], which yields an estimate of the average length scale (usually given at redshift z = 0) that a DM particle would have travelled from the time
of production tprod until today (t0) — had it not been gravitationally trapped. It can be computed via
to |—1
Afs = f dt . (3.15)
J a(t)
tprod ^^^
A free-streaming horizon below something like 0.01 Mpc indicates a scenario that — in terms of structure formation — cannot be distinguished from cold Dark Matter (CDM), although the border is arbitrary to some extend (but reasonably well motivated to be about one order of magnitude below the border to HDM). A value of 0.1 Mpc, the typical scale of dwarf
galaxies, conventionally indicates the cross-over from models that suppress power on small „___
scales (compared to CDM) in a way that is still in accordance with observational data from HDM models that suppress too much power. As already mentioned, in the literature this
regime is mostly referred to by the very unfortunate term warm DM, which seems to indicate (jiJ a thermal spectrum. However, as we will see in the following, realistic distribution functions are in many cases (highly) non-thermal, which may be very badly reflected by such an over-simplistic quantity like the FS horizon.
It is essential to bear in mind that Afs is — at best — to be understood as an order-of-magnitude estimator. In fact, depending on the case under consideration, it might have to be modified by an O (1) factor to give more accurate classifications of models. Moreover, by its definition in eq. (3.15), it does not take into account the full spectral form of the DM distribution and will therefore never accurately treat non-thermal distributions. We have nevertheless exemplified this analysis (and its failure) for one selected mass of the scalar for the sake of a comparison, cf. appendix C.
A more robust analysis uses the linear matter power spectrum P(k), which reflects the full spectral information. For every point in the parameter space, we therefore compute the linear power spectrum using the CLASS code [115, 116] for non-cold DM species, and then normalise it to the linear power spectrum of a perfectly cold distribution. This ratio, P(k)/Pcdm (k) is also known as squared transfer function T2 (k), which gives information on how strongly the power (i.e., the number of halos formed) is suppressed compared to a perfectly cold distribution at the scale k. Another advantage is that this quantity can then be compared to limits on the squared transfer function obtained from Lyman-a data, displayed by a limiting squared transfer function T2m(k).
Limiting functions (k) are usually obtained assuming a thermal DM distribution (in that case truly warm DM). Accordingly, it is for principle reasons difficult to compare a given model to the observational limits. Ideally, one would re-analyse the Lyman-a data using the exact spectral shape of the DM distribution and the respective mass of the sterile neutrino. However, this is impossible given the virtually infinite number of possible distribution functions. Nonetheless, if T2(k) > T2m(k) for all k, we can safely categorise the model as in agreement with the Lyman-a bound applied. If, on the other hand, T2(k) < 7l2m(k) for all k, the model clearly has to be discarded.
There is a problem with this procedure, though. In some cases, T2(k) > 7l2m(k) might only hold true for a certain interval of wave numbers k, but not for the entire range. For example, T2(k) > T[2m(k) could be true for all k smaller than a particular value but not C | anymore for k > This complication originates from the very fact that the limiting transfer functions are derived from thermal spectra, while realistic spectra can be highly non-thermal, which results in a difference in the exact evolution of the squared transfer function around its cutoff (e.g., the slope may be different). h^
In order to still use the limiting squared transfer functions, which are very good estimators accounting for observational relevance, we have developed the following procedure:
1. First, we compute the half-mode k^2, i.e., the wave number at which the squared transfer function has dropped to 1/2:
k1/2 ^ T2(kV2) = 1/2 . (3.16)
2. Subsequently, we check whether T2(k) > T2m(k) is fulfilled all k < k1/2. If the condition is met within this range of small wave numbers (i.e., larger length scales), we consider the model as being consistent with a given bound. This is clearly an approximate method, since we intrinsically disregard some relevant information (namely the power spectrum below the half-mode). Furthermore, the value of 1/2 is, ultimately, pure u> choice and one could equally well justify other values. However, it is clear that the qq border has to be somewhere between the two extreme cases of a certain squared transfer functions being completely untouched by the Lyman-a bounds (restrictive view) or by at least in small parts being allowed by the bound (conservation view), which is why a value of 1/2 appears to be a fair compromise. We have in addition shown that our results are, in fact, quite robust with respect to this arbitrary choice: due to the power spectra considered still having a slope somewhat similar to that of a thermal relic, even a choice of 0.05 instead of 0.5 would not alter our results very significantly, cf. appendix D. We are currently working on more advanced analyses that will result into even more reliable procedures.
This procedure is in fact rather simple, as can be seen from the cartoon-like illustration in figure 2. For the Lyman-a bounds, we use the squared transfer functions derived from thermal spectra with thermal masses of rnum = 2.0 keV (rnum = 3.3 keV) in a conservative (more restrictive) scenario. Note that only the restrictive bound is displayed in figure 2, as the purpose of this figure is only to illustrate the principles behind the procedure, while later on we will display both bounds for completeness. The values quoted are motivated in [24], where the authors also provide an analytical fit formula for the transfer function of a given thermally distributed species of mass rnnm, which is assumed to make up all the DM.
^ 0.6 ST
g 0.4 0.2
y' I 1 ¿r ^t*»*10**T I \Allowed \
High-power region consistent with ' Lyman-a bound Half-mode: \ », \ \ V \ \ » \ t • \ \ -\ \ \ * \ \ \ \
T 2№I/2)=1/2 » \ t V I * ♦ 1 1 \ * V -\ \ 1 \ X \
Lyman-a bound (thermal relic with a mass of 3.3 keV) \ 1 \ V\ \ V v
k [h/Mpc]
^ 0.6 ST
-----_ \ \"" \ V J - \ \ \ ^ \ > \ \ \ \ . \ \ \ \ \ \ \ \ \ \ Half-mode: High-power region inconsistent with ^----- Lyman-a bound \ Forbidden w\ \l \
T 2№I/2)=1/2 K\........................... » \ \ 1 v \ i \ \
Lyman-a bound \ \ \ V \
(thermal relic with \ 4
a mass of 3.3 keV)
k [Ä/Mpc]
Figure 2. Illustration of how we classify certain squared transfer functions as consistent (left) or inconsistent (right) with a certain bound.
Bounds from the observed amount of dark radiation. The combination of the momentum distribution function of the sterile neutrino and its mass is also constrained by the observation of the cosmic microwave background (CMB) and by the light element abundances related to big bang nucleosynthesis (BBN). The corresponding observables allow to constrain the effective number Neff of neutrinos and its deviation from the SM value of 3.046 [117]. The effective number of neutrinos is a measure of the radiation present in anything else than photons at a given epoch. It is therefore a measure of the Universe's expansion rate, which in turn leaves its imprint on both the CMB and the abundances of the elements produced during BBN. A similar analysis was shown in [84, figure 9], however, with a small error in the numerical computation of the contribution to ANej from the non-cold sterile neutrinos. While this error led to an overestimation of their impact, these bounds would in any case only be relevant in the region that is already excluded by the requirement of DM not being classified as "hot". Of course the precise numbers have to be computed numerically, but it can easily be understood that the bounds from structure formation and from Neff must be closely related, since they both critically depend on how long DM stays ultra-relativistic.
Note that, apart from possibly impacting the region in the parameter space where the DM is produced extremely late, the dark radiation bound could in fact also be relevant for very small DM masses. However, as we have seen, any value of mN below 0.5 keV is already excluded by the TG bound.
Model-dependent bounds on the most minimal particle physics setting. So far, we have fixed the mass of the sterile neutrino by matching its number density to the observed DM energy density. In the most minimal setting, we could simultaneously demand that the mass of the sterile neutrino is generated by a VEV of the singlet scalar S, i.e., mN = y (S). This would allow us to use bounds from perturbative unitarity and bounds derived on the scalar mixing angle from its contribution to the W-mass. In principle, there are also direct collider bounds [118], but they are off our plots and do in practice not constrain the relevant region of the parameter space.
INJ O M CTi
Following [118], we can bound the VEV of the scalar by:
(S >>V ¿ms . (3.17)
Using (S> = ^ enables us to conclude an upper bound on the Yukawa coupling:
mN /16n , ,
y <-\ . (3.18)
ms V 3
A singlet scalar mixing with the Higgs via its VEV can also yield a radiative correction to the W-boson mass. Combining eqs. (8), (9), and (11) from ref. [118], we can deduce an upper limit on the Higgs portal coupling:
A < Amax = y sinmax(2a)|mS mhl , (3.19) Nj
2vEWmw
where vEW is the VEV of the SM-Higgs and a is related to the mixing of the scalars. |—
This completes our list of bounds which can possibly be relevant. As we will see, apart from the Lyman-a bound, the most relevant constraint arises from the TG bound. In the most minimal setting, collider-related constraints can become relevant, too, however these are strongly model-dependent and in fact very easy to circumvent. Finally, the overclosure C ) and dark radiation bounds play no role in practice.
4 Results
In this section, we present our main results for sterile neutrino Dark Matter. However, before o discussing the sterile neutrinos themselves, it is very useful to understand the evolution of the singlet scalar in the early Universe, which will make many of our results much easier to grasp.
So let us, just for a moment, assume that the scalar is stable and understand what happens to it. As we had already mentioned, the scalar can either freeze in or freeze out in the early Universe, depending on its interaction strength A. We have illustrated all the different cases in figure 3, where we display the ratio between interaction rates rint of the singlet scalar and the Universe's expansion rate H, as functions of the time parameter r = m^/T. The interaction rates rint are computed using all diagrams available in the different regimes, cf. section 2, which depending on the processes may receive contributions from different scattering or decay diagrams. As can be concluded from table 1, all diagrams are proportional to A, which is why the evolution of the ratio rint/H with r looks identical for every plot in figure 3, up to a rescaling. However, it is an important information whether rint/H < 1, in which case the scalar freezes in, or whether rint/H > 1 and falls off only later on, in which case the scalar first thermalises and then freezes out.
Freeze-in. Let us start with the freeze-in case. This regime does depend somewhat sensitively on the scalar mass, which is why we have in figure 3 depicted the two cases of A = 10-9 (left panel: for ms = 30, 60GeV) and A = 10-8 (central panel: for ms = 65, 100, 500, 1000 GeV). In all cases, we can clearly see that the scalars are always far from
Scalar freezes in
I A=10-9
Figure 3. Interaction rates rint compared to the Universe's expansion rate H, as functions of the time parameter r = mh/T. We illustrate two freeze-in cases — depending on the scalar mass we show either A = 10-9 (left) or A = 10-8 (center) — as well as one freeze-out case, with A = 10-5 ( right). Notably, in the latter, all scalars undergo a cold freeze-out a temperatures within [mS/20, mS/4], which explicitly proves that it is a good approximation to consider the scalar to be at rest at the point of freeze-out [70], contrary to the claim made in [119].
equilibrium.8 We can furthermore see the change in the interaction rates due to the new diagrams available after the EWPT. In particular for very small scalar masses, ms = 30 or 60 GeV, we can see that the interaction rate jumps up by some two orders of magnitude. This is because for such small masses the Higgs decay h ^ SS is available, which is the dominant contribution as noted in ref. [88]. For larger scalar masses, in turn, there is hardly any difference visible (if at all, the interaction rates seem to decrease a little); however, the main reason is that larger scalar masses are thermally suppressed at this late stage, as most clearly visible for the largest scalar masses of ms = 500 or 1000 GeV. But even for medium-scale masses of ms = 65 or 100 GeV, a small drop in the rate is visible. This is in parts due to kinematics, since the W-, Z-, and Higgs bosons as well as the top quark are already thermally suppressed at this stage, cf. diagrams in table 1, but it happens also because of the different structures of the diagrams involved. In any case, the freeze-in of the scalar ceases at a temperature that is roughly equal to its mass, and indeed all interaction rates start to drop shortly after (or even before) T = ms is reached. A larger coupling A translates into a larger number density, as expected.
Freeze-out. For the larger coupling, A = 10-5, scalars of all masses equilibrate. As visible in the right panel of figure 3, all interaction rates are larger than the expansion rate already for very early times. As long as the scalars are equilibrated, the actual interaction rate does not matter much. However, a boost in the rate — as due to the EWPT for ms = 30 and 60 GeV — can "delay" the freeze-out. This is particularly visible for a ms = 60 GeV, which freezes out at r « 41 (or at T « ms/20), and thus even later than the scalar with ms = 30 GeV, for which freeze-out happens at r « 30 (or T « ms/7). Thus, as we will later see to be correct, we can expect a much larger number density of scalars (and hence of sterile
8In the central plot in figure 3, one can also see the curves for ms = 30, 60 GeV in light colours, which cases apparently can thermalise after a while. This would then lead to an "overabundance", at least for this part of the parameter space, which will later also be visible in our main plots — the sterile neutrinos would need to be so light in order to compensate for the massive number density produced that they are already excluded by several bounds.
neutrinos) for ms = 30 GeV than for 60 GeV. Similarly, the case of ms = 65 GeV should yield a larger abundance than the one for ms = 100 GeV, which also turns out to be correct.
Sterile neutrinos. Let us now come to sterile neutrinos as DM. In order to include as much information as possible in the figures, without however making them too busy, we illustrate the allowed ranges and all relevant constraints as "spaghetti plots" displayed in figures 4, 7, 10. In these plots, we depict the regions in the A-y parameter space where the correct DM abundance is met (where only part of the abundance is produced) by the dark coloured solid lines (by the lightly coloured bands), for different values of the sterile neutrino mass: mN = 2, 7.1, 20, 50, 100 keV. To obtain these regions, we have for each combination of the parameters (A, y, ms) computed the resulting distribution function along the lines described in section 3, and then integrated the result to obtain the DM abundance according to eq. (3.14).
In addition, we display the TG phase space bound, as explained in section 3.2. The overclosure bound (gray area) marks the part of the parameter space where we would obtain overclosure of the Universe even for the minimum sterile neutrino mass of 0.5 keV.
We also display the model-dependent bounds, which only hold in the most minimal |—
setting, cf. section 3.2. While they may not even exist in more involved settings, they do |_
serve the purpose of representing the maximal effect collider-related bounds could possibly have. To indicate some numbers, the effect of the unitarity bound, eq. (3.18), for a scalar mass 100 GeV is that the Yukawa coupling y has to be smaller than about 8.2 ■ 10-8 [2.9■ 10-7, 8.2 ■ 10-7] for mN = 2keV [7.1 keV, 20keV], in very good agreement with figures 4, 7, 10. For a scalar mass of 1000 GeV, these bounds get stronger by a factor of roughly 10. The W-boson mass correction, eq. (3.19), is somewhat more subtle, since the scalar mixing angle a depends on the Higgs portal A which in turn influences the DM abundance and is thus related non-trivially to the Yukawa coupling y and mN. As discussed, these bounds do only »—^ exist in the most minimal model, which is why we have marked them by thin dashed lines (see areas left of the TG bound and down in the lower right corners). For some cases, they
are even off the plot. (jJ
The main bound however arises from structure formation, where we take into account the two Lyman-a bounds derived in ref. [120] and perform the half-mode analysis described in section 3.2 to determine whether a given distribution function is consistent with the data, or not. We have for each scalar mass ms numerically computed the linear power spectrum for each pair (A,y). Reproducing the correct abundance results in a condition on the sterile neutrino mass mN, so that every point with the correct abundance can be characterised by a point (ms, A, y, mN) in the parameter space. Depending on which Lyman-a bound we used, we have in our plots marked the following regions:
• forbidden: if the upper half of the squared transfer function is forbidden by both the conservative and restrictive Lyman-a bounds, the respective point is red.
• constrained: if the upper half of the squared transfer function is only forbidden by the restrictive Lyman-a bound but allowed by the conservative one, the respective point is purple.
• allowed: if the upper half of the squared transfer function is allowed by both the conservative and restrictive Lyman-a bounds, the respective point is blue.
This colour code is slightly reminiscent of the historically grown terms "hot", "warm", and "cold" DM, however, we would like to stress once more that the distributions obtained are
Figure 4. Abundance and constraints for small scalar masses, ms < m^/2. The crosses mark the points displayed explicitly in figure 5.
h-> h->
non-thermal, and one thus cannot associate any temperatures with them; see appendix C for an explicit counterexample showing that the classification drawn from thermal spectra is bound to fail for non-thermal distributions. Nevertheless it is correct to say that, roughly, the red regime is associated with rather high momenta compared to the DM mass, while the blue one tends to feature smaller momenta. By the light colours, we have indicated regions with a sizable but insufficient abundance for a given mass. Technically, these scenarios correspond to a slightly larger mass which yields the same effect as replacing a fraction of, say 10% of the Dark Matter by a perfectly cold Dark Matter component. We have checked that this simplified procedure yields correct results up to the level of accuracy limiting analyses using Boltzmann equations anyway.
A final comment on the unavoidable but tiny contribution by Dodelson-Widrow production at temperatures of a few 100 MeV: as had been shown in ref. [60], this modification does not only hardly contribute to the abundance for sterile neutrino masses of 4keV and higher, also the modification of the actual spectrum is completely negligible. Thus, in our plots, there would be no effect visible even if we did take into account the largest modification allowed by data. This may only be different for a sterile neutrino mass of 2keV, but even then the error by not considering the modification is of about 10%, which is in practice still negligible. For even smaller masses, the correction could be sizable, but in such cases a mixed scenario would be excluded by structure formation anyway.
4.1 Very light scalars: mS < mh/2
Let us begin with the case of very light scalar masses, i.e., ms < m^/2, displayed in figure 4. As discussed in section 2, this leaves us with an interval of roughly 5 GeV < ms < 62 GeV. Here, the lower limit arises from the temperature where the relevant interactions are frozen out and/or not efficient anymore, while the upper limit is obtained from m^ — 125 GeV. Here, depending on the coupling A, we may enter different regimes. First of all, for very small couplings A, we will be in the FIMP region, i.e., in the left part of the plots. Freeze-in generically is most efficient at temperatures T around the mass ms of the FIMP, i.e., in
regime III (see figure 1 and table 1). Of course the production of scalars has started in regime I, at high T, but in any case — even with all diagrams of regime III at work — the interaction remains feeble. Once the scalar is produced, it decays with a certain lifetime proportional to y-2 into sterile neutrinos. Thus, for y large, the decays happen very close to the freeze-in of S and thus the resulting DM particles will be more strongly redshifted and "cooler". For too small y, in turn, the scalars remain present in the Universe and decay only very late, thus injecting too much energy into the sterile neutrinos and rendering them HDM — and thus excluded. Notably, as can be seen in the plots, the actual mass of the scalar within the allowed interval does not make much of a difference. The decay rate is proportional to y2ms, such that a larger scalar mass ms should translate into more red-shift, both because the correct abundance is obtained earlier and and on top of that the decay happens faster. However, for a fixed y, the unboosted decay rate also becomes smaller, such that the lifetime ^^ effectively increases, leaving less time for the steriles to redshift. To a first approximation, these two effects compensate each other, but of course there will be corrections due to scalars of different masses being boosted more or less (and thus having bigger or smaller lifetmes in the cosmic frame). Apart from these minor effects, the FIMP regions for ms = 30GeV and ms = 60 GeV appear rather similar, up to a shift in the Higgs portal A.
On the other hand, for larger values of A, the scalar can enter thermal equilibrium and then freeze out, see the lower right regions of the plots. Depending on the exact combination of (A, y), the decay into sterile neutrinos happens while the scalar is in equilibrium (horizontal part of the isoabundance-lines — which is independent of A from a certain value on), after freeze-out (vertical part — which is similar to freeze-in in the sense that the scalars cool
down before their decay and thus "forget" about their cosmological history for a long enough lifetime or small enough coupling y), or in between (kink-region — where a strong dependence on both A and y is present). For this part of the plot, the question of whether there is at all any allowed region depends strongly on ms. This observation can be understood by realising that the scalar freezes out earlier and thus with a much larger abundance for ms = 30 GeV, which is due to the dependence of the interaction rates on the scalar masses. Hence, given many more scalars than for ms = 60 GeV, a too large abundance can only be avoided for sufficiently small sterile neutrino masses — which is why only one strip is visible in the lower right corner of the left plot, and this one is forbidden by structure formation due to the small DM mass favouring a "hotter" spectrum, i.e., with a stronger tendency for large momenta.
Two representative example evolutions of the yield can be seen in the top row of figure 5, both for the scalar freezing in (left) or out (right), which correspond to the crosses marked in the right figure 4. If the scalar was stable (gray dashed line), it would compare trivially to the hypothetical equilibrium curve (gray dotted line): as a FIMP (left), it would gradually be produced but never reach a thermal abundance before freezing in, while as WIMP (right) it would quickly thermalise before freezing out. Once the decay into sterile neutrinos is switched on, the true abundance of the scalar (black solid line) again increases initially but ultimately decreases to zero. Note that, in both plots, one can see the increase in the interaction rate due to the EWPT, cf. figure 3. The sterile neutrino abundance (red solid line), in turn, increases until the point where ultimately all scalars available have decayed.
Let us make a quick comparison to the results obtained in ref. [88]. Looking at the top two plots in their figure 2, which feature the same scalar masses as our top two plots in figure 5, our yields look qualitatively very similar except for a missing "bump" in the yield of the scalar for temperatures around 10GeV, which appears in ref. [88] but not in our result. While we cannot explain this feature, we do however have good agreement with the overall
Figure 5. Example evolutions of the yield (top row) and sterile neutrino distributions (bottom row) for a scalar with a mass of 60GeV undergoing freeze-in (left) or freeze-out (right) before decaying into sterile neutrinos, for the two points marked in the right figure 4.
abundance obtained in [88],9 up to some 20% difference, which can be attributed to applying different numerical procedures and also due to ref. [88] using some approximations such as rate equations assuming a (suppressed) thermal shape.
The evolution of the resulting sterile neutrino distributions for the same two points is displayed in the bottom row of figure 5, with all spectra being functions of £ = ()1/3T , cf. the second eq. (3.9). We can see a certain deformation in both curves, although the one on the left (for the scalar being a FIMP), it is probably a small bump rather than an actual double peak. In both cases, these shapes come from two different production phases. In the FIMP case (left), the two phases correspond to the production before and after the
9To see this, one has to adjust the pair (A, y) to the case from [88], though, taking into account the different normalisation of A in their eq. (5) compared to our eq. (2.3).
0.0 0,
cons. Ly-a \ res. Ly-a
- mS=60GeV 'A \
\\ •
A=10"8'83 A \
- y=10-800 \
mN=7.11keV \ \
\ \ « \ X \ \ \ %
cons. Ly-a \ \ res. Ly-a
mS=60GeV A=10"6-86 y=10-869 mN=20keV
5 10 k [h/Mpc]
0.0^ 0.5
5 10 k [h/Mpc]
Figure 6. Example squared transfer functions for a scalar with a mass of 60 GeV undergoing freeze-in (left) or freeze-out (right) before decaying into sterile neutrinos, for the two points marked in the right figure 4.
EWPT, whereas the visible bump for larger £ for the WIMP case (right) rather comes from production while the scalar is still in equilibrium and after its freeze-out, respectively.10
Finally, in figure 6, we display the squared transfer functions for the two example points. Glancing at the right figure 4, we can see that both points lie inside "purple" regions, i.e., they are constrained but not excluded by the Lyman-a bounds. According to the procedure detailed in section 3.2, this would imply that the parts of the squared transfer functions with k < k1/2 should both be consistent with the conservative Lyman-a bound, but inconsistent with the restrictive one. This is just what we can see in figure 6. It is also visible that, in particular for the FIMP case depicted on the left, the slope of the functions is different from that of the bounds, for which the distributions were thermal by construction. Furthermore, note that the sterile neutrino mass on the left plot, 7.1 keV, is significantly different from the mass of 2 keV corresponding to the conservative bound, which originates from the nonthermal shape of the scalar-decay produced DM spectrum. Obviously, according to the left figure 6, one would classify the case displayed as "warmer" than the restrictive bound corresponding to a mass of 3.3 keV, even though the true DM mass is even larger than that. This is one more reflection of conclusions which do hold for thermal spectra being invalidated once the shape of the momentum distribution function is more complicated. Thus, any "translation" of the structure formation properties to a hypothetical spectrum of thermal shape, as attempted for DW-produced sterile neutrinos on several occasions in the literature (e.g. refs. [26, 99, 120-123]), must be treated with extreme care and should not be regarded as a precision statement. The latter statement is in agreement with ref. [124], where the authors conclude that there are differences between resonantly produced sterile neutrinos and thermal models in the mildly non-linear regime at high redshifts.
We can also see that, although both of the points that we had explicitly illustrated above fall into the same category in what concerns structure formation, the FIMP case tends to be slightly "warmer", i.e., closer to being excluded than the WIMP example. This may come as a small surprise, given that the average momentum of the WIMP-produced sterile
10Note that, in the spectrum on the right, a tiny bump stemming from the changes during the EWPT is nevertheless visible at around £ ~ 0.03, thereby implying a third momentum scale. However it is so tiny that it does not have a real influence on the spectrum.
Figure 7. Abundance and constraints for intermediate scalar masses, mh/2 < mS < mh. The crosses mark the points displayed explicitly in figure 8.
neutrino is nearly two-times larger than the FIMP-produced one.11 However, when dividing the average momenta by the sterile neutrino masses, hence looking at the average velocities, the factor of nearly two remains present — just that in this case the FIMP-produced particle yields the larger value.
Looking back at figure 4, to get a more global picture, it is visible in the freeze-in (left) regions that for larger Yukawa couplings y the DM particles get "colder" (i.e., shifted to smaller momenta). This is due to the parent particles decaying earlier, thus leaving more time for the DM particles to redshift. Of course, this effect is more pronounced when the DM particles are heavier, see upper left corners of the plots. Interestingly, the unitarity bound would cut off part of that otherwise allowed parameter space (but only if mN = y (S)). For the freeze-out (right) regions, instead, the case of small Yukawa coupling corresponds to a very late decay, always far after the freeze-out. This leaves the DM particles no time to cool down, so that they will generically be threatened by structure formation. The second limit of large A keeps the scalars in equilibrium for a very long time, so that practically all sterile neutrinos are produced while the scalar is still equilibrated. That leads to a "colder" spectrum, because the DM particles have more time to redshift. The turnover, corresponding to the double peak structure observed in ref. [84] for the first time, can be somewhere in between. Depending on the (non-trivial) interplay between the different parameters and in particular on the relative strength of the two peaks, it can be ruled in or out by structure formation.
4.2 Light scalars: m^/2 < ms <
The next case to be discussed, cf. figure 7, is that of "light" scalars, in the sense that their mass is still less than the mass m^ of the Higgs boson, but larger than half the Higgs mass. This corresponds to regime II in table 1 for low temperatures T, while the production nevertheless starts at high T and thus in regime I, cf. figure 1. Here, the main characteristic is that, at the EWPT, all of a sudden lots of new channels open up. For the freeze-in case, i.e. small A,
2 O H*
11 Note that, while one has to be slightly careful when converting the quantity £ to a physical momentum, cf. section 3.1.1, it can still be used to compare the relative average momenta of two given spectra.
this means that it is easier to produce scalars at a temperature just above their mass, which will result in an increase in their abundance and thus also in that of sterile neutrinos — however this increase is much less dramatic than for the very light scalars, and hence hardly visible in the plots, due to the Higgs decay into two scalars not anymore being accessible in this mass range. For the freeze-out case, i.e. large A, this means that the scalars will be kept in equilibrium for a longer time in regime II, resulting into a large overall number of sterile neutrinos produced while the scalar is still equilibrated.
The spaghetti plots for two different scalar masses, ms = 65 and 100 GeV, are depicted in figure 7. At first sight, these plots appear qualitatively similar to those depicted in figure 4, except for an overall shift towards larger values of A. However, an interesting observation is that, for the FIMP case, the spectrum seems to become "warmer" when going from ms = 65 C | to 100 GeV, while the trend was opposite when going from ms = 30 to 60 GeV. The reason is that, for the case at hand, Higgs decay does not play a role in the production of the scalar. Instead, many scalars are already produced rather early (see upper left panel of figure 8), such that they can indeed decay earlier (and are less boosted), so that in this case a smaller scalar mass results into smaller initial sterile neutrino velocities and this effect dominates the production.
As already anticipated, the lighter scalar (ms = 65 GeV) freezes out earlier, namely at T ~ 17GeV, than the heavier one (ms = 100 GeV), which does so at T ~ 15 GeV, cf. figure 3. This may look surprising at first, since usually heavier particles tend to freeze out earlier due to the thermal suppression for non-relativistic particles. But the thermal suppression is similar for both cases, e-65/100 ~ 0.5, so that the actual interactions are decisive. Already at T ~ 50 GeV, all heavy bosons and the top are thermally suppressed, so |— that they cannot easily be produced anymore. However, a scalar with a mass of 100 GeV can even at rest annihilate into both W+W- and Z0Z0, while a lighter scalar cannot. Thus, the interaction rate of the heavier scalar is considerably larger, thereby keeping it in equilibrium significantly longer. It follows that the number density of the lighter scalar, and thus of the o sterile neutrinos originating from its decay, is higher and hence must be compensated by a smaller sterile neutrino mass — which is the origin of only small masses being allowed in the freeze-out region of the left figure 7.
Again, two example points are depicted in figure 8, with the order of the plots as in figure 5. It is clearly visible that, in the yields of the scalars, the effect of the EWPT is much less dramatic than for the very light scalars. In the resulting sterile neutrino spectra, the effect of the EWPT is thus much less pronounced, although little bumps can be spotted on the left sides of both spectra in the bottom row of the figure. Qualitatively, the course of the DM production is similar to case of very light scalars: again, if the scalar freezes in, it gradually decays into sterile neutrinos; if it freezes out, then it may decay either in or out of equilibrium, which a gradual transition between the two regimes. Note that, for the masses displayed, the comparison to ref. [88] is in fact much more difficult. While the same scalar masses are used in that reference, we strongly suspect that paper to have a small error in the degrees of freedom of the Higgs field above the EWPT. While the Higgs field does possess four physical components in the unbroken phase, in ref. [88], only one seems to have been used. In other words, at high T, while the gauge boson contributions are switched off, the contributions of the would-be Goldstone bosons are not switched on, as one would need to do. We will confirm this claim later in section 4.3, where it results in a factor of four between our results and those obtained in ref. [88]. However, this factor is only close to four for very heavy scalar masses (even when taking into account the difference in the normalisation of A,
10-4 0.001 0.010 0.100 1 f
0.2 0.5 1 2 r=mH/T 5
10-2 Distribution Function ___ —^(0-1.55 10-2
mS=65GeV /
A=10-7'82 /
10-4 y=10-700 / 10.0 10-4
mN=7.1keV / 0.91\
10-6 0.29 \ \ 1 10-6
10-8 \ \ 1 10-8
10-10 - 1 I \ |l \ 1 III 10-10
1 1 11
10-12 \ I | ||| 10-12
r=mHIT
50 100
Distribution Function mS=65GeV A=10-422 (0-15.9 —\__422
y=10-865 ^^ mN=7.1keV y^ 14.0\
0.50 N\ \ \ 1 1 t 1 1
r=0.19 " i 1 1 i 1 i 1 '. » I < I i i
10 100
10-4 0.001 0.010 0.100 1 f
10 100
Figure 8. Example evolutions of the yield (top row) and sterile neutrino distributions (bottom row) for a scalar with a mass of 65 GeV undergoing freeze-in (left) or freeze-out (right) before decaying into sterile neutrinos, for the two points marked in the left figure 7.
which also happens to be a factor of four), where the production almost entirely happens in regime I. But for the comparatively light scalars treated in this subsection, the production is spread out over regimes I and II, cf. figure 1, such that a discrepancy by less than a factor of four but no perfect agreement can be expected. This is confirmed by our numerical results, making us confident that our treatment is the correct one for high temperatures.
The squared transfer functions (i.e., the ratios of the resulting linear power spectra to the CDM case) are depicted in figure 9. As to be anticipated from the left figure 7, the FIMP point is located in a blue region, and it should thus be perfectly allowed by all bounds. The WIMP point, in contrast, is in a region that is entirely red — and it should thus correspond to DM that is by far too hot. Indeed, the squared transfer functions displayed in figure 9 confirm this conclusion. The whole curve of the FIMP case is in the allowed region, whereas the whole curve of the WIMP case is forbidden by far, even by the conservative bound. For
conS' Ly-a \ \ res. Ly-a
mS=65GeV \ \ \
A=10-7'82 \ \ \
- y=10-7'00 \ \ \ ^
mN=7'11 keV V
res. Ly-a
mS=65GeV : A=10-4'22 0.2 y=10-8'65
mN=7.11keV
5 10 k [h/Mpc]
0.0 0.5
5 10 k [h/Mpc]
Figure 9. Example squared transfer functions for a scalar with a mass of 65 GeV undergoing freeze-in (left) or freeze-out (right) before decaying into sterile neutrinos, for the two points marked in the left figure 7.
Figure 10. Abundance and constraints for large scalar masses, m^ < ms. The crosses mark the points displayed explicitly in figure 11.
such clear cases, an elaborate analysis would not even be necessary. Again we observe that, although both points feature a sterile neutrino mass of 7.1 keV, the predictions for cosmic structure formation are different by far. This clearly illustrates that the DM mass itself is not at all a good way to discriminate between two different DM spectra, even if viewed at one and the same temperature.
4.3 Heavy scalars: mh < mS
Finally, the scalars may also be heavier than the Higgs, see figure 10, which is the case that had already been discussed in our earlier reference [84]. In this case, given that both freeze-in and freeze-out are correlated to the mass of the scalar, most of the production of scalars happens above the EWPT — we basically remain in regime I in figure 1 during the whole production. Only for somewhat light scalar masses, such that ms/20 < Tewpt, a
small amount may be produced after the EWPT. However, for most of the production, only the simple 4-scalar interaction in the top row of table 1 plays a role. In particular, one can greatly simplify the equations involved for (mh/ms)2 ^ 1, see ref. [84] for details. As we can see from the central and right panels of figure 3, it is in fact a good approximation to assume that the production of sterile neutrinos from very heavy scalars happens entirely before the EWPT — thereby providing a clear justification of the assumptions made previously in [84].
Understanding the results for this case is somewhat less subtle than for the previous two. Comparing the two plots depicted in figure 10, one can see that there are no overly drastic changes. The reason is that, as argued in ref. [84], there is a trade-off between a heavier scalar and an earlier decay: a larger scalar mass means that the production of scalars ceases at higher temperatures. In the freeze-in case, this indicates a lower overall production, which is what the vertical lines in that regime slightly shift to the right when going from ms = 500 GeV to 1000 GeV. However, given that no new production channels open up while staying in regime I, the change is rather mild. In the other limiting case, where the scalar freezes out, a larger scalar mass means a larger abundance for the case where the scalar
decays mainly after freeze-out (i.e., the vertical lines on the bottom right of the plots). Thus, these vertical lines should slightly shift to the left to compensate for that, which they indeed do, as visible in the plots. For the horizontal parts of the lines, where the scalar decays while in equilibrium, however, freezing out earlier means less sterile neutrinos unless this is compensated by a larger decay rate: thus, these lines should shift to slightly larger values of the Yukawa coupling y when going from ms = 500 GeV to 1000 GeV, which is just what is visible in the plot.
Two representative example evolutions of the yield can be seen in the upper row of figure 11, both for the scalar freezing in (left) or out (right). Just as before, in the FIMP case, the scalar is gradually produced but never reaches a thermal abundance before freezing in (while constantly decaying), whereas for the WIMP-case it quickly thermalises before freezing out. Also here, the decay can happen earlier or later, depending on the value of C ) y. As already hinted in section 4.2, we expect a factor four difference compared to the results obtained in ref. [88], as DM production happens nearly completely in regime I. This approximation is best for very large scalar masses, however, the most extreme case discussed in [88] involves a scalar with ms = 500 GeV. Taking the same couplings as in the reference and taking into account the different normalisation of A, the discrepancy obtained numerically is already larger than a factor of three, and should converge to four for even larger scalar masses (which we however cannot explicitly check, given the absence of this case in [88]). Yet, given that no distribution functions had been computed in ref. [88], their results can be rectified by simply including the missing factor.
The corresponding evolutions of the distribution functions are depicted on the bottom row of figure 11, with similar characteristics as for smaller masses. One big difference, though, is that the DM production is finished much earlier, which comes from heavier scalars decaying faster. This allows more time to redshift, but it is compensated to some extend by a larger initial velocity.
In what concerns structure formation, however, these shifts do not change very much. In the freeze-out case, while the sterile neutrinos are produced by a heavier scalar and are thus faster, they are also produced earlier (since the decay rate of the scalar is proportional to its mass), so that hardly any change in the colouring can be seen for the freeze-in case when comparing the two plots in figure 10. For the case of the scalar freezing out, only the in-equilibrium decays are phenomenologically relevant, as the late decays in any case
0.05 0.10 0.50 1 r=mH/T 5
Distribution Function 1 <0=4.25
mS=500GeV
A=10"740
y=10"750 \
mN=7.1keV 0.55 \ 6.83
0.10 * \ 111 \ \ In
I \ \ 11
10"4 0.001 0.010 0.100 f
10-7 -
0.05 0.10 0.50 1 5 10 r=mHIT
Distribution Function (f\~3 03
mS=500GeV
A=10"430
y=10"809 S \
mN=20keV
'' 2.84 \ 47.2 \
\ 11 0.08 \ 1 111 4 1 I 1 1
1 I 1 1 1 1 1 1 I r=0.025 '"••. 1 111
Figure 11. Example evolutions of the yield (top row) and sterile neutrino distributions (bottom row) for a scalar with a mass of 500GeV undergoing freeze-in (left) or freeze-out (right) before decaying into sterile neutrinos, for the two points marked in the left figure 10.
produce too "hot" Dark Matter. However, for scalars decaying while in equilibrium, a larger mass translates into larger Yukawa couplings and thus into earlier decays so that, again, the resulting sterile neutrinos have a larger momentum at production but they have also more time to redshift, such that the prediction for structure formation remains basically unaffected — hence the similarity between both plots. Note that, for large scalar masses, none of the "kink" regions on the bottom right of the plot is allowed, which would correspond to a double peak structure in the sterile neutrino's momentum distribution function (due to nearly equal contributions from the decays before and after freeze-out of the scalar) — such distributions do exist but they are already peaking into the region of the restrictive Lyman-a bounds. This result is consistent with the findings from ref. [84], which analysed just that case, and it is also backed up by the example squared transfer functions, cf. figure 12.
cons. Ly-a
mS=500GeV A=1Ü-7-40 y=10-750 mN=7.11 keV
0.0 L 0.5
res. Ly-a
5 10 k [h/Mpc]
0.0 0.5
cons. Ly-a \ \ res. Ly-a
\ * \ -
- mS=500GeV \ \ -
'is 1 \ .
A=10"430 \ \ \
- y=10-809 \ \ \
mN=21keV .........\
5 10 k [h/Mpc]
Figure 12. Example squared transfer functions for a scalar with a mass of 500 GeV undergoing freeze-in (left) or freeze-out (right) before decaying into sterile neutrinos, for the two points marked in the left figure 10.
5 Conclusions & outlook
In this work, we have completed the investigation of sterile neutrino Dark Matter production from scalar decays. We have shown how to perform the full computation using Boltzmann equations on the level of momentum distribution functions. While this would in principle amount to numerically solving a coupled system of partial integro-differential equations of a rather uncommon (and possibly unknown) type, we have argued why it is allowed to decouple the equations, how to transform the partial into an ordinary integro-differential equation, and which steps to take to finally solve it using an algorithm based on a discretisation in one variable. All that comprised pioneering work, which had previously not been available without drastic assumptions or simplifications. While our findings are of course somewhat specific to the examples studied, our methods carry over to more general settings such as decays of non-scalar particles or multiple generations of sterile neutrinos: any reader inclined to investigate such scenarios can adopt our methods and will be able to apply the same strategies to crack the equations involved.
The importance of investigating sterile neutrino Dark Matter on the level of momentum distribution functions lies in these quantities carrying all the information about the Dark Matter model at hand. Not only can they be used to extract aggregate information such as the Dark Matter number or energy densities, but having the distribution functions is vital to derive predictions for and apply constraints from cosmic structure formation. Many previous works tried to go around this by solving oversimplified rate equations and using a simplistic estimate of the free-streaming horizon to conclude about whether or not a setting is consistent. However, we have shown that this procedure can fail once non-thermal distribution functions are involved. In these cases, it is just not a very reliable strategy. We have instead proposed to compute the linear power spectrum (from which we derive the squared transfer function), for which task dedicated tools are available, and to use the half-mode analysis developed in this work to obtain a fairly good estimate of the compatibility of a certain distribution function with data [125].
Having investigated the production of sterile neutrino Dark Matter from scalar decay in full generality, the final question is where to go from here. While we have shown how to obtain a reliable constraint from Lyman-a data using the half-mode analysis presented
in section 3.2, in principle, one could refine the analysis presented, as we intend to do, and/or study further halo properties of settings like the one presented, as done e.g. for mixed Dark Matter settings [126, 127]. The distribution functions obtained from scalar decay are conceptually not very different, however, as we have seen, they are highly non-thermal in both shape and number of momentum scales, so that even a description by a superposition of two thermally-shaped distributions may in fact not be sufficient. Various constraints can be applied, with maybe the most generic apart from the Lyman-a forest being dwarf satellite galaxy counts [128].12
The main limitation of our method is that, currently, we are using only the linear power spectrum [129], while more considerable differences of decay-produced Dark Matter to thermal spectra may be visible in the non-linear regime. Ideally, one would perform and |
N-body simulation for each distribution presented here, however, this is clearly not doable ^^ given that these computations are numerically very expensive. Thus, a better strategy is to first obtain a pre-selection of cases which could be potentially interesting for a full N-body simulation, while discarding those which are clearly ruled out, or in practice no different from cold Dark Matter. One way to do such an analysis is to estimate the effect of the non-linearities, by the extended Press-Schechter approach [130, 131], modified by incorporating a collapse model that describes how the free-streaming Dark Matter particles get bound
to form structures [128, 132]. This is the next step to be done, which we will leave for ___
future work [125]. Given that present-day data is already sufficient to yield vastly different k » constraints on different sterile neutrino Dark Matter production mechanisms [70] it can be expected that — with more detailed computations — we may be able to clearly distinguish
various production mechanisms by observational data.
Acknowledgments
We are grateful to M. Volpp for very useful discussions concerning the numerical solution of integro-differential equations, and to N. Menci, A. Schneider, and M. Viel for sharing their insights on cosmic structure formation. We would furthermore like to thank A. Adulprav-itchai and M.A. Schmidt for giving us detailed information on their earlier work. We also QO received valuable comments on our manuscript from M.A. Schmidt and A. Schneider. AM acknowledges partial support by the Micron Technology Foundation, Inc. AM furthermore acknowledges partial support by the European Union through the FP7 Marie Curie Actions ITN INVISIBLES (PITN-GA-2011-289442) and by the Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreements No. 690575 (InvisiblesPlus RISE) and No. 674896 (Elusives ITN). MT acknowledges support by Studienstiftung des deutschen Volkes as well as support by the IMPRS-EPP.
A Details on the Boltzmann equation
The purpose of this appendix is to discuss some technical details related to the Boltzmann equations used in the text. Apart from reporting the explicit forms of the collision terms, we will also show how to perform a convenient transformation of variables, which enables us to deal with considerably simpler equations.
12Note that this has nothing to do with the missing satellite problem. On the contrary, for non-cold Dark Matter, it could possibly happen that not enough satellites are produced, instead of producing too many. This yields quite a strong bound: at least the number of observed satellites has to be met, since satellites cannot be produced from larger structures.
A.1 Collision terms
Let us first present the explicit versions of the collision terms. Note that we will only very briefly describe the basic form of the Boltzmann equation, as information on this part can be found in many textbooks, such as refs. [133, 134]. Constructing explicit forms for the collision terms is in fact somewhat trivial, although one of course has to be careful to include all relevant factors. Example derivations can be found in the appendix of ref. [84], however, we will for illustration nevertheless present one explicit derivation, while we only list the results for all other cases.
As explained in the main text, we have two decisive structures, namely "2 ^ 2" (scattering) and "1 ^ 2" (decay) processes. The question is how to extract their forms from a general collision term. The most general form possible for a collision term describing the reaction ^ + a + b + ... o a + P + ... is [134]:
/ ] = I dPadPb . . . dPadps ... x (2n)V4)(p + pa + p6 + ... - pa - pHi - ...) x |M|2
x [// ... (1 ± fa)(1 ± fb) ... (1 ± /^) - // ../ (1 ± fa)(1 ± /) ...],
with p being the 4-momentum and Ep the energy of the particle ^ under consideration. The symbol |M|2 denotes the (initial and final) spin-averaged matrix element, which also contains ^—^ symmetry factors for identical final and initial states,13 as well as multiplicity factors for reactions involving multiple particles per process. Otherwise, the standard definitions apply: the phase space element is dPX = 2Edf(27r)3, for a particle X with internal degrees of
x |M|2/q(|p + p'|,T) - /s(p,T)/s(p',T)] . (A.3)
13This is in accordance with the conventions from ref. [134]. The reason behind symmetry factors also appearing for initial state is that, to arrive at rate equations, we have to integrate over all phase space elements, for both initial and final state particles. Thus, to avoid double counting, the symmetry factors have to be included here. This is different from, e.g., the situation at a particle collider, where the initial state is prepared in a certain way, but never integrated over.
freedom, momentum pX, and energy EX = ^JmX + pX with mass mX
To give one concrete example, we will now explicitly derive the collision term describing the decay of a SM-like Higgs into two singlet scalars S, [/s](p, T). The following
quantities label the 4- and 3-momenta, as well as the absolute value of the latter:
• P, P, P (= |p|): for either one of the scalars,
• P', p', p' (= |p'|): for the remaining scalar, 00
• q, q, q (= |q|): for the Higgs boson.
We furthermore abbreviate Ep. = ^m2 + p2. Starting from eq. (A.1), we obtain:
cLss [/s](p,T) = 2ES / p^ (¿pjhf <2n)4i<4)(p+p'- q)
X |M|2[/hq(?,T) - /s(p,T)/s(p',T)] , (A.2)
where |M|2 = 16A2v2, cf. eq. (A.21). The 3-dimensional ¿-function eliminates q:
1 t d3p'
[/s](p,T) = ^ J (2n)24Ep(E(p+p/| + E(S - Elp+P'l)
/h P P
Due to the remaining ¿-function, the term in parentheses can be cast into:
1 f d3p' r.Si o
cLss [/s ](p,T) = 2ES / (2n)24i?S; Ep+p/| ä(es + - ^p+p^
x |M|2 [/eq(p,T)/eq(p/,T) - fs(p,T)/s(p',T)] . (A.4)
Using spherical coordinates, d3p' = 2np/2dp'd(cosa), where a is the angle between p and p', we can rewrite the ¿-function to obtain:
S(ES + Ep - ßhp+p'l) = ¿(cos a - cos ao) |p+pp'"="° . (A.5)
Here, a0 is the value of a such that /(cos a0) = 0. Using these results, we find
CLss[/s](P,T) = 2Es j (^Epy / d(cOS a) ^(cOS a - COS ao)
-p J v2n) -^p'.
0 p -1
-p-p- 1 d(cos a) ¿(cos a — cos an) ^^^
'-"-' I"1
= 1 if cos a0 E [-1,1] ,
= 0 otherwise. | ^
x |M|2 /q(p,T)/eV,T) - /s(p,T)fs(p',T)] . (A.6) ^
The result of the integration d(cos a) ¿(cos a — cos ao) can be encoded in the limits of the integration over p', i.e., by restricting that integral to those values of p' for which the
cos a-integral yields one. So we need to solve |_1
(ES + ES )2 — mh — P2 — p'2 6
± 1 = cos ao =--—-—:--(A.7)
for p'. The solutions are C )
(mh — 2m2)p ± mh\/(m| — 4mS)(mSS + p2) p' =-^- for cos ao = +1, (A.8)
t — (mh — 2m^)p ± mhy/(m' — 4ms)(ms + P2)
p =-—2- for cos ao = —1. (A.9)
As can readily be seen, these amount not to four but to only two distinct values:
(mh — 2m2)p + mh\ (mh — 4m|)(m| + p2)
p1 =-2m-• (A-10)
(mh — 2m2)p — mh\j (mh — 4m| )(m| + p2)
(A.11)
Employing arguments on continuity and limits of cos a0 for p' ^ to and p' ^ 0, these two have to be the boundaries of the p'-integral, so the final form of this collision term is:
p1 ' '
CS^ss[/s](p, T) = -^p^J pdp- |M|2[/Sq(p,T)/Sq(p',T) — /s(p,T)/s(p',T)] . (A.12)
Having seen one concrete example, we will now list the full set of collision terms needed to reproduce our results:
• 2 ^ 2-scattering processes. All scattering processes are of the form CLss [fs ](P,T )= (A.13)
g2 r p'2dp' fcos amax w I 4m2
j d(cos ah/1 — j
16 Jm2 + p2(2n)3 0 A/m2 + p'2 V-i V S(p,p''cosa, m)
,2 _ s + 2mh
|Mss^hh|2 = Mhh^ss |2 = 32A2 —-f , (A.16)
Vs — mh J
Mss^f = IM^ss|2 = 8A2m2 ( ,^ ' 2r2 , (A.17)
(s — mh) + mhr h
|mss^w+w-12 = |mw+wss|2 = -gA2mw^+^miprw ' (A.18)
i »> i2 i »> i2 8 ,2 s2 — 4mZ s + 12mZ in,
|Ms's-^zzI2 = = gA (s _„,£ + „,hrf •
x |Mss^«(p,p', cosa)|2(/!>,T) /eq(p',T) — fs(p, T) fs(p', T}) ,
where ii = hh, ii, W+WZZ and f^q(p) = exp (— jJp2 + ms/T) is the (would-be) equilibrium distribution of the scalar S, which is of Boltzmann-shape by virtue of the principle of detailed balance. Furthermore, s is the square of the centre-of-mass energy, explicitly given by:
s(p,p', cos a, ms) = 2(ms + y (ms + p2)(ms + p'2) — pp' cos a) • (A.14)
In addition, p is the momentum of the second scalar participating in the process and a is the angle between p and p', whose maximum value is given by cos amax = min { max{cos aim, —1}, 1}, defined by 4m2 = s(p,p', cos aim, m», ms). This upper integration boundary excludes all values of cos a for which the integral becomes imaginary due to the square root contained in eq. (A.13).
The spin-averaged matrix elements including a factor of 2 in all of the following because ——^ two scalars are annihilated or produced, respectively, in each process, are given by (assuming CP-invariance):
Mss^mI2 = IM^ss |2 = 32A2, (A.15)
• 1 ^ 2-decay processes. Several different decays have to be taken into account. Starting with the decay of a SM-Higgs into two singlet scalars, the corresponding collision term
14Here, 0 denotes any component of the SM-Higgs doublet
15The hh o SS result is only to leading order in A. All others are full tree-level results. These values also include appropriate factors to account for identical particles in the initial or final state.
is given by
CLss [/S](p,t) = mh—ss|2
+ p2 p'dp'
x/ (/eq(p, t ) /ev, t ) - /s (p, t ) /S (p', t »,
-'p'min Jm2S + p'2
(A.20)
/ I ^ik ^_ (^l2 _ 2^12 I / ^ik C l( ^l2 _ 2^12
with boundaries p'min = |-— I and p'max = -— , where =
^(mh — 4mS)(mS + p2), as well as the matrix element including a factor 2 because of annihilation/production of two scalars:
Mh—ss |2 = Mss—h|2 = 16A2v2. (A.21)
Furthermore, there are two collision terms related to the decay of a singlet scalar S into two sterile neutrinos N. The first is the one used in the Boltzmann equation for S,
CS^nn[/s](p,T) = — -rs^NN/s(p,T), (A.22)
VmS + p2 2
with the decay widths rs—NN = y2ms/(16n). The second version is the one used in the Boltzmann equation for N,
CN-NN[/s](p, T) = msrp^ / dp' p /s(p',T) , (A.23)
p / A/m2 + p'2
P'min,N V ^
with p'min n = |p — ip |. Note that the two collision terms Cj— NN and C|—NN appear to be somewhat different, although one may very naively expect one to be just the negative of each other. However, we should keep in mind that we are working on the level of momentum distribution functions, which implies that the collision terms look rather different depending on whether or not the desired distribution function is integrated over.
In these equations it is understood that p, p', and T are substituted in favour of £, £', and r, as specified in eq. (3.9).
A.2 Transformation of variables
As stated in the main text, we perform a transformation of variables in order to bring the Liouville operator L = Jt — Hpdp into a more convenient form. To see how this results into eq. (3.10), consider a general transformation into new variables r and £:
p} - {? = , (A-)
These new variables can be inserted into the Liouville operator L:
- dr d d TT , .. /dr d d \ . . .
L = dr dr + dtdde— Hp(r, dpddr + dpdd^J. (A.25)
To simplify L, we need to get rid of one of the two differential operators. In other words, the new variables would be most useful if we could choose them in such a way that they transform the partial differential equation into an effective ordinary differential equation. The first step is to demand that r does not depend on p, which eliminates one term:
S dr d dt dr
| — f
Next, we demand
I = Hp(r.e) |
This is a rather simple partial differential equation. Fixing the initial condition
£(p,to) = Co(p),
where £o is some arbitrary C 1-function, it has the simple solution
e(p,t) = eo
(A.26)
(A.27)
(A.28)
(A.2g)
This implies that, if we fulfill the requirements that r only depends on t and the dependence of e on p and t is given by eq. (A.29), the Liouville operator in terms of the new coordinates will have the simple form
L = (A-30)
We can make our life even easier by a smart choice for the functions r(t) and £0. Exploiting the one-to-one correspondence between temperature T and time t, one possible choice is:
and e =
1 a(t) To a(t(To))
ffs(To)\
(A.31)
for some reference mass mo and some reference temperature To, both of which we choose to equal the Higgs mass:
mo = To = mh • (A.32)
For the last equality in eq. (A.31), we have used the fact that the comoving entropy density s is constant,
s(T)a(T)3 = ^gs(T) T3 a3(T) = const., (A.33)
which allows to relate the scale factor a(T) to the effective number gs(T) of relativistic entropy degrees of freedom. Eq. (A.33) can also be used to derive the time-temperature relation
dT = WTm+1
dt V 3gs(T)
Plugging this into the Liouville operator from eq. (A.30) yields its final form,
L rH ( 3gs +0 <9r •
(A.34)
(A.35)
CTi —
This completes the proof of eq. (3.10).
B The freeze-out of the Higgs boson
We argued that we do not have to solve a system of Boltzmann equations for all species in the early Universe, but only for the scalar singlet S and the sterile neutrino N, since we rely on the SM particles sourcing the production of S being in thermal equilibrium. This assumption is for sure good for ms » mH, but we should assess its quality if we want to proceed to smaller ms. The smaller ms the later (in cosmic time) the scalar will be produced in general, and hence the less reliable the assumption of particles like the Higgs or gauge bosons being in equilibrium could be. In fact, all species relevant for sourcing S in the broken phase (W±, Z, h, t) have similar masses and are therefore expected to decouple from the thermal plasma at a similar time.
Still, if we restrict ourselves to ms > 30 GeV, it is enough to know whether the source particles are still in equilibrium at the corresponding temperature. Even if they are in O equilibrium much longer, all the particles mentioned will by then have disappeared due to their masses resulting in strong Boltzmann-suppressions. Since we have not found any ^^ detailed source to cross-check this assumption, we have assessed it ourselves with the following analysis.
Let us assume that the Higgs boson h and the top quark t decouple from the plasma |—1 before the gauge-bosons, since their mass is larger by a factor of O (1). We therefore looked at the following system of coupled Boltzmann equations:
L/h = ^h^SM'SM7[/h] + ^hh^SM'SM7[/h] + ^h^iit/^ /i] ,
L = Csm'sm [/h] + Ctt^SM'SM7 [/h ] + Ch^ti[/h, /i] , (B.1)
where SM' denotes all SM degrees of freedom except for the t and the h. The matrix elements ^^ going into all the collision terms were computed at tree-level only. The Higgs decay width was taken from ref. [112].
Solving them on the level of rate equations, we find that in this system, both the Higgs and the top closely track their equilibrium abundance all the way down to temperatures of about 1 GeV, where their abundances are exponentially suppressed already. When switching off inverse decays and taking into account only two-to-two-processes, the top would freeze out around 5 GeV while the Higgs would freeze out around 3 GeV. Note that the massive gauge bosons can also be produced via inverse decays, which is another argument to assume that they will decouple only after the Higgs and the top.
C Failure of the free-streaming horizon as measure for non-thermal Dark
Matter
Here, we would like to discuss one important point which we had stressed on several occasions throughout the manuscript. Given that we are dealing with highly non-thermal DM spectra, we cannot expect paradigms developed for thermal relics to carry over to this much more general case. In particular, a non-thermal spectrum cannot be described very well by a single number, such as a temperature, an average momentum, or a velocity. Given that, we have in fact no reason to believe that a quantity based on an average velocity, such as the freestreaming (FS) horizon Afs as defined in eq. (3.15), should give us any countable information on the DM spectrum. Yet, it is incorrectly used in many occasions in the literature.
In order to clearly illustrate how the FS horizon fails, we compare two versions of an example spaghetti plot, which are depicted in figure 13. Here, on the left panel, we
Figure 13. Comparison of half-mode analysis (left; figure identical to the right figure 7, except for the model-dependent bounds not being displayed here for simplicity) with the computation of the free-streaming horizon (right).
can see the plot for ms = 100 GeV how we obtained it in section 4.2 (this figure is basically identical to the right figure 7, apart from the missing collider-related bounds which we skipped here in order not to distract the reader). The identical patch of the parameter space is displayed on the right panel of figure 13, however, this time with an analysis based on the FS horizon, classifying the different points as "cold", "warm", or "hot" — even though, as already explained, these categories do not really suit any non-thermal spectrum. Comparing both plots, one can immedately see that the analysis based on the FS horizon16 is much more pessimistic than the result based on the half-mode analysis described in section 3.2. While there is no perfect correspondence of the colours red, purple, and blue between the two plots, it is nevertheless evident that some points which are not even constrained by current data in the left plot, seem to be excluded completely when looking at the right plot.
This is particularly true for the two points marked in the plots:
A : mN = 7.1 keV, A = 10-7 77, y = 10-650, B: mN = 20 keV, A = 10-425, y = 10-847.
Obviously, both these points are unconstrained (blue) in the half-mode analysis, while point A would be classified as "hot" (red) in the analyses based on the FS horizon, and even point B would still be labeled as "warm" (purple). This second classification is based on the results obtained for the FS horizon, which turn out to be:
A: Afs = 0.117MPc ^ "hot", B : Afs = 0.089 MPc ^ "warm"
which perfectly matches the classification from the right panel of figure 13.
16This we computed following the numerical computation of the FS horizon as reported in ref. [88], just with a small numerical error in gs which we have corrected in this version — the resulting plot would be in between the "numerical" and "analytic" versions of Afs used in [88].
k [h/Mpc]
Figure 14. Squared transfer functions for the points A and B marked in figure 13. Even though point A (point B) would be completely discarded (strongly constrained) in an analysis based on the FS horizon, both points are in fact in full agreement with data.
But how can we be sure that this classification is insufficient? The simplest way is to confront the DM distribution functions with the actual data, which we can do by explicitly displaying the corresponding squared transfer functions in comparison to the Lyman-a data, as done in figure 14. Looking at the curves for both points A & B, we can immediately see that both of them are not at all constrained by the data. Thus, both these points are, in fact, even indistinguishable from cold DM, from a structure formation point of view.
Given that this is truly obvious for the two example cases, it serves as a clear example for the FS horizon analysis leading to a conclusion that would be completely incorrect (namely discarding point A all along, while viewing point B as borderline case). Thus, as should now be obvious, the free-streaming horizon is not at all a suitable measure when applied to non-thermal DM distributions.
D Robustness of the halfmode analysis
At first sight, the definition of the halfmode in eq. (3.16) might seem just as arbitrary as the free-streaming boundary between "cold" and "warm" ("warm" and "hot") at 0.01 Mpc (0.1 Mpc). Even though the transfer function falls off steeply around kx/2, it is a priori unclear at which value of the squared transfer function the discrimination becomes indeed negligible. While a more fundamental analysis of structure formation (e.g. by rederiving Lyman-a bounds for non-thermal spectra) is beyond the scope of this paper and is projected for future work, we can nonetheless test the robustness of our analysis against changes in the threshold given in eq. (3.16). More specifically, we can make our analysis more restrictive by comparing not only wavenumbers smaller than k\/2 but wavenumbers smaller than some kx with x < 0.5, which translates to k\/2 < kx. In order to demsonstrate that the analysis is rather robust even against large changes in this threshold, we have reanalysed the case of a scalar with ms — 60 GeV with a kx — ko .05. Figure shows that even this rather drastic change in the threshold power (by one whole order of magnitude) inflicts rather mild
Figure 15. Compatability of the scalar decay model with mS = 60GeV for different threshold wave numbers. The left panel contains the analysis with kx = k1/2 (as used throughout section 4; the plot displayed here is basically identical to the right figure 4, except for the model-dependent bounds not being displayed to enable a better comparison), while the right panel displays kx = ko.05. The comparison shows that the changes are minor, but they will be finally specified in an upcoming work of us aiming to compared several advances methods to derive bounds from structure formation.
changes on the results. This would be dramatically different when changing the values put in as boundaries in a free-streaming analysis (cf. Figure 13).
Even though being approximate itself, our method of incorporating bounds from structure formation, using Lyman-a data derived for thermal spectra, clearly proves more reliable than the free-streaming approach.
References
[1] Planck collaboration, P.A.R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13 [arXiv:1502.01589] [inSPIRE].
[2] K.G. Begeman, A.H. Broeils and R.H. Sanders, Extended rotation curves of spiral galaxies: dark haloes and modified dynamics, Mon. Not. Roy. Astron. Soc. 249 (1991) 523 [inSPIRE].
[3] D. Merritt, The distribution of dark matter in the Coma cluster, Astrophys. J. 313 (1987) 121.
[4] D. Clowe et al., A direct empirical proof of the existence of dark matter, Astrophys. J. 648 (2006) L109 [astro-ph/0608407] [inSPIRE].
[5] C. Lage and G.R. Farrar, The bullet cluster is not a cosmological anomaly, JCAP 02 (2015) 038 [arXiv:1406.6703] [inSPIRE].
[6] G. Bertone, D. Hooper and J. Silk, Particle dark matter: evidence, candidates and constraints, Phys. Rept. 405 (2005) 279 [hep-ph/0404175] [inSPIRE].
[7] J.R. Primack, Dark matter and structure formation, astro-ph/9707285 [inSPIRE].
[8] V. Springel et al., Simulating the joint evolution of quasars, galaxies and their large-scale distribution, Nature 435 (2005) 629 [astro-ph/0504097] [inSPIRE].
[9] B.W. Lee and S. Weinberg, Cosmological lower bound on heavy neutrino masses, Phys. Rev. Lett. 39 (1977) 165 [inSPIRE].
J. Bernstein, L.S. Brown and G. Feinberg, The cosmological heavy neutrino problem revisited, Phys. Rev. D 32 (1985) 3261 [inSPIRE].
K. Abazajian, E.R. Switzer, S. Dodelson, K. Heitmann and S. Habib, Nonlinear cosmological matter power spectrum with massive neutrinos: the halo model, Phys. Rev. D 71 (2005) 043507 [astro-ph/0411552] [inSPIRE].
R. de Putter et al., New neutrino mass bounds from Sloan Digital Sky Survey III data release 8 photometric luminous galaxies, Astrophys. J. 761 (2012) 12 [arXiv:1201.1909] [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].
B. Moore et al., Dark matter substructure within galactic halos, Astrophys. J. 524 (1999) L19 [astro-ph/9907411] [inSPIRE].
A. Klypin, I. Karachentsev, D. Makarov and O. Nasonova, Abundance of field galaxies, Mon. Not. Roy. Astron. Soc. 454 (2015) 1798 [arXiv:1405.4523] [inSPIRE].
M. Boylan-Kolchin, J.S. Bullock and M. Kaplinghat, Too big to fail? The puzzling darkness of massive Milky Way subhaloes, Mon. Not. Roy. Astron. Soc. 415 (2011) L40
[arXiv:1103.0007] [inSPIRE]. |—
E. Papastergis, R. Giovanelli, M.P. Haynes and F. Shankar, Is there a "too big to fail" |— problem in the field?, Astron. Astrophys. 574 (2015) A113 [arXiv:1407.4665] [inSPIRE].
J. Dubinski and R.G. Carlberg, The structure of cold dark matter halos, Astrophys. J. 378 (1991) 496 [inSPIRE].
B. Moore, Evidence against dissipationless dark matter from observations of galaxy haloes, O
Nature 370 (1994) 629 [inSPIRE].
M. Schaller et al., Baryon effects on the internal structure of ACDM haloes in the EAGLE simulations, Mon. Not. Roy. Astron. Soc. 451 (2015) 1247 [arXiv:1409.8617] [inSPIRE].
T.K. Chan et al., The impact of baryonic physics on the structure of dark matter haloes: the view from the FIRE cosmological simulations, Mon. Not. Roy. Astron. Soc. 454 (2015) 2981 C 3
[arXiv:1507.02282] [inSPIRE].
M.A. Henson, D.J. Barnes, S.T. Kay, I.G. McCarthy and J. Schaye, The impact of baryons on massive galaxy clusters: halo structure and cluster mass estimates, arXiv:1607.08550 [inSPIRE].
V.K. Narayanan, D.N. Spergel, R. Dave and C.-P. Ma, Constraints on the 'mass of warm dark matter particles and the shape of the linear power spectrum from the Lya forest, Astrophys. J. 543 (2000) L103 [astro-ph/0005095] [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].
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].
M. Viel, G.D. Becker, J.S. Bolton and M.G. Haehnelt, Warm dark matter as a solution to the small scale crisis: new constraints from high redshift Lyman-a forest data, Phys. Rev. D 88 (2013) 043502 [arXiv:1306.2314] [inSPIRE].
N. Menci, F. Fiore and A. Lamastra, The evolution of active galactic nuclei in warm dark matter cosmology, Astrophys. J. 766 (2013) 110 [arXiv:1302.2000] [inSPIRE].
N. Menci, N.G. Sanchez, M. Castellano and A. Grazian, Constraining the warm dark matter particle mass through ultra-deep UV luminosity functions at z = 2, Astrophys. J. 818 (2016) 90 [arXiv:1601.01820] [inSPIRE].
[31 [32 [33 [34 [35 [36 [37 [38
[40 [41 [42
[43 [44 [45 [46 [47 [48
N. Menci, A. Grazian, M. Castellano and N.G. Sanchez, A stringent limit on the 'warm dark matter particle masses from the abundance of z = 6 galaxies in the Hubble frontier fields, Astrophys. J. 825 (2016) L1 [arXiv:1606.02530] [inSPIRE].
N. Yoshida, A. Sokasian, L. Hernquist and V. Springel, Early structure formation and reionization in a warm dark matter cosmology, Astrophys. J. 591 (2003) L1 [astro-ph/0303622] [inSPIRE].
M.R. Lovell et al., The properties of warm dark matter haloes, Mon. Not. Roy. Astron. Soc. 439 (2014) 300 [arXiv:1308.1399] [inSPIRE].
S. Bose et al., Substructure and galaxy formation in the Copernicus Complexio warm dark matter simulations, arXiv:1604.07409 [inSPIRE].
R. Adhikari et al., A white paper on keV sterile neutrino dark matter, arXiv:1602.04816 [inSPIRE].
N.A. Ky and N.T.H. Van, Scalar sextet in the 331 model with right-handed neutrinos, Phys. Rev. D 72 (2005) 115017 [hep-ph/0512096] [inSPIRE].
A.G. Dias, C.A. de S. Pires and P.S. Rodrigues da Silva, Naturally light right-handed,
neutrinos in a 3-3-1 'model, Phys. Lett. B 628 (2005) 85 [hep-ph/0508186] [inSPIRE]. |—
M. Shaposhnikov, A possible symmetry of the vMSM, Nucl. Phys. B 763 (2007) 49 |—1
[hep-ph/0605047] [inSPIRE].
D. Cogollo, H. Diniz and C.A. de S. Pires, KeV right-handed neutrinos from type-II seesaw mechanism in a 3-3-1 model, Phys. Lett. B 677 (2009) 338 [arXiv:0903.0370] [inSPIRE].
A.G. Dias, C.A. de S. Pires and P.S. Rodrigues da Silva, The left-right SU(3)l x SU(3)r x U(1)x model with light, keV and heavy neutrinos, | I
Phys. Rev. D 82 (2010) 035013 [arXiv:1003.3260] [inSPIRE].
M. Lindner, A. Merle and V. Niro, Soft Le-LM-LT flavour symmetry breaking and sterile neutrino keV dark matter, JCAP 01 (2011) 034 [Erratum ibid. 07 (2014) E01] [arXiv:1011.4950] [inSPIRE].
A. Kusenko, F. Takahashi and T.T. Yanagida, Dark matter from split seesaw, Phys. Lett. B 693 (2010) 144 [arXiv:1006.1731] [inSPIRE].
A. Merle and V. Niro, Deriving models for keV sterile neutrino dark matter with the Froggatt-Nielsen mechanism, JCAP 07 (2011) 023 [arXiv:1105.5136] [inSPIRE].
T. Araki and Y.F. Li, Qe flavor symmetry model for the extension of the minimal standard model by three right-handed sterile neutrinos, Phys. Rev. D 85 (2012) 065016 [arXiv:1112.5819] [inSPIRE].
A. Adulpravitchai and R. Takahashi, A4 flavor models in split seesaw mechanism, JHEP 09 (2011) 127 [arXiv:1107.3829] [inSPIRE].
H. Zhang, Light sterile neutrino in the minimal extended seesaw, Phys. Lett. B 714 (2012) 262 [arXiv:1110.6838] [inSPIRE].
D.J. Robinson and Y. Tsai, KeV warm dark matter and composite neutrinos, JHEP 08 (2012) 161 [arXiv:1205.0569] [inSPIRE].
N.E. Mavromatos and A. Pilaftsis, Anomalous Majorana neutrino masses from torsionful quantum gravity, Phys. Rev. D 86 (2012) 124038 [arXiv:1209.6387] [inSPIRE].
P.S.B. Dev and A. Pilaftsis, Minimal radiative neutrino mass mechanism for inverse seesaw models, Phys. Rev. D 86 (2012) 113001 [arXiv:1209.4051] [inSPIRE].
R. Takahashi, Separate seesaw and its applications to dark matter and baryogenesis, Prog. Theor. Exp. Phys. 2013 (2013) 063B04 [arXiv:1303.0108] [inSPIRE].
[49 [50 [51 [52 [53
[54 [55 [56 [57 [58 [59 [60
[61 [62 [63 [64 [65 [66
D. Borah and R. Adhikari, Common radiative origin of active and sterile neutrino masses, Phys. Lett. B 729 (2014) 143 [arXiv:1310.5419] [inSPIRE].
A. Merle, keV neutrino model building, Int. J. Mod. Phys. D 22 (2013) 1330020 [arXiv:1302.2625] [inSPIRE].
D.J. Robinson and Y. Tsai, Dynamical framework for KeV Dirac neutrino warm dark matter, Phys. Rev. D 90 (2014) 045030 [arXiv:1404.7118] [inSPIRE].
E. Bulbul et al., Detection of an unidentified emission line in the stacked X-ray spectrum of 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].
HlTOMl collaboration, F.A. Aharonian et al., Hitomi constraints on the 3.5keV line in the Ci
Perseus galaxy cluster, arXiv:1607.07420 [inSPIRE].
J.P. Conlon, F. Day, N. Jennings, S. Krippendorf and M. Rummel, Consistency of Hitomi, XMM-Newton and Chandra 3.5keV data from Perseus, arXiv:1608.01684 [inSPIRE].
P. Langacker, On the cosmological production of light sterile-neutrinos, Report UPR-0401T (1989) [inSPIRE].
U. Seljak, A. Makarov, P. McDonald and H. Trac, Can sterile neutrinos be the dark matter?, Phys. Rev. Lett. 97 (2006) 191303 [astro-ph/0602430] [inSPIRE].
K. Enqvist, K. Kainulainen and J. Maalampi, Resonant neutrino transitions and, nucleosynthesis, Phys. Lett. B 249 (1990) 531 [inSPIRE].
X.-D. Shi and G.M. Fuller, New dark matter candidate: nonthermal sterile neutrinos, Phys. Rev. Lett. 82 (1999) 2832 [astro-ph/9810076] [inSPIRE].
K. Abazajian, G.M. Fuller and M. Patel, Sterile neutrino hot, warm and cold dark matter, Phys. Rev. D 64 (2001) 023501 [astro-ph/0101524] [inSPIRE].
M. Laine and M. Shaposhnikov, Sterile neutrino dark matter as a consequence of vMSM-induced lepton asymmetry, JCAP 06 (2008) 031 [arXiv:0804.4543] [inSPIRE].
C.T. Kishimoto and G.M. Fuller, Lepton number-driven sterile neutrino production in the early universe, Phys. Rev. D 78 (2008) 023524 [arXiv:0802.3377] [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].
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].
J. Ghiglieri and M. Laine, Improved determination of sterile neutrino dark matter spectrum, JHEP 11 (2015) 171 [arXiv:1506.06752] [inSPIRE].
S. Dodelson and L.M. Widrow, Sterile-neutrinos as dark matter, Phys. Rev. Lett. 72 (1994) 17 [hep-ph/9303287] [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. Merle, A. Schneider and M. Totzauer, Dodelson-Widrow production of sterile neutrino dark matter with non-trivial initial abundance, JCAP 04 (2016) 003 [arXiv:1512.05369] ( J
[inSPIRE].
[70 [71 [72 [73
[74 [75 [76 [77 [78 [79 [80 [81 [82 [83 [84 [85 [86 [87 [88 [89
T. Venumadhav, F.-Y. Cyr-Racine, K.N. Abazajian and C.M. Hirata, Sterile neutrino dark matter: weak interactions in the strong coupling epoch, Phys. Rev. D 94 (2016) 043515 [arXiv:1507.06655] [inSPIRE].
A. Merle and A. Schneider, Production of sterile neutrino dark matter and the 3.5keV line, Phys. Lett. B 749 (2015) 283 [arXiv:1409.6311] [inSPIRE].
S. Horiuchi et al., Properties of resonantly produced sterile neutrino dark matter subhaloes, Mon. Not. Roy. Astron. Soc. 456 (2016) 4346 [arXiv:1512.04548] [inSPIRE].
A. Schneider, Astrophysical constraints on resonantly produced sterile neutrino dark matter, JCAP 04 (2016) 059 [arXiv:1601.07553] [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].
A.V. Patwardhan, G.M. Fuller, C.T. Kishimoto and A. Kusenko, Diluted equilibrium sterile neutrino dark matter, Phys. Rev. D 92 (2015) 103509 [arXiv:1507.01977] [inSPIRE]. |—
S.F. King and A. Merle, Warm dark matter from keVins, JCAP 08 (2012) 016 |—\
[arXiv:1205.0551] [inSPIRE].
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 o [arXiv:0912.0390] [inSPIRE]. |_
F. Bezrukov and D. Gorbunov, Relic gravity waves and 7keV dark matter from a GeV scale inflaton, Phys. Lett. B 736 (2014) 494 [arXiv:1403.4638] [inSPIRE].
A. Kusenko, Sterile neutrinos, dark matter and the pulsar velocities in models with a Higgs
singlet, Phys. Rev. Lett. 97 (2006) 241301 [hep-ph/0609081] [inSPIRE]. CD
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].
A. Kusenko, Sterile neutrinos: the dark side of the light fermions, Phys. Rept. 481 (2009) 1 [arXiv:0906.2968] [inSPIRE].
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].
A. Merle and M. Totzauer, keV sterile neutrino dark matter from singlet scalar decays: basic concepts and subtle features, JCAP 06 (2015) 011 [arXiv:1502.01011] [inSPIRE].
M. Klasen and C.E. Yaguna, Warm and cold fermionic dark matter via freeze-in, JCAP 11 (2013) 039 [arXiv:1309.2777] [inSPIRE].
Z. Kang, Upgrading sterile neutrino dark matter to FImP using scale invariance, Eur. Phys. J. C 75 (2015) 471 [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].
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].
P. Humbert, M. Lindner and J. Smirnov, The inverse seesaw in conformal electro-weak symmetry breaking and phenomenological consequences, JHEP 06 (2015) 035 [arXiv:1503.03066] [inSPIRE].
[92 [93 [94 [95 [96 [97 [98 [99 100 101 102
S. Yaser Ayazi, S.M. Firouzabadi and S.P. Zakeri, Freeze-in production of fermionic dark matter with pseudo-scalar and phenomenological aspects, J. Phys. G 43 (2016) 095006 [arXiv:1511.07736] [inSPIRE].
J. McDonald, Warm dark matter via ultra-violet freeze-in: reheating temperature and non-thermal distribution for fermionic Higgs portal dark matter, JCAP 08 (2016) 035 [arXiv:1512.06422] [inSPIRE].
A. Adulpravitchai and M.A. Schmidt, Sterile neutrino dark matter production in the neutrino-phillic two Higgs doublet model, JHEP 12 (2015) 023 [arXiv:1507.05694] [inSPIRE].
B. Shakya, Sterile neutrino dark matter from freeze-in,
Mod. Phys. Lett. A 31 (2016) 1630005 [arXiv:1512.02751] [inSPIRE].
L. Heurtier and D. Teresi, Dark matter and observable lepton flavour violation, |
arXiv:1607.01798 [inSPIRE].
K. Kaneta, Z. Kang and H.-S. Lee, Right-handed neutrino dark matter under the B — L gauge interaction, arXiv:1606.09317 [inSPIRE].
D. Boyanovsky, Clustering properties of a sterile neutrino dark matter candidate,
Phys. Rev. D 78 (2008) 103505 [arXiv:0807.0646] [inSPIRE]. |—1
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].
A. Biswas and A. Gupta, Freeze-in production of sterile neutrino dark matter in U(1)b_l model, JCAP 09 (2016) 044 [arXiv:1607.01469] [inSPIRE].
A. Abada, G. Arcadi and M. Lucente, Dark matter in the minimal inverse seesaw mechanism, o JCAP 10 (2014) 001 [arXiv:1406.6556] [inSPIRE]. |_
L. Lello and D. Boyanovsky, Cosmological implications of light sterile neutrinos produced after the QCD phase transition, Phys. Rev. D 91 (2015) 063502 jarXiv:1411.2690] [inSPIRE].
L. Lello and D. Boyanovsky, The case for mixed dark matter from sterile neutrinos,
JCAP 06 (2016) 011 [arXiv:1508.04077] [inSPIRE]. O
S. Nurmi, T. Tenkanen and K. Tuominen, Inflationary imprints on dark matter, (jU
JCAP 11 (2015) 001 [arXiv:1506.04048] [inSPIRE].
K. Kainulainen, S. Nurmi, T. Tenkanen, K. Tuominen and V. Vaskonen, Isocurvature constraints on portal couplings, JCAP 06 (2016) 022 [arXiv:1601.07733] [inSPIRE].
M. Drewes and J.U. Kang, Sterile neutrino dark matter production from scalar decay in a thermal bath, JHEP 05 (2016) 051 [arXiv:1510.05646] [inSPIRE].
M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, Observational constraints on decoupled hidden sectors, Phys. Rev. D 94 (2016) 063506 [arXiv:1604.02401] [inSPIRE].
N. Sekiya, N.Y. Yamasaki and K. Mitsuda, A search for a keV signature of radiatively decaying dark matter with Suzaku XIS observations of the X-ray diffuse background, Publ. Astron. Soc. Jap. 68 (2016) S31 [arXiv:1504.02826] [inSPIRE].
M. Frigerio, T. Hambye and E. Masso, Sub-GeV dark matter as pseudo-Nambu-Goldstone bosons from the seesaw scale, Phys. Rev. X 1 (2011) 021026 [arXiv:1107.4564] [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].
O. Wantz and E.P.S. Shellard, Axion cosmology revisited, Phys. Rev. D 82 (2010) 123508 [arXiv:0910.1066] [inSPIRE].
L. Shampine and M. Reichelt, The MATLAB ODE suite, SIAM J. Sa. Comp. 18 (1997) 1.
111 112
120 121 122
L. Shampine, M. Reichelt and J. Kierzenka, Solving Index-1 DAEs in MATLAB and Simulink, SIAM Rev. 41 (1999) 538.
Particle Data Group collaboration, K.A. Olive et al., Review of particle physics, Chin. Phys. C 38 (2014) 090001 [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].
J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS). Part IV: Efficient implementation of non-cold relics, JCAP 09 (2011) 032 [arXiv:1104.2935] [inSPIRE].
G. Mangano et al., Relic neutrino decoupling including flavor oscillations, Nucl. Phys. B 729 (2005) 221 [hep-ph/0506164] [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].
F. Bezrukov and D. Gorbunov, Applicability of approximations used in calculations of the spectrum of dark matter particles produced in particle decays, Phys. Rev. D 93 (2016) 063502 [arXiv:1412.1341] [inSPIRE]. CD
K. Markovic and M. Viel, Lyman-a forest and cosmic 'weak lensing in a warm dark matter | I
universe, Publ. Astron. Soc. Austral. 31 (2014) e006 [arXiv:1311.5223] [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].
L.A. Popa and D. Tonoiu, Subdominant dark matter sterile neutrino resonant production in the light of Planck, JCAP 09 (2015) 066 [arXiv:1501.06355] [inSPIRE]. (jJ
J. Baur, N. Palanque-Delabrouille, C. Yeche, C. Magneville and M. Viel, Lyman-a forests cool QQ warm dark matter, in SDSS-IV Collaboration Meeting, Madrid Spain, 20-23 Jul 2015 [JCAP 08 (2016) 012] [arXiv:1512.01981] [inSPIRE].
B. Bozek et al., Resonant sterile neutrino dark matter in the local and high-z universe, Mon. Not. Roy. Astron. Soc. 459 (2016) 1489 [arXiv:1512.04544] [inSPIRE].
N. Menci, A. Schneider and M. Viel, private communication.
A.V. Maccio, O. Ruchayskiy, A. Boyarsky and J.C. Muñoz-Cuartas, The inner structure of haloes in cold+warm dark matter models, Mon. Not. Roy. Astron. Soc. 428 (2013) 882 [arXiv:1202.2858] [inSPIRE].
D. Anderhalden, A. Schneider, A.V. Maccioo, J. Diemand and G. Bertone, Hints on the nature of dark matter from the properties of Milky Way satellites, JCAP 03 (2013) 014 [arXiv:1212.2967] [inSPIRE].
A. Schneider, Structure formation with suppressed small-scale perturbations, Mon. Not. Roy. Astron. Soc. 451 (2015) 3117 [arXiv:1412.2133] [inSPIRE].
C.-P. Ma and E. Bertschinger, Cosmological perturbation theory in the synchronous and conformal Newtonian gauges, Astrophys. J. 455 (1995) 7 [astro-ph/9506072] [inSPIRE].
W.H. Press and P. Schechter, Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation, Astrophys. J. 187 (1974) 425 [inSPIRE].
[131] R.K. Sheth and G. Tormen, Large-scale bias and the peak background split, Mon. Not. Roy. Astron. Soc. 308 (1999) 119 [astro-ph/9901122]
[132] A. Schneider, R.E. Smith and D. Reed, Halo mass function and the free streaming scale, Mon. Not. Roy. Astron. Soc. 433 (2013) 1573 [arXiv:1303.0839] [inSPIRE].
[133] J. Bernstein, Kinetic theory in the expanding universe, Cambridge University Press, Cambridge U.K. (1988) [inSPIRE].
[134] E.W. Kolb and M.S. Turner, The early universe, Westview Press (1994).
