A&A 382, 84-91 (2002)
DOI: 10.1051/0004-6361:20011620

Density profiles in a spherical infall model
with non-radial motions

N. Hiotelis[*]

Lyceum of Naxos, Chora of Naxos, Naxos 84300, Greece

Received 3 September 2001 / Accepted 5 November 2001

Abstract
A generalized version of the Spherical Infall Model (SIM) is used to study the effect of angular momentum on the final density profile of a spherical structure. The numerical method presented is able to handle a variety of initial density profiles (scale or not scale free) and no assumption of self-similar evolution is required. The realistic initial overdensity profiles used are derived by a CDM power spectrum. We show that the amount of angular momentum and the initial overdensity profile affect the slope of the final density profile at the inner regions. Thus, a larger amount of angular momentum or shallower initial overdensity profiles lead to shallower final density profiles at the inner regions. On the other hand, the slope at the outer regions is not affected by the amount of angular momentum and has an almost constant value equal to that predicted in the radial collapse case.

Key words: galaxies: formation - galaxies: halos - galaxies: structure - methods: numerical


1 Introduction

It is likely that dark matter halos are formed by the evolution of small density perturbations in the early Universe. The matter contained in a perturbed region progressively detaches from the general flow and after reaching a radius of maximum expansion it collapses to form an individual structure. The most simple case is when this region is spherical, isolated and undergoing a radial collapse. This is the spherical infall model (hereafter SIM). SIM has been extensively discussed in the literature (Gunn & Gott 1972; Gott 1975; Gunn 1977; Fillmore & Goldreich 1984; Bertschinger 1985; Hoffman & Shaham 1985, hereafter HS; White & Zaritsky 1992).

The final density profile after a collisionless evolution of the matter depends on its initial density profile and the underlying cosmology. Self-similarity solutions (Fillmore & Goldreich 1984; Bertschinger 1985, HS) show that a power-law initial density profile relaxes to a final density profile given by $\rho(r)\propto r^{-\alpha}$ with $\alpha\geq 2$. Furthermore, recent numerical studies that relax the assumption of self-similarity, also give final density profiles steeper than r-2. Lokas & Hoffman (2000b) found values of $\alpha$ in the range 2 to 2.3.

The density profiles of galactic halos do not seem to follow power laws. Numerical studies (Quinn et al. 1986; Frenk et al. 1988; Dubinski & Galberg 1991; Crone et al. 1994; Navarro et al. 1997; Cole & Lacey 1996; Huss et al. 1999; Fukushinge & Makino 1997; Moore et al. 1998; Jing & Suto 2000) showed that the profile of relaxed halos steepens monotonically with radius. The logarithmic slope $\alpha=-\frac{{\rm d}\ln\rho}{{\rm d}\ln r}$ is less than 2 near the center and larger than 2 near the virial radius of the system. The value of $\alpha$ near the center of the halo is not yet known. Navarro et al. (1997) claimed $\alpha=1$while Kravtsov et al. (1998) initially claimed $\alpha\sim 0.7$ but in their revised conclusions (Klypin et al. 2000) they argue that the inner slope varies from 1 to 1.5. Moore et al. (1998) found a slope $\alpha=1.5$ at the inner regions of their N-body systems.
In this paper we study the final density profiles predicted by the SIM, when non-radial motions are included. Studies concerning the role of non-radial motions have been presented by Ryden & Gunn (1987), Ryden (1988), Gurevich & Zybin (1988a, 1988b), Avila-Reese et al. (1998), White & Zaritsky (1992), Sikivie et al. (1997) and recently by Nusser (2001).
In Sect. 2 we discuss the SIM and the associated problems. In Sect. 3 the description of the numerical method as well as the way the angular momentum is included is given. The initial conditions and the results are given in Sect. 4 and are summarized in Sect. 5.

2 The spherical infall model (SIM)

SIM is based on the physical process described in Gunn (1977) and in Zaroubi & Hoffman (1993): in an expanding spherical region the maximum expansion radius $\zeta$(apapsis) of a shell is a monotonic increasing function of its initial radius $x_{\rm i}$ and is given by the relation:


 \begin{displaymath}%
