"Non-cold" dark matter at small scales: a general approach

TL;DR

This paper introduces a flexible analytical formula for the power spectrum of non-thermal dark matter models, enabling better modeling of small-scale structure formation and providing preliminary astrophysical constraints.

Contribution

It presents a new fitting formula that captures a wide range of non-thermal dark matter models, improving upon the traditional warm dark matter approach.

Findings

The formula accurately reproduces clustering signals of various non-thermal DM models.

N-body simulations using the formula reveal detailed nonlinear structure statistics.

Preliminary constraints from satellite counts and Lyman-alpha forest data are obtained.

Abstract

Structure formation at small cosmological scales provides an important frontier for dark matter (DM) research. Scenarios with small DM particle masses, large momenta or hidden interactions tend to suppress the gravitational clustering at small scales. The details of this suppression depend on the DM particle nature, allowing for a direct link between DM models and astrophysical observations. However, most of the astrophysical constraints obtained so far refer to a very specific shape of the power suppression, corresponding to thermal warm dark matter (WDM), i.e., candidates with a Fermi-Dirac or Bose-Einstein momentum distribution. In this work we introduce a new analytical fitting formula for the power spectrum, which is simple yet flexible enough to reproduce the clustering signal of large classes of non-thermal DM models, which are not at all adequately described by the oversimplified notion of WDM. We show that the formula is able to fully cover the parameter space of sterile neutrinos (whether resonantly produced or from particle decay), mixed cold and warm models, fuzzy dark matter, as well as other models suggested by effective theory of structure formation (ETHOS). Based on this fitting formula, we perform a large suite of N-body simulations and we extract important nonlinear statistics, such as the matter power spectrum and the halo mass function. Finally, we present first preliminary astrophysical constraints, based on linear theory, from both the number of Milky Way satellites and the Lyman-$\alpha$ forest. This paper is a first step towards a general and comprehensive modeling of small-scale departures from the standard cold DM model.

Figures

