<figure><img src="image_1.png" /></figure>

# Fundamental Physics with the Hubble Frontier Fields: Constraining Dark Matter Models with the Abundance of Extremely Faint and Distant Galaxies

N. Menci<fnref n="1" />, A. Merle<fnref n="2" />, M. Totzauer<fnref n="2" />, A. Schneider<fnref n="3" />, A. Grazian<fnref n="1" />, M. Castellano<fnref n="1" />, and N. G. Sanchez<fnref n="4" />

## Abstract

We show that the measured abundance of ultra-faint lensed galaxies at $z \approx 6$ in the Hubble Frontier Fields (HFF) provides stringent constraints on the parameter space of (i) dark matter models based on keV sterile neutrinos; (ii) “fuzzy” wavelike dark matter models, based on Bose–Einstein condensates of ultra-light particles. For the case of sterile neutrinos, we consider two production mechanisms: resonant production through mixing with active neutrinos and the decay of scalar particles. For the former model, we derive constraints for the combination of sterile neutrino mass $m_\nu$ and mixing parameter $\sin^2(2\theta)$ which provide the tightest lower bounds on the mixing angle (and hence on the lepton asymmetry) derived so far by methods independent of baryonic physics. For the latter we compute the allowed combinations of the scalar mass, its coupling to the Higgs field, and the Yukawa coupling of scalar to sterile neutrinos. We compare our results to independent existing astrophysical bounds on sterile neutrinos in the same mass range. For the case of “fuzzy” dark matter, we show that the observed number density $\approx 1$/Mpc<fnref n="3" /> of high-redshift galaxies in the HFF sets a lower limit $m_\psi \ge 8 \cdot 10^{-22}$ eV (at the $3-\sigma$ confidence level) on the particle mass, a result that strongly disfavors wavelike bosonic dark matter as a viable model for structure formation. We discuss the impact on our results of uncertainties due to systematics in the selection of highly magnified, faint galaxies at high redshift.

**Key words:** dark matter – galaxies: abundances – galaxies: formation

## 1. Introduction

Understanding the nature of the dark matter (DM) component of the universe constitutes a key issue in fundamental physics and in cosmology. Over the last two decades, the formation and the growth of cosmic structures has progressively adopted the cold dark matter (CDM) paradigm as a baseline (Peebles 1982; Blumenthal et al. 1984). This envisages DM particles to be characterized by thermal velocities small enough to produce negligible free streaming on scales relevant to structure formation. Typically, this corresponds either to assuming DM particles to be massive ($m_X > 0.1$ GeV) or to be constituted by condensates of light axions (with mass $\sim 10^{-5}-10^{-1}$ eV). Basic motivations for such a scenario are its simple properties (the corresponding power spectrum has a self-similar power law behavior in the whole range of mass scales involved in galaxy formation) and the fact that particles with mass and cross sections characteristic of the weak scales (weakly interacting massive particles, WIMPs, with masses $m_X \sim 100$ GeV) produce approximately the correct abundance of DM when they freeze-out of equilibrium in the early universe. On the other hand, extensive studies of structure formation showed that the CDM model provides an excellent baseline to explain the properties of large-scale structures and of galaxies on a huge range of mass scales ranging from $M \approx 10^{16}$ to $M \approx 10^9 M_\odot$ (see, e.g., Diemand & Moore 2011, for a review).

However, to date, both direct (see, e.g., Aprile et al. 2012, 2016; Akerib et al. 2014) and indirect (see, e.g., Adriani et al. 2013; Ackermann et al. 2015) CDM detection experiments have failed to provide a definite confirmation of such a scenario. Also, no evidence for CDM candidates with mass $10^2-10^4$ GeV has been

found in experiments at the LHC (e.g., Aad et al. 2013), while experiments aimed to detect axions as DM components have produced no evidence in the explored portion of the parameter space (see Graham et al. 2015; Marsh 2016). On the structure formation side, several critical issues are affecting the CDM scenario at the mass scales of dwarf galaxies ($M \approx 10^7-10^9 M_\odot$). These are all connected to the excess of power in the CDM power spectrum at such scales compared to a variety of observations, and include: the excess DM density of the inner regions of dwarf galaxies compared to the observed cored profiles (see de Vega et al. 2014), the over-abundance of faint dwarfs around our Galaxy and in the Local Group (see, e.g., Lovell et al. 2012) as well as in the field (Maccio et al. 2012; Menci et al. 2012; Papastergis et al. 2015), the excess of massive satellite DM halos with virial velocities $V_{vir} \ge 20$ km s$^{-1}$ relative to the number of observed bright dwarf galaxies (Boylan-Kolchin et al. 2011), and—most of all—the over-prediction of the abundance of field dwarfs with $V_{vir} \approx 40-60$ km s$^{-1}$ (Klypin et al. 2015).

While a refined treatment of baryonic effects entering galaxy formation (in particular feedback from supernovae) has been proposed to solve or at least alleviate the above problems (see, e.g., Governato et al. 2012, 2015; Di Cintio et al. 2014), the combination of astrophysical issues with the lack of direct or indirect detection of candidate particles has stimulated interest toward different DM scenarios, characterized by power spectra with suppressed amplitude at small mass scales ($M \lesssim 10^8-10^9 M_\odot$) with respect to the CDM case. In fact, several groups have started to investigate galaxy formation in a number of alternative models, such as self-interacting DM (Rocha et al. 2013; Vogelsberger et al. 2014), decaying DM (Wang et al. 2014), late-forming DM (Agarwal et al. 2015), atomic

<footnote n="1">¹INAF—Osservatorio Astronomico di Roma, via di Frascati 33, I-00040 Monte Porzio Catone, Italy</footnote>

<footnote n="2">²Max-Plank-Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, D-80805 München, Germany</footnote>

<footnote n="3">³Institute for Astronomy, Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093, Zurich, Switzerland</footnote>

<footnote n="4">⁴LERMA, CNRS UMR 8112, 61, Observatoire de Paris PSL, Sorbonne Universités, UPMC Univ. Paris 6, 61 Avenue de l'Observatoire, F-75014 Paris, France
*Received 2016 November 18; revised 2016 December 22; accepted 2017 January 3; published 2017 February 8*</footnote>
DM (Cyr-Racine & Sigurdson 2013), and DM interacting with dark radiation (Buckley et al. 2014; Chu & Dasgupta 2014).

In this framework, a specific focus has been given to warm dark matter (WDM) scenarios, which assume DM to be composed by particles with masses $m_X$ in the keV range. Their larger thermal velocities (corresponding to larger free-streaming lengths) suppress structure formation at scales $M = 10^7-10^9 M_\odot$, depending on the exact value of $m_X$ (since a thermalized species has no memory of the details of its production). While WDM candidates may result from the freeze-out of particles initially in thermal equilibrium in the early universe (like, e.g., gravitinos, see Steffen 2006 for a review), a similar suppression at these scales can be obtained by a variety of models featuring particles in the keV mass range with *non-thermal* spectra, like sterile neutrinos. Note that, in the case of non-thermal spectra, the production mechanism is essential in determining the suppression of the power spectra with respect to CDM. The shape of the power spectra are—although even that only qualitatively—somewhat similar to thermal WDM cases only for specific regions of the parameter space of the assumed production model, since a non-thermal spectrum cannot be associated with a temperature straightforwardly. Thus, to provide accurate limits to the mass of a non-thermal candidate, a detailed exploration of the parameter space of the selected models has to be performed. In this work, we tackle this task focusing on models of keV sterile neutrinos, which have received particular interest in the literature in recent years (Adhikari et al. 2016). This is due to both solid fundamental physics motivations (right-handed neutrinos constitute a natural extension of the Standard Model to provide mass terms for active neutrinos, see Merle 2013) and to the fact that such particles constitute the simplest candidates (see, e.g., Abazajian 2014) for a DM interpretation of the potential X-ray line in stacked observations of galaxy clusters and in the Perseus cluster (Boyarsky et al. 2014; Bulbul et al. 2014). In fact, in the presence of a tiny admixtures $\sin^2(2\theta)$ with active neutrinos, the decay of sterile neutrinos can result in photon emission at energies close to $m_\nu/2$. The non-detection of such a line in galaxies, dwarf galaxies, or the Milky Way (see, e.g., Sekiya et al. 2015; Adhikari et al. 2016; Jeltema & Profumo 2016; Riemer-Sorensen 2016, for an extended discussion) yields effective *upper* limits on the mass of sterile neutrinos for each value of the mixing angle in the range $10^{-13} \le \sin^2(2\theta) \le 10^{-9}$ (Canetti et al. 2013; Ng et al. 2015; Tamura et al. 2015), which rule out sterile neutrino models based on the non-resonant production mechanism of Dodelson & Widrow (1994).<fnref n="5" />

However, general *lower* mass bounds (arising from phase space density of nearby dwarf galaxies) are rather loose, yielding a model-independent limit $m_\nu \ge 0.4$ keV (Boyarsky et al. 2009), and leave a major portion of the parameter space largely unconstrained, in scenarios in which sterile neutrinos are resonantly produced (RP models) in the presence of a lepton asymmetry $L$ in the background medium (as proposed by Shi & Fuller 1999). They can, however, be significantly tightened when the information available from the production mechanism is taken into account (Merle & Schneider 2015; Schneider 2016), thereby strongly pushing such scenarios.

Other settings, where sterile neutrinos are produced by decays of, e.g., scalar particles (see a few concrete examples in Shaposhnikov & Tkachev 2006; Kusenko 2006; Petraki & Kusenko 2008; Klasen & Yaguna 2013; Merle et al. 2014; Merle & Totzauer 2015; Adulpravitchai & Schmidt 2015; Shakya 2016; König et al. 2016), are largely unconstrained by present X-ray observations, since they could in principle operate even without active-sterile mixing. Note that the corresponding parent in such scalar decay (SD) models must themselves be coupled to the Standard Model, typically via the Higgs sector. Thus, the key free parameters in such settings are the Higgs portal coupling and the Yukawa coupling between the scalar and the sterile neutrinos.<fnref n="6" /> For more complicated spectra, they can even look entirely different from thermal cases and exhibit qualitatively new features such as more than one momentum scale, long tails, or similar.