\zeta=g(x_{\rm i})=\frac{1+\Delta_{\rm i}(x_{\rm i})}{1-\Omega^{-1}
_{\rm i}+\Delta_{\rm i}(x_{\rm i})}x_{\rm i},
\\
\end{displaymath} (1)

where $\Omega_{\rm i}$ is the initial value of the density parameter of the Universe and $\Delta_{\rm i}$ is the relative excess of mass inside the sphere of radius $x_{\rm i}$, given by


 \begin{displaymath}%
\Delta_{\rm i}(x_{\rm i})=\frac{M(x_{\rm i})-M_{\rm b}(x_{\...
..._{{\rm i}}}\int_0^{x_{\rm i}}x^2\delta_{\rm i}(x){\rm d}x.
\\
\end{displaymath} (2)

In (2), M is the mass of the spherical region, $M_ {\rm b}$ is the mass of the unperturbed Universe and $\delta_{\rm i}$ is the spherically symmetric perturbation of the density field ( $\delta_{\rm i}(x)=\frac{\rho(x)-\rho_{{\rm b},{\rm i}}}{\rho_{{\rm b},{\rm i}}}$where $\rho$ is the density and $\rho_{{\rm b},{\rm i}}$ is the constant density of the homogeneous Universe at the initial conditions). The time spent for a shell to reach its above turnaround radius is:

 \begin{displaymath}%
t_{{\rm ta}}=\frac{1+\Delta_{\rm i}(x_{\rm i})}{2H_{\rm i}\...
...{-1}_{\rm i}+\Delta_{\rm i}(x_{\rm i})]^{\frac{3}{2}}}\pi,
\\
\end{displaymath} (3)

where $H_{\rm i}$ is the value of Hubble's constant at the epoch of the initial conditions. The above equations are valid for bound shells, so the condition $1-\Omega^{-1}_{\rm i}+\Delta_{\rm i}(x_{\rm i})>0$ is satisfied. A shell, after reaching its turnaround radius, collapses, re-expands to a new (smaller) turnaround radius and so on. The limiting value of this radius, after a large number of such oscillations, is the final radius of the shell, corresponding to the relaxed state of the system. In any reasonable potential a shell spends most of its orbital time near the maximum radius (Gunn 1977). So, if shell crossing -occurring during the collapse stage- had no dynamical consequences, the final distribution of mass could be approximated by the distribution resulting if every shell stopped at its maximum expansion radius. This should lead to a "turnaround density profile'' $\rho_{\rm ta}$of the form


 \begin{displaymath}%
\rho_{\rm ta}(\zeta)=\left(\frac{x_{\rm i}}{\zeta}\right)^2...
...)
\left(\frac{{\rm d}\zeta}{{\rm d}x_{\rm i}}\right)^{-1}.
\\
\end{displaymath} (4)

In deriving (4) the conservation of mass is used ( $M(x_{\rm i})=M(\zeta)$). This is an important relation since this distribution of mass is used as the initial one in SIM. SIM assumes that the collapse is gentle enough. This means that the orbital period of the inner shell is much smaller than the collapse time of the outer shells (Zaroubi & Hoffman 1993). This implies that the radial action $\oint v(r){\rm d}r$, where v is the radial velocity, is an adiabatic invariant of the inner shell. As the outer shells collapse, the potential changes slowly and because of the above adiabatic invariant, the inner shell shrinks. The collapse factor depends on the time the mass of the outer shells (passing momentarily) spends inside the inner shell. Consider a shell with apapsis $\zeta$ and initial radius $x_{\rm i}$. The mass inside radius $\zeta$ is a sum of two components. The first one, (permanent component, $M_{\rm p}$), is due to the shells with apapsis smaller than $\zeta$ and the second (additional mass, $M_{\rm add}$) is the contribution of the outer shells passing momentarily through the shell $\zeta$. Because of the mass conservation, the permanent component is given by the following relation

 \begin{displaymath}%
M_{\rm p}(\zeta)=M(x_{\rm i})=\frac{4}{3}\pi\rho_{{\rm b},{\rm i}}x^3_{{\rm i}}[1+\Delta_{\rm i}(x_{\rm i})].
\\
\end{displaymath} (5)

The additional component is:

 \begin{displaymath}%
M_{\rm add}(\zeta)=\int_{\zeta}^R P_{\zeta}(x)
\frac{{\rm d}M(x)}{{\rm d}x}{\rm d}x.
\\
\end{displaymath} (6)