a of Refs. [[32](#bib.bib32), [33](#bib.bib33)], and we are also grateful for Nicola Menci for useful discussions on the small-scale part of the transfer functions. 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). RM and MV are supported by the INFN PD51 grant INDA

Tables

Table 1. Table 1 : Each { α , β , γ } 𝛼 𝛽 𝛾 \{\alpha,\beta,\gamma\} -combination corresponds to a different nCDM model. The first 27 models are associated to the combinations of { α , β , γ } 𝛼 𝛽 𝛾 \{\alpha,\beta,\gamma\} that constitute the 3 × 3 × 3 3 3 3 3\times 3\times 3 non-regular grid in the parameter space, while the last 28 models refer to the additional triplets that we have considered. In the last column we list the corresponding values of k 1 / 2 subscript 𝑘 1 2 k_{1/2} .

	$\beta$	$\gamma$	$\alpha$ (Mpc/$h$)	$k_{1/2}$ ($h$/Mpc)
nCDM1	1.5	$-10.0$	0.0075	14.323
nCDM2	1.5	$-10.0$	0.0048	22.463
nCDM3	1.5	$-10.0$	0.0035	30.914
nCDM4	1.5	$-5.0$	0.012	14.323
nCDM5	1.5	$-5.0$	0.0077	22.463
nCDM6	1.5	$-5.0$	0.0056	30.914
nCDM7	1.5	$-1.0$	0.039	14.323
nCDM8	1.5	$-1.0$	0.025	22.463
nCDM9	1.5	$-1.0$	0.018	30.914
nCDM10	2.0	$-10.0$	0.013	14.323
nCDM11	2.0	$-10.0$	0.0084	22.463
nCDM12	2.0	$-10.0$	0.0061	30.914
nCDM13	2.0	$-5.0$	0.019	14.323
nCDM14	2.0	$-5.0$	0.012	22.463
nCDM16	2.0	$-1.0$	0.045	14.323
nCDM15	2.0	$-5.0$	0.0087	30.914
nCDM17	2.0	$-1.0$	0.029	22.463
nCDM18	2.0	$-1.0$	0.021	30.914
nCDM19	2.5	$-10.0$	0.018	14.323
nCDM20	2.5	$-10.0$	0.012	22.463
nCDM21	2.5	$-10.0$	0.0085	30.914
nCDM22	2.5	$-5.0$	0.024	14.323
nCDM23	2.5	$-5.0$	0.016	22.463
nCDM24	2.5	$-5.0$	0.011	30.914
nCDM25	2.5	$-1.0$	0.049	14.323
nCDM26	2.5	$-1.0$	0.031	22.463
nCDM27	2.5	$-1.0$	0.023	30.914

Table 2. Table 2 : Number of subhalos (with mass M sub ≥ 10 8 ​ M ⊙ / h subscript 𝑀 sub superscript 10 8 subscript 𝑀 direct-product ℎ M_{\rm sub}\geq 10^{8}~{}M_{\odot}/h ) within the virial radius of a halo with mass M halo = 1.7 ⋅ 10 12 ​ M ⊙ / h subscript 𝑀 halo ⋅ 1.7 superscript 10 12 subscript 𝑀 direct-product ℎ M_{\rm halo}=1.7\cdot 10^{12}~{}M_{\odot}/h . Each of the 55 models corresponds to a different { α , β , γ } 𝛼 𝛽 𝛾 \{\alpha,\beta,\gamma\} -combination, according to Tab.  1 . Models that predict a number of subhalos consistent with observations (i.e., predicting at least as many subhalos as the number of observed MW satellites, N sat = 63 subscript 𝑁 sat 63 N_{\rm sat}=63 ) are written in bold-face, while those not surviving the reality check are not bold-faced.

	$N_{\mathrm{sub}}$
nCDM1	39
nCDM2	78
nCDM3	105
nCDM4	25
nCDM5	68
nCDM6	104
nCDM7	18
nCDM8	59
nCDM9	103
nCDM10	41
nCDM11	78

Table 3. Table 3 : Here we list the 55 models that we have studied, each of them with its corresponding δ ​ A 𝛿 𝐴 \delta A , namely the estimator of the small-scale power suppression associated to it. A model is excluded (at 95% C.L.) if δ ​ A > δ ​ A REF , 1 𝛿 𝐴 𝛿 subscript 𝐴 REF 1 \delta A>\delta A_{\rm REF,1} , i.e., if it shows a power suppression ≳ 38 % greater-than-or-equivalent-to absent percent 38 \gtrsim 38\% with respect to the Λ Λ \Lambda CDM power spectrum. Each of the models corresponds to a different { α , β , γ } 𝛼 𝛽 𝛾 \{\alpha,\beta,\gamma\} -combination, according to Tab.  1 . Accepted models are highlighted in bold-face.

	$\delta A$
nCDM1	0.61
nCDM2	0.45
nCDM3	0.34
nCDM4	0.63
nCDM5	0.45
nCDM6	0.32
nCDM7	0.64
nCDM8	0.45
nCDM9	0.31
nCDM10	0.61
nCDM11	0.45

Table 4. Table 4 : Here we list 20 { α , β , γ } 𝛼 𝛽 𝛾 \{\alpha,\beta,\gamma\} -triplets, with the corresponding values of k 1 / 2 subscript 𝑘 1 2 k_{1/2} , which represent the real model examples presented in Sec.  3 , split into five groups. For each case, we have confronted the fit (the real model) with both halo counting, N sub fit subscript superscript 𝑁 fit sub N^{\rm fit}_{\rm sub} ( N sub true subscript superscript 𝑁 true sub N^{\rm true}_{\rm sub} ), and the Lyman- α 𝛼 \alpha forest, δ ​ A fit 𝛿 subscript 𝐴 fit \delta A_{\rm fit} ( δ ​ A true 𝛿 subscript 𝐴 true \delta A_{\rm true} ), where for each case we also indicate the percentage (in purple for better visibility) by which the value predicted from the fit differs from the “true” value predicted by the model point. Bold-faced numbers indicate that both the restrictive and conservative bounds are met, while Italic numbers indicate that only the conservative bound was met. For each case we have indicated whether the conclusion drawn from the fit – i.e., whether or not a certain choice of parameters is allowed by the data – does ( ✓ ✓ {\large\color[rgb]{0.1,0.8,0.1}\checkmark} ) or does not ( × ) agree with the one drawn from the data for the real model. Looking at the observables, it is visible that the values predicted by the fits and the real models differ by a few per cent at most and, indeed, the very few red crosses in the table indicate that our fitting function is remarkably powerful in reproducing the correct conclusion. Thus, in most cases, it is sufficient to fit a model with our general formula and check whether the fitted point is in the allowed region.

	$\alpha$	$\beta$	$\gamma$	$k_{1/2}\left[h/\mathrm{Mpc}\right]$	$N_{\mathrm{sub}}^{\mathrm{fit}}$ ($N_{\mathrm{sub}}^{\mathrm{true}}$) [%]	Agree?	$\delta A_{\mathrm{fit}}$ ($\delta A_{\mathrm{true}}$) [%]	Agree?
	$0.025$	$2.3$	$-2.6$	17.276	38 (39) [$-2.6\%$]	$✓$	0.555 (0.571) [$-2.8\%$]	$✓$
RP	$0.071$	$2.3$	$-1.0$	9.828	15 (14) [$+7.1\%$]	$✓$	0.743 (0.754) [$-1.5\%$]	$✓$
neutrinos	$0.038$	$2.3$	$-4.4$	8.604	5 (5) [$\pm 0.0\%$]	$✓$	0.799 (0.810) [$-1.4\%$]	$✓$
	$0.035$	$2.1$	$-1.5$	15.073	35 (37) [$-5.4\%$]	$✓$	0.599 (0.613) [$-2.3\%$]	$✓$
Neutrinos	$0.016$	$2.6$	$-8.1$	19.012	38 (42) [$-9.5\%$]	$✓$	0.521 (0.535) [$-2.6\%$]	$✓$
from	$0.011$	$2.7$	$-8.5$	28.647	91 (97) [$-6.2\%$]	$✓$	0.339 (0.360) [$-5.8\%$]	$✓$
particle	$0.019$	$2.5$	$-6.9$	16.478	27 (28) [$-3.6\%$]	$✓$	0.582 (0.576) [$+1.0\%$]	$✓$
decay	$0.011$	$2.7$	$-9.8$	26.31	79 (87) [$-9.2\%$]	$✓$	0.375 (0.390) [$-3.8\%$]	×
	$0.16$	$3.2$	$-0.4$	6.743	9 (9) [$\pm 0.0\%$]	$✓$	0.823 (0.834) [$-1.3\%$]	$✓$
Mixed	$0.20$	$3.7$	$-0.18$	7.931	28 (27) [$+3.7\%$]	$✓$	0.738 (0.752) [$-1.9\%$]	$✓$
models	$0.21$	$3.7$	$-0.1$	11.36	60 (62) [$-3.2\%$]	$✓$	0.596 (0.610) [$-2.3\%$]	$✓$
	$0.21$	$3.4$	$-0.053$	33.251	110 (114) [$-3.5\%$]	$✓$	0.365 (0.377) [$-3.2\%$]	$✓$
	$0.054$	$5.4$	$-2.3$	13.116	8 (9) [$-11.1\%$]	$✓$	0.691 (0.708) [$-2.4\%$]	$✓$
Fuzzy	$0.040$	$5.4$	$-2.1$	18.106	21 (23) [$-8.7\%$]	$✓$	0.543 (0.565) [$-3.9\%$]	$✓$
DM	$0.030$	$5.5$	$-1.9$	25.016	56 (60) [$-6.7\%$]	$✓$	0.376 (0.399) [$-5.8\%$]	×
	$0.022$	$5.6$	$-1.7$	34.590	121 (126) [$-4.0\%$]	$✓$	0.228 (0.250) [$-8.8\%$]	$✓$
	$0.0072$	$1.1$	$-9.9$	7.274	18 (19) [$-5.3\%$]	$✓$	0.780 (0.788) [$-1.0\%$]	$✓$
ETHOS	$0.013$	$2.1$	$-9.3$	16.880	36 (39) [$-7.7\%$]	$✓$	0.568 (0.581) [$-2.2\%$]	$✓$
models	$0.014$	$2.9$	$-10.0$	21.584	50 (53) [$-5.7\%$]	$✓$	0.463 (0.477) [$-2.9\%$]	$✓$
	$0.016$	$3.4$	$-9.3$	23.045	53 (56) [$-5.4\%$]	$✓$	0.430 (0.439) [$-2.1\%$]	$✓$

Equations

$$T^2(k)\equiv \frac{P_{\mathrm{nCDM}}(k)}{P_{\mathrm{CDM}}(k)},$$

$$T(k)=[1+(\alpha k{)}^{2\nu }{]}^{-5/\nu },$$

$$\begin{array}{rl}\alpha & =0.24{\left(\frac{m_x/T_x}{1\mathrm{keV}/{\mathrm{T}}_{\nu }}\right)}^{-0.83}{\left(\frac{{\omega }_x}{0.25\cdot (0.7{)}^2}\right)}^{-0.16}\mathrm{Mpc}\\ & =0.049{\left(\frac{m_x}{1\mathrm{keV}}\right)}^{-1.11}{\left(\frac{{\Omega }_x}{0.25}\right)}^{0.11}{\left(\frac{h}{0.7}\right)}^{1.22}{h}^{-1}\mathrm{Mpc},\end{array}$$

$$T(k)=[1+(\alpha k{)}^{\beta }{]}^{\gamma },$$

$$k_{1/2}=\frac{1}{\alpha }{\left[{\left(\frac{1}{\sqrt{2}}\right)}^{1/\gamma }-1\right]}^{1/\beta }.$$

$$\begin{array}{r}{k}_{1/2}^{\prime }=14.323h/\mathrm{Mpc}(\text{if thermal}:⟷m_x^{\prime }=2\mathrm{keV}),\\ k_{1/2}^{\prime \prime }=22.463h/\mathrm{Mpc}(\text{if thermal}:⟷m_x^{\prime \prime }=3\mathrm{keV}),\\ k_{1/2}^{\prime \prime \prime }=30.914h/\mathrm{Mpc}(\text{if thermal}:⟷m_x^{\prime \prime \prime }=4\mathrm{keV}).\end{array}$$

$$T(k)=[1+(\alpha k{)}^{\beta }{]}^{\gamma }\to T(k)=[1+(\alpha k{)}^{\beta }{]}^{\gamma }+\delta .$$

$$T(k)=[1+(\alpha k{)}^{\beta }{]}^{\gamma }\cdot {\cos }^2(\delta k)$$

$$\frac{\mathrm{d}N}{\mathrm{d}{M}_{\mathrm{sub}}}=\frac{1}{44.5}\frac{1}{6{\pi }^2}\frac{M_{\mathrm{halo}}}{M_{\mathrm{sub}}^2}\frac{P(1/R_{\mathrm{sub}})}{R_{\mathrm{sub}}^3\sqrt{2\pi (S_{\mathrm{sub}}-S_{\mathrm{halo}})}},$$

$$S_i=\frac{1}{2{\pi }^2}\underset{0}{\overset{1/R_i}{\int }}\mathrm{d}k{k}^{2}P(k);M_i=\frac{4\pi }{3}{\Omega }_m{\rho }_c(c{R}_i{)}^3;c=2.5;$$

$$\begin{array}{cccccccccc}\alpha \le 0.067\mathrm{Mpc}/h(95\%C.L.)& & & & & & & & & \text{requiring}{N}_{\mathrm{sub}}\ge 63,\text{hfill}\text{hfill}\text{hfill}\text{hfill}\text{hfill}\\ \alpha \le 0.061\mathrm{Mpc}/h(95\%C.L.)& & & & & & & & & \text{requiring}{N}_{\mathrm{sub}}\ge 57,\end{array}$$

$$r(k)=\frac{P_{1\mathrm{D}}(k)}{P_{1\mathrm{D}}^{\Lambda \mathrm{CDM}}(k)},$$

$$P_{1\mathrm{D}}(k)=\frac{1}{2\pi }\underset{k}{\overset{\infty }{\int }}\mathrm{d}{k}^{\prime }{k}^{\prime }P(k^{\prime }),$$

$$\delta A\equiv \frac{A_{\mathrm{\Lambda CDM}}-A}{A_{\mathrm{\Lambda CDM}}},$$

$$A=\underset{k_{\min }}{\overset{k_{\max }}{\int }}\mathrm{d}kr(k),$$

$$\begin{array}{cccccccccc}\alpha \le 0.058\mathrm{Mpc}/h(95\%C.L.),& & & & & & & & & \text{conservative analysis,}\text{hfill}\text{hfill}\text{hfill}\text{hfill}\text{hfill}\\ \alpha \le 0.044\mathrm{Mpc}/h(95\%C.L.),& & & & & & & & & \text{non-conservative analysis.}\end{array}$$

$${T\left(k\right)|}_{k\alpha \gg 1}\simeq {\alpha }^{\beta \gamma }{k}^{\beta \gamma }.$$

$${\alpha }^{\beta \gamma }{k}^{\beta \gamma }\to {\tilde{\alpha }}^{\beta \gamma }{k}^{\beta \gamma }={\alpha }^{\beta \gamma }{x}^{\beta \gamma }{k}^{\beta \gamma }={\alpha }^{\beta \gamma }{\left(k^{\ln x/\ln k}\right)}^{\beta \gamma }{k}^{\beta \gamma }={\alpha }^{\beta \gamma }{k}^{\beta \gamma (1+\ln x/\ln k)},$$

$$\gamma \to \tilde{\gamma }(k)=\gamma (1+\ln x/\ln k),$$

$$\tilde{\gamma }(k)\simeq const.$$

$$\frac{\mathrm{d}n}{\mathrm{d}M}=\frac{1}{12{\pi }^2}\frac{\overline{\rho }}{M^2}\nu f(\nu )\frac{P(1/R)}{{\delta }_c^2{R}^3},$$

Full text

1 Introduction

Cosmic microwave background and large-scale structure data have provided the consistent picture that the present universe is mainly composed by a cosmological constant, denoted by $\Lambda$, and by cold dark matter (CDM). Whereas this standard cosmological model, or $\Lambda$CDM model, provides very convincing predictions for the large-scale structures, it still exhibits some limits at very small sub-galactic scales (for a review see [1]), where tensions exist between predictions and observations. In fact, assuming the $\Lambda$CDM model, cosmological $N$-body simulations predict too many dwarf galaxies (missing satellite problem [2, 3]) and too much dark matter (DM) in the innermost regions of galaxies (cusp-core problem [4]) with respect to observations. The situation both under the theoretical/numerical and observational point of view is not entirely clear, and there are continuous efforts to exploit astrophysical consequences of DM both in (spiral) galaxies [5, 6], dwarf galaxies of our Milky Way (MW), or, even more recently, tidal streams [7]. Moreover, the dynamical properties of massive MW satellites are not reproduced in simulations (too-big-to-fail problem [8, 9]). These small-scale problems could also be alleviated or solved by invoking baryon physics, still not perfectly understood and/or fully implemented in cosmological simulations. For instance, it is known that photo-evaporation from UV sources during the reionisation period pushes gas out from small halos, preventing star formation and reducing the number of observed substructures [10]. Moreover, supernova feedback may be able to make the inner parts of halo density profiles significantly shallower [11]. However, these baryonic feedback effects are currently implemented as semi-analytical models into hydrodynamical simulations, making it difficult to obtain results with full predictive power.

Another possible solution to the small-scale crisis of the $\Lambda$CDM model is to modify the nature of DM, by going beyond the standard CDM paradigm. In fact, despite huge efforts both in particle physics and cosmology, the nature and composition of the “dark sector” are still unknown. From an astrophysical perspective, DM candidates can be classified according to either a velocity dispersion or a pressure term which counteracts the effects of gravity at small scales. Below this characteristic scale of a given DM candidate, density fluctuations are erased and gravitational clustering is thus suppressed. The velocity dispersion of CDM particles is by definition so small to be totally negligible for cosmological structure formation.

Many “non-cold” DM (nCDM) candidates, well motivated by particle physics theories, such as sterile neutrinos [12, 13, 14, 15, 16, 17] or axion-like particles [18, 19, 20, 21, 22], have been proposed in order to give a better description of the structure formation and distribution at small scales, with respect to the $\Lambda$CDM model. In addition, there are other non-standard hypotheses potentially able to induce a small-scale suppression in the matter power spectrum: a mixed (cold and warm) DM fluid [23, 24, 25, 26], DM particles coupled to dark energy [27, 28, 29] or to some relativistic fluid [30, 31], or Self-Interacting Dark Matter (SIDM) [32, 33]. Different scenarios lead to different shapes in the suppression of the power spectrum. However, most of the constraints from structure formation data which have been published so far, refer to a very specific shape of the small-scale power suppression, corresponding to the case of thermal warm dark matter (WDM), i.e., candidates with a Fermi-Dirac momentum distribution [34, 35, 36, 37]. Nonetheless, most of viable nCDM candidates do not have thermally distributed momenta (i.e., they feature truly non-thermal distribution functions), which may lead to non-trivial suppressions in their power spectra. In particular, given that non-thermal distributions may feature more than one characteristic momentum scale, these settings may in fact provide a whole new approach to resolving the small scale issues of CDM.

The suppression of gravitational clustering in nCDM models is usually parametrised through the so-called transfer function $T(k)$, i.e., the square root of the ratio of the matter power spectrum predicted by the given model with respect to that in the presence of CDM only, for fixed cosmological parameters. In this work we introduce a new, more general analytical fitting formula for the transfer function. We show that it accurately describes the behaviour of the most popular DM scenarios, being able to reproduce a large variety of shapes with only three free parameters. Finally, we present the first constraints from structure formation on its free parameters.

This paper is organised as follows: in Sec. 2 we introduce and motivate our general approach; in Sec. 3 we show that the new formula is an useful tool for fitting the behaviour of the most viable DM models provided by particle physics; in Sec. 4 we discuss the results of the cosmological simulations that we have performed in order to investigate different parametrisations of the new transfer function; in Sec. 5 we present the first preliminary astrophysical constraints on its free parameters; finally, in Sec. 6, we apply these constraints to the theoretical particle physics models previously discussed, in order to compare the predictions in terms of structure formation by our fits with the actual transfer functions and check the accuracy of our method.

2 Method and parametrisation

The suppression of gravitational clustering in nCDM scenarios can be described through the transfer function $T(k)$, given by:

[TABLE]

where $P_{\rm CDM}$ and $P_{\rm nCDM}$ are the power spectra of the CDM and the nCDM models, respectively.

For the special case of thermal WDM, the transfer function can be well approximated by the fitting function [38]:

[TABLE]

where $\alpha$ is the only free parameter and $\nu=1.12$. Therefore, constraints on the mass of the WDM candidate translate into bounds on $\alpha$, by the following formula [24]:

[TABLE]

where the subscripts $x$ and $\nu$ refer to WDM and active-neutrino properties, respectively, and the second equation holds only in the case of thermal relics.

Let us now generalize Eq. (2.2) and write down the following fitting formula111Note that equivalent fitting functions have already been used in, e.g., Refs. [39, 40, 41]. However, they have only been applied to special cases, and its general applicability had not been recognised to our knowledge.

[TABLE]

which is a function of three free parameters: $\alpha$, $\beta$, and $\gamma$. In this paper we show that the simple function given by Eq. (2.4) is generic enough to describe the majority of nCDM models from the literature.

As a next step, let us define the characteristic half-mode scale $k_{1/2}$, obtained by setting $T^{2}\equiv 1/2$. Using Eq. (2.4), we therefore have:

[TABLE]

Whereas through Eqs. (2.2) and (2.3) we had a one-to-one correspondence between $m_{x}$ and $\alpha$, constraints on the DM mass are now, by Eqs. (2.4) and (2.5), mapped to 3-dimensional surfaces in the $\{\alpha,\beta,\gamma\}$ space. In other words, given a value of $k_{1/2}$, which corresponds to a certain (thermal) WDM mass, we can easily compute the corresponding surface in the 3-dimensional parameter space from Eq. (2.5) – but this information alone is not yet sufficient to decide about the validity of the point under consideration. In Fig. 1 we plot the three surfaces associated to the $k_{1/2}$ values listed below:

[TABLE]

We also build up a $3\times 3\times 3$ non-regular grid that brackets the volume of the parameter space between the blue upper surface (corresponding to $k^{\prime}_{1/2}$) and the red lower surface (corresponding to $k^{\prime\prime\prime}_{1/2}$) in Fig. 1. This is done by taking all the possible combinations of the two triplets $\beta=\{1.5,2,2.5\}$ and $\gamma=\{-1,-5,-10\}$ and computing the corresponding values of $\alpha$, by plugging the values of $k_{1/2}$ listed in Eq. (2.6) into Eq. (2.5).

Thus, we obtain a table with 27 combinations of $\{\alpha,\beta,\gamma\}$ that sample the volume of the parameter space associated to thermal WDM masses between 2 and 4 keV, each of them corresponding to a different nCDM model. The models are listed in Tab. 1, with labels from 1 to 27, and are represented in Fig. 1 by the black points. In order to investigate even shallower and steeper transfer functions, we have furthermore considered 28 additional points in the parameter space, corresponding to the $\{\alpha,\beta,\gamma\}$-combinations listed in Tab. 1, with labels from 28 to 55, represented by the red points in Fig. 1. Note that, although the points marked in Fig. 1 may appear to be somewhat sparcely distributed at first sight, they in fact cover a large fraction of the relevant parameter space. The reason for this lies in a quasi-degeneracy between the two parameters $\alpha$ and $\gamma$, which we discuss in detail in Appendix A.

For each of the models listed in Tab. 1, we have computed the corresponding transfer function by using Eq. (2.4). We plot them in the left panel of Fig. 2, where the green and blue dashed lines represent the “old” transfer functions, i.e. computed through Eq. (2.2), for $m^{\prime}_{x}=2~{}\rm{keV}$ and $m^{\prime\prime\prime}_{x}=4~{}\rm{keV}$, respectively. Fig. 2 gives an overview over the large variety of transfer functions investigated in this paper.

Let us now qualitatively describe the role of the different parameters in the generalised fit for the transfer function, Eq. (2.4). The value of $\alpha$ gives the general scale of suppression, i.e., it is the most important parameter for setting the position of $k_{1/2}$. The parameters $\beta$ and $\gamma$ are responsible for the slope of the transfer function before and after the half-mode scale $k_{1/2}$, respectively. The parameter $\beta$ has to be positive in order to have meaningful transfer functions, since $\beta<0$ gives a $T(k)$ that differs from 1 at large scales. The larger is $\beta$, the flatter is the transfer function before $k_{1/2}$. Analogously, the larger the absolute value of $\gamma$, the sharper is the cut-off.

3 Connection to particle physics models

The purpose of this section is to see to which extent the suggested 3-parameter fitting formula is able to match the transfer functions from different nCDM models. We will focus on sterile neutrinos from resonant production, sterile neutrinos from particle decay production, mixed (cold plus warm) DM models, ultra-light axions, and another class of models suggested by the effective theory of structure formation (ETHOS).

3.1 Sterile neutrinos by resonant production

Given that keV sterile neutrinos generically mix with the active-neutrino sector, it is a natural idea to use this mixing to produce sterile neutrino DM in the early universe. While it is nowadays known that the production by non-resonant transitions (“Dodelson-Widrow mechanism”) [42, 43, 44]222Note that, contrary to previous statements in the literature [43, 45], non-resonantly produced sterile neutrinos also feature a non-thermal distribution [44], rather than a suppressed thermal spectrum. However, while this slightly changes the published numerical values of the bounds on this setting, the basic conclusion of non-resonant production being excluded remains valid (in fact, it is even made stronger) [44]. is incompatible with structure formation [34, 12], a suitable lepton number asymmetry in the early universe (whose origin is not necessarily clear, though) can resonantly enhance the active-sterile transitions (“Shi-Fuller mechanism”) and yield spectra that are more likely to be in agreement with data [46, 47, 48, 49, 50, 51] (note, however, that also this mechanism is restricted to a small successful region in the parameter space [52, 12]).

As can be seen from Refs. [48, 50, 51], the distribution functions resulting from resonant production can be highly non-thermal: typically, they feature one or more narrow peaks on top of a continuous spectrum, see Fig. 3. This figure shows different momentum distribution functions (left) and the corresponding transfer functions (right), for a few example values of the sterile neutrino mass $m_{N}$, the active-sterile mixing angle $\sin^{2}(2\theta)$, and the lepton asymmetry $\mathcal{L}_{\mu}$.333Due to the current technical limitations of the software developed in conjunction with Ref. [50], the lepton asymmetry can only be placed in the muon sector, if the package sterile-dm is used. However, as shown in Ref. [51], the results would not be altered dramatically if the lepton asymmetry was present in another sector. Note that the green curve actually features $\mathcal{L}_{\mu}\equiv 0$, i.e., a case of non-resonant production. Compared to the red curve, one can see that in this case a larger mixing angle is required to meet the correct abundance, and also the spectrum is different from the resonant cases. The plots in the right panel illustrate the corresponding transfer functions (solid lines), along with the fits obtained from Eq. (2.4) using a least squares approach. As can be seen already by eye (and is confirmed by a goodness-of-fit test), our general transfer function, Eq. (2.4), provides excellent fits to these cases, with parameter values well within the parameter space probed in this work, cf. Fig. 1. This remains true also for distributions with more than one scale, which can e.g. be seen for the orange curve in the plots. Indeed, the trick of fitting the transfer function gets us rid of many redundancies.

Given that the distributions we tested (more than shown here) are very representative for resonant production of sterile neutrino DM, we conclude that our fit function describes this class of models very well.

3.2 Sterile neutrinos from particle decays

Another potential production mechanism for sterile neutrino DM relies on the decay of a hypothetical parent particle in the early universe, whose properties (in terms of the momentum distribution) translate into those of the resulting keV sterile neutrino. A very simple case discussed frequently in the literature is that of a singlet scalar particle which may via its interactions with the Standard Model Higgs boson either thermalise (and thus be equilibrated and ultimately freeze-out) or not (and thus freeze-in), see Refs. [53, 54, 55, 56, 57, 13] for very detailed treatments. Other possibilities for parent particles, like pions [58] or electrically charged scalars [59], do not exhibit any behaviour that would be qualitatively different.

For decay production, the resulting distribution functions are highly non-thermal, with spectra not only having a shape very different from a thermal one but also featuring, in general, two distinct characteristic scales [56, 13] (or even three if a subdominant subsequent Dodelson-Widrow modification is taken into account [44]). Four example distributions are depicted in the left panel of Fig. 4, for different values of the sterile neutrino and decaying scalar masses, along with the two parameters $\lambda$ (Higgs portal) and $y$ (Yukawa coupling), which shape the distributions – see Refs. [56, 13] for details. The plots in the right panel illustrate the corresponding transfer functions (solid lines), along with the fits obtained from Eq. (2.4) using a least squares approach. As for the case of resonantly produced sterile neutrinos, the general parametrisation of Eq. (2.4), provides an excellent fit with parameter values well within the range studied in this paper (cf. Fig. 1). This remains true also for distributions with more than one scale, which can be seen for the green curve in the plots. Indeed, the trick of fitting the transfer function gets us rid of many redundancies.

Given that the distributions we tested (more than shown here) are very representative for decay production of DM, independently of the details of the particle physics setting under consideration, we conclude that our fit function describes this class of models very well.

3.3 Mixed models

In principle, the DM sector may consist of a complicated mixture of different DM particles [60]. Here we study a toy model that assumes the presence of both a cold and a warm (thermal) component. This simple model, dubbed mixed DM, leads to a large variety of shapes in the transfer function, therefore providing an ideal test for the parametrisation used in this paper (i.e. the fitting formula of Eq. (2.4)).

Mixed DM is characterised by two parameters: the mass of the WDM component and the fraction $f$ of the warm to the total DM abundance, i.e., $f=\Omega_{\rm WDM}/\Omega_{\rm total}$, where $\Omega_{\rm total}=\Omega_{\rm WDM}+\Omega_{\rm CDM}$ denotes the total DM abundance in the universe. The fraction $f$ parametrises the example settings illustrated in Fig. 5, where we depict both the power spectra (left) and the transfer functions (right). For mixed DM, it had been pointed out in Ref. [60] that a non-zero plateau can be present in the transfer function for large $k$, corresponding to the remaining CDM component once the reduction of small scales by the warm component has died off. However, although our fitting function, Eq. (2.4), does not formally feature a plateau, it provides a very good fit to the majority of mixed DM cases, not only by eye but also when performing a goodness-of-fit test.

3.4 Fuzzy dark matter

A conceptually different class of DM candidates that also affects the small scales of structure formation is the so-called Fuzzy DM. [18, 19, 20]. This type of DM consists of (initially) condensed scalar particles with tiny masses, $\sim 10^{-22}$ eV, such that their wave-like nature becomes relevant at astrophysical scales. These particles are assumed to have no self-interactions, quite like axions [20], which could modify the picture if sufficiently strong [62]. In the absence of such interactions, the class of Fuzzy DM models is conveniently described by a single parameter, namely the DM mass $m_{22}\equiv m_{\psi}/10^{-22}$ eV, where $m_{\psi}$ denotes the actual physical particle mass. However, note that strong constraints exist on these scenarios, e.g. upper limits $m_{22}<1.5$ from the kinematics of dwarf galaxies [63, 64], which were recently complemented by strong lower limits from the abundance of high-$z$ galaxies [65], $m_{22}>10$, superseeding earlier limits from their luminosity functions ($m_{22}>1.2$ [66]) by nearly one order of magnitude. Even more recently, the Intergalactic Medium (IGM) have provided the tightest limits on the mass of $m_{22}>20$ (2$\sigma$ C.L.) for a very conservative analysis of high redshift data [21] (while $m_{22}>37.5$ is obtained for a less conservative scenario where some priors on the IGM thermal history are assumed).

In Fig. 6, we show a few example power spectra (left) and transfer functions (right), associated with different values of $m_{22}$. The point that we want to illustrate is that, even though the known oscillations are present in the fuzzy DM power spectra, their transfer functions are still well described by our general parametrisation from Eq. (2.4), simply because the oscillations appear only at the smallest scales or, equivalently, at large values of $k$. To achieve a good fit, we have simply cut off the oscillations, which are very suppressed and therefore negligible for most applications, including the MW satellite counting and the Lyman-$\alpha$ forest methods.

Using this approximation, we can see that our parametrisation from Eq. (2.4) in fact provides a very good fit to the fuzzy DM transfer functions, cf. right panel of Fig. 6. We will see in Sec. 6 how well this strategy truly works, when both the actual models and our fits are subjected to a reality test.

3.5 Effective theory of structure formation (ETHOS)

In order to further demonstrate the variability of our approach, we should compare it to that of ETHOS [32, 33], which consists of an attempt to formulate an effective theory of cosmic structure formation, to map virtually any particle physics model to the constraints from astrophysics and cosmology. While to some extent this work is done in the same spirit as ours, we should point out the differences in our approach:

Langrangian-based vs. transfer function-based approach:

While we are simply trying to parametrise the transfer function for nCDM, cf. Eq. (2.4), the approach of ETHOS is to start directly from the particle physics Lagrangian. While at first this seems like a clear disadvantage of our strategy, since a particle physicist would need to compute the matter power spectrum (e.g. using CLASS [61]) before being able to apply our results to their model, we would like to point out that the mapping of DM models into the quantities relevant for cosmic structure formation is injective. Put in simpler words, many models that look quite different from the particle physics point of view will yield precisely the same predictions for structure formation (this can be seen easily for thermal examples, e.g., by rescaling both temperature and DM mass). With our approach it is therefore possible to summarise most DM models in a much simpler framework. In fact, it may be a too big effort to start from any possible particle physics Lagrangian, when the key point for comparison lies in the transfer function.

We thus consider our approach to be the most economic. 2. 2.

Validity for small scales, i.e., for large $k$:

As we have pointed out e.g. in Sec. 3.4, our approach cannot capture the very smallest scales, which are suppressed in the transfer function due to dividing by the CDM power spectrum. However, given that this is in any case the part of the spectrum with the smallest power, we do not expect big effects of this region, unless we find an observable that is truly sensitive to very small scales. The ones we are using (i.e., satellite counts and Lyman-$\alpha$) are not sensitive to that extreme region. For instance, as visible in the left panel of Fig. 7, oscillations in the transfer functions seem very prominent when plotted in log-scales, although they are in fact unimportant for most aspects of structure formation.

We thus consider our approach to be safe as long as the regions for (very) large $k$ play no role. 3. 3.

Model coverage:

By construction, our approach of fitting the transfer function is much simpler and therefore less versatile than a method starting from the particle Lagrangian. However, we want to point that the original fitting function of Eq. 2.4 can be easily generalised if necessary. For example, a plateau in the transfer function can be described by adding one more parameter $\delta$, i.e.

[TABLE]

In principle, oscillationary patterns could also be included, for example by simply adding a cosine function of the form

[TABLE]

where $\delta$ is an additional free parameter.

We thus consider our approach to be easily extendable to basically cover the same range of models as ETHOS does.

Given that, we can try to fit some of the transfer functions obtained by ETHOS, and for this purpose we take the ones given in Fig. 1a of Ref. [32] as example (see that reference for the details on the data chosen).444Note that Ref. [32] uses an alternative definition of the transfer function of $T_{\rm ETHOS}(k)\equiv P(k)_{\rm{nCDM}}/P(k)_{\rm{CDM}}$, instead of the more common definition of $T^{2}(k)\equiv P(k)_{\rm{nCDM}}/P(k)_{\rm{CDM}}$, which is the one we are using. This difference in definitions is what creates the seeming difference between the right panel of our Fig. 7 and Fig. 1a of Ref. [32]. This will allow us to investigate whether our approach yields a similar result as the more detailed ETHOS analysis. We have again fitted the transfer functions with Eq. (2.4), however, note that this time we have only fitted the part for small $k$, i.e., left of the first “oscillation”. In Sec. 6 we will show that regarding the number of Milky-Way satellites as well as the power spectrum from the Lyman-$\alpha$ forest, there is hardly any difference between our simplifying fit without oscillations and the full transfer function from ETHOS.

4 Results from ${\boldsymbol{N}}$-body simulations

In order to perform cosmological $N$-body simulations and to put constraints on the nCDM models listed in Tab. 1, we have modified the numerical code 2LPT [68] – which generates initial conditions for running cosmological simulations – by implementing the new, general transfer function: now the code takes as inputs $\{\alpha,\beta,\gamma\}$ instead of the thermal WDM mass, and it computes the corresponding $T(k)$ from Eq. (2.4).

In the right panel of Fig. 2 we plot the matter power spectra associated with the combinations of $\{\alpha,\beta,\gamma\}$ listed in Tab. 1, computed at redshift $z=99$ by using the modified version of 2LPT, plus the matter power spectrum associated to a purely CDM universe (the dashed line in the plot). The dotted line refers to the matter power spectrum computed, in linear theory, with the numerical Boltzmann solver CLASS [61] in a $\Lambda$CDM universe at $z=99$ as well. All the relevant cosmological parameters are set to their Planck values [69]: $\Omega_{m}=0.317$, $\Omega_{\Lambda}=0.683$, $h=0.67$, $\sigma_{8}=0.8$.

We have used these snapshots as initial conditions for running 55+1 DM-only simulations (5123 particles in a 20 Mpc$/h$ box) with the GADGET-3 code, a modified version of the publicly available GADGET-2 code [70, 71]. On top of these simulations, we have used a Friends-of-Friends (FoF) algorithm [72] with the standard linking length $b=0.2$, in order to extract the DM halos. We also run the SUBFIND code [73] for searching for the substructures bound to each main FoF group.

The suite of simulations used here has the following purposes: $i)$ present a first assessment of non-linearities for the models discussed; $ii)$ address quantitatively how the mass functions based on linear theory predictions compare to the actual results of the $N$-body simulations; $iii)$ assess whether non-linearities in the matter power spectrum could affect the conclusions regarding the derived constraints from our approximate Lyman-$\alpha$ method. In the following, we will see that the constraints given are primarily based on the linear theory and the simulations are mainly used as a first cross-check that the results are indeed not altered by non-linearities and/or by a more accurate modeling of either the number of MW satellites or the IGM structures.

