text:
# DARK-HALO CUSP: ASYMPTOTIC CONVERGENCE

## 1 INTRODUCTION

A relatively robust universal shape for the density profile of dark-matter halos has been seen in cosmological N-body simulations of dissipationless hierarchical clustering from Gaussian initial fluctuations. It can be approximated by the functional form

$$\rho (r)=\rho _\mathrm{s}\left(\frac{r}{r_\mathrm{s}}\right)^{\alpha _{\mathrm{in}}}\left(1+\frac{r}{r_\mathrm{s}}\right)^{\alpha _{\mathrm{in}}\alpha _{\mathrm{out}}},$$
(1)
where $`r_\mathrm{s}`$ and $`\rho _\mathrm{s}`$ are characteristic inner radius and density respectively.111a useful generalized functional form has been proposed by Zhao (1996) This density profile is characterized by an inner “cusp” $`r^{\alpha _{\mathrm{in}}}`$ and a continuous steepening through a bend near $`r_\mathrm{s}`$ towards $`r^3`$ near the virial radius $`R_\mathrm{v}`$ (defined by a fixed mean overdensity $`\mathrm{\Delta }_\mathrm{v}`$ above the universal mean, with $`\mathrm{\Delta }_\mathrm{v}=180`$ to 340, depending on time and the cosmological model). Navarro, Frenk & White (1995; 1996; 1997, hereafter NFW) found eq. (1) with $`\alpha _{\mathrm{in}}1`$ and $`\alpha _{\mathrm{out}}3`$ to be a good fit to halos in simulations over the radius range $`(0.011)R_\mathrm{v}`$, for a wide range of halo masses and for a range of hierarchical cosmological scenarios with different power spectra of initial fluctuations. Cole & Lacey (1996) came to a similar conclusion for self-similar scenarios with power-law power spectra, $`P_kk^n`$ with $`n=0,1,2`$, in an Einstein-deSitter cosmology. High-resolution simulations of a few individual halos in a cosmological environment (Moore et al. 1998; Ghigna et al. 2000; Klypin et al. 2001) found that the typical asymptotic cusp profile at $`rr_\mathrm{s}`$ is sometimes somewhat steeper, closer to $`\alpha _{\mathrm{in}}1.5`$. A careful convergence analysis by Power et al. (2002), who explored the robustness to numerical errors, found for the standard $`\mathrm{\Lambda }`$CDM cosmology that $`\alpha _{\mathrm{in}}`$ reaches a slope shallower than 1.2 at their innermost resolved point of $`r0.005R_\mathrm{v}`$. Thus, the robust result of the simulations is that dark halos have inner cusps with a characteristic slope $`1\alpha _{\mathrm{in}}1.5`$. We seek a basic theoretical understanding of the origin of these cusps.

An even more intriguing puzzle is introduced by observations which indicate that at least in some cases the actual inner halo density profiles are close to flat cores, with $`\alpha _{\mathrm{in}}0`$. This has been detected directly by rotation curves in low surface brightness (LSB) galaxies, whose centers are dominated by their dark halos (van den Bosch et al. 2000; de Block et al. 2001; Marchesini et al. 2002). Cores have also been argued to exist in normal disk galaxies, where more involved modeling is required (Salucci & Burkert 2000; Salucci 2001; Borriello & Salucci 2001). The presence of a core seems to introduce a severe challenge to the CDM cosmological paradigm. We mention in the conclusion section several mechanisms that have been proposed for the origin of such cores. In particular, attempts to turn a cusp into a core by direct stellar feedback effects in the present halos, which looked promising at a first sight (Navarro, Eke & Frenk 1996), seem not to work (e.g., Geyer & Burkert 2001; Gnedin & Zhao 2002).

We find it useful to first try to understand in simple basic terms the possible origin of the universal cusp in the gravitational N-body simulations of cold dark matter. This may provide us with a tool for addressing the formation and survival of flat cores by other mechanisms, in particular by baryonic feedback processes within the hierarchical CDM framework.

Several mechanisms have been studied in the context of dark-halo profiles in dissipationless simulations. The outer slope of $`r^3`$ (and steeper) may possibly be explained in terms of violent relaxation (e.g., Barnes & Hernquist 1991; Pearce, Thomas & Couchman 1993 and references therein). We note in general that any finite system would tend to have a steep density fall off at large radii due to diffusion of particles outwards. Secondary spherical infall is expected to produce a profile closer to $`\rho r^2`$, which may explain the behavior in the intermediate regions of the halo, but is too steep to explain the flatter inner cusp (Lokas & Hoffman 2000 and references therein). Thus, none of the above mechanisms seem to provide a natural explanation for the characteristic cusp of $`\alpha \stackrel{>}{}1`$.

The clear impression from the cosmological N-body simulations of hierarchical clustering scenarios is that halos are largely built up by a sequence of mergers of smaller structures.222though cusps may also emerge in other scenarions where mergers are suppressed, see a discussion in §8. In a typical merger, a bound satellite halo spirals into the center of the larger halo due to gravity and dynamical friction. The satellite gradually transfers mass into the host halo due to tidal stripping or by eventually melting into the halo inner region. This process is likely to have an important effect in shaping up the density profile. Indeed, Syer & White (1998), Nusser & Sheth (1999) and Subramanian, Cen & Ostriker (2000, hereafter SCO) argued, using certain simple models and simulations of a sequence of mergers, that the buildup by mergers may naturally lead to a stable profile. However, they find their predicted profile to be quite sensitive to the power spectrum of fluctuations and to allow an inner slope of $`\alpha <1`$, in conflict with the robust result of the cosmological simulations. In fact, when trying to repeat the Syer & White analysis using their simplified modeling of the stripping process but with higher resolution, we find that in the long run the profile does not really converge to a stable cusp but rather continues to steepen slowly towards $`\alpha _{\mathrm{in}}=3`$ (see §8). Either way, it seems that something is not adequate in the simplified model adopted to describe the mass transfer from the satellite to the halo.

We re-visit the buildup of halo profile by merging satellites and gain an encouraging new insight. We add two important new ingredients to the tidal effects. In another paper (Dekel, Devor & Hetzroni 2003, hereafter DDH) we argue that for a flat mean density profile with $`\alpha 1`$ the tidal effects on typical satellites induce three-dimensional compression with no local mass deposit, which results in a rapid steepening of the inner profile to $`\alpha >1`$. In this paper, we derive a useful prescription for tidal mass transfer at $`\alpha >1`$, and obtain higher deposit efficiency at higher $`\alpha `$. We then show that this tends to flatten steep profiles with large $`\alpha `$ and thus slowly leads to an asymptotic fixed point at a certain $`\alpha =\alpha _{\mathrm{as}}\stackrel{>}{}1`$.

In §2 we address the compression at $`\alpha <1`$ and the resultant steepening (DDH). In §3 we derive a simple mass-transfer prescription in the range $`\alpha >1`$, based on merger N-body simulations and toy-model understanding. In §4 we explain why this prescription should lead to an asymptotic cusp. In §5 we use linear perturbation analysis to compute the asymptotic slope for satellites of a given mass, and in §6 we extend the analysis to a cosmological distribution of masses. In §7 we demonstrate this process via a semi-analytic simulation of a cosmological sequence of mergers. In §8 we discuss our results.

## 2 CORE TO CUSP BY TIDAL COMPRESSION

This introductory section describes the robust transition of a core ($`\alpha _{\mathrm{in}}1`$) to a cusp ($`\alpha _{\mathrm{in}}>1`$), which is the main theme of another paper, DDH. We summarize it here for completeness and as a background to the independent analysis described in the following sections.

A useful quantity in describing the tidal forces exerted by a halo of mass profile $`M(r)`$ is the local logarithmic slope of its mean density profile $`\overline{\rho }(r)M(r)/r^3`$,

$$\alpha (r)\frac{d\mathrm{ln}\overline{\rho }}{d\mathrm{ln}r},$$
(2)
such that locally $`\overline{\rho }r^\alpha `$. We assume that $`\alpha `$ is either constant or monotonically increasing as a function of $`r`$, with values in the range $`0\alpha 3`$. The extreme values of $`\alpha =0`$ and 3 correspond to a constant-density halo and a point mass respectively. Note that if the profile inside $`r`$ is a power law, then the local and mean density profiles have the same logarithmic slope. In general, they are related via $`\rho (r)=[1\alpha (r)/3]\overline{\rho }(r)`$, but they do not necessarily have the same slope at a given $`r`$. The slope of $`\rho (r)`$ is equal to or larger than the slope of $`\overline{\rho }(r)`$. The following analysis refers to $`\alpha `$ as the slope of $`\overline{\rho }(r)`$.

We then consider a satellite of mass $`m_\mathrm{v}M_\mathrm{v}`$, moving under the gravity excreted by the halo, when its center of mass is at position $`𝐫`$ as measured from the halo center. The tidal acceleration exerted by the halo mass distribution on a satellite particle at position vector $`\mathrm{}`$$`\mathrm{}`$$`\mathrm{}`$ relative to the satellite center of mass is obtained by transforming the gravitational attraction exerted by the halo on the particle into the (non-rotating) rest frame of the accelerated satellite. In the tidal limit $`\mathrm{}r`$, to first order in $`\mathrm{}/r`$, using Cartesian coordinates about the satellite center, where $`\mathrm{}\mathrm{}\mathrm{}=(\mathrm{}_1,\mathrm{}_2,\mathrm{}_3)`$ and $`\mathrm{}\mathrm{}\mathrm{}_1`$ lies along $`𝐫`$, we obtain

$$𝐅_\mathrm{t}=\frac{GM(r)}{r^3}\left([\alpha (r)1]\mathrm{}\mathrm{}\mathrm{}_1\mathrm{}\mathrm{}\mathrm{}_2\mathrm{}\mathrm{}\mathrm{}_3\right).$$
(3)
The components perpendicular to the line connecting the centers of mass are always of compression towards the satellite center, with an amplitude that does not explicitly depend on $`\alpha `$. The maximum radial tidal force outwards is obtained along the line connecting the centers of mass. In the limit where the tides are exerted by a point-mass halo, $`\alpha =3`$, the pull outwards is maximal. For flatter halo slopes, the tidal stretching becomes weaker in proportion to $`(\alpha 1)`$, until it vanishes at $`\alpha =1`$ and reverses direction into compression for $`\alpha <1`$. Thus, while for $`\alpha >1`$ there is always a tidal component pulling outwards, for $`\alpha <1`$ the tidal forces are of compression everywhere in the satellite, and in the limit of a core with $`\alpha =0`$ the tides induce symmetric compression in all directions.

We argue in DDH that this critical transition at $`\alpha 1`$ is the main source for the origin of a cusp steeper than $`r^1`$. The idea is that if the local tidal mass transfer from the satellite to the halo stops when the satellite’s orbit has decayed into a core region where $`\alpha (r)1`$, the satellite would continue to sink in due to dynamical friction without further mass loss until it settles in the halo center. This would inevitably cause a general steepening of the core profile towards $`\alpha >1`$. In DDH we show in different ways that, indeed, no tidal transfer of mass from the satellite to the halo is expected in a region where $`\alpha 1`$. We address this point in simple cases using analytic approximations in the impulse and adiabatic limits and then demonstrate the anticipated effects using merger N-body simulations.