In (6), R stands for the radius of the system (the apapsis of the outer shell) and the distribution of mass M(x) is given by (4). $P_{\zeta}(x)$ is the probability of finding the shell with apapsis x inside radius $\zeta$, calculated as the ratio of the time the outer shell (with apapsis x) spends inside radius $\zeta$ to its period. In the general case of non-radial collapse this ratio is given by the relation $P_{\zeta}(x)=I(\zeta)/I(x)$ where


\begin{displaymath}%
I(r)=\int_{x_{\rm p}}^{\rm r}\frac{{\rm d}n}{v_x(n)},
\end{displaymath} (7)

where $x_{\rm p}$ is the pericenter of the shell with apapsis x and vx(n) is the radial velocity of the shell with apapsis x as it passes from radius n. If the collapse is radial then $x_{\rm p}=0$. After the calculation of $M_{\rm add}$ the collapse factor $f(x_{\rm i})$ of a shell with initial radius $x_{\rm i}$ and apapsis $\zeta=g(x_{\rm i})$ is given by

 \begin{displaymath}%
f(x_{\rm i})=\frac{M_{\rm p}(\zeta)}{M_{\rm p}(\zeta)+M_{\rm add}(\zeta)}\cdot
\end{displaymath} (8)

The final radius of the shell is $x=f(x_{\rm i})\zeta$ and mass conservation leads to the following final density profile


 \begin{displaymath}%
\rho(x)=
\rho_{\rm ta}(\zeta)f^{-3}(x_i)
\left[1+ \frac{{\rm d}\,\ln f(x_i)}{{\rm d}\,\ln g(x_i)} \right]^{-1}.
\end{displaymath} (9)

The potential energy of the system $W_{\rm f}$ in the relaxed state is related to its total energy by the virial theorem $W_{\rm f}=2E$. The energy of the system in a radial collapse is that at the turnaround epoch when all shells are assumed to have zero velocities at the same time. Thus $W_{\rm f}=2W_{\rm ta}$. A simple collapse factor satisfying this requirement is f=0.5 for every shell, leading to similarity solutions (a final density profile parallel in log-log space to that of the turnaround epoch). However, the collapse factor is not constant. N-body simulations (e.g. Voglis et al. 1995, hereafter VHH) show that f is an increasing function of the initial radius (or the turnaround radius) of the shell and its form is related to the initial profile of the density perturbation.

In a radial collapse case Eq. (8) gives $f\rightarrow 0$as $x_{\rm i}\rightarrow 0$ resulting in very condensed central regions with very steep density profiles (Lokas 2000a). As it is shown by Lokas & Hoffman (2000b) the inner slope of the density profile is between 2 and 2.3 even for low initial density peaks, far from the values derived from N-body simulations that range from 1 to 1.5 for the inner regions. On the other hand, the collapse factor in N-body simulations is hardly less than 0.05, even in the very central regions (VHH). This difference is easy to understand because no radial collapse exists in N-body simulations. During the expansion and the early stage of collapse particles acquire random velocities that prevent them from penetrating the inner regions. The consequence is the reduction of $M_{\rm add}$ in (8) that leads to larger values of f. In this sense the approximations based on constant f may be closer to the results of N-body simulations. The fit of the NFW profile is characteristic using an approximation given by del Popolo et al. (2000). However, in such an approximation the final density profile depends only on the "turnaround density profile'' and not on the amount of angular momentum acquired by the system during its expansion phase.

3 The numerical method and the angular momentum

The calculation of the collapse factor requires the evaluation of the integral in (6). Changing the variables from the turnaround radius to the initial one, this is written:

 \begin{displaymath}%
M_{\rm add}(\zeta)=4\pi\rho_{{\rm b}, {\rm i}}
\int_{x_{\rm...
...})[1+\delta_{\rm i}(x'_{\rm i})]x_{\rm i}'^2{\rm d}x'_{\rm i},
\end{displaymath} (10)

where $P_{x_{\rm i}}(x'_{\rm i})=I(x_{\rm i})/I(x'_{\rm i})$ with

 \begin{displaymath}%
I(r)=\int_{x'_{\rm p}}^r\frac{1}{v_{{\rm g}(x'_{\rm i})}(g(n))}\frac{{\rm d}g(n)}{{\rm d}n}{\rm d}n,\\
\end{displaymath} (11)

