format: text

text:
# The Dynamical Evolution of Substructure

## 1 Introduction

In the context of Cold Dark Matter dominated (CDM), hierarchical clustering cosmology, large galaxies are formed from the merging and accretion of small, less massive progenitors (clumps). The exact role played by this substructure in the formation of galaxies is still unclear, and there remain unresolved several very important and fundamental issues that are closely related to the role of these clumps.

The angular momentum content of galaxies is a fundamental parameter. Cosmological N-body simulations have verified the proposal that tidal torques between proto-galaxies generate angular momentum (Peebles 1969), and have shown that the dimensionless spin parameter $`\lambda J|E|^{1/2}G^1M^{5/2}`$ has a well-defined mean value of $`0.06`$, independent of the details of the power spectrum of density fluctuations (Efstathiou & Jones 1979; Aarseth & Fall 1980; Barnes & Efstathiou 1987; Warren et al. 1992; Steinmetz & Bartelmann 1995; Cole & Lacey 1996). Since the spin parameter $`\lambda `$ is much smaller for present day ellipticals – indeed equal to the value of $`0.06`$ of a typical dark halo – than is the value of $`\lambda `$ for spirals, while the effective radius for spirals is comparable to that of ellipticals at a given luminosity (Fall 1980), it is important to know how the ellipticals succeed in losing angular momentum, compared to the expectation of smooth collapse and spin-up within a dark halo (Fall & Efstathiou 1980) while still dissipating the same amount of binding energy as do proto-spirals. On the other hand, there exists an angular momentum problem in that while disks are formed in fully self-consistent hierarchical-clustering simulations, the angular momentum transport inherent in the merging process produces disks that are too centrally-concentrated and contain too little angular momentum compared with observed spirals (e.g. Navarro & Benz 1991; Evrard, Summers & Davis 1994; Vedel, Hellsten & Sommer-Larsen 1994; Navarro & Steinmetz 1997). Indeed standard analytic and semi-analytic models of the formation of spiral galaxies show that extended disks, as observed, form through detailed angular momentum conservation, with proto-disk material retaining the same value of specific angular momentum as the dark halo (e.g. White & Rees 1978; Fall & Efstathiou 1980; Gunn 1982; Jones & Wyse 1983; Dalcanton, Spergel & Summers 1997; Mo, Mao & White 1998; Zhang & Wyse 2000, Silk 2001). Hence models must retain angular momentum for spirals, while losing it for ellipticals. Mergers and the associated gravitational torques and angular momentum transport clearly play a role. However, detailed merger simulations of two disk galaxies have shown that there are difficulties with the simplest merger picture of how ellipticals form (Bendo & Barnes 2000; Cretton, Naab, Rix & Burkert 2001), in that equal mass mergers seem to produce too much kinematic misalignment, while in unequal mass mergers the strong signatures of the more massive disk survive.

What are the effects of substructure such as gas clouds, stellar (globular?) clusters, dwarf galaxies or sub-halos which fall later into the galactic-scale potential well? Do they significantly change the angular momentum distribution and further change the final galaxy classification? An analytical approach will have difficulties in dealing with this problem, although the secondary infall model is successful in describing the smooth post-collapse halo profile (Gunn & Gott 1972; Gunn 1977; Filmore & Goldreich 1984; Hoffman & Shaham 1985; Bertschinger & Watts 1988; Zaroubi, Naim & Hoffman 1996). However, numerical simulations have the capability of treating the formation and the non-linear evolution of dark halos in a cosmological context (Quinn, Salmon & Zurek 1986; Zurek, Quinn & Salmon 1988; Quinn & Zurek 1988; Warren et al. 1992; Dubinski & Carlberg 1991). These studies show that the angular momentum and binding energy are indeed redistributed, in an orderly manner, during the relaxation of the halos (Quinn & Zurek 1988). Further, during the merging process inherent in hierarchical clustering, baryonic substructure can lose angular momentum to the dark halo (Frenk et al. 1985; Zurek, Quinn & Salmon 1988; Barnes 1988; Navarro & Steinmetz 1997; Navarro, Frenk & White 1995;). Zurek, Quinn & Salmon (1988) suggest that whether the baryonic clumps are gas clouds (and hence dissipative) or stellar could be an important criterion in determining whether the forming galaxy will become a spiral or an elliptical.

Proposed solutions to the problem with disks, within the context of CDM cosmologies, have included delaying the onset of disk formation until after the merging process is essentially complete, by appeal to suitable ‘feedback’ from the first stars (Weil, Eke & Efstathiou 1998), but this solution must address the old age for disk stars in the solar neighborhood (e.g. Edvardsson et al. 1993; Wyse 1999) and in the outer parts of the M31 disk (Ferguson & Johnson 2001), in addition to the observation of apparently fully-formed disks at redshifts above unity (Brinchmann & Ellis 2000) with a Tully-Fisher relation that is essentially indistinguishable from that of present-day disks, apart from passive evolution. The density and mass of the substructure are clearly important, and a major motivation for the present work is to quantify the effects of substructure of a range of each. Further, when is substructure disrupted and why?

Another issue we seek to clarify is how much of the angular momentum change seen in the cosmological simulations is due to the boundary conditions, namely that the system under study is not closed. It is hard to separate the effects of the real angular momentum transfer from clumps to halo from the angular momentum gain by the accretion of material or the torques applied by the neighboring environment or the later infall of dark matter. Thus it is worthwhile to study, as here, the evolution of angular momentum in isolated systems. In the isolated N-body simulations by Barnes (1988) of the encounter of two galaxies, it is found that the dark halo can extract some angular momentum from other components.

In this paper, we shall extend such N-body studies of the dynamical evolution to analyse the evolution of a system of clumps in a larger, smooth, isolated system. Our simulations have relevance beyond a forming proto-galaxy. We may address the dynamical processes occurring in a virialized system with substructure, such as a cluster of galaxies or a single galaxy and its retinue of dwarf companion galaxies. The relative importance of the effects of the dissipationless collapse of the overall system, tidal disruption of clumps, dynamical friction of clumps, clump-clump heating and merging can possibly be related to the formation of cD galaxies in clusters of galaxies, to the formation of central density cusps in dark halos, to the triaxial nature of dark halos, to the formation of substructure composed of disrupted debris in the outer halo regions of galaxies, etc. Thus we initiated a consistent comparison among these dynamical processes in the context of galaxy formation.

In section 2, we present the parameters of our N-body simulations. Section 3 contains a discussion of our expectations for the dynamics of the substructure based on analytic calculations and comparisons with the N-body simulations, while section 4 discusses the angular momentum content and transport. The internal kinematics of the different components are discussed in Section 5, while their morphology is discussed in Section 6. Section 7 summarizes our conclusions. The more technical aspects are given in the Appendices.

## 2 The N-body Models

### 2.1 The Technique

We adopt a numerical N-body code based on the hierarchical tree algorithm (Barnes & Hut 1986; Hernquist 1987). Our simulations study the dynamical evolution of a number of clumps embedded in a smooth dark halo. Due to the collisionless property of the particles used in our N-body simulation, it should be emphasized that the clump can represent either stellar or dark matter depending on our interpretation of the physical entity. The initial dark halo profile is modelled with either the Hernquist profile (Hernquist 1990) or the Plummer profile (Plummer 1911), the details of which are given in Appendix A. The Hernquist profile, with $`\rho r^1`$ at the centre and $`\rho r^4`$ at large radius, was introduced as being, when projected, close to the observed surface brightness profiles of elliptical galaxies, and is close to the mass profile seen in hierarchical clustering simulations (the NFW profile; Navarro, Frenk & White 1997), which has a cuspy centre, $`\rho r^1`$ (though flatter than the $`r^{1.5}`$ behaviour of high-resolution simulations; Ghigna et al. 2000), and $`\rho r^3`$ at large radius. The Plummer profile has close to a harmonic potential in the core (only slowly varying density) and $`\rho r^5`$ at large radius. In all realisations, each clump is simulated with a Plummer profile.

We ran simulations of systems in initial equilibrium plus out-of-equilibrium systems that underwent an initial collapse. In order to study the angular momentum behaviour of the clumps and the dark halo, we require that the system contain enough angular momentum, consistent with the cosmological initial conditions expected, i.e. $`\lambda 0.1`$ (e.g. Efstathiou & Jones 1979; Barnes & Efstathiou 1987). First, we build the dark halo with no streaming velocity. Then, we add some angular momentum about the $`z`$ direction, which is defined to be the axis of the halo. The orbital rotation velocity added to each particle is the component projected on the equatorial plane of a constant fraction, $`\beta `$, of the circular velocity where the particle is located, i.e., $`V_{rot}=\beta V_c(r)cos\theta `$, where $`\theta `$ is the angle between $`\widehat{𝐫}`$ and $`\widehat{𝐳}`$. We further stretch the spherical profile to an oblate profile with ellipticity $`ϵ_0=0.53`$ (Warren et al. 1992). Assuming, if the system is virialized, that the flattening is only caused by the addition of rotation, we adopt the empirical relation $`\beta =\sqrt{ϵ_0/5}`$, which may be seen by noting the scaling of the rotation parameter $`v_0/\sigma \sqrt{ϵ_0}`$ for isotropic systems flattened by rotation (Binney & Tremaine 1987). With this choice of $`\beta `$, the initial spin parameter of the halo in the virialized simulations (with the virial ratio $`2T/W1`$, where $`T`$ is the total kinetic energy and $`W`$ is the total potential energy) is $`\lambda 0.16`$. The collapse simulation models (with the virial ratio $`2T/W0.1`$) have $`\lambda 0.07`$. In order to achieve an initial condition suitable for pre-collapse, i.e. with small $`2T/W`$, we decrease the total velocity of the dark halo particles and the velocity of each clump by a factor of $`\sqrt{10}`$ for $`2T/W0.1`$.

No internal rotation is added to the clumps, which are spherical. The distribution of clumps is chosen to follow the dark halo profile with similar treatment of the flattening and rotation velocity. For the simulation models with the same number of clumps, we require the orbital characteristics for the clumps be the same with the exception of models B and E (see Table 1).

### 2.2 The Parameters and Units