## 3 MASS-TRANSFER PRESCRIPTION AT $`\alpha >1`$$`\alpha >1`$$`\alpha >1`$

Once satellites continue to merge with a halo of an inner cusp $`\alpha >1`$, the final slope is determined by the details of the tidal mass transfer from the satellite to the halo. We analyze this process via a quite general toy model prescription, which we justify and calibrate using N-body simulations.

### 3.1 Toy Model Prescription

We consider a halo of mean density profile $`\overline{\rho }(r)M(r)/r^3`$ with a corresponding logarithmic slope profile $`\alpha (r)`$, and a merging satellite of initial mass profile $`m(\mathrm{})`$ with a corresponding mean density profile $`\overline{\sigma }(\mathrm{})m(\mathrm{})/\mathrm{}^3`$. Let $`m_\mathrm{f}(r)`$ describe the final distribution of stripped satellite mass in spheres about the halo center. By equating $`m_\mathrm{f}(r)`$ and $`m(\mathrm{})`$ we obtain a one to one correspondence between $`\mathrm{}`$ and $`r`$. We then define the initial-density ratio

$$\psi [\alpha (r)]=\frac{\overline{\rho }(r)}{\overline{\sigma }[\mathrm{}(r)]}.$$
(4)
This is an operational definition independent of any model assumption. We argue below that $`\psi (\alpha )`$ is a relatively robust function quite insensitive to the specific nature of the merger, with the distinct property that it is a monotonically decreasing function of $`\alpha `$, and with values significantly below unity in the range $`\alpha >1`$. The robustness of $`\psi (\alpha )`$ is indicated in §3.2 by a preliminary set of merger simulations, to be confirmed and refined by a more complete suite of simulations in an associated paper.

In order to obtain a qualitative feeling for why $`\psi (\alpha )`$ should decrease with $`\alpha `$ we appeal to a very simplified toy model. We assume that the halo is spherical and fixed during the merger, and that the satellite is spherical and is tidally stripped shell by shell outside a momentary tidal radius. The stripped material from a given satellite shell, identified by its original radius $`\mathrm{}`$, is assumed to be deposited on average in a halo radius $`r`$.

The condition $`\psi =1`$ refers to a crude resonance condition, $`\overline{\sigma }(\mathrm{})=\overline{\rho }(r)`$, where the orbital period at $`\mathrm{}`$ within the (unperturbed) satellite equals the orbital period of the satellite within the halo at $`r`$. This is commonly assumed to approximate the momentary tidal radius (e.g., Syer & White 1998; Klypin 1999a). However, this condition ignores the inevitable $`\alpha `$ dependence of the stripping process, and it fails to address the structural changes of the satellite before stripping and the difference between where the stripping occurs and where the mass is actually deposited.

We demonstrated in DDH that the tidal mass transfer becomes weak at $`\alpha <1`$, which means that $`\psi `$ is higher there compared to where $`\alpha >1`$. If we define the tidal radius $`\mathrm{}_\mathrm{t}`$ by the Lagrangian point where the net force in the satellite rest frame vanishes along the line connecting the centers of mass of halo and satellite, we obtain for a satellite on a circular orbit

$$\alpha (r)\overline{\rho }(r)=\overline{\sigma }(\mathrm{}_\mathrm{t})$$
(5)
and for a satellite on a radial orbit

$$[\alpha (r)1]\overline{\rho }(r)=\overline{\sigma }(\mathrm{}_\mathrm{t}).$$
(6)
If the satellite structure inside the tidal radius is assumed to remain fixed, and if the stripped material outside $`\mathrm{}_\mathrm{t}`$ is assumed to be effectively deposited at the $`r`$ where it is stripped, we obtain $`\psi (\alpha )\alpha ^1`$ and $`(\alpha 1)^1`$ respectively. This is the kind of decrease in $`\psi (\alpha )`$ expected due to the $`\alpha `$ dependence of the stripping efficiency.

Another source of $`\alpha `$ dependence is due to the systematic difference between the radii of stripping and deposit. Along a typical orbit, the satellite distance from the halo center oscillates periodically between the radii of apocenter $`r_\mathrm{a}`$ and pericenter $`r_\mathrm{p}`$, while their amplitudes gradually decay due to dynamical friction. Ghigna et al. (1998) studied the distribution of satellite orbits in a high-resolution N-body simulation of a cluster emerging from a CDM cosmological background and found that the median ratio $`r_\mathrm{a}/r_\mathrm{p}`$ is 6:1, with radial orbits common and circular orbits rare, and with a distribution of eccentricities quite independent of $`r`$. They also demonstrated that the tidal radii of the satellites are consistent with being determined near pericenter, approximately under the general resonance condition $`\overline{\rho }(r_\mathrm{p})\overline{\sigma }(\mathrm{}_\mathrm{t})`$. Particles that escape from the satellite can be assumed, on average, to continue on an orbit about the halo center with apocenter and pericenter radii “frozen” at their values near the time of escape, suffering no further decay due to dynamical friction. Since a particle spends most of its time near the apocenter of its orbit, we can say that the escapers are effectively “deposited” near the apocenter radius valid at the time of escape. We can thus assume that the ratio of stripping radius ($`\stackrel{>}{}r_\mathrm{p}`$) and deposit radius ($`\stackrel{<}{}r_\mathrm{a}`$) is typically $`ϵ\stackrel{>}{}r_\mathrm{p}/r_\mathrm{a}1/6`$. If we ignore the structural evolution of the satellite inside the tidal radius as well as the $`\alpha `$ dependence of the stripping efficiency addressed above, we obtain straightforwardly a decrease as a function of $`\alpha `$ of the sort $`\psi (\alpha )=ϵ^\alpha `$, where $`\alpha `$ is the effective slope of the mean density profile between $`r_\mathrm{p}`$ and $`r_\mathrm{a}`$.

A third source of $`\alpha `$ dependence is the distortion of the satellite inside the tidal radius, with an $`\alpha `$-dependent stretching along one axis and an $`\alpha `$-independent compression along the others. The resultant decrease in satellite density before stripping has been seen in simulations (e.g. Klypin et al. 1999a, Fig. 6; Hayashi et al. 2002). At higher $`\alpha `$ values the stretching is stronger, which should make the stripping even more efficient there.

The above three effects provide qualitative hints for the expected properties of $`\psi (\alpha )`$, and in particular for it being a decreasing function of $`\alpha `$ with values below unity. We turn to N-body simulations in order to quantify these predictions.

### 3.2 Merger Simulations

In order to evaluate $`\psi (\alpha )`$ at $`\alpha >1`$ and test its robust properties, we ran several N-body simulations of isolated mergers between a large halo and a satellite halo of mass ratio $`m/M=0.1`$. Similar simulations were used in DDH to study the inner core region and they are also described there using several additional figures which are not essential for our current purpose. We use the Tree code by Mihos & Hernquist (1996 and references therein) but with dark-matter halos only (no gaseous disks). In our default simulations, the large host halo is represented by $`10^5`$ equal-mass particles and the satellite by $`10^4`$ particles. We re-ran one case with $`0.55\times 10^6`$ particles and the force resolution higher by a factor $`\sqrt{5}`$, and confirmed that the results were practically identical. The simulation units are: length $`3.5`$kpc, mass $`5.6\times 10^{10}M_{}`$, and time $`13.06`$Myr. The force softening length is 0.08 units, i.e. $`0.28`$kpc

The initial halo density profile, as measured in the unperturbed initial conditions, is a truncated isothermal sphere with a flat core,

$$\rho (r)=\frac{\rho _\mathrm{s}e^{(r/r_\mathrm{t})^2}}{1+(r/r_\mathrm{s})^2},$$
(7)
with $`\rho _\mathrm{s}=10.36`$, $`r_\mathrm{s}=1`$ and $`r_\mathrm{t}=10`$ The density at the characteristic radius $`r_\mathrm{s}`$ is thus $`\rho (r_\mathrm{s})=5.13`$. The internal velocities are constructed to fulfill the isotropic Jeans equation which ensures an equilibrium configuration as discussed in Mihos & Hernquist (1996). When run in isolation, the halo profile has been tested to be very stable for many dynamical times. The logarithmic slope $`\alpha (r)`$ of the mean density profile spans the range of interest between $`\alpha =0`$ and 3, and its variation as a function of radius can be described to a good approximation by $`\alpha (r)1.73\mathrm{log}r+0.67`$ throughout the range $`0.3\alpha 2.9`$. The initial halo density profile can be seen in Fig. 3 below. It resembles the generalized NFW profile of eq. (1) with a core of $`\alpha 1`$. We chose to start with an inner core in order to also simulate the rapid steepening to a cusp, but the outer profile, where $`\alpha >1`$, is quite generic.

The satellite initial density profile is fit by a Hernquist profile,

$$\sigma (\mathrm{})=\frac{\sigma _\mathrm{s}}{(\mathrm{}/\mathrm{}_\mathrm{s})[1+(\mathrm{}/\mathrm{}_\mathrm{s})]^3},$$
(8)
with the default choice $`\sigma _\mathrm{s}=19.2`$ and $`\mathrm{}_\mathrm{s}=1`$ defining our typical “compact” satellite. The initial satellite density profile can also be seen in Fig. 3. In the inner region it convergence to the NFW profile, $`\alpha _{\mathrm{in}}=1`$. If we fit this Hernquist profile with an NFW profile by matching the characteristic radii where the local logarithmic slope is $`2`$, we find that the radius corresponding to the characteristic radius of the NFW profile is $`\mathrm{}=\mathrm{}_\mathrm{s}/2`$. At this radius the density is $`16\sigma _\mathrm{s}/27`$, which for our default satellite is $`\sigma (\mathrm{}_\mathrm{s}/2)=11.38`$.

The satellite parameters were chosen to roughly mimic a typical compact satellite according to the distribution of halo properties in the $`\mathrm{\Lambda }`$CDM scenario. In a hierarchical clustering scenario with a fluctuation power spectrum $`P_kk^n`$, the halo characteristic radii and densities scale like $`\mathrm{}_\mathrm{s}/r_\mathrm{s}m^{(1+\nu )/3}`$ and $`\sigma _\mathrm{s}/\rho _\mathrm{s}m^\nu `$, with $`\nu =(3+n)/2`$. For $`\mathrm{\Lambda }`$CDM on galactic scales one has $`\nu 0.33`$. This is confirmed by cosmological simulations (NFW; Bullock et al. 2001a), where the typical halo profile is NFW. The halos of lower masses are thus typically more compact; they tend to have lower characteristic radii and higher corresponding mean densities within these radii (corresponding to higher virial concentration parameters). For CDM halos of mass ratio $`m/M=1/10`$, the typical ratio of characteristic radii and corresponding densities are expected to be roughly $`0.4`$ and $`2.1`$. The corresponding ratios in the initial conditions of our simulations are approximately $`0.5`$ and $`2.2`$, providing a reasonable match. It may also be interesting to note that the satellite mean density interior to $`\mathrm{}_\mathrm{s}/2`$ (or $`\mathrm{}_\mathrm{s}`$) is about $`3.9`$ (or $`1.09`$) times the halo mean density interior to $`r_\mathrm{s}`$.

