On the Structure of Dark Matter Halos at the Damping Scale of the Power Spectrum with and without Relict Velocities
We report a series of high–resolution cosmological N-body simulations designed to explore the formation and properties of dark matter halos with masses close to the damping scale of the primordial power spectrum of density fluctuations. We further investigate the effect that the addition of a random component, , into the particle velocity field has on the structure of halos. We adopted as a fiducial model the warm dark matter (WDM) cosmology with a non–thermal sterile neutrino mass of 0.5 keV. The filtering mass corresponds then to . Halos of masses close to were simulated with several million of particles. The results show that, on one hand, the inner density slope of these halos (at radii the virial radius R) is systematically steeper than the one corresponding to the Navarro-Frenk-White (NFW) fit or to the cold dark matter counterpart. On the other hand, the overall density profile (radii larger than 0.02 R) is less curved and less concentrated than the NFW fit, with an outer slope shallower than -3. For simulations with , the inner halo density profiles flatten significantly at radii smaller than 2–3 kpc (). A constant density core is not detected in our simulations, with the exception of one halo for which the flat core radius is . Nevertheless, if “cored” density profiles are used to fit the halo profiles, the inferred core radii are , in rough agreement with theoretical predictions based on phase–space constrains, and on dynamical models of warm gravitational collapse. A reduction of by a factor of 3 produces a modest decrease in core radii, less than a factor of 1.5. We discuss the extension of our results into several contexts, for example, to the structure of the cold DM micro–halos at the damping scale of this model.
Subject headings:cosmology:dark matter — galaxies:halos — methods:N–body simulations
The nature of dark matter (DM) is one of the most intriguing and fundamental problems in cosmology and particle physics. The standard hypothesis assumes that dark matter is made of non–baryonic collisionless elemental particles that become non–relativistic very early in the history of the Universe (cold). This minimal scenario, named Cold Dark Matter (CDM), has successfully explained the observed structure of the universe at large scales, like the two–point correlation function of galaxies and the Cosmic Microwave Background Radiation (CMBR) anisotropies (see for recent results Springel et al., 2006; Spergel et al., 2006). Confrontation of model predictions with obervations turns out to be more complicated at galactic scales, because non–linear dynamics and baryonic processes may distort considerably the underlaying DM distribution. Thus, the predictions of the CDM scenario and its variants at the scale of galaxies, are an active subject of study.
Two of the most controversial CDM predictions are: (i) the large abundance of subhalos in galaxy–sized halos (Kauffmann et al., 1993), and (ii) the cuspy inner density profile of dark halos (e.g., Navarro et al., 2004; Diemand et al., 2004). Based on comparison with observations, it has been suggested that both predictions may indicate a flaw of the CDM scenario (Klypin et al., 1999a; Moore et al., 1999a; Moore, 1994; de Blok et al., 2001, Gentile et al. 2005,2007). However, these comparisons might be biased by astrophysical processes in action during galaxy assembly and evolution, as for example the inhibition of star formation in small subhalos (Bullock, Kravtsov, & Weinberg, 2000; Benson et al., 2003; Governato et al., 2007) or the halo core expansion due to energy or angular momentum transfer from dark/baryonic structures (Ma & Boylan-Kolchin, 2004; El-Zant et al., 2004; Weinberg & Katz, 2007), but see also (Colín et al., 2006; Ceverino & Klypin, 2007; Sellwood, 2007). It has also been showed that the disagreements may be a consequence of systematics in the observational inferences, (e.g., Simon & Geha, 2007; Rhee et al., 2004; Hayashi & Navarro, 2006; Valenzuela et al., 2007).
It is also possible that slight modifications to the CDM particle properties could solve or ameloriate the mentioned potential problems if they persist (e.g., Spergel & Steinhardt, 2000; Colín et al., 2002). As the precision of observations and the control of systematics improve, the confrontation with model predictions opens a valuable window for constraining the dark matter properties. Among the “slight” modifications with the CDM scenario is the introduction of warm particles, instead of cold ones (Warm Dark Matter, hereafter WDM). WDM implies two more degrees of freedom in comparison to CDM: (i) a damping (filtering) of the power spectrum at some intermediate scale, , due to the relativistic free streaming, and (ii) some primordial random velocity in the dark particles, (for CDM, the mass corresponding to , , is comparable to planet masses, e.g., Diemand et al., 2005; Profumo et al., 2006, and ).
Previous numerical studies have shown that WDM can be very effective in reducing the amount of (sub)structure below the filtering mass (e.g., Colín, Avila-Reese, & Valenzuela, 2000; Avila-Reese et al., 2001; Bode et al., 2001; Knebe et al., 2002). In addition, both the peculiar dynamical formation history of these halos and its random can also impose an upper limit on the phase space density, potentially producing an observable core of constant density (Hogan & Dalcanton, 2000; Avila-Reese et al., 2001). If has a thermal origin its amplitude is linked directly to the mass of the WDM particle (Hogan & Dalcanton, 2000); notice, however, that the amplitude of the random velocities may depend on other physical factors not directly related to the particle mass. This is the case, for instance, of gravitinos produced non–thermally by late decays of the Next to Lighest Supersymmetric Particle (NLSP) (Feng, Rajaraman & Takayama, 2003; Strigari et al., 2007, see for a recent review Steffen 2006). Thus, an exploration of the effect of a random velocity component independently of the dark matter particle mass seems to be necessary.
Predictions of the core radius in WDM halos have been computed assuming a King profile (Hogan & Dalcanton, 2000) or a subclass of the Zhao (1996) profiles (Strigari et al., 2006, herafter S2006). It is still unknown which estimates give the more accurate value, yet these predictions are necessary for comparison with observations (see e.g., S2006; Gilmore et al., 2007). Unfortunately, the predicted core radii for masses of the most popular WDM particle, the sterile neutrino, allowed by observational constraints (, see §5 for references) are below the resolved scales in current simulations.
The structure of WDM halos of scales below the filtering radius might be different from their CDM counterparts not only in the central parts but also in their overall mass distribution. This is somehow expected because the assembly history of these halos is different from the hierarchical one. Besides, they form later and have concentrations lower than those ones derived for CDM halos (Avila-Reese et al., 2001; Bode et al., 2001; Knebe et al., 2002). However, it is controversial whether the shape of the density profile differs systematically from the corresponding CDM one (see for different results e.g., Huss et al., 1999; Moore et al., 1999b; Knebe et al., 2003). On the other hand, WDM halos of mass close to (with set to 0) can be thought of as scaled-up versions of the first microhalos in a CDM cosmology. The earliest collapse of CDM microhalos is a subject of considerable current interest (e.g., Diemand et al., 2005; Gao et al., 2005).
In this paper we explore by means of numerical simulations the two questions mentioned above: the overall structure of dark halos with masses close or just below the cutoff mass in the power spectrum of fluctuations, and the inner density profile of these halos when a random velocity is added to the particles. For our numerical study we use the truncated power spectrum corresponding to a non–thermal sterile neutrino of keV (). We initially neglect the particle random velocity (=0), and later we will consider two values of that cover the range of velocities corresponding to thermal and non–thermal keV WDM particles. It is important to remark that our goal is not to study a specific WDM model but to explore in general the effects on halos of the power spectrum truncation and the addition of random velocities to the particles.
The structure of the paper is as follows. In §2 we describe the cosmological model that we use for our investigation : a WDM model with a filtering radius at the scale of Milky Way size halos. Halos of these scales were resimulated with higher resolutions, first without adding the corresponding random velocity, and then with adding this velocity component to the particles. Details of the numerical simulations carried out in this paper are given in §3. The results from our different simulations are presented in §4. Section §5 is devoted to a discussion of the results and their implications. Finally, in §6 we present the main conclusions of the paper.
2. The cosmological model
The general cosmological background that we use for our numerical simulations corresponds to the popular flat low–density model with , and (the Hubble constant in units of ).
For the experiments designed to explore the structure of dark matter halos with masses close to , we adopt an initial power spectrum corresponding to a non–thermal sterile neutrino of keV. Even if this WDM model is ruled out by observations (see §5 for references), it is however adequate for the purposes stated in the Introduction. As is shown below, the filtering mass corresponding to this WDM particle is of the order of Milky Way–size halos, namely the halos that we are able to follow with high resolution in a cosmological simulation. The high resolution of the simulations avoids that the formation of halos with a mass scale near to will be dominated by discreteness effects (see §3.1). Besides, for the resolution that we attain, we expect to resolve the inner regions where a flattening in the halo is predicted density profile is predicted for the case when is introduced.
We use here the transfer function for the non–thermal sterile neutrino derived in Abazajian (2006a). The WDM power spectrum is then given by
and is the CDM power spectrum given by
Klypin & Holtzman (1997). This fit is
in excellent agreement, at the scales of interest, with the power spectrum
obtained with linger, which is contained inside the
where , , , , , and . The power spectrum is normalized to , a value close to that estimated from the third-release of WMAP (Spergel et al., 2006). Here is the rms of mass fluctations estimated with the top-hat window of radius .
As in Avila-Reese et al. (2001), we have defined the free–streaming (damping) wavenumber, , as the for which the WDM transfer function decreased to 0.5, and compute the corresponding filtering mass in the linear power spectrum as
where is the present–day mean density of the universe. The filtering wavelength of a non–thermal sterile neutrino of is , which corresponds to a filtering mass . The random component was linearly added to the velocities calculated with the Zel’dovich approximation at the onset of the simulation.
A number of Milky Way–size halos are simulated with the WDM power spectrum given by eq. (1) and with =0. In this way we isolate the effects of the power spectrum filtering on the structure of the halos. Later, the same halos are resimulated but adding random velocities to the particles. We approximate the shape of the particle phase space distribution function (DF) with the corresponding thermal one; i.e, we use a thermal equilibrium Fermi-Dirac DF. This is a good approximation for a non–thermal sterile neutrino of 0.5 keV (Abazajian, 2006a). Further, as a first approximation, we assume that the amplitude of the DF of our non–thermal sterile neutrinos is equal to the one corresponding to their thermal partners with =0.5 keV (Hogan & Dalcanton, 2000). With this assumption for DF shape and amplitude, the neutrino random velocity at , the initial redshift in our simulations, is .
The particle velocity dispersion in the linear regime decreases adiabatically with the expansion. The convention is to define in terms of the value extrapolated to the present epoch, ; therefore, for our case km/s. We should note that the rms velocity amplitudes of sterile neutrinos and thermal particles of the same mass can not be the same, but larger because the former become non–relativistic later than the latter. A rough calculation shows that for the 0.5 keV sterile neutrino, should be approximately three times higher than for a 0.5 keV thermal WDM particle, i.e. km/s. We resimulate halos with both values of the random velocity, and 0.3 km/s, having in mind that our main goal is not to study a particular WDM model but to explore the structure of the dark halos with different initial conditions .
3. Numerical Simulations
A series of high resolution simulations of Milky Way size galactic halos are performed using the Adaptive Refinement Tree (ART) N-body code (Kravtsov, Klypin, & Khokhlov, 1997) in its multiple mass scheme (Klypin et al., 2001). The ART code achieves high spatial resolution by refining the root grid in high-density regions with an automated refinement algorithm.
In all of our experiments we first run a low–mass resolution (LMR) simulation
of a box of Mpc on a side and or
particles. We then select a spherical region centered on a Milky
Way size halo of radius times the virial
We started the simulations at because the power at the Nyquist frequency at this redshift is in the linear regime. Note that spurious noise can influence the evolution of simulations that include if they start too early. The noise might be particularly important when the generated have amplitudes comparable to those of the Zel’dovich peculiar velocities. In order to show this, we started a simulation at and evolved it up to . At the velocities are on average 2.5 times greater than at . Figure 1 shows the measured initial power spectrum of the simulation started at (solid line) and the power spectrum for the simulation started at and measured at (squares). The latter simulation developed spurious power at frequencies higher than about 0.6 Mpc. We run another simulation started at , but with no addition of a relic velocity component. In this case we did not detect the evolution of spurious power by . In order to make sure that this spurious noise does not appear when the simulation is started at we repeated the experiment but now the onset of the simulation is set at and the power spectra are measured at . Unlike the previous case, no differences between the power spectra were detected. In other words, as far as the initial power spectrum is concerned it does not matter if the simulation is started at o .
Concern may arise about the structure of halos simulated in a relatively small Mpc box, specially in a WDM cosmology, where there is a scale below which the power spectrum exponentially drops to zero. Avila-Reese et al. (2001) discussed this potential issue and concluded that in order to be confident about the simulated halo structure, a box size greater than the filtering length should be used. For our simulations is 2.5 larger than . In any case, we also experimented with other box sizes, namely, 15Mpc and 20Mpc (not shown in the Tables), for the case, and found results similar to those reported here for the 10Mpc box.
The bound density maxima (BDM) group finding algorithm (Klypin et al., 1999b), or a variant of it (Kravtsov et al., 2004), is used to locate the halos in the simulations and to generate their density profiles. The BDM finds positions of local maxima in the density field smoothed at the scale of interest and applies physically motivated criteria to test whether a group of particles is a gravitationally bound halo.
Aside from those halos shown in Table 1, for the halo D with =0 we have also run a very high resolution simulation with about 31 million particles in the high resolution zone. This halo was taken from the same 10 Mpc box run but with particles in the LMR mode. As far as we know this is the highest resolution simulation of a halo run in a WDM cosmology. The same halo was also simulated with less resolution for a convergence test. The parameters of the sequence of halos D are resumed in Table 2.
The simulations presented here differ in several aspects from previous WDM simulations. First, it should be emphasized that our aim rather than discussing a specific WDM model is to explore the influence of the truncation of the power spectrum and/or the addition of random velocities upon the structure of dark halos of masses close to the truncation scale. For this aim we need to simulate (a) halos with very high resolution, and (b) halos with masses close to . The halos simulated in Avila-Reese et al. (2001) had several times less particles than the best–resolved halos presented here and the aims in that paper focused in exploring general halo properties for a concrete WDM model. Other papers aimed to study the properties of WDM halos (Bode et al., 2001; Knebe et al., 2002; Busha et al., 2006) focused more in the statistical aspects than in details of the inner halo structure; therefore, the halos in these papers had resolutions much lower than those attained here. The properties of the WDM halos simulated here are in general agreement with previous findings; for example, their concentrations are systematically lower (Avila-Reese et al., 2001; Eke et al., 2001; Bode et al., 2001) and they form later (Knebe et al., 2002; Busha et al., 2006) than the corresponding CDM halos.
3.1. Discreteness effects
One of the motivations of this paper is to investigate the structure of well resolved halos with masses close or below the damping (truncation) scale in the power spectrum, . The origin of these halos is controversial. Halos with masses close to (truncation halos) could be formed by a quasi–monolithic collapse of filaments of size (e.g., Avila-Reese et al., 2001). They could also be just the result of an incomplete collapse, highly deviated from the spherical–symmetric case, of originally larger structures assembled hierarchically (Busha et al., 2006). On the other hand, it has been suggested that halos with masses considerably less than form by fragmentation of the shrinking filaments of size (e.g., Valinia et al., 1997; Avila-Reese et al., 2001; Bode et al., 2001; Götz & Sommer-Larsen, 2003; Knebe et al., 2003). However, it is also known that the filaments in Hot Dark Matter simulations that start from a cubic lattice break up into regularly spaced clumps, which reflect the initial grid pattern. Therefore, some of these halos seen in WDM simulations could be spurious, product of discreteness effects. Recently, Wang & White (2007) have shown that this artifact is present even for a glass–like initial particle load (White, 1996).
As Wang & White (2007) show, halos of masses smaller than a given effective fraction of , which depends on the resolution of the simulation, will be spurious. We selected the WDM model (§2) and the number of particles in the simulations in such a way that the halos studied here can neiher be fake nor affected by discreteness effects. Following Wang & White (2007), in our model with effective number of particles, only structures with masses lower than 1.3 are candidate to be spurious. For the models with and effective number of particles, the masses of structures triggered by the initial grid spacing are and 1.04 , respectively. For our WDM model, the first halos to form have masses , which are far from the mentioned effective resolution limits.
We first present the results of our WDM simulations with . In this way, we explore the inner structure of halos close to the truncation mass , which could be scaled-up versions of the first CDM microhalos (of masses for a neutralino mass of 100 GeV). Afterwards we present the results of the same simulations but introducing random velocities to the particles with two different amplitudes, and 0.3 km/s (see §2). The main goal of the last simulations is to explore the predicted flattening in the halo inner density profile produced by the addition of a random component to the particle velocities(see section 1). Table 1 resumes the main properties of the resimulated halos, which were selected to have masses close to the truncation mass of the initial power spectrum.
4.1. The Structure of Halos at the scale of Damping
Figures 2 and 3 show the spherically averaged density profiles measured for the halos of masses around the power–spectrum filtering mass, , and with =0. The first panel of Fig. 2 shows halos A256 and B256 (the latter was shifted by in the log); the second and third panels show halos C512 and E512, respectively. In Fig. 3, the same halo but simulated with 4 different resolutions is shown. In both Figs., the thin dashed lines are the best Navarro–White–Frenk (NFW, Navarro, Frenk, & White, 1997) fits to the showed density profiles. In the lower panels, the residuals of the measured density profile and the NFW fit are plotted with the same line coding as in the corresponding upper panels.
Halos A256 and B256 in the left panel and C512 and E512 in the second and third panels of Fig. 2, were simulated with formal spatial resolutions and 0.152kpc, respectively (see Table 1). Previous convergence studies for CDM halos simulated with the ART code have shown that the innermost halo density is reliable only for radii larger than four times and containing more than 200 particles (Klypin et al., 2001). For all the density profiles shown in Figs. 2 and 3, the innermost plotted point corresponds to radii larger than by and they contain more than 200 particles. The convergence analysis that we have carried out for our WDM halos suggests that instead of , the innermost radius should be close to , .
Figure 3 compares the density profiles of halo D, which was re–simulated with four different resolutions separated each by a factor of eight in the particle mass (see Table 2). As can be seen, convergence is achieved at about . The arrows in Figs. 2 and 3 indicate for the corresponding simulations. In Fig.3 , the solid–line arrow is for the highest–resolution simulation ( effective number of particles), while the dashed–line arrow is for the simulation; the NFW fit is shown only for these two cases.
The inner density profiles of all halos simulated here are systematically steeper than the corresponding fitted NFW law. The slopes of the profiles at the virial radii, R, span a range from to . For comparison, the slope of the density profile of a typical LCDM halo of at R is (a NFW profile was used with the corresponding concentration given by Bullock et al., 2001, and re-scaled to ). At radii smaller than R, the slopes tend to become shallower but the halos are still denser and slopes steeper than the corresponding NFW fit up to the resolution limits (, see Figs. 2 and 3).
The overall density profile shapes of our halos are also somehow different from the NFW function. For radii larger than R, the profiles tend to be in general slightly less curved than the NFW model. This is why the residuals shown in Figs. 2 and 3 indicate a systematical defect at intermedium radii and then an excess at the outer radii. The outer slopes are .
In summary, the density profiles of the halos simulated here, with masses close or below the filtering mass, have a shape slightly flatter than the NFW law for radii R and slopes significantly steeper at radii R. Nevertheless, each profile is different. The profile of halo E512 has minimal deviations from the NFW function, while the profile of halo D512 significantly deviates from this function.
Table 1 details the main properties of the halos studied here. In columns 6 to 14 are reported respectively the virial mass, M, the maximum circular velocity, , the c and NFW concentration parameters, the average density within 1% the virial radius, and three core radii estimated by different criteria (the latter quantities apply only to halos simulated with , see §§4.2). The NFW and c concentrations are defined respectively as the ratios between R and the NFW scale radius, and between R and the radius where 1/5th of M is contained (Avila-Reese et al., 1999). As found in previous results, halos of masses below the truncation mass in the power spectrum tend to be less concentrated than CDM halos of similar masses (Avila-Reese et al., 2001; Bode et al., 2001; Knebe et al., 2002). We have simulated some CDM () halos of masses and measured NFW and c concentrations around 8–15 and 6.0–11.0, respectively, to be compared with the values given in Table 1 for the =0 cases.
4.2. The inner structure of halos simulated with rms velocity added to the particles
We have rerun the halos presented in the previous sub–section but now introducing a random velocity component, . As explained in §2, the particle DF used corresponds to a Fermi–Dirac function. Regarding the amplitude, we use two values: and 0.3 km/s. These values are for WDM particles of thermal origin and for non–thermal sterile neutrinos, respectively, in both cases with =0.5 keV. Our goal is to explore whether the inner structure of the halos becomes affected significantly or not by adding .
We first present results for the case km/s, and then explore how the inner halo structures change when is increased from 0.1 to 0.3 km/s, a more appropiate value for the 0.5 keV sterile neutrino model used to generate the initial power spectrum of the simulations (see §2). For halo D, we have simulated with the case. Unfortunately, due to limitations in our computational resources, it was not possible to run the simulation with particles for .
In Fig. 4 we present spherically–averaged profiles of different properties for all the simulated halos without and with added. In the first and third columns two halos (A and B, and Db and D, respectively) are presented but the latter ones are down shifted for clarity. The upper panels of Fig. 4 show the density profiles of our halos without (black solid lines), and with km/s (red dot–dashed line) and 0.3 km/s (magenta dashed line). In the medium panels are plotted the corresponding three–dimensional velocity dispersion profiles, , and the lower panels show the corresponding coarse–grained phase–space density profiles, . As in Fig. 2, the arrows indicate the strong resolution limit radius of the simulations. For halo D (third column), the simulations with were carried out with particles and using two different random seeds for the particle DF calculation (see below).
To some degree, all the simulated halos have been affected in their inner regions by the injection of initial random velocities to their particles. For the less resolved halos A, B (left panels), the deviations at of the inner density profiles from the profiles obtained in the simulations with =0 are yet marginal, but in the expected direction. For the halos C, D, Db, and E simulated with particles, the deviations down to are significant: the density profiles systematically flatten with respect to the corresponding =0 cases. The radii at which the density profiles of halos with start to deviate (flatten) from the ones without , are , well above the resolution limit of the simulations.
In a second series of experiments, we have resimulated halos C, Db and E ( particles) with three times larger, i.e. km/s. The profiles corresponding to these simulations are plotted in Fig. 4 with dashed (magenta) lines. The flattening of the inner density profiles is clearly more pronounced than in the simulations with km/s.
The innermost density profiles actually vary from halo to halo. Again, halo E is the less affected not only by the damping in the power spectrum, but also by the injection of , and halo D is the most affected by both effects (lower dot–dashed curve in the corresponding panel). The latter actually shows a “true” flat core already at . Since the core of this halo is too different with respect to the other ones, we decided to explore if large difference can be explained by a rare fluctuation in the random procedure of particle velocity assignment. Thus, the same halo D512 was resimulated with a different seed in the random number generator used to draw the particle velocities. The upper curves in the third column panels of Fig. 4 correspond to the profiles for this halo, called D512b. The inner density profile of this halo is not too different from the profiles of the other halos, though it remains as the flatest one among all the simulated halos with km/s.
In Fig. 5 we attempt to fit different functions to the density profiles of halos C, D, and E ( particles) with and 0.3 km/s. A general function to describe density profiles of cosmic objects was proposed by Zhao (1996):
The NFW profile corresponds to . This function does not provide a good description for the profiles of our halos with , in particular in the inner regions. We have fitted the halo profiles to the NFW in order to obtain an estimate of the c concentration given in Table 1.
Strigari et al. (2006) suggested a “cored” density profile in order to derive constraints on the size of a possible shallow core in the halo of the Fornax dwarf spheroidal galaxy (Goerdt et al., 2006; Sánchez-Salcedo et al., 2006). This profile is described by eq. (5) with , and they define the core radius, r as the radius where the inner slope, , reaches the value of . Thus, . The dot–dashed curves in Fig. 5 show the best fits using the S2006 profile. As can be seen, this profile does not describe well the density profiles of the WDM halos in the simulations. Note that characterizes the sharpness of the change in logarithmic slope. As already seen in Figs. 2 and 4, the profiles of our halos tend to be less curved than the usual NFW profile. Therefore, values of smaller than 1 should be used instead of larger than 1. We have obtained a reasonable description of our WDM profiles with . The (magenta) dashed curves in Fig. 5 are the best fits with these profiles. The core radius, as defined above, is in this case . Finally, we have also tried fits to the S2006 function but taking into account only the central halo regions, up to kpc. The fits are shown in Fig. 5 with (green) dotted curves. The fit is specially good for the halo D512 and those with (0)=0.3 km/s. Columns 12 to 14 in Table 1 report the values of r obtained with the three different fits.
As to the 3–D velocity dispersion profiles of the simulated halos, they do not differ significantly between the cases with and without , the exception being halo D. In the innermost regions, is similar or higher for halos with than for those without . The larger differences are for halo D, which after introducing random velocities to the particles produce a relatively hot core.
Finally, from the measured density and dispersion velocity profiles, we calculate the coarse–grained phase–space density profiles, . As seen in the lower panels of Fig. 4, excepting the innermost regions, the profiles are well described by a power law with close to that obtained for CDM halos by Taylor & Navarro (2001). For the inner regions, the profile of the halos simulated without tends to steepen, specially in halos D512 and C512. The opposite happens for the halos simulated with adding , the inner profile tends to be flatter as is larger.
5.1. Robustness of the results
The halos studied here are on one hand among the first virialized structures to form in our simulations, and on the other their assembling process started relatively late in the universe, between and (for discussions about the mass assembling process of halos with masses close to the damping scale in the power spectrum see e.g., Moore et al., 2001; Avila-Reese et al., 2001; Bode et al., 2001; Knebe et al., 2002; Busha et al., 2006). Because of the late collapse, one might argue that the halos studied here are not relaxed and therefore is not suprising to have the deviations from the NFW profile, as reported in §4.1 for the experiments with =0. We have followed the evolution of some of our halos for more than a Hubble time (scale factor ), finding a negligible evolution in the density profiles since . For example, halo D512 was run until (18.4 Gyr). The density profile of the halo at this epoch is practically the same as at (13.7 Gyr, see Fig. 4). As discussed in previos studies (see the references above), the collapse of halos of scales close to the damping scale seems to be quasi–monolithic (though highly non spherical). Thus, in regions that remain relatively isolated, as in the case of the halos selected for our study, halos suffer low mass accretion, and their structures remain almost unaltered since the initial collapse.
We checked that the effects on the structure of halos reported in §4 are systematic by running for a WDM model –not shown in Table 1– the corresponding CDM simulation, using the same random phases and changing only the initial power spectrum. We found that the density profile of the CDM halo is well fitted by the NFW function, while the corresponding WDM halo present the systematic deviation already seen in Fig. 2. Figure 6 compares the density and circular velocities profiles for the halo in question in its two versions, CDM and WDM without random velocities. Since the WDM halo ends up at with a slightly lower mass than the CDM halo we correct the profiles so as to make the comparison at a fixed mass (). Notice, in particular, that the inner density profile of the WDM halo is indeed steeper than their CDM counterpart.
Concerning the resolution limit in our simulations, based on the convergence study carried out for halo D (see Fig. 2), we find that a strong limit is , but resolution might be still acceptable for radii slightly larger than , a value suggested previously for CDM halos in “equilibrium” simulated with the ART code (Klypin et al., 2001). With a resolution limit at , our simulations allow to resolve the inner structure of halos down to kpc for the runs and down to kpc for the run. These radii correspond respectively to % and 0.25% virial radius in our halos.
Finally, it is important to recall that the masses of the the WDM halos analyzed in our simulations, are well above from the mass scale affected by discreteness effects, like the spurious formation of structures and substructures due to the initial grid pattern (see §3.1).
5.2. Do soft cores form in WDM halos?
Early structure formation studies based on a WDM cosmology considered particles originated in thermal equilibrium. For this case, both and depend only on the particle mass . The smaller , the larger and . Controlled numerical simulations of isolated halos showed that in order to produce “observable” soft cores, the amplitude of should be several times higher than the values corresponding to thermal WDM particles of masses keV (Avila-Reese et al., 2001). Particle masses smaller than keV are not allowed by the constraints on satellite galaxy abundances as well as by the Ly– power spectrum alone or combined with CMBR and large scale structure data. The Ly– power spectrum is the strongest of the constraints. For the non–thermal sterile neutrino, it places a limit on it mass at keV (Seljak et al., 2005; Viel et al., 2006; Abazajian, 2006b). For thermal WDM particles, the observational constraints give a limit of keV (Narayanan et al., 2000; Viel et al., 2005; Abazajian, 2006b), while a distinct analysis, using different simulations provide a stronger limit, keV (Seljak et al., 2006).
We may estimate the expected flat core radii of WDM halos for thermal particles in the mass range keV by using the approximation given in Avila-Reese et al. (2001, their eq. 13). This approximation is based on the monolithic collapse of halos with non–negligible particle random velocities before the collapse. For thermal WDM particles of masses 2 and 0.5 keV, the of a perturbation are and 1.9, respectively, while the corresponding to these masses are 0.1 and 0.015 km/s. Therefore, the expected core radii are and 210 pc.
The value of for a non–thermal sterile neutrino of 0.5 keV is approximately three times larger than the corresponding to the thermal particle of the same mass. Therefore, for this case pc. So, the resolutions that we may attain in our simulations of WDM halos for the =0.5 keV sterile neutrino are already close to these estimates of the flat core radii.
Recently, alternative particle models, like super–WIMPS, were proposed. The dark particles in these models may acquire random velocities non–thermally, for example, through the decay process of NLSP particles (e.g., charged Sleptons into gravitinos, Feng, Rajaraman & Takayama, 2003, and see for more references Steffen 2006). In these cases, does not depend directly on the damping scale of the linear power spectrum. However, an extra parameter is introduced, the decaying epoch. Strategies to constraint this parameter using astrophysical observations have been proposed (Feng, Rajaraman & Takayama, 2003; Strigari et al., 2007).
Some cosmological models with super–WIMP particles may be in agreement with constraints based on structure formation, specially the Ly- forest, and still allow for relatively large velocity dispersions, able to ameloriate the potential problems of CDM at small scales. This is the case of the so called Meta–dark matter models that consider the late decay of neutralinos into gravitinos. These models preserve a power spectrum similar to CDM models and at the same time set a phase space limit in the innermost structure of dark halos due to the injection of random velocities to the particles (Strigari et al., 2007; Kaplinghat, 2005). However, as mentioned in the Introduction, it was still an open question of how efficient the introduction of is for producing significant effects on the inner structure of simulated dark halos. One of the goals of the present paper was just to expore this question by means of numerical simulations able to resolve the WDM halos down to R.
The results presented in §4.2 show definitively that the addition of flattens the inner density profiles of WDM halos, and more as is higher. Differences in density profiles of halos simulated with and those with , start to be evident at radii of 2–3 kpc() for our halos. The halo average inner density measured at 0.01R, decreases on average by factors 1.5 and 2.5 for the models with and 0.3 km/s, respectively (column 11 in Table 1). Nevertheless, we notice that the concentration of the halos seems not to change significantly by the introduction of .
In general, the NFW function does not describe well the inner density profiles of our simulated halos; their profiles are much shallower than as can be seen in Fig. 5. However, is the size of the random velocity effect as theoretically expected? The theoretical predictions are based on the existence of an upper limit in the fine–grained phase–space density due to the collisionless nature and finite relict velocity dispersion of particles. This upper limit, , implies that the halo density profile must saturate and form a constant–density core (Hogan & Dalcanton, 2000).
For the random velocities used in this study ( and 0.3 km/s, corresponding respectively to thermal and non–thermal 0.5 keV sterile neutrinos), and Mpc/(km/s). According to Hogan & Dalcanton (2000), the core radius produced by the phase–space packing scales as , where is the asymptotic circular velocity for the assumed non–singular isothermal sphere (their eq. 18). This implies that more massive halos have smaller, more tightly bound cores. Applying the same equation for keV and , the core radius for the thermal particle ( Mpc/(km/s)) is 85 pc, while for the sterile neutrino ( Mpc/(km/s)) is 450 pc.
In any case, our results can only marginally test these estimates. The resolution limit in our simulations with non–zero is in between and 1.7 kpc. If the inner density profile response to the introduction of is gradual, one expect to see yet some effects at these radii, and this happen to be the case.
A way to attempt to infer (extrapolate) the sizes of possible flat cores in the simulated halos is by fitting the measured density profiles to an analytical function that implies a flat core. Results of these fits were presented in §4.2 using the Zhao (1996) profile with as well as the one suggested by S2006. The latter function gives a poor description of the overall measured density profiles, which tend to be significantly less curved than the analytical model (see Fig. 5). The obtained (overestimated) values for r (, see §4.2) are reported in column 13 of Table 3. When the S2006 function is fitted to only the inner 30 kpc, the fits improve and the estimated core radii become smaller by roughly a factor of two (column 14). However, even for this case the core radii seems to be upper limits, with the exception of halo D512.
We have found that the WDM profiles are better described by the Zhao function with . The best fits to the measured profiles give extrapolated core radii times smaller than the S2006 profile fitted to only the inner 30 kpc. Should we have the sufficient resolution to resolve the flat cores, their radii would lie in between r and . The only halo for which the flat core is patent at our resolution limit is D512; a visual inspection shows that the core radius is close to 1 kpc.
Our results show that the different estimates of the core radius increase by less than a factor of 1.5 from the simulations with (0)=0.1 km/s to the ones with (0)=0.3 km/s. The amount of this increase is less than the one we would predict using the monolithic collapse approximation of Avila-Reese et al. (2001); according to this approximation, the core radius of halos formed at the same time depends linearly on the injected at the maximum expansion of the perturbation, . We have estimated the predicted values of r for our halos by using this approximation. From the simulations, we find that the redshifts of maximum expansion of halos C, D, and E are roughly and 1.8, respectively; these redshifts are practically the same for the different values of . The calculated r for (0) = 0.1 (0.3) km/s are then 221 (663), 240 (720), and 209 (627) pc, respectively. Thus, in general these predictions give core radii just in between r and (see Table 1), though the dependence on is much more pronounced than for the (extrapolated) core radii estimated from the fits to our halos.
5.3. How is the structure of halos at the damping scale?
The simulations carried out in this paper allowed us also to explore the structure of dark matter halos formed from perturbations at the scale of damping of the linear power spectrum. The cutoff in the power spectrum used here corresponds to a relatively large mass, M. Therefore, the formation of the first structures in this model, namely those structures with masses close to , happens relatively late. We speculate that the formation process of the truncation halos is generic. If this is true, then the structure of the late–formed truncation halos simulated here with several millions of particles (=0) should be similar to the structure of early–formed truncation (micro)halos in models with much smaller filtering masses than the one used here, for example in the CDM models. If this is the case, then our results may enrich the discussion about the formation and structure of the first microhalos in CDM models (Earth–mass scales).
We have found a clear systematic trend in the density profiles of the simulated truncation halos: they are significantly steeper than in the inner regions, , and lie below the best NFW fits in the intermediate region. CDM halos with an inner slope steeper than have been reported in other contexts: recently merged group– and cluster–size halos (Knebe et al., 2002; Tasitsiomi et al., 2004) or microhalos formed at the scale of CDM power–spectrum damping (Diemand et al., 2005). It has been argued that the recent major merger is the one to blame for the steepening of the density profile while subsequent secondary infall modifies the external region. For the halos at the damping scale, rather than a major merger, the dynamical situation corresponds to a fast (cuasi–monolithic) collapse. However, in both cases, the process is dynamically violent. We have checked that the obtained density profiles do not correspond to a transient configuration. As mentioned in §5.1, for halo D the profile remains almost unchanged until .
Extrapolating our results to the damping scale of CDM, we can speculate that CDM microhalos may be significantly steeper than the NFW profile. This implies that the possible contribution of surviving microhalos to the ray flux originated by neutralino annihilation, might be comparable to the central flux from host halos (Diemand et al., 2006). Another implication of our results could be related to the buildup of the inner density profile of dark halos in general. Dehnen (2005) and Kazantzidis et al. (2006) argued that the assembly of halo inner density profiles happens very early in the history of the universe and specifically that the inner slope of the cuspier progenitor survives up to the final halo. If the density profiles in our simulated halos are representative of objects formed at the damping of power scales in general, then, according to the mentioned studies, the central slope of present day dark matter halos should be much steeper than the NFW profile.
In order to verify our results, a more systematic study halo structure at the scale of damping is required, exploring any possible dependence with the shape of the cutoff and the power spectrum slope at the scale of damping. Currently the only studies discussing a similar situation report different results for the scale of galaxy clusters(Moore et al., 1999b) and the microhalos (Diemand et al., 2005). It is unclear if this suggests a slope dependence on the profile of the smallest dark matter halos with the power spectrum slope.
We have studied by means of cosmological N–body simulations the structure of dark halos formed in the context of a WDM model corresponding to a non–thermal sterile neutrino particle of mass =0.5 keV. The first series of simulations did not include the injection of a random velocity component, , to the particles and were aimed at exploring the structure of the halos formed from perturbations at the damping scale in the linear power spectrum ( for the concrete WDM model studied here). The second series of simulations included a component of a (i) a =0.5 keV thermal neutrino ( km/s), and one that roughly corresponds (ii) to a =0.5 keV non–thermal sterile neutrino ( km/s). This latter model was used to generate the initial power spectrum. These simulations were aimed at exploring the effect that a random velocity component has on the inner structure of halos; in particular, at dilucidating whether constant–density cores are produced or not in the simulations. The results of our study lead us to the following two main conclusions:
The structure of halos formed from perturbations of scales close to , and resolved with up to more than 16 millions of particles (with =0), is peculiar: the inner density profile () is systematically steeper than the best corresponding NFW fit (and the respective CDM counterpart), and the overall density profile () tends to be less curved than the best NFW fit; the outer profile slope is never steeper than . According to our tests, these differences with respect to the structure of halos assembled hierarchically, can hardly be attributed to a peculiar dynamical state of the halos simulated here.
The effect of adding to the particles produces a significant flattening of the inner density profile ( corresponding to ) of the simulated halos. The different estimated (extrapolated) sizes of the nearly constant–density cores are of the order of the theoretical predictions, which give values below our resolution limit. For the halo masses simulated here, M, the flat core radii estimated from different fittings are between kpc. An increase in (0) from 0.1 to 0.3 km/s produces an increase in the extrapolated core radii of a factor 1.5 or less. For one of our simulations (halo D512), the presence of a nearly constant–density core, of radius kpc, is already revealed at the resolution limit; the same halo simulated with a different random velocity seed is less flattened.
Although the simulations presented here refer to a concrete WDM model, they can be interpreted within a wide range of contexts. The density profile of dark halos with masses close to the truncation mass in the linear power spectrum is systematically different from the NFW profile; in particular, the inner regions tend to be steeper. These could have important implications in the context of CDM models if a significant fraction of microhalos formed at the free–streaming CDM scales () have survived until the present epoch. In this case, the predicted ray flux from the neutralino annihilation in the center of these cuspy microhalos might be comparable to the central flux from host halos. On the other hand, the fact that the microhalos are so cuspy, could have some interesting implications in the building up of the next hierarchies of the halo assembling, as well as in the inner structure of the larger halos.
Regarding the effects of random velocity injection to the particles, our results show evidence of significant inner flattening of the halo density profile at our resolution radii. These resolutions are not enough to test directly the predicted core radii by phase–space constraints (e.g., Hogan & Dalcanton, 2000) or by dynamical models of gravitational collapse with initial random velocities (e.g., Avila-Reese et al., 2001; Bode et al., 2001). However, the inner extrapolations of the best–fit models to our simulated halos are consistent with such predictions.
- This radius is defined as the radius at which the mean overdensity is equal to the predicted by a top-hat spherical collapse which is 337 for our selected cosmological model.
- The random velocity component was linearly added to the velocities correspondingt to the Zeldovich approximation at the onset of the simulation.
- footnotetext: All the halos presented in this Table were resimulated from a 10Mpc box simulation.
- footnotetext: Note that the NFW function does not provide a good fit to the density profiles of halos with . However, for completeness, we report here the value of c obtained from the fit.
- Abazajian, K. 2006a, PhRvD, 73, 063506
- Abazajian, K. 2006b, PhRvD, 73, 063513
- Avila-Reese, V., Firmani, C., Klypin, A., & Kravtsov, A.V. 1999, MNRAS, 310, 527
- Avila-Reese, V. Colín, P., Valenzuela, O., D’Onghia, E., & Firmani, C. 2001, ApJ, 559, 516
- Benson, A.J., Frenk, C.S., Lacey, C.G., Baugh, C.M. & Cole, S. 2002, 333, 177
- Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93
- Bullock, J. S., Kravtsov, A.V., & Weinberg, D.H. 2000, ApJ, 539, 517
- Bullock, J. S., Kolatt, T. S., Sigad, Y., Somerville, R. S., Kravtsov, A. V., Klypin, A. A., Primack, J. R., & Dekel, A. 2001, MNRAS, 321, 559
- Busha, M. T., Evrard, A. E., & Adams, F. C. 2007, ApJ, 665, 1
- Cembranos, J. A., Feng, J. L., Rajaraman, A., & Takayama, F. 2005, Physical Review Letters, 95, 181301
- Ceverino, D., & Klypin, A. 2007, ArXiv Astrophysics e-prints, arXiv:astro-ph/0703544
- Colín, P., Avila-Reese, V., & Valenzuela, O. 2000, ApJ, 542, 622
- Colín, P., Avila-Reese, V., Valenzuela, O., & Firmani, C. 2002, ApJ, 581, 777
- Colín, P., Valenzuela, O., & Klypin, A. 2006, ApJ, 644, 687
- de Blok, W. J. G., McGaugh, S. S., Bosma, A., & Rubin, V. C. 2001, ApJ, 552, L23
- Dehnen, W. 2005, MNRAS, 360, 892
- Diemand, J., Moore, B., & Stadel, J. 2004, MNRAS, 353, 624
- Diemand, J., Moore, B., & Stadel, J. 2005, Nature, 433, 389
- Diemand, J., Zemp, M., Moore, B., Stadel, J., & Carollo, C.M. 2005, MNRAS, 364, 665
- Diemand, J., Kuhlen, M., & Madau, P. 2006, ApJ, 649, 1
- Diemand, J., Kuhlen, ??? & Madau, P. 2006, ApJaccepted (astro-ph/0611370)
- Dubinski, J., & Carlberg, R.G. 1991, ApJ, 378, 496
- Eke, V. R., Navarro, J. F., & Steinmetz, M. 2001, ApJ, 554, 114
- El-Zant, A. A., Hoffman, Y., Primack, J., Combes, F., & Shlosman, I. 2004, ApJ, 607, L75
- J. L. Feng, A. Rajaraman, & F. Takayama 2003, PhRvL. 91, 011302
- Gao, L., White, S. D. M., Jenkins, A., Frenk, C. S., & Springel, V. 2005, MNRAS, 363, 379
- Gentile, G., Burkert, A., Salucci, P., Klein, U., & Walter, F. 2005, ApJ, 634, L145
- Gentile, G., Tonini, C., & Salucci, P. 2007, A&A, 467, 925
- Gilmore, G., Wilkinson, M.I., Wyse, R., Kleyna, J.T., Koch, A., & Evans, N.W. 2007, ApJ, accepted (astro-ph/0703308)
- Goerdt, T., Moore, B., Read, J. I., Stadel, J., & Zemp, M. 2006, MNRAS, 368, 1073
- Governato, F., Willman, B., Mayer, L., Brooks, A., Stinson, G., Valenzuela, O., Wadsley, J., & Quinn, T. 2007, MNRAS, 374, 1479
- Götz, M., & Sommer-Larsen, J. 2003, Ap&SS, 284, 341
- Hayashi, E., Navarro, J. F., Taylor, J. E., Stadel, J., & Quinn, T. 2003, ApJ, 584, 541
- Hayashi, E., & Navarro, J. F. 2006, MNRAS, 373, 1117
- Hogan, C.J., & Dalcanton, J.J. 2000, PhRvD, 62, 063511
- Huss, A., Jain, B., & Steinmetz, M. 1999, ApJ, 517, 64
- Kaplinghat, M. 2005, PhRvD, 72, 063510
- Kamionkowski, M., & Liddle, A. R. 2000, PhRvL, 84, 4525
- Kauffmann , G., White, S. D. M., & Guideroni, B. 1993, MNRAS, 264, 201
- Kazantzidis, S., Zentner, A. R., & Kravtsov, A. V. 2006, ApJ, 641, 647
- Klypin, A., & Holtzman, J. 1997, preprint (astro-ph/9712217)
- Klypin, A.A., Kravtsov, A.V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82
- Klypin, A., Gottlöber, S., Kravtsov, A.V., & Khokhlov, A.M. 1999, ApJ, 516, 530
- Klypin, A., Kravtsov, A.V., Bullock, J.S., & Primack, J.R. 2001, ApJ, 554, 903
- Knebe, A., Devriendt, J. E. G., Mahmood, A., & Silk, J. 2002, MNRAS, 329, 813
- Knebe, A., Devriendt, J. E. G., Gibson, B. K., & Silk, J. 2003, MNRAS, 345, 1285
- Kravtsov, A.V., Klypin, A., & Khokhlov, A.M. 1997, ApJS, 111, 73
- Kravtsov, A.V., Berlind, A.A., Wechsler, R.H., Klypin, A.A., Gottlöber, S. Allgood, B., & Primack, J.R. 2004, ApJ, 609, 35
- Ma, C.-P., & Boylan-Kolchin, M. 2004, PhRvL, 93, 021301
- Moore, B. 1994, Nature, 370, 629
- Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J., Tozzi, P. 1999, ApJ, 524, L19
- Moore, B., Quinn, T., Governato, F., Stadel, J., & Lake, G.1999, MNRAS, 310, 1147
- Moore, B., Quinn, T., Governato, F., Stadel, J., & Lake, G. 1999, MNRAS, 310, 1147
- Moore, B., Calcáneo-Roldán,C., Stadel, J. Quinn, T., Lake, G., Ghigna, S., & Governato, F. 2001, Phys. Rev. D, 64, 063508
- Narayanan, V.K., Spergel, D.N., Davé, R., Ma, C.-P. 2000, ApJ, 543, 103
- Navarro, J.F., Frenk, C.S., & White, S.D.M. 1997, ApJ, 490, 493
- Navarro, J.F., Hayashi, E., Power, C., Jenkins, A.R., Frenk, C.S., White, S.D.M., Springel, V., Stadel, J., & Quinn, T.R. 2004, MNRAS, 349, 1039
- Oda, T., Totani, T., & Nagashima, M. 2005, ApJ, 633, L65
- Profumo, S., Sigurdson, K., & Kamionkowski, M. 2006, PhRvL, 97, 031301
- Rhee, G., Valenzuela, O., Klypin, A., Holtzman, J., & Moorthy, B. 2004, ApJ, 617, 1059
- Sánchez-Salcedo, F. J., Reyes-Iturbide, J., & Hernandez, X. 2006, MNRAS, 370, 1829
- Seljak, U., et al. 2005, Phys. Rev. D, 71, 103515
- Seljak, U., et al. 2006, Phys. Rev. D, 73, 063513
- Sellwood, J. A. 2007, ArXiv e-prints, 704, arXiv:0704.1047
- Simon, J.D., Bolatto, A.D., Leroy, A., Blitz, L., Gates, E.L. 2005, ApJ, 621, 757
- Simon, J. D., & Geha, M. 2007, ArXiv e-prints, 706, arXiv:0706.0516
- Spergel, D. N., & Steinhardt, P. J. 2000, PhRvL, 84, 3760
- Spergel, D.N., et al. 2006, ApJ, submitted (astro-ph/0603449)
- Springel, V., Frenk, C. S., & White, S. D. M. 2006, Nature, 440, 1137
- Stasielak, J., Biermann, P.L., & Kusenko, A. 2007, ApJ, 654, 290
- Steen, F. D. 2006, JCAP, 9, 1
- Steffen, F. D. 2006, Journal of Cosmology and Astro-Particle Physics, 9, 1
- Strigari, L.E., Bullock, J.S., Kaplinghat, M., Kravtsov, A.V., Gnedin, O.Y., Abazajian, K., & Klypin, A. A. 2006, ApJ, 6502, 306 (S2006)
- Strigari, L. E., Kaplinghat, M., & Bullock, J. S. 2007, Phys. Rev. D, 75, 061303
- Tasitsiomi, A., Kravtsov, A. V., Gottlöber, S., & Klypin, A. A. 2004, ApJ, 607, 125
- Tremaine, S. & Gunn, J.E. 1979, PhRvL, 42, 407
- Valenzuela, O., Rhee, G., Klypin, A., Governato, F., Stinson, G., Quinn, T., & Wadsley, J. 2007, ApJ,657, 773
- Valinia, A., Shapiro, P.R., Martel, H., Vishniac, E.T. 1997, ApJ, 479, 46
- Viel, M., Lesgourgues, J., Haehnelt, M.G., Matarrese, S., & Riotto, A. 2005, PhRvD, 71, 063534
- Viel, M., Lesgourgues, J., Haehnelt, M. G., Matarrese, S., & Riotto, A. 2006, PhRvL, 97, 071301
- Wang, J., & White, S. D. M. 2007, ArXiv Astrophysics e-prints, arXiv:astro-ph/0702575
- Weinberg, M. D., & Katz, N. 2007, MNRAS, 1499
- White S. D. M., 1996, in “Cosmology and Large Scale Structure Formation and Evolution of Galaxies”, Schaeffer R., Silk J., Spiro M., Zinn-Justin J., eds.
- Zhao, H. S. 1996, MNRAS, 278, 488