4.1 Non-linear matter power spectra

In Fig. 8 we plot the ratios between the nCDM non-linear matter power spectra with respect to the $\Lambda$CDM power spectrum. We show the power spectra of all nCDM models at redshift 5 (left panel) and redshift 2 (right panel). Additionally, we illustrate the thermal WDM cases with 2 and 4 keV, which are in agreement with a similar study made by Ref. [40]. The differences between the models gradually decrease when going to smaller redshifts. Below redshift 2, the small-scale power enhancement from the non-linear structure evolution starts to dominate the signal from different nCDM models, resulting in a progressive shift of the corresponding half-mode scales towards larger values of $k$. This evolution is the reason why the Lyman-alpha data from the highest redshift bins provide the strongest limits on the nature of DM.

4.2 Halo mass functions

In Fig. 9 we plot the ratios between the nCDM halo mass functions with respect to the $\Lambda$CDM halo mass function, at redshift 1 (left) and 0 (right). Note that all the nCDM models predict a lower abundance of halos with masses between $10^{8}$ and $10^{9}$ $M_{\odot}/h$: for some of the $\{\alpha,\beta,\gamma\}$-triplets the number of halos included in this mass range is expected to be 80$\%$ smaller with respect to $\Lambda$CDM predictions.

The visible upturn at low masses is not physical, but a consequence of the phenomenon of artificial clumping, due the limited resolution of our simulations, which affects models with suppressed initial power spectra [74, 75, 76]. This is why it is crucial to have good theoretical predictions for the low-mass behaviour of the halo (and subhalo) mass functions. This issue, essential for the analyses presented in Sec. 5.1, is discussed in more detail in the Appendix B. The oscillatory pattern which characterises the region corresponding to masses $\gtrsim 10^{10}~{}M_{\odot}/h$ is due to the cosmic variance, since the size of our simulated box does not allow to have enough large halos to have statistically fully meaningful results.