We simulated three cases of initial merger orbits: a radial orbit, a circular orbit, and a typical elongated orbit with $`r_\mathrm{p}/r_\mathrm{a}1/6`$ initially. The unperturbed satellite is put initially at $`r=20`$, i.e. at about $`2r_\mathrm{t}`$, where the initial circular period is about $`230`$Myr. In the circular and radial cases the magnitude of the initial satellite velocity was set to equal the circular velocity of the halo at that radius, namely a bound orbit with an orbital kinetic energy that equals half the absolute value of the total energy. For the elongated orbit the initial tangential velocity was one half of the circular velocity at $`r=20`$. Each merger has been followed until the satellite’s bound core has practically settled at the halo center.

Fig. 1 shows the satellite mass distribution in 6 snapshots during the elongated merger, projected onto the orbital plane (while DDH show analogous plots for the radial and circular mergers). The satellite oscillates about the halo center through a sequence of pericenter and apocenter passages. The shape of the orbit remains roughly constant while it shrinks in scale; when measuring the ratio of pericenter to the following apocenter it is $`r_\mathrm{p}/r_\mathrm{a}1/3.5`$, and when measuring the ratio of pericenter to the preceding apocenter it is $`r_\mathrm{p}/r_\mathrm{a}1/3.5`$, but these ratios remain roughly the same for all detectable pericenters. The oscillations decay due to dynamical friction until the satellite becomes confined to the halo core after $`125`$Myr and 5 pericenter passages. By the first pericenter, the satellite is already stretched and stripped along its orbit, while it is temporarily shrunk in the perpendicular direction (as expected in DDH, see §2). This is followed by a re-bounce and significant mass loss about the following apocenter. The visual impression confirms the notion that the particles that are torn away near a pericenter radius continue on orbits that reflect on average the satellite orbit at the time of stripping, while the bound remnant continues to sink into smaller radii due to dynamical friction. For example, in the snapshots corresponding to the second and third apocenters we clearly see a large amount of stripped satellite material spread about the location of the previous apocenter. We can crudely say that mass is stripped near pericenter and is practically “deposited” about the following apocenter radius. The final distribution of satellite mass extends quite smoothly about the halo center in a puffy oblate ellipsoid that looks quite symmetric in the orbital plane.

Fig. 2 describes the time evolution of the satellite spherical mass profile during the elongated merger by showing the mean radii of concentric spherical shells about the satellite maximum-density center, each encompassing a given fraction of the satellite mass (i.e., not necessarily the same population of particles at different times). Pericenter passages are identified five times. Overall contraction of bound satellite shells seems to start roughly at these times, as expected both in the impulse and adiabatic limits (see DDH). Each major contraction is followed by a re-bounce as the satellite moves towards apocenter, which results in overall expansion and stripping of the outer shells. This seems like a manifestation of the expected delayed stripping due to the energy pumped into the satellite by the inwards impulse in the core, as well as the adiabatic stretching and stripping outside the tidal radius. Once the satellite becomes confined to the inner halo where $`\alpha <1`$, after $`125`$Myr and 5 apocenters, there is no apparent overall shell expansion anymore, indicating that stripping has stopped, as expected. Similar effects are noticed in DDH for the radial and circular mergers.

Fig. 3 addresses the cusp formation straightforwardly by showing the density profile of the halo before and after the merger on a radial orbit. In DDH we also show similar plots for the circular and elongated merger orbits, with very similar results. For the standard, compact satellite, the left panel demonstrates the inevitable steepening of the profile in the core region, starting near the core boundary where $`\alpha 1`$. The slight depletion of final density in the region near $`\alpha 1`$, compared to the slight increment in halo density at larger radii, is consistent with no mass transfer in this region while the orbit of the remaining satellite continues to decay into smaller radii.

How stable is this cusp under further mergers with similar satellites? To test this, we performed a series of mergers following each other. The main progenitor halo in each merger is taken to be the outcome of the previous merger, except that its total mass is scaled down back to the original halo mass (such that all the density profiles can be directly compared). We achieve this by letting each simulation involve the same numbers of $`N`$ halo particles and $`n`$ satellite particles (each of a fixed mass $`m`$), where the halo particles are selected at random from the $`N+n`$ particles of the halo produced in the previous merger. In these tandem mergers we used $`N=20,000`$ and $`n=2,000`$. The halo profile in the range of interest is found to remain stable for at least several hundred Myr when the halo is run in isolation. The satellite starts with the same profile as the original compact satellite. It is put at the same distance of $`r=20`$ as before but in an elongated orbit of a random spatial orientation. The initial satellite velocity is tangential to the line connecting the centers of mass of halo and satellite, with an amplitude in units of the circular velocity at $`r=20`$ chosen at random in the range $`v/v_\mathrm{c}=0.20.6`$ (compared to the typical case simulated before, where $`v/v_\mathrm{c}=0.49`$ led to peri/apocenter ratio of roughly 1:6). Each merger was followed for 325Myr before the following merger started. The density profiles after each of the first 5 mergers are shown in the middle panel of Fig. 3. The first merger reproduces a cusp very similar to the cusp produced when the merger was simulated with $`N=100,000`$ (and with $`N=500,000`$), indicating that even $`N=20,000`$ is adequate for crude results in the regime of interest here. We see that inside the initial core radius the cusp is stable at a slope of $`\alpha 1.5`$, with no systematic tendency to deviate from the power-law profile resulting already after the first merger. At the radius where $`\alpha \stackrel{>}{}1.5`$ the profile flattens, and the amplitude becomes lower partly because of the renormalization of the total mass between every two mergers. Our main goal in the current paper is to try to understand the origin of the stable cusp as indicated by these simulations of tandem mergers.

How much puffing is needed in order to prevent damage by merging satellites to the halo core? In order to obtain a first clue, we have performed simulations similar to those described above, except that the compact satellites of CDM have been replaced by more puffy satellites of the same total mass. In the Hernquist profile, eq. (8), where the default compact satellites had $`\sigma _\mathrm{s}=19.2`$ and $`\mathrm{}_\mathrm{s}=1`$, the parameters of the puffy satellites are now $`\sigma _\mathrm{s}^{}=\sigma _\mathrm{s}/8`$ and $`\mathrm{}_\mathrm{s}^{}=2\mathrm{}_\mathrm{s}`$. With this choice, the density at $`\mathrm{}=1/2`$, the characteristic radius of the original compact satellite, is scaled down by a factor of $`2.3`$, to just below the density of the host halo at its characteristic radius. This corresponds to a reduction by a factor of $`2.8`$ in the mean density interior to $`\mathrm{}=1/2`$. The right panel of Fig. 3 shows the effect of such a merger with a puffy satellite on the halo density profile in the case of a radial radial merger (DDH also show similar results for other merger orbits). When the satellite is puffy, we see that almost all the satellite mass is stripped before the satellite orbit decays to the $`\alpha 1`$ zone, and as a result the halo core is practically unaffected. We learn that a modest reduction in the initial satellite inner density is enough for preventing the cusp formation seen in the compact-satellite case.

### 3.3 Measuring $`\psi (\alpha )`$$`\psi (\alpha )`$$`\psi (\alpha )`$ in the Simulations

The measurement of $`\psi (\alpha )`$ in the simulations is straightforward (and free of any model assumption). We simply measure the mass profile of the stripped satellite mass about the halo center at the final time, $`m_\mathrm{f}(r)`$, and equate it with the initial satellite mass profile $`m(\mathrm{})`$ to obtain the deposit relation $`\mathrm{}(r)`$. Then $`\psi `$ is evaluated from eq. (4) at any desired $`r`$, and is expressed as $`\psi (\alpha )`$ given the halo slope profile $`\alpha (r)`$.

Fig. 4 shows the measured $`\psi (\alpha )`$ in three different merger simulations, where the orbits are very different from each other: circular, radial and elongated 1:6. All the curves show a similar general behavior, consistent with our qualitative expectations. First, the value of $`\psi `$ is significantly below unity, demonstrating that the mass transfer is more efficient than predicted by the naive $`\psi =1`$ model. Second, the function $`\psi (\alpha )`$ is monotonically decreasing with increasing $`\alpha `$ in the range $`1<\alpha <3`$. The general behavior can be crudely approximated by several different functions, such as the power law $`\psi (\alpha )0.5/\alpha `$, or the exponential $`\psi (\alpha )0.5^\alpha `$. The scatter about these fits is roughly $`\pm 0.1`$. The relative robustness of $`\psi (\alpha )`$ to the nature of the merger orbit indicates that eq. (4) may serve as a useful approximate recipe for tidal mass transfer in a general merger.

In order to evaluate the sensitivity of our results to the resolution of the simulation, we ran a case identical to the merger on an elongated orbit described above, but now with 5 times more particles ($`N=0.55\times 10^6`$) and a softening length smaller by a factor of $`5^{1/3}`$ accordingly. The results are found to be practically identical in all respects, namely the decay rate of the satellite radius within the halo, the mass loss from the satellite, the final halo density profile, and the derived $`\psi (\alpha )`$. This indicates that the resolution of our simulations is adequate in the range of radii of relevance, well inside the core/cusp and down to below $`0.2r_\mathrm{s}`$.

In the current simulations we have only tested in a limited way the relative robustness to the merger orbit. The simulated halo and satellite are of a typical mass ratio and their profiles relate to each other in general accordance with the expected average scaling in the $`\mathrm{\Lambda }`$CDM cosmology, which makes the obtained $`\psi (\alpha )`$ a sensible first guess. However, the robustness of $`\psi (\alpha )`$ to variations about this mass ratio and these profiles should be confirmed and refined using a more complete suite of merger simulations. We adopt below $`\psi (\alpha )=0.5/\alpha `$ as an illustrative example for obtaining numerical results but the following analysis is valid for a quite general $`\psi (\alpha )`$.

## 4 AN ASYMPTOTIC PROFILE: QUALITATIVE

We saw in DDH that when the halo profile is flat, $`\alpha 1`$, tidal compression causes rapid steepening to $`\alpha >1`$. In this section we investigate the development of the profile due to a sequence of mergers between similar halos where $`\alpha `$ is of order unity or larger. We show that if the tidal stripping is described by a condition similar to eq. (4), with $`\psi (\alpha )`$ monotonically decreasing rapidly enough, then the profile evolves slowly towards an asymptotic stable power law $`r^{\alpha _{\mathrm{as}}}`$, with $`\alpha _{\mathrm{as}}`$ larger than unity, which we can evaluate once $`\psi (\alpha )`$ is given. The evolution to this asymptotic profile is through a sequence of profiles which crudely resemble the generalized NFW shape, and with an inner cusp of slope slightly larger than unity, $`1<\alpha <\alpha _{\mathrm{as}}`$. While the value of $`\alpha _{\mathrm{as}}`$ may depend on the exact shape of $`\psi (\alpha )`$, the following analysis showing the convergence to an asymptotic profile is not sensitive to it as long as $`\psi (\alpha )`$ is decreasing rapidly enough.

### 4.1 Homologous Halos

We consider the buildup of the halo of mass $`M`$ by a sequence of mergers with satellites of masses $`mM`$ drawn from a cosmological distribution of halos (similar to Syre & White 1998). In our search for a self-similar evolution, we assume that the halo and satellite before the encounter are homologous, meaning that their unperturbed mean density profiles are scaled version of each other,