The parameter values of interest for our simulations are given in Table 1. The total number of particles in our simulations ranged from 20,000 to 120,000, with smooth dark halo particles and clump particles having an equal mass. The total number of clumps, $`n_c`$, varied from 5 to 80 and the total clump mass fraction, $`f`$, varied from 5% to 100%. These ranges were chosen to study the effect of clumpiness on the outcome of the simulations. The core radius of a clump (a Plummer sphere) is set to achieve the chosen density contrast, $`\rho _{cl}/\rho _0`$, the ratio of the mean density within the half mass radius of the clump to that of the whole system (dark matter plus the clumps). We vary the density parameter $`\rho _{cl}/\rho _{DM}`$ (with $`\rho _{cl}/\rho _0=(1f)\rho _{cl}/\rho _{DM}`$) over the range of 3 to 64 in our simulations, to study the effects of global tides in disrupting the clumps. A value of three for the density contrast is the minimum expected for a system to withstand global tides (the Roche criterion).

The total mass of the system is normalized to be unity. The unit of distance is the core size of the dark halo ($`a=1`$, see Appendix A) and we further set the gravitation constant $`G=1`$.

With these units and as detailed in Appendix A, the mean crossing time for a Hernquist halo is 5.3, while that for a Plummer halo is 2.1, and the simulations are run for $`25`$ crossing times, or approximately a Hubble time in physical units.

It is convenient to express some of our parameters, such as the clump mass fraction $`f`$, the number of clumps $`n_c`$, density contrast of each clump $`\rho _{cl}/\rho _0`$, and crossing time $`t_{\frac{1}{2}hc}`$, in terms of dimensionless quantities. From our definition of the density contrast of the clump, the ratio of the half mass radius of a clump to that of the whole system can be written as

$$r_{\frac{1}{2}c}/r_{\frac{1}{2}h}=f^{1/3}n_c^{1/3}\left(\frac{\rho _{cl}}{\rho _0}\right)^{1/3}.$$
(1)
The ratio of the one-dimensional internal velocity dispersion of the clump and the halo is $`\sigma _c/\sigma _h=\frac{v_{\frac{1}{2}c}}{v_{\frac{1}{2}h}}=\frac{r_{\frac{1}{2}c}}{r_{\frac{1}{2}h}}\frac{t_{\frac{1}{2}hc}}{t_{\frac{1}{2}cc}}=\left(\frac{r_{\frac{1}{2}c}}{r_{\frac{1}{2}h}}\right)\left(\frac{\rho _{cl}}{\rho _0}\right)^{1/2}`$. This can be further written as

$$\sigma _c/\sigma _h=f^{1/3}n_c^{1/3}\left(\frac{\rho _{cl}}{\rho _0}\right)^{1/6}.$$
(2)
### 2.3 Numerical Effects

As demonstrated in Appendix B, while there may well be spurious clump heating effects in our simulations, due to the limited number of particles used, they do not pose a problem for our analysis. The softening-length in the treecode is set at $`r_s=0.5r_{\frac{1}{2}cl}N_{cl}^{1/3}`$, where $`N_{cl}`$ is the number of particles used in each clump. This relatively large value is chosen to suppress the unphysical two-body relaxation (White 1978) which otherwise could occur due to the small number of particles used in each clump (ranging from 150 to 2400).

Figure 1 shows the mass profiles for examples of unperturbed clumps consisting of 150 particles (from simulation A) and 900 particles (from simulation H) at the beginning (solid line) and the end (dashed line) of our simulations. The radius of a given mass fraction does not change much even for the clump with 150 particles, demonstrating that indeed relaxation effects are not important.

### 2.4 The Analysis Method

In this section, we describe the methods we used in calculating the quantities of physical interest, such as the angular momentum, the properties of the debris from the disrupted clumps, and the ellipticity or triaxiality, density and kinematics of the clumps and the halo.

In calculating the angular momentum, one needs to choose the origin carefully. Here we calculate the angular momentum of the different components of the system in two ways, relative to the centre of mass of the whole system and relative to the highest density point of the dark halo. The centre of mass of the whole system can be calculated directly, while the highest density point of the dark halo can be searched for in an iterative way: first we choose the centre of mass of the whole system as the starting point $`C_0`$, and calculate the centre of mass $`C_1`$ of the dark halo particles contained within a radius $`r_1=20`$ relative to the point $`C_0`$. Then we adopt the point $`C_1`$ as the next starting point to repeat the above process until the centre of mass converges; then we decrease the radius $`r_1`$ by half and repeat the first step. After several iterations of these two steps, the highest density point is located when the search radius $`r_1`$ has decreased greatly. Throughout this paper we adopt the highest density point as the origin when calculating other physical quantities such as the density profile, ellipticity and kinematics.

The debris disrupted from each clump is calculated in an iterative way: at some timestep, e.g. initially, for each clump we know those particles still bound to that clump; at the next timestep we check which of those particles remain bound to the clump, and we assign the unbound particles to the debris. We repeat this process at each timestep until the end of the simulation is reached or the clump is completely disrupted. It should be noted that our simulation code is not designed to search for merging events, which would require a more sophisticated algorithm that includes a check on a merging criterion at each timestep of the integration. Thus the ‘debris’ we identify could include some clump particles that are no longer bound to their initial clump due to the fact that this initial clump has merged with another clump, and the particles are bound to the new merger remnant, rather than being truely unbound. In other words, we do not re-assign particles to a new, larger clump should one clump be subsumed in another. We simply keep track of whether the particles remain bound or not to the initial clump. However, our calculation of the debris would be inaccurate only in a very few cases in which the merging process is very efficient.

In addition, when identifying debris, we ignore the possibility that a clump can capture particles that were removed from another clump, but this is not an important process for our simulations which do not include a dissipative component.

The intrinsic triaxiality of the final dark halo or debris component is calculated in a simple way. We calculate the moment of inertia formed by the relevant particles within a chosen radius (e.g. half mass radius). From the eigenvalues of the moment of inertia $`I_1<I_2<I_3`$, we can fit the moment of inertia with a triaxial ellipsoid with the intrinsic axial ratios $`ϵ_1=b/a`$ and $`ϵ_2=c/a`$ where $`abc`$. It can be calculated that $`ϵ_1=\sqrt{\frac{PQP+Q}{PQPQ}}`$ and $`ϵ_2=\sqrt{\frac{PQ+PQ}{PQPQ}}`$ where $`P=I_1/I_2`$ and $`Q=I_1/I_3`$. The intrinsic axial ratios calculated in this way can be slightly overestimated compared to their real values, but provide a consistent comparison among our simulations. The ellipticity of the projected (surface density) images of the final dark halo or debris can be calculated in a more direct way, for a range of viewing angles. We calculate an isodensity contour map by assigning a Gaussian density profile to each particle, with width proportional to the softening length. Then we simply fit the contour with the ellipse to obtain the axial ratio $`b/a`$ and the ellipticity $`ϵ=1b/a`$.

The radial density profile is calculated by dividing the particles into radial bins, each containing equal numbers of particles (about 16). The velocity dispersion and rotation velocity are calculated by the conventional approach.

## 3 The Dynamics of Clumps

We can study the dynamical evolution of clumps that contain enough particles. The study of the dynamics of substructure has applications on both the galaxy size scale, where star clusters, gas clouds and dwarf satellite galaxies play the role of clumps (Fall & Rees 1977), and the galaxy cluster scale, where galaxies play the role of clumps (Dressler 1984). It is also related to the ‘undermerging’ problem that has emerged recently, in that high-resolution CDM cosmological simulations predict too many surviving dark matter satellites around halos of the size of our Galaxy (Klypin et al. 1999; Moore et al. 1999)).

The dynamical processes associated with the evolution of these clumps we shall investigate are global tidal stripping, dynamical friction, and merging and close encounters between clumps. We first derive expectations, based on analytic arguments, for the roles these dynamical processes have played and on which parameters they depend, for the virialized models. A preliminary comparison between these theoretical estimates and our simulations is provided at the end of the section. Snapshots of each model are plotted in Figures 2–7.

### 3.1 Analytic Expectations

We derived analytic expectations for the amplitudes of the various processes we believe are operating in the simulations, to gain an understanding of the results. Many processes are working together, often with the same net result (disruption of the clumps), which would greatly complicate the interpretation of the N-body simulations in the absence of analytic insight.

#### 3.1.1 The Effects of Global Tides

Global tides due to the spatial variation of the underlying large-scale dark halo potential provide an important mechanism for the disruption of substructures as they orbit through the dark halo. Analytical studies of the tidal effect on a clump show that, for a given clump orbit, the disruption efficiency depends on the ratio between the clump density to the mean density of halo within the pericentre of the clump orbit (e.g. King 1962; Binney & Tremaine 1987; Johnston, Hernquist & Bolte 1996). The tidal radius of the clump, during its disruption, is the Roche radius at which its gravity is equal to the tidal force exerted by the dark halo, which roughly scales as

$$r_{tide}\left(\frac{M_c}{M_{tot}}\right)^{1/3}R,$$
(3)
where $`R`$ is the pericentre distance. For our simulations, we expect that the clumps with density contrast $`\left(\rho _{cl}/\rho _0\right)=3`$ should suffer more efficient tidal disruption than those with $`\left(\rho _{cl}/\rho _0\right)=64`$. For a clump in a circular orbit at the half mass radius of the dark halo, the fraction of mass contained within the tidal radius is $`97\%`$ for $`\left(\rho _{cl}/\rho _0\right)=64`$ and $`74\%`$ for $`\left(\rho _{cl}/\rho _0\right)=3`$. Since the clumps in our models are each represented by a Plummer profile, which has a fairly constant density core (see Appendix A), while the dark halo is represented typically by a Hernquist profile, which has a cuspy central density profile, the minimum pericentre distance for the approximately constant density core to be disrupted is $`0.085`$ for $`\left(\rho _{cl}/\rho _0\right)=64`$ and $`0.72`$ for $`\left(\rho _{cl}/\rho _0\right)=3`$.

#### 3.1.2 Merging Effects

The merging cross section between a pair of identical spherical clumps has been studied by Makino & Hut (1997) numerically, using a variety of mass profiles for the clumps. They find that in a cluster of clumps with one-dimensional velocity dispersion $`\sigma _h`$ the number of merging events per unit time per unit volume, $`_m`$, is given by

$$_m=\frac{18}{\sqrt{\pi }}\frac{1}{x^3}n^2r_{cvi}^2\sigma _c_0(x)=\frac{18}{\sqrt{\pi }}\frac{1}{x^4}n^2r_{cvi}^2\sigma _h_0(x),$$
(4)
where $`x=\sigma _h/\sigma _c`$, with $`\sigma _h`$ and $`\sigma _c`$ being the one-dimensional velocity dispersion of the system of clumps (the halo) and the internal velocity dispersion of a clump respectively, $`n`$ is the number density of the clumps within the half mass radius of the dark halo, $`r_{cvi}`$ is the virial radius of the clump, and the dimensionless quantity $`_0(x)=\frac{Ax^2}{x^2+B},`$ with $`A=12`$ and $`B=0.4`$ for clumps with a Plummer profile. Written this way, one can see explicitly a strong dependence of the merging rate on the parameter $`x`$ (the ratio of the internal velocity dispersions of halo and clump).

