Self-Similar Spherical Collapse with Tidal Torque
aa r X i v : . [ a s t r o - ph . C O ] D ec Self-Similar Spherical Collapse with Tidal Torque
Phillip Zukin ∗ and Edmund Bertschinger Department of Physics, MIT, 77 Massachusetts Ave., Cambridge, MA 02139
N-body simulations have revealed a wealth of information about dark matter halos however theirresults are largely empirical. Using analytic means, we attempt to shed light on simulation resultsby generalizing the self-similar secondary infall model to include tidal torque. In this first of twopapers, we describe our halo formation model and compare our results to empirical mass profilesinspired by N-body simulations. Each halo is determined by four parameters. One parameter setsthe mass scale and the other three define how particles within a mass shell are torqued throughoutevolution. We choose torque parameters motivated by tidal torque theory and N-body simulationsand analytically calculate the structure of the halo in different radial regimes. We find that angularmomentum plays an important role in determining the density profile at small radii. For cosmologicalinitial conditions, the density profile on small scales is set by the time rate of change of the angularmomentum of particles as well as the halo mass. On intermediate scales, however, ρ ∝ r − , while ρ ∝ r − close to the virial radius. I. INTRODUCTION
The structure of dark matter halos affects our under-standing of galaxy formation and evolution and has im-plications for dark matter detection. Progress in our un-derstanding of dark matter halos has been made bothnumerically and analytically. Analytic treatments beganwith work by Gunn and Gott; they analyzed how boundmass shells that accrete onto an initially collapsed objectcan explain the morphology of the Coma cluster [1, 2] andelliptical galaxies [3]. This continuous accretion processis known as secondary infall.Secondary infall introduces a characteristic lengthscale: the shell’s turnaround radius r ∗ . This is the ra-dius at which a particular mass shell first turns around.Since the average density is a decreasing function of dis-tance from the collapsed object, mass shells initially far-ther away will turnaround later. This characteristic scaleshould be expected since the radius of a mass shell, likethe radius of a shock wave in the Sedov Taylor solution,can only depend on the initial energy of the shell, thebackground density, and time [4]. By imposing that thestructure of the halo is self-similar – all quantities de-scribing the halo only depend on the background density, r ta (the current turnaround radius), and lengths scaled to r ta – Bertschinger [4] and Fillmore and Goldreich (here-after referred to as FG) [5] were able to relate the asymp-totic slope of the nonlinear density profile to the initiallinear density perturbation.Assuming purely radial orbits, FG analytically showedthat the slope ν of the halo density distribution ρ ∝ r − ν falls in the range 2 < ν < .
25 for r/r ta ≪
1. This devi-ates strongly from N-body simulations which find ν ∼ < ν ∼ . ν ∼ . ∗ [email protected] while orbits in simulations and observed galaxies containtangential components, it is analytically tractable anddoes not suffer from resolution limits. Numerical darkmatter simulations, on the other hand, do not make anysimplifying assumptions and have finite dynamic range.Moreover, it is difficult to draw understanding from theiranalysis and computational resources limit the smallestresolvable radius, since smaller scales require more par-ticles and smaller time steps. It seems natural, then, togeneralize the work done by FG in order to explain thefeatures predicted in simulations and observed in galax-ies. This paper, in particular, investigates how non-radialmotion affects the structure of dark matter halos.Numerous authors have investigated how angular mo-mentum affects the asymptotic density profile. Ry-den and Gunn analyzed the effects of non-radial motioncaused by substructure [13] while others have examinedhow an angular momentum, or a distribution of angularmomenta, assigned to each mass shell at turnaround, af-fects the structure of the halo [14–21]. Note that many ofthese authors do not impose self-similarity. Those thatdo assume that a shell’s angular momentum remains con-stant after turnaround.This paper extends previous work by Nusser [14].Assuming self-similarity, he analytically calculated thestructure of the halo in different radial regimes for shellswith constant angular momentum after turnaround. Hefound that the inclusion of angular momentum allows0 < ν < .
25. According to Hoffman and Shaham [22], ν depends on the effective primordial power spectral in-dex ( d ln P/d ln k ), which varies for different mass halos.For galactic size halos, Nusser’s analytic work predicts ν ∼ .
3, in disagreement with simulation results [6, 8].In order to address this discrepancy, we extend Nusser’swork by including torque. We consistently keep track of aparticle’s angular momentum, allowing it to build up be-fore turnaround because of tidal interactions with neigh-boring protogalaxies [23] and to evolve after turnaroundbecause of nonlinear effects within the halo. Moreover,we compare the predictions of our halo model to simula-tion results.Self-similar secondary infall requires Ω m = 1 since anonvanishing Ω Λ introduces an additional scale. Apply-ing self-similarity to halo formation in the ΛCDM modeltherefore requires approximations and a mapping to ha-los in an Einstein de-Sitter universe. We assume thatthe linear power spectrum and background matter den-sity ¯ ρ m today are equal in both universes so that thestatistics, masses, and length scales of halos found in thetwo models are equivalent. Since the scale factor evolvesdifferently in both universes, the halo assembly historieswill differ.In section II, we define our self-similar system andtorqueing parameters. In section III, we set initial con-ditions and evolve the mass shells before turnaround. Insection IV, we describe evolution after turnaround andthen analyze the asymptotic behavior of the density pro-file at different scales in section V. In section VI, wegive numerical results, discuss the overall structure of thehalo, and compare to N-body simulations. We concludein section VII. II. SELF-SIMILAR DEFINITIONS
Here we explicitly define our self-similar system andderive constraints on the functional form of the mass dis-tribution within a halo and the angular momentum ofparticles in a particular shell.If the infall process is self-similar, then the halo’s ap-pearance does not change once all lengths are scaled tothe current turnaround radius. For our analysis, we de-fine the current turnaround radius as r ta ( t ) ≡ Ct β whereboth C and β are positive constants. The exponent β ,as we will find, depends on the initial perturbation spec-trum.The evolution of a particular mass shell must dependon time t and the shell’s turnaround time t ∗ . More explic-itly, assuming spherical symmetry, we have r = R ( t, t ∗ ).We define a self-similar system as one in which every tra-jectory obeys the following scaling: R (Λ t, Λ t ∗ ) = Λ β R ( t, t ∗ ) (1)where Λ is a constant. The above implies that the tra-jectory of one mass shell with turnaround time t canbe mapped to the trajectory of another mass shell withturnaround time t = Λ t . The exponent β follows since R ( t , t ) /R ( t , t ) = ( t /t ) β .Each shell of a self-similar system must also follow thesame equation of motion. From Newton’s law, the radialequation of motion for a mass shell with angular momen-tum is given by:¨ R ( t, t ∗ ) = − GM (cid:16) R ( t, t ∗ ) , t (cid:17) R ( t, t ∗ ) + L (cid:16) R ( t, t ∗ ) , t, t ∗ (cid:17) R ( t, t ∗ ) (2) where dots denote derivatives with respect to the firstargument, M is the mass of the halo interior to r and L is the angular momentum per unit mass of a particle inthe shell. Note that we enforce the mass to not dependexplicitly on t ∗ , while the angular momentum can. Aswe will show below, this is physically motivated. Fromeq. (1), we find:¨ R (Λ t, Λ t ∗ ) = Λ β − R ( t, t ∗ ) (3)Plugging in eqs. (1) and (3) into eq. (2) and simplifying,we find:¨ R (Λ t, Λ t ∗ ) = − Λ β − GM (cid:16) Λ − β R (Λ t, Λ t ∗ ) , t (cid:17) R (Λ t, Λ t ∗ )+ Λ β − L (cid:16) Λ − β R (Λ t, Λ t ∗ ) , t, t ∗ (cid:17) R (Λ t, Λ t ∗ ) (4)Changing variables from R ( t, t ∗ ) to R ( t, t ∗ ) /Ct β for themass and angular momentum and rewriting eq. (4), wefind:¨ R (Λ t, Λ t ∗ ) = − Λ β − GM (cid:16) R (Λ t, Λ t ∗ ) /C (Λ t ) β , t (cid:17) R (Λ t, Λ t ∗ )+ Λ β − L (cid:16) R (Λ t, Λ t ∗ ) /C (Λ t ) β , t, t ∗ (cid:17) R (Λ t, Λ t ∗ ) (5)Relabeling coordinates and enforcing consistency witheq. (2), we find the following constraints on the func-tional forms of the mass and angular momentum. M (cid:16) R ( t, t ∗ ) /Ct β , t (cid:17) = Λ β − M (cid:16) R ( t, t ∗ ) /Ct β , t/ Λ (cid:17) (6) L (cid:16) R ( t, t ∗ ) , t, t ∗ (cid:17) = Λ β − L (cid:16) R ( t, t ∗ ) /Ct β , t/ Λ , t ∗ / Λ (cid:17) (7)With the above in mind, we define the angular mo-mentum per unit mass L of a particle in a shell at r , andthe density ρ and mass M of the halo as follows. L ( r, t ) = B r ( t ) t f ( λ, t/t ∗ ) (8) ρ ( r, t ) = ρ B ( t ) D ( λ ) (9) M ( r, t ) = 4 π ρ B ( t ) r ( t ) M ( λ ) (10)where λ ≡ r/r ta ( t ) is the radius scaled to the currentturnaround radius and ρ B = 1 / πGt is the backgrounddensity for an Einstein de-Sitter (flat Ω m = 1) universe.Using eq. (1), it is straightforward to show: λ ( t, Λ t ∗ ) = λ ( t/ Λ , t ∗ ) (11)Eq. (11) implies that if one can compute λ ( t, t ∗ ) for aparticular mass shell t ∗ at all times, then one also knowsthe position of all other mass shells, labeled by Λ t ∗ withvarying Λ, at a particular time. This interpretation isvery powerful and will be used later in order to calculatethe mass profile after turnaround.If the mass profile M ( λ ) also depended explicitly on t ∗ , then the mass would not have to grow like the back-ground mass enclosed in the current turnaround radius.This is clearly not physical. Hence we suppressed the ex-plicit dependence on t ∗ . On the other hand, we’ve keptthe dependence on t ∗ in the angular momentum in orderto have this extra freedom. Inspired by tidal torque the-ory and numerical simulations, in eq. (8) we take f tobe: f ( λ, t/t ∗ ) = ( λ − γ if t < t ∗ , ( t/t ∗ ) ̟ +1 − β if t > t ∗ . (12)The constant B sets the amplitude of the angularmomentum at turnaround while γ ( ̟ ) controls howquickly the angular momentum increases before (after)turnaround. Constraints on B , γ , and ̟ will be dis-cussed in later sections.We’ve assumed that the halo is spherically symmetric.While simulated halos are triaxial [24], the descriptionabove is meant to represent an average halo. Since thereare no preferred directions in the universe, it should beexpected that a statistically averaged halo is sphericallysymmetric.In the above, L represents the angular momentum perunit mass of all particles in the shell. We impose thatall particles in the shell have orbital planes that are ran-domly distributed. This implies that the total vectorangular momentum of the mass shell, and hence the to-tal angular momentum of the halo J , vanishes. Hence,while individual particles on a mass shell gain angularmomentum in random directions throughout evolution,on average the mass shell remains spherical. Therefore,like we’ve assumed above, only one radial equation ofmotion is necessary to describe the evolution of the shell.Since our statistically averaged halo has a vanishingtotal angular momentum, this model cannot address thenonzero spin parameters observed in individually simu-lated halos [25, 26]. Nor can it reproduce the nonzerovalue of (cid:10) J (cid:11) expected from cosmological perturbationtheory [27–29]. However, R L dm where dm is the massof a shell, does not vanish for this model. We will use thisquantity, which is a measure of the tangential dispersionin the halo, to constrain our torque parameters. III. BEFORE TURNAROUND
The trajectory of the mass shell after turnaround de-termines the halo mass profile. In order to start inte-grating at turnaround, however, the enclosed mass of thehalo must be known. For the case of purely radial orbits,the enclosed mass at turnaround can be analytically cal-culated [4, 5]. For the case of orbits that have a timevarying angular momentum, we must numerically evolveboth the trajectory and M ( λ ) before turnaround in orderto determine the enclosed mass at turnaround.The trajectory of a mass shell follows from Newton’slaw. We have: d rdt = − GM ( r, t ) r + L ( r, t ) r (13)Rewriting eq. (13) in terms of λ and ξ ≡ log( t/t i ), where t i is the initial time, and plugging in eqs. (8), (10) and(12), we find: d λdξ + (2 β − dλdξ + β ( β − λ = − M ( λ ) λ + B λ − γ − (14)The angular momentum before turnaround was chosenso that eq. (14) does not explicitly depend on ξ . Thisallows for a cleaner perturbative analysis. Since r is anapproximate power law in t at early times, we still havethe freedom to choose a particular torque model inspiredby tidal torque theory. This will be discussed at the endof this section.In order to numerically solve the above equation, onemust know M ( λ ), a function we do not have a priori.Before turnaround, however, the enclosed mass of a par-ticular shell remains constant throughout evolution sinceno shells cross. Taking advantage of this, we relate dr/dt to M by taking a total derivative of eq. (10). We find: (cid:18) drdt (cid:19) M = β rt − (3 β − Ct β − MM ′ (15)In the above, a prime represents a derivative taken withrespect to λ . Taking another derivative of the above withrespect to time, plugging into eq. (13) and simplifying,we find an evolution equation for M : β ( β − λ + (3 β − β − MM ′ − (16)(3 β − M M ′′ ( M ′ ) = − M λ + B λ − γ − Given eqns. (14) and (16), we must now specify initialconditions when λ ≫
1. We assume the following pertur-bative solutions for D ( λ ) and λ ( ξ ) valid at early times. M ( λ ) follows from eq. (10). log H M (cid:144) M Ÿ L n FIG. 1. The variation of model parameter n with halo mass.Larger mass halos map to steeper initial density profiles. D ( λ ) = 1 + δ λ − n + δ λ − p + ... (17) M ( λ ) = λ (cid:18) δ − n λ − n + 3 δ − p λ − p + ... (cid:19) (18) λ ( ξ ) = λ e (2 / − β ) ξ (1 + λ e α ξ + λ e α ξ + ... ) (19)In the above, n characterizes the first order correctionto the background density. It is related to the FG pa-rameter ǫ through n = 3 ǫ . It is also related to the ef-fective power spectral index n eff = d ln P/d ln k through n = n eff + 3 [22]. Since n eff depends on scale and hencehalo mass (Appendix A), we have a relationship between n and halo mass. As Figure 1 shows, larger mass haloshave larger n . This is expected since larger smoothinglengths imply steeper initial density profiles. As in FG,we restrict 0 < n < n > . p characterizes the next order correction to the backgrounddensity caused by angular momentum. Consistency withour perturbative expansion (eqs. 17 through 19) demandsthat we take n < p < n . However, it is straightforwardto generalize to other cases (Section III A). As we showbelow, the constants { δ , λ , λ , α , α } are set by theequations of motion. The constants { δ , λ } are set byboundary conditions. Plugging eq. (18) into eq. (16) andenforcing equality between terms proportional to λ − n we find two possible solutions. β = 23 (cid:18) n (cid:19) (20) β = 23 (cid:18) − n (cid:19) (21) These represent the two solutions to the second orderdifferential equation (eq. 16). For the first case, theturnaround radius grows faster than the Hubble flowwhile for the second it grows slower. Hence, the firstsolution represents the growing mode of the perturba-tion while the second is the decaying mode. Since weare interested in the growth of halos, we will only con-sider the growing mode from now on. Next, imposing p = 2 γ + 4 and enforcing equality between terms propor-tional to λ − p in eq. (16), we find: δ = 9 n B ( p − p − n )(3 n + 2 p ) (22)Comparing eqs. (18) and (22), we see that the correc-tion to the initial mass caused by angular momentum isnegative since p > n . This is expected since the angularmomentum acts against gravity.Next, we find constraints for the parameters in eq. (19).Plugging in eq. (19) into eq. (14) and setting terms linearin δ , δ , λ , and λ equal to each other, we find: α = 23 (23) λ = δ n − λ − n (24) α = p (cid:18) β − (cid:19) (25) λ = δ p − λ − p (26)Eqns. (17) through (26) set the initial conditions foreqs (14) and (16). We evolve λ and M and choose { λ , δ } such that turnaround occurs ( dλ/dξ = − λβ )when λ = 1. The free parameters for evolution beforeturnaround are { n, B, p } . For the zero angular momen-tum case analyzed in FG, the enclosed mass is the sameregardless of n ; M (1) = (3 π/ . Including torque, how-ever, we find that the enclosed mass depends on the pa-rameters B and p . Figure 2 shows a contour plot of16 M (1) / π for n = 1. As expected, the enclosed massat turnaround must be larger than the no-torque casein order to overcome the additional angular momentumbarrier. B sets the amplitude of the angular momen-tum while p controls how the angular momentum growsin comparison to the mass perturbation. Larger B andsmaller p correspond to stronger torques on the massshell. Contour plots with different values of n give thesame features.The above perturbative scheme ensures that the inclu-sion of angular momentum preserves cosmological initialconditions. Analyzing eq. (19) for ξ → −∞ , and usingeqs. (23) and (25), we see that since p > n , the angularmomentum correction is subdominant to the density per-turbation correction. More importantly, at early times,the shell moves with the Hubble flow. Last, in orderto be consistent with cosmological initial conditions, the p B FIG. 2. Contour plot of 16 M (1) / π for n = 1 as a function oftorquing parameters B and p . Smaller p and larger B result inlarger torques on the mass shell. These larger torques requirebigger enclosed masses at turnaround, in order to counteractthe stronger angular momentum barrier. angular momentum of particles within a mass shell mustvanish at early times. Imposing r ∝ t / , and plugginginto eq. (8), we find that L ∝ t (1+ p/n ) / . Hence, for allvalues of p we consider, the angular momentum vanishesat early times and has a value of Br ∗ /t ∗ at turnaround.Note that shells which turn around later have larger an-gular momentum. A. Tidal Torque Theory
According to Hoyle’s tidal torque mechanism, massshells before turnaround gain angular momentumthrough their interactions with the tidal fields of neigh-boring protogalaxies [23]. Peebles claimed that the angu-lar momentum of protogalaxies in an Einstein de-Sitteruniverse grows as t / [27] while Doroshkevich showedthat for non-spherical regions, the angular momentumgrows as t [28, 29]. White confirmed Doroshkevich’s anal-ysis with N-body simulations [28]. Since the net angularmomentum of our model’s halo vanishes, we will insteaduse ˜ σ , defined below, to constrain p and B .˜ σ ≡ Z V L dm | ( r − r ) × ( v − v ) | (27)In the above, we integrate over the Lagrangian volume V L of the halo, r and v are the physical radius and velocityof particles within the halo, r is the center of mass of thehalo and v ≡ v ( r ). As described earlier in this section, λ ≫ (cid:10) ˜ σ (cid:11) using cosmologicallinear perturbation theory and compare to expectationsfrom our model.Using the Zel’dovich approximation [30], assuming aspherical Lagrangian volume with radius R , and workingto first order, we find: (cid:10) ˜ σ (cid:11) M = 6 a ˙ D M x A ( R ) (28)where M is the mass of the halo, a is the scale factor, D is the linear growth factor, dots denote derivativeswith respect to proper time, x max is the lagrangian radiusof the volume, R is the spherical top hat radius for ahalo of mass M and A ( R ) is a time independent functiondefined in Appendix B which has units of length. Notethat the scale factor a and ˙ D are the only quantitieswhich vary with time. For a matter dominated universe,( (cid:10) ˜ σ (cid:11) M ) / ∝ t , just as in the White analysis. This isexpected since the Lagrangian mass is time independent.Next we calculate eq. (27) from the perspective of ourmodel. Using eqns. (8), (10), and (12) and assuming firstorder corrections to M are negligible, we find:˜ σ = Z r max r min B r ta t λ − γ ∂M ( r, t ) ∂r dr = 4 π − γ B ρ B ( t ) r ta ( t ) t h λ − γ max ( t ) − λ − γ min ( t ) i (29)The lower limit of integration sets an effective smoothinglength which we choose to be λ ≫ p = 2 γ + 4and n < p < n , then for the range of n we consider,2 γ < r max . Equating eqs. (28)and (29) and assuming the first order corrections to r max in eq. (19) are negligible, we find p = 2 n and: B = 23 p − n ) M (1) ( n − / A ( R ) R (30)Eqs. (28) and (30) are derived in Appendix B. Note thatthe perturbative analysis, presented above, which is usedto calculate M (1) is not valid for p = 2 n . Redoing theanalysis for this special case, we find: α = 43 (31) δ = 914 B (2 n −
3) + (7 n − n − n − δ (32) λ = 914 λ − n (cid:18) B − δ ( n − (cid:19) (33) log H M (cid:144) M Ÿ L B FIG. 3. The variation of model parameter B with halo mass. B sets the angular momentum of particles at turnaround. For the remainder of this paper we impose p = 2 n , sothat angular momentum grows in accordance with cos-mological perturbation theory, and set n according tothe halo mass. Unfortunately, when comparing to N-body simulations (Section VI A), eq. (30) overestimatesthe angular momentum of particles at turnaround by afactor of 1.5 to 2.3. We discuss possible reasons for thisdiscrepancy in Appendix B. For convenience, B . ( B . )denotes B calculated using eq. (30) with the right handside divided by 1.5 (2.3). As described above, M (1) ineq. (30) depends on B . Therefore, in order to find B , wecalculate B and M (1) iteratively until eq. (30) is satis-fied.The relationship between n and n eff as well as tidaltorque theory implies that { n, p, B } are all set by thehalo mass. Figure 3 shows the variation of B with halomass. A ( R ) increases with halo mass since more power atlarge scales is included; R also increases with halo mass.These two competing effects cause a slight variation in B over seven orders of magnitude in halo mass. IV. AFTER TURNAROUND
Given the enclosed mass found at turnaround, we nowsolve for the trajectory and mass profile after turnaround.For convenience, we redefine the time variable to be ξ ≡ ln( t/t ta ), where t ta is the current turnaround time.The trajectory’s evolution equation after turnaround,with the appropriate torque model (eq. 12), is shown be-low. d λdξ +(2 β − dλdξ + β ( β − λ = − M ( λ ) λ + B λ e ̟ +1 − β ) ξ (34) The torque model after turnaround was chosen to notexplicitly depend on r , since r begins to oscillate ona much faster timescale than the growth of the halo.Nusser [14] and Sikivie et al. [17] focused on the case ̟ = 0. However, as was discussed before, this results indensity profiles steeper than what is predicted by numer-ical simulations.There are a number of dynamical processes that cancause a particle’s angular momentum to evolve afterturnaround. Dynamical friction [31] transfers the an-gular momentum of massive bound objects – like blackholes, globular clusters and merging satellite galaxies –to the background halo. A massive black hole at thecenter of the halo that dominates the potential at smallscales tends to make the velocity dispersion isotropic [32–34]. Bars [35, 36] and supermassive black hole binaries[37, 38] are also expected to perturb the dark mattervelocity distribution. While the torque model proposedafter turnaround is clearly very simplistic and may notaccurately describe some of the above phenomena, it stillallows us to get intuition for how torques acting on massshells change the structure of the halo.Analytically calculating ̟ is difficult since the halo af-ter turnaround is nonlinear. In Appendix C, we showin a simplistic manner how ̟ is sourced by substruc-ture and argue that dark matter dominated substructureshould cause steeper density profiles than baryon domi-nated substructure. In order to properly constrain ̟ , N-body simulations are required. This is beyond the scopeof this work.The initial conditions for eq. (34), enforced in theabove section, are λ ( ξ = 0) = 1 and dλ/dξ ( ξ = 0) = − β .As discussed before, self-similarity implies that all massshells follow the same trajectory λ ( ξ ). Hence, λ ( ξ ) caneither be interpreted as labeling the location of a partic-ular mass shell at different times, or labeling the locationof all mass shells at a particular time. We take advan-tage of the second interpretation in order to calculate themass profile.After turnaround, shells cross since dark matter is col-lisionless. Therefore, the mass interior to a particularshell does not stay constant. However, since λ ( ξ ) spec-ifies the location of all mass shells at a particular time,the mass interior to a given scale is simply the sum of allmass shells interior to it. The mass profile is then givenby [4, 5]: M ( λ ) = 2 n M (1) Z ∞ dξ exp[ − (2 /n ) ξ ] H [ λ − λ ( ξ )]= M (1) X i ( − i − exp[ − (2 /n ) ξ i ] (35)where M (1) is the normalization constant found in theprior section, H [ u ] is the heaviside function, and ξ i is the i th root that satisfies λ ( ξ ) = λ . The above is straightfor-ward to interpret. The roots ξ i label shell crossings at aparticular scale and the exponential factor accounts forthe mass difference between shells that turn around atdifferent times.Since the trajectory and the mass profile depend oneach other, it is necessary to first assume a mass profile,then calculate the trajectory from eq. (34) and resultingmass profile from eq. (35) and repeat until convergenceis reached.The density profile D ( λ ) is straightforward to deriveusing eqs. 9, 10 and 35. We find [4, 5]: D ( λ ) = 13 λ d M dλ = 23 n M (1) λ X i ( − i exp[ − (2 /n ) ξ i ] (cid:18) dλdξ (cid:19) − i (36) V. ASYMPTOTIC BEHAVIOR
Unlike N-body experiments, self-similar systems arenot limited by resolution. One can analytically infer theasymptotic slope of the mass profile close to the origin.FG did this by taking advantage of adiabatic invariance,self-consistently calculating the mass profile, and ana-lyzing the limit of the mass profile as λ →
0. Below,we generalize their analysis to the case of particles withchanging angular momentum. Unlike Nusser [14], we donot restrict our analysis to the case ̟ = 0.Just as in FG, we start by parameterizing the halomass and the variation of the apocenter distance r a . M ( r, t ) = κ ( t ) r α (37) r a r ∗ = (cid:18) tt ∗ (cid:19) q (38)In the above r ∗ is the turnaround radius of a mass shellwhich turns around at t ∗ . It is possible to relate q and α to n by taking advantage of adiabatic invariance. Theequation of motion for the mass shell is: d rdt = − Gκ ( t ) r α − + L ( t ) r (39)At late times, the orbital period is much smaller thanthe time scale for the mass and angular momentum togrow. Taking κ ( t ) and L ( t ) to stay roughly constant overan orbit and integrating the above equation, we find theenergy equation: (cid:18) drdt (cid:19) = 2 Gκ ( t ) α − r α − a − r α − ) − L ( t )( r − − r − a ) (40)The above relationship tells us how the pericenters r p evolve with time. Note that we only consider torquingmodels (eq. 12) which give rise to bound orbits. This restriction on ̟ will be discussed below. Defining y ≡ r p /r a and evaluating the above at r = r p , we find:1 − y α − y − − ≡ A ( y ) = ( α − L ( t )2 Gκ ( t ) r α +1 a ( t ) (41)For ̟ > ( < ) 0, the angular momentum of particles inthe mass shell increases (decreases). This gives rise topericenters that increase (decrease). Hence at late times,the orbit of a mass shell with increasing angular momen-tum will circularize and have y ∼
1, while the orbit ofa mass shell with decreasing angular momentum will be-come more radial, with y ≪
1. With this in mind, wecan now calculate the radial action in order to find how q relates to n and α . The radial action is given by: J = 2 Z r a r p dr (cid:18) drdt (cid:19) = 2 (cid:18) Gκ ( t ) α − (cid:19) / r ( α +1) / a × Z y ( t ) du (cid:2) (1 − u α − ) − A ( y )( u − − (cid:3) / (42)In the above we’ve assumed α >
1. Generalizing tothe case α < α = 1will be addressed later. For y ( t ) ≪
1, the above inte-gral is dominated by the region in which y ( t ) ≪ u ≪ κ ( t ) r α +1 a = const. For y ( t ) ∼
1, theorbit is circular, which implies the radial action vanishesand L ( t ) = Gκ ( t ) r α +1 a . Using eq. (38), and noting that κ ( t ) ∝ t s where s = 3 β − − αβ , we find at late times: q = ( α +1 { ̟ + n [ α (1 + n ) − } if ̟ ≥ n ( α +1) [ α (1 + n ) −
3] if ̟ < ̟ <
0, taking advantage of y ≪ y ( t, t ∗ ) = y ( t/t ∗ ) l , where: l = ( ̟ if α > , ̟/ ( α + 1) if α < . (44)and y r ∗ is the pericenter of a mass shell at turnaround.Constant angular momentum after turnaround corre-sponds to ̟ = 0. This case was addressed analyticallyin [14] and numerically in [17].We next take advantage of the functional form of themass profile. Following FG, we define P ( r/r a , y ) to bethe fraction of time a particle with apocenter distance r a and pericenter yr a , at a particular time t , spends inside r . P ( v, y ) = 0 ( v < y ) P ( v, y ) = I ( v, y ) I (1 , y ) ( y < v ≤ P ( v, y ) = 1 ( v >
1) (45)where I ( v, y ) ≡ (R vy du ((1 − u α − ) − A ( y )( u − − / if α > , R vy du (( u α − − A ( y )( u − − / if α < . (46)We see that the presence of pericenters causes the newcase v < y , which did not exist in the FG analysis. Selfconsistency demands that (cid:18) rr ta (cid:19) α = M ( r, t ) M ( r ta , t ) = Z M ta dM ∗ M ta P (cid:18) rr a ( t, t ∗ ) , y ( t, t ∗ ) (cid:19) (47)where M ∗ is the mass internal to a shell that turns aroundat t ∗ and M ta is the current turnaround mass. The in-tegral assigns a weight to each shell depending on howoften that shell is below the scale r . Noting from eq. (10)that M ∗ = M ta (cid:18) t ∗ t (cid:19) β − , (48)using eq. (38) and transforming integration variables, wefind: (cid:18) rr ta (cid:19) α − k = k Z ∞ r/r ta duu k P ( u, y ( t, t ∗ )) (49)where k = 62 + n (2 − q ) (50)As u increases, the above integral sums over shells withsmaller t ∗ . Since the pericenter of a shell evolves withtime, the second argument of P depends on u . The de-pendence, as we showed, varies with torque model (signof ̟ ); hence we’ve kept the dependence on u implicit.Next we analyze the above for certain regimes of r/r ta ,and certain torquing models, in order to constrain therelationship between α and k . A. r/r ta ≪ y , ̟ < For ̟ <
0, particles lose angular momentum over time.When probing scales r/r ta ≪ y , mass shells with t ∗ ≪ t ta only contribute. As a result, y ( t, t ∗ ) ≪
1. Using eq.(44), we then find: y ( t, t ∗ ) = y (cid:18) tt ∗ (cid:19) l = y (cid:18) rur ta (cid:19) δ (51)where δ ≡ l/ ( q − β ). For bound mass shells, q − β < δ >
0, the first argument of P in eq. (49)increases while the second decreases as we sum over shellsthat have turned around at earlier and earlier times ( u →∞ ). For r/r ta ≪ y , mass shells which most recentlyturned around do not contribute to the mass inside r/r ta since we are probing scales below their pericenters. Massshells only begin to contribute when the two argument of P are roughly equal to each other. This occurs around: u = y ≡ (cid:0) y ( r/r ta ) δ (cid:1) / (1+ δ ) (52)Hence, we can replace the lower limit of integration ineq. (49) with y . We next want to calculate the behaviorof eq. (49) close to y in order to determine whether theintegrand is dominated by mass shells around y or massshells that have turned around at much earlier times.The first step is to calculate the behavior of P ( u, y ) for u ≈ y . We find: P ( u, y ) ∝ u / (1 − y/u ) / × ( y / if α > ,y − α/ if α < . (53)Given the above, we evaluate the indefinite integral ineq. (49), noting that y is a function of u (eq. 51). For u ∼ y , we find: Z duu k P u, y (cid:18) rur ta (cid:19) δ ! ∝ ( u/y − / ( y − k if α > ,y / − k − α/ if α < . (54)Now comes the heart of the argument. Following thelogic in FG, if we keep u/y fixed and the integrand blowsup as r/r ta →
0, then the left hand side of eq. (49) mustdiverge in the same way as the right hand side shown ineq. (54). Therefore, using eq. (52): α − k = ( δ (1 − k ) / (1 + δ ) if α > ,δ (3 / − k − α/ δ ) if α < . (55)Otherwise, if the right hand side converges, then the inte-grand will not depend on r/r ta as r/r ta →
0. Therefore,the left hand side cannot depend on r/r ta either, whichimplies α = k . Solving the above system of equationsfor α given eqs. (43) and (50) and making sure the solu-tion is consistent, (ie: using eq. (55) only if the integranddiverges), we find:For n ≤ α = 1 + n − p (1 + n ) + 9 n̟ ( n̟ − n̟k = 1 + n + 3 n̟ − p (1 + n ) + 9 n̟ ( n̟ − n̟ (4 + n ) q = 1 + n − n̟ − p (1 + n ) + 9 n̟ ( n̟ − n For n ≥ α = k = 31 + n , q = 0 (56)The above solutions are continuous at n = 2. Moreover,taking the no-torque limit ( ̟ →
0) for n ≤ n ≥
2, which is consistent with analyticand numeric results from [14, 17]. Taking the limit, ̟ →−∞ reproduces the FG solution, as expected, since theshell loses its angular momentum instantly. The solutionfor n ≥ ̟ . This is because the mass isdominated by shells with turnaround time t ∗ ≪ t ta whichhave effectively no angular momentum. In other words,for ̟ <
0, the solution should only depend on torquingparameters when the mass is dominated by shells thathave turned around recently. B. r/r ta ≪ y , ̟ > For ̟ >
0, the angular momentum of particles increasewith time. As mentioned above, when probing scales r/r ta ≪ y , mass shells with t ∗ ≪ t ta only contribute.As a result, y ( t, t ∗ ) ∼
1. In other words, the orbits areroughly circular. We can therefore replace the lower limitof integration in eq. (49) with 1 since mass shells will onlystart contributing to the sum when u ∼ y ∼
1. Hence,the right hand side of eq. (49) does not depend on r/r ta ,which implies α = k . Using eq. (43) and (50), we find: α = k = 31 + n − n̟ , q = 2 ̟ , for 0 ≤ n ≤ ̟ = 0, is consistent with theanalysis in the prior subsection. The singularity ̟ =(1 + n ) / n implies q = β . This physically corresponds tothe orbital radius of mass shells increasing at the samerate as the turnaround radius and results in orbits thatare not bound and a cored profile where there are noparticles internal to a particular radius. This breaks theassumption of a power law mass profile (eq. 37); hencewe only consider ̟ < (1 + n ) / n .Unlike the Nusser solution, certain parameters give α >
3, which corresponds to a density profile which con-verges as r →
0. Since the angular momentum acts like aheat source dρ/dr > C. y ≪ r/r ta ≪ In this regime, we are probing scales larger than thepericenters of the most recently turned around massshells. As a result, P ( u, y ) is dominated by the contribu-tion from the integrand when u ≫ y . Therefore: P ( u, y ) ∝ ( u if α > ,u (3 − α ) / if α < . (58)Hence the integral in eq. (49) becomes: Z duu k P ( u, y ( t, t ∗ )) ∝ ( u − k if α > ,u / − k − α/ if α < . (59)Following the logic in the prior section, if the integraldiverges as r/r ta →
0, then we set the exponents on theleft hand side and right hand side equal to each otherso that both sides diverge in the same way. If the inte-gral converges, then the left hand side cannot depend on r/r ta , which implies α = k . Given these arguments, andimposing consistency with the above inequalities on α tofind the appropriate ranges for n , we find: α = 1 , k = 64 + n , q = n − n , for n ≤ α = k = 31 + n , q = 0 , for n ≥ α = 1 and yet, for certain partsof parameters space, eqs. (56), (57) and (60) give α = 1.However, since the solutions are continuous as α → α = 1as well. VI. STRUCTURE OF THE HALO
In this section, we discuss the radial structure of galac-tic size halos and compare directly to numerical N-bodysimulations. Note however, that the mass of a halo is notwell defined when our model is applied to cosmologicalstructure formation since it is unclear how the sphericaltop hat mass which characterizes the halo when it is lin-ear relates to the virial mass which characterizes the halowhen it is nonlinear. For halos today with galactic sizevirial masses, we assume the model parameter n whichcharacterizes the initial density field, is set by a sphericaltop hat mass of 10 M ⊙ . As described in prior sections,0specifying the top hat mass also sets model parameters B and p . Before comparing directly to N-body simulations,we first describe how ̟ influences the halo.Figure 4 shows the mass M ( λ ) and density profiles D ( λ ) for galactic size halos { n = 0 . , p = 2 n, B . =0 . } with varying ̟ . The spikes in the density pro-file are caustics which form at the shell’s turning points.They form because of unphysical initial conditions; we as-sume each shell has zero radial velocity dispersion. Thestructure of the halos naturally break down into threedifferent regions. The dividing points between these re-gions are roughly the virial radius ( r v ) and y r ta , thepericenter of the mass shell which most recently turnedaround.As in dark matter N-body simulations, we associatethe virial radius with r , the radius at which M ( r ) =800 πρ B r / λ ∼ . r > r v , the mass profile flattens and then starts to in-crease. The flattening is equivalent to what is seen onlarge scales in N-body simulations where ρ ∝ r − . Themass profile then starts to increase again because at largeradii, where λ ≫
1, the density is roughly constant, whichimplies
M ∝ λ . For r ∼ r v , it is difficult to make an-alytic predictions for the mass profile because adiabaticinvariance breaks down. In other words, the mass of thehalo and angular momentum of a shell change on thesame time scale as the shell’s orbital period.As discussed in the prior section, for y r ta ≪ r ≪ r v ,we can take advantage of adiabatic invariance to infer thelogarithmic slope of the mass profile. Since this regimeprobes a scale much larger than the pericenters of themass shells, the angular momentum does not affect thedynamics and we recover the FG solution. For our par-ticular choice of n = 0 .
77, this gives an isothermal profilewith ρ ∝ r − . However, since n < r/r ta ≪ y , angular momentum begins toplay a role and the halo starts to exhibit different featuresthan the FG solution. The behavior is very intuitive. Themass of a particular shell does not contribute to the in-ternal mass when probing radii less than the pericenterof that mass shell. Therefore, as one probes radii smallerthan the pericenter of the most recently turned aroundmass shell, one expects a steeper fall off than the FG so-lution, since less mass is enclosed interior to that radius.Moreover, varying ̟ varies the pericenter of mass shellsover time. Increasing (decreasing) angular momentum, ̟ > ( < ) 0, causes the pericenters to increase (decrease)over time. This results in profiles which are steeper (shal-lower) than the no-torque case.It is informative to find how the transition radius y depends on model parameters. Since the mass and angu-lar momentum grow significantly before the first pericen-ter, we can only approximately determine this relation-ship. We assume that the profile is isothermal on large x -6 x -5 x -4 ϖ = − 0.4 ϖ = − 0.2 ϖ = 0 ϖ = 0.2 ϖ = 0.4 λ M ( λ ) rv ∕ rtay λ x x x x D ( λ ) ϖ = − 0.4 ϖ = − 0.2 ϖ = 0 ϖ = 0.2 ϖ = 0.4 y rv ∕ rta FIG. 4. The mass and density profiles for galactic size halos { n = 0 . , p = 2 n, B . = 0 . } with varying ̟ . The valueof ̟ changes how pericenters evolve with time and therebyaffects how many shells at a particular scale contribute tothe internal mass. The above numerically computed profilesmatch analytic predictions. The virial radius ( r v ) and firstpericenter passage ( y r ta ) are labeled for clarity. scales and the halo mass and shell angular momentumare fixed to their turnaround values. For y ≪
1, thetransitional radius y solves the transcendental equation y ln( y ) = − B / M (1). As expected, B → y →
0. As shown in Figure 2, M (1) varies with { B, p } . However, for reasonable parameter values, themass normalization is changed at most by a factor of 2.Therefore, y most strongly depends on B , and p has anegligible effect on the structure of the halo. As seen infigure 4, y should also depends on ̟ . The above ap-proximation neglects this dependence since we assumedthe angular momentum is set to the turnaround value.Figure 5 shows the radius of a mass shell, normal-1 t (cid:144) t * r (cid:144) r * t (cid:144) t * r (cid:144) r * t (cid:144) t * r (cid:144) r * FIG. 5. The radius of a mass shell, normalized to itsturnaround radius, for a galactic size halo { n = 0 . , p =2 n, B . = 0 . } , plotted vs time. The top panel shows ashell with particles that have decreasing angular momentum( ̟ = − . ̟ = 0). The bottompanel shows a shell with particles that have increasing angularmomentum ( ̟ = 0 . - - Λ t * r * d r d t - - Λ t * r * d r d t - - Λ t * r * d r d t FIG. 6. A phase space diagram of a galactic size halo { n =0 . , p = 2 n, B . = 0 . } at the current turnaround time.Velocities are normalized to the turnaround time and radiusof each shell. The top panel shows a halo with ̟ = − .
2. Themiddle panel shows a halo with ̟ = 0. The bottom panelshows a halo with ̟ = 0 . r ∗ , for a galactic size halo { n = 0 . , p = 2 n, B . = 0 . } , as a function of time.In the top panel, particles in the shell lose angular mo-mentum ( ̟ = − . ̟ = 0), while in the bot-tom panel particles in the shell gains angular momentum( ̟ = 0 . ̟ . The period of oscillation is set by the shell’sapocenter r a and the mass internal to r a . Using the adi-abatic invariance relations we found in Section V andassuming Kepler’s third law, we find that the period ofthe orbit P ∝ r a for shells with decreasing angular mo-mentum and P ∝ r a for shells with increasing angularmomentum. Moreover, from eqs. (56) and (60), r a de-creases (increases) with time for ̟ < ( > ) 0. ThoughKepler’s third law doesn’t hold for this system, it stillgives intuition for the above results.Figure 6 shows the phase space diagram for a galacticsize halo { n = 0 . , p = 2 n, B . = 0 . } . In the toppanel, the particles in the shell lose angular momentum( ̟ = − . ̟ = 0), while in the bottom panel theparticles in the shell gain angular momentum ( ̟ = 0 . t ∗ and radius r ∗ .Unlike FG, the presence of angular momentum results incaustics associated with pericenters as well, which can beseen in the lower panel of Figure 4. In addition, since anincreasing angular momentum results in increasing peri-centers, the pericenter caustics are more closely spaced inthe lower panel than in the upper panel. Moreover, theamplitude of the radial velocity is smaller in the lowerpanel because orbits are circularizing. The phase spacecurve appears to intersect itself because we did not plotthe tangential velocity component. In full generality, thedistribution in the phase space ( r, v r , v t = L/r ) for ourmodel is a non-self-intersecting one-dimensional curve.
A. Comparing with N-body Simulations
In this subsection, we compare the density profile ofour model’s halo to empirical fits inspired by N-bodysimulations. We first numerically calculate the densityprofile for a galactic size halo with ̟ = 0 .
12. Thisvalue of ̟ was chosen so that ρ ∝ r − on small scales.We then compute the spherically averaged density in 50spherical shells equally spaced in log r over the range1 . × − < r/r v <
3, and take r v = r (definedabove). This is the same procedure followed with therecent Aquarius simulation [6]. Next, we calculate r − ,the radius where r ρ reaches a maximum. For our halo, r ∕ r −2 ρ r [ ρ − r − ] αE = 0.159B r ∕ r −2 ρ r [ ρ − r − ] αE = 0.159B FIG. 7. Spherically averaged density profile for the secondaryinfall model compared with NFW and Einasto profiles. Thesecondary infall model is calculated for a galactic size halo { n = 0 . , p = 2 n } with ̟ = . B . = 0 .
39 while in the bottom panel we chose B . = 0 . as discussed above, the profile is isothermal over a rangeof r . Moreover, the maximum peaks associated with thecaustics are unphysical. So, we choose a value of r − in the isothermal regime that gives good agreement withthe empirical fits. Changing r − does not change ourinterpretation of the results.In Figure 7 we compare our spherically averaged den-sity profile to NFW and Einasto profiles. We plot r ρ inorder to highlight differences. The NFW profile is givenby [39]: ρ ( r ) = 4 ρ − ( r/r − )(1 + r/r − ) (61)while the Einasto profile is given by:3ln h ρ ( r ) /ρ − i = ( − /α E )[( r/r − ) α E −
1] (62)where ρ − is the density of our halo at r − and α E , knownas the shape parameter, sets the width of the r ρ peak.In the top panel, we use B . = .
39 while in the bottompanel, we use B . = .
26. We choose α E = 0 .
159 sincethis value was used in Figure 3 of Navarro et. al. [6].We see that the secondary infall model works surpris-ingly well. The peaks are a result of the caustics thatarise because of cold radial initial conditions. The firstspike on the right comes from the first apocenter passagewhile the second comes from the first pericenter passage.The location of pericenter is most strongly influenced bythe model parameter B . Hence, the isothermal region issmaller in the top panel than in the bottom panel sinceparticles have less angular momentum at turnaround inthe lower panel than in the upper panel. The parameter B then plays the same role as α E ; it sets the width ofthe isothermal region. If we assume N-body simulationsfaithfully represent dark matter halos, then Figure 7 im-plies that our estimate of B in eq. (30) overestimates theactual value by 1.5 to 2.3. We discuss possible reasonsfor this in Appendix B. VII. DISCUSSION
N-body simulations reveal a wealth of informationabout dark matter halos. Older simulations predict den-sity profiles that are well approximated by an NFW pro-file [39], while more recent simulations find density pro-files that fit better with a modified NFW profile [40] orthe Einasto profile [6]. In an attempt to gain intuitionfor these empirical profiles, we’ve generalized the self-similar secondary infall model to include torque. Thismodel doesn’t suffer from resolution limits and is muchless computationally expensive than a full N-body simu-lation. Moreover, it is analytically tractable. Using thismodel, we were able to analytically calculate the densityprofile for r/r ta ≪ y and y ≪ r/r ta ≪
1. Note thatthe self-similar framework we’ve extended predicts powerlaw mass profiles on small scales. Hence, it is inconsistentwith an Einasto profile.It is clear from our analysis that angular momentumplays an essential role in determining the structure ofthe halo in two important ways. First, the amount ofangular momentum at turnaround ( B ) sets the width ofthe isothermal region. Second, the presence of pericenterssoftens the inner density slope relative to the FG solutionbecause less mass shells contribute to the enclosed mass.Moreover, the interior density profile is sensitive to theway in which particles are torqued after turnaround ( ̟ ).If we assume that ̟ is constant for all halos, then thissecondary infall model predicts steeper interior densityprofiles for larger mass halos. More specifically, if we usethe value of ̟ = 0 .
12 which gave ρ ∝ r − for galactic size halos, then ρ ∝ r − . for a 10 M ⊙ halo and ρ ∝ r − . for a 10 M ⊙ halo. This trend towards steeper interiorslopes for larger mass halos, and hence non-universality,has been noticed in recent numerical simulations [41, 42]as well as more general secondary infall models [43]. Onthe other hand, if we assume that all halos have ρ ∝ r − as r →
0, then ̟ must vary with halo mass. More specif-ically, halos with mass M < M ⊙ must have particleswhich lose angular momentum over time ( ̟ <
0) whilehalos with mass
M > M ⊙ must have particles whichgain angular momentum over time ( ̟ > ̟ must conspire to erase anydependence on initial conditions. A more thorough treat-ment requires the use of N-body simulations, which isbeyond the scope of this paper.It is also possible to predict a dark matter halo’s den-sity distribution if one assumes a mapping between amass shell’s initial radius, when the structure is linear,to its final average radius, which is some fraction of itsturnaround radius. From this scheme, one can also infera velocity dispersion, using the virial theorem. Unfortu-nately, this scheme does not give any information aboutthe halo’s velocity anisotropy. Our self-similar prescrip-tion discussed above, on the other hand, contains all ofthe velocity information. Hence, one can reconstruct thevelocity anisotropy profile given the trajectory of a massshell. The velocity anisotropy is significant since it de-scribes to what degree orbits are radial. Moreover, it canbreak degeneracies between n and ̟ in our halo model.We will discuss the velocity structure of our halo model,including the pseudo-phase-space density profile, in moredetail in Paper 2 of this series. There we will once againcompare our halo predictions to the recent Aquarius sim-ulation results [6].While the above self-similar prescription has its clearadvantages, it’s also unphysical since mass shells atturnaround are radially cold. The same tidal torquemechanisms which cause a tangential velocity dispersion[23], should also give rise to a radial velocity dispersion.For a more physical model, one would need to imposeself-similarity to a phase space description of the haloand include sources of torque as diffusion terms in theBoltzmann equation. This will be the subject of Paper 3of this series.As we’ve shown, the way in which particles are torquedafter turnaround ( ̟ ) influences the interior power law ofthe density profile. One way to source this change inangular momentum is through substructure that is as-pherically distributed throughout the halo. It is reason-able to assume that substructure dominated by baryonstorque halo particles more strongly than substructuredominated by dark matter since baryons can achievehigher densities and hence are not tidally disrupted aseasily. If this is the case, then torques sourced by baryonswould result in a larger value of ̟ than torques sourcedby dark matter. According to the predictions of this sec-ondary infall model, this would lead to less cuspy profiles4(See Appendix C for a more detailed discussion). There-fore a more thorough understanding of ̟ coupled withthis simplified model of halo formation could potentiallyshed light on the Cusp Core problem and thereby possi-bly bridge the gap between simulations and observations. ACKNOWLEDGMENTS
The authors acknowledge support from NASA grantNNG06GG99G and thank Robyn Sanderson and PaulSchechter for many useful conversations.
Appendix A: Calculating n eff The effective primordial power spectral index, n eff , re-lates our model parameter n to the halo mass M . Theeffective index is defined by: n eff ≡ − d ln σ R d ln R − σ R ≡ Z d k (2 π ) P ( k ) W R ( k ) (A2)and W R ( k ) ≡ kR ) − ( kR ) cos( kR )( kR ) . (A3)The scale R is set by the top hat mass of the halo ( M =4 πρ m R / ρ m is the dark matter backgrounddensity today. The power spectrum today, P ( k ), is givenby [44]: P ( k ) = (cid:18) k H Ω m (cid:19) P R ( k ) T ( k ) D ( a = 1) (A4)where k P R ( k )2 π = △ R ( k ) (cid:18) kk (cid:19) n s − , (A5) D is the linear growth factor normalized so that D/a → a → a is the scale factor, T ( k ) is the transferfunction and we choose cosmological parameters derivedfrom WMAP7: △ R ( k ) = 2 . × − , h = 0 . m =0 . n s = 0 . k = 0 .
002 Mpc − [45]. Wecalculate T ( k ) using CMBFAST [46]. Appendix B: Tidal Torque Theory
We first derive eq. (28) using cosmological linear per-turbation theory. Starting with the Zel’dovich approxi-mation [30], we have: r ( q , t ) = a ( t ) (cid:16) q − D ( t ) ∇ φ ( q ) (cid:17) (B1)where r is the physical radius, q is a Lagrangian coordi-nate, and φ , which is time independent, is related to theNewtonian potential Φ through the following: φ = 14 πGρ m Da Φ (B2)where ρ m is the dark matter background density. Thevelocity, ∂ r /∂t , then is given by: v ( q , t ) = H ( t ) r ( q , t ) − a ( t ) ˙ D ( t ) ∇ φ ( q ) (B3)where H ≡ d ln a/dt and dots denote derivatives withrespect to time. Therefore, to first order in φ :( r − r ) × ( v − v ) = − a ˙ D ( q − q ) × ( ∇ φ − ∇ φ )(B4)where ∇ φ ≡ ∇ φ ( q ). Plugging this expression into eq.(27), taking an expectation value, defining x ≡ q − q ,rewriting in terms of the velocity perturbation variable Ψ ≡ − D ∇ φ and using index notation, we find: < ˜ σ > = a ˙ D D ǫ ijk ǫ ilm Z V L d xρ m a x j x l × < (cid:16) Ψ k (0) − Ψ k ( x ) (cid:17)(cid:16) Ψ m (0) − Ψ m ( x ) (cid:17) > C (B5)where ǫ ijk is the Levi-Civita tensor and V L is the La-grangian volume of the halo. We assume V L is sphericalwith radius x max . Note that to linear order in φ , theLagrangian coordinate x is equivalent to a comoving co-ordinate. The subscript C denotes an expectation valuetaken over constrained gaussian fields since we want toaverage only over peaks in the density field.Using the formalism developed in [47], we choose, forsimplicity, to use the zeroth and first order derivatives ofthe smoothed density field to constrain the velocity per-turbation. Our treatment and conventions are identicalto that used in Appendix A of Ma and Bertschinger [48].For more details, please refer to this reference. We find:5 h Ψ i ( x )Ψ j ( x ) i C = (cid:16) σ − ¯ η ( x ) σ (cid:17) ( δ ij − ˆ x i ˆ x j ) + σ − x ¯ η ( x ) σ − (cid:0) ¯ ξ ( x ) − η ( x ) (cid:1) σ ! ˆ x i ˆ x j (B6) h Ψ i ( x )Ψ j (0) i C = (cid:16) γ ( x ) x − ¯ η (0)¯ η ( x ) σ (cid:17) δ ij + (cid:16) dγ ( x ) dx − γ ( x ) x − x ¯ η (0) σ d ¯ η ( x ) dx (cid:17) ˆ x i ˆ x j (B7) h Ψ i (0)Ψ j (0) i C = (cid:16) σ − ¯ η (0) σ (cid:17) δ ij (B8)where: σ ≡ Z d k (2 π ) k − P ( k ) (B9) σ ≡ Z d k (2 π ) P ( k ) W R ( k ) (B10) σ ≡ Z d k (2 π ) k P ( k ) W R ( k ) (B11)¯ η ( x ) ≡ Z d k (2 π ) P ( k ) W R ( k ) j ( kx ) kx (B12)¯ ξ ( x ) ≡ Z d k (2 π ) P ( k ) W R ( k ) j ( kx ) (B13) γ ( x ) ≡ Z d k (2 π ) P ( k ) k − j ( kx ) . (B14) W R ( k ) and P ( k ) are defined in eq. (A3) and (A4) re-spectively and the spherical Bessel functions are j ( x ) = x − sin x and j ( x ) = x − (sin x − x cos x ). Eq. (B7) cor-rects an error in equation (A15) of Ref. [48].Notice from eqs. (B6) through (B8) that the expres-sions separate into terms proportional to δ ij and termsproportional to ˆ x i ˆ x j . The terms proportional to ˆ x i ˆ x j vanish in eq. (B5) because of the antisymmetry of theLevi-Civita tensors. Eq. (B5) then reduces to: (cid:10) ˜ σ (cid:11) M = 2 a ˙ D D Z V L d xρ m a x f ( x, R ) (B15)where f ( x, R ) = 2 σ − γ ( x ) x − ¯ η ( x ) σ + 2¯ η (0)¯ η ( x ) σ − ¯ η (0) σ (B16)Last, defining u = x/x max , we find: (cid:10) ˜ σ (cid:11) M = 6 a ˙ D D M x max Z u f ( ux max , R ) du ≡ a ˙ D M x A ( R ) (B17) In the above A has units of Mpc and x max = R since thescale of the galaxy today is equivalent to its lagrangiansize to linear order in perturbation theory. Note that f /D is time independent.Now, we derive eq. (30). Consistency with the sec-ondary infall model demands that we assume Ω m = 1.Equating the time dependence of eq. (28) and eq. (29)for an Einstein de-Sitter universe ( D = a ) at early times( r max ∝ t / ), we find p = 2 n where p = 2 γ + 4. Giventhis relationship, we now equate eq. (28) to eq. (29) andsolve for B . We find: B = 23 p − n ) a (cid:18) r max r ta (cid:19) n − A ( R ) r ta (B18)where we’ve used eqs. (10) and (18) evaluated at earlytimes to cancel the mass as well as the relationship r max = ax max . Now we evaluate eq. (B18) at early timessince tidal torque theory only applies when the halo islinear. To relate r max at some initial time ( t i ) to r ta today ( t ), we use the conservation of mass.4 π ρ B ( t i ) r ( t i ) = 4 π M (1) ρ B ( t ) r ta ( t ) (B19)Evaluating eq. (B18) at t i with the use of eq. (B19) andnoting that r ta ∝ t β , we find: B = 23 p − n ) M (1) ( n − / A ( R ) r ta ( t ) (B20)As expected, the time dependence of B vanishes. Last,assuming r ta ( t ) = R , we reproduce eq. (30). The quanti-ties n, R, A ( R ) are calculated in an Ω m = 1 universe withthe same background matter density and power spectrumof ΛCDM today. This ensures that the statistics, mass,and size of halos in both universes are equivalent today.As mentioned previously, eq. (30) overestimates the an-gular momentum of particles at turnaround by a factor of1.5 to 2.3. One potential source of error is to assume thelagrangian volume is spherical. Assuming an ellipsoidallagrangian volume with axis ratios 1 : a : b , we find that B is at most reduced by 8% when 0 . < a, b <
1. Thediscrepancy may be caused by not including higher or-der constraints on the smoothed density field. However,comparing B when calculated using zeroth and first orderderivative constraints with B when calculated withoutusing constraints results in only percent level differences;so it seems unlikely that constraints would have a signif-icant effect.N-body simulations use friends-of-friends group findersin order to identify halos and subhalos [49, 50]. Thisalgorithm, however, removes particles that are groupedto neighboring halos and hence neglects a contributionto ˜ σ . Trying to mimick this selection effect, we replacedeq. (B15) with6 (cid:10) ˜ σ (cid:11) M = 2 a ˙ D D Z ∞ d xe − x / R ρ m a x f ( x, R )(B21)where ¯ R ≡ (2 / π ) / R ensures that the mass enclosedwithin the lagrangian volume is equal to the mass of thehalo calculated by the simulation. However, calculating B in this manner leads to overestimating the angularmomentum at turnaround by ∼
4, as opposed to ∼ f ( x, R ) is an increasing function of x ( f ∼ x . near R for 10 M ⊙ halos), most of the con-tribution to A ( R ) comes from close to R . Therefore,while the gaussian cutoff decreases the contribution to B around R , it includes contributions beyond R , leadingto a worse estimate. This highlights that B significantlydepends on the outer parts of the halo. Hence, overesti-mating B by ∼ B is set during the lin-ear regime. Assuming that the shell is dominated bysubstructure at turnaround, nonlinear interactions likedynamical friction and tidal stripping play an importantrole from the time of turnaround to the first pericenterpassage [51, 52]. As time goes on, these effects becomeless important since substructure in the shell becomessubdominant. Including these extra interactions shouldlead to smaller estimates of B at first pericenter passageand hence potentially explain our overestimate, but isbeyond the scope of this work. Appendix C: Evolution After Turnaround
In this Appendix, we use dimensional analysis in orderto gain intuition about ̟ , a parameter that describestidal torque after turnaround. First, consider the timederivative of L , where L is the angular momentum perunit mass of a particle with radius r at time t . dL dt = 2 (cid:16) r ( v · a ) − ( r · a )( r · v ) (cid:17) (C1)In the above, a ( v ) is the acceleration (velocity) of theparticle. We now decompose the acceleration vector intoa radial (ˆ r ) and tangential (ˆ t ) component and use thisbasis to rewrite the velocity vector. a = a r ˆ r + a t ˆ t (C2) v = v r ˆ r + v t ˆ t + v p ˆ p (C3)The direction ˆ p is orthogonal to both ˆ r and ˆ t . Notethat all basis vectors depend on position. Plugging inthe above decomposed vectors into eq. (C1), we find: dL dt = 2 r v t a t . (C4) As expected, changes in L are sourced by deviationsfrom spherical symmetry that create nonzero a t . Now,imagine a spherically symmetric halo, roughly describedby our self-similar infall model with ̟ >
0, with a clumpof mass m in the shell at radius r . We assume that m is small enough so that it does not influence the radialequation of motion of the shell at r . We focus on ̟ > M ⊙ halo to be consistent with the NFW profile(Section VI A).Next, consider averaging dL /dt over an orbital periodand over a spherical shell of radius r , in order to comparethe change in angular momentum sourced by the clumpat r to our model’s prescription for angular momentumevolution. For ̟ >
0, orbits are roughly circular at latetimes. Hence, we assume averaging over an orbital pe-riod is equivalent to evaluating the right hand side of eq.(C4) at roughly the apocenter radius of the shell. Asdescribed in Section II, the orbital planes of all parti-cles in a shell at r are randomly aligned. Therefore weexpect v t averaged over a sphere to vanish. However, ifthere exists an excess mass m , then all the particles willbe pulled slightly in that direction, leading to a nonzeroaverage. We therefore assume v t ∝ P a t where P , theorbital period of the particle, is taken to be a dynamicaltime ( P ∝ ρ ( r ) − / ). Last assuming r ≪ r , we find: (cid:28) dL dt (cid:29) ∝ r m r ρ / . (C5)In order for the secondary infall halo model to be con-sistent, the right hand and left hand side must have thesame time scaling. Assuming m ∝ t µ , associating r and r with their respective apocenters, and noting from eq.(57) that r ∝ t ̟ , r ∝ t ̟ , and ρ ∝ t − ̟ , we find: ̟ = 13 (1 + 2 µ ) (C6)The same scaling relationship holds for r ≫ r . Henceone could imagine substructure in the shell at r sourc-ing a change in L of the shell at r and substructurein the shell at r sourcing a change in L of the shell at r . Therefore, a hierarchy of substructure non-sphericallydistributed, which is subdominant to the monopole con-tribution of the halo, would result in a halo roughly con-sistent with our described secondary infall model.Eq. (C6), which is only valid for µ > − / ̟ >
0, together with eq. (57) relates the steepnessof the inner density profile to the mass loss rate of sub-structure. If the clump does not lose mass ( µ = 0), then ̟ = 1 / ρ ∝ r . If the clump loses mass ( µ < should have steeper density profiles than galaxies7which include baryons. This is expected since baryonsstir particles around more efficiently, causing larger peri- centers and less dense interiors. A more thorough treat-ment that involves constraining ̟ with simulations isbeyond the scope of this paper. [1] J. E. Gunn and J. R. Gott, III, Astrophys. J. , 1(Aug. 1972).[2] J. R. Gott, III, Astrophys. J. , 296 (Oct. 1975).[3] J. E. Gunn, Astrophys. J. , 592 (Dec. 1977).[4] E. Bertschinger, Astrophys. J. , 39 (May 1985).[5] J. A. Fillmore and P. Goldreich, Astrophys. J. , 1(Jun. 1984).[6] J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vo-gelsberger, S. D. M. White, A. Jenkins, C. S. Frenk,and A. Helmi, Mon. Not. R. Astron. Soc. , 21 (Feb.2010), arXiv:0810.1522.[7] A. W. Graham, D. Merritt, B. Moore, J. Diemand,and B. Terzi´c, Astron. J. , 2701 (Dec. 2006),arXiv:astro-ph/0608613.[8] J. Diemand, B. Moore, and J. Stadel,Mon. Not. R. Astron. Soc. , 624 (Sep. 2004),arXiv:astro-ph/0402267.[9] W. J. G. de Blok, in Revista Mexicana de Astronomiay Astrofisica Conference Series , Revista Mexicana deAstronomia y Astrofisica, vol. 27, Vol. 17, edited byV. Avila-Reese, C. Firmani, C. S. Frenk, & C. Allen(2003) pp. 17–18.[10] G. Gentile, P. Salucci, U. Klein, D. Vergani, andP. Kalberla, Mon. Not. R. Astron. Soc. , 903 (Jul.2004), arXiv:astro-ph/0403154.[11] P. Salucci, A. Lapi, C. Tonini, G. Gentile, I. Yegorova,and U. Klein, Mon. Not. R. Astron. Soc. , 41 (Jun.2007), arXiv:astro-ph/0703115.[12] F. Donato, G. Gentile, P. Salucci, C. Frigerio Martins,M. I. Wilkinson, G. Gilmore, E. K. Grebel, A. Koch, andR. Wyse, Mon. Not. R. Astron. Soc. , 1169 (Aug.2009), arXiv:0904.4054 [astro-ph.CO].[13] B. S. Ryden and J. E. Gunn, Astrophys. J. , 15 (Jul.1987).[14] A. Nusser, Mon. Not. R. Astron. Soc. , 1397 (Aug.2001), arXiv:astro-ph/0008217.[15] N. Hiotelis, Astron. Astrophys. , 84 (Jan. 2002),arXiv:astro-ph/0111324.[16] L. L. R. Williams, A. Babul, and J. J. Dal-canton, Astrophys. J. , 18 (Mar. 2004),arXiv:astro-ph/0312002.[17] P. Sikivie, I. I. Tkachev, and Y. Wang, Phys. Rev. D ,1863 (Aug. 1997), arXiv:astro-ph/9609022.[18] A. Del Popolo, Astrophys. J. , 2093 (Jun. 2009),arXiv:0906.4447.[19] S. D. M. White and D. Zaritsky, Astrophys. J. , 1(Jul. 1992).[20] M. Le Delliou and R. N. Henriksen, Astron. Astrophys. , 27 (Sep. 2003), arXiv:astro-ph/0307046.[21] Y. Ascasibar, G. Yepes, S. Gottl¨ober, and V. M¨uller,Mon. Not. R. Astron. Soc. , 1109 (Aug. 2004),arXiv:astro-ph/0312221.[22] Y. Hoffman and J. Shaham, Astrophys. J. , 16 (Oct.1985).[23] F. Hoyle, in Problems of Cosmical Aerodynamics (1951)pp. 195–197. [24] E. Hayashi, J. F. Navarro, and V. Springel,Mon. Not. R. Astron. Soc. , 50 (May 2007),arXiv:astro-ph/0612327.[25] J. Barnes and G. Efstathiou, Astrophys. J. , 575(Aug. 1987).[26] M. Boylan-Kolchin, V. Springel, S. D. M. White, andA. Jenkins, Mon. Not. R. Astron. Soc. , 896 (Jun.2010), arXiv:0911.4484.[27] P. J. E. Peebles, Astrophys. J. , 393 (Feb. 1969).[28] S. D. M. White, Astrophys. J. , 38 (Nov. 1984).[29] A. G. Doroshkevich, Astrofizika , 581 (1970).[30] Ya. B. Zel’dovich, Astron. Astrophys. , 84 (Mar. 1970).[31] S. Chandrasekhar, Astrophys. J. , 255 (Mar. 1943).[32] O. E. Gerhard and J. Binney, Mon. Not. R. Astron. Soc. , 467 (Sep. 1985).[33] D. Merritt and G. D. Quinlan, Astrophys. J. , 625(May 1998), arXiv:astro-ph/9709106.[34] F. Cruz, H. Vel´azquez, and H. Aceves, Revista Mexicanade Astronomia y Astrofisica , 95 (Apr. 2007).[35] A. J. Kalnajs, in Dynamics of Disc Galaxies , edited byB. Sundelius (1991) pp. 323–+.[36] W. Dehnen, Astron. J. , 800 (Feb. 2000),arXiv:astro-ph/9911161.[37] M. Milosavljevi´c and D. Merritt, Astrophys. J. , 860(Oct. 2003), arXiv:astro-ph/0212459.[38] A. Sesana, F. Haardt, and P. Madau, Astrophys. J. ,546 (May 2007), arXiv:astro-ph/0612265.[39] J. F. Navarro, C. S. Frenk, and S. D. M.White, Astrophys. J. , 563 (May 1996),arXiv:astro-ph/9508025.[40] B. Moore, T. Quinn, F. Governato, J. Stadel, andG. Lake, Mon. Not. R. Astron. Soc. , 1147 (Dec.1999), arXiv:astro-ph/9903164.[41] M. Ricotti, A. Pontzen, and M. Viel, Astrophys. J. Lett. , L53 (Jul. 2007), arXiv:0706.0856.[42] R. Cen, F. Dong, P. Bode, and J. P. Ostriker, ArXiv As-trophysics e-prints(Mar. 2004), arXiv:astro-ph/0403352.[43] A. Del Popolo, Mon. Not. R. Astron. Soc., 1312(Sep.2010).[44] M. Takada, E. Komatsu, and T. Futamase, Phys. Rev. D , 083520 (Apr. 2006), arXiv:astro-ph/0512374.[45] E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett,B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. R.Nolta, L. Page, D. N. Spergel, M. Halpern, R. S. Hill,A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S.Tucker, J. L. Weiland, E. Wollack, and E. L. Wright,ArXiv e-prints(Jan. 2010), arXiv:1001.4538.[46] U. Seljak and M. Zaldarriaga, Astrophys. J. , 437(Oct. 1996), arXiv:astro-ph/9603033.[47] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay,Astrophys. J. , 15 (May 1986).[48] C. Ma and E. Bertschinger, Astrophys. J. , 28 (Sep.2004), arXiv:astro-ph/0311049.[49] M. Davis, G. Efstathiou, C. S. Frenk, and S. D. M. White,Astrophys. J. , 371 (May 1985).[50] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow, A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, andS. D. M. White, Mon. Not. R. Astron. Soc. , 1685(Dec. 2008), arXiv:0809.0898.[51] J. Gan, X. Kang, F. C. van den Bosch, and J. Hou, ArXiv e-prints(Jun. 2010), arXiv:1007.0023.[52] A. R. Zentner, A. A. Berlind, J. S. Bullock, A. V.Kravtsov, and R. H. Wechsler, Astrophys. J.624