$\zeta=g(x_{\rm i})$, $x'_{\rm p}=g^{-1}(x_{\rm p})$ where $x_{\rm p}$ is the pericenter of the shell with initial radius $x'_{\rm i}$. The upper limit $x_{\rm b}$ of the integral in (10), is taken to be the initial radius of the sphere that has collapsed at the present epoch. The radial velocity v of a shell with apapsis $x=g(x_{\rm i})$ as it reaches the radius $r=g(r_{\rm i})$ is given by the conservation of the energy of the shell and is:

\begin{displaymath}%
v^2_x(r)=2[\Psi(r)-\varepsilon_x]-\frac{j^2_x}{r^2},\\
\end{displaymath} (12)

where $\Psi$ equals minus the potential $\Phi$, $\varepsilon_x$equals minus the specific energy and jx is the specific angular momentum of the shell. The potential $\Psi$ after the change of variables is given by the expression:

\begin{displaymath}%
\Psi[g(r_{\rm i})]=\frac{GM(x_{\rm b})}{g(x_{\rm b})}+G\int...
...rac{{\rm d}g(x_{\rm i})}{{\rm d}x_{\rm i}}{\rm d}x_{\rm i},\\
\end{displaymath} (13)

where the distribution of mass $M(x_{\rm i})$ is that at the initial conditions. The energy of the shell is calculated by:

\begin{displaymath}%
\varepsilon_x=\Psi[g(x_{\rm i})]-\frac{j^2_x}{2g^2(x_{\rm i})}\cdot
\end{displaymath} (14)

The angular momentum is introduced by the following scheme:

Each shell expands radially from its initial radius $x_{\rm i}$ up to its maximum expansion radius x. At this stage a specific angular momentum jx is added, given by $j_x=\mathcal{L}\sqrt{M(x)x}=\mathcal{L}\sqrt{M(x_{\rm i})g(x_{\rm i})}$, where $\mathcal{L}$ is a constant. This way of introducing angular momentum is consistent with the angular momentum distribution in N-body simulations (e.g. Barnes & Efstathiou 1987) and does not introduce any additional physical scale. It has been used by Avila-Reese et al. (1998) and recently by Nusser (2001).

In this way the apocenter $r_{\rm a}$ of a shell is its turnaround radius x while its pericenter $r_{\rm p}$ is found by the solution of the equation:

\begin{displaymath}%
2r^2[\Psi(r)-\varepsilon_x]-j^2_x=0.
\end{displaymath} (15)

The change of variables described above requires the solution for r of the equation

\begin{displaymath}%
2g^2(r)[\Psi(g(r))-\varepsilon_x]-j^2_x=0,
\end{displaymath} (16)

where the potential $\Psi$ is that of the turnaround epoch.

Nusser (2001) proved the following two important properties:

1.
If the angular momentum is introduced in the above described way in a spherical system with a power law density profile, then all shells have the same eccentricity. In fact, if $\varphi \equiv r_{\rm p}/r_{\rm a}$ then the following equation holds;

\begin{displaymath}G{\mathcal{L}}^2(\varphi^{-2}-1)=\frac{2[\Psi(r_{\rm p})-\Psi(r_{\rm a})]r_{\rm a}}{M(r_{\rm a})}\cdot
\end{displaymath} (17)

The right hand side of the above equation can be expressed in terms of $\varphi$ in the case of a power law density profile and completes the proof.
2.
If the potential evolves adiabatically, given at time t by the relation $\Psi_{\rm t}(r)=k(t)\Psi(r)$ with k(t) a slowly varying function of t and $\Psi$ the potential at the turnaround epoch, and the radial action is indeed invariant, then the eccentricity of every shell remains constant during the evolution.

The radial action can be written in the form:

\begin{displaymath}J_{\rm r}=j\times
\end{displaymath}


\begin{displaymath}\int_{\varphi}^1\left[(\varphi^{-2}-1)
\frac{\Psi(ur_{\rm a})...
...p})-\Psi(r_{\rm a})}
+(1-u^{-2})\right]^{\frac{1}{2}}{\rm d}u.
\end{displaymath} (18)