The number of merging events per unit time, $`R_m`$, expected in our simulations can be estimated by the product of $`_m`$ and the volume within the half mass radius of the halo profile (equal to the half mass radius of the system of clumps). After some algebraic manipulations using equations (1) and (2) given in section 2.2 and the scalings in Appendix A, we obtain

$$R_m=\frac{4}{3}\pi r_{\frac{1}{2}h}^3_m0.47f^2_0(x)\frac{1}{t_{\frac{1}{2}hc}}.$$
(5)
Thus the number of merging events per halo crossing time is $`0.47f^2_0(x)`$, and depends fairly strongly on the clump mass fraction, $`f`$, but with only a weak residual explicit dependence on the parameter $`x`$ (many parameters are interdependent). In the limit of a high relative velocity compared to internal velocity dispersion, $`x^20.4`$, which is approximately valid for all our models, $`_0(x)12`$. This limit does not favour mergers, and results in the number of merging events per halo crossing time being $`0.05`$ for $`f=0.1`$, $`1`$ for $`f=0.4`$ and $`5`$ for $`f=1`$.

We can also estimate the merging timescale for a clump. For a given clump within the half-mass radius of the halo, the probability of it merging with another clump per unit time is $`_m/n`$. Thus the merging time scale is

$$t_m=n/_m=\frac{n_c}{f^2_0(x)}t_{\frac{1}{2}hc}.$$
(6)
It should be emphasized that the internal velocity dispersion ratio of the halo to the clump, $`x`$, itself depends on the clump mass fraction $`f`$, the clump number $`n_c`$ and the density contrast $`\rho _{cl}/\rho _0`$ by equation (2). We can see here that, given the values of the model parameters, the merging time scale is only very weakly dependent on the density contrast, which is instead a crucial factor for the efficiency of the global tidal stripping. Merging between clumps is an important effect for large values of the clump mass fraction, $`f`$, and for small numbers of clumps, $`n_c`$. In our models with $`f0.1`$, the merging effect is very small with $`t_m50t_{\frac{1}{2}hc}`$, while for larger $`f>0.4`$, merging can be important, and the minimum merging timescale is $`t_m1.7t_{\frac{1}{2}hc}`$ (model I).

The merging rate calculated above is the mean quantity at the initial half mass radius. Further, it does not take into account other dynamical effects, such as global tidal effects and close encounters, that can greatly decrease the calculated merging rate (Makino & Hut 1997). The effects of global tides can decrease the merging cross section by tending to tear apart a pair of clumps which, if isolated, would have merged. Similarly, the local tides from a third clump can also act to tear apart a merging pair of clumps. Furthermore, the subsequent nonlinear evolution in the presence of other more important dynamical processes can make the estimated merging timescale change quickly with time. We will see the limitation of this analytic estimate when we compare with the simulations below.

#### 3.1.3 Clump-Clump Heating

Disruption of clumps can also be caused by close encounters with other clumps. Similarly to the discussion given in section 2.3.2, we can estimate this clump disruption timescale as follows: For a given clump, the heating rate is

$$\dot{E}=\frac{4\sqrt{\pi }G^2m_c^3n\overline{r_c^2}}{3g\sigma _hr_{\frac{1}{2}c}^2},$$
(7)
where $`n`$ is the mean clump number density within the half mass radius of the dark halo, and all other parameters are as defined above. As discussed in Appendix B2, we adopt $`g3`$ and $`\overline{r_c^2}2r_{\frac{1}{2}c}^2`$. Thus the clump disruption timescale due to clump-clump heating is given by:

$$t_{cc}=0.03f^{4/3}n_c^{1/3}\left(\rho _{cl}/\rho _0\right)^{1/3}t_{\frac{1}{2}hc}.$$
(8)
Thus clump-clump heating is important for a large clump mass fraction, $`f`$, small clump number, $`n_c`$, and low density contrast, $`\rho _{cl}/\rho _0`$. We have $`t_{cc}11t_{\frac{1}{2}hc}`$ for model A, and $`t_{cc}1.5t_{\frac{1}{2}hc}`$ for model B.

However, we should keep in mind that the above calculation could overestimate the effects of clump-clump heating. The choice of $`\overline{r_c^2}`$ is uncertain. Further, the impulse approximation holds only for $`\sigma _c/\sigma _h1`$ if the impact parameter $`b`$ is chosen to be $`r_{\frac{1}{2}c}`$. If the encounter time is too long compared with the clump internal crossing time, particles in the clump may adiabatically respond to the encounter, hence the encounter could leave no net effects (Binney & Tremaine 1987). This could possibly explain the apparent inconsistency between our analytic estimates for the clump-clump merging effect and the clump-clump heating effect (equations (6) and (8)), in that they both become more important with large clump mass fraction and small numbers of clumps. From equation (2), we can see that $`\sigma _c/\sigma _h`$ also tends to be large for a large clump mass fraction and a small number of clumps (provided that $`\sigma _c/\sigma _h1`$). When $`\sigma _c/\sigma _h`$ increases above $`0.5`$, our prediction of the clump-clump heating effect in equation (8) becomes invalid, and the clump-clump heating effect will decrease quickly while the merging effect will increase continuously.

The above arguments can be understood in terms of three different regimes depending on the value of the encounter speed, basically equivalent to the role of $`\sigma _h`$. (1) In the case that the encounter speed $`V`$ is very high, the two clumps will not merge and the clumps are just heated somewhat during the fast encounter. (2) In the case that the encounter speed is just slightly reduced, it can be envisioned that the chance for merging is slightly increased, and the heating effect is also slightly enhanced since the duration of the encounter is slightly longer. (3) In the case that the encounter velocity is very low, the merging rate is greatly increased, and the heating effect can be inhibited, since if the duration of the encounter is long enough compared to the internal crossing timescale in the clump, the clump just adiabatically responds during the encounter and returns to the initial state after the encounter.

#### 3.1.4 Dynamical Friction

Dynamical friction can drive clumps to the centre of the dark halo, where they can be tidally disrupted more easily. The slowing down of the clumps also can increase their merging cross-section. At the same time this process can extract angular momentum from the clumps, and transfer it to the dark halo. Since Chandrasekhar (1943) introduced the concept of dynamical friction, namely that an object moving through an infinite and homogeneous medium made of small mass particles suffers a drag force, there have been many studies on this subject by numerical simulations (e.g. White 1978, 1983; Lin & Tremaine 1983; Bontekoe & van Albada 1987; Zaritsky & White 1988; Hernquist & Weinberg 1989; van den Bosch et al. 1999) and analytical methods (e.g. Tremaine 1981; Tremaine & Weinberg 1984; Weinberg 1989; Maoz 1993; Domínguez-Tenreiro & Gómez-Flechoso 1998; Colpi, Mayer & Governato 1999; Tsuchiya & Shimada 2000). An overview of past work is given by Cora, Muzzio & Vergne (1997). In spite of the difficulties encountered in the study of dynamical friction, many authors find Chandrasekhar’s formula is a remarkably good approximation (e.g. Velazquez & White 1999).

Adopting Chandrasekhar’s formula the deceleration of a clump is then

$$\frac{d𝐯_M}{dt}=\frac{𝐯_M}{t_{df}},$$
(9)
where

$$t_{df}^1=16\pi ^2\mathrm{ln}\mathrm{\Lambda }G^2m_p(M_c+m_p)\frac{_0^{v_M}f(v)v^2𝑑v}{v_M^3},$$
(10)
and $`m_p`$ is the mass of a background halo particle, $`M_c`$ is the mass of the clump, $`v_M`$ is the velocity of the clump, $`f(v)`$ is the phase-space number density of the background halo medium, and $`\mathrm{\Lambda }`$ is the ratio between the maximum and minimum impact parameters $`b_{max}/b_{min}`$. Using the values at the half mass radius, we can immediately see that the characteristic dynamical friction time,

$$t_{df}\frac{n_c}{f(1f)}t_{\frac{1}{2}hc}\frac{1}{M_c(1f)}t_{\frac{1}{2}hc}$$
(11)
scales inversely with the mass of a clump ($`f/n_c`$, with the total mass normalized to unity as here).

In applying these formulae to our simulations with virialized initial conditions, $`2T/W1`$, and for which the dark halo is taken to follow the Hernquist profile, we have chosen to use the phase space distribution function given by Hernquist (1990), hence ignoring the angular momentum dependence of the distribution function. The phase-space density of the dark halo is then as detailed in Appendix C. For a clump initially on a circular orbit at the half mass radius, the dynamical friction time calculated as in Appendix C is

$$t_{df,circ}=\frac{0.16n_c}{f(1f)}t_{\frac{1}{2}hc}.$$

Analyses of the orbital eccentricities of substructure in spherical potentials have found that for isotropic distribution functions the typical orbital eccentricity is $`0.6`$ (van den Bosch et al. 1999). With this value the dynamical friction time is decreased by a factor of up to 2, and the typical dynamical friction time for our computed models should then be

$$t_{df}0.08\frac{n_c}{f(1f)}t_{\frac{1}{2}hc}.$$
(12)
Thus dynamical friction is more important for a large clump mass fraction (provided $`f<0.5`$) and a small number of clumps, giving a large mass for each clump. This dependence on $`f`$ and $`n_c`$ in a general sense is consistent with that of clump-clump merging in all parameter ranges and that of clump-clump heating in some restricted parameter ranges. For our models with $`n_c=80`$, dynamical friction is not important, with $`t_{df}>20t_{\frac{1}{2}hc}`$. For our models with $`n_c=20`$, dynamical friction is not important for $`f=0.1`$, with $`t_{df}18t_{\frac{1}{2}hc}`$, but for $`f=0.4`$ it becomes important, with $`t_{df}6.7t_{\frac{1}{2}hc}`$. For our models with $`n_c=5`$, dynamical friction is important, with $`t_{df}8.4t_{\frac{1}{2}hc}`$ for $`f=0.05`$ and $`t_{df}1.7t_{\frac{1}{2}hc}`$ for $`f=0.4`$.

It should be emphasized that efficient dynamical friction leads to significant angular momentum loss from the clumps, and thus can drive clumps to the centre, where global tidal effects are stronger. Thus the net efficiency of tidal stripping by the global potential includes a dependence on the clump mass fraction and number of clumps similar to that of dynamical friction.

### 3.2 Comparison with Simulations