$$\overline{\sigma }(\mathrm{})=\varrho \overline{\rho }(\mathrm{}/\lambda ),\varrho =m/\lambda ^3,$$
(9)
where $`m`$ stands for the satellite-halo mass ratio $`m/M`$ (we use hereafter $`M=1`$). If we refer, for example, to the functional form of eq. (1), then $`\varrho =\sigma _\mathrm{s}/\rho _\mathrm{s}`$ and $`\lambda =\mathrm{}_\mathrm{s}/r_\mathrm{s}`$. The cosmological N-body simulations of CDM as well as of power-law power spectra (NFW; Bullock et al. 2001a) show that halos of lower masses tend to have lower characteristic radii and higher corresponding mean densities within these radii (corresponding to higher virial concentration parameters). The scaling of the averages of these quantities, at any given time, is found to be

$$\lambda =m^{(1+\nu )/3},\varrho =m^\nu ,$$
(10)
with $`\nu 0.33`$ for $`\mathrm{\Lambda }`$CDM. Given a power-law power spectrum $`P(k)k^n`$, a simple scaling argument based on linear theory predicts this behavior with $`\nu =(3+n)/2`$. Thus, $`\nu 0.33`$ corresponds to $`n2.3`$, as expected for a $`\mathrm{\Lambda }`$CDM power spectrum at the relevant scales.

### 4.2 Convergence to an Asymptotic Profile

The origin of an asymptotic slope can be qualitatively understood in simple terms via the illustration in Fig. 5. We show a schematic mean density profile of a halo and a homologous satellite scaled accordingly. As in the toy-model interpretation of eq. (4), we assume that satellite shell $`\mathrm{}`$ is being deposited at halo radius $`r`$ where the slope is $`\alpha `$, with a given decreasing function $`\psi (\alpha )`$. To visualize this deposit process, we connect by arrows two points on the log-log curve $`\overline{\sigma }(\mathrm{})`$ with the corresponding points on the log-log curve $`\overline{\rho }(r)`$, illustrating both the (horizontal) distance ratio $`r/\mathrm{}`$ and the (vertical) density ratio $`\psi `$. The simplified condition $`\psi =1`$ (as adopted, e.g., by Syre & White 1998) leads to a correspondence via horizontal arrows. In this case, the satellite slope at $`\mathrm{}`$ is always steeper than the halo slope at $`r`$, naturally leading to a steepening of the halo profile. The mean profile evolves in this case towards the steepest possible power law, $`\alpha =3`$. On the other hand, when $`\psi `$ is properly decreasing with $`\alpha `$, the relation between $`\mathrm{}`$ and $`r`$ is now illustrated by the solid arrows, where the vertical displacement is a growing function of $`r`$ and $`\alpha `$. This is expected to weaken the steepening and possibly turn it into flattening of the profile, because the slope at $`\mathrm{}`$ is now closer to the slope at $`r`$, and may be even flatter. Another way to understand this flattening effect is by noticing that due to the decrease of $`\psi (\alpha )`$ more satellite mass is deposited at outer halo radii.

We can be a bit more quantitative as follows. With $`mM`$, we assume that the original halo is not affected by tides, and that its mean density profile evolves only by the addition of the satellite mass wherever it is being deposited, namely,

$$\overline{\rho }_{\mathrm{final}}(r)=\overline{\rho }(r)+\overline{\sigma }[\mathrm{}(r)]\frac{\mathrm{}(r)^3}{r^3},$$
(11)
where $`\mathrm{}(r)`$ is defined by $`m_\mathrm{f}(r)=m(\mathrm{})`$. Given that $`\alpha (r)`$ is minus the logarithmic derivative of $`\overline{\rho }(r)`$, it is straightforward to show that the change of $`\alpha `$ at $`r`$ due to the merger is

$$\mathrm{\Delta }\alpha (r)=\frac{F(r)}{1+F(r)}\frac{\mathrm{d}F(r)}{\mathrm{d}r},F(r)\frac{\overline{\sigma }(\mathrm{})}{\overline{\rho }(r)}\frac{\mathrm{}^3}{r^3}.$$
(12)
Note that the first term in $`F(r)`$ is $`1/\psi (r)`$, and that $`F(r)=m(\mathrm{})/M(r)`$.

Every power law is a self-similar solution, where $`\mathrm{\Delta }\alpha (r)=0`$. This is because when $`\alpha (r)`$ is constant, so is $`\psi (r)`$, namely $`\overline{\sigma }(\mathrm{})/\overline{\rho }(r)`$ is constant, implying that $`\mathrm{}/r=\mathrm{const}.`$, and therefore the derivative $`\mathrm{d}F/\mathrm{d}r`$ vanishes. However, $`\mathrm{\Delta }\alpha (r)=0`$ does not guarantee a stable solution to which the process would converge.

When $`\alpha `$ is varying with $`r`$, as long as $`\alpha `$ is small and $`\psi `$ is close to unity, the ratio $`\mathrm{}/r`$ is a decreasing function of $`r`$. Thus, if $`\psi [\alpha (r)]=`$const., the $`\mathrm{}^3/r^3`$ term determines a positive $`\mathrm{\Delta }\alpha (r)`$, namely continuous steepening towards $`\alpha =3`$. On the other hand, when $`\psi (\alpha )`$ is properly decreasing with $`\alpha `$, $`1/\psi [\alpha (r)]`$ is an increasing function of $`r`$, an increase which can balance the decrease of $`\mathrm{}^3/r^3`$ and produce a fixed point of $`\mathrm{\Delta }\alpha =0`$ at some asymptotic slope $`\alpha _{\mathrm{as}}`$. This is a stable solution, where $`\mathrm{\Delta }\alpha >0`$ at $`\alpha <\alpha _{\mathrm{as}}`$, and $`\mathrm{\Delta }\alpha <0`$ at $`\alpha >\alpha _{\mathrm{as}}`$.

### 4.3 Linear Perturbation Analysis: Summary

A rigorous perturbation analysis, for a given $`\psi (\alpha )`$ and $`m/M`$, is described in the following section. We provide here a summary of this analysis and its main result so that the reader who is not interested in the mathematical details can skip §5, which is rather technical.

The linear analysis is based on the following two assumptions, which are argued to be of quite general validity:

* Near the asymptotic solution the slope $`\alpha (r)`$ is varying slowly, namely, the first and higher derivatives of $`\alpha (r)`$ are much smaller than $`\alpha (r)`$. This allows us to write $`\alpha (r)\alpha _0+ϵ\delta (r)`$, and examine the behavior of $`\mathrm{\Delta }\alpha `$ at first order in the small parameter $`ϵ`$.
* The (very small) varying part of $`\alpha (r)`$ is approximately linear in $`\mathrm{ln}r`$, namely, the second and higher derivatives of $`\alpha (\mathrm{ln}r)`$ are much smaller than $`\mathrm{d}\alpha /\mathrm{d}\mathrm{ln}r`$.

This analysis leads to an explicit expression for $`\mathrm{\Delta }\alpha `$ at $`r`$:

$`\mathrm{\Delta }\alpha (r)`$ $``$ $`{\displaystyle \frac{\alpha ^{}(r)r}{1+\lambda ^3\varrho ^{3/\alpha }\psi (\alpha )^{13/\alpha }}}`$
$`\times {\displaystyle \frac{1}{\alpha ^2}}\left(\alpha (\alpha 3){\displaystyle \frac{\psi ^{}(\alpha )}{\psi (\alpha )}}+3\mathrm{ln}[\varrho \psi (\alpha )]\right).`$
Here $`\alpha =\alpha (r)`$, $`\alpha ^{}(r)\mathrm{d}\alpha /\mathrm{d}r`$ and $`\psi ^{}(\alpha )\mathrm{d}\psi /\mathrm{d}\alpha `$. Thus, the sign of $`\mathrm{\Delta }\alpha (r)`$ is the sign of the expression in big brackets, and the asymptotic slope $`\alpha _{\mathrm{as}}`$ is the solution of the equation

$$\alpha (\alpha 3)\frac{\psi ^{}(\alpha )}{\psi (\alpha )}+3\mathrm{ln}[\varrho \psi (\alpha )]=0.$$
(14)
The factor in front of the big parentheses in eq. (4.3) determines the rate of convergence to $`\alpha _{\mathrm{as}}`$; note that the rate slows down as $`\alpha `$ approaches a constant and $`\alpha ^{}0`$. Recall that $`\psi (\alpha )`$ is given, and that the factors $`\varrho `$ and $`\lambda `$ are functions of the mass ratio in the merger and the cosmological power spectrum, eq. (10).

The lines in Fig. 6 show the asymptotic slope $`\alpha _{\mathrm{as}}`$ as predicted by the linear theory, eq. (14), for a sequence of mergers with the same mass ratio $`m/M`$, as a function of $`m/M`$. The assumed mass-transfer recipe in both cases is $`\psi (\alpha )=0.5/\alpha `$. The middle curve is for $`\mathrm{\Lambda }`$CDM, $`n=2.33`$, and the other curves are for $`n=3`$ and $`n=1.67`$.

The discussion so far referred to a sequence of mergers with the same mass ratio. In §6 we will generalize the analysis to the case of a realistic cosmological distribution of mass ratios.

## 5 LINEAR ANALYSIS OF THE ASYMPTOTIC PROFILE

This section is rather technical, and can be skipped by the non-practitioner reader without hurting the flow of the discussion.

### 5.1 The General Case

The assumptions and analysis become more transparent if we switch to logarithmic variables,

$$\stackrel{~}{r}\mathrm{ln}r,\stackrel{~}{\rho }\mathrm{ln}\overline{\rho }.$$
(15)
Hereafter we use a tilde to mark logarithmic quantities in general. Then eq. (11) becomes

$$\stackrel{~}{\rho }_{\mathrm{final}}(\stackrel{~}{r})=\stackrel{~}{\rho }(\stackrel{~}{r})+\mathrm{ln}[1+e^{\stackrel{~}{\sigma }(\stackrel{~}{\mathrm{}})\stackrel{~}{\rho }(\stackrel{~}{r})+3(\stackrel{~}{\mathrm{}}\stackrel{~}{r})}].$$
(16)
Using the log variables, the local slope is simply (minus) the first derivative of the density profile,

$$\alpha (\stackrel{~}{r})=\mathrm{d}\stackrel{~}{\rho }(\stackrel{~}{r})/\mathrm{d}\stackrel{~}{r}.$$
(17)
Then, from eq. (16), the change in the slope in the merger, $`\alpha _{\mathrm{final}}(\stackrel{~}{r})\alpha (\stackrel{~}{r})`$, is