In the case of a power law density profile the quantity $\frac{\Psi(ur_{\rm a})-\Psi(r_{\rm a})}{\Psi(r_{\rm p})-\Psi(r_{\rm a})}$ is written in terms of $\varphi$ and u. Since $J_{\rm r}$ is constant, then $\varphi$ is constant.
Nusser (2001) used these properties to estimate the asymptotic behavior of the density profile near and far from the center of a system with a power-law initial density profile. Unfortunately these properties do not hold for more realistic density profiles. Similar power-law density profiles have been used by Sikivie et al. (1997) who used a CDM power spectrum to estimate an "effective'' exponent for this power-law on the galactic scale. Sikivie et al. use a self-similar evolution of the system in order to calculate its properties at the relaxed state. Unfortunately, the self-similar evolution is not valid for more realistic initial density profiles (non power-law profiles). The numerical method presented in our study is more general. It is able to deal with various initial density profiles (scale or not scale free) and the assumption of self-similarity is not required. Our results (presented in Sect. 4) have derived for realistic initial density profiles that have a finite value of the density perturbation at the location of the peak. Therefore, these results could be more reliable at least regarding the final state of the central region of the system. Moreover, it is shown that taking into account the angular momentum, final density profiles are well fitted by two power law density profiles with slopes less than 2at the central regions of the systems and larger than 2 at the outer regions. This class of final density profiles is consistent with the results of N-body simulations (e.g. Subramanian et al. 2000).


  \begin{figure}
\par\includegraphics[width=8.6cm,clip]{hiotfg1.eps}\end{figure} Figure 1: Overdensity profiles versus the distance x from a $3\sigma $ peak. Dotted curve corresponds to a smoothing length of $1.2~{\rm Mpc}$ while the dashed and solid curves correspond to smoothing lengths of 0.6 and $0.3~{\rm Mpc}$ respectively. The values of $<\delta (x)>$ and x are both normalized to the present.
Open with DEXTER

4 The initial density profile and the results

The averaged overdensity profile $<\delta (x)>$ at distance xfrom a $n\sigma=n\xi^{1/2}(0)$ extremum of a smoothed density is given in Bardeen et al. (1986, hereafter BBKS) by the equation:

\begin{displaymath}%
<\delta(x)> =\frac{n\xi(x)}{\xi(0)^{1/2}}
\end{displaymath}


\begin{displaymath}-\frac{\theta(n\gamma,\gamma)}{\gamma(1-{\gamma}^2)}
[{\gamma}^2\xi(x)+\frac{R^2_*}{3}\nabla^2\xi(x)]/{\xi(0)}^{1/2},
\end{displaymath} (19)

where $\gamma\equiv I(4)/[I(2)I(6)]^{1/2}$ and $R_*\equiv[3I(4)/I(6)]^{1/2}$ with $I(l)=\int_0^{\infty}k^lP(k){\rm d}k.$

In the above relations $\xi$ is the correlation function and Pthe power spectrum. These are related by:

\begin{displaymath}%
\xi(x)=\frac{1}{2{\pi}^2x}\int_0^\infty P(k)k{\rm sin}(kx){\rm d}k.
\end{displaymath} (20)

The function $\theta(n\gamma,\gamma)$ is given by the relation

\begin{displaymath}%
\theta(n\gamma,\gamma)=\frac{3(1-{\gamma}^2)+(1.216-0.9{\ga...
....45+(\frac{n\gamma}{2})^{2}]^{\frac{1}{2}}+
\frac{n\gamma}{2}}
\end{displaymath} (21)

in the range $1< n\gamma < 3$. The rms mass excess, $\frac{\Delta
M}{M}$, within a sphere of radius x, is given by:

\begin{displaymath}%
\sigma_{x}=\frac{1}{(2\pi^{2})^{\frac{1}{2}}}\frac{3}{x^3}
...
...\sin kx-kx\cos
kx)^{2}}{k^4}{\rm d}k\right]^{\frac{1}{2}}\cdot
\end{displaymath} (22)

We used the spectrum calculated by BBKS for a CDM-dominated Universe with $\Omega=1$ and h=0.5. This is given by the equation

\begin{displaymath}%
P_{\rm CDM}=Ak^{-1}[\ln(1+4.164k)]^2[G(k)]^{-\frac{1}{2}}
\end{displaymath} (23)