It should be noted that the above analysis is mainly applicable to virialized models. For the pre-collapse models, the disruption of clumps due to the collapse process itself dominates. For the Plummer dark halo profile, the details of the calculation of dynamical friction from section 3.1.4 are not applicable, but the scaling should be the same as in the Hernquist models. Thus, for ease of comparison, the dynamical friction time is normalized to the crossing time scale at the half mass radius, as are all timescales.

As noted earlier, all models with the same number of clumps have the same orbital characteristics for the clumps, with the exception of models B and E. Thus it should be the different choices of our three free parameters – the density contrast $`\rho _{cl}/\rho _0`$, the clump mass fraction $`f`$ and the number of clumps $`n_c`$ – that are responsible for differences in the evolution of the different models.

The most important uncertainty in the comparison of the theoretical predictions with our simulation results comes from two factors that enter the determination of the theoretical clump-clump heating timescale, namely the normalization quantity, $`\overline{r_c^2}`$, and the use of the impulse approximation even for relatively large values of $`\sigma _c/\sigma _h`$. However, there are also uncertainties in the analysis of our simulations. Our code is not designed to follow the merging events that occur throughout the simulations, and as mentioned above, the ‘debris’ includes both particles genuinely removed from clumps and orbiting freely in the global potential and the remnants of merging clumps. Thus in our analysis, both heating and merging produce debris. As a rough estimate of the mean disruption timescale of a clump in the simulations, we simply use the time required for the disruption of half of the clump mass.

Table 2 lists, for the virialized models only, all the timescales, in units of $`t_{\frac{1}{2}hc}`$, calculated from our analytic expressions, together with the disruption time of half of the mass in clumps measured from our simulations. The reader should note that the relative importance of the processes has not be derived from the simulations explicitly. Fig. 8 shows the dependence of various processes on the mass of the clump, and also on the mass fraction in clumps, $`f`$, the latter denoted by the symbol size.

The discussion of each virialized model is given in turn below:

(1) For model A, since $`\sigma _c/\sigma _h=0.21`$, the impulse approximation is valid in our calculation of the clump-clump heating timescale $`t_{cc}`$. Indeed, the dominant dynamical process is clump-clump heating and the analytic prediction is quite consistent with the simulation result.

(2) For model B, again clump-clump heating is the dominant dynamical process. However our analytic estimate of $`t_{cc}`$ is a factor of 3 below the simulation result. This can be understood since $`\sigma _c/\sigma _h=0.31`$ and thus we are in the regime of marginal applicability of the impulse approximation, and we may have underestimated $`t_{cc}`$.

(3) For model C, we can see $`\sigma _c/\sigma _h=0.46`$ and thus our underestimate of $`t_{cc}`$ could be significant, possibly by as much as the factor of 30 needed for consistency with the simulation result. It is also possible that $`t_{merge}`$ is underestimated by the factor of 2 discrepancy with the simulation result. The dominant processes appear to be merging and clump-clump heating.

(4) For model G, $`\sigma _c/\sigma _h=0.34`$, and again we find our apparent underestimate of $`t_{cc}`$ is about a factor of 3. The dominant process is apparently dynamical friction or possibly clump-clump heating. Note that dynamical friction acts to drive clumps to the centre, where they will be tidally disrupted more easily.

(5) For model H, $`\sigma _c/\sigma _h=0.5`$, and again we have an apparent underestimate of $`t_{cc}`$ by about a factor of 20, and an underestimate of $`t_{merge}`$ by a factor of 2. The dominant process appears to be dynamical friction.

(6) For model I, $`\sigma _c/\sigma _h=0.74`$, sufficiently close to unity that our underestimate of $`t_{cc}`$ could be very significant. We can simply assume that there is no clump-clump heating for this case (adiabatic encounters). The dominant physical process is then merging of the clumps. As can be seen from Figure 2, in this model all the clumps finally merge into two large systems. Since in this case the ‘debris’ is mainly in these two merging remnants, we do not calculate the disruption time in the simulation.

(7) For model M, $`\sigma _c/\sigma _h=0.34`$, and we again assume our underestimate of $`t_{cc}`$ is about a factor of 3. As noted above, our analytic estimates of the dynamical friction timescale are not applicable for Plummer profiles, but as a rough estimate we can adopt the value from the equivalent Hernquist profile model G. The dominant processes are clump-clump heating and possibly dynamical friction.

(8) For model N, as above, we adopt the dynamical friction timescale from the corresponding Hernquist profile model H. The velocity dispersion ratio is $`\sigma _c/\sigma _h=0.5`$, and again the anaytic expressions have apparently underestimed $`t_{cc}`$ by about a factor of 20, and $`t_{merge}`$ by about a factor of 2. The dominant process is dynamical friction.

(9) For model O, again the ratio of velocity dispersion is sufficiently close to unity ($`\sigma _c/\sigma _h=0.74`$) that we assume there is no clump-clump heating. The only possible dominant physical process is merging of clumps, and the analytic estimate of $`t_{merge}`$ is apparently underestimated by a factor of 4.

(10) For model S, $`\sigma _c/\sigma _h=0.20`$. The dominant disruption processes are global tidal effects and clump-clump heating.

(11) For model T, $`\sigma _c/\sigma _h=0.30`$. Again we assume our underestimate of $`t_{cc}`$ is about a factor of 3. The dominant disruption processes are global tidal effects and clump-clump heating.

(12) For model W, $`\sigma _c/\sigma _h=0.43`$, and the analytic underestimate of $`t_{cc}`$ is at least a factor of 10. The only dominant process is dynamical friction.

(13) For model X, $`\sigma _c/\sigma _h=0.79`$, again sufficiently high that we ignore the clump-clump heating effect. Again we assume our underestimate of $`t_{merge}`$ is about a factor of 2. The possible dominant processes are dynamical friction and merging.

From the above, we can see that the modification to our analytic expression for the clump-clump heating timescale $`t_{cc}`$ should be performed in a consistent way: for $`\sigma _c/\sigma _h<0.20`$, no modification is needed; for $`\sigma _c/\sigma _h0.3`$, we multiply by a factor of 3; for $`\sigma _c/\sigma _h0.4`$, we multiply by a factor of 10; for $`\sigma _c/\sigma _h0.5`$, we multiply by a a factor of 30; for $`\sigma _c/\sigma _h>0.75`$, we ignore the clump-clump heating effect (i.e. make the timescale infinite).

From the above we can see that the merging timescale $`t_{merge}`$ should be multiplied by a factor of 2 for all values of the velocity-dispersion ratio. From the comparison of the analytic theoretical predictions with our simulation results, we can see that with some reasonable and consistent modifications to our theoretical estimate of the timescale for clump-clump heating and merging, we can reach rough consistency with the estimates measured from our simulations.

This allows us to develop an understanding of the dependence of the various processes of global tidal stripping, clump-clump heating, clump-clump merging and dynamical friction on the number of clumps, mass fraction of clumps and density contrast of a clump in the halo. In general, a decrease in the number of clumps and an increase in the clump mass fraction can increase the velocity dispersion ratio of clump to halo, and can further increase directly or indirectly the efficiency of all four processes, if the ratio is much smaller than the order of unity. For larger values of this velocity dispersion ratio, the efficiency of the clump-clump heating process drops quickly, while the merging process becomes more important, and the dynamical friction becomes more important as long as clump mass fraction is less than 50%. A decrease in the density contrast between clump and halo can enhance the disruptive effects of global tides, and of clump-clump heating (if the velocity dispersion ratio of clump to halo is much less than unity).

## 4 The Angular Momentum

All the simulations have good conservation of total angular momentum and of total energy. The total angular momentum is conserved relative to the initial origin and to the centre of the total mass. When significant substructure is involved, as here, the calculation of angular momentum can be complicated by two factors. The first is that the smooth dark matter component and the clump component do not always share the same centre of mass. The second is that the centre of mass of the smooth dark matter component does not always correspond to its highest density point – which is usually taken to be the centre of the relevant system, the galaxy or the cluster of galaxies – due to the fact that the dark halo density profile to be simulated does not decrease with radius fast enough, e.g. for Hernquist models, $`\rho r^4`$, and the simulation can only allow a limited number of particles to be distributed to infinity. The first factor is possibly non-trivial, since though using many clumps can make their centre of mass more closely coincide with that of the smooth component, there do exist physical environments, such as clusters of galaxies, in which indeed a small number of clumps are embedded in a smooth component. The second factor however is artificial; the few particles in the outermost regions have a distribution that is more extended and asymmetric, and thus have more weight in the determination of the centre of mass. This problem could possibly be avoided by making the density profile decrease more rapidly with radius, e.g., applying a cutoff radius. Indeed, our simulations show that the coincidence between the centre of mass of the dark halo and the highest density point is better for the Plummer profile than for the Hernquist profile.

Table 3 lists the different measures of the angular momentum contents of our computed models, with the exception of the models with a clump mass fraction $`f=1`$. The quantity $`\delta _{CL}`$, in Column 7, is defined as

$$\delta _{CL}=\frac{_{k=1}^{n_c}|J_{zk}^iJ_{zk}^f|}{_{k=1}^{n_c}|J_{zk}^i|},$$
(13)
where $`J_{zk}^i`$ and $`J_{zk}^f`$ are the initial and final angular momentum, calculated relative to the highest density point of the dark matter, of the $`k`$th clump. It should be noted that any direct comparison between the Plummer models and the Hernquist models should use the halo crossing time as the time unit; the final times here are fixed in these units for each of the Plummer models and the Hernquist models.

We can see from the entries in this Table that shifting the origin from the centre of mass of the total system, to the highest density point of the dark matter, can produce apparent non-conservation of the total angular momentum; this arises since the highest density point of the dark matter defines a non-inertial frame.

A real physical change in angular momentum with respect to the centre of mass of the total system occurs for each of the clumps and the dark matter (compare the changes with the estimated errors). This ensures that the angular momentum exchange between the clump component and DM component is real (Barnes 1988). However, from the observational point of view, or in cosmological simulations, the angular momentum is rather calculated relative to the highest density point of the DM. Thus the question we want to ask is whether the angular momentum of the clumps changes significantly when measured in this way. This remains non-trivial since, for example, if a clump is extremely far from the centre, and thus it contains most of the angular momentum of the clumpy component, no matter how efficient is the angular momentum transfer from the remaining clumps to the dark matter, the angular momentum change of the total clump system will be very small. Thus we must address not only the total angular momentum change of the clumps but also the angular momentum change of each individual clump.

