Velocity Structure of Self-Similar Spherically Collapsed Halos

# Velocity Structure of Self-Similar Spherically Collapsed Halos

Phillip Zukin    Edmund Bertschinger Department of Physics, MIT, 77 Massachusetts Ave., Cambridge, MA 02139
###### Abstract

Using a generalized self-similar secondary infall model, which accounts for tidal torques acting on the halo, we analyze the velocity profiles of halos in order to gain intuition for N-body simulation results. We analytically calculate the asymptotic behavior of the internal radial and tangential kinetic energy profiles in different radial regimes. We then numerically compute the velocity anisotropy and pseudo-phase-space density profiles and compare them to recent N-body simulations. For cosmological initial conditions, we find both numerically and analytically that the anisotropy profile asymptotes at small radii to a constant set by model parameters. It rises on intermediate scales as the velocity dispersion becomes more radially dominated and then drops off at radii larger than the virial radius where the radial velocity dispersion vanishes in our model. The pseudo-phase-space density is universal on intermediate and large scales. However, its asymptotic slope on small scales depends on the halo mass and on how mass shells are torqued after turnaround. The results largely confirm N-body simulations but show some differences that are likely due to our assumption of a one-dimensional phase space manifold.

###### pacs:
04.50.Kd,04.50.-h,95.36.+x

## I Introduction

Recent N-body simulations have revealed a wealth of information about the velocity structure of halos Navarro et al. (2010); Ludlow et al. (2010); Vogelsberger et al. (2010). However, simulations have finite dynamic range. Moreover, it is difficult to draw understanding from their analysis, and computational resources limit the smallest resolvable radius, since probing smaller scales require using more particles and smaller time steps. Hence, it seems natural to use analytic techniques, which do not suffer from resolution limits, to analyze the velocity distributions of halos.

Numerous authors have analytically investigated the density profiles of halos. Work began with Gunn and Gott where they analyzed the continuous accretion of mass shells onto an initial overdensity Gunn and Gott (1972); Gott (1975); Gunn (1977). This process is known as secondary infall. By imposing that the mass accretion is self-similar, Fillmore and Goldreich Fillmore and Goldreich (1984) and Bertschinger Bertschinger (1985), assuming purely radial orbits, were able to analytically calculate the asymptotic slope of the density profile. Since then, there have been numerous extensions, some which do not assume self-similarity, that take into account non-radial motions Ryden and Gunn (1987); Nusser (2001); Hiotelis (2002); Williams et al. (2004); Sikivie et al. (1997); Del Popolo (2009); White and Zaritsky (1992); Le Delliou and Henriksen (2003); Ascasibar et al. (2004). Those works that do not impose self-similarity can only infer information about the velocity dispersion using the virial theorem. Hence they cannot predict a halo’s velocity anisotropy. Those works that do impose self-similarity focus only on the asymptotic slopes of density profiles.

In this paper, we analytically and numerically analyze the velocity structure of halos using an extended self-similar secondary infall model Zukin and Bertschinger (2010). The work which introduced this extended infall model will hereafter be referred to as Paper I. We then compare the predictions of our halo model to simulation results, focusing on the velocity anisotropy Binney and Tremaine (2008) and pseudo-phase-space density profiles Taylor and Navarro (2001); Navarro et al. (2010); Ludlow et al. (2010).

Density profiles do not uniquely determine a self gravitating system. In order to more fully characterize dark matter halos, one needs to probe their phase space distributions. The velocity anisotropy and pseudo-phase-space density profiles are thereby useful since they complement density profiles by revealing additional information about the phase space structure of the halo.

In section II we summarize our generalized secondary infall model and discuss how to numerically calculate the radial and tangential velocity dispersions in the halo. In section III we analytically calculate the asymptotic behavior of the radial and tangential kinetic energy profiles on small and intermediate scales. In section IV we compare our numerically calculated anisotropy and pseudo-phase-space density profiles to recent N-body results and conclude in section V.

## Ii Self-Similar Model

Here we first summarize the self-similar halo formation model developed in Paper I. The model is characterized by four parameters which are described below.