5 Comparison with observations

In this section we present the first constraints on $\{\alpha,\beta,\gamma\}$ from structure formation data, based on linear theory. We have constrained the parameters by using MW satellite counts and Lyman-$\alpha$ forest data, which constitute two powerful independent methods for testing the “non-coldness” of DM. Satellite counts rely on the simple fact that any nCDM model must predict a number of substructures within the MW virial radius not smaller than the actual number of MW satellites that we observe. The Lyman-$\alpha$ forest data analysis instead provides information on the matter power spectrum along our line of sight, down to scales corresponding to $k\sim 10~{}h/$Mpc, currently giving the most stringent limits on the masses of thermal WDM candidates, which is why they can also be expected to be strong for non-thermal cases.

5.1 Constraints from Milky Way satellite counts

Assuming the standard $\Lambda$CDM model, cosmological $N$-body simulations predict too many dwarf galaxies within the MW virial radius, with respect to observations. Therefore, MW satellite counts represent a very useful tool for constraining DM properties (see, e.g., Refs. [77, 78, 79, 80]). We can estimate the number of MW satellites following the approach of [81, 14], i.e., by multiplying the 15 ultra-faint satellites from SDSS by a factor 3.5, in order to account for the limited sky coverage of the survey, and by finally summing the 11 MW classical satellites. We obtain $N_{\rm sat}=63$ as an estimate of the number of observed satellites within the MW viral radius.