The parameter $`\delta _{CL}`$ is useful in quantifying the mean angular momentum change of an individual clump. Figures 9 and 10 show the binding energy and angular momentum for each individual clump. It should be noted that in our simulations in which the smooth (dark matter) component has a Plummer profile, the initial direction of the angular momentum vector of the clumpy component is in the opposite direction to that of the DM component. This is simply a result of the random assignment of the initial velocity, leading to negative values of the angular momentum for (several of) the clumps at large radius. From the plots of models M and P (Figure 9 and Figure 10 respectively) we can see this fluctuation leads to a lower initial specific angular momentum than in the Hernquist models. Thus the total angular momentum of the clumps can be changed drastically compared to its initial value, even if the change of angular momentum of an individual clump is slight. Thus we conclude that $`\delta _{CL}`$ is a more robust parameter to measure changes in angular momentum than is the percentage change of the total angular momentum of the clumps, especially for the Plummer models.

From the entries in Table 3 we can see that $`\delta _{CL}`$ for all the collapse models is much larger than most of the virialized models. By checking individual clumps, as illustrated in Figures 9 and 10, we find that clumps close to the centre have significant angular momentum loss, while the clumps at large distance can have significant angular momentum change – gains and losses – in collapse cases.

For the virialized models A, S, B, G, W and T, $`\delta _{CL}`$ is small, while for the remaining virialized models $`\delta _{CL}`$ is large. We list the dynamical friction time and the time for half of the mass in clumps to be disrupted, in units of $`t_{\frac{1}{2}hc}`$, for these models, in increasing order of the value of $`\delta _{CL}`$, in Table 5. We can see that there are two factors that determine the angular momentum change $`\delta _{CL}`$: the dynamical friction time and the disruption time. Generally, long dynamical friction times and short disruption times produce small values of $`\delta _{CL}`$. As discussed in section 3.1.4, dynamical friction acts to drive clumps to the centre, with loss of angular momentum of the clumps. On the other hand, efficient disruption of a clump, especially due to processes unrelated to dynamical friction such as clump-clump heating or global tidal effects, helps to decrease the mass of the clump, which slows down the dynamical friction process, and thus can decrease the angular momentum loss. In the extreme case, a clump will no longer suffer dynamical friction after being disrupted completely. Thus the angular momentum change is most significant for those massive clumps moving slowly, close to the centre, that are dense enough to withstand disruption by tides.

For the collapse models, the large value of $`\delta _{CL}`$ is caused by not only dynamical friction, but also by dynamical mixing during the post-collapse relaxation, which can redistribute both the binding energy and the angular momentum (Quinn & Zurek 1988). We plot the Lindblad diagram for the smooth DM component for models K and D in Figure 11. The most bound $`90\%`$ of the DM particles are divided into nine equal size bins, sorted by the binding energy of each particle at the end of the simulation. We plot the vector angular momentum $`J_z=|\mathrm{\Sigma }j_z|`$, and the scalar angular momentum $`J_s=\mathrm{\Sigma }|j|`$, for each bin. We can see that for all the collapse models, regardless of the clump mass fraction or number of clumps, they show the same behaviour, in that inner, more-bound material loses angular momentum $`J_z`$, while outer, less-bound material gains angular momentum $`J_z`$. The scalar angular momentum $`J_s`$ increases for each bin. This result is consistent with previous results for the redistribution of angular momentum in protogalaxies (Quinn & Zurek 1988). The clumps close to the centre are disrupted quickly during the collapse, but still suffer significant angular momentum loss, both while still bound in clumps and after disruption. This material becomes more bound to the halo. For those clumps that survive the collapse, subsequent disruption and dynamical friction effects occur, as discussed above for the case of an initially virialized system. In addition, the final dark matter distribution is triaxial (see section 6.2 below) and exerts a torque that helps transfer angular momentum between clumps and the dark matter, and between the inner and outer regions (resulting in the angular momentum loss from the inner debris).

In the context of galaxy formation, the quantity of interest is the total angular momentum change, over a Hubble time, of the baryonic component (here given by the clumps). However as can be seen from Table 3, the total angular momentum $`J_z`$ of all the clumps does not change as much as would be expected from the parameter $`\delta _{CL}`$, except for model Y and all the Plummer models. For the Plummer models, as we discussed above, the large change in the total angular momentum of the clumpy component is due to the particular random choice of the initial (low) specific angular momentum of this component. As an upper limit, we may adopt the parameter $`\delta _{CL}`$ as the estimate of the angular momentum change for the clump system. Our simulations show that the total angular momentum of the clumpy component does not change significantly in a Hubble time for a small clump mass fraction and a large number of clumps, but does change significantly (quantified by $`\delta _{CL}`$ in Table 3 being greater than $`30\%`$) in models with few and massive clumps; compare model G with model H, or model S with model T, or model B with model X. Such a high amplitude of angular momentum loss from proto-disk material cannot allow a spiral to form that matches observation. Thus we conclude that spiral galaxies can be formed through the accretion of many, but small, clumps, that conserve angular momentum in the process, but cannot be formed in an environment that contains only a few, but large, clumps (cf. Silk & Wyse 1993).

## 5 The Kinematics

In all of the collapse models the kinematic properties for the dark halo component are similar, regardless of the different parameter values. From Figure 12, we can see that after collapse the radial velocity dispersion $`\sigma _r`$ increases at all radii, but the increase is larger in the central regions than at large radius. The other components of the velocity dispersion tensor, $`\sigma _\varphi `$ and $`\sigma _\theta `$, increase less than does $`\sigma _r`$, with $`\sigma _\varphi `$ slightly larger than $`\sigma _\theta `$, and they do not increase significantly at large radius. Thus the final ratios $`\sigma _\theta /\sigma _r`$ and $`\sigma _\varphi /\sigma _r`$ decrease with increasing radius. However, at the centre, the velocity ellipsoid approaches isotropy, i.e. $`\sigma _\theta \sigma _\varphi \sigma _r`$. The rotation velocity, $`V_{rot}`$, increases, while $`V_{rot}/\sigma _r`$ decreases. $`V_{rot}/\sigma _\varphi `$ decreases during the collapse in the central regions but is fairly constant in time at large radius.

The Lindblad diagram for the dark halo component is shown in Figure 11, and shows that for the most bound energy bin the scalar angular momentum $`J_s=\mathrm{\Sigma }|j_z|`$ increases, despite the fact that the amplitude of the vector angular momentum $`J_z=|\mathrm{\Sigma }j_z|`$ decreases. The different behaviours of the scalar and the vector angular momenta for the most bound energy bin can be explained by noting that individual particles can gain angular momentum in addition to kinetic energy, but can be deflected isotropically in direction during the post-collapse relaxation process, thus the vector angular momentum of all the particles in that energy bin decreases. As for the virialized models, there are no significant changes for the dark halo component.

The kinematic properties of the clump debris for each model are somewhat different, but regardless of the parameter values, the velocity ellipsoid properties are similar to those of the dark halo component, especially in the central regions. Figure 13 shows the same quantities for the debris as does Figure 12 for its dark halo. In the outer regions, the debris is largely still composed of many streams formed from the disrupted clumps. To illustrate, Figure 14 plots $`v_x`$ and $`v_y`$ at the final time for the clump particles within a box of diameter of $`3,`$ located at the coordinate (4,4,0) (beyond the half-mass radius) for model U. The different symbols indicate membership of different initial individual clumps; in this example the box contains several streams, which are composed of eight different disrupted clumps.

## 6 The Morphology

### 6.1 The Density Profile

The density profiles of the disrupted clumps for all the collapse simulations follow the profiles of the dark matter, the surface density profile of which can be best fit by the $`R^{1/4}`$ law. This result is consistent with the early studies of van Albada (1982), which demonstrated that no matter how clumpy the pre-collapse initial condition is, the final profile is consistent with the de Vaucouleurs $`R^{1/4}`$ law, provided the collapse factor is large enough (the virial ratio $`2T/W0.1`$). However, in the virialized simulations, the density profile of the debris can be different from that of the dark matter. Figures 15 and 16 show the density profile and surface density profile for the debris (normalized to their values at the debris half-mass radius). The solid lines correspond to the collapse cases, while the dashed lines correspond to the virialized cases. The density profile for model H shows a peak at some distance from the centre, due to the presence of a merger remnant there. We can see that, once normalized, the profiles of the debris in all the models are very similar. At the central regions, the density profile of the virialized models is slightly shallower than that of the collapse models. At the normalized distance $`r=0.1`$, i.e. at one-tenth of the half-mass radius, $`\rho r^2`$ for collapse models while $`\rho r^{1.7}`$ for virialized models. At the very centre, due to the different resolutions for our different models, we can only infer that the central density profile is cuspy.

The models with the Plummer profile for the dark halo behave somewhat differently. Figure 17 shows the final density profiles of the debris and of the dark matter for models M and N. These demonstrate that the central density of the debris can be larger than that of the dark halo component, and the central density profile of the debris is cuspy and does not have a (approximately) constant density core, in contrast to the dark halo profile. The initial and final density profiles of each component for models M and N are shown in Figure 18. Since the distribution of the clumpy components initially follows the same Plummer profile as does the dark matter, the initial density profile including all components is approximately a Plummer profile. The final density profile, including all components, shows a cuspy structure at the central region, for example at $`r0.1`$, model M has central profile $`\rho r^{1.1}`$ and model N has central profile $`\rho r^{1.7}`$. A larger clump mass fraction (model N) leads to a steeper central density profile and a higher normalization, compared to the results with a lower clump mass fraction (model M). This is because in the former case dense and massive clumps can move to the centre through dynamical friction without being significantly disrupted, and thus can contribute to the mass and density profile at the centre. This result has the application, in CDM cosmology, that the merging into a lower density parent halo of a number of higher density subhalos can produce a final cuspy and dense halo. Furthermore, how cuspy and dense the final halo is depends mainly on the number of subhalos in the parent halo and the density contrast of a subhalo to the parent halo, both of which are related to the density fluctuation power spectrum index, $`n`$, for $`P(k)k^n`$. This provides an explanation for the central cuspy density profiles of halos found in cosmological numerical simulations, as the consequence of the self-limiting merging/accretion of small, dense halos, as suggested by Syer and White (1998).