where G(k) is the following polynomial of k.
$\displaystyle %
G(k)=192.9+1340k+1.599\times10^5k^2$      
$\displaystyle +1.78\times10^5k^3+3.995\times10^6k^4.$     (24)

The above spectrum is smoothed on various scales according to the relation:

\begin{displaymath}%
P(k)=P_{\rm CDM}{\rm e}^{-(\frac{k}{k_{\rm c}})^2}.
\end{displaymath} (25)

Regarding the smoothing of the above spectrum three cases are examined. In case A the spectrum is smoothed on a scale ${k_{\rm c}}^{-1}=1.2$ Mpc, in case B on a scale ${k_{\rm c}}^{-1}=0.6$ Mpc and in case C on ${k_{\rm c}}^{-1}=0.3$ Mpc. The first scale length corresponds to a mass $0.5\times10^{12}\,M_{\odot}$ the second to a mass $6.25\times10^{10}\,M_{\odot}$ and the third to $7.8\times10^{9}\,M_{\odot}$. The constant of proportionality A is chosen by the condition $\sigma_{8}=0.7$, (the rms mass excess within 8 Mpc to be 0.7). Then $<\delta (x)>$ is calculated for n=3 and plotted in Fig. 1. The dotted curve corresponds to the case A, the dashed curve to the case B, while the solid one to the case C. We note that the values presented in this figure are normalized to the present.

In the case of linear growth of the overdensity the following hold (Gunn & Gott 1972): a shell with initial velocity equal to the Hubble flow and an initial comoving radius x has expanded up to a maximum radius $r_{\rm max}$ in a time $t_{\rm ta}$ (given by (3), $\Omega_{\rm i}=1$ in our case). This maximum radius is given by the relation

 \begin{displaymath}%
r_{\rm max}=x/\Delta,
\end{displaymath} (26)

while the collapse time of the shell, $t_{\rm c}$, is related to the age of the Universe, t0, by

 \begin{displaymath}%
t_{\rm c}=\frac{3\pi}{2}t_0{\Delta}^{-3/2}.
\end{displaymath} (27)

However calculating

 \begin{displaymath}%
\Delta=\frac{3}{{x}^3}\int_0^{x}<\delta(u)>u^2{\rm d}u,
\end{displaymath} (28)

and using (27) and (28) the collapse time of a shell and its turnaround radius are found. In fact the condition $t_{\rm c}/t_0=1$ gives the value of $\Delta$ for the shell that collapses today. Then x is found by solving numerically (28) and $r_{\rm max}$ is calculated by (26). In case A the mass inside the shell that collapses today is about $9.8\times10^{12}\,M_{\odot}$ and the radius of maximum expansion is 1130 Kpc. For the case B the mass is $4.3\times10^{12}\,M_{\odot}$and the radius of maximum expansion is 870 Kpc while for the case C the values are $2.1\times10^{12}\,M_{\odot}$ and 685 Kpc respectively. The values of x resulting from the numerical solution of (28) are 3.18, 2.45 and 1.93 Mpc respectively. Finally, the initial conditions at redshift $z_{\rm i}$can be derived by dividing both x and $<\delta (x)>$ by $1+z_{\rm i}$. Using $1+z_{\rm i}=1000$ the value of $x_{\rm b}$ in (10) are 3.18, 2.45 and 1.93 Kpc for the cases A, B and C respectively.

The amount of the angular momentum in the system is adjusted by the value of $\mathcal{L}$, (see Sect. 3), and is measured by the value of the dimensionless spin parameter

\begin{displaymath}%
\lambda\equiv\frac{L\vert E\vert^{\frac{1}{2}}}{GM^{\frac{5}{2}}}
\end{displaymath} (29)

where L, E and M are the total angular momentum, the total energy and the total mass of the system respectively. The mean value of $\lambda $ resulting from N-body simulations (e.g. Efstathiou & Jones 1979; Barnes & Efstathiou 1987) seems to be about 0.05. In our calculations we used different values of $\lambda $ with a maximum of 0.12. The maximum value of $\lambda $ corresponds to $\mathcal{L}=0.26$.

The resulting density profiles are fitted by a two-power law curve of the form

