{'fraction_non_alphanumeric': 0.04485130816315476, 'fraction_numerical': 0.018417203835111858, 'mean_word_length': 4.589826402906742, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Apparent Power Law in CDM Halo Pseudo Phase Space Density Profiles

We investigate the apparent power-law scaling of the pseudo phase space density (PPSD) in CDM halos. We study fluid collapse, using the close analogy between the gas entropy and the PPSD in the fluid approximation. Our hydrodynamic calculations allow for a precise evaluation of logarithmic derivatives. For scale-free initial conditions, entropy is a power law in Lagrangian (mass) coordinates, but not in Eulerian (radial) coordinates. The deviation from a radial power law arises from incomplete hydrostatic equilibrium (HSE), linked to bulk inflow and mass accretion, and the convergence to the asymptotic central power-law slope is very slow. For more realistic collapse, entropy is not a power law with either radius or mass due to deviations from HSE and scale-dependent initial conditions. Instead, it is a slowly rolling power law that appears approximately linear on a log-log plot. Our fluid calculations recover PPSD power-law slopes and residual amplitudes similar to N-body simulations, indicating that deviations from a power law are not numerical artefacts. In addition, we find that realistic collapse is not self-similar: scale lengths such as the shock radius and the turnaround radius are not power-law functions of time. We therefore argue that the apparent power-law PPSD cannot be used to make detailed dynamical inferences or extrapolate halo profiles inward, and that it does not indicate any hidden integrals of motion. We also suggest that the apparent agreement between the PPSD and the asymptotic Bertschinger slope is purely coincidental.


INTRODUCTION
What determines the final state of a large collection of particles interacting under Newtonian gravity? Despite decades of effort, a satisfactory answer to this question remains elusive. The difficulty is due primarily to the long range, unshielded nature of gravity. While gaseous systems with shortrange interactions quickly relax to a Maxwellian equilibrium, self-gravitating systems have a collisional relaxation time that scales with the number of particles, making two-body relaxation essentially negligible on galactic scales. Instead, collective processes such as phase-mixing and violent relaxation govern the approach to equilibrium. The fact that relaxation is incomplete is clear because the results of N-body simulations depend on initial conditions, such as the initial virial ratio (van Albada 1982). Indeed, the only distribution function consistent with complete relaxation results in a singular isothermal profile with infinite mass (Lynden-Bell 1967). The entire notion of thermodynamic equilibrium in a self-gravitating system is suspect, since at fixed mass and energy one can always arbitrarily increase the entropy of a system by increasing its central concentration (Lynden-Bell & Wood 1968;Tremaine et al. 1986). These problems arise directly from the fact that gravitational systems are nonextensive and non-additive 1 . Similar difficulties appear in other systems with long-range interactions for which standard Boltzmann-Gibbs statistics do not apply; examples include 2D hydrodynamic systems, 2D elastic systems, and charged or dipolar systems (Campa et al. 2009) Despite these theoretical difficulties, there is now overwhelming evidence from N-body simulations that structure formation results in nearly universal CDM halo profiles (see Frenk & White (2012) for a review); the most famous of these features is the density profile, as described by either the NFW (Navarro et al. 1997) or Einasto (Navarro et al. 2004) fitting formulae. The overall shape of the density profile is independent of the initial fluctuation spectrum, of halo mass and formation epoch, and of cosmological parameters. Indeed, similar profiles even arise in the absence of hierarchical growth (Huss et al. 1999;Moore et al. 1999;Wang & White 2009). The logarithmic slope of the density profile continually steepens from ∼ −1 or even shallower in the innermost regions to ∼ −3 at the halo outskirts. These universal density profiles in turn imply universal circular velocity profiles. In the original NFW fit, the only weak scale dependence is encapsulated by the halo concentration c = r200/r−2 (where r−2 is the radius at which d(logρ)/d(logr) = −2), which has been shown to reflect halo formation time.
Most work to date has focused on explaining the origin of the dark matter density profile, which for instance has been attributed to tidal disruption of substructure (Syer & White 1998;Subramanian et al. 2000;Dekel et al. 2003a,b), the shape of the matter power spectrum (Nusser 2001), a modulation of the accretion rate (Lu et al. 2006), the adiabatic contraction of the peaks of Gaussian random fields (Dalal et al. 2010), and conservation of orbital actions (Pontzen & Governato 2013). At this point, a consensus explanation seems unlikely (Frenk & White 2012). However, there are other intriguing regularities that provide insights into halo structure. Halos clearly show common trends in orbital anisotropy β(r) = 1 − σ 2 θ /σ 2 r : they are quasi-isotropic (β ≈ 0) at the centre and radially anisotropic outwards, with β ≈ 0.25 at r−2 and β ≈ 0.5 farther out Ludlow et al. 2010). Within r−2, the anisotropy parameter β(r) is correlated with the logarithmic slope of the density profile γ(r) via β(r) ≈ −0.15 − 0.2γ(r) (Hansen & Moore 2006). Similarly, the cumulative distribution of specific angular momentum j can be fit by a universal function M (< j) ∝ j (with a flattening at large j), albeit with large scatter (Bullock et al. 2001;Bett et al. 2010).
Most remarkably, the spherically averaged quantity Q(r) ≡ ρ/σ 3 follows an approximate power law Q(r) ∝ r −α , where α ≈ 1.9, over three orders of magnitude in radius over the entire relaxed halo profile (Taylor & Navarro 2001;Navarro et al. 2010;Ludlow et al. 2010). This relationship holds despite the fact that neither ρ(r) nor σ(r) are power laws; in fact, these quantities vary substantially from halo to halo. Like the NFW profile, Q(r) is robust to mergers and whether structure formation is hierarchical Wang & White 2009), and its slope only depends very weakly on the power spectrum (Knollmann et al. 2008), so its origin is likely unrelated to any of these. As r → 0, Q(r) asymptotically approaches a similar power law to the one predicted by the 1D model of self-similar collapse in Bertschinger (1985), which otherwise disagrees with simulations since it predicts incorrect asymptotic power-law slopes for the density and velocity dispersion profiles. While Q(r) has units of phase space density, it is in fact a ratio of moments of the distribution function 3 , and hence is often dubbed the 'pseudo phase space density'. A very similar power law holds if one substitutes the radial velocity dispersion for the total velocity dispersion (Dehnen & McLaughlin 2005;Navarro et al. 2010).
The origin of this remarkable power-law relationship is not known, and understanding it could pay considerable div-idends. In contrast to empirical fitting formulae such as the NFW and Einasto profiles, a scale-free power law emerges naturally in a similarity solution, generally as the result of conserved quantities. A true power-law Q(r) would therefore be a 'fundamental' universal feature of CDM halos, in the sense that its power-law index should be calculable from the equations and boundary conditions governing halo collapse. If this is the case, the relationship could be used when numerical resolution peters out. For instance, density profiles must differ from a simple Einasto form if a power-law Q(r) holds at all radii (Ma et al. 2009), but current numerical simulations cannot test this claim (Ludlow et al. 2011). In addition, many of the present universal features of equilibrium CDM halos would follow as a consequence of a power-law Q(r). In particular, the Jeans equation, which essentially describes hydrostatic equilibrium, has three free functions: ρ(r), σ 2 r (r), and β(r). Q(r) relates ρ(r) and σ 2 r (r) via an effective 'equation of state', so only β(r) needs to be specified to close the system of equations. Assuming orbital isotropy (β(r) ≈ 0), as appropriate for halo centres, the observed power-law Q(r) ∝ r −α with α ≈ 1.9 is the only solution of equation 1 that yields a physical (non-negative and non-truncated) density profile, and this profile agrees with the NFW fit (Taylor & Navarro 2001;Hansen 2004;Dehnen & McLaughlin 2005). Additional assumptions are needed for the more general anisotropic case, but if one assumes β(r) ∝ γ(r), as seen in numerical simulations, then one can recover density, velocity dispersion, and anisotropy profiles in good agreement with cosmological simulations, independent of the assumed slope of the β(r) − γ(r) relation (Dehnen & McLaughlin 2005).
A power law pseudo phase space density also has some observational support. From member galaxy kinematics and from lensing, X-ray, and kinematical mass profiles, a powerlaw Q(r) in agreement with simulations has been inferred for galaxy cluster Abell 2142 (Munari et al. 2014). Measurements of the mass density ρ(r) and the radial stellar velocity dispersion σr, * (r) of ∼ 2000 SDSS elliptical galaxies yield a power-law Q(r) ∝ r −χ , with χ = 1.860 ± 0.035, in remarkable agreement with values for simulated CDM halos (Chae 2014;Chae et al. 2014). This is particularly striking because density profiles for ellipticals are not universal and differ sharply from halo density profiles, largely because the stellar components of ellipticals are influenced by dissipative processes (e.g., gas cooling to form discs or the merger of disc galaxies to form a bulge). The fact that Q(r) follows a power law with a similar slope for elliptical galaxies, galaxy clusters, and dark matter halos reinforces the view that Q(r) is a more universal and fundamental feature of gravitational collapse than ρ(r).
The apparent power-law nature of Q(r) can therefore deliver powerful insights into the universality of CDM halo structure. However, a distinction must be drawn between approximate and exact power laws. Exact power laws often have deep physical significance; the associated scale invariance generally arises from conserved integrals of motion (e.g., energy conservation in Taylor-Sedov explosions). In the case of 'Type I similarity', the (rational) power-law exponent follows from simple dimensional analysis. In 'Type II similarity', the integrals of motion and the (irrational) power-law exponent can only be obtained by solving an eigenvalue problem; this type of self-similarity arises from invariance under the renormalisation group (Goldenfeld 1992;Barenblatt 1996). On the other hand, the power-law nature of Q(r) has been inferred simply by fitting noisy data on log-log plots, and it is notoriously easy to obtain spurious power-law relationships in this manner (e.g., see Clauset et al. (2009) for the plethora of pitfalls associated with inferring power-law distributions). State of the art simulations show ∼ 20% deviations from the best-fitting powerlaw Q(r), which moreover varies systematically with radius. Given the current numerical limitations, it is difficult to assess the significance of these deviations. Thus, simple toy models that can be solved to high precision (and benchmarked by comparison to exact self-similar calculations) are valuable. In this paper, we show that simple 1D calculations of fluid collapse satisfy these criteria. While these models do not possess all of the complex degrees of freedom encompassed by cosmological N-body simulations, they are very illustrative. We show that the entropy profiles from fluid calculations are remarkably similar to the pseudo phase space density profiles of CDM halos in N-body simulations; in particular, they have similar straight-line fits on log-log plots. However, these profiles deviate from true power laws for simple physical reasons, even for idealised collapse with scalefree initial conditions. We therefore argue that the power-law fits to Q(r) do not have deep physical significance. This paper is organised as follows. In §2, we discuss the applicability of the fluid approximation to CDM halos. In §3, we describe our method for simulating realistic halo collapse and extracting pseudo phase space density profiles via the fluid approximation. We also consider halo collapse with scale-free initial conditions and compare our simulation with an analytic calculation. We present our results in §4 and conclude in §5.