It is well-known that the surface density profile of a cD galaxy in its outer regions is systematically above the $`R^{1/4}`$ law that fits the more central parts. Recent N-body simulations of groups of 50 galaxies studied by Garijo, Athanassoula & García-Gómez (1997) show that a giant central galaxy can be formed if there is a initial seed galaxy at the centre. Further, these authors find that their virialized models produce a giant central galaxy with such a cD envelope. Previous studies of the formation of cD galaxies, using non-cosmological N-body simulations of isolated clusters of galaxies, (Richstone & Malumuth 1983; Bode et al. 1994) concentrated on identifying the relative roles of stripping, merging and collapse. Can a cD galaxy be formed by the galaxies that were disrupted at the cluster centre or by galaxies that merged with the central galaxy? Does the formation of a cD galaxy happen before, or after, the collapse of the cluster? Dubinski (1998) explored the formation of brightest cluster galaxies through the merging of several massive galaxies, in the context of hierarchical clustering cosmologies, but still failed to create the characteristic extended envelope of cD galaxies. One might expect that a larger clump mass fraction and a small number of clumps can lead to efficient merging or disruption of galaxies, which are plausibly important processes to form cD galaxies. Though our simulations do only explore a limited parameter space, they include a simulation with only five clumps with a mass fraction of 40%, and also one with only twenty clumps with a mass fraction of unity. Neither formed a cD. All our giant central galaxies can be fit by an $`R^{1/4}`$ law very well, over a range of more than 4 orders of magnitude in surface density. We speculate that collapse and small-scale structure alone cannot produce cD galaxies. It is plausible that the continuous later infall of large galaxies could provide a mechanism to form cD galaxies.

### 6.2 Triaxiality

CDM-dominated cosmological N-body simulations show the protohalos are usually triaxial (Warren et al. 1992; Dubinski & Carlberg 1991) with the tendency to be prolate, and with the mean ratio of minor to major axes, $`c/a`$, approximately 0.5 at small radius. The study of Tremblay & Merritt (1995; Merritt & Tremblay 1996) of the intrinsic shapes of elliptical galaxies found that bright ellipticals ($`M_B<20`$) are systematically rounder than faint ellipticals, with the ratio of minor to major axes, $`c/a`$, peaked at 0.75 for bright ellipticals, while for faint ellipticals $`c/a`$ is peaked at 0.65. Small ellipticals also tend to be flattened by rotation (Davies et al. 1983; Rix, Carollo & Freeman 1999) rather than by anisotropic stellar velocity dispersions. Studies of the collapse of a proto-galaxy with some initial triaxiality have explored the possibility that angular momentum can be transferred from the inner halo to the outer parts (Subramanian 1988; Curir, Diaferio & de Felice 1993; Curir & Diaferio 1994). Our simulations study the evolution of a system of clumps embedded in a smooth dark halo with an initial oblate profile, and thus we do not address the question of the role of any initial triaxiality, but rather investigate what parameters determine the triaxiality of the final system.

We first investigate the level of flattening of various components that can be due to rotation, by comparing the ellipticity with the ratio of rotation to velocity dispersion. The ellipticity, $`ϵ`$, is calculated from the projected isodensity contour of the simulation particles from sixteen different viewing angles. The central line-of-sight velocity dispersion, $`\sigma _0`$, and rotational streaming velocity, $`v_0`$, are also calculated for these viewing angles, and the ratio $`v_0/\sigma _0`$ formed. In Figure 19 we plot the location of each viewing angle on the ($`v_0/\sigma _0`$,$`ϵ`$) plane for models J, G, U and K. Since the morphology of the debris generally follows that of dark halo, here we just show the dark halo measured in this way. From Figure 19 we can see for the collapse cases the value $`v_0/\sigma _0`$ is decreased significantly and the system becomes anisotropic. For model U, with low density contrast between the clumps and the halo, the ellipticity increases, while for model K, with large clump mass fraction, the ellipticity decreases. The virialized cases are actually initially marginally out of equilbrium, and generally evolve such that the ellipticity and $`v_0/\sigma _0`$ are decreased, to make the system rounder and more isotropic (remember these simulations are constructed to be initially close to isotropic and flattened by rotation).

The three-dimensional shapes can be calculated from the moments of inertia, giving the triaxiality for both dark halo and debris. As we discussed in section 2.4, this technique for calculating the axial ratio systematically underestimates the flattening. For example, the initial flattening of the dark halo in all simulations is $`c/a0.53`$, while the calculation from the moment of inertia gives $`c/a0.63`$. However, since we are interested more in trends in the dependence of final triaxiality on the initial conditions, such a systematic overestimate is not a problem. Table 6 summarizes the values of $`b/a`$ and $`c/a`$ in the simulations, with as usual $`a>b>c`$ where a, b and c are calculated from the eigenvalues of the moment of inertia tensor within the half mass radius as described in section 2.4.

From Table 6, and Fig. 20 we can see that:

(a) all the collapse models show some level of triaxiality for both the dark halo component and the debris component. Cases with a large clump mass fraction ($`f=0.4`$) are slightly rounder than the cases with small clump mass fraction ($`f=0.1`$). For the completely clumpy cases ($`f=1`$) the final shape of the debris is very spherical, with the exception of the collapse model with a initial Plummer profile (model R). The triaxiality for the low density contrast cases, such as models U and V, is more significant than that of high density contrast cases, such as models J and K. The triaxiality for models with a small number of clumps, such as models Y and Z, is less significant than that for those with a large number of clumps, such as models J and K, or D and E. Also we should note that, for our choice of initial configuration, the clump mass fraction is weakly correlated with the value of the virial ratio, $`2T/W`$, for the collapse models. Actually the real collapse factor for all the collapse models is almost a constant. For the pure clump models with $`f=1`$, calculating a modified virial ratio (denoted $`2T^{}/W^{}`$) by excluding the clump internal kinetic energy and potential energy, we have $`2T^{}/W^{}=0.13`$. So, for our models, it is the level of clumpiness, rather than the collapse factor, that determines the triaxiality of the final systems. Large values of the clump density contrast, large clump mass fractions, and small numbers of (massive) clumps all indicate that clumpiness is important; a higher level of ‘clumpiness’ makes the final system rounder. Clearly, the collapse factor will also affect the triaxiality of the final system.

(b) for the virialized cases, the final shape of the dark matter component is oblate with $`c/a>0.8`$, slightly rounder than its initial configuration; the shape of the debris component is close to oblate, with $`b/a>0.9`$ and $`c/a>0.7`$.

As we can see from the virialized cases, which is an extreme case of approximately zero collapse factor, the final system tends to be oblate or round. This result implies that the generation of triaxiality requires thin box orbits, or radial orbits, that get very close to the centre. In this case the collapse is mostly radial and creates a lot of triaxiality; more clumpiness disrupts these radial orbits. This suggests that giant ellipticals, which tend to be rounder than small ellipticals, may be formed in a more clumpy environment than are small ellipticals, but this inference is complicated by the fact that the small ellipticals tend to be oblate, with significant dissipation being implicated (e.g. Wyse & Jones 1984; Rix et al. 1999).

### 6.3 The Substructure in the Outer Halo

The existence of a population of intracluster stars in the core region of the Virgo cluster, detected through planetary nebulae (Méndez et al. 1997) and red giant stars (Ferguson, Tanvir & von Hippel 1998) indicates that the stripping of stars from cluster galaxies, through the multiple disruption processes that we discussed above, is a reality during the formation of clusters of galaxies. In addition, M87, the giant elliptical in the Virgo cluster, exhibits a broad, diffuse and extended ‘fan’ in surface brightness at a projected distance of around 100 kpc from its centre. The accretion of a small stellar satellite with $`0.1\%`$ of the mass of M87, with small impact parameter, can allow the formation of a surface brightness feature like this (Weil, Bland-Hawthorn & Malin 1997). Even in our Galaxy, large overdensities of A-stars (Yanny et al. 2000) and RR Lyrae stars (Ivezic et al. 2000) are observed in the halo at distance of up to 60 kpc from the Galactic centre, through the analysis of Sloan Digital Sky Survey imaging data (Yanny et al. 2000). This structure may be tidal debris from the Sagittarius dwarf spheroidal galaxy (Ibata et al. 2001). The observation of substructure in the outer halo can be useful in constraining the formation history of galaxies. From Figures 21 and 22 we can see that in the outer halo region of our simulations there are indeed prominent streams from the disrupted clumps, seen in surface density. These unrelaxed streams can account for the over-density above the $`R^{1/4}`$ law seen at the very outer regions in Figure 16.

A detailed study of the substructure is a mathematical challenge needing a serious treatment, but as a first step we adopt a simple approach to explore its presence and evolution. To quantify the lumpiness of substructure as it might be observed, the concept of the filling factor of disrupted clumps - the fraction of sky containing one or more streams - was introduced by Tremaine (1993) and used by many authors subsequently in studying the evolution of tidally disrupted stellar systems (Johnston 1998; Helmi & White 1999). The shortcoming of the filling factor as a measure of streaming is that it can become significantly larger than unity as streams wrap around the sky.

Here we introduce a different measure of inhomogenities in the debris, restricting our analysis to the spatial part of phase space (the complications in velocity space are illustrated in Figure 14). We define the quantity $`\eta `$, defined as the fraction of the debris which is in regions with overdensity above some critical value, or $`\delta n/\overline{n}>\delta _n`$, restricting the calculation to regions beyond a critical radius $`r_{cr}`$ projected onto the celestial sphere. At a given time $`t`$, the critical radius $`r_{cr}`$ is chosen to be the radius where $`t_c(r_c)=c_tt=2\pi r/v_c(r)`$, where $`v_c(r)`$ is the circular velocity at radius $`r`$, and $`c_t`$ is a free parameter of the order of unity. This critical radius can be regarded as the division between the inner, smooth debris where mixing is complete, and the outer regions where mixing is incomplete. Though the assumption is over-simplified, a visual check of Figure 21 justifies a choice of $`c_t2`$. For collapse cases, we have approximately $`r_{cr}=0.7,2.3,4.3`$ respectively for $`t=14,56,126`$, while for virialized cases $`r_{cr}=0.6,2.1,4.1`$, for times $`t=14,56,126`$.

It should be noted that we do not have any ‘baryons’ that were not initially in clumps, so that all the unbound clump particles are ‘debris’; in the interpretation of the evolution of substructure in a galaxy like the Milky Way, the entire field halo would be debris, and we are quantifying the structure in that debris.

We simply divide the sky into a number of cells of equal solid angle. The average number of debris particles in each cell is $`\overline{n}`$. If the cell has more than $`(1+\delta _n)\overline{n}`$ debris particles, then it is counted as a region of overdensity. If $`\overline{n}16`$, the overdensity of $`\delta _n=0.75`$ can be used to indicate the existence of streaming structure in that cell, since the density fluctuation is then greater than 3$`\sigma `$, with $`\sigma =\sqrt{\overline{n}}`$. Thus an accurate calculation of $`\eta `$ requires that both the number of cells and the number of debris particles be large enough. However, this condition is not satisfied in many of our simulations. Thus here we can only give a rough estimate of $`\eta `$ for the simulations with large clump mass fraction, using $`\delta _n=0.75`$ and $`c_t=2`$.