\begin{displaymath}%
\rho_{\rm fit}(r)=\frac{\rho_{\rm c}}{(\frac{r}{r_{\rm s}})^{\beta}(1+\frac{r}{r_{\rm s}})^{\mu}}
\end{displaymath} (30)

where the fitting parameters $\rho_{\rm c}$, $r_{\rm s}$, $\beta$ and $\mu$are calculated finding the minimum of the sum

\begin{displaymath}%
S=\sum_{i=1}^{NP}\left(\frac{\log\rho_{\rm SIM}(r_i)-\log\rho_{\rm fit}(r_i)}{\log\rho_{\rm SIM}(r_i)}\right)^2
\end{displaymath} (31)

where $\rho_{\rm SIM}$ are the predictions of SIM. NP is the number of points where the density is found. We used 100 points equally spaced on a log scale. The estimation of the above fitting parameters is done using the unconstrained minimizing subroutine ZXMWD of IMSL mathematical library. The quality of fit is very good as can be seen in the following three figures.


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg2.eps}\end{figure} Figure 2: Case A. Final density profiles derived for three different values of the spin parameter $\lambda $. Solid curves: SIM predictions. Dotted curves: the fits by a two-power law. From the bottom of the figure the curves correspond to $\lambda = 0.09, 0.05$ and 0.0 (radial collapse) respectively.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg3.eps}\end{figure} Figure 3: Case B. As in Fig. 2. Solid curves: density profiles derived from the SIM. From the higher to the lower curve the values of the spin parameter are 0.0, 0.05, 0.09 and 0.12respectively. Dotted curves: The fits of the solid curves by a two-power law density profile.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg4.eps}\end{figure} Figure 4: Case C. As in Fig. 2. From the higher to the lower solid line the values of $\lambda $ are 0.0, 0.05, 0.09 and 0.12 respectively. Dotted lines are the fits of the solid lines by a two-power law density profile.
Open with DEXTER

In Fig. 2 the final density profiles for the case A are shown for three different values of $\lambda $. The radial collapse ($\lambda=0$) corresponds to the higher solid curve. The intermediate curve corresponds to $\lambda=0.05$ while the lower one to $\lambda=0.09$. Dotted lines are the fits by the above described two-power law density profile. Distances are normalized to the virial radius $r_{\rm vir}$ which is 448 kpc while densities are normalized to the critical density $\rho_{\rm crit}$. The virial radius is the radius of the sphere with mean density $\approx$178times the present density of the Universe $\rho_{\rm crit}$ (Cole & Lacey 1996). The virial mass of the system (the mass contained inside the virial radius) is about $4.7\times10^{12}\,M_{\odot}$. It is characteristic that the efficiency of the angular momentum leads to shallower density profiles in the inner regions of the system. Additionally at the outer regions the density profile does not change even for the maximum amount of angular momentum used. We note that for larger values of the spin parameters the density profile becomes unrealistic (increasing at the inner regions). As it will be shown below this is a consequence of the shallow initial profile of this case.

The results for the case B are presented in Fig. 3. The values of the spin parameter are 0., 0.05, 0.09, and 0.12. The virial radius of the system is about 347 ${\rm Kpcs}$ and contains a mass of about $2.16\times10^{12}\,M_{\odot}$.

The case C is presented in Fig. 4 for the same values of $\lambda $ as in the case B. The virial radius of this system is 273 Kpc and its virial mass $1.06\times10^{12}\,M_{\odot}$.

Note that in the cases B and C the profile of the density decreases even for larger values of $\lambda $ than used in case A, because of their steeper initial density profiles.

In the following three figures the collapse factors for each case are presented. Figure 5 shows the collapse factor f of mass M on a logarithmic scale. The role of angular momentum is clear. Larger values of $\lambda $ lead to smaller values of f and consequently to shallower density profiles as shown in Figs. 2-4. Figure 6 refers to case B and Fig. 7 to case C. It is clearly shown in the above three figures that the collapse factor at the outer regions of the system is not affected by the amount of the angular momentum and it is almost the same as that of the radial collapse case. The efficiency of angular momentum in creating shallow density profiles depends on the initial density profile. This can be shown in the following three Figs where the slope of the two-power law fit versus radius is plotted. This is given by the relation