THE FLUID APPROXIMATION
In this paper, we use the 'fluid approximation' for dark matter. This technique utilises the fact that the zeroth, first, and second moments of the collisionless Boltzmann equation resemble fluid equations for the conservation of mass, momentum, and energy. The main difference with a true fluid that has scalar pressure p is that collisionless systems are supported against gravity by gradients in the (generally anisotropic) stress tensor ρσ 2 ij . The fluid approximation has been used in stellar dynamics to study the gravothermal catastrophe (Larson 1970;Lynden-Bell & Eggleton 1980), giving results in good agreement with numerical Fokker-Planck and N-body calculations. More recently, it has been used to study the collapse of collisionless (Teyssier et al. 1997;Subramanian et al. 2000;Lapi & Cavaliere 2011) and self-interacting (Ahn & Shapiro 2005) CDM halos. Given the same assumptions as the classic particle-based similarity solutions of Fillmore & Goldreich (1984) and Bertschinger (1985) (1D radial infall and scale-free initial conditions in an Einstein de-Sitter universe), fluid calculations correctly recover the same density and velocity profiles.
The solution procedure is as follows. Initially, the cold fluid obeys the standard parametric relations for turnaround and free-fall onto an overdensity (Peebles 1980). The hyperbolic fluid equations undergo a discontinuous 'shock' at the radius where shell-crossing would occur, and the singlestream flow becomes multi-stream. This can be treated by the standard Rankine-Hugoniot shock jump conditions. In the particle-based approach, this corresponds to the first caustic, and relaxation (i.e., phase-mixing) occurs in a thicker region, rather than instantaneously, due to successive shell-crossings. The shock thermalizes the fluid, converting ordered bulk motion into random motions.
The assumption that the infinite BBGKY hierarchy can be truncated at its second moment deserves justification. It amounts to the assumption that vr has a symmetric distribution, and in particular that it is skewless 4 : (vr −vr) 3 = 0. An asymmetric distribution implies net radial energy flow through the system, and in general skewness is zero except near caustics. In the particle approach, this can be seen directly from the symmetric distribution of radial velocities in phase diagrams. The hypothesis that the radial action is an adiabatic invariant (Fillmore & Goldreich 1984) is equivalent to our assumption that the fluid does not experience shell-crossing, since in both cases the enclosed mass within a shell is constant on orbital timescales (Teyssier et al. 1997).
For our purposes, the fluid approximation has a number of advantages. Self-similar calculations are straightforward in the fluid approximation; agreement with or departure from collisionless N-body results can shed further insights. Indeed, the case of pure fluid collapse is an interesting problem in its own right and can be tested against hydrodynamic simulations. Furthermore, since we focus on the pseudo phase space density Q(r), the fluid approximation allows us to exploit a very interesting and useful analogy with the gas entropy kg = 3kBTgρ −2/3 g (µmp) −1 = σ 2 g ρ −2/3 g . Simulations of galaxy clusters have shown that, outside the central core (where gas entropy is often flat 5 ), the analogous quantity kDM = σ 2 DM ρ −2/3 DM matches the gas entropy in both its power-law radial scaling (kDM ∝ kg ∝ r 1.2 ) and in its normalisation, once gas bulk motions are taken into account (Faltenbacher et al. 2007). A similar power-law entropy profile kg ∝ r 1.21±0.39 is inferred from X-ray observations (Cavagnolo et al. 2009). A quick heuristic way to understand the power-law gas entropy profile follows from the accretion history of the cluster (Tozzi & Norman 2001;Voit et al. 2003). The gas entropy is a power-law function of mass, vin is the infall velocity at the shock and M ∝ t η with 0.9 η 1.9 for a Ωm = 0.3 flat cosmology with the relevant power spectrum indices −2 n −1 (Voit et al. 2003). This estimate agrees with our more detailed calculations (for example, see Figure 7). A running power-law scaling of kg with radius then follows by demanding hydrostatic equilibrium. Using the fluid analogy -that dark matter entropy is generated at 'shocks' (i.e., caustics) via phase-mixing and is thereafter conserved -we can construct pseudo phase space density profiles for dark matter halos and compare against simulations using Q(r) ∝ k −3/2 An important caveat is in order. The fluid model does not give a clear-cut prescription for the velocity anisotropy β(r); in particular, it does not specify the directionality of the random motions generated at the shock (of course, in a true fluid, velocities are isotropic). Tangential velocity dispersions arise from asphericity, tidal torques, and the radial orbit instability (ROI). The ROI is particularly interesting: it arises because purely radial orbits are unstable to precessional motion, leading to a quadrupolar bar-like instability. This in turn generates tangential velocity dispersions that are ultimately stabilising (Antonov 1973;Polyachenko & Shukhman 1981;Merritt & Aguilar 1985;Palmer & Papaloizou 1987). The ROI clearly plays an important role in shaping halo density profiles; in simulations of monolithic collapse where only radial forces are evolved, the density profile is close to isothermal. However, once the ROI is allowed to operate, orbits near the centre are isotropised, the central cusp softens, and density profiles transform into an NFW form (Huss et al. 1999;MacMillan et al. 2006). Unfortunately, there is currently no successful analytic theory for the ROI (see Merritt (1999); Binney & Tremaine (2008) for reviews) or how it saturates (Adams et al. 2007). To incorporate these effects, semi-analytic models have inserted tangential velocity dispersions and/or angular momentum by hand, generally near turnaround (Nusser 2001;Barnes et al. 2005;Lu et al. 2006;Zukin & Bertschinger 2010); suitably tuned, this method reproduces the NFW profile. Fortunately, for our goal of a deeper understanding of the power-law nature of Q(r), a detailed understanding of the ROI is unnecessary. A power-law Q(r) is remarkably robust and arises (both in ρ/σ 3 and ρ/σ 3 r ) in purely radial collapse (Bertschinger 1985), in isotropic fluid collapse, and in N-body simulations where the ROI develops naturally (MacMillan et al. 2006) or where tangential velocities are manipulated by hand Lapi & Cavaliere 2011). It is therefore unlikely to be sensitive to the details of the ROI.
Even if we want to understand the density profile, the detailed shape of β(r) may not be necessary. From a phaseplane analysis of the Jeans equation and its derivatives, Dehnen & McLaughlin (2005) have shown that, if we assume a power-law Q(r) and that β(r) ∝ γ(r) ≡ d logρ/d logr (as seen in simulations; Hansen & Moore (2006)), then we only need to know β0 = β(r = 0) to construct an analytic density profile in excellent agreement with simulations. This result is independent of the assumed slope of the β − γ relation. The velocity distribution is close to isotropic near halo centres (β0 ≈ 0, with some scatter) in simulations; this can be understood as the effect of efficient violent relaxation during the rapid collapse phase (Lu et al. 2006), and it is related to the slope of the initial density perturbation (Vogelsberger & White 2011). Also, for the purposes of understanding the central cusp (and extrapolating inward beyond the limits of numerical resolution), the details of β(r) at large radii are irrelevant.
In this paper, we only consider spherical symmetry and monolithic collapse. Although adding triaxiality does affect the density and anisotropy profiles (Lithwick & Dalal 2011), Q(r) is virtually unaffected, remaining a single power law with an index that depends very weakly on the similarity index of the initial conditions (Vogelsberger & White 2011). This perhaps hints that entropy stratification is governed primarily by the gravitational potential, which has greater spherical symmetry than the density profile. The hierarchical nature of structure formation does not appear to be important in determining the structure of CDM halos. As noted above, key features (e.g., NFW-like density profiles and power-law Q(r) profiles) already arise in monolithic collapse simulations. Furthermore, simulations of equal mass merger events -the most violent possible -show that mixing is remarkably inefficient and that the overall shapes of the progenitor Q(r) profiles (and the coarse-grained phase space densityf (x, v)) are remarkably well-preserved (Vass et al. 2009), regardless of the merger geometry or the number of mergers.