Another proposed solution to the small-scale problems in galaxy formation is based on Bose condensates of ultra-light (pseudo) scalar field DM with mass $m_\psi \approx 10^{-22}$ eV (Hu et al. 2000; Marsh & Silk 2014; Schive et al. 2014), which may be in the form of axions arising in string theory (Arvanitaki et al. 2010) or other extensions of the Standard Model of particle physics (see Kim 1987 for a historic review). In these settings, DM particles can be described by a single coherent wave function with a single free parameter $m_\psi$, the DM particle mass. The quantum pressure arising from the de Broglie wavelength produces a steep suppression of the transfer function below the corresponding Jeans length (Khlopov et al. 1985; Hu et al. 2000), making this scenario a viable alternative solution to the small-scale problems in galaxy formations (see Woo & Chiueh 2009; Mielke & Perez 2009; Chavanis 2011; Marsh & Silk 2014; Bozek et al. 2015; Harko & Lobo 2015; Madarassy & Toth 2015; Martinez-Medina et al. 2015). In fact, the soliton solution allows for cored inner density profiles in dwarf galaxies, while the substructures in DM halos, arising from fine-scale, large-amplitude cellular interference, would yield a suppressed abundance of satellite compared to the CDM case (Du et al. 2016). For such models, often referred to as “fuzzy” DM, the parameter space corresponding to the different possible power spectra is rather simple, since it depends only on the particle mass $m_\psi$. This is expected to be $m_\psi \lesssim 1.2 \times 10^{-22}$ eV (see, e.g., Marsh & Pop 2015) to resolve the cusp-core problem (without recourse to baryon feedback, or other astrophysical effects), but values up to $m_\psi \approx 4 \times 10^{-22}$ eV have been considered in previous works (see Schive et al. 2016).

The abundance of low-mass cosmic structures provides an important key to either restrict the range of allowed sterile neutrino production models, or rule out currently allowed DM scenarios. In the case of thermal relic WDM particles, the one-to-one correspondence between the WDM particle mass and the suppression in the power spectrum at small scales has allowed us to derive limits on $m_X$ by comparing the predictions from N-body WDM simulations or semi-analytic models with the abundance of observed ultra-faint satellites. On this basis, different authors have derived limits ranging from $m_X \ge 1.5$ keV (Lovell et al. 2012) to

<footnote n="5">⁵ Note that, contrary to the claims tracing back to the early references (Dodelson & Widrow 1994; Colombi et al. 1996), the spectrum resulting from non-resonant production is *non-thermal* (Merle et al. 2016), rather than being proportional to a Fermi-Dirac distribution multiplied by a suppression factor.</footnote>

<footnote n="6">⁶ Again, non-thermal DM models can suppress so much power at large scales that they are ruled out by the observed cosmic structure. However, the term “hot dark matter” should be used with care also in this case, since a non-thermal spectrum cannot be associated with a temperature straightforwardly.</footnote>
$m_X \ge 1.8$ keV (Horiuchi et al. 2014), $m_X \ge 2$ keV (Kennedy et al. 2014), and $m_X \ge 2.3$ keV (Polisensky & Ricotti 2011); relevant constraints have also been obtained for the parameter space of RP sterile neutrino models (Schneider 2016), in that latter case also taking into account the actual shape of the distribution functions (and, thus, not using a WDM approximation to a non-thermal case). Note that, however, such a method is appreciably sensitive to the assumed completeness corrections (see the discussion in Abazajian et al. 2011; Schultz et al. 2014), to the treatment of sub-halo stripping, and to the assumed values for the DM mass of the host halo and of the satellites. At higher redshifts, $z \approx 6$, a limit $m_X \gtrsim 1$ keV has been derived from the UV luminosity functions of faint galaxies ($M_{UV} \approx -16$) in Schultz et al. (2014); a similar approach by Corasaniti et al. (2016) yields $m_X \gtrsim 1.5$ keV. The same method has been applied to fuzzy DM models in Schive et al. (2016) (see also Corasaniti et al. 2016), deriving a consistency of such a model with the observed galaxy abundances for the whole particle mass range, $10^{-22} \le m_\psi/eV \le 4 \times 10^{-22}$. Since these approaches are based on the comparison between the observed luminosity functions and the predicted mass function of DM halos in different WDM models, the delicate issue in these methods is their dependence on the physics of baryons, determining the mass-to-light ratio of faint galaxies. However, to a lesser extent, uncertainties in the baryonic physics also affect (Viel et al. 2013; Garzilli et al. 2015) the tighter constraints achieved so far $m_X \ge 3.3$ keV for WDM thermal relics, derived by comparing small scale structure in the Ly$\alpha$ forest of high-resolution ($z > 4$) quasar spectra with hydrodynamical N-body simulations (Viel et al. 2013); for a generalization of this method to sterile neutrinos models, see Schneider (2016).

To derive robust limits on alternative DM scenarios, it is crucial to bypass the uncertainties related to the physics of baryons involved in galaxy formation, in order to get rid of the degeneracies between the effect of baryons and of the DM power spectrum in suppressing the number of observed low-mass structures. With this aim, Menci et al. (2016a, 2016b) have exploited the downturn of the halo mass distribution $\phi(M, z)$ in models with suppressed power spectra, which yields a maximum number density $\bar{\phi}$ of DM halos in the cumulative mass distribution which in turn depends on the adopted DM model. Since luminous galaxies cannot outnumber DM halos, an observed galaxy density $\phi_{\text{obs}} > \bar{\phi}$ would rule out the adopted DM model independently of the baryonic processes determining the luminous properties of galaxies. Such a method, first applied to lensed galaxies at $z=10$ by Pacucci et al. (2013) and to galaxies at $z=7$ in the Hubble Deep Field by Lapi & Danese (2015), has acquired an increased potential with the first results of the Hubble Frontier Field (HFF) programme. By exploiting the magnification power of gravitational lensing produced by foreground clusters, HFF has enabled the detection of galaxies fainter than the detection limit of the Hubble Deep Field at $z \ge 6$ (see, e.g., Atek et al. 2015; Ishigaki et al. 2015; Laporte et al. 2015; Castellano et al. 2016), to reach unprecedented faint magnitudes $M_{UV} = -12.5$ in the measurement of the luminosity function of galaxies at $z=6$ (Livermore et al. 2016). The large number density $\phi_{\text{obs}} \ge 1.3$ Mpc$^{-3}$ (at the $2-\sigma$ confidence level) corresponding to the observed luminosity function has allowed us to set a robust lower limit $m_X \ge 2.5$ keV (at $2-\sigma$) to the mass of thermal relic WDM particles. This constitutes the tightest

constraint derived so far on thermal WDM candidates independent of the baryon physics involved in galaxy formation.

Given the potential and the robustness of this method, it is only natural to apply it to different sterile neutrino production models, which offer a natural framework for physically motivated DM candidates as discussed above. In this paper, by requiring the maximum number density attained in a given DM model to be larger than the observed HFF value, we will provide unprecedented constraints on the parameter space for different sterile neutrino production mechanisms. In particular, the combinations $m_\nu - \sin^2(2\theta)$ defining the different RP models will be explored to significantly restrict the allowed regions in combination with existing bounds, while for SD models we will set limits on the different combinations of free parameters, i.e., the scalar mass, the Higgs portal coupling, and the Yukawa coupling between the scalar and the sterile neutrinos. We shall then turn to apply the method to the fuzzy DM scenario to provide lower limits on the ultra-light pseudo-axion mass $m_\psi$ which will turn out to rule out this class of models.

The paper is organized as follows. In Section 2 we provide a description of the method. In Section 3 we briefly summarize the basic features of the DM models we are probing in this work, i.e., the RP and the SD production models for sterile neutrinos, and the fuzzy DM model. In Section 4 we present and discuss the results, while Section 5 is devoted to our conclusions.

## 2. Method

To derive constraints on the parameter space of the sterile neutrino (RP and SD) DM models and of the fuzzy DM model, we compute the maximum number density of DM halos $\bar{\phi}$ expected at redshift $z=6$ for each point of the parameter space, and compare it to the observed number density $\phi_{\text{obs}}$ of galaxies at the same redshift obtained from the galaxy luminosity function measured by Livermore et al. (2016). Since the observed galaxies cannot outnumber their host DM halos, the condition $\bar{\phi} \ge \phi_{\text{obs}}$ determines the set of parameters admitted for each DM model.

The method is based on the drop of the differential halo mass function $d\phi/dM$ at small masses in models where the power spectrum is strongly suppressed with respect to CDM at masses $M \lesssim 10^7-10^9 M_\odot$. As a consequence, the corresponding cumulative halo mass function $\int_M^\infty dm (d\phi/dm)$ saturates to a maximum value $\bar{\phi}$ when the integral is extended down to progressively smaller values of $M$. This provides a maximum value for the number density of DM halos regardless of the underlying mass-luminosity relation, and—hence—is completely independent of the baryon physics entering galaxy formation.

In the following, we describe our computation of the observed number density $\phi_{\text{obs}}$ and or the maximum predicted number density $\bar{\phi}$, in turn.

### 2.1. The Galaxy Number Density at $z=6$ from Observed Luminosity Functions