\begin{displaymath}%
\alpha(r)=-\frac{{\rm d}\ln\rho_{\rm fit}(r)}{{\rm d}\ln
r}=\beta+\mu\frac{\frac{r}{r_{\rm s}}}{1+\frac{r}{r_{\rm s}}}\cdot
\end{displaymath} (32)

Figure 8 corresponds to case A. The higher line is the slope resulting after a radial collapse ( $\lambda=0.0$) where the values of $\alpha$ in the interval $0.001r_{\rm vir}$ to $r_{\rm vir}$ are in the range 2. to 2.25. The intermediate line corresponds to $\lambda=0.05$ while the lower line corresponds to $\lambda=0.09$. Figure 9 shows the slope for the case B. The lines, from the higher to the lower, correspond to $\lambda=0.0, 0.05, 0.09$ and 0.12 respectively. The results of case C are shown in Fig. 10 for the same values of $\lambda $ as in case B. The values of $\alpha$ for the three cases at $r=0.001r_{\rm vir}$ are clearly shown in the above figures. At $r=0.01r_{\rm vir}$ and for $\lambda=0.05$ the values of $\alpha$ are 1.76, 1.82 and 1.83for the cases A, B and C respectively. A similar trend for the inner regions - smaller a for shallower initial density profile - is also clear for all values of $\lambda $. This trend is reversed at the outer regions. The slopes at $r=r_{\rm vir}$ are approximately 2.25, 2.20 and 2.15 for the three cases respectively. It is also clear from Fig. 10 that the radial collapse case leads to an almost exact power law profile. In this case the slope at $r=0.01r_{\rm vir}$ is 2.12 while at $r=r_{\rm vir}$ is 2.15.


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg5.eps}\end{figure} Figure 5: The collapse factor of mass M versus M on a logarithmic scale for case A. From the lower line the values of $\lambda $ are 0.0, 0.05, 0.09 respectively.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg6.eps}\end{figure} Figure 6: As in Fig. 5 but for case B. The values of $\lambda $, from the lower to the higher line are 0.0, 0.05, 0.09, and 0.12 respectively.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg7.eps}\end{figure} Figure 7: As in Fig. 5 but for the case C.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg8.eps}\end{figure} Figure 8: The slope of the final density profile versus radius for case A. Distance is normalized to the virial radius. From the lower to the higher line the values of $\lambda $ are 0.09, 0.05and 0.0 respectively.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg9.eps}\end{figure} Figure 9: As in Fig. 8 but for the case B. From the lower to the higher line the values of $\lambda $ are 0.12, 0.09, 0.05and 0.0 respectively.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{hiotfg10.eps}\end{figure} Figure 10: As in Fig. 8 but for the case C. From the lower to the higher line the values of $\lambda $ are 0.12, 0.09, 0.05and 0.0 respectively.
Open with DEXTER

5 Conclusions

The predictions of the SIM presented in this paper are summarized as follows:

1.
Radial collapse does not lead to power law final density profiles. However the difference in the slopes between the inner and the outer regions of the system is not large. Decreasing the smoothing scale of the power spectrum (leading to steeper initial profiles) this difference becomes smaller, leading to an almost power law for steep enough initial density profiles. The slopes in the radial collapse case are, in agreement with theoretical predictions (HS, Bertschinger 1985), in the range 2 to 2.25.
2.
Angular momentum leads to shallower inner density profiles. The inner slope depends on the amount of the angular momentum, measured in our results by the value of the spin parameter, and on the form of the initial density profile. Angular momentum becomes more efficient, in decreasing $\alpha$, for shallower initial density profiles.
3.
The slope of density profiles does not change significantly at the outer regions of the system even in cases where a large amount of angular momentum is assigned to the system. At $r=r_{\rm vir}$ the slope is approximately that of the radial collapse case.
We note that the above results are limited by a large number of assumptions, by the specific underlying cosmology and the particular form of the power spectrum used. However, they show systematic trends that could help us to better understand the relation between the initial conditions and the final density profiles. If things go the way described above, then the results of N-body simulations could be approximated by adding angular momentum to a case where the radial collapse results in a r-3density profile at the outer regions. However, the role of different parameters of the problem is under study.

Acknowledgements
Thanks to the Empirikion Foundation for its support

References

 
Copyright ESO 2002