Self-similar Fluid Collapse
In an Einstein-de Sitter (EdS) universe with scale-free initial conditions, the evolution of density perturbations is selfsimilar. As we demonstrate below, the fluid equations then reduce from a set of PDEs to a set of ODEs that can easily be solved to very high precision. This is useful for two reasons: (i) We can test the reliability of our hydrodynamic code for initial conditions where the exact solution is known; (ii) We can study the entropy profile kg(r) for self-similar collapse. If the entropy profile is a power law, then we can study how it changes if exact self-similarity is broken; if it is not a power law, then kg(r) is unlikely to be a power law in the more realistic case! Our treatment closely follows that of Bertschinger (1985), who considered radial self-similar fluid collapse. The main new element in our analysis is the focus on the entropy profile kg(r).
As in Fillmore & Goldreich (1984), scale-free initial density perturbations can be characterised by where Mi is the unperturbed mass, M0 is a reference mass, and > 0 is a constant. The positive mass excess implies that each radial shell is bound and collapses after reaching a maximum turnaround radius rta, which scales with time according to Here and throughout this paper, α denotes the power-law index of characteristic scale radii as functions of time (note that α is constant for self-similar collapse). We can derive equation 3 by dimensional analysis as follows. First, note that the initial conditions of equation 2 imply that the specific energy δE ∝ GδM/r ∝ M 1− /M 1/3 ∝ M 2/3− , where E =ṙ 2 /2 − GM/r is the energy integral of the equation of motion for each shell. The maximum radius each shell reaches satisfies δE = GM/rta, so rta ∝ M/M 2/3− = M 1/3+ . The only characteristic time scale is the free-fall shock radius of the infalling matter scales identically by selfsimilarity.
Since there are no other characteristic scales in the problem, the behaviour of the system is self-similar once rescaled to the turnaround radius or the shock radius. In analogy with Bertschinger (1985), we non-dimensionalize the problem by factoring out the time-dependent background densityρ(t) and a characteristic scale radius rs(t): The usual choice for the scale radius is the turnaround radius rta(t), but in this paper we will use the shock radius for rs(t). Thus, λ refers to the radius normalised to the radius of the fluid shock. Note that the ratio rs/rta is constant for self-similar collapse, but this is no longer true for more realistic initial conditions ( §3.3). Using equations 4-8, the fluid equations in spherical symmetry become: where the non-dimensionalized velocity, density, pressure, and mass profiles are functions of λ and primes denote total derivatives with respect to λ. Note that M (λ) is the mass interior to a fluid element at radius λ.
The above system of ODEs can be solved numerically. We describe the procedure using the conventional turnaround radius non-dimensionalization. We assume that at the initial time ti the radial velocity field is pure unperturbed Hubble flow. The cold pre-shock fluid is pressureless (P = 0), so there are only 3 variables: D, V , and M . The initial conditions for these profiles at turnaround are V = 0, M = (3π/4) 2 , and D = M/(3 + 1) (Fillmore & Goldreich 1984). The equations can be integrated inward from the turnaround radius to the non-dimensionalized shock radius λs = rs/rta, whereupon the shock jump conditions apply, and the integration can continue inward. Note that αλs is the non-dimensionalized shock speed. The unknown shock eigenvalue λs is varied until the boundary conditions are satisfied. M (0) = 0 guarantees that the mass profile is non-singular -i.e., there is no central point mass (such as a black hole). V (0) = 0 immediately follows from this requirement, since a nonzero V (0) in a collisional fluid with no shell-crossing leads to the development of a mass singularity.
Throughout this paper, we use an adiabatic index γ = 5/3. The fact that fluid calculations with γ = 5/3 match calculations of collisionless infall was first noted by Bertschinger (1985). He remarked that one might expect a γ = 3 fluid calculation, corresponding to a gas with only one translational degree of freedom, to agree with the 1D collisionless calculation for purely radial infall. This is not the case, since the fluid equations as written above assume an isotropic 3D pressure and use a 3D mass density. Some authors (Subramanian et al. 2000;Lapi & Cavaliere 2011) have subsequently formulated the self-similar calculation in 1D with γ = 3 (i.e., using a radial pressure and density), but we use the conventional formulation.
Equations 9-12 give the mass and entropy integrals (Bertschinger 1985 where we have defined the fluid entropy kg = P D −γ . The mass and entropy integrals should be conserved exactly, so they are useful for checking the fidelity of our numerical results.

Fluid Collapse Simulation
To simulate dark matter halo collapse, we use a 1D spherically symmetric Lagrangian hydrodynamics/gravity code, similar to Thoul & Weinberg (1995). We initialise concentric fluid shells with a density profile ρi(r) and a radial velocity profile vi(r) corresponding to pure unperturbed Hubble flow at a fiducial initial redshift zi. The self-gravitating system then evolves according to the fluid equations; at the radius where shell-crossing would occur, the fluid undergoes a discontinuous shock, which we treat with an artificial viscosity technique as in Thoul & Weinberg (1995). We allow the system to evolve to a final redshift z f , which we take as z f = 0 throughout. We choose M0 = 10 12 M for our reference mass.
We perform several tests of our code. First, we initialise N = 1000 shells in an EdS universe with a Gaussian initial overdensity profile δi(r) = [ρi(r) −ρi]/ρi = δi(0)e −r 2 /r 2 i , where ri is a scale radius, following Thoul & Weinberg (1995). Note that for r ri, δM → const, so δM/M ∝ M −1 and this is equivalent to an = 1 mass perturbation. The resulting non-dimensionalized shell trajectories and fluid variable profiles (Figure 1) agree with the results from Figures  3 and 4 of Thoul & Weinberg (1995) and with the analytic self-similar solution for = 1 in Bertschinger (1985). We also reproduce Figures 6 and 7 of Thoul & Weinberg (1995) for collapse from initial conditions generated by an n = −2 power-law power spectrum in an EdS universe. More generally, we have compared our hydrodynamic code to selfsimilar calculations ( §3.1) for arbitrary , finding excellent agreement (for example, see Figure 6). The integrals of motion (equations 18, 19) are conserved to high precision.  Thoul & Weinberg (1995) for self-similar fluid collapse with a Gaussian initial overdensity profile δ i (r) = δ i (0)e −r 2 /r 2 i , where r i is a scale radius. Top: all shells follow the same trajectory in nondimensionalized coordinates. Bottom: non-dimensionalized velocity, pressure, mass, and density profiles at two snapshots in time (blue and green lines). These profiles agree with our analytic results.

Realistic Initial Conditions
We now describe our method for simulating dark matter halo collapse with realistic initial conditions in a ΛCDM universe with {ΩΛ, Ωm} = {0.721, 0.279}. To generate realistic initial conditions, we follow Lu et al. (2006) by using linear perturbation theory to construct initial overdensity profiles that correspond to realistic halo mass accretion histories. In particular, we use the universal form for halo mass accretion histories found in Wechsler et al. (2002) as a fit to numerical simulations: where the parameter zc is the redshift at which the mass accretion rate d(logM )/d(loga) falls below a critical value S = 2 and effectively transitions from 'fast' to 'slow' accretion. A spherical shell collapses when its average linear overdensity reaches a critical value δc ≈ 1.686, so for a is the linear growth factor (fit by Carroll et al. 1992) and z(M ) is given by inverting equation 20. Furthermore, the mean mass density at redshift zi isρ(zi) = ρcrit,0Ωm(1+zi) 3 , where ρcrit,0 is the critical density of the universe at z = 0. We choose N = 4000 equally spaced shells for our simulation; an example of the initial overdensity profile generated for a M0 = 10 12 M halo with zc = 3 and zi = 500 is shown in Figure 2. Our conclusions are insensitive to the values of M0, zc, and zi, as long as zi >> zc so that the linear perturbation technique is valid.
In Figure 3, we show the (M ) profiles for zc = 0, 1, and 3, where (M ) ≡ −d(log[δM/M ])/d(logM ) (motivated by δM/M ∝ M − in the scale-free case) and M is evaluated at z = 0. Though ≈ 1/6 in the inner regions of the collapsed halos (as we expect from the n eff ∼ −2.5 power-law index of the power spectrum on galactic scales), (M ) varies substantially over each halo profile. In this paper, we choose zc = 3 as our fiducial case; this value is characteristic of halos with M halo ∼ 10 12 M at z = 0 (e.g., see Figure 11 of Wechsler et al. (2002)). Note that zc = 3 is similar to the value chosen in Shapiro et al. (2006), where a comparable spherical collapse model is applied. We obtain very similar results for halos with different values of zc.
To simulate halo collapse in a ΛCDM universe, we model the effective gravitational acceleration of the jth mass shell as M is the enclosed mass at z = 0 normalised to the mass enclosed at the shock radius. We choose zc = 3 (blue) as our fiducial case throughout this paper.
where rj is the radius of the shell and Mj is the mass it encloses. The initial velocity profile is purely radial and corresponds to the unperturbed Hubble flow at redshift zi. Figure  4 shows the z = 0 density and circular velocity profiles for collapse from the realistic initial conditions of Figure 2. Our results are fit reasonably well by the corresponding NFW profiles, which shows that the universality in halo density and circular velocity profiles can be understood as the result of realistic initial conditions 6 . The NFW concentration parameter that fits our results (c = 16) is large for a M0 = 10 12 M halo (see Ludlow et al. 2016). We therefore explore the effects of varying the adiabatic index γ in order to crudely model velocity anisotropy, which is a crucial effect in shaping density profiles with lower values of c (Lu et al. 2006). We find that fluid collapse with γ = 3 -i.e., a radial velocity dispersion -yields a much more reasonable concentration parameter c ∼ 6. Varying γ shellby-shell can in principle account for velocity anisotropy; the large value we obtain for c is therefore an artefact of our fluid model. Since our focus is on the entropy profile, we sidestep this issue in this paper.
The density profile in Figure 4 differs from an NFW profile in detail, since our 1D fluid calculation only approximates 3D collisionless dynamics; however, the inner and outer slopes of the density profile are consistent with the NFW values. Moreover, we will see that the entropy profiles produced by our 1D model correspond to pseudo phase space density profiles that closely resemble those found in Nbody simulations. This suggests that the detailed behaviour of the density profile is largely irrelevant for understanding whether the pseudo phase space density is a power law.  As discussed above, we can construct the fluid analogue of the pseudo phase space density using Q(r) ∝ k −3/2 g . Figure 5 shows the z = 0 pseudo phase space density profiles that correspond to the three initial conditions from Figure  3. The profiles agree reasonably well with the power law favoured by N-body simulations over the dynamic range shown; Shapiro et al. (2006) reach a similar result and conclude that the universality in Q(r) profiles can be understood as the result of realistic initial conditions, even in the context of simple spherical collapse models. In particular, they argue that a power-law fit to Q(r) very similar to those obtained from N-body simulations results from a realistic time-varying mass accretion history in a ΛCDM universe. However, we will find that the Q(r) profiles produced by our spherical collapse model deviate significantly from power laws over the dynamic range relevant to N-body simulations. In particular, the logarithmic derivatives of the curves in Figure 5 are not constant functions (see Figure 8). Thus, by carefully examining the logarithmic derivatives of our pseudo phase space density profiles, we arrive at different conclusions than Shapiro et al. (2006).

Self-similar Fluid Collapse
We first examine the entropy profiles for collapse from scalefree initial conditions in an EdS cosmology, which results in self-similar collapse ( §3.1). These calculations provide much of the physical insight required to understand more realistic entropy profiles, which we consider in §4.2. We study scalefree collapse with = 1/3 (green lines) and 1/6 (red lines) in detail in Figures 6-9. We choose these values for detailed study because they correspond to the central (M → 0) and Figure 5. Pseudo phase space density profiles Q(λ) ∝ k −3/2 g at z = 0 for three halos that collapse from realistic initial conditions at z i = 500 (λ is the radius normalised to the fluid shock). The profiles seem to follow a universal power law that agrees reasonably well with the fit from N-body simulations. However, these profiles are not true power laws; see Figure 8.
average values of (M ) for our fiducial realistic profile with zc = 3 (Figure 3; note that = 1/6 is more representative for lower values of zc). Just as for more realistic initial conditions, scale-free collapse yields power-law-like pseudo phase space density (PPSD) profiles with promising straight-line fits on log-log plots. We now list our main results.
Entropy is a power-law function of mass, but not of radius. From equation 19, we expect the gas entropy to be an exact power-law function of mass, since entropy is a Lagrangian conserved quantity in an adiabatic fluid. In particular, for scale-free collapse in an EdS universe, where n = 1, 2, 3 for planar, cylindrical, and spherical perturbations (see, e.g., Chuzhoy & Nusser 2000); for n = 3, this reproduces equation 19. The lower panel of Figure 7 shows that our fluid simulation recovers this result. The fact that this integral of motion is conserved to high precision demonstrates the accuracy of our hydrodynamic simulation. However, entropy is not a power-law function of radius. In Figure 6, we show the logarithmic entropy derivative d(logkg)/d(logλ) from our analytic calculation ( §3.1) and our hydrodynamic simulation ( §3.2) in an EdS cosmology. The close agreement between the profiles lends confidence to the the fidelity of our simulation. Thus, we only present simulation results in Figures 7 and 8. The key takeaway from Figure 6 is that d(logkg)/d(logλ) is not constant, so the gas entropy is not a power-law function of radius, even for self-similar collapse! Departures from hydrostatic equilibrium (HSE) are responsible for the deviation from a radial power law. Using kg = P D −γ and equation 11, we obtain where λ is the radius scaled to the fluid shock and α is the power-law index of characteristic scale radii as functions of time (equation 3). Equation 24 implies that, if the halo is in hydrostatic equilibrium (V = 0), the entropy profile is an exact power law. Otherwise, for a nonzero radial velocity V (λ), the entropy only attains a constant logarithmic slope asymptotically 7 , as λ → 0, recalling the boundary condition V (0) = 0. Deviations from hydrostatic equilibrium therefore cause the power-law slope of the entropy profile to run systematically with radius. It is easy to see that violations to hydrostatic equilibrium should be sufficient to cause the radial power-law index of the entropy profile, α λ ≡ d(logkg)/d(logλ), to 'roll'. For α λ to be constant, we require V αλ. However, at the shock, V ∼ O(αλ), where αλ is the non-dimensionalized shock speed. V slowly declines throughout the halo as fluid elements are decelerated by the excess pressure gradient, until V (0) = 0. We note that none of the fluid variables for self-similar collapse are power-law functions of radius, even though they might appear to be straight lines on a loglog plot (e.g., see Figure 1). These quantities only asymptotically approach power laws as λ → 0. In this limit, V → V0λ, D ∝ λ δ , P ∝ λ η , and M ∝ λ δ+η , where V0 = {4(6 − 1)/(45 ), 0}, δ = {−3(3 + 2)/(3 + 7), −9 /(3 + 1)}, and η = {0, 2(1 − 6 )/(3 + 1)} for { 1/6, > 1/6} respectively 8 (Chuzhoy & Nusser 2000).
Convergence to HSE -and to the asymptotic power-law entropy slope -is slow. We can understand this by examining the velocity profiles. In Figure 9, we plot the quantity |log[1 − V /(αλ)]| versus λ (in the realistic case, α is determined as a function of λ by using the realistic (λ) profile in equation 3). For the entropy profile to approach a power law, this quantity must approach 0 (for V → 0) or a constant (for V → V0λ). Figure 9 shows that the velocity profiles approach their asymptotic values very slowly. In addition, the convergence to HSE is faster for larger values of , corresponding to more sharply peaked initial profiles. The slow approach to hydrostatic equilibrium, and the trend with , are consistent with the lower panel of Figure 8. In particular, the entropy profiles in Figure 8 approach their asymptotic central values (shown by the dashed lines) very slowly (i.e., at small λ), and the convergence is faster (i.e., it occurs at larger λ) for larger values of .
Note that the close similarity between the d(logkg)/d(logλ) profiles for = 1/3 and 1/6 is coincidental. The asymptotic slope for = 1/6 is steeper than that for = 1/3 (as shown by the red and green dashed lines in Figure 8), but the = 1/6 profile approaches its asymptotic value more slowly due to the slower convergence to HSE. In the radial range we study, these two effects happen to cancel. The accidental nature of this cancellation is highlighted by the differing d(logkg)/d(logλ) profiles for = 2/3 (purple) and = 1 (light blue) in the bottom panel of Figure 8.
The convergence rate is linked to the accretion rate.
7 A power law is also obtained if V ∝ λ. In practice, this only occurs in a neighbourhood near λ = 0, where the boundary condition V (0) = 0 guarantees that the linear term in a Taylor expansion is a good approximation. 8 Note the dichotomy between 1/6, where the central density continuously increases with time, and > 1/6, where the central density is constant (which can be seen from V 0 = 0). = 1/6 Figure 6. Logarithmic derivatives of the z = 0 entropy profiles from our analytic calculation (dashed lines) and our simulation (solid lines) for collapse from scale-free initial conditions δM/M ∝ M − , with = 1/3 (top) and = 1/6 (bottom). λ is the radius normalised to the fluid shock. Note that the physical shock radius becomes too small to resolve using our analytic method for log(λ) −2.0.
Deviations from hydrostatic equilibrium are related to the mass accretion rate. Given sufficient time, an isolated selfgravitating object will settle into perfect hydrostatic equilibrium. Mass accretion acts as a forcing function that promotes mass inflow (V = 0) and deepening of the potential well; the slower the mass accretion rate, the weaker the nature of this gravitational forcing. Since M ∝ t 2/(3 ) , larger values of correspond to smaller mass accretion rates, and thus smaller departures from HSE. This leads to faster convergence to asymptotic power-law entropy slopes for larger values of , as we noted above. These results lead to the strong suspicion that, barring some fortuitous cancellation, entropy profiles for more realistic scenarios should not be power-law functions of radius. We now examine the more realistic case.

Realistic Initial Conditions
We consider the realistic initial conditions and cosmology described in §3.3. There are now two features that break self-similarity. First, the transition from an EdS to a ΛCDM cosmology removes the scale-free nature of Hubble expansion (for instance, the expansion rate and the mean density are no longer power-law functions of time). In addition, as we have seen in §3.3, the mass excess δM/M is no longer a power law; instead, it is a rolling power law with an index (M ) (see Figure 3). As noted above, we choose zc = 3 as our fiducial case, although we obtain very similar results for initial conditions with different values of zc. Our results are as follows.
Entropy is no longer a power-law function of mass. Since entropy is a power-law function of mass for self-similar collapse, we first check whether this remains true for more realistic initial conditions. Figure 7 shows that kg(M ) is not a power law. Instead, the power-law index of the entropy profile as a function of mass, αM ≡ d(logkg)/d(logM ), rolls with mass in both an EdS and a ΛCDM cosmology. Since the results are similar in either cosmology, this behaviour must arise from the scale-dependent initial conditions. Indeed, we can largely reproduce the realistic αM profile by substituting the instantaneous value of α( ) (equation 3) into equation 23, which gives for spherical collapse with γ = 5/3. This relation can be understood intuitively from the argument in §2: the gas entropy kg ∝ T /ρ 2/3 ∝ v 2 in /ρ 2/3 ∝ (GM/R)/ρ 2/3 ∝ (M t) 2/3 ∝ M 2/3(1+3 /2) ∝ M +2/3 , where we have used M ∝ t 2/(3 ) in the penultimate step. Equation 25 is an approximate equality because α is only well defined for scale-free collapse. However, this relation seems quite robust; the αM profile predicted using the realistic (M ) profile in equation 25 matches our simulation to better than 15% over the entire collapsed halo (see Figure 7). In fact, this approximation largely reproduces the shape of the realistic αM profile and differs mainly in its normalisation.
From equation 25, we see that the slowly rolling d(logkg)/d(logM ) profile is directly related to (M ). The latter rolls slowly with mass ( Figure 3) because the excess mass δM is an integral quantity that cannot change abruptly 9 if δM M , as required for linear theory to be valid 10 . Thus, kg ∝ M ζ , where the mean index ζ ≈ 1, with a maximum ∼ 20% deviation on either end (see Figure 10).
Entropy is not a power-law function of radius. The previous results lead us to believe that the realistic entropy profile should not be a radial power law; Figure 8 shows that the radial power-law index of the entropy profile rolls as a function of λ, as we expect. The departure from a power law is somewhat larger than for self-similar collapse, since two effects now contribute: deviations from hydrostatic equilib- rium, and the fact that the collapse is no longer scale free. This result strongly suggests that the PPSD is not a power law for realistic CDM halos.
Interestingly, Figure 9 shows that the relative departure from hydrostatic equilibrium is smaller for realistic collapse than for self-similar collapse with = 1/3 or 1/6. This makes sense in light of the fact that departures from HSE are governed by the accretion rate. The current accretion rate is low for the realistic case (for zc = 3, it is more similar to the accretion rate for an = 3/5 perturbation; see Figure 3), so we expect smaller departures from HSE. Hydrostatic equilibrium is therefore a very good approximation deep within realistic halos. This explains, to the extent that ≈ const in these inner regions (see Figure 3), why Q(r) is an approximate power law. However, despite the fact that hydrostatic equilibrium is a good approximation deep within the halo, the realistic profile deviates from hydrostatic equilibrium by ∼ 10% at log(λ) = −0.1, with the deviation increas- ing sharply near the shock 11 . These values are comparable to departures from the equilibrium Jeans equation found in N-body simulations (e.g., Austin et al. (2005)).
Cosmology makes little difference. The upper panel of Figure 7 shows that scale-free collapse in a ΛCDM cosmology results in kg(M ) profiles that deviate from power laws in the outer mass shells (which are accreted at low redshifts), as one might expect. However, the magnitude of this deviation is small compared to the variation in the realistic entropy profile. Thus, the departure from self-similarity is dominated by the departure from scale-invariant initial conditions. Indeed, in Figure 8, the power-law index d(logkg)/d(logλ) is very similar in both an EdS and a ΛCDM cosmology -the choice of cosmology makes little difference. We can under-11 The size of the deviation at log(λ) = −0.1 depends weakly on the scheme used to smooth the ∇P term in the hydrostatic equilibrium equation, but the conclusion that the realistic profile deviates more strongly than the scale-free profiles is unaffected. . |log[1 − V /(αλ)]| for collapse from realistic initial conditions in a ΛCDM cosmology (blue) and from our analytic calculation of scale-free collapse with = 1 (light blue, dot-dashed), = 2/3 (purple, dot-dashed), = 1/3 (green), and = 1/6 (red).
stand this by noting that every overdense shell that eventually collapses behaves like a closed (Ω > 1) universe, so ΩΛ contributes a subdominant amount to the total energy of each bound shell.
Deviations from power-law entropy profiles are similar to N-body results. The deviations of our kg(M ) and kg(λ) profiles from their respective best-fitting power laws are shown in Figure 10; the amplitude of our realistic residual profile is in good agreement with the amplitudes of PPSD power-law residuals found in simulations Ludlow et al. 2010). Note that, unlike 3D N-body simulations, our idealised fluid model is much less subject to issues of numerical resolution, and it resolves entropy profiles to high precision. Our results imply that the 10-20% deviation from a power-law PPSD found in simulations is a real effect. Indeed, plotting the logarithmic slopes produced by our fluid model (Figures 7 and 8) as histograms and fitting with Gaussian distributions yields standard deviations of 14.4% for d(logkg)/d(logM ) and 12.3% for d(logkg)/d(logλ). The characteristic pattern of residuals in Figure 10 stems from the progressive steepening (flattening) of the entropy profile with mass (radius) (Figures 7 and 8, respectively). Thus, the best-fitting power-law profile for kg(M ) is too steep in the outer regions and too shallow in the inner regions, and vice versa for kg(λ).
We note that the power-law residuals for kg(λ) are very similar for scale-free and realistic initial conditions, even though the departure from a power law is determined entirely by bulk inflow in the scale-free case. For realistic collapse, smaller deviations from HSE ( Figure 9) compensate for the lack of scale-invariant initial conditions.
Rough agreement with the Bertschinger power-law index is a red herring. We regard the fact that the PPSD power-law fits are very close to the Bertschinger (1985) prediction Q(r) ∝ r −1.875 (since Q ∝ k −3/2 g , this corresponds to kg ∝ r 5/4 ) to be purely coincidental and of no real physical significance. The Bertschinger index applies to the secondary infall model ( = 1), which is well outside the range of (M )  Figure 10. Entropy profiles as a function of mass (top) and radius (bottom), along with best-fitting power laws (dot-dashed lines) and residuals, for collapse from realistic initial conditions in a ΛCDM cosmology (blue). Analytic results for scale-free collapse with = 1/3 (green) and = 1/6 (red) are shown in the bottom panel. λ is the radius normalised to the fluid shock.
for a realistic overdensity profile (see Figure 3). Moreover, the Bertschinger result is an asymptotic power-law index. As we have seen, convergence to asymptotic slopes is very slow and occurs at small radii that are inaccessible to N-body simulations. The logarithmic entropy slope for realistic collapse has a broad range of values ( Figure 8); the Bertschinger value just happens to lie in the middle of this range. Realistic collape is not self-similar. The realistic case does not have any 'hidden' self-similar or scale-free behaviour. Neither kg(M ) nor kg(λ) are power laws. Moreover, the lack of self-similarity can be seen from the fact that scale radii such as the turnaround radius rta(t) and the shock radius rs(t) are not power-law functions of time; Figure 11 shows the running power-law behaviour. The powerlaw scaling r(t) ∝ t α is a requirement for self-similar solutions (Barenblatt 1996), indicating the presence of a conserved quantity A ≡ r(t)/t α ; the exponent α can either be derived by dimensional analysis (Type I similarity; for example, energy conservation in a Taylor-Sedov explosion) or by solving an eigenvalue problem (Type II similarity; for exam-ple, the implosion of a spherical shock wave, as in Guderley (1942)). The fact that rta(t) and rs(t) are not power laws strengthens the case that realistic collapse does not display any 'intermediate asymptotics' (Barenblatt 1996), which allow convergence to a self-similar solution from initial conditions that are not scale free.