$`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ $`=`$ $`{\displaystyle \frac{1}{1+e^{\stackrel{~}{\rho }(\stackrel{~}{r})\stackrel{~}{\sigma }(\stackrel{~}{\mathrm{}})+3(\stackrel{~}{r}\stackrel{~}{\mathrm{}})}}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\stackrel{~}{r}}}[\stackrel{~}{\sigma }(\stackrel{~}{\mathrm{}})\stackrel{~}{\rho }(\stackrel{~}{r})+3(\stackrel{~}{\mathrm{}}\stackrel{~}{r})]`$ (18)
$`=`$ $`{\displaystyle \frac{1}{1+e^{\stackrel{~}{\rho }(\stackrel{~}{r})\stackrel{~}{\sigma }(\stackrel{~}{\mathrm{}})+3(\stackrel{~}{r}\stackrel{~}{\mathrm{}})}}}\left[\alpha (\stackrel{~}{\mathrm{}}\stackrel{~}{\lambda }){\displaystyle \frac{d\stackrel{~}{\mathrm{}}}{d\stackrel{~}{r}}}\alpha (\stackrel{~}{r})3\left({\displaystyle \frac{d\stackrel{~}{\mathrm{}}}{d\stackrel{~}{r}}}1\right)\right].`$
In the second equation we have assumed that the profiles of the halo and satellite are homologous, eq. (9).

In order to make the above expression for $`\mathrm{\Delta }\alpha `$ useful, we need to know the function $`\stackrel{~}{\mathrm{}}(\stackrel{~}{r})`$, and especially its derivative $`\mathrm{d}\stackrel{~}{\mathrm{}}/\mathrm{d}\stackrel{~}{r}`$. This is provided by the mass-transfer prescription, eq. (4), which in the log variables takes the form

$$\stackrel{~}{\rho }(\stackrel{~}{r})\stackrel{~}{\sigma }(\stackrel{~}{\mathrm{}})=\stackrel{~}{\psi }(\alpha ),$$
(19)
with the function $`\psi (\alpha )`$ given, and $`\alpha =\alpha (\stackrel{~}{r})`$.

Substituting in eq. (19) the scaling relation for homologous halo and satellite, eq. (9), we obtain the equation for mass transfer in the homologous case:

$$\stackrel{~}{\rho }(\stackrel{~}{r})\stackrel{~}{\rho }(\stackrel{~}{\mathrm{}}\stackrel{~}{\lambda })=\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha ).$$
(20)
Once the scaling factors $`\lambda `$ and $`\varrho `$ are given, and the functions $`\psi (\alpha )`$ and $`\alpha (r)`$ are known, this equation should allow us to evaluate the desired $`\stackrel{~}{\mathrm{}}(\stackrel{~}{r})`$, to be used in eq. (18).

### 5.2 First-Order Approximations

For any pure power law, $`\alpha (\stackrel{~}{r})=\mathrm{const}.`$, we have $`\mathrm{\Delta }\alpha (\stackrel{~}{r})=0`$ for all $`\stackrel{~}{r}`$’s. This is because then $`\psi `$ is constant, meaning that $`\stackrel{~}{\rho }\stackrel{~}{\sigma }`$ is constant, so $`\stackrel{~}{\mathrm{}}\stackrel{~}{r}`$ is constant, and then the derivative vanishes in eq. (18). However, not all these solutions are stable. In order to find a stable solution, we examine $`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ for small perturbations of order $`ϵ`$ about $`\alpha =\alpha _0`$, and see how the first order in $`ϵ`$ behaves.

We first wish to evaluate $`\mathrm{d}\stackrel{~}{\mathrm{}}/\mathrm{d}\stackrel{~}{r}`$ at first order in $`ϵ`$. Based on the definition of $`\alpha `$, the homologous mass-transfer condition, eq. (20), can be written as

$$_{\stackrel{~}{r}}^{\stackrel{~}{\mathrm{}}\stackrel{~}{\lambda }}ds\alpha (s)=\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha ).$$
(21)
Into this equation we substitute the perturbed quantities

$`\alpha (\stackrel{~}{r})`$ $`=`$ $`\alpha _0+ϵ\delta (\stackrel{~}{r}),`$ (22)
$`\stackrel{~}{\mathrm{}}(\stackrel{~}{r})`$ $`=`$ $`\stackrel{~}{\mathrm{}}_0(\stackrel{~}{r})+ϵ\stackrel{~}{\mathrm{}}_1(\stackrel{~}{r}),`$ (23)
where $`ϵ`$ is small, $`\delta `$ is of the order of $`\alpha `$, and $`\mathrm{}_1`$ is of the order of $`\mathrm{}`$. We now collect terms of the same order. The zero’s order term yields

$$[\stackrel{~}{\mathrm{}}_0(\stackrel{~}{r})\stackrel{~}{\lambda }\stackrel{~}{r}]\alpha _0=\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha ),$$
(24)
$$\stackrel{~}{\mathrm{}}_0(\stackrel{~}{r})=[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha _0)]/\alpha _0+\stackrel{~}{\lambda }+\stackrel{~}{r}.$$
(25)
The first order term gives

$$_{\stackrel{~}{r}}^{\stackrel{~}{\mathrm{}}_0\stackrel{~}{\lambda }}ds\delta (s)+\stackrel{~}{\mathrm{}}_1(\stackrel{~}{r})\alpha _0=\frac{\psi ^{}(\alpha _0)}{\psi (\alpha _0)}\delta (\stackrel{~}{r}),$$
(26)
$$\stackrel{~}{\mathrm{}}_1(\stackrel{~}{r})=\frac{1}{\alpha _0}\left[\frac{\psi ^{}(\alpha _0)}{\psi (\alpha _0)}\delta (\stackrel{~}{r})_{\stackrel{~}{r}}^{\stackrel{~}{\mathrm{}}_0\stackrel{~}{\lambda }}ds\delta (s)\right],$$
(27)
where $`\psi ^{}(\alpha )\mathrm{d}\psi /\mathrm{d}\alpha `$. Using these results we can now calculate $`\mathrm{d}\stackrel{~}{\mathrm{}}(\stackrel{~}{r})/\mathrm{d}\stackrel{~}{r}`$ to first order:

$$\frac{d\stackrel{~}{\mathrm{}}(\stackrel{~}{r})}{d\stackrel{~}{r}}=1+\frac{ϵ}{\alpha _0}\left[\frac{\psi ^{}(\alpha _0)}{\psi (\alpha _0)}\delta ^{}(\stackrel{~}{r})+\delta (\stackrel{~}{r})\delta (\stackrel{~}{\mathrm{}}_0\stackrel{~}{\lambda })\right],$$
(28)
where $`\delta ^{}(\stackrel{~}{r})\mathrm{d}\delta /\mathrm{d}\mathrm{ln}r`$.

When we plug this result into eq. (18), we obtain after some algebra

$`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ $`=`$ $`{\displaystyle \frac{ϵ}{1+e^{\stackrel{~}{\rho }(\stackrel{~}{r})\stackrel{~}{\sigma }(\stackrel{~}{\mathrm{}}_0)+3(\stackrel{~}{r}\stackrel{~}{\mathrm{}}_0)}}}`$
$`\times \left({\displaystyle \frac{\psi ^{}(\alpha _0)}{\psi (\alpha _0)}}\delta ^{}(\stackrel{~}{r}){\displaystyle \frac{3}{\alpha _0}}\left[{\displaystyle \frac{\psi ^{}(\alpha _0)}{\psi (\alpha _0)}}\delta ^{}(\stackrel{~}{r})+\delta (\stackrel{~}{r})\delta (\stackrel{~}{\mathrm{}}_0\stackrel{~}{\lambda })\right]\right)+𝒪(ϵ^2).`$
We can further simplify the above expression by substituting $`\mathrm{}_0`$ with its explicit value taken from eq. (25). Using eq. (20) and the definition of $`\lambda `$, we obtain

$`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ $`=`$ $`{\displaystyle \frac{ϵ}{1+\lambda ^3\varrho ^{3/\alpha _0}\psi (\alpha _0)^{13/\alpha _0}}}`$
$`\times \left\{{\displaystyle \frac{\psi ^{}(\alpha _0)}{\psi (\alpha _0)}}\delta ^{}(\stackrel{~}{r}){\displaystyle \frac{3}{\alpha _0}}\left[{\displaystyle \frac{\psi ^{}(\alpha _0)}{\psi (\alpha _0)}}\delta ^{}(\stackrel{~}{r})+\delta (\stackrel{~}{r})\delta \left(\stackrel{~}{r}+[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha _0)]/\alpha _0\right)\right]\right\}+𝒪(ϵ^2).`$
Since the above expression is first-order in $`ϵ`$, we may safely replace $`\alpha _0`$ by $`\alpha (\stackrel{~}{r})`$ because the difference between the two is also first order, eq. (22). Additionally, eq. (22) implies the identities

$`ϵ\delta ^{}(\stackrel{~}{r})`$ $`=`$ $`\alpha ^{}(\stackrel{~}{r}),`$ (31)
$`ϵ\delta (\stackrel{~}{r})ϵ\delta \left(\stackrel{~}{r}+[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha _0)]/\alpha _0\right)`$ $`=`$ $`\alpha (\stackrel{~}{r})\alpha \left(\stackrel{~}{r}+[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha _0)]/\alpha _0\right),`$ (32)
which allow us to rewrite eq. (5.2) in terms of $`\alpha (\stackrel{~}{r})`$ and its first derivative $`\alpha ^{}(\stackrel{~}{r})`$ \[with no other explicit functions of $`\stackrel{~}{r}`$ such as $`\delta (\stackrel{~}{r})`$ or $`\delta ^{}(\stackrel{~}{r})`$\]:

$`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ $`=`$ $`{\displaystyle \frac{\alpha ^{}(\stackrel{~}{r})}{1+\lambda ^3\varrho ^{3/\alpha }\psi (\alpha )^{13/\alpha }}}`$
$`\times \left[{\displaystyle \frac{\psi ^{}(\alpha )}{\psi (\alpha )}}{\displaystyle \frac{3}{\alpha }}\left({\displaystyle \frac{\psi ^{}(\alpha )}{\psi (\alpha )}}{\displaystyle \frac{\alpha \left(\stackrel{~}{r}+[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha )]/\alpha \right)\alpha }{\alpha ^{}(\stackrel{~}{r})}}\right)\right]+𝒪(ϵ^2).`$
Here, and in what follows, $`\alpha `$ stands for $`\alpha (\stackrel{~}{r})`$. This notation hides the explicit dependence of $`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ on $`ϵ`$, but for the calculation to be self consistent $`ϵ`$ has to be on the order of

$$ϵ\frac{\alpha \left(\stackrel{~}{r}+[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha )]/\alpha \right)\alpha }{\alpha }.$$
(34)
We thus require the right-hand side of eq. (34) to be much smaller than unity.

In principle, the sign of $`\mathrm{\Delta }\alpha `$ should allow an analysis of the evolution of $`\alpha `$ into a fixed point. However, the third term in eq. (5.2) is still problematic; it involves the function $`\alpha ()`$ evaluated recursively at a point which depends on $`\alpha `$ itself. In order to make this term useful, one must make further assumptions regarding the functional form of $`\alpha ()`$.

### 5.3 The Linear Regime

The most natural assumption for the functional form of $`\alpha (\stackrel{~}{r})`$ is that it is linear in $`\stackrel{~}{r}`$. As seen in Fig. 3 of DDH, this is clearly a good approximation for the profile used in our merger N-body simulation, and it is also valid for the family of profiles described in eq. (1). In other words, we assume that the second and higher derivatives of $`\alpha (\stackrel{~}{r})`$ are much smaller than the first derivative — at least in the interval $`[\stackrel{~}{r},\stackrel{~}{\mathrm{}_0}\stackrel{~}{\lambda }]`$. Under this assumption we may approximate the above problematic term by