In this model, the universe is initially composed of a linear spherically symmetric density perturbation with mass shells that move approximately with the hubble flow. Because of the central overdensity, mass shells eventually stop their radially outward motion and turn around. The radius at which a mass shell first turns around, or its first apocenter, is known as the turnaround radius. Since the average density is a decreasing function from the central overdensity, mass shells initially farther away will turnaround later. The halo grows by continuously accreting mass shells. Mass shells are labeled by their turnaround time or their turnaround radius .

Model parameter characterizes how quickly the initial linear density field falls off with radius (). It is related to the effective primordial power spectral index () through Hoffman and Shaham (1985). Since depends on scale, is set by the halo mass. As in Paper I, we restrict our attention to so that the initial density field decreases with radius while the excess mass increases with radius. Note however that corresponds to objects larger than galaxy clusters today. Model parameter also sets the growth of the current turnaround radius: where Zukin and Bertschinger (2010).

Self-similarity imposes that at time , the angular momentum per unit mass of a particle in a shell at , and the density and mass of the halo have the following functional forms Zukin and Bertschinger (2010).

 L(r,t) = Br2ta(t)tf(λ,t/t∗) (1) ρ(r,t) = ρB(t)D(λ) (2) M(r,t) = 4π3ρB(t)r3ta(t)M(λ) (3)