Table 7 summarizes the values of $`\eta `$ at different times for our simulations, with the conclusion that for our models $`\eta `$ is more or less a constant, with a value varying only from about $`10\%`$ to $`20\%`$. By checking the number of distinct overdense regions on the celestial sphere used in calculating $`\eta `$, and also by a visual check of the snapshots of the models for 5, 20, and 80 clumps, we can see that, not surprisingly, the number of streams correlates with the number of clumps. The disruption of clumps is more efficient for collapse cases, resulting in streams that are more radially aligned compared to those formed in our virialized models. The observations of streams in phase space in the outer halos of galaxies are thus potentially important in constraining the formation history of galaxies, especially in determining the abundance of accreted substructure of different sizes and densities, through comparison with simulations.

## 7 Summary

In this paper, we use N-body simulations to study the evolution of clumps, or substructure, embedded in non-dissipative, isolated dark halos. We ran 26 different combinations of 4 independent free parameters, namely the clump mass fraction, the number of clumps, the density contrast of clump to halo and the initial virial ratio. We also investigated two different density profiles for the smooth dark halo component, the cuspy Hernquist profile and the (approximately) constant-density core Plummer profile. For the initial virialized systems, we derived analytic estimates for the importances of various dynamical processes that occured in our simulations, such as global tidal stripping, clump-clump heating, clump-clump merging and dynamical friction, including explicit dependence on the 4 independent free parameters. These analytic estimates were used to gain physical insight, and a comparison between our numerical results and our analytical estimate showed general agreement.

With an increase in the mass fraction of clumps and a decrease in the number of clumps, the velocity dispersion ratio of a clump to the halo increases, and all the dynamical processes become more efficient in a general sense, provided this velocity dispersion ratio is still much less than one. As this velocity dispersion ratio increases close to unity, the heating effect by close encounters between clumps becomes less important, while other dynamical processes continue to become more important.

We find that due to the angular momentum transport seen, spiral galaxies cannot be formed in a system with a large clump mass fraction and a small number of clumps. Our simulations show that the angular momentum of the clumpy component does not change significantly in a Hubble time for a small mass fraction in the clumps and a large number of clumps. This result confirms the standard formation scenario for disk galaxies, in that the merging history of spirals is restricted to the accretion of many, small systems, still conserving their specific angular momentum. This is similar in some aspects to the model of Vitvitska et al. (2001), who find that one may also end up with ellipticals with rapidly spinning halos, especially if only a small number of massive clumps are involved in the final merging events.

The final density profile, starting with clumps embedded in a smooth Plummer component, shows a rather steep inner part, depending on the clump mass fraction. Thus our experiments support the explanation for the origin of the central cusps in the density profiles of halos in cosmological numerical simulation as the consequence of self-limiting merging of small dense halos, as suggested by Syer and White (1998). The surface density profile of the debris formed from the disrupted clumps can be best fit by the de Vaucouleurs law for all our simulations. Thus our simulations do not produce cD galaxies, a deficiency which suggests that the formation of cD galaxies reauires conditions outside the parameter ranges covered by our simulations. Plausibly the formation of the outer envelope characteristic of cDs requires some later infall of large galaxies into the central regions of the cluster, possibly enhanced by dynamical friction.

The iso-density contours of both the dark halo and debris are consistently triaxial after the collapse, even though the initial configurations are oblate. The triaxiality of the final system is found to depend on the clumpiness of the initial system. A high level of clumpiness, such as having only a few, but massive, clumps, helps to produce a round system. This further suggests that giant ellipticals, which tend to be rounder than small ellipticals, may be formed in a more clumpy environment than small ellipticals.

We introduced a new measure of ‘streams’ in the final system, and found that persistent signatures of disrupted clumps in the outer halo can, at an overdensity of $`>0.75`$, cover 10% to 20% of the sky. The properties of these streams depend on the initial conditions in an intuitive way, in that the number of streams correlates with the number of initial clumps, and their orientation is more radial in collapse simulations.

Some extensions of our work that would further our understanding of the role of substructure in the formation of galaxies and of cluster of galaxies include N-body simulations with (1) a mass spectrum of clumps rather than identical clumps; (2) less restricted coverage of parameter space, in particular to larger collapse factors and (3) explicit incorporation of a merging criterion into the calculation at each time step.

We acknowledge support from NSF grant AST-9804706 (RFGW) and the PPARC visitor programme. RFGW thanks all in the Astronomy group at St Andrews for providing a pleasant and stimulating environment.

## Appendix A Parameters and Units of the Simulations

The density profile, mass profile, corresponding potential, and the total potential energy of a Hernquist sphere are:

$`\rho (r)`$ $`=`$ $`{\displaystyle \frac{M}{2\pi }}{\displaystyle \frac{a}{r(a+r)^3}},`$ (A1)
$`M(r)`$ $`=`$ $`{\displaystyle \frac{Mr^2}{(r+a)^2}},`$ (A2)
$`\mathrm{\Psi }(r)`$ $`=`$ $`{\displaystyle \frac{GM}{a+r}},`$ (A3)
$`\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{GM^2}{6a}},`$ (A4)
where $`a`$ is the core radius.

The density profile, mass profile, corresponding potential and the total potential energy of a Plummer sphere are:

$`\rho (r)`$ $`=`$ $`{\displaystyle \frac{3M}{4\pi a^3}}\left(1+{\displaystyle \frac{r^2}{a^2}}\right)^{5/2},`$ (A5)
$`M(r)`$ $`=`$ $`{\displaystyle \frac{Mr^3}{(r^2+a^2)^{3/2}}},`$ (A6)
$`\mathrm{\Psi }(r)`$ $`=`$ $`{\displaystyle \frac{GM}{\sqrt{r^2+a^2}}},`$ (A7)
$`\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{3\pi }{32}}{\displaystyle \frac{GM^2}{a}},`$ (A8)
where $`a`$ is the core radius.

The half mass radius for the Hernquist profile is $`r_{\frac{1}{2}}=\left(1+\sqrt{2}\right)a,`$ while for the Plummer profile, $`r_{\frac{1}{2}}1.3a`$. The mean crossing time of the dark halo, given by the ratio of the half mass radius, $`r_{\frac{1}{2}h}`$, to the circular velocity at the half mass radius, $`v_{\frac{1}{2}h}`$, is, for the Hernquist dark halo, $`t_{\frac{1}{2}hc}=r_{\frac{1}{2}h}/v_{\frac{1}{2}h}=\sqrt{2}\left(1+\sqrt{2}\right)^{1.5}=5.3`$, while for the Plummer dark halo $`t_{\frac{1}{2}hc}=2.1`$. The Hernquist models are allowed to run for about $`t_{sim}126`$, i.e. approximately 23 times the mean halo crossing time while for Plummer models, the simulation is stopped after $`t_{sim}56`$, i.e. approximately 27 times the mean halo crossing time. For a typical dark halo of a galaxy, the mean crossing time is about 0.5 Gyr, for $`v_{\frac{1}{2}h}200kms^1`$ and $`r_{\frac{1}{2}h}100kpc`$. Thus a typical simulation duration of 25 mean crossing times is approximately a Hubble time. The timescale for a dark halo typical of a cluster of galaxies is approximately the same.

There are some simple and convenient relations among the mean quantities of a given mass profile. For many reasonable spherical mass profiles, including both the Plummer and Hernquist spheres, there is an empirical relation between the half mass radius $`r_{\frac{1}{2}}`$ and the virial radius $`r_{vi}`$, $`r_{\frac{1}{2}}0.8r_{vi}`$, where the virial radius is defined by the potential energy $`|\mathrm{\Phi }|\frac{GM^2}{2r_{vi}}`$ (Spitzer 1987). In the case of an isotropic velocity dispersion tensor, we can use the mean one-dimensional velocity dispersion $`\sigma `$, defined by $`|\frac{\mathrm{\Phi }}{2}|\frac{3}{2}M\sigma ^2`$, instead of the circular velocity at the half mass radius radius $`v_{\frac{1}{2}}`$. These are related by $`v_{\frac{1}{2}}=\sqrt{\frac{6}{1.6}}\sigma 2\sigma `$.

## Appendix B Numerical Effects

We first discuss the effects of two-body relaxation, for which the two bodies can be two particles both from a single clump, or two particles with one from a clump and the other from the dark halo, or two particles both from the dark halo. Then we discuss the heating effect of a clump by dark halo particles.

### B.1 Two-body Relaxation

First we just consider the case that the clump is isolated, i.e. the two bodies under consideration are both from a single clump. We know the evaporation timescale is usually 100 times the relaxation timescale (Binney & Tremaine 1987) for a given Maxwellian system; this sets the timescale over which a clump could suffer evaporative mass loss if the two-body relaxation time is too short compared to the simulation time. The two-body relaxation time in a single clump is (Sellwood 1987)

$$t_{relax}\frac{0.14N_{cl}}{\mathrm{ln}\left(\frac{r_c}{3r_s}\right)}t_{\frac{1}{2}cc,}$$

where the $`N_{cl}`$ is the number of particles in a single clump, $`r_s`$ is the softening-length and $`t_{\frac{1}{2}cc}`$ is the internal crossing time within a clump. The softening-length in the treecode is set at $`r_s=0.5r_{\frac{1}{2}cl}N_{cl}^{1/3}`$; this relatively large value is chosen to suppress the unphysical two-body relaxation (White 1978), which otherwise could occur due to the small number of particles used in each clump (ranging from 150 to 2400).

In terms of the halo crossing time, this two-body relaxation time can be written as

$$t_{relax}\frac{0.14N_{cl}}{\mathrm{ln}\left(\frac{r_c}{3r_s}\right)}\frac{t_{\frac{1}{2}hc}}{\sqrt{\rho _{cl}/\rho _0}}.$$

Thus the two-body relaxation time in a single clump is shorter, in terms of the halo crossing time, for larger values of the density contrast $`\rho _{cl}/\rho _0`$. This poses a limitation for us in exploring the parameters for clumps of large density contrast, in addition to the limitation from the CPU time used, which is determined by the total number of particles. For our computed models with $`N_{cl}=300`$ and $`\rho _{cl}/\rho _0=64`$, we have $`t_{relax}3.5t_{\frac{1}{2}hc}`$. In the worst case, with $`N_{cl}=150`$ and $`\rho _{cl}/\rho _0=64`$, we have $`t_{relax}2t_{\frac{1}{2}hc}`$. Since our simulation time is about $`25t_{\frac{1}{2}hc}`$, generally our clumps should suffer somewhat from two-body relaxation effects, which can lead to some internal manifestation of collisional behaviour such as mass segregation (to which we are insensitive since we have adopted equal mass particles). However, we do not expect there to be significant mass loss due to evaporation effects, not only because the estimated evaporation time scale is very long compared to our simulation time, but also due to the fact that adoption of a softening-length can inhibit evaporation. Furthermore, we are not interested, in this study, in the evolution of the internal structure of the clump, and thus our chosen parameter values are acceptable.