We can now compare $N_{\rm sat}$ with the number of subhalos $N_{\rm sub}$ predicted by our models, in order to extract some constraints on $\{\alpha,\beta,\gamma\}$. According to [14, 76], we can use the following formula for estimating the number of substructures for a given model:

[TABLE]

where $M_{\rm sub}$ and $S_{\rm sub}$ are the mass and the variance of a given subhalo, $M_{\rm halo}$ and $S_{\rm halo}$ are the mass and the variance of the main halo, defined as follows:

[TABLE]

with $P(k)$ being the linear power spectrum of a given model, computed at redshift $z=0$. In Appendix B we explicitly show that the mass functions for the grid of simulations performed in the present work are in good agreement with the theoretical mass function formalism outlined above.

Under the assumption of a MW halo mass $M_{\rm halo}=1.7\cdot 10^{12}~{}M_{\odot}/h$ [82] and by considering subhalos with masses $M_{\rm sub}\geq 10^{8}~{}M_{\odot}/h$, we can obtain the number of subhalos $N_{\rm sub}$ predicted by our models, by simply integrating Eq. (5.1). The results are reported in Tab. 2, where the models highlighted in bold-face are those in agreement with the number of observed satellites, i.e. with $N_{\rm sub}\geq 63$.

In Fig. 10 we show a 3-dimensional contour plot in the $\{\alpha,\beta,\gamma\}$-space, where each triplet is associated with a different model: the orange contour represents the upper bound on the region of the parameter space which is in agreement with MW satellite counts, according to the method that we have just outlined. Hence, in the right panel of Fig. 10, all the $\{\alpha,\beta,\gamma\}$-combinations which sample the orange volume correspond to models that predict a number of substructures at least equal to $N_{\rm sat}=63$. In the left panel, instead, we plot the allowed volume of the parameter space whether we require the number of subhalos predicted by each nCDM model to be equal or larger with respect to a more conservative estimate for the number of MW satellites, $N_{\rm sat}=57$. This number has been chosen in order to account for a ten percent sampling variance in the number of satellites at a given MW halo mass.

