Small scale clustering of late forming dark matter
Shankar Agarwal, Pier Stefano Corasaniti, Subinoy Das, Yann Rasera
SSmall scale clustering of late forming dark matter
S. Agarwal , P.-S. Corasaniti , S. Das , Y. Rasera Laboratoire Univers et Th´eories (LUTh), UMR 8102,CNRS, Observatoire de Paris, Universit´e Paris Diderot,5 Place Jules Janssen, 92190 Meudon, France and Indian Institute of Astrophysics, 560034 Bangalore, India (Dated: November 10, 2018)We perform a study of the nonlinear clustering of matter in the late-forming dark matter (LFDM)scenario in which dark matter results from the transition of a nonminimally coupled scalar field fromradiation to collisionless matter. A distinct feature of this model is the presence of a damped oscil-latory cutoff in the linear matter power spectrum at small scales. We use a suite of high-resolutionN-body simulations to study the imprints of LFDM on the nonlinear matter power spectrum, thehalo mass and velocity functions and the halo density profiles. The model largely satisfies high-redshift matter power spectrum constraints from Lyman- α forest measurements, while it predictssuppressed abundance of low-mass halos ( ∼ − h − M (cid:12) ) at all redshifts compared to a vanillaΛCDM model. The analysis of the LFDM halo velocity function shows a better agreement thanthe ΛCDM prediction with the observed abundance of low-velocity galaxies in the local volume.Halos with mass M (cid:38) h − M (cid:12) show minor departures of the density profiles from ΛCDMexpectations, while smaller-mass halos are less dense, consistent with the fact that they form laterthan their ΛCDM counterparts. I. INTRODUCTION
The cold dark matter (CDM) scenario has beentremendously successful in reproducing observations ofthe distribution of matter on cosmic scales [1–5]. In spiteof this remarkable success the origin of this invisible com-ponent is still not known. This is because cosmologicalobservations shed no light on the particle physics natureof DM and only suggest that DM consists of an approxi-mately pressureless component that clusters gravitation-ally since a few e-foldings before matter-radiation equal-ity.The leading high-energy physics candidates to DM areweakly interacting massive particles (WIMPs). In prin-ciple these can be detected in underground laboratoryexperiments with a signal characterized by an annualmodulation pattern. This has been thought to be a dis-tinctive signature of DM direct detection against pos-sible contamination from radioactive background noise.In fact the annual modulation is commonly considered aconsequence of the orientation of the Earth’s orbit alongor against the flow of DM particles in the Galactic haloduring the annual motion of the Earth around the Sun(for a review see Ref. [6]). A number of underground ex-periments such as DAMA/LIBRA [7], CoGeNT [8] andCRESSST-II [9] have claimed the detection of DM par-ticles with an annual modulation pattern. In contrast,the CDMS-II [10], XENON10 [11] and more recentlyLUX [12] collaborations have reported none. Searches forelectroweak WIMPs at the Large Hadron Collider havealso given negative results so far. More puzzling is thefact that the detections from the underground experi-ments constrain different regions of the WIMPs’ param-eter space [13]. Thus, it is possible that such contrastingoutcomes may be hiding a much richer physics in theweak scale dark sector than previously thought or indi- cate the need of a completely new dark matter paradigmbeyond weak scale WIMPs.The small-scale clustering of matter in the Universe hasemerged as the new arena to test the nature of DM. Thediscovery of a number of anomalies at small scales hascast doubts on the validity of the CDM hypothesis. Thecore-vs-cusp problem [14, 15] and the missing satelliteproblem (see e.g. [16, 17]) have motivated the study ofalternative scenarios beyond the CDM paradigm. More-over, recent studies of the dynamical properties of themost luminous Milky Way dwarf spheroidal satellitegalaxies have pointed to a new anomaly, the so-calledtoo-big-to-fail problem [18, 19]. Baryonic processes ingalaxy formation have been invoked as the natural solu-tion to these discrepancies (see e.g. Refs. [20–30]). Asan example, baryon feedback and observational incom-pleteness can account for the missing satellite problem(see e.g. Refs. [31, 32]). Similarly, statistical variationsin the predictions of subhalo abundances from N-bodysimulations combined with the uncertain value of theMilky Way virial mass may considerably alleviate thetoo-big-to-fail problem [33, 34]. However, it is still un-clear whether baryonic feedback models can provide aunique self-consistent explanation to the entirety of DManomalies (see e.g. Refs. [35–37]). Furthermore, analy-ses of dwarf galaxies in the local field have shown thattheir abundance and properties differ significantly fromthe ΛCDM predictions [38–43]. This is more difficult toreconcile with the CDM paradigm in terms of baryonfeedback models given the fact that such isolated sys-tems, contrary to satellites, undergo less complex pro-cesses. Therefore, it cannot be a priori excluded thatsuch anomalies in the small-scale clustering of matter are(also) related to the unknown nature of DM (for a reviewsee Ref. [44]) or a hint of broken scale invariance of theinflationary power spectrum at small scales [45–47]. a r X i v : . [ a s t r o - ph . C O ] A ug Among the scenarios alternative to CDM, the existenceof a warm dark matter (WDM) component has been mo-tivated by particle physics models of sterile neutrinos (seee.g. Refs. [48–52]). Differently from the CDM, WDMparticles with mass of a few keV free-stream on a scale (cid:46)
100 kpc due to velocity dispersion [47, 53]. This causesa characteristic suppression of the linear matter densitypower spectrum at small scales and alters the proper-ties of DM halos. Comparison with observations thusprovides bounds on the DM particle mass. As an exam-ple, it was pointed out in Ref. [54] that the solution tothe core-vs-cusp problem requires WDM particles withmass m WDM ∼ . . (cid:46) m WDM [keV] (cid:46) α forest measurements [58–60]. Forinstance, the recent analysis by Viel et al. [61] indicatesa lower limit m WDM (cid:38) . σ confidence level.Even assuming a particle mass within such a bound (e.g. m WDM = 4 keV) the authors of Ref. [62] have shown thatthe dynamical properties of galactic WDM subhalos areindistinguishable from those of CDM, thus leaving thetoo-big-to-fail problem unsolved. Indeed, if the small-scale anomalies are entirely due to the particle physicsnature of dark matter, then requiring their simultaneoussolution poses strong constraints on DM models (see e.g.Ref. [63]).Another class of models alternative to the standardCDM hypothesis has developed around the possibilitythat DM particles have significant self-interactions [64].This scenario, also known as self-interacting dark matter(SIDM), arises in the context of particle physics mod-els of the DM sector (see e.g. Refs. [65–67]). A keyprediction of these models is the fact that DM par-ticles have velocity-dependent self-scattering cross sec-tions. These self-interactions, which can be assimilatedto gravitational scalar forces, alter the nonlinear gravi-tational collapse, leaving characteristic signatures in thesmall-scale clustering of matter. N-body simulations ofSIDM models have been presented in numerous studies.Recent analysis (see e.g. [68–70]) has shown that in anarrow range of the SIDM model parameter space thedynamics of subhalos in parent halos with mass aboutthat of the Milky Way reproduces that of the dwarf-spheroidal galaxies in the Milky Way. At the same time,the subhalos possess core profiles O (1 kpc), thus simul-taneously solving the core-cusp problem and the too-big-to-fail problem. On the other hand, in SIDM models thesubhalo mass function remains identical to that of CDM,thus leaving the missing satellite problem unsolved. Ithas been shown in Refs. [71–73] that such a class of mod-els can provide a successful solution if an additional inter-action of the DM component with radiation is considered.A different family of models consists of scenarios in which DM is the result of a decay or a phase transitionprocess. For instance, the authors of Ref. [74] have ex-plored the possibility that neutral dark matter particlesare the result of the decay of charged particles coupledto the photon-baryon plasma before recombination, whileRef. [75] has proposed a scenario in which early dark mat-ter particles decay into a less massive species and a mass-less noninteracting one. In the latter case the resultingspecies eventually disrupt the formation of low-mass ha-los as investigated in Refs. [76, 77], while still satisfyingbounds from Lyman- α observations [78].Here, we consider a scenario in which the DM is conse-quence of a phase transition in the dynamics of a scalarfield. In particular, we focus on late-forming dark mat-ter (LFDM) [79], inspired by particle physics models ofneutrino dark energy [80] which aim to provide a unifieddescription of dark matter and dark energy. In LFDMmodels, DM particles are the result of a phase transitionin the equation of state of a scalar field from radiation( w ∼ /
3) to matter ( w ∼
0) (see Ref. [81] for a differentrealization of this scenario). A distinctive characteristicof LFDM models is the presence of a damped oscilla-tory tail at the small scales of the linear matter powerspectrum. However, contrary to other nonstandard DMmodels, the range of scales where such a distinctive fea-ture occurs is set by the epoch of the phase transition(rather than the value of the DM particle mass or the am-plitude of the self-interaction cross section). The earlierthe transition, the smaller the scale where the suppres-sion of power occurs. In this regard, models of ultralightaxion (ULA) dark matter can also be seen as a form oflate-forming dark matter where the axion field transitionsfrom a vacuum state ( w ∼ −
1) to matter ( w ∼
0) whenthe field mass becomes comparable to the Hubble scale(see e.g. Ref. [82]). ULA dark matter alters the linearcosmic structure formation; thus it is constrained by ob-servations of the large-scale structures as recently shownin [83]. This occurs in LFDM only if the transition toDM state occurs after recombination, otherwise LFDMexclusively affects the linear matter power spectrum atsmall scales, thus leaving potentially observable featuresin the nonlinear structure formation. Given the similar-ities between the ULA and LFDM models, it is reason-able to expect that in some viable region of the ULAparameter space, the two scenarios may share the samephenomenology of the matter clustering at small scales,which is yet to be studied. In this paper we present afirst study of the nonlinear structure formation of LFDMand assess the viability of this scenario using a series ofhigh-resolution N-body simulations.The paper is organized as follows: In Sec. II we recallthe main features of the LFDM model. In Sec. III wedescribe the numerical methods, the generation of ini-tial conditions and N-body simulation characteristics. InSec. IV we analyze the nonlinear matter power spectrum.In Sec. V we describe the imprints on the halo mass func-tion. In Sec. VI we present the results from the halo den-sity profiles. In Sec. VII we probe the circular velocitydistribution, and we conclude in Sec. VIII.
II. COSMOLOGY OF LATE-FORMING DARKMATTER
In the following, we will briefly review the main fea-tures of the LFDM scenario that are relevant to the studyof cosmic structure formation. We refer the reader to theoriginal paper by Das and Weiner [79] for a detailed de-scription of the model.In this scenario, DM is the result of the dynamical evo-lution of a scalar field coupled to a thermal bath of rela-tivistic particles (e.g. eV sterile neutrinos). At high tem-perature (early times) the field is trapped in a metastablestate and behaves as a dark radiation fluid with an equa-tion of state w ∼ /
3. As temperature drops due to thecosmic expansion, a new lower energy minimum appearsand the field rolls into the lower energy state while oscil-lating around the true minimum of the scalar potential,thus behaving as a collisionless DM component ( w ∼ z t .The later the phase transition, the larger the scales wherethe power spectrum exhibits a damped oscillatory tail.The absence of these features in the large-scale powerspectrum of galaxies imposes constraints on the redshiftof the transition. In particular, a recent analysis indi-cates that z t > [87]. Since after the phase transitionLFDM behaves as a collisionless component, the linearcosmic structure formation is indistinguishable from thatof the standard ΛCDM model. However, in order to pro-duce sufficient dark matter abundance, an excess of scalardark radiation (in addition to that of neutrinos and pho-tons) is needed at early times. As shown in Ref. [87], asmall excess contributing to ∆ N eff ∼ . − .
1, thus wellwithin the CMB constraints [88], is sufficient to give theright amount of dark matter particles.As an example, in the left panel of Fig. 1 we plotthe cosmic microwave background (CMB) temperatureanisotropy power spectrum for a flat Λ-LFDM modelwith cosmological constant against a ΛCDM model bestfit to Planck data [88, 89]. In this specific LFDM modelthe phase transition occurs at z t = 1 . × , and the cos-mological parameters are set to the values of the standardΛCDM model. We have computed the CMB and matterpower spectrum of the LFDM model using a modified version of the camb code [90] to solve the dynamics ofthe coupled linear perturbation equations associated withthe LFDM component. We may notice that even with-out a model parameter optimization requiring a MonteCarlo likelihood analysis of available data, the choice ofthe Λ-LFDM model parameters provides a good fit tothe Planck measurements (the small differences arise es-sentially from the small excess of dark radiation densitybefore matter-radiation equality). In the right panel ofFig. 1, we plot the corresponding linear matter powerspectra at redshift z = 0 against the luminous red galaxy(LRG) power spectrum from the Sloan Digital Sky Sur-vey seventh data release (SDSS DR7) [91]. As we cansee the spectra are identical on the large scales, whiledifferences are only present at k (cid:38)
10 h Mpc − due toLFDM damped oscillations on nonlinear scales. The ear-lier the phase transition (i.e. higher z t ), the smaller thepower spectrum cutoff scale, and the damped part of thespectrum shifts towards larger wave numbers.It is worth mentioning that the linear LFDM powerspectrum considered here resembles that of the modelsstudied in Refs. [72, 73]. Nonetheless, there are a few no-ticeable differences. In Ref. [72] the authors have consid-ered nonstandard models with DM-photon interactionscharacterized by a linear power spectrum that in therange 10 < k [ h Mpc − ] <
100 has a damped oscillatorypattern with a slope different from that of the LFDM case(see Fig. 1 in Ref. [72]), while in Ref. [73] the amplitude ofthe damped oscillations is significantly more suppressedthan in our model realization (see Fig. 1 in Ref. [73]).Certain realizations of mixed cold and warm dark mat-ter models such as those considered in Refs. [52, 92] alsohave spectra that approximate the envelope of the LFDMspectrum; however, they have different power distribu-tion in the range 10 < k [ h Mpc − ] < m WDM = 1 .
465 keV characterized by a cutoff scale sim-ilar to that of LFDM, and of a broken-scale-invariance(BSI) inflationary model [93, 94] as the one considered inRef. [45] with highly tuned parameters such as to havethe same cutoff scale and a slope of power suppressionsimilar to that of the LFDM model. We can see thatthe spectra still carry differences that are likely to alterthe nonlinear structure formation. In particular, in thecase of the BSI spectrum, the model predicts an excessof power compared to the ΛCDM case before the cutoff; The primordial spectrum for BSI models [93, 94] has a universalform and depends only on two parameters, the cutoff scale andthe amplitude of the power suppression p . In Refs. [46, 47] it wasassumed p = 4, while to closely reproduce the LFDM spectrumwe set p = 100. FIG. 1: Left panel: CMB temperature anisotropy power spectrum for a flat ΛCDM model (black solid line) best fit to Planckdata [88, 89] (filled squares) and a flat Λ-LFDM model (red dot-dashed line) with identical cosmological parameters. Rightpanel: Linear matter density power spectra for the ΛCDM and Λ-LFDM models shown in the left panel. Data points correspondto LRG power spectrum from SDSS-DR7 [91]. For illustrative purposes only we also show the linear power spectrum for aWDM model with thermal relic particle mass m WDM = 1 .
465 keV (blue dotted line) and a broken-scale-invariant inflationarymodel (cyan dotted line). thus it is likely to be highly constrained by high-redshiftLyman- α measurements. This is not the case for theWDM model, which differs from the LFDM beyond thecutoff scale, where the power exponentially drops sev-eral orders of magnitude below the envelope of LFDMdamped oscillations. Even considering a slightly largerthermal mass particle, this will only shift the exponen-tial tail of the WDM spectrum to higher wave numbers,while still remaining more suppressed than the LFDMprediction due to the shallower slope of the LFDM en-velope. Hence, it is reasonable to expect that such dif-ferences may manifest with different predictions of theabundance and properties of low-mass halos. It is forthese very reasons that, despite the numerous numeri-cal studies of the nonlinear clustering of WDM models,it is opportune to investigate the phenomenology of theLFDM scenario and the like. III. N-BODY SIMULATIONS
We run a series of high-resolution N-body simulationsusing ramses [95], an adaptive mesh refinement (AMR)code with a tree-based data structure in which particlesare evolved using a particle-mesh solver and the Poissonequation is solved using a multigrid method [96].We generate the initial conditions with the code mp-grafic [97], which uses the Zel’dovich approximation.We set the initial redshift such that the standard devia- tion of the initial density field smoothed on the scale ofthe coarse grid ∆ coarse x is σ (∆ coarse x ) = 0 .
02. With thischoice, the initial redshift of the simulations is sufficientlylarge to suppress spurious effects due to transients [98].We assume a flat Λ-LFDM model with phase transi-tion redshift z t = 1 . × and characterized by thelinear power spectrum shown in Fig. 1. The cosmologi-cal parameters are set to the following values: Ω m = 0 . . σ = 0 . n s = 0 .
96 and Ω b = 0 . . − Mpc) volume with 512 particles, a sim-ilar volume with 1024 particles, an intermediate boxof (64 h − Mpc) with 512 particles and a large box of(110 h − Mpc) volume with 2048 particles. This sim-ulation suite enables us to control numerical systematicerrors due to mass resolution and volume effects. Thecharacteristics of the simulations are summarized in Ta-ble I. We can see that for the highest mass resolution runsthe spatial resolution of the simulations at the level of thecoarse grid varies from 27 h − kpc to 54 h − kpc. However,the actual resolution of the simulations is much higherdue to the AMR scheme used in the ramses code. Forinstance, in the case of the (27 . − Mpc) volume simu-lation with 1024 particles we have up to six refinementlevels of the coarse grid which are triggered during therun; thus the densest objects are resolved with a spatialresolution of ≈ . − kpc. L (h − Mpc) N p z ini m p (h − M (cid:12) ) ∆ coarse x (h − kpc)27 .
320 1 . × .
324 1 . ×
303 1 . ×
320 1 . × N p is the number of N-body parti-cles, z ini is the initial redshift of the simulation, m p the massresolution, and ∆ coarse x is the spatial resolution of the coarsegrid. In addition to the LFDM model simulations, we haverun a (27 . − Mpc) volume with 512 particles and onewith 1024 particles for a flat ΛCDM model with identicalcosmological parameters (with z ini = 424 and z ini = 502,respectively) and the same phase of the initial conditions,which we use to evaluate differences with respect to thesmall-scale clustering of Λ-LFDM.The simulations were run on the ADA supercomputerof the Institute for Development and Resources in In-tensive Scientific Computing (IDRIS). In particular, theLFDM simulation of (110 h − Mpc) volume with 2048 particles was run on 3000 Intel Sandy Bridge E5-4650processors for a total running time of 3 × hours. IV. NONLINEAR MATTER POWERSPECTRUM
We compute the matter power spectrum with the code powergrid [97]. This computes the power spectrumby performing a Fourier transform of the density fieldin band powers ∆ k = 2 π/L , where L is the simulationbox length. We correct for the smoothing effect due tothe cloud-in-cell (CIC) algorithm used to estimate thedensity field from the particle distribution. To be con-servative, we restrict the scope to wave numbers belowthe Nyquist frequency of the coarse grid.In the left panel of Fig. 2, we plot the nonlinear powerspectra at z = 0 , z = 0 case,with similar trend at higher redshifts. First, we maynotice that the spectra of the large (solid line) and in-termediate (dotted line) volumes are fairly in agreementat k < − , while in the same range the am-plitude of the spectrum of the smaller simulation box(dashed line) is slightly lower. This is due to finite vol-ume effects as well as the choice of the initial phase (seee.g. Refs. [99–101]). Volume effects are negligible atlarge k ; in particular for k >
13 h Mpc − , the spec-trum of the (27 . − Mpc) volume simulation with 512 particles and that of the (110 h − Mpc) box differ byless than 2%. On the other hand, at larger wave num-bers (near the Nyquist frequency of the simulations), wecan see the systematic suppression of power due to massresolution errors for the lower resolution runs compared to the higher one. At small scales, this is the domi-nant source of numerical uncertainty; in particular, inthe range 10 < k [h Mpc − ] <
50, the spectra of the(27 . − Mpc) volume simulation with 512 and thatwith 1024 particles differ by a few percent in the low-end limit and up to 50% in the high end. However, sincewe are interested in the relative difference between theLFDM and ΛCDM spectra, we expect (and show in theright panel of Fig. 2) the ratio to be less affected by massresolution errors.We plot in the right panel of Fig. 2 the relative differ-ence of the nonlinear matter power spectrum of the Λ-LFDM model with respect to the ΛCDM case for k (cid:38)
10h Mpc − at z = 0 , , , . . − Mpc) volume simulations. At z = 5 .
5, the runs with 512 and 1024 particles differ from less than a percent at k = 10 h Mpc − and increase up to 4% at the high endof the interval. These differences due to mass resolutionare clearly subdominant compared to the relative differ-ences between the models. Volume effects in this rangeof wave numbers are expected to be negligible. Hence,we are confident the differences seen between Λ-LFDMand ΛCDM are not due to numerical artifacts. Noticethat any signature of the LFDM damped oscillations inthe initial linear power spectrum has been completelyerased.Constraints on the high-redshift matter power spec-trum (3 < z < .
4) inferred from the recent Lyman- α measurements [61] indicate that the nonlinear powerspectrum for WDM thermal relic particles of mass m wdm > . σ confidence level) may deviatefrom ΛCDM by no more than ∼
5% at k = 10 h Mpc − (as can be inferred from Fig. 1 in Ref. [61]). As shownin the right panel of Fig. 2, the Λ-LFDM model consid-ered here is consistent with these bounds for z (cid:46)
4, whilemodels with slightly higher z t will largely evade theseconstraints in the entire high-redshift range. V. HALO MASS FUNCTIONA. Numerical convergence analysis
We detect DM halos using the code pFoF [102] basedon the friend-of-friend algorithm [103], which identifieshalos as group of particles characterized by an intraparti-cle distance smaller than a given linking length parame-ter b . In our analysis we set this to the standard value b = 0 . FIG. 2: Left panel: Nonlinear matter power spectrum of the Λ-LFDM model from the simulations of box length L = 110 h − Mpcwith 2048 particles (solid line), L = 64 h − Mpc with 512 particles (dotted line), and L = 27 . − Mpc with 512 (short dashedline) and 1024 (long dashed line) particles at z = 0 , z = 0 , , , . . − Mpc) volume with 512 (black solid line) and 1024 (red dotted line) particles. ysis of the Λ-LFDM model presented here and refer thereader to Ref. [107] for a dedicated study on the re-moval of spurious halos. Our spurious halo identificationmethod relies on the fact that since spurious halos orig-inate from the fragmentation of the density field, theseare highly nonspherical groups of particles characterizedby extreme values of the spin and shape parameters whileexhibiting large deviation from the virial theorem. Thus,in the range of masses where mass resolution effects aresubdominant, spurious halos can be removed by simplyretaining those halos that approximately satisfy the virialcondition and checking that the distribution of the spinand shape parameters of the remaining halos have well-behaved tails.We compute the mass function as dnd ln M = 1 L N ∆ ln M , (1)where L is the simulation box length and N is the numberof halos in a mass bin of size ∆ ln M . Throughout thispaper, log and ln denote base-10 and base- e logarithms,respectively.In Fig. 3 we plot the mass function of the Λ-LFDMmodel at z = 0 from simulations with volumes of(110 h − Mpc) with 2048 particles (blue filled squares),(64 h − Mpc) with 512 particles (brown open circles),(27 . − Mpc) with 512 (green open squares) and1024 (red filled circles) particles in bins of size ∆ M/M =0 .
1. In all four cases we can see that the high-mass endof the mass function is characterized by a large level of scatter due to finite volume effects. Mass resolution ef-fects are relevant in the low-mass end; this can be seenin the case of the intermediate simulation box. for whichthe mass function deviates at more than a 5% level for
M < × h − M (cid:12) from the higher-resolution runs.Likewise, the comparison between the mass functions ofthe (27 . − Mpc) volume simulation with 512 parti-cles and that with 1024 shows deviation at more than a5% level for M < × h − M (cid:12) . It is worth noticingthe drop of the mass function from the (27 . − Mpc) volume run with 1024 at low masses; this is consistentwith the small-scale cutoff in the linear matter powerspectrum. Overall, we find convergence at the 5% level inthe mass range 5 × (cid:46) M [h − M (cid:12) ] (cid:46) for the massfunctions from the (27 . − Mpc) and (110 h − Mpc) volume simulations. B. Redshift evolution and cosmological imprints
In Fig. 4 we plot the redshift evolution of the Λ-LFDMmass function from the (110 h − Mpc) volume simulationat z = 0 , , , . z = 0 and extrapolated tohigher redshifts (see the Appendix).We can see that despite the suppressed amplitude ofthe linear matter power spectrum at small scales (seeright panel in Fig. 1), the halo formation proceeds as in FIG. 3: Halo mass function of the Λ-LFDM model at z = 0in bins of size ∆ M/M = 0 . − Mpc) with 2048 particles (blue filled squares),(64 h − Mpc) with 512 particles (brown open circles), and(27 . − Mpc) with 512 (green open squares) and 1024 (red filled circles) particles for halos with at least 100 particles.Error bars are given by Poisson errors.FIG. 4: Redshift evolution of the Λ-LFDM halo mass func-tion at z = 0 (blue open squares), 1 (red solid circles), 3 (ma-genta open triangles), 5 . − Mpc) volume simulation. Error bars aregiven by Poisson errors. The solid lines show the mass func-tion assuming the Sheth-Tormen multiplicity function withparameters calibrated at z = 0. a hierarchical bottom-up scenario, and the abundance ofhalos increases from high to low redshifts as a functionof halo mass. At z = 0, the ST mass function provides agood fit to the N-body measurements to better than 10%in the mass range 10 < M [h − M (cid:12) ] < × . Forlarger masses, deviations from the ST fit are dominatedby finite volume errors. In contrast, at lower masses, theST overestimates the abundance of low-mass halos withdeviations up to 30%. At redshifts z = 1 , z = 9, the low-mass endis still overestimated, while the abundance of halos withmass M (cid:38) h − M (cid:12) is underestimated. This trend isa manifestation of departure from the universality of thehalo mass function already found in standard CDM cos-mologies (see e.g. Ref. [109–111]). Nevertheless, the factthat ST systematically overestimates the abundance oflow-mass halos ( M (cid:46) h − M (cid:12) ) suggests a substantialdeparture of the nonlinear gravitational dynamics fromthat of the standard CDM ellipsoidal collapse model en-coded in the ST multiplicity function (see the Appendix).To see this more clearly we compare the multiplicity func-tion of the Λ-LFDM model to that of the ΛCDM. For thiswe use data from the (27 . − Mpc) volume simulationwith 512 particles. As with the case of the relative differ-ence between the matter power spectra presented earlier,the ratio of the multiplicity functions is less affected bymass resolution errors.In the left panel of Fig. 5, we plot the ratio of the mul-tiplicity functions at z = 0 , , , . M/M = 0 . z = 0 to z = 9, as volume effects are largest at thehigh-mass end. We can see that independently of theredshift, for M (cid:38) h − M (cid:12) the multiplicity functionsof the two simulated cosmologies are identical to numeri-cal precision. We may also notice that in this mass rangethe ratio of the ST multiplicity functions calibrated tothe data at z = 0 for the two cosmological models re-produces fairly well the numerical results. This is notthe case at lower masses, where ST predictions largelydeviate from the N-body results. This suggests that theuse of standard mass function formulas developed in theframework of the CDM paradigm to infer predictions fornonstandard CDM models characterized by a cutoff inthe linear matter power spectrum at small scales maylead to significant misestimates of the halo abundancesat low mass.Let us now turn to the comparison of the halo abun-dance between ΛCDM and Λ-LFDM. In the right panelof Fig. 5, we plot the mass function of ΛCDM halos rel-ative to that of the Λ-LFDM at z = 0 , , , . M/M = 0 . M (cid:38) h − M (cid:12) which is consistentwith the fact that in this range the multiplicity functions FIG. 5: Left panel: Ratio of the multiplicity function of the ΛCDM model to that of the Λ-LFDM from the (27 . − Mpc) volume simulations with 512 particles in mass bins of size ∆ M/M = 0 . z = 0 to 9 as in Fig. 4. Error bars are given by the propagation of Poisson errors. The different lines (solid, short dashed,dotted, long-dashed, dot-dashed) correspond to the ratio at different redshifts (0 to 9, respectively) obtained assuming theSheth-Tormen multiplicity function with parameters calibrated at z = 0 for the ΛCDM and Λ-LFDM models. Right panel:Ratio of the ΛCDM mass function to that of the Λ-LFDM model in mass bins of size ∆ M/M = 0 . . − Mpc) box simulations, as in the left panel. of the two models are identical. Furthermore, in thisrange also the variance of the smoothed linear densityfield predicted by the two cosmologies is nearly identical,thus leading to a similar halo abundance well reproducedby the ST mass function. In contrast, we can see thatfor M (cid:46) h − M (cid:12) , halos in ΛCDM are up to a fac-tor of 3 more abundant than in Λ-LFDM. This is dueto the differences in the multiplicity function shown inthe left panel of Fig. 5 and amplified by the fact that onthese scales the smoothed linear density field of the Λ-LFDM model is suppressed compared to the ΛCDM case.Also notice that there is no redshift evolution of the rel-ative abundance between z = 3 and 0, while a modestevolution occurs at higher redshifts for halos with mass M ∼ h − M (cid:12) . This is consistent with the fact thatin Λ-LFDM, small-mass halos form later than in ΛCDM,thus leading to differences that can be tested throughhigh-redshift universe observations as already exploredin the case of WDM models in Ref. [112]. The lowerabundance of low-mass halos in LFDM models comparedto ΛCDM is a direct consequence of the suppression ofpower at small scale in the linear matter power spectrumbeyond the cutoff scale. Since in the LFDM scenario thelocation of the cutoff depends on the epoch of dark mat-ter formation, we can expect that DM forming earlierthan that in the model considered here will decrease thediscrepancy with respect to the ΛCDM model, while a later formation will increase it. Thus, measurements ofthe abundance of dwarf galaxies may provide direct con-straints on the epoch of dark matter formation in LFDMmodels. VI. HALO DENSITY PROFILE
We now focus on the density profile of DM halos. InFig. 6 we plot the relative difference of the ΛCDM aver-aged halo density profile at z = 0 with respect to thatof the Λ-LFDM case for halos in three different massbins. Within each mass bin, individual halo profiles arestacked to obtain the averaged profile. To facilitate stack-ing, radius is normalized by r corresponding to an en-closed overdensity ∆ = 178 relative to the cosmic meanmatter density. To control numerical resolution effects,we plot results from the (27 . − Mpc) box simulationswith 512 (blue lines) and 1024 (red lines) particles.Moreover, we limit the curves up to the physical scaleassociated with the size of the most refined cell of thesimulations. We can see that the lower-resolution simula-tions systematically underestimate the profile differenceswith respect to the higher-resolution runs with deviationsup to ∼
3% level.We find that density profiles for halos with mass
M > h − M (cid:12) are nearly identical in the two cosmologies FIG. 6: Relative difference of the averaged halo density profileat z = 0 in ΛCDM with respect to the Λ-LFDM model, as afunction of radius normalized by r . The averaged profilesare obtained by stacking the halo profiles within mass bins M = 10 − h − M (cid:12) (dashed line), 10 − h − M (cid:12) (dotted line) and 10 − h − M (cid:12) (solid line) for the(27 . − Mpc) box simulations with 512 (blue lines) and1024 (red lines) particles. for r/r (cid:38) .
1. At smaller radii, deviations do notexceed the 10%, level with ΛCDM halos being on averagedenser than in Λ-LFDM. In the case of halos with massin the range M = 10 − h − M (cid:12) , we find largerdeviations. In particular, in the CDM case the averagedprofile at radii r/r (cid:46) . ∼
20% denser than theLFDM counterpart. As already noticed in Ref. [73], thisis consistent with the fact that the two models only differfor the amplitude of the initial power spectrum at smallscales. In fact, as shown by the analysis of the massfunction, the suppression of power at small scales in theΛ-LFDM model causes small-mass halos to form laterthan in ΛCDM, thus resulting in lower density profiles.Halos of higher mass, on the other hand, assemble nearlyat the same time in the two cosmologies, thus leading tonearly identical profiles. If the trend shown in Fig. 6 atsmall radii is extrapolated to lower-mass halos, then theseresults do not exclude the possibility that halos with mass
M < h − M (cid:12) may have cored profiles; however, weare not able to address this point at the moment, sinceour simulations do not possess the required resolution. VII. CIRCULAR VELOCITY DISTRIBUTION
The abundance of dark matter halos can be probedthrough measurements of the galaxy luminosity functionprovided prior knowledge of the relation between galaxyluminosity and halo mass exists. However, this is hard to predict, since it depends on the galaxy formation processitself. In contrast, the velocity function, i.e. the numberdensity of galaxies with a given circular velocity, can bepredicted more reliably and tested against observations.It is indeed the comparison of the circular velocity ofthe satellites of the Local Group against N-body simula-tion results that originally pointed to the missing satel-lite problem [16, 17]. Another advantage of the velocityfunction is that it probes the abundance of all virializedstructures. In the past few years, analyses of the ve-locity function inferred from a number of surveys of thelocal volume have pointed to suppressed abundances oflow-velocity field galaxies compared to expectations fromΛCDM simulations [38, 39, 43]. These findings have beenrecently confirmed by the measurements of the circularvelocity distribution in galaxies within 10 Mpc from theLocal Volume catalog [41]. The circular velocity distribu-tion is indeed a promising tool to constrain DM proper-ties; however, it is worth mentioning that any analysis us-ing DM-only simulations can introduce a mass-dependentbias, since the measured circular velocity is that of thebaryon component inside dark matter halos.Here, we compare the measurements of the circu-lar velocity distribution from Ref. [41] against esti-mates from our ΛCDM and LFDM halo catalogs of the(27 . − Mpc) volume simulations with 1024 particles.To derive predictions from DM-only simulations that canbe compared to observations, we follow the procedure de-scribed in Ref. [41]. More specifically, we correct the esti-mated halo circular velocities to account for the effect ofbaryons using the prescription described in Ref. [41] andbased on the analysis by Ref. [113]. Then, we multiplythe estimated velocity functions by a factor 1 .
25 to ac-count for the contribution of subhalos not included in ournumerical catalogs. As pointed out in Ref. [41], this fac-tor corresponds to the fraction of subhalos in ΛCDM forvelocities V (cid:46)
200 km s − . In the LFDM scenario, thefraction of subhalos may be even smaller due to the powersuppression at small scales; thus the choice adopted hereis a rather conservative one.In Fig. 7, we plot the differential circular velocity func-tion dN/d log V for the ΛCDM (red triangles) and Λ-LFDM (blue squares) models. The dotted straight lineis the fitting function to the ΛCDM simulation given byEq. (5) in Ref. [41], while the dashed line includes thebaryonic corrections. Since the cosmological parametersconsidered here are marginally different from Ref. [41],we normalize these fitting functions to the circular veloc-ity function obtained from our ΛCDM simulation. Thecyan solid lines enclose the region of the observed veloc-ity function (with a fitting function given by Eq. (12) inRef. [41]) to within 15% error. As we can see, the LFDMmodel predicts a lower abundance of low-velocity galaxiesin better agreement with observations than ΛCDM. Fur-thermore, since the suppression is related to the transi-tion redshift z t of DM formation (the later the transition,the larger the suppression of low-mass halo abundances),it is possible that setting z t to a value slightly smaller0 FIG. 7: Differential circular velocity function for the(27 . − Mpc) box simulations with 1024 particles for theΛCDM (red triangles) and Λ-LFDM (blue squares) models.Error bars are given by Poisson errors. The dotted straightline is the fitting function to the ΛCDM simulation given byEq. (5) in Ref. [41], while the dashed line includes the bary-onic corrections. The cyan solid lines denote the region of theobserved velocity function to within 15% error from Ref. [41]. than that assumed here may reproduce the observed ve-locity distribution up to statistical errors. VIII. CONCLUSIONS
The distribution of matter at small scales may provideinsights on the nature of dark matter particles and poten-tially constrain scenarios alternative to the standard colddark matter paradigm. In this work we have performeda first study of the nonlinear clustering of late-formingdark matter using a suite of high-resolution N-body sim-ulations. The main phenomenological feature of this sce-nario is the presence of damped acoustic oscillations atlarge wave numbers in the linear matter power spectrum.The scale where such a power suppression occurs de-pends on the epoch of the phase transition of a scalarfield coupled to relativistic particles that generates theDM particles. Without requiring an extreme fine-tuningof the microscopic model parameters, this may occur be-fore matter-radiation equality and imprint a pattern ofdamped oscillations at k (cid:38)
10 h Mpc − in the initialconditions. Because of this, LFDM models are indistin-guishable from standard CDM on large scales, while theyare constrained by probes of the small-scale clustering ofmatter.From the analysis of N-body simulations, we haveshown that LFDM could be a viable alternative to CDM, as it provides predictions of the nonlinear mat-ter power spectrum largely consistent with bounds fromhigh-redshift Lyman- α power spectrum measurements.It should be noted that the Lyman- α constraints haveso far been obtained for a restricted class of DM models(such as WDM) with a sharp cutoff in the linear powerspectrum. As LFDM models have spectra with con-siderably shallower slopes than WDM, it is very likelythat LFDM models that resolve ΛCDM inconsistencieson galactic and subgalactic scales would also be in betteragreement with Lyman- α measurements.We have evaluated the halo abundance; comparisonwith predictions from the standard ΛCDM model showsthat small-mass LFDM halos (10 < M [h − M (cid:12) ] < )are less abundant than their CDM counterparts. We havealso evaluated the corresponding differences in the circu-lar velocity function at z = 0 and shown that the LFDMscenario is in better agreement with observations of thenumber density of low-velocity galaxies than ΛCDM. Thenumber density of low-mass halos varies with redshift,with low-mass halo abundance up to ∼ z >
3. This indicates that halos in this massrange form later than in ΛCDM, thus leading to lowerdensities in the internal part of the halos as confirmedby the study of density profiles. On the other hand, theassembly of more massive halos occurs as in ΛCDM, andwe find no significant differences in their density profilesbetween the two cosmologies.These results confirm previous numerical studies ofnonstandard DM models characterized by initial linearmatter power spectra similar to that of LFDM. SinceLFDM behaves as a collisionless component, the lack ofcored profiles in Milky Way–like halos seems inherentlyrelated to the fact that in order to satisfy the Lyman- α constraints, the suppression of power in the initial powerspectrum must occur on wave number k >
10 h Mpc − .This suggests that only by including DM self-interactionswould it be possible to develop cored profiles. On theother hand, one should note that baryonic feedback doesplay a role at these scales. Given the differences in theabundance and assembly of small-mass halos in Λ-LFDMand ΛCDM, it will be of interest to investigate the dy-namics of baryons in the LFDM scenario. This may leavedistinct observational features in the formation of starsand galaxies which warrant further investigation. Acknowledgments
We would like to thank Matteo Viel and Jes´us Zavalafor useful comments and suggestions. This work wasgranted access to HPC resources of IDRIS through timeallocations made by GENCI (Grand ´Equipement Na-tional de Calcul Intensif) on the machine ADA No.x2014042287. We acknowledge support from the DIMACAV of the Region ˆIle-de-France. The research lead-ing to these results has received funding from the Euro-pean Research Council under the European Community’s1Seventh Framework Programme (FP7/2007-2013 GrantAgreement No. 279954).
Appendix: MULTIPLICITY FUNCTION
In the Press-Schechter approach [114] the halo massfunction can be written as dndM = ρ m M d ln σ − dM f ( σ ) , (A.1)where ρ m is the mean cosmic matter density and σ is theroot-mean-square fluctuation of the linear matter densityfield smoothed on a scale R ( M ) enclosing a mass M with σ ( M ) = 12 π (cid:90) dkk P ( k ) W [ k, R ( M )] , (A.2)where P ( k ) is the linear DM power spectrum and W ( k, R ) is the Fourier transform of the smoothing func-tion in real space (here we consider a top hat filter).The function f ( σ ) is the so-called “multiplicity function”which carries all information on the nonlinear gravita-tional processes that lead to the formation of halos.In the framework of the excursion set theory [115], thiscan be computed from the distribution of random walksfirst crossing a collapse density threshold that encodesthe nonlinear gravitational dynamics of matter collapse.It was shown in Ref. [116] that the multiplicity functionfor uncorrelated walks with a threshold barrier motivatedby the ellipsoidal collapse model is well approximated by a the Sheth-Tormen function (apart from a fudge factorcoefficient necessary to recover agreement with the N-body results) given by [108] f ST ( σ ) = A (cid:114) aπ δ c σ (cid:34) (cid:18) δ c √ aσ (cid:19) − p (cid:35) e − aδ c / σ , (A.3)where δ c is the linearly extrapolated spherical collapsedensity threshold. Since after the transition LFDM be-haves as a collisionless component, the spherical collapsedynamics is identical to that of CDM . At z = 0 we have δ c = 1 . δ c ( z ) by numerically solving the spherical collapseequations as in Ref. [111].We calibrate the ST parameters A , a and p to the nu-merical mass function at z = 0 from the (27 . − Mpc) volume simulation of the ΛCDM and Λ-LFDM model.We find A = 0 .
145 and a = 0 .
695 for both models, while p = 0 .
15 for ΛCDM and p = 0 . p is determined by theshape of the mass function at the low-mass end, whichin turn depends on the form of the ellipsoidal collapsethreshold (see also Refs. [117, 118] for an explicit rela-tion between the mass function and the parameters ofan ellipsoidal-collapse-inspired barrier). Hence, the largedeviations from unity of the ratio of the multiplicity func-tions at low masses shown in the left panel of Fig. 5 areclearly indicative of a departure of the nonlinear collapsedynamics of the Λ-LFDM model from that of the ΛCDM. This might not be the case for the ellipsoidal collapse modelwhich explicitly depends on the mass of the collapsing object, while the spherical collapse is independent of the mass. [1] G. Efstathiou et al. , Mon. Not. R. Astron. Soc. ,L29 (2002)[2] D. N. Spergel et al. , Astrophys. J. Suppl. Ser. , 175(2003)[3] M. Tegmark et al. , Phys. Rev. D. , 103501 (2004)[4] D. Clowe et al. , Astrophys. J. , L109 (2006)[5] P. A. R. Ade et al. (Planck Collaboration), Astron. &Astroph. , A1 (2011)[6] K. Freese, M. Lisanti, C. Savage, Rev. Mod. Phys. 85,1561 (2013)[7] R. Bernabei et al. , Eur. Phys. J. C , 39 (2010)[8] C. E. Aalseth et al. , Phys. Rev. Lett. , 141301 (2011)[9] G. Angloher et al. , Eur. Phys. J. C , 1971 (2012)[10] Z. Ahmed et al. , Phys. Rev. Lett. , 131302 (2011)[11] J. Angle et al. , Phys. Rev. Lett. , 051301 (2011)[12] D.S. Akerib et al. (LUX Collaboration), Phys. Rev.Lett. , 091303 (2014)[13] D. Hooper, J. Cosmol. Astropart. Phys. , 035 (2013)[14] B. Moore, Nature (London) , 629 (1994)[15] R. Kuzio de Naray, T. Kaufmann, Mon. Not. R. Astron.Soc. , 3617 (2011)[16] A. Klypin, A. V. Kravtsov, O. Valenzuela, F. Prada,Astrophys. J. , 82 (1999)[17] B. Moore et al. , Astrophys. J. , L19 (1999) [18] M. Boylan-Kolchin, J. S. Bullock, M. Kaplinghat, Mon.Not. R. Astron. Soc. , L40 (2011)[19] M. Boylan-Kolchin, J. S. Bullock, M. Kaplinghat, Mon.Not. R. Astron. Soc. , 1203 (2012)[20] J. S. Bullock, A. V. Kravtsov, D. H. Weinberg, Astro-phys. J. , 517 (2000)[21] A. J. Benson et al. , Mon. Not. R. Astron. Soc. , 177(2002)[22] R. S. Somerville, Astrophys. J. , L23 (2002)[23] M. Ricotti, N. Y. Gnedin, M. J. Shull, Astrophys. J. , 49 (2002)[24] M. Ricotti, N. Y. Gnedin, Astrophys. J. , 259 (2005)[25] J. I. Read, A. P. Pontzen, M. Viel, Mon. Not. R. Astron.Soc. , 885 (2006)[26] N. I. Libeskind et al. , Mon. Not. R. Astron. Soc. ,16 (2007)[27] A. V. Maccio et al. , Mon. Not. R. Astron. Soc. ,1995 (2010)[28] A. S. Font et al. , Mon. Not. R. Astron. Soc. , 1260(2011)[29] A. Pontzen, F. Governato, Mon. Not. R. Astron. Soc. , 3464 (2012)[30] R. Teyssier, A. Pontzen, Y. Dubois, J. I. Read, Mon.Not. R. Astron. Soc. , 3068 (2013) [31] S. Koposov et al. , Astrophys. J. , 279 (2008)[32] Q. Guo et al. , Mon. Not. R. Astron. Soc. , 101 (2011)[33] J. Wang, C. S. Frenk, J. F. Navarro, L. Gao, T. Sawala,Mon. Not. R. Astron. Soc. , 2715 (2012)[34] C. W. Purcell, A. R. Zentner, J. Cosmol. Astropart.Phys. , 007 (2012)[35] J. Pa˜ n arrubia, A. Pontzen, M. G. Walker, S. E. Ko-posov, Astrophys. J. , L42 (2012)[36] I. Ferrero et al. , Mon. Not. R. Astron. Soc. , 2817(2012)[37] S. Garrison-Kimmel, M. Rocha, M. Boylan-Kolchin, J.S. Bullock, J. Lally , Mon. Not. R. Astron. Soc. ,3539 (2013)[38] J. Zavala et al. , Astrophys. J. , 1779 (2009)[39] E. Papastergis, A. M. Martin, R. Giovanelli, M. P.Haynes, Astrophys. J. , 38 (2011)[40] E. N. Kirby, J. S. Bullock, M. Boylan-Kolchin, M.Kaplinghat, J. G. Cohen, Mon. Not. R. Astron. Soc. , 1015 (2014)[41] A. Klypin, I. Karachentsev, D. Makarov, O. Nasonova,arXiv:1405.4523[42] S. Garrison-Kimmel, M. Boylan-Kolchin, J. S. Bullock,E. N. Kirby1, Mon. Not. R. Astron. Soc. , 222 (2014)[43] E. Papastergis, R. Giovanelli, M. P. Haynes, F. Shankar,Astron. & Astrophys. , A113 (2015)[44] D. H. Weinberg et al. , arXiv:1306.0913[45] M. Kamionkowski, A. R. Liddle, Phys. Rev. Lett. ,4525 (2000)[46] A. R. Zentner, J. S. Bullock, Phys. Rev. D , 043003(2002)[47] A. R. Zentner, J. S. Bullock, Astrophys. J. , 49(2003)[48] S. Dodelson, L. M. Widrow, Phys. Rev. Lett. , 17(1994)[49] G. M. Fuller, A. Kusenko, I. Mocioiu, S. Pascoli, Phys.Rev. D , 103002 (2003)[50] K. Abazajian, Phys. Rev. D , 063513 (2006)[51] A. Boyarsky, J. Lesgourgues, O. Ruchayskiy, M. Viel,Phys. Rev. Lett. , 201304 (2009)[52] A. Boyarsky, J. Lesgourgues, O. Ruchayskiy, M. Viel,J. Cosmol. Astropart. Phys. , 012 (2009)[53] D. Boyanovsky, H. J. de Vega, N. G. Sanchez, Phys.Rev. D , 063546 (2008)[54] A. V. Maccio et al. , Mon. Not. R. Astron. Soc. ,1105 (2012)[55] M. R. Lovell et al. , Mon. Not. R. Astron. Soc. , 2318(2012)[56] M. R. Lovell et al. , Mon. Not. R. Astron. Soc. , 300(2014)[57] K. N. Abazajian, Phys. Rev. Lett. , 161303 (2014)[58] M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese,A. Riotto, Phys. Rev. D , 063534 (2005)[59] U. Seljak, A. Makarov, P. McDonald, H. Trac, Phys.Rev. Lett. , 191303 (2006)[60] M. Viel et al. , Phys. Rev. Lett. , 041304 (2008)[61] M. Viel, G. D. Becker, J. S. Bolton, M. G. Haehnelt,Phys. Rev. D , 043502 (2013)[62] A. Schneider, D. Anderhalden, A. V. Macci`o, J. Die-mand, Mon. Not. R. Astron. Soc. , L6 (2014)[63] D. Anderhalden et al. , J. Cosmol. Astropart. Phys. ,014 (2013)[64] D. N. Spergel, P. J. Steinhardt, Phys. Rev. Lett. ,3760 (2000)[65] J. L. Feng, M. Kaplinghat, H. Tu, H.-B. Yu, J. Cosmol. Astropart. Phys. , 004 (2009)[66] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, N.Weiner, Phys. Rev. D , 015014 (2009)[67] J. M. Cline, Z. Liu, G. D. Moore, W. Xue, Phys. Rev.D , 015023 (2014)[68] M. Vogelsberger, J. Zavala, A. Loeb, Mon. Not. R. As-tron. Soc. , 3740 (2012)[69] J. Zavala, M. Vogelsberger, M. G. Walker, Mon. Not.R. Astron. Soc. , L20 (2013)[70] M. Vogelsberger, J. Zavala, S. Christine, A. Jenkins,Mon. Not. R. Astron. Soc. , 3684 (2014)[71] A. Kamada, N. Yoshida, K. Kohri, T. Takahashi, J.Cosmol. Astropart. Phys. , 008 (2013)[72] C. Bœhm et al. , Mon. Not. R. Astron. Soc. , L31(2014)[73] M. R. Buckley, J. Zavala, F. Y. Cyr-Racine, K. Sigurd-son, M. Vogelsberger, Phys. Rev. D , 043524 (2014)[74] K. Sigurdson, M. Kamionkowski, Phys. Rev. Lett. ,171302 (2004)[75] P. H. G. Annika, C. E. Moody, M. Kamionkowski, Phys.Rev. D , 103501 (2010)[76] M.-Y. Wang, A. R. Zentner, Phys. Rev. D , 043514(2012)[77] M.-Y. Wang et al. , Mon. Not. R. Astron. Soc. , 614(2014)[78] M.-Y. Wang et al. , Phys. Rev. D , 123515 (2013)[79] S. Das, N. Weiner, Phys. Rev. D , 123511 (2011)[80] R. Fardon, A. E. Nelson, N. Weiner, J. High EnergyPhys. , 042 (2006)[81] S. Das, Journ. Phys.: Conf. Series , 012011 (2012)[82] D. J. E. Marsh, P. G. Ferreira, Phys. Rev. D , 103528(2010)[83] R. Hlozek, D. Grin, D. J. E. Marsh, P. Ferreira, Phys.Rev. D , 103512 (2015)[84] C. Boehm, A. Riazuelo, S. H. Hansen, R. Schaeffer,Phys. Rev. D , 083505 (2002)[85] F.-Y. Cyr-Racine, K. Sigurdson, Phys. Rev. D ,103515 (2013)[86] R. J. Wilkinson, J. Lesgourgues and C. Boehm, J. Cos-mol. Astropart. Phys. , 026 (2014)[87] A. Sarkar, S. Das, S.K. Sethi, J. Cosmol. Astropart.Phys. , 004 (2015)[88] Planck Collaboration: P. A. R. Ade et al. , Astron. As-trophys. , A16 (2014)[89] Planck Collaboration: P. A. R. Ade et al. , Astron. As-trophys. , A15 (2014)[90] A. Lewis, A. Challinor, A. Lasenby, Astrophys. J.
473 (2000)[91] B. A. Reid et al. , Mon. Not. R. Astron. Soc. o z-Cuartas, Mon. Not. R. Astron. Soc. , 882 (2013)[93] A. A. Staboninsky, JETP Lett. , 489 (1992)[94] J. Lesgourgues, D. Polarski, A. A. Starobinsky, Mon.Not. Roy. Astron. Soc. , 769 (1998)[95] R. Teyssier, Astron. & Astrophys. , 337 (2002)[96] T. Guillet, R. Teyssier, J. Comput. Phys. , 4756(2011)[97] S. Prunet et al. , Astrophys. J. Suppl. Ser. , 179(2008)[98] M. Crocce, S. Pueblas, R. Scoccimarro, Mon. Not. R.Astron. Soc. , 369 (2006)[99] K. Heitmann et al. , Astrophys. J. , 104 (2010)[100] C. Orban, Phys. Rev. D , 023509 (2014) [101] Y. Rasera et al. , Mon. Not. R. Astron. Soc. , 1420(2014)[102] F. Roy, V. Bouillot, Y. Rasera, Astron. & Astrophys. , A13 (2014)[103] M. Davis, G. Efstathiou, C. S. Frenk, S. D. M. White,Astrophys. J. , 371 (1985)[104] M. Gotz, J. Sommer-Larsen, Astrophys. Space Sci. ,415 (2002)[105] M. Gotz, J. Sommer-Larsen, Astrophys. Space Sci. ,341 (2003)[106] J. Wang, S. D. M. White, Mon. Not. R. Astron. Soc. , 93 (2007)[107] S. Agarwal, P.-S. Corasaniti, Phys. Rev. D , 123509(2015)[108] R.K. Sheth, G. Tormen, Mon. Not. R. Astron. Soc. ,119 (1999)[109] J. Tinker et al. , Astrophys. J.
709 (2008)[110] M. Crocce, P. Fosalba, F. J. Castander, E. Gazta˜ n aga, Mon. Not. R. Astron. Soc. , 1353 (2010)[111] J. Courtin et al. , Mon. Not. R. Astron. Soc. , 1911(2011)[112] C. Schultz, J. O˜ n orbe, K. N. Abazajian, J. S. Bullock,Mon. Not. R. Astron. Soc. , 1597 (2014)[113] S. Trujillo-Gomez, A. Klypin, J. Primack, A. J. Ro-manowsky, Astrophys. J. , 16 (2011)[114] W.H. Press, P. Schechter, Astrophys. J. , 425 (1974)[115] J. R. Bond, S. Cole, G. Efstathiou and G. Kaiser, As-trophys. J. , 440 (1991).[116] R. K. Sheth, H. J. Mo and G. Tormen, Mont. Not. R.Astron. Soc. , 1 (2001)[117] P.S. Corasaniti, I. Achitouv, Phys. Rev. Lett. ,241302 (2011)[118] A. Lapi, P. Salucci, L. Danese, Astrophys. J.772