A&A 458, 31-37 (2006)
DOI: 10.1051/0004-6361:20065551

Semi-analytical predictions of density profiles for $\Lambda$CDM haloes

N. Hiotelis

1st Experimental Lyceum of Athens, Ipitou 15, Plaka, 10557 Athens, Greece

Received 5 May 2006 / Accepted 29 June 2006

Abstract
We combine the results of merger trees realizations, predicted by the extended Press-Schechter theory, to the assumption of "stable clustering'' in order to calculate density profiles of dark matter haloes. Our results show that: 1) Haloes of different masses have different concentrations. More specifically, concentration is a decreasing function of the mass. The relation between concentration and virial mass predicted by our results is in good agreement with the predictions of large cosmological N-body simulations. 2) The slope of the density profile of dark matter haloes is flatter at the inner regions and steepens outwards. At a given fraction of the virial radius the slope (defined as minus the derivative of the logarithm of the density with respect to the logarithm of the radial distance) is a decreasing function of the virial mass of the halo. At 0.01 of the virial radius $R_{\rm vir}$ the slope ranges from 0.4 to 1.5 for masses in the range $5\times 10^9~{h}^{-1}~M _{\odot}$ to $10^{14}~{h}^{-1}~M_{\odot}$. At distance $R_{\rm vir}$ the values of the slope are in the range 2.5 to 5.5. We note that the small values of the inner slope predicted in this paper are closer to the predictions of observations than the results of N-body simulations. 3) Comparing haloes of the same present day mass we found that: a) Haloes with large rate of recent mass increase, show flatter outer density profiles than those with small recent mass increase. b) The concentration becomes larger for increasing recent mass growth rate.

Key words: galaxies: formation - galaxies: halos - galaxies: structure - methods: numerical - methods: analytical - cosmology: dark matter

1 Introduction

It is likely that structures in the Universe grow from small initially Gaussian density perturbations that progressively detach from the general expansion, reach a maximum radius and then collapse to form bound objects. Larger haloes are formed hierarchically by mergers between smaller ones.

The above scenario of formation is usually studied by two different kinds of methods. The first kind is the N-body simulations that are able to follow the evolution of a large number of particles under the influence of the mutual gravity, from initial conditions to the present epoch. The second kind consists of semi-analytical methods. Among them, the Press-Schechter (PS) approach and its extensions (EPS) are of great interest since they allow us to compute mass functions (Press & Schechter 1974; Bond et al. 1991) to approximate merging histories (Lacey & Cole 1993, LC93 hereafter; Bower 1991; Sheth & Lemson 1999b) and to estimate the spatial clustering of dark matter haloes (Mo & White 1996; Catelan et al. 1998; Sheth & Lemson 1999a).

One of the interesting problems related to the formation of dark matter haloes is that of their density profiles. Numerical simulations show that a good model for the density profile  $\rho_{\rm M}$ of a halo is that given by the formula,

$$\rho_{\rm M}(r)=\frac{\rho_{\rm c}}{(r/r_{\rm s})^{\kappa}[1+(r/r_{\rm s})^{\lambda}]^{\mu}}$$ (1)

where $\rho_{\rm c},\kappa, \lambda, \mu$ and $r_{\rm s}$ are constants (in the sense that they do not depend on the radial distance rfrom the center of the system). The logarithmic slope $\gamma$, is then given by the relation