By marginalising over $\beta$ and $\gamma$ we obtain the following limits on $\alpha$:

[TABLE]

which would correspond, in the old one-to-one parametrisation, to a thermal WDM particle with mass $m_{\rm{WDM}}\approx 2~{}\rm{keV}$ (see Eq. (2.2) and Eq. (2.3)). These limits are less constraining than the latest constraints from structure formation data: as expected, modeling the power suppression with three free parameters weakens the constraints on $k_{1/2}$. Within our general approach, due to the dependence of $\alpha$ on $\beta$ and $\gamma$, even lighter DM candidates may provide suppressed power spectra in agreement with MW satellite counts.

Looking at Eq. (5.3), it may seem surprising that the constraint on $\alpha$ strengthens when we use a weaker rejection criterion for the nCDM models, i.e. a smaller value for $N_{\rm sat}$. The reason is that, by relaxing the limit on the number of substructures, we accept a larger number of $\{\alpha,\beta,\gamma\}$-triplets characterised by very small values of $\alpha$. This is a straightforward consequence of the shape of the volume of the parameter space shown in Fig. 10. By accepting all those models which predict $57\leq N_{\rm sub}<63$, we are slightly shifting towards higher values of $\alpha$ the whole orange 3-dimensional surface in the plot. As one can see by comparing the two panels of Fig. 10, due to the geometry of the orange contour, this shift mainly concerns $\{\alpha,\beta,\gamma\}$-combinations with $\alpha$ very close to zero. Therefore, by imposing a smaller value for $N_{\rm sat}$ and marginalising over $\beta$ and $\gamma$, small values of $\alpha$ increase their contribution to its probability distribution with respect to the high-value tail of the distribution, which instead is only minimally affected by the choice of a more conservative value for $N_{\rm sat}$. We are therefore pushing the peak of the probability distribution of $\alpha$ towards zero, yielding to a stronger upper bound on it. However, at this approximate level of analysis, the difference between the two bounds reported in Eq. (5.3) is practically negligible, as visible in Fig. 10.

5.2 Constraints from the Lyman-$\alpha$ forest data

The Lyman-$\alpha$ forest (namely the Lyman-$\alpha$ absorption produced by intervening intergalactic neutral hydrogen in the spectra of distant quasars) provides a powerful tool for constraining small-scale DM properties. It is indeed well established that the Lyman-$\alpha$ absorption is produced by the inhomogeneous distribution of the Intergalactic Medium (IGM) along different line of sights to distant bright sources [83]. Therefore, the Lyman-$\alpha$ forest has been used widely as a probe of the matter power spectrum on scales $0.5~{}{\rm Mpc}/h\,<\lambda<100~{}{\rm Mpc}/h\,$ [34, 35, 37].

Instead of computing absolute bounds, which can only be obtained through a full statistical analysis, our goal is to investigate the deviations of our nCDM models with respect to a thermal WDM reference model, i.e., $m_{\rm WDM}=3.5$ keV, which is one of the most updated constraints on WDM candidates (at 95% C.L.), obtained through a recent comprehensive analysis of Lyman-$\alpha$ forest data [37].