DISCUSSION
Our bottom line is easy to summarise: the pseudo phase space density profile in CDM halos is only a pseudo power law. The 10-20% deviation from a power-law PPSD found in N-body simulations is a real effect with physical origins, including deviations from hydrostatic equilibrium within collapsed halos and scale-dependent initial conditions. The fact that the PPSD is intrinsically different from a power law is somewhat obvious based on the results for scale-free fluid collapse: even scale-free initial conditions produce entropy profiles that are not power laws, so we should not expect a power law in the realistic case.
We find that the analogous Lagrangian quantity kg(M ) is robust to variations in the velocity profile and instead depends primarily on the initial conditions for halo collapse. Even in a ΛCDM cosmology, scale-free initial conditions produce kg(M ) profiles that are very close to power laws (while they are exact power laws in an EdS cosmology). The realistic kg(M ) profile is mainly determined by the function (M ) ≡ −d(log[δM/M ])/d(logM ) (motivated by δM/M ∝ M − in the self-similar case), which specifies the initial conditions. In particular, there is an intriguing linear relation between (M ) and the logarithmic slope of the entropy profile: d(logkg)/d(logM ) = (M ) + 2/3 (equation 25). This relation agrees with our realistic simulation to better than 15% over the entire collapsed halo.
An obvious caveat to our conclusions is that we only simulate fluid collapse in 1D, which differs in detail from collisionless collapse in 3D. At the very least, the fluid approximation ignores the effect of velocity anisotropy. It is conceivable that velocity anisotropy somehow conspires to create a true power-law PPSD, rather than a pseudo power law 12 . While we cannot rule out this possibility, we regard it as unlikely. In the absence of a plausible mechanism for accomplishing this result (which has not been proposed), remarkable fine-tuning would be required. Our fluid model gives power-law slopes ( Figure 5) and residuals ( Figure 10) consistent with those from collisionless N-body CDM simulations . Likewise, we find deviations from hydrostatic equilibrium that are comparable to deviations from the equilibrium Jeans equation found in N-body simulations (Austin et al. 2005).
The fact that the PPSD is not an exact power law, and that characteristic scales such as the turnaround radius and the shock radius are not power-law functions of time (a fundamental requirement for self-similarity), has several consequences. It implies that the PPSD does not have an elevated status relative to the density or velocity profiles, which show 12 For instance, velocity anisotropy clearly affects the density profile, reducing the halo concentration relative to fluid models. r ta r s Figure 11. d(logrta)/d(logt) (blue) and d(logrs)/d(logt) (green) from our simulation of halo collapse with realistic initial conditions in a ΛCDM universe. rta(t) and rs(t) are not power-law functions of time, which implies that realistic halo collapse does not converge to a self-similar solution.
approximate -but not exact -universality 13 . It also implies that we cannot use a strict power-law fit to make detailed dynamical inferences about halo structure, or to extrapolate halo profiles inward, beyond the limits of numerical resolution. Most importantly, since self-similar behaviour is always associated with conserved quantities, the PPSD is unlikely to be a smoking gun indicator of some hitherto unknown integral of motion (such an unknown conserved quantity has been invoked to explain the apparent universality in CDM halo structure, e.g., Pontzen & Governato (2013)). Instead, the PPSD is a rolling power law that runs very slowly and appears approximately linear on a log-log plot because it is a ratio of correlated quantities. This interpretation is disappointing, but it would not be the first time that an apparent power-law -thought to be of deep physical significance -was later found to be a coincidence. One such example is the two-point galaxy correlation function ξ(r). From the 1970s to 2005, galaxy surveys found a power-law ξ(r) that showed no new features in the transition from linear to non-linear scales, or even from nonlinear scales to collapsed and virialized groups and clusters. This power law is difficult to achieve in the standard cosmological framework without delicate fine-tuning (Berlind & Weinberg 2002). SDSS eventually detected a statistically significant deviation from a power law (Zehavi et al. 2004(Zehavi et al. , 2005; today, the power-law correlation function is thought of as a cosmic coincidence, arising in a limited luminosity and redshift range that happened to coincide with the scope of former low-z galaxy surveys (Watson et al. 2011). We suggest that a similar fate may lie in store for the PPSD.