For the observed number density $\phi_{\text{obs}}$, we take the value derived integrating the galaxy luminosity function at $z=6$ down to the faintest bin $M_{UV} = -12.5$ (Livermore
et al. (2016). The luminosity function was estimated from objects in the Abell 2744 and MACS 0416 cluster fields, selected on the basis of their photometric redshift. The effective volume sampled by the observations was computed from simulations of $5 \times 10^5$ mock galaxies with realistic colors and shapes, using all lensing models available for the two fields in the MAST archive. The simulated sources are added to the images and analyzed with the same procedure used for real sources, including photometric redshift analysis, in order to estimate the flux-dependent selection completeness. The resulting fiducial luminosity function with the corresponding 1-$\sigma$ uncertainties in each magnitude bin is estimated on the basis of the median magnification for each galaxy in the sample and is reported in Figure 10 of Livermore et al. (2016). From this we have derived the observed cumulative number density $\phi_{\text{obs}}$ (and its confidence levels) through a Monte Carlo procedure (Menci et al. 2016a). We extracted random values $\Phi_{\text{random}}(M_{\text{UV}})$ of the luminosity function in each magnitude bin according to a Gaussian distribution with variance given by the error bar in Livermore et al. (2016). Thus, for each simulation we produced a new realization of the luminosity function at $z=6$. From this, a cumulative number density $\phi_{\text{random}}$ was derived by summing up the values of $\Phi_{\text{random}}(M_{\text{UV}})$ in all the observed magnitude bins in the range $-22.5 \le M_{\text{UV}} \le -12.5$. We carried out $N_{\text{sim}} = 10^7$ simulations to compute the probability distribution function (PDF) of the cumulative number density $\phi_{\text{random}}$. We obtain a median value $\log \phi_{\text{obs}}/\text{Mpc}^{-3} = 0.54$, while from the relevant percentiles of the PDF we derive lower bounds 0.26, 0.01, and $-0.32$ at the 1, 2, and 3-$\sigma$ confidence levels, respectively. We have checked that the median value of the differential luminosity function $\Phi_{\text{random}}$ obtained from our simulations is consistent (within 3%) with the best fit value of the luminosity function obtained by Livermore et al. (2016). Note that such a procedure assumes that each magnitude bin is uncorrelated with the adjacent ones (as indeed is done in both the Livermore 2016 and the Bouwens 2016 analysis). In case of correlated bins larger errorbars are expected for blank field observations (see the discussion in Castellano et al. 2010), but a quantitative estimate of such an effect for measurements involving lensing magnification is lacking; this would constitute an interesting improvement over the present treatment.

Note that the measurements of the luminosity functions derived by Livermore et al. (2016) are particularly delicate at the faint end where large lensing magnifications are involved. Indeed, Bouwens et al. (2016b) adopted a different estimate of the impact of lensing magnification, finding not only a lower median value for the number density of galaxies at $M_{\text{UV}} = -12.5$ compared to Livermore et al. (2016), but also larger errorbars. In fact, assuming the luminosity functions in Bouwens et al. (2016) we find for the maximum number density a median value $\log \phi_{\text{obs}}/\text{Mpc}^{-3} = -0.25$, with lower bounds $-0.47, -0.62$, and $-0.9$ (at the 1, 2, and 3-$\sigma$ confidence levels), yielding looser limits on the parameter space of DM models. We discuss the impact of assuming such values for $\phi_{\text{obs}}$ in Section 5. In the same section we present a critical discussion of the systematics associated with the measurements of highly magnified sources.

## 2.2. The Maximum Number Density of Halos in Dark Matter Models with Suppressed Power Spectra

The computation of the differential halo mass function $d\phi/d(\log M)$ in sterile neutrino models is based on the standard procedure described and tested against N-body simulations in Schneider et al. (2012, 2013), Benson et al. (2013), and Angulo et al. (2013). Our computation has been tested against simulations in Menci et al. (2016b) and presented in Menci et al. (2016a). Here we summarize the key points of the computation and we refer to the above papers for details.

The key quantity entering the computations is the variance of the linear power spectrum $P(k)$ of DM perturbations (in terms of the wave-number $k = 2\pi/r$). Its dependence on the spatial scale $r$ of perturbations is given by

$$ \frac{d(\ln \sigma^2)}{d(\ln r)} = -\frac{1}{2\pi^2 \sigma^2(r)} \frac{P(1/r)}{r^3}. \quad (1) $$

Here we have used a sharp-$k$ form (i.e., a top-hat sphere in Fourier space) for the window function $W(kr)$ relating the variance to the power spectrum $\sigma^2(M) = \int dk k^2 P(k)W(kr)/(2\pi^2)$—see e.g., Schneider et al. (2013) for other common choices. The normalization $c$ entering the relation between the halo mass $M = 4\pi \bar{\rho}(cr)^3/3$ and the filter scale $r$ must be calibrated through simulations (here, $\bar{\rho}$ is the background density of the universe). All studies in the literature yield values for $c$ in the range $c = 2.5-2.7$ (see, e.g., Angulo et al. 2013; Benson et al. 2013; Schneider et al. 2013). The effect of such an uncertainty will be considered in deriving the constraints presented in Section 4.

For the sterile neutrino RP and SD models, the power spectrum is computed directly by solving the Boltzmann equation after computing the distribution function for all points of the parameter space, as described in detail in Section 3. Then, the differential halo mass function (per unit log $M$) based on the extended Press & Schechter approach (Bardeen et al. 1986; Benson et al. 2013; Schneider et al. 2013) reads

$$ \frac{d\phi}{d(\ln M)} = \frac{1}{6} \frac{\bar{\rho}}{M} f(\nu) \frac{d(\ln \sigma^2)}{d(\ln r)}. \quad (2) $$

Here, $\nu \equiv \delta_c^2(t)/\sigma^2$ depends on the linearly extrapolated density for collapse in the spherical model, $\delta_c = 1.686/D(t)$, and $D(t)$ is the growth factor of DM perturbations. We conservatively assume a spherical collapse model, for which $f(\nu) = \sqrt{2\nu/\pi} \exp(-\nu/2)$. Assuming an ellipsoidal collapse model would yield a lower halo mass function at the low-mass end and—hence—even tighter constraints on the DM particle mass.

For sterile neutrino DM, yielding power spectra suppressed at small scales compared to CDM, the resulting differential mass functions, see Equation (2), are characterized by a maximum value at masses close to the “half-mode” mass (Schneider et al. 2012; Angulo et al. 2013; Benson et al. 2013), the mass scale at which the spectrum is suppressed by 1/2 compared to CDM. This function depends strongly on the sterile neutrino mass; for RP models it also depends on the lepton asymmetry assumed and, hence, on the resulting mixing angle $\theta$; typical power spectra in such models yield half-mode masses ranging from $M_{\text{hm}} \approx 10^{10} M_{\odot}$ to $M_{\text{hm}} \approx 10^8 M_{\odot}$. Correspondingly, the cumulative mass functions saturate to a
maximum value $\bar{\phi}(z) \approx \phi(M_{\text{hm}}, z)$, defining the maximum number density of DM halos associated with the considered power spectrum.

A similar behavior characterizes the halo mass functions in the fuzzy DM case. Dedicated N-body simulations (Schive et al. 2016) yield for the differential mass function the form

$$ \frac{d \phi}{d(\ln M)} = \frac{d \phi}{d(\ln M)} \left|_{\text{CDM}} \left[ 1 + \left( \frac{M}{M_0} \right)^{-1.1} \right]^{-1.2} \right. , \quad (3) $$

where $|d\phi/d(\ln M)|_{\text{CDM}}$ is the halo mass function in the CDM scenario computed after Equations (1) and (2), assuming for $P(k)$ the CDM form given in Bardeen et al. (1986). The auxiliary mass scale $M_0 = 1.6 \times 10^{10} (m_{\psi}/10^{-22} \text{eV})^{-4/3} M_\odot$, determining the suppression of the halo mass function compared to the CDM case, depends on the fuzzy DM candidate mass, and it plays a role analogous to the half-mode mass scale for sterile neutrino models.<fnref n="7" /> For $m_{\psi} \ge 1.5 \times 10^{-22} \text{eV}$ the uncertainties in the above expression for the halo mass functions are below 10% (Schive et al. 2016) and will be considered when comparing with observed galaxy number densities.

For each DM model considered, the method described above allows us to exclude the region of the parameter space yielding $\phi_{\text{obs}} > \bar{\phi}$. In the next section, we briefly recall the properties of the parameter space we explore for the different DM models we consider.

### 3. The Dark Matter Models: Power Spectra and the Parameter Space

#### 3.1. Sterile Neutrinos: Resonant Production from Mixing with Active Neutrinos

The minimal setup for sterile neutrino DM is the production via mixing with one or several active neutrino flavors. Active neutrinos are weakly interacting and are therefore in thermal equilibrium with other Standard Model particles in the early universe (i.e., at temperatures above the MeV range). During that epoch, the sterile neutrino abundance builds up gradually via occasional oscillations from the active to the sterile sector. Depending on mixing angle and particle mass, this freeze-in production can lead to the right amount of sterile neutrino DM (Dodelson & Widrow 1994) and is usually referred to as the Dodelson-Widrow (DW) mechanism or non-resonant production.<fnref n="8" /> Recent investigations based on combined limits from structure formation and X-ray observations have ruled out this production mechanism as the dominant contributor to the DM sector (Seljak et al. 2006; Viel et al. 2006; Horiuchi et al. 2014; Merle et al. 2016).

However, it was noticed early on that the active-sterile oscillation may be enhanced by a resonance (Shi & Fuller 1999), provided there exists a significant lepton asymmetry $L$ in the early universe. Such a resonance allows for significantly smaller mixing angles $\theta$, relaxing the tight limits from X-ray observations. Furthermore, the resulting particle momentum distributions may be colder and, therefore, in better agreement with structure formation (Abazajian

et al. 2001). The Shi & Fuller (SF) or resonant production mechanism of sterile neutrino DM (Shi & Fuller 1999), first mentioned in Enqvist et al. (1990), also plays an important part in the framework of the Neutrino Minimal Standard Model (νMSM), which attempts to simultaneously solve the problems of DM, non-zero neutrino masses, and baryon asymmetry by solely adding three additional sterile neutrino flavors to the Standard Model (see e.g., Asaka et al. 2005; Asaka & Shaposhnikov 2005; Canetti et al. 2013), at the cost of introducing a considerable fine tuning to produce a suitable lepton asymmetry. Note that in the νMSM, there is an upper limit for the allowed lepton asymmetry (see Canetti et al. 2013; Laine & Shaposhnikov 2008; Shaposhnikov 2008) and hence, there is a lower limit on the mixing angle even in the resonant case if the production mechanism is to explain the total DM abundance.

Since, for any given sterile neutrino mass, the mixing angle is related to the adopted lepton asymmetry $L$, in this work we describe the parameter space of RP sterile neutrino models in terms of combinations of sterile neutrino masses $m_\nu$ and mixing amplitudes $\sin^2(2\theta)$. Each one of such combinations corresponds to a different momentum distribution, which strongly differs from a generic Fermi-Dirac form (Abazajian et al. 2001). In many cases such non-thermal momentum distributions are characterized by colder mean particle momenta and by a larger range of different momenta compared to the case of non-resonant production. In terms of structure formation, this may result in density perturbations which are more gradually suppressed and survive down to smaller scales. This is, however, only true for certain parts of the parameter space, while other parts show similar (or even stronger) suppression of matter perturbations compared to the case of non-resonant sterile neutrino DM (for more details see Schneider 2016).

In this work the sterile neutrino momentum distributions of RP are computed with the public code sterile-dm of Venumadhav et al. (2016), which is an extension of earlier works (Abazajian et al. 2001; Kishimoto & Fuller 2008; Abazajian 2014). The computation is based on the Boltzmann equation and includes detailed calculations of the lepton asymmetry around the quark-hadron transition. Independent calculations by Ghiglieri & Laine (2015) give similar results. To obtain the power spectra, we use the publicly available Boltzmann solver CLASS (Blas et al. 2011; Lesgourgues & Tram 2011).

#### 3.2. Sterile Neutrinos: Production from SD

Production from SD is described by a generic model that invokes one real scalar singlet $S$ and (at least) one sterile neutrino $N$ beyond the Standard Model. The scalar singlet couples to the SM Higgs doublet $\Phi$ via a Higgs portal,

$$ \mathcal{L} \supset 2\lambda(\Phi^\dagger\Phi)S^2, \quad (4) $$

where $\lambda$ is a dimensionless coupling. The interaction between the scalar and the sterile neutrino is encoded in

$$ \mathcal{L} \supset -\frac{y}{2} S \overline{N^c} N, \quad (5) $$

where $y$ is a Yukawa-type coupling. In the most general case, the complete model Lagrangian should contain terms for the mixing between active and sterile neutrinos. We will, however, neglect these terms, since non-resonant active-sterile mixing cannot contribute significantly to the production of sterile

<footnote n="7">⁷ Instead of using the fitting function given by Schive et al. (2016), one could also apply Equation (2) to the case of fuzzy DM. We use Equation (3) for simplicity and because this guarantees full consistency with previous work on fuzzy DM.</footnote>

<footnote n="8">⁸ Note that the original proposal of this mechanism was prior to DW in Langacker (1989); however, DW were the first to link it to DM.</footnote>
neutrinos when X-ray limits and limits of structure formation
are taken into account (Merle et al. 2016). Note that the
assumption of having only one right-handed neutrino is not
very restrictive, since most statements can easily be generalized
to the case of several right-handed states (König et al. 2016).

The free parameters of the SD model are the following.

1. The Higgs portal coupling $\lambda$, which determines the production rate and the kinematics of the scalar from the SM degrees of freedom of the Higgs doublet.

2. The Yukawa coupling $y$, which enters into the decay rate of the scalar and hence controls how fast the scalar decays into sterile neutrinos.

3. The mass of the scalar singlet, $m_S$, which determines which channels contribute to the production of scalars and thereby finally to the abundance of sterile neutrinos. If $m_S < m_h/2$, where $m_h$ is the mass of the physical Higgs boson, thermal Higgses can directly decay into scalars which increases the production drastically compared to those cases where $S$ can only be produced from pairwise annihilation of SM degrees of freedom. For a detailed discussion of the individual regimes, see Section 2 of König et al. (2016).

4. The mass of the sterile neutrino, denoted $m_\nu$. This will strongly influence the effects on cosmological structure formation. Nonetheless we want to stress again that the SD mechanism is highly non-thermal and that no premature conclusion on the sterile neutrino being “too hot” or “quasi-cold” should be made.

In this work, we will treat $\lambda$, $y$, and $m_S$ as free parameters. For each triple of ($\lambda$, $y$, $m_S$), we fix the mass of the sterile neutrino by requiring it to reproduce the observed relic DM abundance. We use the best-fit values published by the Planck collaboration (Ade et al. 2016). We have scanned the remaining three-dimensional parameter space by first solving the homogeneous and isotropic Boltzmann equation governing the production of sterile neutrinos in the early universe (König et al. 2016). The resulting momentum distribution functions and the value of $m_\nu$ inferred from the relic abundance constraint then serves as an input to compute the linear power spectra, again using the CLASS code (Blas et al. 2011; Lesgourgues & Tram 2011).

Let us very briefly discuss the interplay between the Higgs portal and the Yukawa coupling, in order to make the results presented later on easier to digest: for small Higgs portal couplings, the scalar itself is produced by freeze-in and is always strongly suppressed compared to its would-be equilibrium abundance. Physically, this means that back-reactions of scalars into SM particles can be neglected completely. In this case, the relic abundance of sterile neutrinos (and hence the mass $m_\nu$) are independent of the Yukawa coupling $y$ for a fixed pair ($m_S, \lambda$). Nonetheless, the Yukawa coupling controls the production time of the sterile neutrinos, which is one of the key factors for structure formation. When $\lambda$ is large enough to equilibrate the scalars, they will be subject to the well-known dynamics of freeze-out. In this regime, sterile neutrinos can be produced from scalars in equilibrium and from those decaying after freeze-out. Accordingly, the number density of steriles and thereby their mass $m_\nu$ can strongly depend on $y$ even for fixed ($m_S, \lambda$). Again, we refer the reader to König et al. (2016) for a more detailed discussion.

### 3.3. Fuzzy DM

Fuzzy DM models assume the DM to be composed of a non-relativistic Bose-Einstein condensate, so that the uncertainty principle counters gravity below a Jeans scale corresponding to the de Broglie wavelength of the ground state. In this case, the suppression of small-scale structures and the formation of galactic cores in dwarf galaxies is in fact entirely due to the uncertainty principle, which counteracts gravity below the Jeans scale, corresponding to a mass scale $M_f = 10^7 M_\odot m_{22}^{-3/2}$ (Marsh & Silk 2014), where $m_{22} \equiv m_\psi/10^{-22}$ eV. In such models, the DM mass $m_\psi$ ultimately determines all the relevant DM physical scales in structure formation, since it determines the scale below which an increase in momentum opposes any attempt to confine the particle any further. For example, the inner halo density profiles in such models are well described by the stable soliton solution of the Schrödinger-Poisson equation, which extends up to a core radius of $r_\psi = 1.6 m_{22}^{-1} (M/10^9 M_\odot)^{-1/3}$ kpc (Schive et al. 2014); at larger radii, the properties of fuzzy DM halos are indistinguishable from CDM. Recent works comparing the observed stellar-kinematical data of dwarf spheroidal galaxies to the density profiles produced in fuzzy DM scenarios derived upper limits for $r_\psi$ which translate into constraints on the DM particle mass $m_\psi \lesssim 1.5 \cdot 10^{-22}$ (Marsh & Pop 2015; Calabrese & Spergel 2016).

Since the DM particle mass is the only free parameter in fuzzy DM models, comparing the mass distribution of collapsed DM halos derived by Schive et al. (2016), Equation (3), with the observed abundance of high redshift galaxies as described in Section 2, provides straightforward constraints on $m_\psi$ which can be compared with the existing bounds mentioned above.

## 4. Results

We now proceed to present the constraints on the parameter space of the different models presented above. In all cases we adopted the Planck values for the cosmological parameters (Ade et al. 2016).

### 4.1. Resonant Production of Sterile Neutrinos

In the case of resonantly produced sterile neutrino DM, we choose the free parameters to be the mass, $m_\nu$, and the mixing amplitude $\sin^2(2\theta)$. For each combination of such quantities, the lepton number $L$ is fixed to the value required to yield the right DM abundance. We first investigate the effect of varying the mixing angle for a fixed sterile neutrino mass by focusing on the case $m_\nu = 7.1$ keV, corresponding to a sterile neutrino whose decay could be at the origin of the potential 3.5 keV line in X-ray spectra of clusters. For such a case, the spectra computed as presented in Section 3 yield the cumulative halo mass functions shown in Figure 1 (left panel) for different values of $\sin^2(2\theta)$ and compared with the observed number density of galaxies with $M_{\text{UV}} \le -12.5$ in the HFF. Note that the lines in Figure 1 constitute upper limits with respect to the theoretical uncertainties discussed in Section 2. At small masses, they saturate to a maximum number density $\bar{\phi}$, which is plotted in the right panel of Figure 1 as a function of $\sin^2(2\theta)$. When compared with the observed number density of luminous galaxies from the HFF (the upper shaded areas), requiring the number density of DM halos to be larger than the observed abundance $\bar{\phi} \ge \phi_{\text{obs}}$ restricts the mixing angle in the range $2 \times 10^{-11} \le \sin^2(2\theta) \le 10^{-9}$ (at the $2-\sigma$ confidence level).
<figure><img src="image_2.png" /></figure>

**Figure 1.** Left panel: the cumulative mass functions computed at $z=6$ for RP sterile neutrino models with $m_ν = 7.1$ keV for different values of the mixing angle $θ$ shown by the labels on the right. The shaded areas correspond to the observed number density of galaxies with $M_{UV} \le -12.5$ within the $1-σ$, $2-σ$, and $3-σ$ confidence levels. Right panel: the maximum value $\bar{\phi}$ (including theoretical uncertainties) of the predicted number density of DM halos at $z=6$ for the case with $m_ν = 7.1$ keV as a function of the mixing angle $θ$. The upper shaded areas represent the observed number densities of galaxies with $M_{UV} \le -12.5$ within the $1-σ$, $2-σ$, and $3-σ$ confidence levels. The vertical filled area corresponds to the range of values of $\sin^2(2θ)$, consistent with the tentative line signal 3.5 keV in X-ray spectra (Boyarsky et al. 2014; Bulbul et al. 2014), with the hatched vertical area corresponding to present upper limits on $\sin^2(2θ)$ from the absence of such a line in the spectra of the Milky Way, and of dwarf galaxies from several authors, as given in Riemer-Sorensen (2016) and references therein.

When combined with existing upper bounds from X-ray observations (hatched region in the right-hand panel of Figure 1; see also Watson et al. 2012; Horiuchi et al. 2014), this restricts the range of allowed values to the small interval $2 \times 10^{-11} \le \sin^2(2θ) \le 10^{-10}$.

To constrain different combinations of the free parameters and to compare with previous results, we explore the whole range of free parameters using a grid of values for both $m_ν$ and $\sin^2(2θ)$. After computing the corresponding power spectra (Section 3), the condition $\bar{\phi} \ge \phi_{\text{obs}}$ leads to the exclusion region in the plane $m_ν-\sin^2(2θ)$ shown in Figure 2.

Note that the exclusion plot is characterized by two excluded regions. The upper region is bounded by the non-resonant DW limit $L=0$ (the upper green curve in Figure 2) and is characterized by non-thermal spectra with resonant peaks at low momenta in the momentum distribution. This leads to a strong suppression of the power spectrum yielding the upper exclusion region. For lower values of $\sin^2(2θ)$, the corresponding increase of the lepton asymmetry (see Section 3.1) shifts the amplitude and position of the resonance peaks toward larger momenta so that the overall spectrum becomes cooler. The corresponding predicted number density of galaxies $\bar{\phi}$ becomes large enough to be consistent with the observed values $\phi_{\text{obs}}$, and results in the allowed region of Figure 2. When the lepton number is further increased, i.e., for even smaller values of $\sin^2(2θ)$, the resonance peak in the momentum distribution shifts to sufficiently high momenta to make the spectrum warmer again, thus yielding the lower exclusion region.

We exclude all models with a sterile neutrino mass below $m_ν \le 5$ keV and also large parts of the parameter space above. The region around the 3.5 keV line ($m_ν = 7.1$ keV) is still allowed, but restricted to a narrow range of values for $-11.3 \le \log \sin^2(2θ) \le -10$. Our lower bounds on the mixing angle extend over a wide range of sterile neutrino

<figure><img src="image_3.png" /></figure>

**Figure 2.** Constraints on the RP sterile neutrino parameter space from our method are represented as exclusion regions, with $3-σ$ and $2-σ$ limits represented by darker and lighter colors. Our constraints are compared with upper bounds from X-ray observations (as reported in Riemer-Sorensen 2016) of the Milky Way and dwarf galaxies. Gray areas are excluded by current limits on the abundance of DM (Ade et al. 2016). The green line corresponds to the non-resonant DW case with vanishing lepton asymmetry $L=0$. We also show as a dashed line the constraint obtained in Schneider (2016) from the abundance of satellites in the Milky Way. The tentative line signal at 7.1 keV (Boyarsky et al. 2014; Bulbul et al. 2014) is shown by the point with errorbars.

masses up to $m_ν \approx 11$ keV, where the mixing angle is constrained to the interval $-12 \le \log \sin^2(2θ) \le -11.4$. Our results constitute the tightest lower bounds on the mixing angle (and hence on the lepton asymmetry) derived so far by methods independent of baryonic physics. For a large range of neutrino masses, our results are very close to the limits obtained in
Schneider (2016) from the local abundance of Milky Way satellites, using a counting based on the extended Press & Schechter approach, calibrated to match the total abundance of satellites measured in test cases through N-body simulations. The similarity of these results is encouraging, given the different method adopted and, most of all, the large difference in cosmic times between the observations used to constrain the models. However, it must be noted that estimates of the Milky Way satellites depend on the assumed upper limit for the DM halo of the Milky Way (3 · 10¹² M☉), on the assumed isotropic distribution of satellites, on the treatment of sub-halo stripping, and on the assumed lower limit for the DM mass satellites (10⁸ M☉; for more details on the approach, see Schneider 2015). Lovell et al. (2016) adopted a similar approach based on satellite abundances computed through a semi-analytic computation and found somewhat looser limits, leaving sin²(2θ) unconstrained above $m_ν \sim 8$ keV.

On the other hand, our results are less stringent compared to the method based on the small-scale structure in the Lyα forest of high-resolution (z > 4) quasar spectra with hydrodynamical N-body simulations (Viel et al. 2013). This method has been directly applied to warm, thermal DM (yielding a lower limit $m_X \ge 3.3$ keV), and later extended to sterile neutrinos (Schneider 2016) by considering the values and the slopes of the suppression of the one-dimensional power spectrum consistent with the observed Lyα structures. However, such a method depends on baryonic physics, and in particular on the thermal history of the intergalactic medium. Indeed, when assumptions concerning this quantity are relaxed, the Lyα limit on the thermal relic mass is reduced by about 1 keV (Viel et al. 2013; Garzilli et al. 2015), and becomes comparable to the constraint from the high-redshift galaxies (see Menci et al. 2016a).

## 4.2. SD Production of Sterile Neutrinos

As discussed in detail in Section 3, in this case the parameter space is three-dimensional, since it includes the mass of the scalar $m_S$, the Higgs portal coupling $\lambda$ and the Yukawa coupling with the scalar $y$. The sterile neutrino mass $m_ν$ is related to a combination of the three free parameters by demanding that the model reproduce the observed relic abundance of dark matter. Before performing a complete exploration of the parameter space, we first show a comparison between the model cumulative halo distributions and the observed number density of galaxies in the illustrative case of a sterile neutrino with $m_ν = 7.1$ keV (the candidate origin of the potential 3.5 keV line) in the limit of small Higgs portal coupling $\lambda \ll 10^{-6}$. In such a freeze-in limit (see Section 3), the scalar mass $m_S$ and the Yukawa coupling $y$ will be related for any given value of $m_ν$, thus reducing the parameter space to be explored. In Figure 3 (left panel) we show the cumulative halo distributions in such a regime for a scalar mass of $m_S = 60$ GeV and for different values of the Yukawa coupling $y$, and we compare them to the observed number density of high-redshift HFF galaxies. The strong dependence of the maximum predicted number density $\bar{\phi}$ on $y$ is highlighted in the right-hand panel of Figure 3, which shows that present data allow us to set a constraint $y \ge 9 \times 10^{-9}$ at the 2-σ confidence level.

We then extend our exploration to cover the whole parameter space of SD production models for sterile neutrinos. With this aim, we consider a grid of $\lambda$ and $y$ values for six different values

of the scalar mass $m_S$/GeV = 60, 65, 100, 170, 500, 1000. The latter are chosen to cover the three different regimes where different production channels are dominant for the production of scalars from their coupling to the SM particles (see König et al. 2016, and Section 3): the light scalar case $m_S \le m_h/2$, the intermediate case $m_h/2 \le m_S \le m_h$, and the heavy scalar case $m_S \ge m_h$ (in terms of the Higgs mass $m_h = 125$ GeV). These different regimes have a rather strong impact on the distribution function of the scalar which directly translates into that of the sterile neutrino and, hence, on the resulting power spectrum.

For each value of $m_S$, we compute the power spectrum corresponding to each point in the $\lambda$-$y$ plane as described in Section 3. Then we represent in Figure 4, left panel, the regions of the parameter space consistent with the galaxy number densities measured in the HFF ($\bar{\phi} \ge \phi_{\text{obs}}$). These regions clearly split into a freeze-out (for $\lambda \ge 10^{-6}$) and freeze-in (for $\lambda \ll 10^{-6}$) family. For the freeze-out family, decreasing the scalar mass $m_S$ leads to a tighter bound on $y$, while yielding an approximate lower bound of $\lambda \gtrsim 10^{-5.2}$ for the Higgs portal coupling. For the freeze-in family, decreasing the scalar mass $m_S$ pushes the admitted values of $\lambda$ to progressively smaller values, while providing progressively stronger limits on $y$. Note that, moving from $m_S = 65$ GeV to $m_S = 60$ GeV, the allowed freeze-in region shifts to appreciably smaller values of $y$ (actually out of the plot window). This is because a scalar $S$ with 60 GeV has a mass $m_S < m_h/2$. This opens up an entirely new production channel, namely Higgs decay $h \to SS$, which actually dominates the production rate. In this case the scalars (and hence the sterile neutrinos) are produced with a much larger number density once the Higgs decay sets in. Thus, for the same combination of $\lambda$-$y$, this larger number density requires a smaller sterile neutrino mass, yielding larger velocities for the DM particles and suppressing the abundance of low-mass halos.

Finally, the right panel of Figure 4 displays an explicit comparison of our high-z HFF bounds for the case $m_S = 100$ GeV with a version of Lyα bounds derived by König et al. (2016). The latter authors have introduced an approximate method of generalizing Lyα bounds to non-thermal cases, by computing the ratio of power spectra of the case under consideration to CDM and demanding the “upper” (i.e., high-power) part of the spectrum to be allowed when compared to a certain bound. Using the thermal relic mass bounds of 2.2 and 3.3 keV, the blue (consistent with both bounds), purple (only consistent with the first bound), and red (inconsistent with both bounds) regions displayed in the right panel of Figure 4 can be derived. Our HFF-bound, in turn, is displayed by the green line. As one can clearly see, the constraints from HFF galaxies basically track the Lyα bounds: although the latter are slightly stronger, they may possibly suffer from unknown systematics, given that the distribution of the intergalactic material is observed and not that of the DM. However, the fact that both types of bounds are basically tracking each other adds significant credibility to both constraints.

## 4.3. Fuzzy DM

The large observed number density of high-redshift galaxies turns out to provide particularly strong constraints on fuzzy DM, i.e., scenarios based on wavelike DM composed by ultra-light bosons. In this case, we compute the cumulative halo mass function directly from the formula in Equation (3) derived
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left panel: the cumulative mass functions computed at $z=6$ for SD sterile neutrino models with $m_S = 60$ GeV and $m_ν = 7.1$ keV, for different values of the Yukawa coupling $y$ (shown by the labels on the right) in the freeze-in limit. The shaded areas correspond to the observed number density of galaxies with $M_{UV} \le -12.5$ within the $1-\sigma$, $2-\sigma$, and $3-\sigma$ confidence levels. Right panel: the maximum value $\bar{\phi}$ (including the theoretical uncertainties) of the predicted number density of DM halos at $z=6$ for the case with $m_ν = 7.1$ keV as a function of the Yukawa coupling $y$. The upper shaded areas represent the observed number density of galaxies with $M_{UV} \le -12.5$ within the $1-\sigma$, $2-\sigma$, and $3-\sigma$ confidence levels.</figcaption></figure>

by Schive et al. (2016) from dedicated N-body simulations. This depends on the single free parameter constituted by the particle mass $m_ψ$, so that the exploration of the parameter space is straightforward. In the left panel of Figure 5 we show the cumulative halo mass function for different values of the DM particle mass (in units of $10^{-22}$ eV), where we have conservatively increased by 10% the model number densities derived from Equation (3) to account for theoretical uncertainties (see Schive et al. 2016). The strong suppression in the number of low-mass halos compared to the CDM case yields a lower limit $m_ψ \ge 10^{-21}$ eV for the DM particle mass (right panel of Figure 5) at the $3-\sigma$ confidence level. This tightens by a large factor the previous bounds on $m_ψ$, derived from the luminosity function of high-redshift galaxies in the Hubble Deep Field ($m_ψ \ge 1.2 \times 10^{-22}$ eV) given by Schive et al. (2016; see also Corasaniti et al. 2016). This is due to the increase (by more than one order of magnitude) in the number density of galaxies achieved by HFF observations of faint galaxies due to the magnification of the foreground clusters, which allowed to push the detection limit at $z=6$ from $M_{UV} = -15$ to $M_{UV} = -12.5$.

Our results constitute the tightest constraint on fuzzy DM particles derived so far, and have a strong impact for the whole class of models based on fuzzy DM. In fact, even allowing for observational and numerical uncertainties, all results in the literature indicate that the mass of fuzzy DM particles should be in the range $m_ψ = (1-5.6) \cdot 10^{-22}$ eV to provide solitonic cores matching the observed density profile of nearby dwarf galaxies (Lora et al. 2012; Schive et al. 2014; Marsh & Pop 2015; Calabrese & Spergel 2016; Gonzáles-Morales et al. 2016). This is inconsistent at more than the $3-\sigma$ confidence level with our lower limits, strongly disfavoring such a scenario as a viable solution to the cusp-core problem of dwarf galaxies.

## 5. Conclusions and Discussion

We have shown that the HFF measurements of the abundance of ultra-faint lensed galaxies at $z \approx 6$ down to faint magnitudes $M_{UV} \approx -12$ have profound implications on the nature of DM, in particular for DM models based on keV sterile neutrinos and on “fuzzy” DM models considered in the present work.

For the case of resonant production of sterile neutrinos, the high galaxy number densities $\sim 1/M_{\text{pc}}^3$ measured at $z=6$ by (Livermore et al. 2016) provide stringent limits on the combination of sterile neutrino mass $m_ν$ and mixing parameter $\sin^2(2\theta)$. In particular, our method provides lower limits on the values of $\sin^2(2\theta)$ for different values of $m_ν$ which are complementary to the upper limits derived by the non-observations of decay lines in the X-ray spectra of dwarf galaxies (Sekiya et al. 2015; Adhikari et al. 2016; Jeltema & Profumo 2016; Riemer-Sorensen 2016), thus restricting the allowed portion of the parameter space of such models to a narrow region; e.g., the region around the potential 3.5 keV line ($m_ν = 7.1$ keV) is still allowed, but restricted to a narrow range of values for $-11.4 \le \log \sin^2(2\theta) \le -10.2$. The constraints we obtain on $\sin^2(2\theta)$ for different values of $m_ν$ are close to those obtained from the abundance of low-mass satellites of the Milky Way (Schneider 2016) for a wide range of values $m_ν \lesssim 10$ keV. This supports the reliability of methods based on structure formation as probes for DM models with suppressed power spectra, especially considering that the two methods are applied to galaxies observed at cosmic times differing by $\sim 13$ Gyr. However, the approach in the present paper is not affected by the uncertainties affecting methods based on the abundance of Milky Way satellites, like the assumed mass of galaxy halos and sub-halos, or the assumed isotropic distribution of satellites. On the other hand, for large sterile neutrino masses of $m_ν \approx 10$ keV, the abundance of high-redshift galaxies seems to provide tighter
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left panel: the constraints on the parameter space of SD production model for sterile neutrinos. The colored regions represent the allowed region ($\bar{\phi} \ge \phi_{\text{obs}}$) in the log $\lambda$–log y plane for each value of the scalar mass $m_S$ (shown by different colors according to the labels on the right). For each scalar mass (and hence for each color) lighter and darker tonalities represent the 1-$\sigma$ and 2-$\sigma$ confidence levels, respectively. Right panel: for the case $m_S = 100$ GeV, we compare our bounds from HFF galaxies derived in the present work (green thick lines) with a version of $\text{Ly}\alpha$ bounds derived by König et al. (2016): the latter are represented as red (forbidden cases), blue (allowed cases), and purple (constrained cases) segments in the lines of constant $m_\nu$ values in the log $\lambda$–log y plane.</figcaption></figure>

limits on $\sin^2(2\theta)$ compared to the abundance of Milky Way satellites. This is related to the larger sensitivity of the former method to the shape of the power spectrum at large wavenumber $k$. In turn, this is due to the different dependence on the variance $\sigma^2(M)$ in the mass function in Equation (2) compared to the mass distribution of satellite sub-halos, and to the minimal mass assumed for satellite halos in the Milky Way ($M \ge 10^8 M_\odot$), which excludes tiny halos as hosts of satellite galaxies, thus reducing the sensitivity of such a method with respect to very large $k$-modes. Regarding the limits from $\text{Ly}\alpha$ absorption lines of distant quasar spectra, the present constraints are less stringent than those obtained by Schneider (2016). It is however important to notice that our method is entirely independent of the baryon physics entering the $\text{Ly}\alpha$ method, which could affect the thermal state of the intergalactic medium (see Viel et al. 2013; Garzilli et al. 2015).

For the case of sterile neutrino production via SD, our method provides constraints on the combination of scalar mass $m_S$, Higgs portal coupling $\lambda$, and Yukawa coupling $y$, which are very close to those derived in König et al. (2016)—as shown in Figure 4 for the case $m_S = 100$ GeV. In the freeze-in limit (see Section 3) of small Higgs portal coupling $\lambda \ll 10^{-6}$, where the Yukawa coupling constitutes the leading quantity in the determining the model properties, this is due to the strong dependence of the predicted maximum halo number density $\bar{\phi}$ on $y$ (shown in the right-hand panel of Figure 3 for the case of sterile neutrino mass $m_\nu = 7$ keV). When compared to the observed high-redshift galaxy abundance this allows us to set a constraint $y \ge 9 \times 10^9$ at the 2-$\sigma$ confidence level.

As for the fuzzy DM, the results from our method have a relevant impact for the whole class of such models. In fact, extending the analysis by Schive et al. (2016) to compare with the recent measurements of galaxy abundance from the HFFs yields unprecedented limits on the mass of candidate DM particles of $m_\psi \ge 1.2 \times 10^{-21}$ eV at the 2-$\sigma$ confidence level and of $m_\psi \ge 8 \times 10^{-22}$ eV at 3-$\sigma$. These values are one order of magnitude larger than those required to yield core sizes of

dwarf galaxies large enough to match the observations (Marsh & Pop 2015). Indeed, such limits for the fuzzy DM mass would correspond to an upper limit $M_I \lesssim 3 \cdot 10^5 M_\odot$ (see Section 3.3) for the mass scale associated with the de Broglie wavelength of wavelike DM, pushing the scale where the “fuzzy” nature of such DM shows up to values much smaller than those involved in galaxy formation. This result is consistent with the independent constraints derived from the phase space density of local galaxies by de Vega & Sanchez (2014).

The bounds presented here are based solely on the statistical properties of DM halos. Indeed, tighter limits can be obtained by adopting a physical description of the baryonic processes involved in galaxy formation, as shown in Menci et al. (2016b). This is the approach taken, e.g., by Corasaniti et al. (2016) who compared the observed luminosity function to theoretical luminosity functions derived from the abundance matching technique applied to high-resolution N-body simulations. The uncertainties associated with baryon physics introduced in such approach can be suppressed by complementing the statistical analysis with the study of relations among galaxy properties that are specific to each DM scenario and, as such, carry additional information on the nature of DM that is potentially testable with observations.

While our method is robust and independent of the baryon physics entering galaxy formation, and does not rely on any estimate for the DM mass of the observed galaxies, the above results are based on the abundance of $z=6$ galaxies measured in the HFF by Livermore et al. (2016). The number density derived by those authors is based on 167 galaxies at $z \ge 6$, and it is thus robust from the statistical point of view; in any case, the statistical uncertainties have been considered in computing the confidence levels of our results. There are however subtle systematic effects related to the estimation of the survey volume, i.e., the variance of the lensing magnification maps of HFF clusters and the physical sizes of faint, high-$z$ galaxies which are strongly magnified by the cluster potential well.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left panel: the cumulative mass functions computed at $z=6$ for fuzzy DM models with different values of the wavelike DM particle mass $m_{22} = m_\psi/10^{-22}$ eV. The shaded areas correspond to the observed number density of galaxies with $M_{\text{UV}} \le -12.5$ within the $1-\sigma$, $2-\sigma$, and $3-\sigma$ confidence levels. Right panel: the maximum value $\bar{\phi}$ (including the theoretical uncertainties) of the predicted number density of DM halos $\phi$ at $z=6$ as a function of $m_{22}$. The upper shaded areas represent the observed number density of galaxies with $M_{\text{UV}} \le -12.5$ within the $1-\sigma$, $2-\sigma$, and $3-\sigma$ confidence levels.</figcaption></figure>

These effects have been critically analyzed by Bouwens et al. (2016b), who showed that they could result in luminosity functions flatter than those derived by Livermore et al. (2016) at $M_{\text{UV}} \ge -14$. As discussed in Section 2.1, the different estimates of the magnification uncertainties obtained by Bouwens et al. (2016) would yield looser limits on the parameter space of DM models. For example, an observational number density $\log \phi_{\text{obs}}/\text{Mpc}^{-3} = -0.62$ corresponding to the luminosity function by Bouwens et al. (at the $2-\sigma$ confidence level, see Section 2.1) would lead to $\log \sin^2(2\theta) \ge -10.2$ for $m_\nu = 5$ keV, to $\log \sin^2(2\theta) \ge -11.5$ for $m_\nu = 7$ keV, and to $\log \sin^2(2\theta) \ge -12.2$ for $m_\nu = 10$ keV in the case of resonantly produced sterile neutrinos, thus lowering the lower bounds in Figure 2. For the SD models of sterile neutrinos, the lower limit on $y$ in Figure 4 would be decreased to $y \ge 5 \times 10^{-8}$ (for the case $m_\nu = 7$ keV), while for the case of fuzzy DM, adopting the Bouwens et al. (2016) luminosity functions would yield $m_\psi \ge 5 \times 10^{-22}$ eV at the $2-\sigma$ confidence level.

Thus, on the observational side, the first step to improve the results presented in this paper consists of a deeper understanding of the systematics associated with the lensing observations of faint, high-redshift galaxies. As for the lensing magnification, with typical values $\mu > 10$ and as large as ~50–100 for the faintest galaxies, the analysis in Livermore et al. (2016) shows that large differences are found in the magnification estimates of individual galaxies when different lensing models are assumed. This is in agreement with the recent analysis by Meneghetti et al. (2016), who compared the performances of several lensing models on artificial images mimicking the depth and resolution of HFF data. They find that the largest uncertainties in the lens models are found near substructures and around the cluster critical lines, concluding that uncertainties in the magnification estimate are growing with the magnification value itself (30% of magnification uncertainties at $\mu \ge 10$). However, while the magnification of each individual galaxy could be subject to a large variance, the

aggregated information of the abundance of high-$z$ galaxies is less subject to such systematics. The analysis in Livermore et al. (2016) shows that the different lensing models introduce minimal changes in the best-fit luminosity functions, even at the faintest absolute magnitudes $M_{\text{UV}} = -12.5$ analyzed at $z \sim 6$, with a systematic uncertainty on the slope $\alpha$ of the luminosity functions at $z \sim 6-7$ below 2%. Although such an uncertainty would not change appreciably the results presented in this work, the proper procedure to be adopted in deriving the variance in the luminosity function due to the magnification is still a matter of debate (see Bouwens et al. 2016b).

The other systematic effect affecting the luminosity functions measured by Livermore et al. (2016) is related to the size distribution of high-$z$ galaxies. In fact, the detection of faint galaxies is strongly affected by their surface brightness distributions, with compact galaxies more easily detected compared to the extended, low surface brightness ones (Grazian et al. 2011). In the simulations carried out by Livermore et al. (2016), a normal distribution of half-light radii $r_h$ was assumed with a peak at 500 pc. While similar (or even higher) number densities are obtained for values $r_h \gtrsim 100$ pc, for smaller values a significant suppression of the faint-end logarithmic slope $\alpha$ of luminosity function is found (up to 10% in the case $r_h = 40$ pc, Bouwens et al. 2016a). While Bouwens et al. find that the typical values of $r_h$ at $z=6$ are of the order of ~25–80 pc for $M_{\text{UV}} \ge -15$, larger values are usually obtained in the literature: extrapolating the size distribution obtained by Kawamata et al. (2015) for $M_{\text{UV}} \sim -20$ galaxies at $z \sim 6$ to fainter magnitudes (assuming a size–luminosity relation of $r_h \propto L^{0.5}$ Grazian et al. 2012), the typical sizes of galaxies with $M_{\text{UV}} = -12.5$ range from ~40 pc to ~100 pc with a log-normal distribution. Laporte et al. (2016) analyzed the size distribution of $z \sim 7$ Lyman-break galaxies, finding a typical value of 250 pc for galaxies with $-19 \le M_{\text{UV}} \le -17$. Thus, while the completeness correction adopted by Livermore et al. (2016) is consistent with most of the extrapolations in the literature, assessing the actual size distribution of ultra-faint
galaxies at $z \approx 6$ would constitute an important improvement for a precise determination of the number density of such galaxies, so as to lower this source of systematics below the present value of ≈10%. In the future we plan to analyze the impact of the size-luminosity relation on the HFF luminosity function based on detailed simulations, using the same technique adopted in Grazian et al. (2011, 2012) to comply with the delicate issues (i.e., distortions due to the strong lensing amplification) involved in the measurement of sizes for such faint, noisy galaxies.

Thus, although the measurement by Livermore et al. (2016) constitutes a state-of-the-art achievement, the analysis of the HFF observations is open to several advancements. This is true not only for the most magnified objects (as discussed above), but also for those moderate magnifications ($\mu \lesssim 10$), for which the largest source of uncertainty is constituted by the accuracy of photometric redshifts. In this respect, improvements can be readily obtained through the combination of different procedures (e.g., Castellano et al. 2016), an approach which is known to reduce systematic effects of photo-z estimates (Dahlen et al. 2013). Indeed, deep MUSE observations (Vanzella et al. 2016) are already able to provide the spectroscopic redshifts for faint $z \ge 6$ galaxies in the HFF pointings, improving both the lensing magnification maps and the photometric redshifts. Finally, the analysis in Livermore et al. (2016) is based at present only on the first two fields of the HFF survey: the inclusion of the remaining four strong lensing clusters (Lotz et al. 2016) will reduce both statistical uncertainties and mitigate possible cosmic variance effects. In addition, it is not to be disregarded that the six available parallel HFF pointings at a depth comparable to the HUDF will improve the determination of the logarithmic faint-end slope of the luminosity function that can work as a valuable baseline for interpreting the abundance of fainter lensed sources. In a few years from now a significant leap will be made possible by the availability of deep James Webb Space Telescope imaging. In particular, the capability of reaching 30.5 AB (at $S/N = 5$) in deep NIRCam fields (e.g., Finkelstein et al. 2015) will improve by 1.5 mag the depth of current HFF imaging, reaching absolute magnitudes of $M_{UV} \approx -11$, and will yield five-times larger samples of high-redshift galaxies (Laporte et al. 2015) while significantly improving photometric selections through the availability of rest-frame optical photometry of high-z sources.

We thank the referee for constructive comments that helped to improve the paper.

## References

* Aad, G., Abajyan, T., Abbott, B., et al. 2013, JHEP, 04, 075

* Abazajian, K., Fuller, G. M., & Patel, M. 2001, PhRvD, 64, 023501

* Abazajian, K. N. 2014, PhRvL, 112, 161303

* Abazajian, K. N., Calabrese, E., Cooray, A., et al. 2011, APh, 35, 177

* Ackermann, M., Ajello, M., Albert, A., et al. 2015, PhRvL, 115, 231301

* Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2016, A&A, 594, A13

* Adhikari, R., Agostini, M., Anh Ky, N., et al. 2017, JCAP, 01, 025

* Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2013, PhRvL, 111, 081102

* Adulpravitchai, A., & Schmidt, M. A. 2015, JHEP, 01, 006

* Agarwal, S., Corasaniti, P. S., Das, S., & Rasera, Y. 2015, PhRvD, 92, 063502

* Akerib, D., Araujo, H. M., Bai, X., et al. 2014, PhRvL, 112, 091303

* Angulo, R. E., Hahn, O., & Abel, T. 2013, MNRAS, 434, 3337

* Aprile, E., Alfonsi, M., Arisaka, K., et al. 2012, PhRvL, 109, 181301

* Aprile, E., Arisaka, K., Arneodo, F., et al. 2016, JCAP, 1604, 027

* Arvanitaki, A., Dimopoulos, S., Dubovsky, S., Kaloper, N., & March-Russell, J. 2010, PhRvD, 81, 123530

* Asaka, T., Blanchet, F., & Shaposhnikov, M. 2005, PhLB, 631, 151

* Asaka, T., & Shaposhnikov, M. 2005, PhLB, 620, 17

* Atek, H., Richard, J., Kneib, J.-P., et al. 2015, ApJ, 800, 18

* Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15

* Benson, A., Farahi, A., Cole, S., et al. 2013, MNRAS, 428, 1774

* Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 1107, 034

* Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Natur, 311, 517

* Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2016a, arXiv:1608.00966

* Bouwens, R. J., Oesch, P. A., Illingworth, G. D., Ellis, R. S., & Stefanon, M. 2016b, arXiv:1610.00283

* Boyarsky, A., Ruchayskiy, O., & Iakubovskyi, D. 2009, JCAP, 0903, 005

* Boyarsky, A., Ruchayskiy, O., Iakubovskyi, D., & Franse, J. 2014, PhRvL, 113, 251301

* Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40

* Bozek, B., Marsh, D. J. E., Silk, J., & Wyse, R. F. G. 2015, MNRAS, 450, 209

* Buckley, M. R., Zavala, J., Cyr-Racine, F.-Y., Sigurdson, K., & Vogelsberger, M. 2014, PhRvD, 90, 043524

* Bulbul, E., Markevitch, M., Foster, A., et al. 2014, ApJ, 789, 13

* Calabrese, E., & Spergel, D. N. 2016, MNRAS, 460, 4397

* Canetti, L., Drewes, M., Frossard, T., & Shaposhnikov, M. 2013, PhRvD, 87, 093006

* Castellano, M., Fontana, A., Paris, D., et al. 2010, A&A, 524, A28

* Castellano, M., Amorín, R., Merlin, E., et al. 2016, A&A, 590, A31

* Chavanis, P.-H. 2011, PhRvD, 84, 063518

* Chu, X., & Dasgupta, B. 2014, PhRvL, 113, 161301

* Colombi, S., Dodelson, S., & Widrow, L. M. 1996, ApJ, 458, 1

* Corasaniti, P. S., Agarwal, S., Marsh, D. J. E., & Das, S. 2016, arXiv:1611.05892

* Cyr-Racine, F.-Y., & Sigurdson, K. 2013, PhRvD, 87, 103515

* Dahlen, T., Mobasher, B., Faber, S. M., et al. 2013, ApJ, 775, 93

* de Vega, H. J., Salucci, P., & Sanchez, N. G. 2014, MNRAS, 442, 2717

* de Vega, H. J., & Sanchez, N. G. 2014, arXiv:1401.1214

* Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2014, MNRAS, 441, 2986

* Diemand, J., & Moore, B. 2011, ASL, 4, 297

* Dodelson, S., & Widrow, L. M. 1994, PhRvL, 72, 17

* Du, X., Behrens, C., & Niemeyer, J. C. 2016, arXiv:1608.02575

* Enqvist, K., Kainulainen, K., & Maalampi, J. 1990, PhLB, 249, 531

* Finkelstein, S. L., Dunlop, J., Le Fevre, O., & Wilkins, S. 2015, arXiv:1512.04530

* Garzilli, A., Boyarsky, A., & Ruchayskiy, O. 2015, arXiv:1510.07006

* Ghiglieri, J., & Laine, M. 2015, JHEP, 11, 171

* Gonzáles-Morales, A. X., Marsh, D. J. E., Peñarrubia, J., & Ureña-López, L. 2016, arXiv:1609.05856

* Governato, F., Weisz, D., Pontzen, A., et al. 2015, MNRAS, 448, 792

* Governato, F., Zolotov, A., Pontzen, A., et al. 2012, MNRAS, 422, 1231

* Graham, P. W., Irastorza, I. G., Lamoreaux, S. K., Lindner, A., & van Bibber, K. A. 2015, ARNPS, 65, 485

* Grazian, A., Castellano, M., Fontana, A., et al. 2012, A&A, 547, A51

* Grazian, A., Castellano, M., Koekemoer, A. M., et al. 2011, A&A, 532, A33

* Harko, T., & Lobo, F. S. N. 2015, PhRvD, 92, 043011

* Horiuchi, S., Humphrey, P. J., Onorbe, J., et al. 2014, PhRvD, 89, 025017

* Hu, W., Barkana, R., & Gruzinov, A. 2000, PhRvL, 85, 1158

* Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2015, ApJ, 799, 12

* Jeltema, T. E., & Profumo, S. 2016, MNRAS, 458, 3592

* Kawamata, R., Ishigaki, M., Shimasaku, K., Oguri, M., & Ouchi, M. 2015, ApJ, 804, 103

* Kennedy, R., Frenk, C., Cole, S., & Benson, A. 2014, MNRAS, 442, 2487

* Khlopov, M., Malomed, B. A., & Zeldovich, I. B. 1985, MNRAS, 215, 575

* Kim, J. E. 1987, PhR, 150, 1

* Kishimoto, C. T., & Fuller, G. M. 2008, PhRvD, 78, 023524

* Klasen, M., & Yaguna, C. E. 2013, JCAP, 1311, 039

* Klypin, A., Karachentsev, I., Makarov, D., & Nasonova, O. 2015, MNRAS, 454, 1798

* König, J., Merle, A., & Totzauer, M. 2016, JCAP, 1611, 038

* Kusenko, A. 2006, PhRvL, 97, 241301

* Laine, M., & Shaposhnikov, M. 2008, JCAP, 0806, 031

* Langacker, P. 1989, Pennsylvania Univ. Internal report n. upr-0401T

* Lapi, A., & Danese, L. 2015, JCAP, 1509, 003

* Laporte, N., Infante, L., Troncoso Iribarren, P., et al. 2016, ApJ, 820, 98

* Laporte, N., Streblyanska, K., Kim, S., et al. 2015, A&A, 575, A92

* Lesgourgues, J., & Tram, T. 2011, JCAP, 1109, 032
Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2016, arXiv:1604.06799

Lora, V., Just, A., Sánchez-Salcedo, F. J., & Grebel, E. K. 2012, ApJ, 757, 87

Lotz, J. M., Koekemoer, A., Coe, D., et al. 2016, arXiv:1605.06567

Lovell, M. R., Bose, S., Boyarsky, A., et al. 2016, MNRAS, 461, 60

Lovell, M. R., Eke, V., Frenk, C. S., et al. 2012, MNRAS, 420, 2318

Maccio, A. V., Paduroiu, S., Anderhalden, D., Schneider, A., & Moore, B. 2012, MNRAS, 424, 1105

Madarassy, E. J. M., & Toth, V. T. 2015, PhRvD, 91, 044041

Marsh, D. J. E. 2016, PhR, 643, 1

Marsh, D. J. E., & Pop, A.-R. 2015, MNRAS, 451, 2479

Marsh, D. J. E., & Silk, J. 2014, MNRAS, 437, 2652

Martinez-Medina, L. A., Robles, V. H., & Matos, T. 2015, PhRvD, 91, 023519

Menci, N., Fiore, F., & Lamastra, A. 2012, MNRAS, 421, 2384

Menci, N., Grazian, A., Castellano, M., & Sanchez, N. G. 2016a, ApJL, 825, L1

Menci, N., Sanchez, N. G., Castellano, M., & Grazian, A. 2016b, ApJ, 818, 90

Meneghetti, M., Natarajan, P., Coe, D., et al. 2016, arXiv:1606.04548

Merle, A. 2013, IMPD, 22, 1330020

Merle, A., Niro, V., & Schmidt, D. 2014, JCAP, 1403, 028

Merle, A., & Schneider, A. 2015, PhLB, 749, 283

Merle, A., Schneider, A., & Totzauer, M. 2016, JCAP, 1604, 003

Merle, A., & Totzauer, M. 2015, JCAP, 1506, 011

Mielke, E. W., & Perez, J. A. V. 2009, PhLB, 671, 174

Ng, K. C. Y., Horiuchi, S., Gaskins, J. M., Smith, M., & Preece, R. 2015, PhRvD, 92, 043503

Pacucci, F., Mesinger, A., & Haiman, Z. 2013, MNRAS, 435, L53

Papastergis, E., Giovanelli, R., Haynes, M. P., & Shankar, F. 2015, APh, 574, A113

Peebles, P. J. E. 1982, ApJL, 263, L1

Petraki, K., & Kusenko, A. 2008, PhRvD, 77, 065014

Polisensky, E., & Ricotti, M. 2011, PhRvD, 83, 043506

Riemer-Sorensen, S. 2016, A&A, 590, A71

Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013, MNRAS, 430, 81

Schive, H.-Y., Chiueh, T., & Broadhurst, T. 2014, NatPh, 10, 496

Schive, H.-Y., Chiueh, T., Broadhurst, T., & Huang, K.-W. 2016, ApJ, 818 89

Schneider, A. 2015, MNRAS, 451, 3117

Schneider, A. 2016, JCAP, 1604, 059

Schneider, A., Smith, R. E., Maccio, A. V., & Moore, B. 2012, MNRAS, 424, 684

Schneider, A., Smith, R. E., & Reed, D. 2013, MNRAS, 433, 1573

Schultz, C., Onorbe, J., Abazajian, K. N., & Bullock, J. S. 2014, MNRAS, 442, 1597

Sekiya, N., Yamasaki, N. Y., & Mitsuda, K. 2016, PASJ, 68, 31

Seljak, U., Makarov, A., McDonald, P., & Trac, H. 2006, PhRvL, 97, 191303

Shakya, B. 2016, MPLA, 31, 1630005

Shaposhnikov, M. 2008, JHEP, 08, 008

Shaposhnikov, M., & Tkachev, I. 2006, PhLB, 639, 414

Shi, X.-D., & Fuller, G. M. 1999, PhRvL, 82, 2832

Steffen, F. D. 2006, JCAP, 9, 001

Tamura, T., Iizuka, R., Maeda, Y., Mitsuda, K., & Yamasaki, N. Y. 2015, PASJ, 67, 23

Vanzella, E., Balestra, I., Gronke, M., et al. 2016, arXiv:1607.03112

Venumadhav, T., Cyr-Racine, F.-Y., Abazajian, K. N., & Hirata, C. M. 2016, PhRvD, 94, 043515

Viel, M., Becker, G. D., Bolton, J. S., & Haehnelt, M. G. 2013, PhRvD, 88, 043502

Viel, M., Lesgourgues, J., Haehnelt, M. G., Matarrese, S., & Riotto, A. 2006, PhRvL, 97, 071301

Vogelsberger, M., Zavala, J., Simpson, C., & Jenkins, A. 2014, MNRAS, 444, 3684

Wang, M.-Y., Peter, A. H. G., Strigari, L. E., et al. 2014, MNRAS, 445, 614

Watson, C. R., Li, Z.-Y., & Polley, N. K. 2012, JCAP, 1203, 018

Woo, T.-P., & Chiueh, T. 2009, ApJ, 697, 850