In order to constrain our models with Lyman-$\alpha$ forest data, we slightly modify the method developed in Ref. [14]. We parametrise the deviation of a model with respect to $\Lambda$CDM by the ratio

[TABLE]

where $P_{1\rm{D}}(k)$ is the 1D power spectrum of the model that we are considering, computed by the following integral on the 3D matter power spectrum:

[TABLE]

where $P(k^{\prime})$ is the 3D linear matter power spectrum, computed at redshift $z=0$.

We are now able to determine whether a model deviates more or less from $\Lambda$CDM, with respect to the thermal WDM reference model that we have chosen, by adopting the following criterion: a model is excluded (at 95% C.L.) if it is characterised by a larger power suppression with respect to the reference model. In order to quantify the suppression in the power spectra, we define the following estimator:

[TABLE]

where $A$ is the integral of $r(k)$ over the typical range of scales probed by Lyman-$\alpha$ observations ($0.5~{}h/{\rm Mpc}<k<20~{}h/{\rm Mpc}$ for the MIKE/HIRES+XQ-100 combined dataset used in [37]), i.e.,

[TABLE]

such that $A_{\rm\Lambda CDM}\equiv k_{\rm max}-k_{\rm min}$ by construction.

Analogously, by plugging the power spectrum of the thermal WDM reference model into Eqs. (5.4) and (5.6), we find $\delta A_{\rm REF,1}=0.38$, which is an estimate of the small-scale power suppression with respect to $\Lambda$CDM for models that are excluded at 95% C.L. by the Lyman-$\alpha$ forest data. In Tab. 3 we list the 55 models that we have studied, each of them with its corresponding $\delta A$: a given model is excluded (at 95% C.L.) if $\delta A>\delta A_{\rm REF,1}$, i.e., if it shows a power suppression $\gtrsim 38\%$ with respect to the $\Lambda$CDM power spectrum. Accepted models are highlighted in bold-face.

Let us now stress that the constraint on the thermal WDM mass associated to $\delta A_{\rm REF,1}$ (i.e. $m_{\rm WDM}=3.5$ keV) has been obtained under very conservative assumptions on the thermal history of the universe (see [37] for details). By modifying these assumptions, the lower limit on thermal WDM masses strengthens to $m_{\rm WDM}=5.3$ keV (at 95% C.L.), which represents the tightest bound from Lyman-$\alpha$ forest data up to date. By taking this limit as reference, we find indeed a corresponding small-scale suppression $\delta A_{\rm REF,2}=0.21$, with respect to the $\Lambda$CDM power spectrum.

We note that the physical observable for the Lyman-$\alpha$ forest data is the flux power spectrum $P_{\rm{F}}(k,z)$ and not the 1D or 3D linear matter power. However, two different key aspects of the Lyman-$\alpha$ forest physics suggest that the analysis presented could be also quantitatively correct. Firstly, the relation between linear matter and flux power can be modelled as $P_{\rm{F}}=b^{2}(k)P_{3\rm{D}}(k)$, with a bias factor $b^{2}(k)$ which differs only very little between $\Lambda$CDM and $\Lambda$nCDM models, at least for models reasonably close to the standard case (see e.g. [84, 24]): this motivates the application of Eq. (5.4) to flux power spectra as well. Secondly, the area criterion is justified by the fact that IGM peculiar velocities (typically $<100$ km/s) tend to redistribute the small-scale power within a relatively wide range of wave-numbers in the probed region of the flux power [85]. Thereby, the derived bounds should be robustly checked against full hydrodynamic simulations which can provide a forward modeling of the flux power spectrum, as done in [37]. We thus regard the approximate method presented here (based on the linear theory) as a first quantitative step towards a more comprehensive analysis.

In Fig. 11 we show a 3-dimensional contour plot in the $\{\alpha,\beta,\gamma\}$-space, which represents the region of the parameter space that contains models in agreement with Lyman-$\alpha$ forest data. The left panel refers to the conservative analysis, whereas the right panel refers to the “non-conservative” case. The red contours represent the 68% C.L. limit on the $\{\alpha,\beta,\gamma\}$-combinations, while the blue and green contours represent the 95% and 98% C.L. limits, respectively. All those models associated to $\{\alpha,\beta,\gamma\}$-triplets placed outside of the 3-dimensional coloured region are therefore excluded at 98% C.L. by our analysis. By marginalising over $\beta$ and $\gamma$ we obtain the following limits on $\alpha$:

[TABLE]

These limits would correspond, in the old one-to-one parametrisation, to a thermal WDM particle with a mass of $m_{\rm{WDM}}\approx 3$ keV, see Eqs. (2.2) and (2.3). Even at this approximate level of analysis, the Lyman-$\alpha$ forest generally tend to provide more stringent constraints than MW satellite counts. Note that, however, this is not true for all models: e.g., the point nCDM35 is allowed by the Lyman-$\alpha$ forest, while it is excluded by satellite counts. This is because, after all, the two methods probe slightly different scales, and thus are in reality complementary when it comes to constraining DM models. Note that this is a qualitative difference of non-thermal settings compared to thermal WDM: indeed, using different methods to constrain DM models is paramount to obtain a clear picture of what is allowed and what is not. Still, the limit shown is weaker compared to the most updated constraints on thermal WDM masses. As we discussed before, this is due to the new general parametrisation of $T(k)$. With our approach, thanks to the mutual dependence among $\alpha$, $\beta$, and $\gamma$, it is possible to model nCDM scenarios with non-trivially suppressed power spectra. Therefore, models with shallower transfer functions may be found to be in agreement with Lyman-$\alpha$ forest data even if the corresponding DM candidate mass lies below the current constraints for thermal WDM masses, given that those constraints refer to a very specific shape of the power suppression (i.e., a very specific $\{\alpha,\beta,\gamma\}$-combination).

As a double-check, we applied our method to the non-linear power spectra extracted from the simulations, finding consistent results with respect to the linear analysis. All the models rejected when comparing their linear power spectra are also rejected when comparing the non-linear ones.

6 Constraints on particle physics models

The goal of this section is to compare predictions in terms of MW satellite counts and power in the range probed by the Lyman-$\alpha$ forest by using the fits and the true transfer functions of some particle physics models defined in Sec. 3, thus to put our fitting formula, Eq. (2.4), to the reality check.

Note that we do not aim to give a full account of the validity (or invalidity) of the different nCDM models presented in Sec. 3, since in any case a few example points will not be able to give us a clear answer. Instead, we would like to find out whether the $\{\alpha,\beta,\gamma\}$-fit to a certain DM setting would lead us to the same conclusion about its validity when confronted with halo counting and Lyman-$\alpha$ bounds, while we do not care very much at this stage about whether a certain point in the model parameter space is now marked as “excluded” or “allowed” – we only want to check whether our conclusion about the points under consideration changes when we look at the fits instead of looking at the actual points.

To do so, we depict in Tab. 4 first of all the fit parameters for the example model points discussed in Sec. 3, which should be rather simple to grasp. For example, the second line for RP corresponds to the green curves in Fig. 3, while the first line for SD corresponds to the blue curves in Fig. 4. For each point, we have computed the number of satellites for both the fit and the corresponding transfer function of the real model, $N^{\rm fit}_{\rm sub}$ ($N^{\rm true}_{\rm sub}$), as well as the Lyman-$\alpha$ estimator, $\delta A_{\rm fit}$ ($\delta A_{\rm true}$) – with the difference of the fit to the real point indicated by the percentages in square brackets – which are in both cases matched to the respective conservative and non-conservative observational constraints. In all cases, no matter if the resulting number corresponds to a real model or to a fitted point, we use bold-faced scripts/Italic scripts/Roman scripts to indicate that a certain number is in agreement with both the conservative and non-conservative bounds/only with the conservative bound/with none of the bounds. As can be seen from the distribution of bold-faced or Italic numbers, many of the example points shown here are not in agreement with the bounds. However, what we are interested in is whether or not the fitted points would have brought us to the same conclusion. Indeed, this is the case for the vast majority of cases. In fact, given that the predictions from the fits deviate from the real model predictions only by a few per cent at most, we would expect agreement of the conclusions drawn from both versions of the transfer function (i.e., fitted and “true”) in all cases up to a few per cent.555Two fails in fourty comparisons, i.e., an empirical failure rate of 5% seems to support this prediction rather well – even though, of course, we have not selected the examples shown completely arbitrarily, but rather we have picked them to somehow reflect some of the variation possible for the different production mechanisms.