Next, we turn to consideration of two-body relaxation in the dark halo, i.e. the two bodies participating in the two-body relaxation are both from the dark halo. Similarly to the expressions above, we have

$$t_{relax,h}\frac{0.14N_h}{\mathrm{ln}\left(\frac{r_h}{3r_s}\right)}t_{\frac{1}{2}hc}$$

where $`N_h`$ is the number of particles in the dark halo. Typically $`t_{relax,h}700t_{\frac{1}{2}hc}`$. Thus the two-body relaxation effect in the halo can be ignored.

In addition to the possible evaporation of clump particles due to two-body heating between two particles within the same clump, there is another heating process we should consider, caused by two-body heating between a halo particle and a clump particle (Carlberg 1994; van Kampen 1995). The disruption time due to this heating process is

$$t_{hpcp}\frac{N_{cl}}{4\mathrm{ln}N_h}\frac{r_{\frac{1}{2}h}}{r_{\frac{1}{2}c}}t_{\frac{1}{2}hc}.$$

For our simulations, this disruption time is much longer than the total integration time; for a typical model, $`t_{hpcp}130t_{\frac{1}{2}hc}`$ and we can safely ignore this heating effect too.

### B.2 Particle-Clump Heating Effect

Numerical limitations resulting from a small number of particles can lead not only to unwelcome collisional effects as discussed above, but also can lead to spurious heating effects within clumps. While in our simulations we have both halo particles and clump particles of the same mass, we shall consider the more general case of unequal particles, and indeed of a heavier halo particle. As a (heavy) dark halo particle passes by or through a clump, it can tidally heat the clump, and the clump can eventually be disrupted (Carlberg 1994). One can roughly estimate the magnitude of this effect by using the impulse approximation, provided the duration of the encounter is short compared to the internal crossing time of the clump (Binney & Tremaine 1987). The binding energy of the clump is given by $`E_b=|\mathrm{\Phi }/2|=\frac{Gm_c^2}{4r_{cvi}}`$. For a given clump, the heating rate is

$$\dot{E}=\frac{4\sqrt{\pi }G^2m_p^2n_pm_c\overline{r_c^2}}{3\sigma _hb_{min}^2},$$
(B1)
where $`m_p`$ is the mass of one dark halo particle, $`n_p`$ is the mean number density of the dark halo particles within the half mass radius of the dark halo, $`m_c`$ is the mass of the clump, $`\overline{r_c^2}`$ is the mean-square radius of the clump, $`\sigma _h`$ is the one-dimensional velocity dispersion of the dark halo and $`b_{min}`$ is the minimum impact parameter. The size of the single dark halo particle can be set to the value of the softening length $`r_s`$. Since $`r_sr_{\frac{1}{2}c}`$, we can choose $`b_{min}r_{\frac{1}{2}c}`$ to include head-on encounters, thus modifying the above equation by multiplying a factor of $`g3`$ (Binney & Tremaine 1987). The appropriate value of $`\overline{r_c^2}`$ is somewhat uncertain, since, for example, for the Plummer profile $`\overline{r_c^2}`$ is divergent. The heating rate however remains finite. This result can be understood in the following way: the heating by such a close encounter causes the particles in the clump to gain a velocity increment proportional to their projected distance from the centre of the clump, and thus while the inner particles are essentially unaffected, the outermost particles, though few in number, can gain most of the kinetic energy, resulting in their removal from the clump. Thus in the calculation of $`\overline{r_c^2}`$, we should ignore those particles lying too far from the clump centre. Here we adopt $`\overline{r_c^2}`$ to be $`2r_{\frac{1}{2}c}^2`$. After some algebraic manipulations using the relations in section 2.2 and Appendix A, we have that the clump disruption timescale, due to halo particle-clump heating, is

$$t_{pc}=0.03N_{tot}f^{2/3}n_c^{2/3}\left(\rho _{cl}/\rho _0\right)^{1/3}t_{\frac{1}{2}hc},$$
(B2)
where $`N_{tot}`$ is the total number of particles used in the simulation. For the worst case (model S) with $`f=0.1`$, $`\rho _{cl}/\rho _0=2.7`$, $`n_c=20`$ and $`N_{tot}60000`$, we have $`t_{pc}73t_{\frac{1}{2}hc}`$. So again this particle-clump heating effect can be ignored in our simulations. As discussed in the text (section 3.1.3), the physical clump-clump heating effect is much more important.

## Appendix C Dynamical Friction

We estimate the expected amplitude of dynamical friction under the conditions we simulated by adopting, for Hernquist profile halos, the phase space distribution function he derived, which neglects angular momentum. This is

$$f(E)=\frac{(1f)M_{tot}}{8\sqrt{2}\pi ^3a^3v_g^3m_p}K_0(q),$$
(C1)
where

$$K_0(q)=\frac{1}{\left(1q^2\right)^{5/2}}\times \left[3\mathrm{sin}^1q+q\left(1q^2\right)^{1/2}\left(12q^2\right)\left(8q^48q^23\right)\right],$$
(C2)
where $`a`$ is the core radius of the halo, $`v_g=\sqrt{GM_{tot}/a}`$, $`q=\sqrt{E}/v_g`$, and $`E=\frac{1}{2}v^2+\mathrm{\Psi }(r)`$ is the specific energy of a halo particle. Define $`xr/a`$. Then we have

$$v^2=\frac{2}{x+1}v_g^22v_g^2q^2.$$
(C3)
Since $`f(E)`$ is non-negative for all $`v<v_e`$, where $`v_e=\sqrt{2\mathrm{\Psi }(r)}`$ is the escape velocity, and if the clump is moving with velocity $`v_M`$ less than the escape velocity at radius r, we can simply substitute for $`f(E)`$ in equation (10). With $`M_cm_p`$ (symbols defined in section 3.1.4), we have

$$_0^{v_M}f(v)v^2𝑑v=\frac{(1f)M_{tot}}{4\pi ^3a^3}f_0(x,v_M)$$
(C4)
with

$$f_0(x,v_M)=_{q_{v=v_M}}^{q_{v=0}}\left(\frac{1}{1+x}q^2\right)^{1/2}K_0(q)q𝑑q$$
(C5)
where $`q_{v=0}=\left(1+x\right)^{1/2}`$ and $`q_{v=v_M}=\sqrt{\frac{1}{1+x}\frac{v_M^2}{2v_g^2}}`$. Thus the characteristic dynamical friction timescale becomes

$$t_{df}^1=\frac{4\mathrm{ln}\mathrm{\Lambda }f(1f)}{\pi n_c}\frac{v_g^3}{v_M^3}\frac{v_g}{a}f_0(x,v_M)$$
(C6)
To obtain the exact dynamical friction time from the above equations, for given initial conditions, one needs to integrate over time for the time-dependent orbit (e.g. Cora, Muzzio & Vergne 1997). While it is generally assumed that orbits are quickly circularized by dynamical friction, and thus that the Chandrasekhar formula for circular orbits may be applied, this is not always a good assumpotion (van den Bosch et al. 1999). However, there are some empirical relations derived from numerical studies (see a short summary in Colpi et al. 1999) that relate the dynamical friction time calculated numerically for a particle on a given initial orbit, i.e. $`E`$ and $`J`$, to the prediction using the circular orbit approximation, giving $`t_{df}=t_cϵ^\alpha `$, where $`ϵ=J(E)/J_c(E)`$, and $`J_c(E)`$ is the angular momentum of a particle on a circular orbit with the same binding energy. The parameter $`ϵ`$ is simply dependent on the orbit eccentricity only, while $`\alpha `$ can range from 0.4 to 0.8. For a typical orbital eccentricity of $`0.6`$, the circular orbit approximation provides an overestimate by a factor of around 2. We then assume that a clump is initially moving on a circular orbit with $`v_M=v_c=v_g\sqrt{x}/(1+x)`$ (Tsuchiya & Shimada 2000), and modify the estimates of the dynamical friction timescales by this factor of two.

Then from equation (9) in the main text, the specific angular momentum loss rate should be equal to the torque caused by the dynamical friction at radius $`r`$, and we have

$$\frac{d\left(rv_c\right)}{dt}=\frac{v_c}{t_{df}}r.$$
(C7)
It can be further written as

$$\frac{dr}{dt}=\frac{2r}{t_{df}}\frac{1}{\frac{\mathrm{ln}M}{\mathrm{ln}r}+1}.$$
(C8)
With some algebraic manipulation, the dynamical friction time can be expressed as

$$t_{df}=\frac{\pi n_c}{4\mathrm{ln}\mathrm{\Lambda }f(1f)}\frac{a}{v_g}B(x_0)$$
(C9)
where

$$B(x_0)=_0^{x_0}\frac{(3+x)x^{1/2}}{2(1+x)^4}\frac{1}{f_0(x)}𝑑x,$$
(C10)
and $`x_0`$ is the initial radius. The appropriate value of $`\mathrm{\Lambda }`$ can be calculated as (White 1976)

$$\mathrm{ln}\mathrm{\Lambda }=\frac{1}{M_c^2}_0^{b_{max}}D^3𝑑D\left[_D^{\mathrm{}}\frac{M_c(r)dr}{r^2\left(r^2D^2\right)^{1/2}}\right]^2,$$
(C11)
where $`M_c(r)`$ is the mass profile of a clump, here given by the Plummer law. The maximum impact parameter is chosen to be the half mass radius of the dark halo, $`b_{max}=r_{\frac{1}{2}h}`$. Thus

$$\mathrm{ln}\mathrm{\Lambda }=\frac{1}{2}\mathrm{ln}\left(\frac{r_{\frac{1}{2}h}^2}{a_c^2}+1\right)+\frac{1}{2}\left(\frac{1}{\frac{r_{\frac{1}{2}h}^2}{a_c^2}+1}1\right)$$
(C12)
where $`a_c`$ is the core radius of the clump Plummer profile. For our computed models, typically $`\mathrm{ln}\mathrm{\Lambda }3`$.

The above expressions are used to derive the expressions in section 3.1.4 of the text.