where is the radius scaled to the current turnaround radius, is the background density for an Einstein de-Sitter (flat ) universe, and is a constant. Inspired by tidal-torque theory and numerical simulations, we take to be:

 f(λ,t/t∗)={λ(4−p)/2if tt∗. (4)

Model parameter , defined above, sets how quickly angular momentum builds up before turnaround while sets the amplitude of angular momentum at turnaround. In Paper I, using cosmological linear perturbation theory, we constrained and so that the angular momentum of particles before turnaround evolves as tidal torque theory predicts Hoyle (1951); Peebles (1969); White (1984); Doroshkevich (1970). Conveniently both and are set by the halo mass. However, after comparing to density profiles from N-body simulations, we found that our expression for derived from linear theory overestimates the actual value. Hence, for the rest of this paper, the notation () signifies using a value of divided by 1.5 (2.3).

Model parameter , defined above, sets how quickly the angular momentum of particles grows after turnaround. This parameter is difficult to constrain analytically since the halo is nonlinear after turnaround. However, in Paper I we showed that a nonzero can be sourced by substructure. Moreover, influences the density profile at small scales since it controls how the pericenters of shells evolve over time.

The trajectory of a shell after turnaround contains all of the velocity information in the halo. The trajectory’s evolution equation, which follows from Newton’s law, is given by:

 d2λdξ2+(2β−1)dλdξ+β(β−1)λ=−29M(λ)λ2+B2λ3e2(ϖ+1−2β)ξ (5)

where and is the current turnaround time. The initial conditions for eq. (5) are and . Calculating requires evolving both the shell’s trajectory and before turnaround Zukin and Bertschinger (2010). Because of self-similarity, the trajectory can either be interpreted as labeling the location of a particular mass shell at different times, or labeling the location of all mass shells at a particular time. We take advantage of the second interpretation in order to numerically calculate the velocity profiles.

Inside the shell that is currently at its second apocenter, multiple shells exist at all radii. This can be seen from Figure 5 of Bertschinger (1985), which plots the location of all shells at a particular time. Therefore the expectation value of a quantity , for example the radial velocity, at radius and time is the value of for each shell at weighted by each shell’s mass. We find:

 ⟨h(r,t)⟩=∫Mta0dM∗Mtah(t,t∗)δD(λ−λ(ξ))∫Mta0dM∗MtaδD(λ−λ(ξ)) (6)

where is the current turnaround mass, represents the value of for the shell with turnaround time , is the mass of the shell with turnaround time , and is the dirac delta function which picks out all shells at .

In order to numerically calculate , we must relate to , the computed trajectory of the shell which turns around at . Using self-similarity, we find:

 vt(t,t∗) = Brtate(1+ϖ−2β)ξ∗1λ(ξ∗) (7) vr(t,t∗) = rtate−βξ∗ddξ[eβξλ(ξ)]|ξ=ξ∗ (8)

where and () is the tangential (radial) velocity. Using eqs. (6)-(8), and taking advantage of the delta function, the tangential () and radial () velocity dispersions become:

 σ2t(r,t) ≡ ⟨v2t(r,t)⟩ (9) = r2tat2∑ie(4−7β+2ϖ)ξiλ−2i|dλ/dξ|−1i∑ie(2−3β)ξi|dλ/dξ|−1i σ2r(r,t) ≡ ⟨v2r(r,t)⟩−⟨vr(r,t)⟩2 = r2tat2∑ie(2−5β)ξi[d(eβξλ)/dξ]2i|dλ/dξ|−1i∑ie(2−3β)ξi|dλ/dξ|−1i − r2tat2⎛⎝∑ie(2−4β)ξi[d(eβξλ)/dξ]i|dλ/dξ|−1i∑ie(2−3β)ξi|dλ/dξ|−1i⎞⎠2

where is the th root that satisfies . In the above, we’ve imposed since our model assumes that the orbital planes of particles in a given shell are oriented in random directions. Note that inside the shell that is currently at its second apocenter, interference between multiple shells forces to quickly go to zero.

## Iii Asymptotic Behavior

Here, using techniques developed in Fillmore and Goldreich (1984), we analytically calculate the logarithmic slope of the tangential and radial kinetic energy in two different radial regimes. We accomplish this by taking advantage of adiabatic invariance and self consistently calculating the total radial and tangential kinetic energy profiles of the halo.

We start by parametrizing the halo mass, radial kinetic energy , tangential kinetic energy and the variation of the apocenter distance .

 M(r,t) = κ(t)rα (11) Kr(r,t) = κr(t)rαr (12) Kt(r,t) = κt(t)rαt (13) rar∗ = (tt∗)q (14)

In the above is the turnaround radius of a mass shell that turns around at . As was shown in Paper I, adiabatic invariance allows us to relate and to . At late times, the orbital period is much smaller than the time scale for the mass and angular momentum to grow. Integrating Newton’s equation and assuming and change little over an orbit, we find:

 (drdt)2=2Gκ(t)α−1(rα−1a−rα−1)−L2(t)(r−2−r−2a). (15)

The above relationship tell us how the pericenters evolve with time. Defining and evaluating the above at , we find:

 1−yα−1y−2−1≡A(y)=(α−1)L2(t)2Gκ(t)rα+1a(t) (16)

In Paper I, by analyzing the radial action, we found that when , . Therefore, for , eq. (16) implies that will decrease over time. However, for , and eq. (16) imply that will increase over time and will at one point violate . Since we only consider bound orbits, the constraint holds. At late times, as the angular momentum continues to increase for , , orbits become approximately circular, the radial action vanishes, and . Hence halos with will have orbits that become more radial over time () while halos with will have orbits that become more circular over time (). The above insight leads to the following constraint.

 q=⎧⎨⎩1α+1{2ϖ+23n[α(1+n)−3]}%if$ϖ≥0$23n(α+1)[α(1+n)−3]if ϖ<0 (17)

For the specific case, , taking advantage of , the adiabatic invariance arguments above, and eqs. (1) and (4), we can rewrite eqn. (16) in the form , where:

 l={ϖif α>1,2ϖ/(α+1)if α<1. (18)

and is the pericenter of a mass shell at turnaround. The special case will be addressed later.

We next calculate the kinetic energy profiles. After a few orbits, shells oscillate at a much higher frequency than the growth rate of the halo. When calculating the internal mass profile, this allows us to weight each mass shell based on how much time it spends interior to a certain scale Fillmore and Goldreich (1984). Likewise, when calculating the total internal kinetic energy profile, we can weight each mass shell by both a time-averaged (or ) and a factor that accounts for how often the shell lies interior to a certain scale. For a derivation, please see the Appendix.

Using eqs. (15) and (16), the kinetic energy weighting at time for a mass shell with apocenter distance , pericenter , below , is:

 Pi(u,y) = 0(u1) (19)

where

 It(u,y) = 12(Br2∗t∗ra)2(tt∗)2ϖ∫uydvv2f(v,y) (20) Ir(u,y) ∝ r2tat2(rarta)α−1∫uyf(v,y)dv (21) I(u,y) = ∫uydvf(v,y) (22)

and

 f(v,y)≡{((1−vα−1)−A(y)(v−2−1))1/2if α>1,((vα−1−1)+A(y)(v−2−1))1/2if α<1. (23)

The index is used for shorthand to represent either the radial or tangential direction and the dependence of and on is implicit. Eq. (21) is not an equality since . Moreover, the proportionality constant varies for different radial regimes in the halo.

Similar to the treatment in Paper I, self consistency demands that:

 (rrta)αi = Ki(r,t)Ki(rta,t) ∝ ∫Mta0dM∗Mtat2r2taPi(rra(t,t∗),y(t,t∗))

where is the mass of a shell that turned around at and is the current turnaround mass. The above is not an equality, even for the tangential kinetic energy, because of a proportionality constant, similar to , that is not included. See the Appendix for details. Its numerically computed value does not affect the asymptotic slopes . Noting from eq. (3) that

 dM∗dlnt∗=(3β−2)Mta(tt∗)3β−2, (25)

using eq. (14) and transforming integration variables to , we find:

where . As increases, the above integral sums over shells with smaller . Since the pericenter of a shell evolves with time, the second argument of depends on . The dependence varies with torque model (sign of ); hence we’ve kept the dependence on implicit. Next we analyze the above for certain regimes of , and certain torquing models, in order to constrain the relationship between and .

### iii.1 r/rta≪y0, ϖ<0

For , particles lose angular momentum over time. When probing scales , mass shells with only contribute. As a result, . Using eq. (18), we then find:

 y(t,t∗)=y0(tt∗)l=y0(rurta)δ (27)

where and . For bound mass shells, . Therefore, since , the first argument of in eq. (26) increases while the second decreases as we sum over shells that have turned around at earlier and earlier times (). For , mass shells which most recently turned around do not contribute to the kinetic energy inside since we are probing scales below their pericenters. Mass shells only begin to contribute when the two argument of are roughly equal to each other. This occurs around:

 u=y1≡(y0(r/rta)δ)1/(1+δ) (28)

Hence, we can replace the lower limit of integration in eq. (26) with . We next want to calculate the behavior of eq. (26) close to in order to determine whether the integrand is dominated by mass shells around or mass shells that have turned around at much earlier times. The first step is to calculate the behavior of for . We find:

 Pt(u,y) ∝ r2tat2(rurta)ltu1/2(1−y/u)1/2 (29) × {y−3/2if α>1,y−1−α/2if α<1. Pr(u,y) ∝ r2tat2(rurta)lru3/2(1−y/u)3/2 (30) × {y−1/2if α>1,y−1+α/2if α<1.

where and . Given the above, we evaluate the indefinite integral in eq. (26), noting that is a function of (eq. 27). For , we find:

 t2r2ta∫duu1+kPt(u,y0(rurta)δ) ∝ (u/y1−1)3/2(rrta)lt×{y−1−k−lt1if α>1,y−1/2−k−lt−α/21if α<1. t2r2ta∫duu1+kPr(u,y0(rurta)δ) ∝ (u/y1−1)5/2(rrta)lr×{y1−k−lr1if α>1,y1/2−k−lr+α/21if α<1.

Following the logic in Paper I, if we keep fixed and the integrand blows up as , then the left hand side of eq. (26) must diverge in the same way as the right hand side shown in eqs. (III.1) and (III.1). Therefore, using eq. (28):

 αt−k−lt = {−δ(1+k+lt)/(1+δ)if α>1,−δ(1/2+k+lt+α/2)(1+δ)if α<1. (33) αr−k−lr = {δ(1−k−lr)/(1+δ)if α>1,δ(1/2−k−lr+α/2)(1+δ)if α<1.

Otherwise, if the integrand converges, then the right hand side is proportional to . Therefore, the left hand side must also have the same scaling, which implies . Solving the above system of equations for simplifies dramatically since we have already solved for in Paper I. Rewritten below for convenience, we found:

 For n ≤ 2: α = 1+n−√(1+n)2+9nϖ(nϖ−2)3nϖ k = 1+n+3nϖ−√(1+n)2+9nϖ(nϖ−2)nϖ(4+n) q = 1+n−3nϖ−√(1+n)2+9nϖ(nϖ−2)3n For n ≥ 2: α = k=31+n,q=0. (35)

Using eq. (III.1) to solve for and making sure the solution is consistent, (ie: using eqs. (33) and (III.1) only if the integrand diverges as ), we find:

 For n ≤ 2: αt = αr=12−9nϖ−2√(1+n)2+9nϖ(nϖ−2)1+n+√(1+n)2+9nϖ(nϖ−2) For n ≥ 2: αt = (4+n)(3nϖ+2n−10)2(1+n)(3nϖ−n−4) αr = 5−n1+n. (36)

The above solutions are continuous at . Taking the no-torque limit (), we find for all . Assuming virial equilibrium, one would predict . Hence, in the no-torque limit, both the radial and tangential kinetic energy are virialized. However, when , only the radial kinetic energy for is virialized. All other profiles are out of virial equilibrium because they are dominated by shells which recently turned around and hence have not had time to virialize. Since all collapsed objects today have , this model predicts unvirialized halos when particles lose angular momentum after turnaround.

### iii.2 r/rta≪y0, ϖ>0

For , the angular momentum of particles increase with time. As mentioned above, when probing scales , mass shells with only contribute. As a result, . In other words, the orbits are roughly circular. We can therefore replace the lower limit of integration in eq. (26) with 1 since mass shells will only start contributing to the sum when . Hence, the right hand side of eq. (26) is proportional to , which implies . Using the results from Paper 1, reproduced below for convenience,

 α=k=31+n−3nϖ,q=2ϖ,for0≤n≤3, (37)

we find:

 αt=αr=5−n+3nϖ1+n−3nϖ,% for0≤n≤3. (38)

The no-torque case, , is consistent with the analysis in the prior subsection. The singularity , as discussed in Paper I, corresponds to orbits that are not bound. Hence we only consider . Eq. (38) shows that the halo, for , is in virial equilibrium (). This is expected since the velocity profiles are dominated by mass shells that have turned around at .

### iii.3 y0≪r/rta≪1

In this regime, we are probing scales larger than the pericenters of the most recently turned around mass shells. As a result, is dominated by the contribution from the integrand when . Therefore:

 Pt(u,y) ∝ r2tat2(rurta)lt×{y−1if α>1,y−(α+1)/2if α<1. (39) Pr(u,y) ∝ r2tat2(rurta)lr×{uif α>1,u(α+1)/2if α<1.

Plugging in the above into eq. (26), using the results of Paper I shown below for convenience,

 α = 1,k=64+n,q=n−23n,% forn≤2 α = k=31+n,q=0,% forn≥2

and utilizing the same divergence and convergence arguments above, we find:

 For n ≤ 2: αt = 0,αr=1 For n ≥ 2: αt = {0if ϖ<5−n3n,5−n−3nϖ1+nif 5−n3n≤ϖ<1+n3n. αr = 5−n1+n. (42)

The above solutions are continuous at . The upper limit on enforces that orbits are bound (). For , results in unbound orbits and hence is not considered. The radial kinetic energy follows the same profile expected from virial equilibrium (), even though recently turned around mass shells dominate the kinetic energy for . We believe this is a result of angular momentum not playing a dynamical role at these scales. Using this logic, and as eq. (III.3) reveals, this is consistent with the the tangential kinetic energy not being virialized (). Taking the limit of eq. (III.1) as , we recover the same expressions as eq. (III.3) for . This is expected since in this limit, particles lose their angular momentum instantly, resulting in purely radial orbits. We do not expect to recover the same expressions for since the tangential kinetic energy vanishes in this limit.

Eqs. (III.1), (37) and (III.3) determine how quickly the mass inside a fixed radius grows as a function of time ( defined in eq. 11). When , apocenters of mass shells settle down to a constant fraction of their turnaround radii, leading to a constant mass inside a fixed radius. For cosmologically relevant structures (), this occurs on small scales for . When , inward migration leads to an increasing mass inside a fixed radius. When , outward migration leads to a decreasing mass inside a fixed radius.

This section assumed and yet, for certain parts of parameter space, eqs. (III.1), (37) and (III.3) give . However, since the solutions are continuous as from the left and right, then the results hold for as well.

## Iv Comparison with N-body Simulations

In this section, using the analytic results derived above to gain intuition, we first analyze how influences the anisotropy and pseudo-phase-space density profiles and then compare our numerically computed profiles to recent N-body simulations of galactic size halos Navarro et al. (2010).

As described in Paper I, the mass of a halo is not well defined when our model is applied to cosmological structure formation since it is unclear how the spherical top hat mass which characterizes the halo when it is linear relates to the virial mass which characterizes the halo when it is nonlinear. For halos today with galactic size virial masses, we assume the model parameter which characterizes the initial density field, is set by a spherical top hat mass of . Specifying the top hat mass also sets model parameters and . For explicit expressions used to calculate model parameters , , and , please see Paper I.

When analyzing the influence of , we use model parameters: , , and . When comparing to N-body simulations, we use model parameters , , , and . This value of ensures on small scales and, as shown in Paper I, this range in gives good agreement with the Einasto and NFW profiles Navarro et al. (1996).

For the N-body comparisons, we average , , and in 50 spherical shells equally spaced in over the range , where satisfies . This is the same procedure followed with the recent Aquarius simulation Navarro et al. (2010). We also calculate , the radius where reaches a maximum. This radius, as well as the virial radius , are commonly referred to in simulation papers. As discussed in Paper I, the density profile is isothermal for our halo over a range of . Moreover, the maximum peaks associated with the caustics are unphysical. So, we choose a value of in the isothermal regime that gives good agreement with empirical density profiles. Changing does not change our interpretation of the results. We find = 0.07 (0.05) for (). For reference, we find the dimensionless radius of first pericenter passage () to be 0.042 (0.026) for ().

As mentioned previously, N-body simulations have finite dynamic range. The innermost radius where the simulation results can be trusted is set by the total number of particles Power et al. (2003). The recent Aquarius simulations characterize their innermost radius based on the convergence of the circular velocity, at a particular radius, for the same halo simulated at different resolutions Navarro et al. (2010). The notation corresponds to the smallest radius such that the circular velocity has converged to or better at larger radii. When these radii are showed in the figures, we use the values quoted in Table 2 of Navarro et al. (2010) for halo Aq-A-2 ( and ) since all six halos were simulated at this resolution.

### iv.1 Anisotropy Profile

Here we analyze the velocity anisotropy for galactic size halos, where the tangential and radial velocity dispersions are defined in eqs. (9) and (II) respectively. Based on the analysis in Section III, we expect to asymptote to a constant for since and to increase for since . Moreover, for radii larger than the first shell crossing (), since only one shell contributes to the dispersion. Hence, in this radial range, .

In the top panel of Figure 1 we plot the velocity anisotropy for galactic size halos with varying . In the bottom panel, we plot the smoothed velocity anisotropy for model parameters that give good agreement with density profiles from simulated galactic size halos. The downward spikes in both panels are caustics which exist because of unphysical radially cold initial conditions. In both panels, as analytically predicted, the velocity anisotropy asymptotes at small radii, increases at intermediate radii, and then drops off near the virial radius.

The top panel shows that affects the radius of first pericenter passage (), the amplitude of close to the virial radius, as well as the asymptotic value of at small radii. This behavior is intuitive since smaller values of give rise to halos populated with less circularized orbits at a given radius. Note however that the envelope of the anisotropy profile begins to increase and become more radially dominated for , contradicting our analytic analysis. More specifically, for , for , and starts to increase at when it should start increasing at , according to Section III. This is a result of assumptions used to calculate breaking down. This is more apparent for since the orbital period is longer. However, as , the assumptions become more valid.

In the bottom panel, we show how model parameter affects the velocity anisotropy. As discussed in Paper I, smaller leads to orbits that take longer to circularize and density profiles with a larger isothermal region (smaller ). The bottom panel should be compared to Figures 9 and 10 of Navarro et al. (2010). Though our model cannot address structure outside , the graphs are qualitatively very similar. The width of the peak predicted in our model agrees with results from N-body simulations. This should be expected since the parameter was chosen so that the width of the isothermal region in the density profiles agree. However, our model over predicts the velocity anisotropy close to and under predicts the velocity anisotropy at small radii. In other words, at large scales the halo is populated with too many radial orbits while on small scales the halo is populated with too many circular orbits.

This trend is most clearly seen in Figure 2. There we plot the local velocity anisotropy versus the logarithmic slope of the density profile for a galactic size halo with as well as a universal relationship relating these two quantities that was derived by Hansen Moore Hansen and Moore (2006). The open circles correspond to while the filled circles correspond to . This figure should be compared to Figure 11 of Navarro et al. (2010). In the Aquarius simulation paper, the Hansen Moore prediction agrees well with N-body results for . However, in our Figure 2, while there is a clear trend between the local velocity anisotropy and the logarithmic slope of the density profile, that trend does not match the derived relationship. Note though that our model, just as the Aquarius simulation claimed, does show deviations from the Hansen Moore trend for . In our model, this deviation is caused by a vanishing radial velocity dispersion. For simulated halos, other effects like non-sphericity or non-self-similarity may also play a role.

The self-similar model’s inability to match the amplitude of the velocity anisotropy seen in N-body simulations reveals a weakness in the model. Clearly, it is unphysical for all particles in a particular shell to have the same amplitude of angular momentum and the same radial velocity. In reality, a given shell should have a radial velocity dispersion and should have a distribution of angular momentum that evolves with time. This possibility will be discussed again in Section V.

### iv.2 Pseudo-Phase-Space Density Profile

Here we analyze the pseudo-phase-space density profiles and for galactic size halos, where . Taylor and Navarro claimed that the pseudo-phase-space density roughly follows the power law for all halos Taylor and Navarro (2001). Surprisingly, this power law matches predictions made by Bertschinger for purely radial self-similar collapse onto a spherical top hat perturbation Bertschinger (1985). Taylor and Navarro’s claim has been verified numerically Navarro et al. (2010); Rasia et al. (2004); Dehnen and McLaughlin (2005); Faltenbacher et al. (2007); Vass et al. (2009); Wang and White (2009), however recently the highest resolution simulations have seen evidence for departures from this power law near their innermost resolved radii Ludlow et al. (2010).

Based on the analysis in Section III, we expect the power law exponent to depend on for . With model parameters which give for galactic size halos, the extended secondary infall model predicts . This is expected for a virialized halo () with . For , the model predicts if the radial velocity dispersion dominates and if the tangential velocity dispersion dominates.

In the left panels of Figure 3, we plot and for galactic size halos with varying . In the right panels, we plot the smoothed pseudo-phase-space densities, with the radius scaled by , for model parameters that give good agreement with density profiles from simulated galactic size halos. In addition, we overlay the radial top hat solution. Scaling the radius to causes the first pericenter () of both models to roughly overlap, leading to less difference in the amplitude of the pseudo-phase-space density at small radii.

The left panels show that the asymptotic slopes vary with . The numerically computed slopes match analytic predictions. The panels for blow up at radii close to the virial radius since vanishes. The right panels should be compared to Figure 13 of Navarro et al. (2010). The extended secondary infall model predicts that simulations of galactic size halos should see significant deviations from Taylor and Navarro’s claim at when analyzing and when analyzing . Looking at the residuals in Figure 13 of Navarro et al. (2010), this prediction seems plausible. If higher resolution simulations do not show deviations from Taylor and Navarro’s claim, then this secondary infall model would be proven incorrect since the model cannot consistently reproduce both the density and velocity profiles of simulated halos.

As shown in Section III, for cosmological initial conditions (), , and have power laws that are independent of initial conditions and torqueing parameters in the regime . This implies that the pseudo-phase-space density is universal on these scales. This universality on intermediate scales may have played a role in Taylor and Navarro’s initial claim.

Figure 4 shows the difference of the pseudo-phase-space density power law exponent from the radial top hat solution, on small scales, as a function of model parameters and . The range in corresponds to . The range in ensures that all orbits are bound. According to the extended secondary infall model, positive is necessary for in order to have on small scales Zukin and Bertschinger (2010). If all halos have on small scales, then halos with will have while halos with will have pseudo-phase-space density exponents that vary with halo mass. If on the other hand, is constant for all halos, then as Figure 4 shows, the power law will vary with mass.

## V Discussion

N-body simulations have revealed a wealth of information about the velocity profiles of dark matter halos. In an attempt to gain intuition for their results, we’ve used a generalized self-similar secondary infall model which takes into accounts tidal torques. The model assumes that halos self-similarly accrete radially cold mass shells. Moreover, each shell is composed of particles with the same amplitude of angular momentum. While the model is simplistic, it does not suffer from resolution limits and is much less computationally expensive than a full N-body simulation. Moreover, it is analytically tractable. Using this model we were able to analytically calculate the radial and tangential kinetic energy profiles for and , where is the dimensionless radius of first pericenter passage, is the virial radius, and is the current turnaround radius.

It is clear from our analysis that angular momentum plays a fundamental role in determining the velocity structure of the halo. The amplitude of angular momentum at turnaround sets the transition scale () between different power law behaviors in the tangential and radial kinetic energy profiles. Also, for collapsed objects today (), , the parameter that quantifies how particles are torqued after turnaround, influences the slopes of both the radial and kinetic velocity dispersions at small radii. Moreover, both the amplitude of angular momentum at turnaround and affect the asymptotic value of the velocity anisotropy profile at small radii.

For , the self-similar halo is not virialized on small scales since the radial and tangential kinetic energy is dominated by mass shells which have not had time to virialize. On the other hand, for , the halo is virialized since the dominant mass shell have had time to virialize. As shown in Paper I, requires for . Hence, positive is favored in order to reproduce N-body simulation density profiles. Quantifying requires analyzing N-body simulations and is beyond the scope of this work. However constraining with simulations will provide a test for this extended secondary infall model.

Our model predicts that the pseudo-phase-space density profile is universal on intermediate and large scales. This could potentially play a role behind the claimed universality of the pseudo-phase-space density Taylor and Navarro (2001). Since we do not understand how depends on halo mass, it is impossible to rule out universality on small scales, since can potentially conspire to erase initial conditions. However, if galactic size halos have , then regardless of universality, the model predicts . While hints of deviations from the radial top hat solution have been seen in recent simulations Ludlow et al. (2010), higher resolution simulations are needed to better test the model.

While our self-similar model has its clear advantages, it is also unphysical. First, all particles in a given mass shell have the same radial velocity. This leads to caustics. The same tidal torque mechanisms which cause a tangential velocity dispersion Hoyle (1951), should give rise to a radial velocity dispersion. Second, while qualitatively similar, the comparison of the model’s predicted velocity anisotropy to N-body simulation results reveals that our treatment of angular momentum is too simplistic. The model predicts too many radial orbits at large radii and too many circular orbits at small radii. In reality, each shell is composed of a distribution of angular momentum that evolves with time. In order to properly take these two effects into account, one would need a statistical phase space description of the halo that includes sources of torque. Ma Bertschinger provided such an analysis in the quasilinear regime Ma and Bertschinger (2004). Therefore, a natural extension of this secondary infall model, which could potentially reproduce both position and velocity space information of N-body simulations, would be to generalize Ma Bertschinger’s analysis to the nonlinear regime and impose self-similarity.

###### Acknowledgements.
The authors acknowledge support from NASA grant NNG06GG99G and thank Robyn Sanderson, Paul Schechter, and Mark Vogelsberger for many useful conversations.

## Appendix A Deriving the Consistency Relationship

In this Appendix, we derive eq. (III). Self-similarity imposes that the total radial (or tangential) kinetic energy at time contained within radius is given by:

 Ki(r,t)=Mta(t)r2tat2Ki(λ) (43)

where is used for shorthand to denote the radial or tangential direction. The kinetic energy also obeys the following relationship.

 Ki(r,t)=12∫Mta0dM∗v2i(t,t∗)H[r−R(t,t∗)] (44)

where is the mass, velocity, and radius of a shell at time which turned around at and is the heaviside function. Since, after a short time, shells begin to oscillate on a timescale much shorter than the growth of the halo, we can replace with a time averaged version and the heaviside function with a weighting that takes into account how often the shell is below . More specifically, considering a shell with turnaround time such that , we have:

 ⟨v2i(t,t∗)⟩