$$\gamma(r)\equiv-\frac{{\rm d}\ln[\rho_{\rm M}(r)]}{{\rm d}\ln... ... \mu \frac{(r/r_{\rm s})^{\lambda}}{1+(r/r_{\rm s})^{\lambda}}$$ (2)

Navarro et al. (1997) (NFW hereafter) proposed the values $\kappa=1,\lambda=1,\mu=2$ while Moore et al. (1998) proposed $\kappa=1.5, \lambda=1.5$ and $\mu=1$ and Jing & Suto (2000) $\kappa=1.5,\lambda=1$ and $\mu=1.5$. The predictions of N-body simulations agree that the slope at large distances is $\kappa+\lambda+\mu=3$. In that sense, Reed et al. (2005, R05 hereafter) proposed a model with $\lambda=1$ and $\mu=3-\kappa$, where $\kappa$ depends on the virial mass of the halo to fit the density profiles resulting from their N-body simulations.

We define the concentration c by the common way, $c\equiv R_{\rm vir}/{r_2}$, that is the ratio of the virial radius of the system $R_{\rm vir}$ to the radius r₂, where the above logarithmic slope becomes equal to 2. For the NFW profile, r₂ equals the scale radius $r_{\rm s}$ but this is not true for other values of $\kappa,\lambda$ and $\mu$.

The results of N-body simulations show that concentration is a decreasing function of the mass of the halo. This shows that in a hierarchical clustering scenario, larger haloes are formed later.

As regards the inner asymptotic slope, given by the value of $\kappa$, the results differ not only between them but also with the predictions of observations that report significantly smaller values of $\kappa$ (Sand et al. 2002, 2004).

In this paper, we present density profiles predicted by semi-analytical methods based on the EPS theory. We hope, in this way, to better understand the physics involved during the collapse and relaxation of dark matter haloes. We compare our results to empirical formulas resulted from the predictions of numerical simulations. We find a satisfactory agreement.

The paper is organized as follows: in Sect. 2, basic equations are summarized. In Sect. 3, we present our results while a discussion is given in Sect. 4.

2 Basic equations

According to the hierarchical scenarios of structure formation, a region collapses at time t if its overdensity at that time exceeds some threshold. The linear extrapolation of this threshold up to the present time is called a barrier, B. A likely form of this barrier is:

$$B(S,t)=\sqrt{aS_*}[1+\beta(S/aS_*)^{\gamma}].$$ (3)

In the above equation $\alpha$, $\beta$ and $\gamma$ are constants, $S_*\equiv S_*(t)\equiv \delta^2_c(t)$ where $\delta_{\rm c}(t)$ is the linear extrapolation up to the present day of the initial overdensity of a spherically symmetric region that collapsed at time t. Additionally, $S\equiv \sigma^2(M)$, where $\sigma^2(M)$ is the present day mass dispersion on comoving scale containing mass M. S depends on the assumed power spectrum.

The spherical collapse model has a barrier that does not depend on the mass (e.g. LC93). For this model, the values of the parameters are a=1 and $\beta=0$. The ellipsoidal collapse model (EC) (Sheth & Tormen 1999) has a barrier that depends on the mass (moving barrier). The values of the parameters are a=0.707, $\beta=0.485$, $\gamma=0.615$ and are adopted either from the dynamics of ellipsoidal collapse or from fits to the results of N-body simulations.

Sheth & Tormen (2002) showed that given a mass element - that is a part of a halo of mass M₀ at time t₀ - the probability that at earlier time t this mass element was a part of a smaller halo with mass M, is given by the equation:

$$f(S,t/S_0,t_0){\rm d}S= \frac{1}{\sqrt{2\pi}}\frac{\vert T(S,t/S_... ...ert}{(\Delta S)^{3/2}} \exp\left[-\frac{(\Delta B)^2}{2\Delta S}\right]{\rm d}S$$ (4)

where $\Delta S=S-S_0$ and $\Delta B=B(S,t)-B(S_0,t_0)$ with S=S(M), S₀=S(M₀).

The function T is given by:

$$T(S,t/S_0,t_0)\!=\!B(S,t)-B(S_0,t_0)\!+\! \sum_{n=1}^{5}\frac{[S_0-S]^n}{n!} \frac{{\rm\partial ^n}}{\partial S^n}B(S,t)$$ (5)

Eq. (4) can obviously predict the unconditional mass probability, f(S,t), which is simply the probability that a mass element is at time t a part of a halo of mass M, by setting S₀=0, and B(S₀,t₀)=0. It is easy to note that the quantity Sf(S,t) is a function of the variable $\nu$ alone, where $\nu\equiv \delta_{\rm c}(t)/\sigma(M)$. $\delta_{\rm c}$ and $\sigma$ evolve with time by the same way, thus the quantity Sf(S,t) is independent of time. Setting $2Sf(S,t)=\nu f(\nu)$ one obtains the so-called multiplicity function $f(\nu)$. The multiplicity function is the distribution of first crossings of a barrier $B(\nu)$ (that is why the shape of the barrier influences the form of the multiplicity function), by independent uncorrelated Brownian random walks (Bond et al. 1991). The multiplicity function is related to the comovimg number density of haloes of mass M at time t, N(M,t), by the relation,

$$\nu f(\nu)=\frac{M^2}{\rho_b(t)}N(M,t)\frac{{\rm d}\ln M}{{\rm d}\ln \nu}$$ (6)

that results from the excursion set approach (Bond et al. 1991). In the above equation, $\rho_b(t)$ is the density of the model of the Universe at time t.

Using a barrier of the form of Eq. (3) in the unconditional mass probability, one finds for $f(\nu)$ the expression:

$$f(\nu) \sqrt{2a / \pi}[1+\beta(a {\nu}^2)^{-\gamma}g(\gamma)]$$ =

$$\times\exp\left(-0.5a\nu^2[1+\beta(a\nu^2)^{-\gamma}]^2\right)$$ (7)

where

$$g(\gamma)= \left\vert 1-\gamma +\frac{\gamma (\gamma -1)}{2... ...(\gamma-1)\cdot \cdot \cdot (\gamma-4)}{5!} \right\vert\cdot$$ (8)

Recent comparisons show that the use of EC model improves the agreement between the results of EPS methods and those of N-body simulations. For example Yahagi et al. (2004) show that the multiplicity function resulting from N-body simulations is far from the predictions of spherical model while it shows an excellent agreement with the results of EC model. On the other hand Lin et al. (2003) compared the distribution of formation times of haloes formed in N-body simulations with the formation times of haloes formed in terms of the spherical collapse model of EPS theory. They found that N-body simulations give smaller formation times. Hiotelis & Del Popolo (2006) showed that using the EC model, formation times are shifted to smaller values than those predicted by a spherical collapse model. Thus EC model is much more realistic than the spherical one and therefore we use this model for our calculations.

We assume a number of haloes with the same present day mass M₀ - at present epoch t₀ - and we study their past using merger-trees. This is done by finding their progenitors -haloes that merged and formed the present day haloes- at previous times. The procedure for a single halo is as follows: A new time t<t₀ is chosen. Then a value $\Delta S$ is chosen from the desired distribution given by Eq. (4). The mass M_p of a progenitor is found by solving for M_p the equation $\Delta S=S(M_p)-S(M_0)$. If the mass left to be resolved M₀-M_p is large enough, the above procedure is repeated so a distribution of the progenitors of the halo is created at t. If the mass left to be resolved - that equals to M₀ minus the sum of the masses of its progenitors - is less than a threshold then we proceed to the next time by analyzing with the same procedure the mass of each progenitor. The most massive progenitor at t is considered as the mass of the initial halo at that time.

The above procedure is repeated for a large number of haloes of the same mass and at every time step, the average mass is found. Thus a smooth curve is derived that shows the growth of mass as a function of time (or the scale factor). A complete description of the above numerical method is given in Hiotelis & Del Popolo (2006). The algorithm - known as N-branch merger-tree- is based on the pioneered works of LC93, Somerville & Kollat (1999) and van de Bosch (2002).

In our calculations, we use a flat model for the Universe with present day density parameters $\Omega_{\rm m,0}=0.3$ and $\Omega_{\Lambda,0}\equiv \Lambda/3H_0^2=0.7$. $\Lambda$ is the cosmological constant and H₀ is the present day value of Hubble's constant. We used the value $H_0=100~{\rm hKMs^{-1}~Mpc^{-1}}$ with h=0.7 and a system of units with $m_{\rm unit}=10^{12}~M_{\odot}~h^{-1}$, $r_{\rm unit}=1~h^{-1}~{\rm Mpc}$ and a gravitational constant G=1. At this system of units $H_0/H_{\rm unit}=1.5276.$

As regards the power spectrum, we used the $\Lambda CDM$ form proposed by Smith et al. (1998). The power spectrum is smoothed using the top-hat window function and is normalized for $\sigma_8\equiv\sigma(R=8~h^{-1}~{\rm Mpc})=1$.

3 Results

We studied seven cases with different masses. In our system of units these masses are $5\times 10^{-3}, 1\times 10^{-2}, 0.1, 1., 10, 100$ and 1000 for the cases 1, 2, 3, 4, 5, 6 and 7, respectively. For every case, we produced a number of $N_{\rm res}=50~000$ realizations. This is the initial number of haloes of given mass. At every time t and for every case, we calculate the mean mass of all realizations and we call mass of a halo at time t, the mean mass of all the realizations of every case. This procedure gives smooth mass growth profiles that are in good agreement with the results of N-body simulations (Hiotelis & Del Popolo 2006).

We proceed by calculating density profiles using the assumption of "stable clustering'' (e.g. Nusser & Sheth 1999; Manrique et al. 2003; Hiotelis 2003). According to this assumption the haloes grow inside-out. The accreted mass is deposited at an outer spherical shell without changing the inner density profile. Although this is an ideal situation, it may not be far from reality. For example, Taylor & Babul (2005), among other results of their N-body simulations, state that "The density of the main system changes with time, although in practice this change is mainly confined to the outer parts of the halo''.

The radial extent of a halo is defined by its virial radius $R_{\rm vir}$, that is the radius containing a mass with mean density $\Delta_{\rm vir}$ times the current mean density of the Universe, $\rho_{b}$. $\Delta_{\rm vir}$ is, in general, a function of time that depends on the assumed cosmology (Bryan & Norman 1998). Consequently, the mass $M_{\rm vir}$ contained inside  $R_{\rm vir}$ at scale factor a, satisfies the following equations:

$$M_{\rm vir}(a)=4\pi\int_{0}^{R_{\rm vir}(a)}r^2\rho(r){\rm d}r$$ (9)

$$M_{\rm vir}(a)=\frac{4}{3}\pi\Delta_{\rm vir}(a)\rho_b(a)R^3_{\rm vir}(a).$$ (10)

Differentiating Eq. (9) with respect to a one gets:

$$\rho(a)=\frac{1}{4\pi R^2_{\rm vir}(a)}\frac{{\rm d}M_{\rm vi... ...\left(\frac{{\rm d}R_{\rm vir}(a)}{{\rm d}a}\right)^{-1}\cdot$$ (11)

We define $X(a)\equiv \Delta_{\rm vir}(a)\Omega_{\rm m}(a)H^2(a)$ where $\Omega _{\rm m}(a)$ and H(a) are the density parameter and the Hubble's constant at scale factor a. X(a) depends on the assumed cosmology alone. Differentiating Eq. (10) with respect to a and substituting in Eq. (11) one has,

$$\rho(a)=\frac{3X(a)}{8\pi G}\left[1-\frac{{\rm d}\ln X(a)/{\... ...}\ln a}{{\rm d}\ln M_{\rm vir}(a)/{\rm d}\ln a} \right]^{-1}$$ (12)

which with the equation

$$R_{\rm vir}(a)=\left[\frac{2GM_{\rm vir}(a)}{X(a)}\right]^{1/3}$$ (13)

that results from Eq. (10), gives the density profile of the halo.

Let us assume that $M_{\rm vir}(a)\propto a^{g}$ while $X(a)\propto a^f$. It is interesting to study approximate predictions of the above two equations. First, X(a) depends on the assumed cosmology alone. Since $\Omega_{\rm m}(a)H^2(a)=\Omega_{\rm m0}H^2_0/a^3$, where subscripts 0 stand for the present day values and $\Delta_{\rm vir}$ is a slowly varying function of a, we see that X(a) scales approximately as a^-3 so $f\approx -3$. On the other hand, g varies significantly during the evolution. It is well known that the mass of a forming object grows fast at the early stages of its evolution, while the growth rate drops at later times when its mass grows very slow. Using (12) and (13) one can show that the slope $\gamma$ follows $\gamma \approx 3f/(f-g)$. Obviously, f<0 and g>0 so $\gamma <3$. Hiotelis & Del Popolo (2006) have shown that the values of ${\rm d}\ln M(a)/{\rm d}\ln a$ for values of a close to unity (present time) are in the range 0.4/a - 0.8/a, from small to large haloes. This gives a rough estimation of the slope at the virial radius in the range 2.64 to 2.36, which is in satisfactory agreement with the models based on the results of N-body simulations as those described in the introduction. For $\Delta_{\rm vir}$, we use the model $\Delta_{\rm vir}(a)=178\Omega^{0.45}_{\rm m}(a)$ resulting in $\Delta_{\rm vir}\approx 100$ at present epoch (a=1) (see Solanes et al. 2005).

From the above approximate expression we see that $\gamma$ is a decreasing function of g. Since g is a decreasing function of the scale factor a (the mass growth rate decreases with time), $\gamma$ increases during the build-up of the halo. Thus, the inner region is characterized by smaller values of $\gamma$ than the outer one. Hiotelis & Del Popolo (2006) found that for haloes of the same present day mass, different analytical models predict different mass growth rates. Namely, at the early stages of the evolution, EC model predicts haloes that grow with a larger mass growth rate than those haloes formed by the spherical model. Thus formation times are smaller in the EC model. This result is in the correct direction since it improves the agreement between the results of analytical and numerical methods (Lin et al. 2003). This improvement is also a motive for the study of density profiles that are predicted by EC model. Some of the differences to be expected are clear: The inner density profile should be flatter in the EC model than the spherical model since at their early stages, haloes evolve with larger g for the first model. However, since the spherical model predicts density profiles in good agreement with the NFW profile (Manrique et al. 2003), we expect inner density profiles with logarithmic slope $\gamma$ smaller than 1.

	Figure 1: Density profiles for four of the cases studied. From top to bottom dotted lines correspond to the results of the cases 1, 4, 6 and 7, respectively. Solid lines are fits by the formula of Eq. (1).
Open with DEXTER

	Figure 2: Triangles show the concentrations predicted by the model studied while the solid line shows the predictions of the toy-model of Bullock et al.
Open with DEXTER

In Fig. 1, we show the density profiles of four of the seven cases studied. From top to bottom, dotted lines show our predictions for cases 1, 4, 6 and 7, respectively, while solid lines show the fits of the above dotted curves by models of the form of Eq. (1) where $\rho_{\rm c},r_{\rm s},\kappa, \lambda$ and $\mu$ are fitting parameters. The values of the parameters are found by minimizing the sum $\sum[\rho_{\rm M}(r)-\rho(r)]^2$, where $\rho(r)$ are the values of the density predicted by our results. The minimization is performed using the unconstrained subroutine ZXMWD of IMSL mathematical library. The fits are shown to be very good. We also tried minimizations for $\kappa=1,\lambda=1$ and $\mu=2$ and $\rho_{\rm c},r_{\rm s}$ free parameters corresponding to the NFW model with $\lambda=1,\mu=3-\kappa$ and $\kappa,\rho_{\rm c},r_{\rm s}$ free parameters, corresponding to the R05 model but the respective fits were found less satisfactory. The resulting values of the parameters are used in order to calculate the concentrations of the haloes. Equation (2) yields that the logarithmic slope $\gamma$ equals n at distance r_n given by

$$r_n=\left[\frac{n-k}{\lambda\mu+k-n}\right]^{1/\lambda}r_{\rm s}.$$ (14)

Applying Eq. (14) for n=2, we calculate the concentration c by $c\equiv R_{\rm vir}/{r_2}$. The resulting concentrations are presented as triangles in Fig. 2. The solid line in the same figure shows the predictions of the toy-model of Bullock et al. (2001), constructed so as to give the values of the concentrations that result from the haloes formed in their large cosmological N-body simulations.

	Figure 3: Solid lines show the logarithmic slopes of the density profiles given by Eq. (2) for all seven cases studied. From top to bottom, these lines correspond to cases 1, 2, 3, 4, 5, 6 and 7, respectively. Dashed lines show, with the same order, the predictions of the model of Reed et al. (2005), based on the results of their N-body simulations.
Open with DEXTER

The value of the concentration depends on the mass of the halo, the time and the underlying cosmology. The procedure proposed by Bullock et al. (2001) is as follows: we are interested in the concentration of a halo that has a mass $M_{\rm vir}$ at scale factor a. We first solve for $a_{\rm c}$ the equation

$$\sigma\left[M_*(a_{\rm c})\right]=\sigma\left[FM_{\rm vir}\right]$$ (15)

where F=0.01 and M_* is the typical collapsing mass which at scale factor, a is given by $\sigma[M_*(a)]=1.686D(1)/D(a)$. Then the concentration is found by the relation: $c\equiv K \frac{a}{a_{\rm c}}$ where K is a constant (typically K=4) and D(a) is the value of the growth factor of overdensities given by the linear theory:

$$D(a)=\frac{1}{H^2_0}\frac{Y^{1/2}(a)}{a}\int_0^aY^{-3/2}(u){\rm d}u$$ (16)

where

$$Y(a)=1+\Omega_{\rm m,0}(a^{-1}-1)+\Omega_{\Lambda ,0}(a^2-1)$$ (17)

(Peebles 1980).

We note that for the cosmology used in this paper, the present-day value of M_* is $1.6\times 10^{13}~M_\odot~ h^{-1}$. It is obvious from Fig. 2, that the agreement is very satisfactory.

In Fig. 3, we compare the logarithmic slopes resulting by our calculations to those proposed by R05 to fit the results of N-body simulations. Solid lines are the predictions of the logarithmic slopes predicted by our modes while dashed lines are the predictions of the procedure proposed by the above authors. From top to bottom, the slopes correspond to the cases 1, 2, 3, 4, 5, 6 and 7, respectively. The procedure proposed by R05 is as follows:

For a halo of mass M, we find the exponent $\kappa$ by the formula $\kappa= 1.4-0.08 \log_{10}(M/M_*)$. Then a concentration c_k, defined by $c_{\kappa}=R_{\rm vir}/r_{\rm s}$, is calculated by the empirical formula $c_{k}\simeq 8(M/M_*)^{-0.15}$, that is a satisfactory fit to the results of N-body simulations, and then the formula $\gamma_{\rm M}(\widetilde{r})=\kappa+ \frac{(3-\kappa)c_{\kappa}\widetilde{r}}{1+c_{\kappa}\widetilde{r}}$ is applied. This formula results from Eq. (2) for $\lambda=1$, and $\mu=3-\kappa$, while $\widetilde{r}$ is the radial distance in units of the virial radius. We also note that R05 report a significant scatter in the values of  $c_{\kappa}$ in the results of their N-body simulations, which obviously leads to significant scatter to the values of the logarithmic slope.

There are some interesting features in this figure. First, the density profiles of heavier haloes are flatter. This order is obvious in both our results and those of N-body simulations. Second, the values of inner slope (at $R=0.01R_{\rm vir}$) resulting from our models are significantly smaller than those resulting from N-body simulations while the inverse occurs at the outer slopes (at $R=R_{\rm vir}$).

In Fig. 4, dotted lines show the central density profiles of the cases 3, 4, 5, 6 and 7 while solid ones show the fits by the formula of Eq. (1). This figure is drawn to show the high quality of the fit and to ensure that the logarithmic slope can be estimated with a good accuracy. Figure 5 presents the central slopes of the five more massive haloes that are represented by the cases 3, 4, 5, 6 and 7. It is clear that the asymptotic slope at the center of the system varies in the range of 0.35 to 0.4.

	Figure 4: Fits of the density profiles of inner regions. From top to bottom dotted lines correspond to the cases 3, 4, 5, 6, 7 of our model while solid lines are the fits by the formula of Eq. (1). The quality of the fit is sufficient for a reasonable estimation of the logarithmic slope.
Open with DEXTER

	Figure 5: The logarithmic slopes of the density profiles for the cases of Fig. 4. From top to bottom the lines correspond to cases 3, 4, 5, 6, 7. It is clear that larger haloes have at the same radial distance (measured in units of the virial radius) flatter inner density profiles. Additionally, it is clear that $\gamma$tends asymptotically to values smaller than 0.5.
Open with DEXTER

	Figure 6: The same as in Fig. 5 but for the outer region. From top to bottom, lines correspond to the cases 1, 2, 3, 4, 5, 6 and 7, respectively. It is clear that the slope $\gamma$ at $R=R_{\rm vir}$is a decreasing function of the virial mass of the halo.
Open with DEXTER

	Figure 7: The dependence of the concentration on the mass growth rate of the halo. Triangles represent the concentrations of haloes with large recent mass growth, squares represent the concentrations of all haloes of every case and circles the concentrations of haloes with small recent mass growth.
Open with DEXTER

Figure 8: The dependence on the outer slope on the mass growth rate. Solid lines correspond to case 3, dashed lines to case 4 and dotted lines to case 7. Upper lines correspond to haloes with small recent mass growth while lower lines correspond to haloes with large recent mass growth. It is obvious that haloes which recently increased significantly their masses have flatter outer density profiles.

Open with DEXTER

We have to note here that there are difficulties concerning the calculation of the slope of the density profile at inner and outer regions of a halo. If the derivation of slopes is carried out by finding the logarithmic derivative of the analytic formula designed to fit the density profile, it is clear that this could lead to large inaccuracies. For example, an analytic formula could fit well a large part of the curve but not its inner and/or outer region. Such an example can be seen in Fig. 8 of Tasitsiomi et al. (2004), where analytic formulas are presented along with the results of N-body simulations. It is obvious that, in some cases presented in that figure, analytic formulas are bad estimators of inner and outer slopes. Regarding our fits we note that in cases where the formula given by Eq. (1) is not able to fit the whole curve to a desired high accuracy, we used it to fit the inner or the outer regions separately. Thus, the resulting inner and outer slopes as well as the values of concentration are calculated with a high accuracy.

We additionally examined the role of the mass growth rate to the resulting concentrations and density profiles So, we analyze our samples of haloes by the following procedure: First, for every case we find the value of the scale factor $a_{\rm eq}$ at which about half of the haloes have mass smaller than one half of the present day mass. Thus, we divide the haloes of every case into two groups. The first group, called group S (small), contains haloes whose mass satisfies $M(a_{\rm eq})< M_0/2$ while the second group, called group L (large), contains haloes whose mass satisfies $M(a_{\rm eq})\geq M_0/2$. It is clear that haloes of group S have to show a larger recent mass growth rate than those of group L since the haloes of both groups have the same present day mass. Then, we found the concentration of every group. The results are presented in Fig. 7. Triangles are the concentrations of group S for every case, squares are those of all haloes and circles show the concentrations of group L. It is obvious that, although the differences are not large, haloes of group S show systematically larger concentrations than those of group L.

In Fig. 8, we plot the outer slopes for every group of the cases 3, 4 and 7. The upper solid line at the right of the figure corresponds to the group L of case 3, while the lower solid line corresponds to the group S of the same case. The upper and lower dashed lines correspond to the groups L and S of the case 4, respectively, while the upper and lower dotted lines correspond to the groups L and S of the case 7. The role of the mass growth is obvious in this figure Haloes which recently increased significantly their masses (they belong to groups S) have flatter outer density profiles than haloes of the same mass but which show no recent significant growth of their masses (group L). We note that the results of Fig. 8 are fully consistent with the approximate predictions of the method that are analyzed above. A significant recent increase of mass corresponds to large values of the exponent g and consequently to small values of $\gamma$.

No systematic dependence of the inner slope on the mass growth rate is found.

It has been reported (Ascasibar et al. 2003; Tasitsiomi et al. 2004) that haloes, which experienced a recent merger event, have lower concentrations and steeper inner profiles than more relaxed systems. We note that in N-body simulations a recent merger event probably leads to a redistribution of particles even inside the virial radius while in our model the halo remains always in virial equilibrium. However, a recent merging event does not necessarily correspond to a recent and significant mass growth rate. Thus the difference that appears (a recent merger event leads to lower concentration while a recent significant mass growth rate to a larger concentration) is due to the obvious inability of the "stable clustering'' approach to deal with all kinds of merger events. Virial equilibrium assumed for our approach holds under certain conditions that are related to the masses and the amounts of the intrinsic and orbital angular momentum of merging haloes. Further studies by both numerical and analytical methods are needed to improve our understanding about the physics of merger events. However, the limits of validity of "stable clustering'' approach are not clear yet. Manrique et al. (2003) and Hiotelis (2003) proposed that the "stable clustering'' approach is valid if the fractional increase of halo's mass by a merger is below a given threshold $\Delta_{\rm m}$ (with $\Delta_{\rm m}$ in the range 0.2-0.5). Although for the above values models seem to work well, it is obvious that kinematic parameters (as the amounts of angular momentum) have also to be taken into account.

4 Discussion

It is well known that observations and simulations give different values of the central logarithmic slope $\gamma$, of the density profile of dark haloes. Most of high resolution numerical simulations give central density profiles $\rho(r)\propto r^{-\gamma}$ with $\gamma$ as steep as 1.5 (Moore et al. 1998; Ghigna et al. 2000; Klypin et al. 2001) or as steep as 1 (NFW, Power et al. 2003; Stoehr et al. 2003). On the other hand, observations differ with those of N-body simulations. HI observations of dwarf galaxies (Flores & Primack 1994; Moore 1994) show that the central density is nearly constant. Additionally, the results from observations of larger systems, such as clusters of galaxies show the same difference too. For example, Sand et al. (2002, 2004) estimate a value for $\gamma \approx 0.35$ at the central region of the cluster MS22137-23. Kelson et al. (2002), found similar results for the cluster Abell 2199. Although both observational and numerical methods have been significantly improved during recent years the disagreement remains and it worsens if one takes into account that real systems show an almost flat central density, despite that their dark matter has undergone an adiabatic compression due to the infalling baryons (Blumental et al. 1986).

It is clear that the formation history of a halo is a very complicated process. Its mass growth history consists of periods of gentle increase, corresponding to the accretion of small amounts of the surrounding matter, and periods of sharp increase that correspond to mergers with other haloes. It is obvious that our assumption of stable clustering can hardly hold in cases of mergers between haloes of similar masses. On the other hand, the prediction of smooth density profiles requires a sufficient smooth approximation of mass growth. Smooth mass growth curves are derived by using a large number of haloes of the same mass for every case. We show that inner and outer slopes depend on the present day mass of the halo. Heavier haloes have smaller values of $\gamma$ in both inner and outer regions. We also show that heavier haloes have smaller concentrations than those of less heavy ones. The predicted values of concentration are very close to those predicted by the results of N-body simulations.

Inner and outer slopes predicted by the model studied here are different from those predicted from N-body simulations. Our results indicate that inner slopes depend mainly on the mass of the halo (or equivalently on the power spectrum)and not on the rate of its mass growth. The values of inner slope are close to those derived from observations. Outer slopes do not depend on the present day mass of the halo only but on the rate of its mass growth too. A low recent mass growth rate results in steep outer density profiles. In any case, these values seem to be different from those reported by numerical simulations. Only very massive haloes seem to have slopes at $R=R_{\rm vir}$ smaller than 3. Less massive haloes show outer slopes even larger than 5.

Since the predictions presented in this paper are so close to the results of observations for the central regions, we believe that the simple physical model used is very promising and further improvements could help us to understand better the physics involved in the formation of central regions of dark matter haloes.

Acknowledgements

We are grateful to the anonymous referee for helpful and constructive comments and discussions. We also acknowledge the Empirikion Foundation for its financial support, Dr. M. Vlachogiannis and K. Konte for assistance in manuscript preparation.

References

• Ascasibar, Y., Yepes, G., M $\ddot{{\rm u}}$ller, V., & Gottl $\ddot{{\rm o}}$ber, S. 2003, MNRAS, 346, 731 [NASA ADS] (In the text)
• Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27 [NASA ADS] [CrossRef]
• Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440 [NASA ADS] [CrossRef] (In the text)
• Bower, R. 1991, MNRAS, 248, 332 [NASA ADS] (In the text)
• Bryan, G., & Norman, M. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] (In the text)
• Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [NASA ADS] (In the text)
• Catelan, P., Lucchin, F., Matarrese, S., & Porciani, C. 1998, MNRAS, 297, 692 [NASA ADS] (In the text)
• Flores, R. A., & Primack, J. R. 1994, ApJ, 427, L1 [NASA ADS] [CrossRef] (In the text)
• Ghigna, S., Moore, B., Governato, F., et al. 2000, ApJ, 544, 616 [NASA ADS] [CrossRef] (In the text)
• Hiotelis, N. 2003, MNRAS, 344, 149 [NASA ADS] (In the text)
• Hiotelis, N., & Del Popolo, A. 2006, Ap&SS, 301, 167 [NASA ADS] (In the text)
• Jing, Y. P., & Suto, Y. 2000, ApJ, 529, L69 [NASA ADS] [CrossRef] (In the text)
• Klypin, A. A., Kratsov, A. V., Bullock, J. S., & Primack, J. R. 2001, ApJ, 554, 903 [NASA ADS] [CrossRef] (In the text)
• Kelson, D. D., Zabludoff, A. I., Williams, K. A., et al. 2002, ApJ, 576, 720 [NASA ADS] [CrossRef] (In the text)
• Lacey, C., & Cole, S. 1993, MNRAS, 262, 627 [NASA ADS] (LC93) (In the text)
• Lin, W. P., Jing, Y. P., & Lin, L. MNRAS, 344, 1327
• Manrique, A., Raig, A., Salvador-Sol ${\rm\acute{e}}$, E., & Sanchis, T. 2003, ApJ, 593, 26 [NASA ADS] [CrossRef] (In the text)
• Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347 [NASA ADS] (In the text)
• Moore, B. 1994, Nature, 370, 629 [NASA ADS] [CrossRef] (In the text)
• Moore, B., Governato, F., Quinn, T., Stadel, J., & Lake, G. 1998, ApJ, 499, L5 [NASA ADS] [CrossRef] (In the text)
• Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] (NFW) (In the text)
• Nusser, A., & Sheth, R. K. 1999, MNRAS, 303, 685 [NASA ADS] (In the text)
• Peebles, P. J. E. 1980, The Large-Scale Structure of the Universe (Princeton, NJ: Princeton Univ. Press) (In the text)
• Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14 [NASA ADS] (In the text)
• Press, W., & Schechter, P. 1974, ApJ, 187, 425 [NASA ADS] [CrossRef] (In the text)
• Reed, D., Governato, F., Verde, L., et al. 2005, MNRAS, 357, 82 [NASA ADS] (R05) (In the text)
• Sand, D. J., Treu, T., & Ellis, R. S. 2002, ApJ, 574, L129 [NASA ADS] [CrossRef] (In the text)
• Sand, D. J., Treu, T., Smith, G. P., & Ellis, R. S. 2004, ApJ, 604, 88 [NASA ADS] [CrossRef] (In the text)
• Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [NASA ADS] (In the text)
• Sheth, R. K., & Lemson, G. 1999a, MNRAS, 304, 767 [NASA ADS]
• Sheth, R. K., & Lemson, G. 1999b, MNRAS, 305, 946 [NASA ADS]
• Sheth, R. K., & Tormen, G. 2002, MNRAS, 329, 61 [NASA ADS] (In the text)
• Smith, C. C., Klypin, A., Gross, M. A. K., Primack, J. R., & Holtzman, J. 1998, MNRAS, 297, 910 [NASA ADS] (In the text)
• Somerville, R. S., & Kollat, T. S. 1999, MNRAS, 305, 1 [NASA ADS] (In the text)
• Solanes, J. M., Manrique, A., Gonzalez-Casado, G., & Salvador-Sole, E. 2005, ApJ, 628, 45 [NASA ADS] [CrossRef] (In the text)
• Stoehr, F., White, S. D. M., Springel, V., Tormen, G., & Yoshida, N. 2003, MNRAS, 345, 1313 [NASA ADS] (In the text)
• Tasitsiomi, A., Kratsov, A. V., Gottlober, S., & Klypin, A. A. 2004, ApJ, 607, 125 [NASA ADS] [CrossRef] (In the text)
• Taylor, J. E., & Babul, A. 2005, MNRAS, 364, 515 [NASA ADS] (In the text)
• Yahagi, H., Nagashima, M., & Yoshii, Y. 2004, ApJ, 605, 709 [NASA ADS] [CrossRef] (In the text)
• van den Bosch, F. G. 2002, MNRAS, 331, 98 [NASA ADS]