$$\frac{\alpha \left(\stackrel{~}{r}+[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha )]/\alpha \right)\alpha (\stackrel{~}{r})}{\alpha ^{}(\stackrel{~}{r})}\frac{1}{\alpha }[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha )],$$
(35)
which brings eq. (5.2) to finally become

$`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ $``$ $`{\displaystyle \frac{\alpha ^{}(\stackrel{~}{r})}{1+\lambda ^3\varrho ^{3/\alpha }\psi (\alpha )^{13/\alpha }}}`$ (36)
$`\times {\displaystyle \frac{1}{\alpha ^2}}\left(\alpha (\alpha 3){\displaystyle \frac{\psi ^{}(\alpha )}{\psi (\alpha )}}+3[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha )]\right).`$
In this approximation, the sign of $`\mathrm{\Delta }\alpha (\stackrel{~}{r})`$ is the sign of the expression in big brackets, and the asymptotic fixed point slope $`\alpha _{\mathrm{as}}`$ is obtained by solving an explicit equation of $`\alpha `$:

$$\alpha (\alpha 3)\frac{\psi ^{}(\alpha )}{\psi (\alpha )}+3[\stackrel{~}{\varrho }+\stackrel{~}{\psi }(\alpha )]=0.$$
(37)
The factor in front of the big brackets of eq. (36) affects the rate at which $`\alpha (\stackrel{~}{r})`$ converges to the asymptotic slope. This concludes the derivation of eq. (4.3) and eq. (14) introduced in §4.3 for the case of a given mass ratio $`m/M`$.

## 6 A COSMOLOGICAL DISTRIBUTION OF SATELLITES

The analysis so far referred to mergers with a fixed mass ratio $`m`$. A real sequence of mergers consists of a cosmological mix of satellite masses. The probability for a mass ratio $`m`$ in a merger can be estimated using the Extended Press Schechter scheme (Lacey & Cole 1993, EPS), which has been tested and calibrated using N-body simulations (Lacey & Cole 1994; Somerville et al. 1999, Fig. 2). We have performed Monte Carlo realizations of merger trees based on EPS (using a scheme developed by Yuval Birnboim, to be published elsewhere, based on the scheme of Somerville & Kolatt 1999), and found that the probability distribution of $`m`$ is quite insensitive to the actual power spectrum of initial fluctuations, and it has no obvious time dependence. The fraction of number of mergers with mass ratio $`m<1`$ in the interval $`(m,m+\mathrm{d}m)`$ can be approximated by

$$p(m)\mathrm{d}mm^{3/2}\mathrm{d}m.$$
(38)
The relevant distribution of mass ratio for the buildup of the inner profile should also take into account the time it takes each satellite’s orbit to decay by dynamical friction from the halo virial radius $`R`$ to the central region. The duration of the decay process can be estimated for circular orbits in an isothermal halo (as in BT, eq. 7-26). We obtain

$$t_{\mathrm{fric}}=\frac{1.17}{m\mathrm{ln}(1/m)}\frac{R}{V_{\mathrm{circ}}}\frac{0.18}{m\mathrm{ln}(1/m)}t_\mathrm{i},$$
(39)
where $`t_\mathrm{i}`$ is the Hubble time when the merger starts with the satellite at the halo virial radius, and where an Einstein-deSitter cosmology is assumed at that time. On one hand eq. (39) is an underestimate because we have ignored the weakening effect of mass loss on the dynamical friction. If we adopt the mass-loss recipe $`m(r)M(r)`$ (which can be deduced from the resonance condition for isothermal halo and satellite) we obtain that $`t_{\mathrm{fric}}`$ becomes about twice as long as when mass loss is ignored. On the other hand eq. (39) is an overestimate in the case of eccentric orbits, where dynamical friction is strong already at early stages of the merger. In our merger simulations we find that a satellite on a typical orbit of eccentricity 1:6 decays to the center on a time scale shorter by a factor of $`3`$ compared to the decay time of a similar satellite on a circular orbit. Thus, the two effects roughly balance each other, and we can adopt eq. (39) as a crude approximation. Note that for $`m\stackrel{<}{}1`$ the estimate in eq. (39) breaks down (e.g., the log term drives $`t_{\mathrm{fric}}`$ to grow with $`m`$ for $`m>e^1`$). Since we wish to push the crude analysis to include such mergers despite the breakdown of our approximations there, we simply keep $`t_{\mathrm{fric}}`$ constant in this range.

The satellite is irrelevant for the cusp at $`t_0`$ if $`t_{\mathrm{fric}}(m)>t_0t_\mathrm{i}`$, which therefore defines a minimum relative satellite mass $`m_{\mathrm{min}}`$. For larger $`m`$ we crudely multiply the mass function of eq. (38) by the approximate correction factor $`f_{\mathrm{fric}}(m)=[1t_{\mathrm{fric}}(m)/(t_0t_\mathrm{i})]`$. For an order-of-magnitude estimate, we assume that a typical merger occurs at $`t_\mathrm{i}=t_0/2`$ (e.g., for an isothermal halo this is the typical time of collapse of the inner half mass of the halo). This implies $`m_{\mathrm{min}}0.067`$ and $`f_{\mathrm{fric}}(10.18/[m\mathrm{ln}(1/m)])`$ for $`m_{\mathrm{min}}me^1`$. A lower value for $`t_{\mathrm{fric}}`$, or an earlier $`t_\mathrm{i}/t_0`$, would correspond to lowering the factor 0.18 in $`f_{\mathrm{fric}}`$ and therefore to a smaller $`m_{\mathrm{min}}`$. Since $`m_{\mathrm{min}}0.067`$ is already much smaller than unity, the effect of the actual value of $`m_{\mathrm{min}}`$ on our final result is weak.

Given the corrected probability distribution $`f_{\mathrm{fric}}(m)p(m)`$, we compute the average over $`m`$ of $`\mathrm{\Delta }\alpha (r,m)`$ as given in eq. (4.3). We obtain the asymptotic slope $`\alpha _{\mathrm{as}}`$ for the distribution of satellites by requiring that the average increment vanishes:

$$_{m_{\mathrm{min}}}^1dmf_{\mathrm{fric}}(m)p(m)\mathrm{\Delta }\alpha (r,m)=0,$$
(40)
namely,

$`0`$ $`=`$ $`{\displaystyle _{m_{\mathrm{min}}}^1}dmf_{\mathrm{fric}}(m)p(m){\displaystyle \frac{1}{1+m^{1\nu (13/\alpha )}\psi (\alpha )^{13/\alpha }}}`$ (41)
$`\times \left(\alpha (\alpha 3){\displaystyle \frac{\psi ^{}(\alpha )}{\psi (\alpha )}}+3\mathrm{ln}[m^\nu \psi (\alpha )]\right).`$
Note that we have expressed $`\varrho `$ and $`\lambda `$ in terms of $`m`$ and $`\nu `$, using the scaling eq. (10). The factor $`\alpha ^{}(\stackrel{~}{r})/\alpha ^2`$ has been taken outside the integral and dropped from the equation because it is independent of the stochastic variable $`m`$ of the current merger (though it is affected by the mass ratios in previous mergers throughout the halo history).

Fig. 7 shows the value of the asymptotic profile obtained by solving eq. (41) for the above distribution of mass ratios and with $`\psi (\alpha )=0.5/\alpha `$ as a function of the power index $`n`$ of the fluctuation power spectrum \[$`\nu =(3+n)/2`$\]. Two values are also marked as horizontal bars in Fig. 6. For $`\mathrm{\Lambda }`$CDM, with $`n2.33`$, we obtain $`\alpha _{\mathrm{as}}=1.29`$. As seen in Fig. 6, the value of $`\alpha _{\mathrm{as}}`$ for a distribution of mass ratios is similar to the values obtained from eq. (14) for $`m0.32`$ \[which turns out to be close the average $`m0.33`$ of the distribution $`f_{\mathrm{fric}}(m)p(m)`$\], indicating that a typical merger in terms of its effect on the cusp profile is with $`m1/3`$. This is true for any $`n`$.

The reason for the robustness of the dominance by $`mm1/3`$ mergers and the insensitivity to the actual value chosen for $`m_{\mathrm{min}}`$ is the relative flatness as a function of $`m`$ of the term multiplying the expression in big parenthesis in the integrand, eq. (41). The factor $`f_{\mathrm{fric}}(m)`$ is roughly constant in most of the range, while $`p(m)m^{3/2}`$. The term expressed as a big fraction is roughly proportional to $`m^{1\nu (13/\alpha )}`$ which is about $`m^{0.5}`$ to $`m`$ in the relevant range of $`\nu `$ and $`\alpha `$. Together they introduce a moderately decreasing dependence on $`m`$, which does not become too large even at small $`m`$, as long as $`m_{\mathrm{min}}`$ is not significantly smaller than 0.01, say.

In Fig. 7 we see a gradual variation of $`\alpha _{\mathrm{as}}`$ with $`n`$: for the extreme case of $`n=3`$ we get a minimum slope of $`\alpha _{\mathrm{as}}=0.98`$ and for $`n=1.67`$ we obtain $`\alpha _{\mathrm{as}}=1.68`$. This range is comparable to the variation seen in simulations of a single cosmology, e.g., $`\mathrm{\Lambda }`$CDM. We shall see below (§7) that the measured cusp slopes are actually expected to be somewhat less $`n`$-sensitive than $`\alpha _{\mathrm{as}}`$. A systematic trend of such magnitude is hard to measure in cosmological simulations of different power spectra. There are hints for a similar trend already in the simulations of Cole & Lacey (1996, Fig. 9), and SCO report from their simulations an average of $`\alpha _{\mathrm{in}}=1.3\pm 0.07`$ for $`n=2`$, then $`\alpha _{\mathrm{in}}=1.6\pm 0.09`$ for $`n=1`$ and $`\alpha _{\mathrm{in}}=1.8\pm 0.09`$ for $`n=0`$.

As a side complementary argument, we use a very simplistic toy model to derive a crude upper limit for $`\alpha _{\mathrm{in}}`$. Following SCO, we consider the accumulation of undigested satellites in the halo center. Based on the assertion that the satellites’ inner densities are significantly higher than that of the host halo, one assumes in this model that the satellites accumulating in the halo center keep their original profiles unchanged by tidal effects. The mean halo profile that is built by a cosmological distribution of satellites is then

$$\overline{\rho }(r)=_{m_{\mathrm{min}}}^1dmf_{\mathrm{fric}}(m)p(m)\overline{\sigma }(m;r)$$
(42)
where $`\overline{\sigma }(m;r)`$ is the profile $`\overline{\sigma }(\mathrm{})`$ at $`r=\mathrm{}`$ of a satellite of mass ratio $`m`$, and where the $`m`$ dependence is given by the cosmological scaling for the given power spectrum. This integral is easy to evaluate if we assume that $`f_{\mathrm{fric}}(m)=\mathrm{const}.`$, adopt $`p(m)m^{3/2}`$ as in eq. (38) based on the realizations of the EPS model, and crudely take $`\overline{\sigma }(m;r)`$ to be a step function of height $`\sigma _\mathrm{s}`$ out to a radius $`\mathrm{}_\mathrm{s}`$, both scaled with $`m`$ as in eq. (10). At a given $`r`$, there is a contribution to $`\overline{\rho }(r)`$ from all satellites of $`\mathrm{}_\mathrm{s}>r`$, which translates by eq. (10) to a lower bound for the integration at a minimum mass $`m=r^{3/(1+\nu )}`$. The result, with $`\nu =(3+n)/2`$, is

$$\alpha _{\mathrm{in}}=3(4+n)/(5+n).$$
(43)
SCO argued for a flatter slope, $`\alpha _{\mathrm{in}}=3(3+n)/(5+n)`$, which could have resulted from eq. (42) if one had assumed $`f_{\mathrm{fric}}(m)p(m)m^1`$. Based on the above discussion of dynamical friction we expect $`f_{\mathrm{fric}}(m)`$ to be a moderately increasing function of $`m`$ for $`m>m_{\mathrm{min}}`$. This would flatten the slope predicted in eq. (43) and thus make it an upper bound for the actual cusp slope. An effective behavior of $`f_{\mathrm{fric}}(m)m^{1/3}`$, with $`p(m)m^{3/2}`$, would roughly recover the $`n`$ dependence found in the simulations of SCO. This simplified toy model ignores the additional effects of tidal mass transfer \[which we model in the current paper via $`\psi (\alpha )`$\] and the associated effects of tidal distortion of the satellites’ inner regions, which may go either way. Nevertheless, we notice that eq. (43) indeed overestimates our model prediction for $`\alpha _{\mathrm{as}}`$ for any $`n`$, and thus serves as an upper bound to $`\alpha _{\mathrm{in}}`$.

## 7 TOY SIMULATIONS

In order to test the predictions of the linear analysis of §5 for the asymptotic profile, and to study the evolution towards this asymptotic profile, we perform toy simulations of the profile buildup by a sequence of mergers, where we implement straightforwardly the mass-transfer recipe of eq. (4) with $`\psi (\alpha )=0.5/\alpha `$ in the range $`\alpha >1`$.333Certain modifications to the recipe may be needed near and below $`\alpha 1`$, where our adopted $`\psi (\alpha )`$ is not well established, and where deposit may occur without stripping. For example, if eq. (4) yields a deposit radius $`r<\mathrm{}`$, one may wish to replace it by $`r=\mathrm{}`$, assuming that the surviving satellite supports itself in a finite configuration while its center sinks to the halo center. We start with a halo of a given mean density profile $`\overline{\rho }(r)`$ and the corresponding monotonically increasing $`\alpha (r)`$. We consider a merging satellite of mass ratio $`m/M`$ with a homologous profile properly scaled using eq. (10). To materialize the merger we solve eq. (4) numerically for $`\mathrm{}(r)`$, with the possible modification for deposit without stripping at $`\alpha <1`$. We then assume that the halo profile changes only due to the added stripped satellite material, eq. (11), and compute the exact change in slope $`\mathrm{\Delta }\alpha (r)`$ at any $`r`$ using eq. (12). The following merger is performed with a satellite that is a scaled version of the new halo profile, and so on.

As a first simple test, we enforce at all times a mean density profile obeying the functional form of eq. (1), with $`r_\mathrm{s}`$ and $`\alpha _{\mathrm{out}}`$ fixed at 1 and 3 respectively. The inner slope $`\alpha _{\mathrm{in}}`$ starts at an arbitrary value in the range $`1\alpha _{\mathrm{in}}2`$ and is allowed to change as a result of the merger. We assume here that $`\alpha _{\mathrm{in}}`$ changes by $`\mathrm{\Delta }\alpha (r)`$ as evaluated at some small radius $`r`$. Fig. 8 shows the resultant evolution of $`\alpha _{\mathrm{in}}`$ for $`m/M=0.33`$ and $`r=0.01r_\mathrm{s}`$, when starting alternatively from either $`\alpha _{\mathrm{in}}=1`$ or $`\alpha _{\mathrm{in}}=2`$, and for three different power spectra of fluctuations. We see that, indeed, there is a slow convergence to an asymptotic profile. The obtained asymptotic slopes are close to the predictions of the linear analysis, Fig. 6. The slight deviations of less than 10% are mostly due to the fact that we compute $`\mathrm{\Delta }\alpha `$ at a finite small $`r`$ rather than at $`r0`$.

We next perform a non-parametric toy simulation where we do not enforce a specific functional form and do not constrain the outer slope to remain steep. We only make sure that $`\alpha (r)`$ is monotonic and relatively smooth, as required for the implementation of the deposit scheme. At each step, the mean density profile of the halo, $`\overline{\rho }(r)`$, is stored in an array of shells spaced logarithmically (1000 shells, $`\mathrm{\Delta }\mathrm{ln}r=0.05`$), and spline interpolated into any desired $`r`$. The slope $`\alpha (r)`$ is computed via 4’th order spline interpolation. To make sure that $`\alpha (r)`$ remains monotonic and smooth, $`\alpha (r)`$ is first smoothed with a Gaussian of width $`\mathrm{\Delta }\mathrm{ln}r=0.5`$, and then interpolated inside intervals of $`\mathrm{\Delta }\mathrm{ln}r=2`$ (with an overlap of $`\pm 1`$ with the neighboring intervals) using locally a functional form for $`\alpha (r)`$ corresponding to eq. (1) with all the parameters free. This is a general smoothing procedure that is not sensitive to the use of this specific functional form in the local fit. Our initial profile is again given by eq. (1), with $`r_\mathrm{s}`$ put at the center of the sampled $`\mathrm{ln}r`$ range.

The ultimate test of the linear approximation is shown in Fig. 6, where the results of the non-parametric toy simulations are presented by the symbols on top of the curves describing the linear predictions. We see that the actual sequence of mergers with a fixed mass ratio $`m/M`$ converges to an asymptotic slope very similar to that predicted by the linear approximation, for all values of $`m/M`$ in the range of interest. Significant deviations are obtained only for very small $`m/M`$ values, below any realistic value of $`m_{\mathrm{min}}`$ as implied by dynamical friction time-scales.

Fig. 9 shows how $`\alpha (r)`$ actually evolves in the non-parametric toy simulations where no constraints are applied on the outer slope. When $`\alpha _{\mathrm{in}}`$ is below unity there is a rapid evolution towards $`\alpha 1`$ in the range $`r<r_\mathrm{s}`$. This is partly due to the factor $`1/\alpha ^2`$ driving $`\mathrm{\Delta }\alpha `$ in eq. (12) and partly due to the processes explained in DDH as reflected in the shape adopted for $`\psi (\alpha )`$ at $`\alpha <1`$. A power-law region with $`\alpha 1`$ develops below the radius where $`\mathrm{\Delta }\alpha =0`$, and it quickly extends to smaller radii. The value of $`\alpha `$ in this inner power-law region grows slowly from near unity towards the asymptotic value $`\alpha _{\mathrm{as}}`$. In parallel, since no constraints are applied on the outer slope, it can decrease slowly from $`\alpha _{\mathrm{out}}=3`$ towards $`\alpha _{\mathrm{as}}`$. The shape of the profile during the long interim period resembles a cusp-halo structure qualitatively similar to the NFW-like profiles of halos seen in cosmological simulations. It takes 10 to 20 mergers (of $`m/M=0.33`$), during which the halo mass grows by about two orders of magnitude, for the shape of the profile to strongly deviate from this characteristic shape. In reality, we may expect the outer profile to keep a steep slope in response to other physical processes not modeled here (see §1), and therefore the asymptotic profile to resemble a shape qualitatively similar to eq. (1) with $`\alpha _{\mathrm{in}}=\alpha _{\mathrm{as}}`$.

## 8 DISCUSSION AND CONCLUSIONS

We have addressed the origin of a cusp in dark-matter halos within the hierarchical clustering scenario, based on tidal effects during the halo buildup by merging of relatively compact satellite halos. In Dekel, Devor & Hetzroni (2003) we highlight the rapid steepening of a flatter core of slope $`\alpha 1`$ into a cusp of $`\alpha >1`$ as a result of vanishing tidal mass transfer in such a core. Merger N-body simulations indicate here that the cusp is quite stable under a sequence of such mergers. we have developed a simple prescription for tidal mass transfer in mergers, argued using N-body simulations that it provides a useful approximation, and showed that it naturally leads to a stable cusp with an asymptotic profile of slope $`\alpha _{\mathrm{as}}`$ slightly larger than unity. Our toy model is one way to explain the cusps seen in cosmological N-body simulations of the hierarchical clustering scenario, and it provides a tool for addressing other processes which may explain the observed flat core in some galaxies.

The mass-transfer recipe relates each radius $`\mathrm{}`$ in the initial satellite with a deposit radius $`r`$ in the halo via the ratio of mean densities $`\overline{\rho }(r)/\overline{\sigma }(\mathrm{})=\psi [\alpha (r)]`$. We argued based on qualitative considerations and found in merger N-body simulations that $`\psi (\alpha )`$ has a characteristic shape which is not too sensitive to the details of the merger orbit; $`\psi (\alpha )`$ is a decreasing function of $`\alpha `$ and it obtains values significantly smaller than unity. A crude fit seems to be provided by $`\psi (\alpha )=0.5/\alpha `$ in the range $`1\alpha 3`$. This makes the mass transfer more efficient than what might be inferred from the naive resonance condition, $`\psi =1`$, and gradually more efficient where the halo profile is steeper. Beyond our immediate purpose of studying the formation of a cusp, this mass-transfer recipe provides a useful tool for general studies using semi-analytic models of galaxy formation.

Admittedly, we have only tested in this paper the robustness of $`\psi (\alpha )`$ to the merger orbit type. The choice of typical properties for the simulated halo and satellite, which relate to each other in general agreement with the scaling relation expected in the $`\mathrm{\Lambda }`$CDM cosmology, makes the obtained $`\psi (\alpha )`$ a sensible crude guess. Nevertheless, the robustness of $`\psi (\alpha )`$ to variations about the simulated properties of halo and satellite must be tested using a more comprehensive set of simulations. We have adopted $`\psi (\alpha )=0.5/\alpha `$ as an illustrative example for obtaining actual numerical values for $`\alpha _{\mathrm{as}}`$, but the convergence analysis, which is the heart of this paper, is much more general.

After understanding why a decreasing $`\psi (\alpha )`$ should lead to a stable asymptotic profile, we have performed a linear perturbation analysis of this process, and confirmed its validity by toy simulations of tandem mergers based on the adopted mass-transfer recipe. The asymptotic slope $`\alpha _{\mathrm{as}}`$ is typically somewhat larger than unity. The system develops a cusp of $`\alpha 1`$ which slowly grows towards $`\alpha _{\mathrm{as}}`$ while the outer halo maintains a steeper profile for a long time. For a typical cosmological distribution of merging satellites in the $`\mathrm{\Lambda }`$CDM cosmology ($`n=2.33`$) the asymptotic slope is $`\alpha _{\mathrm{as}}1.3`$. This is close to the typical cusp slope found in cosmological N-body simulations (e.g. Power et al. 2002 and other references in §1). The typical mergers driving this process are with a relatively high mass ratio of $`m/M1/3`$. This means that one should not take our quantitative results too literally because our toy model assumes $`mM`$.

The asymptotic profile obtained in our model has a relatively weak dependence on the power spectrum of fluctuations. In particular, the effect of tidal compression (§2) and the $`1/\alpha ^2`$ factor in $`\mathrm{\Delta }\alpha `$ \[eq. (12)\] drive the profile to $`\alpha \stackrel{>}{}1`$ independently of the fluctuation power spectrum. Also, the cosmological dependence of the ultimate asymptotic slope is weakened by the fact that the profile is determined by mergers of relatively large mass ratio. For values of power-spectrum index in the range $`3<n<1`$ the asymptotic slope is predicted to vary gradually in the range $`1<\alpha _{\mathrm{as}}<2`$. This can be compared to the only slightly weaker trend indicated in the cosmological N-body simulations of SCO.

Syer & White (1998) also addressed the profile resulting from mergers and obtained results somewhat different from ours. They adopted the condition of tidal stripping at resonance, $`\psi =1`$, ignoring the $`\alpha `$ dependence of $`\psi `$. When we substitute the recipe $`\psi =1`$ in eq. (12), we find that $`\mathrm{\Delta }\alpha `$ is always positive, for any merger and at any $`r`$. Indeed, when trying to repeat the Syer & White toy simulations with higher resolution and following more mergers we actually find that the profile does not really converge to a stable cusp but rather continue to steepen slowly towards $`\alpha _{\mathrm{in}}=3`$. Only when using the revised mass-transfer prescription where $`\psi `$ is decreasing with $`\alpha `$ do we obtain convergence to a flatter asymptotic profile.

Syer & White (1998), Nusser & Sheth (1999) and Subramanian, Cen & Ostriker (2000) describe different toy models which predict $`\alpha _{\mathrm{in}}=(9+3n)/(5+n)`$,444Nusser & Sheth also predict an upper bound of $`\alpha _{\mathrm{in}}=(9+3n)/(4+n)`$. recovering the hierarchical-clustering toy model by Peebles (1980, eq. 26.8). While this scaling relation may be in reasonable agreement with the simulations (and with our model predictions) for $`n1`$, it admits values below $`\alpha _{\mathrm{in}}=1`$ for $`n<2`$, in conflict with the findings of the simulations. Our revision of the toy model of SCO yields $`\alpha _{\mathrm{in}}<3(4+n)/(5+n)`$, to be flattened further by the $`m`$ dependence of dynamical friction. The basic model analyzed in the current paper does not admit cores significantly flatter than $`\alpha _{\mathrm{in}}1`$ (that is by non-dissipative processes only) and is thus in better agreement with the simulation results.

The gravitational processes leading to a cusp are modeled in our analysis as tidal effects during the buildup of the halo by a sequence of mergers. This picture is likely to be valid in the CDM hierarchical clustering scenario, where numerous sub-galactic halos exist on all scales and are continuously merging (e.g., Klypin et al. 1999b; Moore et al. 1999a; Springel et al. 2001). This merger picture is confirmed by a careful inspection of high resolution CDM simulations, where a special effort is made to identify the merging clumps which otherwise could have been easily missed (e.g., Wechsler, Dekel et al.  , in preparation). It is therefore clear that cores, independently of how they form, cannot survive in a pure CDM scenario without some modification — they efficiently turn into cusps as a result of the mergers with the relatively compact building blocks.

Nevertheless, a cusp, though somewhat flatter, is reported to be seen also in simulations where the initial fluctuations had less power on small scales, thus suppressing the number of sub-galactic satellites and the associated merger rate (Moore et al. 1999b; Avila-Reese et al. 2001; Bullock, Kravtsov & Colin 2002). Further indications for the generality of cusp formation comes from simulations by Huss, Jain & Steinmetz (1999), who find that cusps also form as a result of collapse from spherical initial density perturbations with suppressed random velocity perturbations, as well as from the simulations of Alvarez, Shapiro & Martel (2002), who find that cusps arise from the gravitational instability and fragmentation of cosmological pancakes. One way to explain this is by noticing that the asymptotic cusp formed in the CDM scenario is driven by mergers with relatively massive satellites (of typical mass ratio 1:3, see Dekel et al. 2002), and realizing that such mergers do happen even when small-scale power is suppressed and when pancakes fragment. It is therefore possible that the cusp is actually driven by mergers to a certain extent also in some of these cases. Our toy model based on mergers would not be directly applicable in cases where a cusp arises while the halo forms by smooth accretion. The gravitational processes involved in such cases somehow mimic a behavior similar to the merger case. We note, for example, that the tidal compression discussed in DDH is expected to amplify density perturbations and possibly make them behave in certain ways like merging clumps.

However, one should accept in general that the halo buildup is a complex gravitational process, whose different aspects can probably be modeled in more than one way, e.g., as a violent relaxation process driven by fluctuations in comparison with a sequence of mergers and substructure accretion.555We encounter an analogous duality, for example, when one manages to explain the origin of galactic spin alternatively via tidal torque theory applied to shells and by summing up the orbital angular momenta in a cosmological sequence of mergers (see Maller, Dekel & Somerville 2002). There is also the intriguing possibility that general statistical considerations in phase space may provide some clue for the origin of a universal halo profile (e.g., Taylor & Navarro 2001). Unfortunately, we currently know of no viable alternative to the tidal effects in mergers as a simple model for the origin of $`\alpha \stackrel{>}{}1`$ cusps in the cosmological simulations. The idealized merger picture provides one possible toy model, approximately valid at least in hierarchical clustering scenarios, within which we understand the origin of the cusp in simple terms. However, the complexity of the problem suggests that this is not the final word on the issue.

Our result implies a necessary condition for the survival of cores in halos independent of their origin — that satellites should be prevented from adding mass to the halo cores. This could be avoided in a CDM scenario if feedback processes manage to puff up the satellites and make them disrupt before they merge with the halo cores (as proposed in DDH). We have not explicitly addressed in this paper the sufficient conditions for the formation of cores, but one can imagine that a sequence of mergers with low-density satellites, where the mass is predominantly deposited outside the inner halo region, would indeed flatten the inner profiles (work in progress). Other processes may also contribute to the development of halo cores. The proposed scenarios include, for example, the disruption of cusps in merging satellites by massive black holes (Merritt & Cruz 2001), the heating by gas clouds spiraling in due to dynamical friction (El-Zant, Shlosman & Hoffman 2002), the angular-momentum transfer from a big temporary rotating bar (Weinberg & Katz 2002), and the delicate resonant reaction of halo-core orbits to the tidal perturbation by the satellite, which could be a strong effect if the dark-matter distribution is much smoother than the current state-of-the-art N-body simulations (Katz & Weinberg 2002). We stress again that, no matter what the origin of the core might be, our analysis implies that such cores could survive only if they are not perturbed by significant mass transfer from merging compact satellites, which implies that small CDM halos must be puffed up before they merge into bigger halos. We also stress that none of these scenarios seem to be capable of explaining the formation of cores in halos on the scales of clusters of galaxies within the CDM scenario.

Supernova feedback effects are probably not strong enough for turning cusps into cores in halos with rotation velocities higher than $`100\mathrm{km}\mathrm{s}^1`$ (Geyer & Burkert 2001; Gnedin & Zhao 2002). However, they are possibly sufficient for the necessary indirect puffing-up of the merging satellites. Gnedin & Zhao (2002) estimate that direct feedback effects may reduce the central satellite densities by a factor of 2 to 6. Based on our simulations with puffed-up satellites, this may be enough by itself to avoid the steepening from a core to a cusp in the framework of the merging scenario. However, if the (so far inconclusive) observational clues for cores in clusters of galaxies are confirmed, supernova feedback is unlikely to provide a viable explanation for their origin. In this case one may search for higher efficiency in supernova feedback either due to microscopic effects such as porosity in a muliphase ISM or due to hypernova from very massive stars (Silk 2002). Alternatively, one may appeal to stronger feedback mechanisms, perhaps associated with radio jets from AGNs. This process may be indicated by an observed correlation between AGN activity and bright galaxies in SDSS (Kauffmann et al. 2003, in preparation), and it may be needed independently in order to explain the missing baryons in big galaxies (Klypin, Zhao & Somerville 2002) and in clusters. Otherwise, such large cores may present a real challenge for the standard hierarchical clustering scenario

The cusp/core problem is only one of the difficulties facing galaxy formation theory within the CDM cosmology. It turns out that other main problems can also be modeled by tidal effects in mergers, and may also be resolved by the inevitable feedback processes. For example, Maller & Dekel (2002) addressed the angular-momentum problem, where simulations including gas produce disks smaller than the galactic disks observed (Navarro & Steinmetz 2000 and references therein; Governato et al. 2002), and with a different internal distribution of angular momentum (Bullock et al. 2001b; van den Bosch, Burkert & Swaters 2001). A toy model has been constructed for the angular-momentum buildup by mergers based on tidal stripping and dynamical friction, which helps us understand the origin of the spin problem as a result of over-cooling in satellites. A simple model of feedback has then been incorporated, motivated by Dekel & Silk (1986). This model can remedy the discrepancies, and in particular explain simultaneously the low baryon fraction and angular-momentum profiles in dwarf disk galaxies.

Various feedback effects may also provide the cure to the missing dwarf problem, where the predicted large number of dwarf halos in CDM can possibly match the observed number of dwarf galaxies only if the mass-to-light ratio in these objects is very high (Klypin et al. 1999b; Moore et al. 1999a; Springel et al. 2001; Kochanek 2001). Bullock, Kravtsov & Weinberg (2000), Somerville (2002) and Tully et al. (2002) appeal to radiative feedback effects which prevent the formation of small dwarfs after cosmological reionization at $`z7`$, Scannapieco, Ferrara & Broadhurst (2000) and Scannapieco & Broadhurst (2001) address the destructive effect of outflows from one galaxy on neighboring protogalaxies via ram pressure, and Dekel & Woo (2002) study the role of supernova feedback in determining the relevant global properties of dwarfs and larger low-surface-brightness galaxies. We note that while the requirements from feedback in explaining the dwarf-galaxy properties and the angular-momentum problem are not too demanding, the solution to the core problem requires that the dark-matter distribution be affected by feedback, which is a non-trivial requirement.

Nevertheless, the successes of such toy models in matching several independent observations indicate that they indeed capture the relevant basic elements of the complex processes involved, and in particular that feedback effects may indeed provide the cure to some or all the main problems of galaxy formation theory within the $`\mathrm{\Lambda }`$CDM cosmology that does so well on larger scales. The alternative solution involving Warm Dark Matter (e.g., Hogan & Dalcanton 2000; Avila-Reese et al. 2001; Bode, Ostriker & Turok 2001) seems to still suffer from the cusp problem, it may still fail to reproduce the angular-momentum profile in galaxies (Bullock, Kravtsov & Colin 2002), and it may be an overkill where the formation of dwarf galaxies is totally suppressed once the inevitable feedback effects are included (Bullock 2001). The speculative alternatives involving self-interacting dark matter are even more problematic (Spergel & Steinhardt 2000; Dave et al. 2001; Hennawi & Ostriker 2002).

## Acknowledgments

We acknowledge stimulating discussions with Oded Ben David, George Blumenthal, Andi Burkert, Doug Lin, Ari Maller and Gary Mamon. This research has been supported by the US-Israel Bi-National Science Foundation grant 98-00217, the German-Israel Science Foundation grant I-629-62.14/1999, and NASA ATP grant NAG5-8218.