Thus, except for a tiny amount of borderline cases, the fitted points always yield the same conclusion as the actual model points. This is the main result of our paper:

Our fitting formula reproduces the true results to a very high degree.

Hence, whenever is it desired to match a nCDM setting to observational data, there is no need to do the whole computation. Instead, it is perfectly sufficient to match the resulting transfer functions to our Eq. (2.4) and to check whether the fitted points are allowed – which is known from the present text already, and which can for sure be improved in future studies. In fact, one could in principle also constrain the parameter space of our fitting formula by elaborate full scale $N$-body simulations, which would to quite some extent get rid of the potential necessity to run a computationally expensive simulation for a quasi-infinite list of models. With such constraints available, our fit could be used to easily identify the interesting regions in the parameter spaces of many types of complicated particle physics models of non-thermal DM, by simply computing linear power spectra.

7 Conclusions

Various dark matter scenarios predict structure formation to be suppressed at small cosmological scales. This suppression can be of different strength and shape depending on the particle nature of the DM candidate. In the past, a lot of effort has gone into investigating the astrophysical consequences of thermal WDM (i.e. candidates with a Fermi-Dirac or Bose-Einstein momentum distribution), which is characterised by a strong suppression of power only depending on the WDM particle mass. While being a good example for a “non-cold” scenario, the case of thermal WDM is far too restrictive to give justice to the large variety of potential nCDM models.

In this paper we develop a new method to accurately capture the clustering of nCDM models, based on a generalised fit of the transfer function (i.e., the square root of the linear power spectrum relative to the case of pure cold DM). The fit, given by Eq. (2.4), is based on a set of 3 free parameters, regulating the scale and the shape of the power suppression in a very flexible way. Our approach is meant to provide a connection between DM model building and astrophysics. It is expected to work efficiently for most of the known nCDM models. We have explicitly shown that this approach covers the full parameter space of sterile neutrinos (both by resonant production and particle decays), mixed (cold and warm) models, fuzzy DM, and models suggested by ETHOS. Due to the mutual dependence among the three free parameters of the new fitting formula, we are now able to disentangle even tiny differences in the shape of the power suppression and thus to distinguish between models with power spectra suppressed at very similar scales, in order to test them individually with cosmic structure formation data.

We have performed a large suite of $N$-body DM-only simulations that conservatively bracket the suppression of power suggested by the constraints obtained up to date from the Lyman-$\alpha$ forest and we have extracted the non-linear matter power spectra and (sub)halo mass functions.

We have given the first astrophysical constraints on the three free parameters which characterise our method, by using two independent astrophysical observables: the number of MW satellites and the Lyman-$\alpha$ forest. From this preliminary analysis it appears that no significant suppression of power can be present at scales larger than $\sim 0.06~{}h^{-1}$ (comoving Mpc) in order to fit both the observables, regardless of the shape of the power spectrum at smaller scales. This result is valid under the framework of abrupt suppressions of power happening at or below the scales $k\sim 1\,h$/Mpc, those typically induced by nCDM models, rendering the more gentle neutrino-induced suppression at larger scales still possible.

Finally we have constrained some of the most popular nCDM scenarios, by fitting the corresponding transfer functions with the new general formula and then testing them with MW satellite counts and Lyman-$\alpha$ forest data, using the methods that we have outlined.

This work represents a first step towards a general and comprehensive modeling of small-scale departures from the standard CDM paradigm. Future developments will include a full statistical analysis of Lyman-$\alpha$ forest data, by performing hydrodynamical simulations in order to extract the flux power spectra for our nCDM scenarios and determine more accurate limits on $\{\alpha,\beta,\gamma\}$, and a weak lensing data analysis, which will provide another independent observable for constraining the parameter space. It is also expected that forthcoming instruments like DESI666http://desi.lbl.gov/ or the Euclid satellite777https://www.euclid-ec.org/ could help in constraining the parameters by exploiting the weak lensing signal and or the clustering of galaxies at small scales, even if it is likely that most of the progress will be made by combining the DESI quasar spectroscopic survey with higher resolution Lyman-$\alpha$ forest spectra. Another interesting possibility, featuring a new observational probe, will be to probe the global intensity mapping signal of the 21 cm transition produced by neutral hydrogen, which has the advantages of being at high redshift where non-linear evolution is less important [86].

Acknowledgments

We would like to thank Francis-Yan Cyr-Racine for providing the data used in Fig. 1a of Refs. [32, 33], and we are also grateful for Nicola Menci for useful discussions on the small-scale part of the transfer functions. 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). RM and MV are supported by the INFN PD51 grant INDARK. MV is also supported by the ERC-StG “cosmoIGM”. AS acknowledges support from the Ambizione grant of the Swiss National Science Foundation (project number PZ00P2$\_$161363). Simulations were performed on the Ulysses SISSA/ICTP supercomputer. MT acknowledges support by Studienstiftung des deutschen Volkes as well as support by the IMPRS-EPP.

Appendix A Quasi-degeneracy between $\boldsymbol{\alpha}$ and $\boldsymbol{\gamma}$

By considering the high-$k$ limit of Eq. (2.4), we can easily notice the quasi-degeneracy between the two parameters $\alpha$ and $\gamma$. For large $k$ (i.e., small scales), such that $(\alpha k)^{\beta}\gg 1$, we have indeed:

[TABLE]

If we now change $\alpha$ to a new value $\alpha\to\tilde{\alpha}=x\cdot\alpha$, where $x$ is some real number, we can absorb this changes into a $k$-dependent change of $\gamma$:

[TABLE]

where we have used the obvious identity $x=k^{\ln x/\ln k}$. Thus, in Eq. (A.1), we can trade a change in $\alpha$ for a $k$-dependent change in $\gamma$:

[TABLE]

which reproduces the effect of the change from $\alpha$ to $\tilde{\alpha}$. Now, this does not seem like a real degeneracy – and mathematically it is not, due to the $k$-dependence of $\tilde{\gamma}$. However, this dependence is only very weak (logarithmic), so that one can, in spite of the $k$-dependence, view $\tilde{\gamma}$ as approximately constant:

[TABLE]

This approximation will be even better if $|\log_{k}x|=|\ln x/\ln k|\ll 1$, due to this quantity only appearing in the sum in Eq. (A.3). Thus, indeed, to a very good approximation there is a degeneracy between $\alpha$ and $\gamma$: a change in one of them can be traded for a change in the other. Hence, by covering several different values of $\gamma$ in Fig. 1, we have information about many different $\alpha$’s at the same time.

Appendix B Comparing simulated halo mass functions with theoretical predictions

We have already highlighted in Sec. 4.2 the low-mass upturn in the halo mass functions, due to artificial clumping. To take this effect into account and subtract the corresponding numerical artefacts, we have estimated the number of subhalos predicted by our models by integrating over Eq. (5.1), which is the theoretical subhalo mass function derived by an extended Press-Schechter approach [87, 88] based on a sharp-$k$ window [89, 90]. This leads to the following expression [75, 76]:

[TABLE]

where $\overline{\rho}$ is the average density of the universe, $\nu=\delta_{c,0}^{2}/S(R)$ is the peak-height of the perturbations (at $z=0$), and $f(\nu)$ is obtained by the excursion-set approach [91]. The relation between the sharp-$k$ filter scale and mass and the variance $S(R)$ are defined in Eq. (5.2). In Fig. 12, we compare example theoretical halo mass function predicted in Eq. (B.1) to the halo mass function extracted from the corresponding $N$-body simulations. The former is plotted as a solid line, the latter as a dotted line. Consistenly with the trend outlined in Fig. 9, below masses of the order of $10^{9}~{}M_{\odot}/h$, artificial clumping induces the upturn in the simulated halo mass functions and hence the discrepancies with the theory which are manifest in Fig. 12, below masses of the order of $10^{9}~{}M_{\odot}/h$.