# Substructure in Dark Halos: Orbital Eccentricities and Dynamical Friction

## 1 Introduction

In the standard cosmological model galaxies and clusters form by hierarchical clustering and merging of small density perturbations that grow by gravitational instability. In this standard picture the mass of the Universe is dominated by dissipationless dark matter which collapses to form dark halos, inside of which the luminous galaxies form. It was once assumed that the previous generation of substructure was erased at each level of the hierarchy (White & Rees 1978). However, high resolution $`N`$-body simulations have recently shown that some substructure is preserved at all levels (Moore, Katz & Lake 1996; Klypin et al. 1997; Brainerd, Goldberg & Villumsen 1998; Moore et al. 1998; Ghigna et al. 1998). This is consistent with observations which reveal substructure in a variety of systems: globular clusters within galaxies, distant satellites and globulars in the halos of galaxies, and galaxies within clusters. This substructure evolves as it is subjected to the forces that try to dissolve it: dynamical friction, tides from the central objects and impulsive collisions with other substructure. The timescales for survival and the nature of the dissolution guide our understanding of the formation processes, most of which depend on the nature of the orbits. Tides strip satellites on elongated orbits. An elongated orbit and a circular one clearly evolve differently owing to dynamical friction, especially if there is a disk involved. The disk heating similarly depends on the orbits of the satellites. The nature of the orbits also effects the nature and persistence of tidal streams at breakup, and the mutual collisions of structures as considered by galaxy harassment.

Since the clustering and merging of halo substructures is one of the cornerstones of the hierarchical structure formation scenario, a comprehensive understanding of their orbital properties is of invaluable importance when seeking to understand the formation and evolution of structure in the Universe. We have found that the properties of orbits within spherical, fully relaxed systems has received little attention and is often misrepresented. Hence, the first goal of this paper is to derive a statistical characterization of the orbits in potential/density distributions that describe galaxies within clusters, globular clusters within galaxies, etc. We find that orbits are far more elongated than typically characterized. Recently, Ghigna et al. (1998) used high resolution $`N`$-body simulations to investigate the orbital properties of halos within clusters. They found that orbits of the subhalos are strongly elongated with a median apo-to-pericenter ratio of approximately $`6`$. We compare our results on equilibrium spherical potentials with the distribution of orbits in their cluster that was simulated in a cosmological context, and show that the orbital eccentricities of the subhalos are consistent with an isothermal halo that is close to isotropic.

In the second part of this paper, we use high resolution $`N`$-body simulations to calculate the dynamical friction on eccentric orbits. Past studies have compared numerical simulations with Chandrasekhar’s dynamical friction formula (e.g., White 1978, 1983; Tremaine 1976, 1981; Lin & Tremaine 1983; Tremaine & Weinberg 1984; Bontekoe & van Albada 1987; Zaritsky & White 1988; Hernquist & Weinberg 1989; Cora, Muzzio & Vergne 1997). Comprehensive overviews of these studies with discussions regarding the local versus global nature of dynamical friction can be found in Zaritsky & White (1988) and Cora et al. (1997). Most of the studies followed the decay of circular or only slightly eccentric orbits. The two exceptions are Bontekoe & van Albada (1987) and Cora et al. (1997). The former examined the orbital decay of a ‘grazing encounter’ of a satellite on an elliptical orbit that grazes a larger galaxy at its pericenter. In this case, dynamical friction occurs only near pericenter. The pericentric radius remains nearly fixed with significant circularization of the orbit in just a few dynamical times (see also Colpi (1998) for an analytical treatment based on linear response theory). Cora et al. followed satellites on eccentric orbits that were completely embedded in a dark halo, but they didn’t discuss the dependence of decay time on orbital eccentricity or the change in orbital eccentricity with time. We use fully self-consistent $`N`$-body simulations with $`50,000`$ halo particles to calculate dynamical friction on eccentric orbits. We show that the timescale for dynamical friction is shorter for more eccentric orbits, but that the dependence is considerably weaker than claimed previously by Lacey & Cole (1993) based on an analytical integration of Chandrasekhar’s formula. In addition, we show that, contrary to common belief, dynamical friction does not lead to circularization. All in all, dynamical friction only leads to a very moderate change in the distribution of orbital eccentricities over time.

In Section 2 we derive the distributions of orbital eccentricities for a number of spherical densities/potentials using both analytical and numerical methods. Section 3 describes our $`N`$-body simulations of dynamical friction on eccentric orbits in an isothermal halo. In Section 4 we discuss a number of applications. Our results and conclusions are presented in Section 5.

## 2 Orbital eccentricities

### 2.1 The singular isothermal sphere

The flat rotation curves observed for spiral galaxies suggest that their dark halos have density profiles that are not too different from isothermal. Hence, we start our investigation with the singular isothermal sphere, whose potential, $`\mathrm{\Phi }`$, and density, $`\rho `$, are given by

$$\mathrm{\Phi }(r)=V_c^2\mathrm{ln}(r),\rho (r)=\frac{V_c^2}{4\pi Gr^2}.$$
(1)   
Here $`V_c`$ is the circular velocity, which is constant with radius.

#### 2.1.1 Analytical method

For non-Keplerian potentials, in which the orbits are not simple ellipses, it is customary to define a generalized orbital eccentricity <sup>6</sup><sup>6</sup>6Throughout this paper we use $`e`$ to denote the orbital eccentricity, $`\mathrm{e}`$ refers to the base of the natural logarithm, and $`E`$ refers to the energy., $`e`$, as:

$$e=\frac{r_+r_{}}{r_++r_{}}.$$
(2)   
Here $`r_{}`$ and $`r_+`$ are the peri- and apocenter respectively. For an orbit with energy $`E`$ and angular momentum $`L`$ in a spherical potential, $`r_{}`$ and $`r_+`$ are the roots for $`r`$ of

$$\frac{1}{r^2}+\frac{2[\mathrm{\Phi }(r)E]}{L^2}=0,$$
(3)   
(Binney & Tremaine 1987). For each energy, the maximum angular momentum is, for a singular isothermal sphere, given by $`L_c(E)=r_c(E)V_c`$. Here $`r_c(E)`$ is the radius of the circular orbit with energy $`E`$ and is given by

$$r_c(E)=\mathrm{exp}\left[\frac{EV_c^2/2}{V_c^2}\right].$$
(4)   
Upon writing $`L=\eta L_c(E)`$ ($`0\eta 1`$)<sup>7</sup><sup>7</sup>7The quantity $`\eta `$ is generally called the orbital circularity , one can rewrite equation (3) for a singular isothermal sphere such that the apo- and pericenter are given by the roots for $`x=r/r_c`$ of

$$\frac{1}{x^2}+\frac{2}{\eta ^2}\mathrm{ln}(x)\frac{1}{\eta ^2}=0.$$
(5)   
As might be expected in this scale free case, the ratio $`r_+/r_{}`$ depends only on the orbital circularity $`\eta `$ and is independent of energy. This dependence is shown in Figure 1.

At a certain radius, the average of any quantity $`S`$ is determined by weighting it by the distribution function (hereafter DF) and integrating over phase space. For the singular isothermal sphere this yields

$$\overline{S}(r)=\frac{4\pi }{r^2\rho (r)}\underset{\mathrm{\Phi }(r)}{\overset{\mathrm{}}{}}𝑑E\underset{0}{\overset{r\sqrt{2(E\mathrm{\Phi })}}{}}f(E,L)S(E,L)\frac{LdL}{\sqrt{2(E\mathrm{\Phi })L^2/r^2}}.$$
(6)   
In what follows we consider the family of quasi-separable DFs

$$f(E,L)=g(E)h_a(\eta ).$$
(7)   
This approach makes the solution of equation (6) analytically tractable. The general properties of spherical galaxies with this family of DFs are discussed in detail by Gerhard (1991, hereafter G91). We adopt a simple parameterization for the function $`h_a(\eta )`$:

$$h_a(\eta )=\{\begin{array}{ccc}\mathrm{tanh}\left(\frac{\eta }{a}\right)/\mathrm{tanh}\left(\frac{1}{a}\right)\hfill & a>0\hfill & \\ 1\hfill & a=0\hfill & \\ \mathrm{tanh}\left(\frac{1\eta }{a}\right)/\mathrm{tanh}\left(\frac{1}{a}\right)\hfill & a<0\hfill & \end{array}$$
(8)   
For $`a=0`$, the DF is isotropic. Radially anisotropic models have $`a<0`$, whereas positive $`a`$ correspond to tangential anisotropy. For a quasi-separable DF of the form (7), and for the eccentricity $`e`$ which depends on $`\eta `$ only, equation (6) yields

$$\overline{e}(r)=\frac{4\pi }{r\rho (r)}\underset{\mathrm{\Phi }(r)}{\overset{\mathrm{}}{}}𝑑Eg(E)L_c(E)\underset{0}{\overset{\eta _{\mathrm{max}}}{}}h_a(\eta )e(\eta )\frac{\eta d\eta }{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}},$$
(9)   
where

$$\eta _{\mathrm{max}}=\frac{r\sqrt{2(E\mathrm{\Phi })}}{L_c(E)}.$$
(10)   
For a singular isothermal sphere, G91 has shown that the energy dependence of the DF is given by

$$g(E)=\frac{\mathrm{e}}{16\pi ^2GV_cu_H}\mathrm{exp}\left[\frac{2E}{V_c^2}\right],$$
(11)   
where

$$u_H=\underset{0}{\overset{\mathrm{}}{}}𝑑u\mathrm{e}^u\underset{0}{\overset{\eta _{\mathrm{max}}}{}}h_a(\eta )\frac{\eta d\eta }{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}}.$$
(12)   
Here $`\eta _{\mathrm{max}}`$ depends on $`u`$ only and is given by

$$\eta _{\mathrm{max}}=\sqrt{2\mathrm{e}}\sqrt{u}\mathrm{e}^u.$$
(13)   
Substitution of (11) and (12) in (9) yields (upon substituting $`u=(E\mathrm{\Phi })/V_c^2`$)

$$\overline{e}(r)=\overline{e}=\frac{1}{u_H}\underset{0}{\overset{\mathrm{}}{}}𝑑u\mathrm{e}^u\underset{0}{\overset{\eta _{\mathrm{max}}}{}}h_a(\eta )e(\eta )\frac{\eta d\eta }{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}}.$$
(14)   
Note that, due to the scale-free nature of the problem, this average is independent of radius.

We have numerically solved this integral as a function of the anisotropy parameter $`a`$. The results are shown in Figure 2. As expected, the average eccentricity decreases going from radial to tangential anisotropy. The average orbital eccentricity of an isotropic, singular isothermal sphere is $`\overline{e}=0.55`$. An orbit with this eccentricity has $`r_+/r_{}=3.5`$.

#### 2.1.2 Numerical method

To examine the actual distribution of eccentricities, rather than just calculate moments of the distribution, we Monte Carlo sample the quasi-separable DF of equation (7) for orbits in a singular isothermal potential and calculate their eccentricities. We provide a detailed description of this method in the Appendix.

The normalized distribution functions of orbital eccentricities for three different values of the anisotropy parameter $`a`$ are shown in Figure 3 (upper panels). The lower panels show the corresponding distributions of apo-to-pericenter ratios. Each distribution is computed from a Monte Carlo simulation with $`10^6`$ orbits. The thin vertical lines in Figure 3 show the 20th (dotted lines), 50th (solid lines), and 80th (dashed lines) percentile points of the distributions. The average eccentricity for the isotropic case ($`a=0`$) computed from the Monte Carlo simulations is $`\overline{e}0.55`$, in excellent agreement with the value determined in Section 2.1.1. About 15 percent of the orbits in the isotropic singular isothermal sphere have apo- to pericenter ratios larger than $`10`$, whereas only $`20`$ percent of the orbits have $`r_+/r_{}<2`$. Note however that these numbers depend strongly on the velocity anisotropy.

### 2.2 Tracer populations

In the previous two sections, we concentrated on the self-consistent case of a singular isothermal halo and the corresponding density distribution that follows from the Poisson equation. Tracer populations, however, do not necessarily follow the self-consistent density distribution. Consider a tracer population in a singular isothermal sphere potential with a density distribution given by

$$\rho _{\mathrm{trace}}(r)=\rho _0\left(\frac{r}{r_0}\right)^\alpha ,$$
(15)   
(the self-consistent case corresponds to $`\alpha =2.0`$). If we consider the same quasi-separable DF as in Section 2.1 (i.e., equation ), the energy dependence of the DF becomes

$$g(E)=\frac{\sqrt{\mathrm{e}}\rho _0r_0^\alpha }{4\pi V_c^3u_{H,\alpha }}\mathrm{exp}\left[\frac{\alpha E}{V_c^2}\right],$$
(16)   
with

$$u_{H,\alpha }=\underset{0}{\overset{\mathrm{}}{}}𝑑u\mathrm{e}^{(1\alpha )u}\underset{0}{\overset{\eta _{\mathrm{max}}}{}}h_a(\eta )\frac{\eta d\eta }{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}}.$$
(17)   
The average eccentricity can be computed by substituting equations (16) and (17) in (9), thereby using the expression of $`\mathrm{\Phi }(r)`$ given by equation (1). In Figure 4 we plot the average eccentricity thus computed as a function of the power-law slope $`\alpha `$. As can be seen, the orbital eccentricities depend only mildly on the slope of the density distribution. Note that for $`\alpha >3`$, the mass within any radius $`r>0`$ is infinite, whereas in the case $`\alpha <3`$, the mass outside any such radius is infinite. These differing properties in the two regimes are apparent in the behavior of the median eccentricity seen in Figure 4, with the intermediate case of $`\alpha =3`$ yielding a minimum. Due to the infinities involved the cases examined may not accurately represent true dark halos. After giving an example of how such infinities cause problems, the subsequent sections we examine other, more realistic, spherical potentials with finite masses.

### 2.3 Stellar hydrodynamics and the virial theorem in the isothermal potential

The difference between infinite and finite samples is evident when comparing the equation of stellar hydrodynamics to the virial theorem. For a sphere with an isotropic DF, the former takes the simple form for the one dimensional velocity dispersion $`\sigma `$ :

$$\frac{d}{dr}\left(\rho \sigma ^2\right)=\rho \frac{d\mathrm{\Phi }}{dr}$$
(18)   
For the tracer population of equation (15) and an isothermal potential, this simplifies to:

$$\sigma ^2=V_c^2\left(\frac{d\mathrm{ln}\rho }{d\mathrm{ln}r}\right)^1=\frac{V_c^2}{\alpha }$$
(19)   
In contrast, the virial theorem states that twice the kinetic energy is equal to the virial. Adopting particle masses of unity yields:

$$v^2=\stackrel{}{F}\stackrel{}{r}=\frac{V_c^2}{r}r$$
(20)   
Since, the expectation value of $`v^2`$ is $`3\sigma ^2`$, this reduces to:

$$\sigma ^2=\frac{V_c^2}{3}$$
(21)   
This difference owes to the assumption of finite versus infinite tracers. The answers match only for $`\alpha =3`$ where the divergence in mass is only logarithmic at both $`r0`$ and $`r\mathrm{}`$. Even the “self-consistent” case of $`\alpha =2`$ has a problem that is pointed out in problem \[4-9\] of Binney & Tremaine (1987). The kinetic energy per particle is $`V_c^2`$ in a model with only circular orbits and $`3V_c^2/2`$ if they are isotropic. Yet, they have the same potential and must satisfy the virial theorem. In Section 4.2, we return to this issue as we examine a case where the equations of stellar hydrodynamics have been used on astrophysical objects where a subsample of a finite number of global tracers has been observed. We now turn to the calculation of the distribution of eccentricities for finite sets of tracers.

### 2.4 Truncated isothermal sphere with core

Given the infinite mass and central singularity of the singular isothermal sphere, dark halos are often modeled as truncated, non-singular, isothermal spheres with a core:

$$\rho (r)=\frac{M}{2\pi ^{3/2}\gamma r_t}\frac{\mathrm{exp}(r^2/r_t^2)}{r^2+r_c^2}.$$
(22)   
Here $`M`$, $`r_t`$, and $`r_c`$ are the mass, truncation radius and core radius, respectively, and $`\gamma `$ is a normalization constant given by

$$\gamma =1\sqrt{\pi }\left(\frac{r_c}{r_t}\right)\mathrm{exp}(r_c^2/r_t^2)\left[1\mathrm{erf}(r_c/r_t)\right].$$
(23)   
Since for this density distribution the DF is not known analytically (i.e., this requires the knowledge of $`\rho (\mathrm{\Phi })`$ in order to solve the Eddington equation), and since this density distribution is no longer scale-free, we have to use a different approach in order to determine the distribution of orbital eccentricities. We employ the method described by Hernquist (1993). We randomly draw positions according to the density distribution. The radial velocity dispersion is computed from the second-order Jeans equations, assuming isotropy, i.e. assuming $`f=f(E)`$:

$$\overline{v_r^2}(r)=\frac{1}{\rho (r)}_r^{\mathrm{}}\rho (r)\frac{d\mathrm{\Phi }}{dr}𝑑r.$$
(24)   
We compute local velocities $`v`$ by drawing randomly a unit vector and then a magnitude from a Gaussian whose second moment is equal to $`\overline{v_r^2}`$, truncated at the local escape speed. Once the six phase-space coordinates are known, the energy and angular momentum are calculated, providing the apo- and pericenter of the orbit by solving equation (3). The resulting distributions of orbital eccentricities are not rigorous, since the velocity field has not been obtained from a stationary DF, but rather only used the second moments. However, as demonstrated by Hernquist (1993), $`N`$-body simulations run from these initial conditions are nearly in equilibrium (see also Section 3.2), which suggests that the eccentricity distributions derived are sufficiently close to the actual equilibrium distributions.

In Figure 5, we plot the 20th (dotted lines), 50th (solid lines), and 80th (dashed lines) percentile points of the distributions of eccentricity (left panel) and apo- to pericenter ratios (right panel), as functions of $`r_t/r_c`$. For $`r_c=r_t`$ the distribution of eccentricities is almost symmetric, with the median equal to $`0.50`$. When $`r_t/r_c`$ increases, the distribution becomes more and more skewed towards high eccentricity orbits; the distribution closely approaches that of the singular isothermal sphere in the limit $`r_t/r_c\mathrm{}`$.

### 2.5 Steeper halo profiles

To examine the dependence of the orbital eccentricities on the actual density distribution of the halos we determine the eccentricity distributions of two well-known models with steeper outer density profiles ($`\rho r^4`$):

$$\rho _J(r)=\frac{M}{4\pi }\frac{a}{r^2(r+a)^2},$$
(25)   
and

$$\rho _H(r)=\frac{M}{2\pi }\frac{a}{r(r+a)^3},$$
(26)   
where $`M`$ is the total mass. These profiles differ only in the steepness of the central cusp: the former one, known as the Jaffe (1983) profile, has a $`r^2`$-cusp, whereas the latter, known as the Hernquist (1990) profile, has a shallower $`r^1`$-cusp. We use the technique described in Section 2.4 to compute the distributions of orbital eccentricities for isotropic DFs $`f(E)`$.

The results are shown in Figure 6, where we plot the normalized eccentricity distributions for the Jaffe and Hernquist spheres. For comparison the results for the singular isothermal sphere with isotropic DF are reproduced as well. The three distributions are remarkably similar (the differences are best appreciated by comparing the thin lines indicating the percentile points). The distributions are progressively skewed toward higher eccentricities in the sequence isothermal, $`\rho _H(r)`$, $`\rho _J(r)`$, but only moderately so.

Navarro, Frenk & White (1995, 1996, 1997) have used cosmological simulations to argue that the outer density profiles of dark halos decline as $`r^3`$. Their profiles are likely created by a variety of numerical artifacts (Moore et al. 1998). However, our results suggest that such debates will not significantly alter the expected distributions of eccentricities; velocity anisotropy is far more important than the details of the density profile.

## 3 Orbital decay of eccentric orbits by dynamical friction

The orbits of substructures within halos change owing to dynamical friction. Chandrasekhar (1943) derived the local deceleration of a body with mass $`M`$ moving through an infinite and homogeneous medium of particles with mass $`m`$. The deceleration is proportional to the mass $`M`$, such that more massive sub-halos sink more rapidly. As long as $`Mm`$ the frictional drag is proportional to the mass density of the medium, but independent of the mass $`m`$ of the constituents. For an isotropic, singular isothermal sphere the deceleration is

$$\frac{d\stackrel{}{v}}{dt}=\frac{GM_s}{r^2}\mathrm{ln}\mathrm{\Lambda }\left(\frac{v}{V_c}\right)^2\left\{\mathrm{erf}\left(\frac{v}{V_c}\right)\frac{2}{\sqrt{\pi }}\left(\frac{v}{V_c}\right)\mathrm{exp}\left[\left(\frac{v}{V_c}\right)^2\right]\right\}\stackrel{}{e}_v$$
(27)   
with $`M_s`$ and $`v`$ the mass and velocity of the object being decelerated, $`r`$ the distance of that object from the center of the halo, $`\mathrm{ln}\mathrm{\Lambda }`$ the Coulomb logarithm, and $`\stackrel{}{e}_v`$ the unit velocity vector (Tremaine 1976; White 1976a).

As mentioned earlier, remarkably little attention has been paid to the effects of dynamical friction on eccentric orbits—our goal for the rest of this section. Our main objective is to use both analytical and numerical tools to investigate the change of orbital eccentricity with time, and the dependence of the dynamical friction timescale on the (intrinsic) eccentricity. Unlike most previous studies, we will not focus on testing Chandrasekhar’s dynamical friction formula, or on studying the exact cause of the deceleration (i.e., local or global), as this has been the topic of discussion in many previous papers.

### 3.1 The time-dependence of orbital eccentricity

We investigate the rate at which orbital eccentricity changes due to dynamical friction. For simplicity, we focus on the evolution of orbital eccentricity in the singular isothermal sphere, for which

$$\frac{de}{dt}=\frac{d\eta }{dt}\frac{de}{d\eta },$$
(28)   
with $`\eta =L/L_c(E)`$ (see Section 2.1). Using equation (4), we find

$$\frac{d\eta }{dt}=\eta \left\{\frac{1}{L}\frac{dL}{dt}\frac{1}{V_c^2}\frac{dE}{dt}\right\}.$$
(29)   
Because of dynamical friction the energy and angular momentum are no longer conserved, and

$$\frac{dE}{dt}=v\frac{dv}{dt}$$
(30)   
and

$$\frac{dL}{dt}=r\frac{dv_{}}{dt}.$$
(31)   
Since the frictional force acts in the direction opposite of the velocity

$$\frac{dv_{}}{dt}=\frac{v_{}}{v}\frac{dv}{dt}.$$
(32)   
and, upon combining equations (28) - (32), we find

$$\frac{de}{dt}=\frac{\eta }{v}\frac{de}{d\eta }\left[1\left(\frac{v}{V_c}\right)^2\right]\frac{dv}{dt}.$$
(33)   
Here $`dv/dt`$ is the frictional deceleration given by equation (27). Since both $`dv/dt<0`$ and $`de/d\eta <0`$ (see Figure 1), we can immediately derive the sign of $`de/dt`$ at apo- and pericenter. At apocenter $`v<V_c`$, such that $`de/dt>0`$, whereas at pericenter $`v>V_c`$ and thus $`de/dt<0`$. This explains the circularization found for ‘grazing encounters’ (Bontekoe & van Albada 1987), as dynamical friction happens only near pericenter. Equation (27) shows that $`dv/dtr^2`$, and the change in eccentricity is thus larger at pericenter than at apocenter. However, the time spent near pericenter is shorter than near apocenter, so that the evolution of the eccentricity can not be determined by inspection. In the next sections, we use numerical simulations.

### 3.2 $`N`$-body simulations

We perform a set of fully self-consistent $`N`$-body simulations with a large number of particles in order to examine the effects of dynamical friction on eccentric orbits in a massive halo. The halo is modeled by a truncated isothermal sphere (equation ), with total mass of unity, a core-radius $`r_c=1`$, and a truncation radius $`r_t=50`$. Scaled to the Milky Way, we adopt a unit of mass of $`10^{12}\mathrm{M}_{}`$ and a unit of length of $`4`$ kpc. With the gravitational constant set to unity, the units of velocity and time are $`1037\mathrm{km}\mathrm{s}^1`$ and $`3.8`$ Myr, respectively.

The initial velocities of the halo particles are set up as described in Section 2.4, following the procedure of Hernquist (1993). Because of this particular method the halo is not necessarily in equilibrium, nor is it expected to be virialized. In order to remove effects of the halo’s virialization on the decay of the orbiting substructure, we first simulate the halo in isolation for 10 Gyrs. At the end, the halo has nicely settled in virial equilibrium. Figure 7 shows the initial density profile compared to that after 10 Gyr. As can be seen, and as already pointed out by Hernquist (1993), the density profile has not changed significantly after 10 Gyr.

We are interested in the effects of dynamical friction on galactic objects that range from $`10^6\mathrm{M}_{}`$ (globular clusters) to $`10^{10}\mathrm{M}_{}`$ (a massive satellite). To simulate dynamical friction on an object of mass $`M_s`$ orbiting in a halo of mass $`M_h`$, we require $`M_sm`$. In our simulations, we insist that $`N10M_h/M_s`$. A simulation of the orbital decay of a globular cluster in a galactic halo of $`10^{12}\mathrm{M}_{}`$ thus requires $`N10^7`$, clearly too large a number for practical purposes. We run our simulations with $`N=50,000`$ particles, roughly an order of magnitude increase over most previous work. The highest resolution, self-consistent simulation aimed at investigating dynamical friction has so far been performed by Hernquist & Weinberg (1989) using $`N=20,000`$ and focusing only on circular orbits. We discuss the influence of the number of halo particles in Section 3.3.4. With $`N=50,000`$, $`M_h=10^{12}\mathrm{M}_{}`$, and the requirement $`N10M_h/M_s`$, our minimum satellite mass is $`2\times 10^8\mathrm{M}_{}`$. A list of the parameters for each of our simulations is presented in Table 1. Models 1, 2, and 3 have an initial apocenter of 160 kpc, and start on an orbit with $`e=0.8`$. The initial pericenter for these orbits is located at $`17.9`$ kpc from the center of the halo, well outside the halo’s core radius ($`r_c=4`$ kpc). Satellites on orbits of lower eccentricity (Models 4, 5, and 6) are started from smaller apocentric radii, such that the initial specific energy of all satellites in all models is equal.

The satellite is initially positioned at $`(x,y,z)=(r_+,0,0)`$, with $`v_x=v_z=0`$ and $`v_y`$ chosen such as to obtain an initial orbit with the desired eccentricity: we randomly draw velocities $`0v_yv_{\mathrm{escape}}`$ and determine the orbital eccentricity of the satellite as described in Section 2.1. We repeat this procedure until we find an eccentricity within one percent of the desired value. At the start of each simulation, the satellite is introduced instantly in the virialized halo potential.

The simulations use PKDGRAV (Dikaiakos & Stadel 1996; Stadel & Quinn 1998), a stable and well-tested, both spatially and temporally adaptive tree code, optimized for massively parallel processors. It uses an open ended variable timestep criteria based upon the local acceleration (Quinn et al. 1997). Forces are computed using terms up to hexadecapole order and a tolerance parameter of $`\theta =0.8`$. The code uses spline kernel softening, for which the forces become completely Newtonian at 2 softening lengths (see Hernquist & Katz 1989 for details). In terms of where the force is 50 percent Newtonian, the equivalent Plummer softening length would be $`0.67`$ times the spline softening length. The softening length of the halo particles is set to $`ϵ=0.05`$, or $`200`$ pc. Particle trajectories are computed using a standard second order symplectic leap-frog integrator, with a maximum time-step $`\mathrm{\Delta }t=1`$ (corresponding to $`3.77`$ Myr in our adopted physical units). Because of the multi-time stepping, some particles are integrated with smaller timesteps. For a typical simulation, approximately $`20`$ percent of the particles are advanced with $`\mathrm{\Delta }t=0.5`$, and $`0.05`$ percent with $`\mathrm{\Delta }t=0.25`$.

The satellite is modeled as a single particle with mass $`M_s`$ and softening length $`ϵ_s`$. Beyond $`2ϵ_s`$ the satellite potential falls as $`r^1`$, so this radius approximates the satellite’s tidal radius, which is mainly determined by conditions at pericenter (King 1962). In principle. we could fix the softening using an appropriate tidal radius. However, the pericentric distance evolves and we would then have to include time evolution of the softening and some mass loss. We opt for a simpler approach and fix the mean density of each satellite:

$$ϵ_s=2.73\mathrm{kpc}\left[\frac{M_s}{10^{10}\mathrm{M}_{}}\right]^{1/3},$$
(34)   
The scaling is set so that a satellite with mass similar to the Large Magellanic Cloud (hereafter LMC) has a softening length comparable to the LMC’s effective radius (de Vaucouleurs & Freeman 1972). This choice is somewhat arbitrary, and we discuss its influence on dynamical friction timescales in Section 3.3.3.

All simulations are run for $`15`$ Gyr on $`2`$ or $`3`$ DEC Alpha processors, each requiring about 48 hours of wallclock time. Energy conservation was typically of the order of one percent over the total length of the simulation.

### 3.3 Results

We determine the eccentricity in two ways. We compute the center of mass of the halo particles and use the galactocentric distance of the satellite ($`r_s`$), and its energy and angular momentum to solve equation (3) for the orbital eccentricity (shown as solid lines in Figures 9 and 10). This eccentricity is only accurate if the potential of the halo has not changed significantly owing to the satellite’s decay. Hence, we also determine the radial turning points of the orbit and compute the approximate eccentricity which we assign to a time midway between the turning points (shown as open circles in Figures 9 and 10).

In Figure 8 we plot the trajectories of the satellites for models 1 to 6. Both the $`x`$$`y`$ (upper panels) and $`x`$$`z`$ projections (smaller lower panels) are shown. The three trajectories plotted on the top vary in satellite mass, whereas those at the bottom vary in their initial orbital eccentricity (see Table 1). The time-dependences of galactocentric radius, eccentricity, energy, and angular momentum for various models are shown in Figures 9 and 10. Energies are scaled by the central potential $`\mathrm{\Phi }_0`$, and angular momenta by their value at $`t=0`$. Eccentricities are only plotted up to $`t_{0.8}`$ (see below), after which the satellite has virtually reached the halo’s center. Models 1 and 2 reveal an almost constant orbital eccentricity. In Models 3 to 9, in which the satellite mass is equal to $`2\times 10^{10}\mathrm{M}_{}`$, the eccentricity reveals a saw-tooth behavior, such that eccentricities decrease near pericenter, and increase near apocenter. This is in perfect agreement with equation (33). It is remarkable that the net effect of $`de/dt`$ is nearly zero: the eccentricity does not change significantly. The alternative definition of eccentricity based on observed turning points (open circles) shows similar results. Small deviations are due to a change of the halo potential induced by the decaying satellite, and the heuristic assignment of a time to the value found from monitoring the radial turning points. The change in energy with time reveals a step-wise behavior, indicating that the pericentric passages dominate the satellite’s energy loss. Note that the energy of the satellites in Models 3 to 9 does not become equal to $`\mathrm{\Phi }_0`$ once the satellite reaches the center of the potential well. This owes to the deposition of energy into the halo particles by the satellite. The details of this process will be examined in a future paper (van den Bosch et al. in preparation).

Because of the elongation of the orbits, the galactocentric distance is not a meaningful parameter to use to characterize the decay times. Instead we use both the energy and the angular momentum. While the energy is well defined, it changes almost stepwise (see Figures 9 and 10). The angular momentum depends on the precise position of the halo’s center which may be poorly determined when the satellite induces an $`m=1`$ mode. We define the following characteristic times: $`t_{0.4}`$, $`t_{0.6}`$, and $`t_{0.8}`$ defined as the time when the satellite’s energy reaches $`40`$, $`60`$, and $`80`$ percent of $`\mathrm{\Phi }_0`$ respectively, and $`t_{1/4}`$, $`t_{1/2}`$, and $`t_{3/4}`$ when the angular momentum is reduced to a quarter, a half and three-quarters of its initial value. These timescales are listed in Table 2. Because of the instantaneous introduction of the satellite in the virialized halo potential the absolute values of these timescales may be off by a few percent. However, we are mainly interested in the variation of the dynamical friction time as function of the orbital eccentricity. We belief that the instantaneous introduction of the satellite will not have a significant influence on this behavior, as this is only a second order effect.

#### 3.3.1 Influence of satellite mass

Models 1, 2, and 3 start on the same initial orbit, but vary in both the mass and size of the satellite. Satellite masses correspond to $`2\times 10^8\mathrm{M}_{}`$, $`2\times 10^9\mathrm{M}_{}`$, and $`2\times 10^{10}\mathrm{M}_{}`$ for Models 1, 2, and 3 respectively. The sizes of the satellites, e.g., their softening lengths, are set by equation (34), such that all satellites have the same mean density. The mass and size of the satellite in model 3 correspond closely to that of the LMC. It is clear from a comparison of Models 1,2, and 3 that dynamical friction by the galactic halo is negligible for satellites with masses $`10^9\mathrm{M}_{}`$. Thus, globular clusters and the dwarf spheroidals in the galactic halo are not expected to have undergone significant changes in their orbital properties induced by dynamical friction by the halo; the only two structures in the galactic halo that have experienced significant amounts of dynamical friction by the halo are the LMC and the SMC, with masses of $`2\times 10^{10}\mathrm{M}_{}`$ and $`2\times 10^9\mathrm{M}_{}`$ respectively (Schommer et al. 1992). Note that we have neglected the action of the disk and bulge, which apply a strong torque on objects passing nearby. This can result in an enhanced decay of the orbit, not taken into account in the simulations presented here.

#### 3.3.2 Influence of orbital eccentricity

Models 3, 4, 5, and 6 differ only in the eccentricity of their initial orbit; satellites on less eccentric orbits are started from smaller apocentric radii, such that the initial energy is the same in each case. The characteristic friction timescales (defined in Section 3.3 and listed in Table 2) as a function of initial eccentricity are shown in Figure 11. The decay time decreases with increasing eccentricity, the exact rate of which depends on the characteristic time employed. The timescales for the most eccentric orbit (Model 3) are a factor $`1.5`$ to $`2`$ shorter than for the circular orbit (Model 6). This exemplifies the importance of a proper treatment of orbital eccentricities, since decay times based on circular orbits alone overestimate the timescales of dynamical friction for typical orbits in a dark halo by up to a factor two.

Using Chandrasekhar’s dynamical friction formula, and integrating over the orbits, Lacey & Cole (1993) found that, for a singular isothermal sphere, the dynamical friction time is proportional to $`\eta ^{0.78}`$ (with $`\eta `$ the orbital circularity defined in Section 2.1.1). The dotted lines in Figure 11 correspond to this dependence, normalized to $`t_{0.6}`$ and $`t_{1/2}`$ for the orbit with an intrinsic eccentricity of $`0.6`$. As can be seen, our results seem to suggest a somewhat weaker dependence of the dynamical friction time on orbital eccentricity. Combining all the different characteristic decay times listed in Table 2, we obtain the best fit to our $`N`$-body results for $`t\eta ^{0.53\pm 0.01}`$. This difference with respect to the results of Lacey & Cole is most likely due to the fact that their computations ignore the global response of the halo to the decaying satellite.

#### 3.3.3 Influence of satellite size

Since the adopted softening length of a satellite, $`ϵ_s`$, is somewhat arbitrary (see Section 3.2), we examine the effect of this choice by running a set of models (3, 7, 8, and 9) that differ only in this quantity. The characteristic friction timescales as a function of $`ϵ_s`$ are shown in Figure 12. In Chandrasekhar’s formalism, the dynamical friction timescale depends on the inverse of the Coulomb logarithm $`\mathrm{ln}\mathrm{\Lambda }`$. This is often approximated by the logarithm of the ratio of the maximum and minimum impact parameters for the satellite. Taking $`b_{\mathrm{max}}=r_t=200`$ kpc and $`b_{\mathrm{min}}=ϵ_s`$, the expected ratios between the timescales of Models 3, 7, 8, and 9 can be determined. The dotted lines in Figure 12 shows these predictions scaled to the characteristic times of Model 3. It is apparent, however, that these simulations are not accurately represented by this simple approximation. This is not surprising as the choices for the minimum and maximum impact parameters are somewhat arbitrary. Because of the radial density profile of the halo, $`b_{\mathrm{max}}`$ is likely to depend on radius. Furthermore, approximating $`b_{\mathrm{min}}`$ by the satellite’s softening length is only appropriate for large values of $`ϵ_s`$; in the limit where the satellite becomes a point mass $`b_{\mathrm{min}}`$ should be taken as $`GM_s/v^2`$, where $`v^2^{1/2}`$ is the rms velocity of the background, a quantity which again depends on the galactocentric distance of the satellite (White 1976b). The predictions in Figure 12 are thus as good as we might have hoped for. The important point we wish to make here is that deviations from the softening lengths given by equation (34), within reasonable bounds, do not significantly modify the characteristic dynamical friction timescales presented here.

#### 3.3.4 Influence of number of halo particles

In order to address the accuracy of our simulations we have repeated Model 4 with both 20,000 and 100,000 halo particles. We henceforth refer to these simulations as Models 10 and 11 respectively (see Table 1). Both these models are first evolved for 10 Gyrs without satellite to let the halo virialize and reach equilibrium. As for the halo with 50,000 particles, the density distribution after 10 Gyrs has not changed significantly.

A comparison of the results of Models 4, 10 and 11 is shown in Figure 13. The solid lines correspond to Model 4, the dotted lines to Model 10 ($`N=20,000`$), and the dashed lines to model 11 ($`N=100,000`$). Clearly our results are robust, with only a negligible dependence on the number of halo particles. The main differences occur at later times, when the satellite has virtually reached the center of the halo. In Models 4 and 11 the satellite creates a core in the halo, which causes $`E/\mathrm{\Phi }_0<1`$ in the center. Model 10 has insufficient particles to resolve this effect, which explains the differences in both the satellite’s energy and eccentricity at later times with respect to models 4 and 11. In Table 2 we list the different characteristic timescales for Model 10 and 11. They are similar to those of Model 4 to an accuracy of $`4`$ percent for Model 10 and $`1`$ percent for Model 11. Finally we note that our main conclusion, that there is no net amount of circularization, holds for different numbers of halo particles.

## 4 Applications

### 4.1 Orbits of Globular clusters

Odenkirchen et al. (1997, hereafter OBGT) used Hipparcos data to determine the proper motions of 15 globular clusters so that all six of their phase-space coordinates are known. Their velocity dispersions show a slight radial anisotropy with $`(\sigma _r,\sigma _\theta ,\sigma _\varphi )=(127\pm 24,116\pm 23,104\pm 20)`$ $`\mathrm{km}\mathrm{s}^1`$. OBGT integrated the orbits using a model for the galactic potential (Allen & Santillan 1991) that approaches an isothermal at large radii. They find a median eccentricity of $`e=0.62`$ and conclude that the globulars are preferentially on orbits of high eccentricity

We can compare their results with the eccentricity distributions expected for a power-law tracer population in a singular isothermal halo. We use the technique described in the Appendix to determine the distribution of orbital eccentricities given a logarithmic slope of the density distribution $`\alpha `$ and an anisotropy parameter $`a`$. We compare the cumulative distribution to OBGT’s sample using the K-S test (e.g., Press et al. 1992). We determine the probabilities that the OBGT sample is drawn randomly from such a distribution using a $`100\times 100`$ grid in $`(\alpha ,a)`$-space with $`\alpha [2,6]`$ and $`a[2,2]`$ and show the contour plot of these probabilities in Figure 14. The contours corresponds to the $`10,20,\mathrm{},80,90`$ percent probability levels, whereby the latter is plotted as a thick contour. Clearly, the slope of the density distribution, $`\alpha `$, is poorly constrained as the dependence is mild (see Section 2.2). However, the velocity anisotropy is well constrained, and we find that a small amount of radial anisotropy is required in order to explain the observed eccentricities of the globular clusters. This is in excellent agreement with the velocity dispersions obtained directly from the data. The solid dot in Figure 14 corresponds to the best fitting model with $`\alpha =3.5`$, the value preferred by Harris (1976) and Zinn (1985). In Figure 15 we plot the cumulative distribution of the eccentricities from the OBGT sample (thin lines) with the cumulative distribution for the model represented by the dot in Figure 14. The K-S test yields a probability of $`94.3`$ percent that the OBGT sample of orbits is drawn randomly from the probability distribution represented by the thick line.

We conclude that the distribution of eccentricities is just what one expects from a population with a mild radial velocity anisotropy; there is no “preference” for high eccentricity orbits as suggested by OBGT. The potential used by OBGT to integrate their orbits deviates significantly from a spherical isothermal in the center, where the disk and bulge dominate. Since the OBGT sample is limited to nearby globulars within approximately 20 kpc of the Sun, the deviations are likely significant. However, the good agreement between the distribution of eccentricities for the OBGT sample, based on the axisymmetric potential of Allen & Santillan (1991), and a spherical isothermal potential suggests that the differences between these two potentials have only a mild influence on the distribution of orbital eccentricities.

### 4.2 Kinematics in $`\omega `$Centauri

Norris et al. (1997) examined the dependence of kinematics on metal abundance in the globular cluster $`\omega `$Centauri (hereafter $`\omega `$Cen) and found that the characteristic velocity dispersion of the most calcium rich stars is $`8\mathrm{km}\mathrm{s}^1`$, while that of the calcium poor stars is $`13\mathrm{km}\mathrm{s}^1`$. The metal rich stars are located closer to the middle where the velocity dispersion is largest and the authors note that there is evidence for rotation in the metal-poor sample (at $`5\mathrm{km}\mathrm{s}^1`$), but not in the metal rich sample. They use all these facts to conclude that “The more metal-rich stars in $`\omega `$Cen are kinematically cooler than the metal-poorer objects”. The metal-rich stars in their sample live in the part of the cluster where the inferred circular velocity is greatest. Hence, they have an average value of $`\stackrel{}{F}\stackrel{}{r}`$ that is greater than for the metal-poor stars. By the virial theorem, their mean kinetic energy must be higher (equation ), yet Norris et al. find just the opposite with the average kinetic energy of a metal-poor star being more than twice that of a metal-rich star.

Since we are not about to abandon the virial theorem, we can only conclude that if the measurements of the dispersions are correct, the kinetic energy per metal-rich star must be at least that of the metal-poor stars implying a much greater kinetic energy in the plane perpendicular to the line of sight than observed along the line of sight. The straightforward way for this to occur is a rotating disk seen face-on. The rotation must be large enough that $`v_{rot}^2+3\sigma ^2`$ is at least as large as seen for the metal-poor stars. This implies that the metal-rich stars have a rotation velocity of $`18\mathrm{km}\mathrm{s}^1`$. The metal-rich component of $`\omega `$Cen must be a rotating disk that is more concentrated than the metal-poor stars. This is exactly the signature that Norris et al. would have ascribed to self-enrichment of the cluster (Morgan & Lake 1989). The rotation signature would be visible as a proper motion of $`0.064^{\prime \prime }`$/century. Rotation of a smaller magnitude has been detected in M22 by Peterson & Cudworth (1994). Note that a face-on disk in $`\omega `$Cen, which is the Galactic globular with the largest projected flattening, implies a triaxial potential.

Norris et al. did not realize that their observations presented this dynamical puzzle. Instead, they believed that the difference could result from the relative radial profiles of the two components as might be seen in the equations of stellar hydrodynamics (equation ); i.e., the density distribution of the metal rich component must fall off more rapidly with radius than for the metal poor component. They used the observations to argue that $`\omega `$Cen was the product of a merger of previous generations of substructure. However, we argue that such a merger would not have the signatures that they see. The metal-poor stars are the overwhelming majority of the stars. Conservation of linear momentum thus implies that the mean radius of the stars that were in the small (metal-rich) object will be greater than those that were in the large (metal-poor) one. Conservation of angular momentum furthermore implies that the rotation velocity of the stars that were in the small (metal-rich) object will be greater than those that were in the large (metal-poor) one. These signatures are exactly opposite of those claimed by Norris et al. to be consistent with the merger model.

### 4.3 Sinking satellites and the heating of galaxy disks

The sinking and subsequent merging of satellites in a galactic halo with an embedded thin disk has been studied by numerous groups (e.g., Quinn & Goodman 1986; Tóth & Ostriker 1992; Quinn, Hernquist & Fullager 1993; Walker, Mihos & Hernquist 1996; Huang & Carlberg 1997). The timescale for merging clearly depends on the eccentricity of its orbit. Once the satellite interacts with the disk, the sinking accelerates. Several of the above studies assumed that circularization owing to dynamical friction by the halo is efficient and examined satellites that started on circular orbits at the edge of the disk. However, we have shown that circularization is largely a myth; satellites will have large eccentricities when they reach the disk. Quinn & Goodman (1986) followed a satellite with a “typical” eccentricity in an isotropic singular isothermal sphere. They derived $`e0.47`$ or $`r_+/r_{}=2.78`$ for this “characteristic” orbit based on some poorly founded arguments. This ratio, however, is significantly smaller than the true median value of $`r_+/r_{}3.5`$; approximately $`63`$ percent of the orbits have $`e>0.47`$. Huang & Carlberg (1997), in an attempt to be as realistic as possible when choosing the initial orbital parameters, started their satellite well outside the disk on an eccentric orbit. However, the eccentricity is only 0.2, clearly too low to be considered a typical orbit.

The eccentricity of the orbit has two effects; more eccentric orbits decay more rapidly in the halo and they touch the disk at an earlier time. The sinking time scales for the satellite caused by its interaction with the disk are more rapid when the difference in velocities between the satellite and the disk stars is smaller; satellites on prograde, circular orbits in the disk decay fastest (see Quinn & Goodman 1986 for a detailed discussion). Thus whereas more eccentric orbits reach the disk sooner, they are less sensitive to the disk interaction because of their high velocities at pericenter. Less eccentric orbits, on the other hand, require a longer time to reach the disk, but once they do their onward decay is rapid (if the orbit is prograde). The exact dependence of the timescales and disk heating on the orbital eccentricities awaits simulations (van den Bosch et al. , in preparation).

We can compare our results to the few satellite orbits determined from proper-motions. Johnston (1998) integrated the orbits of the LMC, Sagittarius, Sculptor, and Ursa Minor in the galactic potential, and found apo- to pericenter ratios of $`2.2`$, $`3.1`$, $`2.4`$, and $`2.0`$ respectively. Using the K-S test, we find that these eccentricities are so small that there is only a $`8.7`$ percent probability that this sub-sample is drawn from the isotropic model of an isothermal sphere. As shown in previous sections the distribution is relatively insensitive to the profiles of the underlying potential and the density distribution of the tracer population. This therefore implies that the velocity distribution is very strongly tangential<sup>8</sup><sup>8</sup>8Even if we adopt the maximum amount of tangential anisotropy allowed by our simple parameterization of $`h_a(\eta )`$ (i.e., $`a\mathrm{}`$) the K-S probability does not exceed 20 percent.. Since we expect that collapsed halos would produce states with dispersions that are preferentially radial, we suspect that the system of galactic satellites has been strongly altered with satellites on more eccentric orbits having been destroyed owing to their faster dynamical friction timescales. In Section 3, we found that more eccentric orbits have smaller dynamical friction timescales, but the effect is only modest (less than a factor two). Furthermore, all the satellites except the Magellanic clouds have masses $`10^9\mathrm{M}_{}`$, and the dynamical friction owing to the halo is almost negligible for these systems. To get the strong effect that is seen, we have to appeal to the Milky Way disk to accelerate the decay and/or tidally disrupt the satellites on the more eccentric orbits. The problem with this solution, however, is that most studies of sinking satellites have shown that they lead to substantial thickening of the disk. Further studies are needed to investigate whether a disk can indeed yield the observed distribution of orbital eccentricities for the (surviving) satellites without being disrupted itself. We are currently using high resolution $`N`$-body simulations to investigate this is detail (van den Bosch et al. , in preparation). Finally, we must note that the sample of satellites with proper motions is small and the proper motions themselves have large errors which bias results toward large transverse motions. Hence, we close this section with the all too common lament of the need for more data with more precision as well as better simulations that include the disk.

### 4.4 Tidal streams

The tidal disruption of satellites orbiting in a galactic halo creates tidal streams (see e.g., McGlynn 1990; Moore & Davis 1994; Oh, Lin & Aarseth 1995; Piatek & Pryor 1995; Johnston, Spergel & Hernquist 1995). These streams are generally long-lived features outlining the satellite’s orbit (Johnston, Hernquist & Bolte 1996). Clearly, a proper understanding of the orbital properties of satellites is of great importance when studying tidal streams and estimating the effect they might have on the statistics of micro-lensing (i.e., Zhao 1998).

An interesting result regarding tidally disrupted satellites was reached by Kroupa (1997) and Klessen & Kroupa (1998). They simulated the tidal disruption of a satellite without dark matter orbiting a halo. After tidal disruption, a stream sighted along its orbit can have a spread in velocities that mimics a bound object with a mass-to-light ratio that can be orders of magnitude larger than the actual, stellar mass-to-light ratio. They show that one can only sight down such a stream if the satellite was on an eccentric orbit. The chance of inferring such a large mass-to-light ratio thus increases with eccentricity. Klessen & Kroupa conclude that the high inferred mass-to-light ratios in observed dwarf spheroidal galaxies could occur from tidal streams (rather than from a dark matter halo surrounding the satellite) if the orbital eccentricities exceed $`0.5`$. We find that $`60`$% of the orbits obey this criterion even without any radial anisotropy (e.g., Figure 3). We can thus not rule out satellites without dark matter based on the distribution of orbits. However, Klessen & Kroupa (1998) neglect to examine the kinematics in detail. The velocity spread owes to a systematic gradient in velocity along the stream. The apocryphal dwarf will appear to rotate rapidly if it is not perfectly aligned with the line of sight, as is clear in the simulations of the Sagittarius dwarf performed by Ibata & Lewis (1998).

### 4.5 Clusters of galaxies

Using very high resolution cosmological $`N`$-body simulations, Ghigna et al. (1998) were able to resolve several hundred subhalos within a rich cluster of galaxies. They examined, amongst others, the orbital properties of these dark matter halos within the larger halo of the cluster, and found that the subhalos followed the same distribution of orbits as the dark matter particles (i.e., those particles in the cluster that are not part of a subhalo). Ghigna et al. report a median apo-to-pericenter ratio of six. We have used the orbital parameters of the subhalos of the cluster analyzed by Ghigna et al. to examine the orbital properties in some more detail. When we only consider subhalos whose apocenters are less than the virial radius of the cluster, $`r_{200}`$, we obtain a sample of $`98`$ halos with a median apo- to pericenter ratio of $`3.98`$. Using the K-S test as in Section 4.1 we obtain a best fit to the distribution of orbital eccentricities for an anisotropy parameter of $`a0.04`$: the virialized region of the cluster is very close to isotropic. Throughout we consider only a tracer population with $`\alpha =2`$, since this is in reasonable agreement with the observed number density of subhalos. Furthermore, the results presented here are very insensitive to the exact value of $`\alpha `$, similar to what we found in Section 4.1. In Figure 16 we plot the cumulative distribution of the orbital eccentricities of the 98 halos with $`r_+<r_{200}`$ (thin line) together with the same distribution for an isothermal sphere with our best fitting anisotropy parameter $`a=0.04`$. The K-S test yields a probability of $`61.3`$ percent that the two data sets are drawn from the same distribution.

When analyzing all subhalos with $`r_+<2r_{200}`$, we obtain a sample of 311 halos with a median apo-to-pericentric ratio of $`4.64`$ and a best fitting value for the anisotropy parameter of $`a=0.27`$. Clearly, the periphery of the cluster, which is not yet virialized, is more radially anisotropic with more orbits on more eccentric orbits.

The $`N`$-body simulations in Section 3.2 are easily scaled to clusters of galaxies, considering the virialized regions of both galaxies and clusters. If we adopt a cluster mass of $`10^{15}\mathrm{M}_{}`$ and take the timescale to be the same as for the Milky Way simulation, the truncation radius becomes $`r_t=2`$ Mpc. In Section 3.3.1 we found that only objects with masses greater than $`0.1`$ percent of the halo mass, $`M_{\mathrm{gal}}10^{12}\mathrm{M}_{}`$ in the cluster, experience significant dynamical friction in a Hubble time. So, only the most massive galaxies experience any orbital decay.

Recently, Moore et al. (1996) showed that high speed encounters combined with global cluster tides—galaxy harassment—causes the morphological transformation of disk systems into spheroids (see also Moore, Lake & Katz 1998). Moore et al. limited themselves to mildly eccentric orbits with $`r_+/r_{}=2`$, but they quoted correctly that this was a low value compared to the typical eccentricity. Their choice was made in order to be conservative and to underestimate the effect of harassment, as the effect increases with orbital eccentricity. They also felt that any larger value would stretch the reader’s credulity as they could refer to no clear presentation of the likely distribution of orbital eccentricities. Our results imply that the effects of harassment were underestimated in their study.

### 4.6 Semi-analytical modeling of galaxy formation

Over the past couple of years several groups have developed semi-analytical models for galaxy formation within the framework of a hierarchical clustering scenario of structure formation (e.g., Kauffmann, White & Guiderdoni 1993; Cole et al. 1994; Heyl et al. 1995; Baugh, Cole & Frenk 1996; Somerville & Primack 1998). The general method of these models is to use the extended Press-Schechter formalism (Bower 1991; Bond et al. 1991; Lacy & Cole 1993) to create merging histories of dark matter halos. Simplified yet physical treatments are subsequently used to describe the evolution of the baryonic gas component in these halos. Using simple recipes for star formation and feedback, coupled to stellar population models, finally allows predictions for galaxies to be made in an observational framework.

A crucial ingredient of these semi-analytical models is the treatment of mergers of galaxies. When two dark halos merge, the fate of their baryonic cores, i.e., the galaxies, depends on a number of factors. First of all, dynamical friction causes the galaxies to spiral to the center of the combined dark halo, thus enhancing the probability that the baryonic cores collide. Secondly, whether or not such a collision results in a merger depends on the ratio of the internal velocity dispersion of the galaxies involved to the encounter velocity. Both depend critically on the masses involved and on the orbital parameters. The dependence on the orbital eccentricities is addressed by Lacey & Cole (1993) who concluded that observations of merger rates and the thinness of galactic disks seem to argue against strongly elongated orbits. This would be problematic in the light of the typical distributions of orbital eccentricities presented here. However, as we pointed out in Section 3.3.2, Lacey & Cole have likely overestimated the dependence of dynamical friction times on orbital eccentricity. Furthermore, as emphasized in Section 4.3 our current understanding of the damaging effect of sinking satellites on thin disks is not well established and may have been overestimated in the past (Huang & Carlberg 1997).

In the actual semi-analytical modeling the merger time-scales are defined by simple scaling laws that depend on the masses only, but that ignore the orbital parameters. The eccentricity distributions presented here, coupled with the $`\eta ^{0.53}`$ dependence of dynamical friction times, may proof helpful in improving the accuracy of the merging timescales in the semi-analytical treatments of galaxy formation. As an illustrative example, we calculate the average dynamical friction time

$$t/t_0=\underset{0}{\overset{1}{}}𝑑\eta \eta ^{0.53}p_a(\eta ),$$
(35)   
where $`t_0`$ is the friction time for a circular orbit and $`p_a(\eta )`$ is the normalized distribution function of orbital circularities in a singular isothermal sphere with anisotropy parameter $`a`$. This average time is plotted as function of $`a`$ in Figure 17 (solid line). For comparison we also plotted the average times obtained by assuming a $`\eta ^{0.78}`$ dependence (dashed line). The average dynamical friction time is typically of the order of $`70`$ to $`80`$ percent of that of the circular orbit. The stronger dependence of Lacey & Cole underestimates the typical friction times by approximately $`10`$ percent.

## 5 Conclusions & Discussion

This paper has presented the distributions of orbital eccentricities in a variety of spherical potentials. In a singular isothermal sphere, the median eccentricity of an orbit is $`e=0.56`$, corresponding to an apo- to pericenter ratio of 3.55. About $`15`$ percent of the orbits have $`r_+/r_{}>10`$, whereas only $`20`$ percent have moderately eccentric orbits with $`r_+/r_{}<2`$. These values depend strongly on the velocity anisotropy of the halo. Collapse is likely to create radially biased velocity anisotropies that skew the distribution to even higher eccentricities. We also examined the distributions of orbital eccentricities of isotropic tracer populations with a power-law density distribution ($`\rho _{\mathrm{tracer}}r^\alpha `$) and found only modest dependence on $`\alpha `$. Due to the unphysical nature of the singular isothermal sphere, we examined more realistic models and found that they differed only slightly from the isothermal case.

We stress that these eccentricity distributions apply only to systems in equilibrium. If a tracer population has not yet fully virialized in the halo’s potential, its orbital eccentricities can be significantly different from the virialized case. Cosmological simulations of galaxy clusters and satellite systems around galaxies show prolonged infall, and recent mergers can produce correlated motions. Hence, care must be taken in applying our results to such systems.

Objects with mass fractions greater than 0.1% experience significant orbital decay owing to dynamical friction. We used high resolution $`N`$-body simulations with $`50,000`$ particles to examine the sinking and (lack of) circularization of eccentric orbits in a truncated, non-singular isothermal halo. We derived, and numerically verified, a formula that describes the change of eccentricity with time; dynamical friction increases the eccentricity of an orbit near pericenter, but decreases it again near apocenter, such that no net amount of circularization occurs. The energy loss owing to dynamical friction is dominated by the deceleration at pericenter resulting in moderately shorter sinking timescales for more eccentric orbits. We find a dependence of the form $`t\eta ^{0.53}`$; the average orbital decay time for an isotropic, isothermal sphere is $`75`$ percent of that of the circular orbit. This dependence is weaker than the predictions of Lacey & Cole (1993), who found $`t\eta ^{0.78}`$ based on analytical integrations of Chandrasekhar’s dynamical friction formula. Since this analytical treatment ignores the global response of the halo to the decaying satellite we believe our results to be more accurate. This relatively weak dependence of decay times on eccentricity, together with the absence of any significant amount of circularization, implies that dynamical friction does not lead to strong changes in the distribution of orbital eccentricities. When scaling the simulations to represent the orbiting of satellites in the galactic halo, we find that the LMC and the SMC are the only objects in the outer Milky Way halo that have experienced significant amounts of energy and angular momentum loss by dynamical friction.

The distribution of orbital eccentricities is important for several physical processes including: timescales for the sinking and destruction of galactic satellites, structure and evolution of streams of tidal debris, harassment in clusters of galaxies, and mass estimates based on the dynamics of the system of globular clusters. The results presented here may proof particularly useful for improving the treatment of galaxy mergers in semi-analytical models of galaxy formation. In Section 4 we showed that the distribution of orbital eccentricities of a subsample of the galactic globular cluster system is consistent with that of a slightly radially anisotropic $`r^{3.5}`$ tracer population in an isothermal potential. A similar result was found for the subhalos orbiting a large cluster of galaxies in a high resolution, cosmological $`N`$-body simulation presented by Ghigna et al. (1998). However, the Milky Way satellites are not consistent with this distribution, but show a bias toward circularity that may have been caused by dynamical friction and/or tidal disruption by the galactic disk. We expect that additional data together with simulations that include the disk will lead to stringent constraints on the formation and evolution of substructure in the Milky Way.

We are grateful to Sebastiano Ghigna for sending us his data in electronic form. We are indebted to Derek Richardson, Jeff Gardner and Thomas Quinn for their help and support with the $`N`$-body simulations, and to the anonymous referee for his helpful comments that improved the paper. FvdB was supported by a Hubble Fellowship, #HF-01102.11-97A, awarded by STScI.     

## Appendix A Monte Carlo method to compute distributions of orbital eccentricities in a singular isothermal potential

To determine the distribution of orbital eccentricities in a spherical potential, one must sample orbits according to the distribution function $`f(E,L)`$. In this appendix, we describe a Monte Carlo method that samples a quasi-separable distribution function (equation ). Once the energy and angular momentum of an orbit are known, the eccentricity is easily determined (Section 2.1). In the following, $``$ is a random number in the interval $`[0,1]`$, $`𝒫`$ denotes a probability distribution, and $`\widehat{f}`$ indicates that the function $`f`$ is normalized.

For a singular isothermal density distribution with a quasi-separable distribution function of the form given by equation (7), the density can be written as

$$\rho (r)=\frac{4\pi }{r}\underset{\mathrm{\Phi }(r)}{\overset{\mathrm{}}{}}𝑑Eg(E)L_c(E)\underset{0}{\overset{\eta _{\mathrm{max}}}{}}h(\eta )\frac{\eta d\eta }{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}}.$$
(A1)   
The joint probability distribution of $`(E,\eta )`$ at a fixed radius $`r_0`$ is therefore

$$𝒫(E,\eta )=\frac{4\pi }{r_0\rho (r_0)}g(E)L_c(E)\frac{\eta h(\eta )}{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}},$$
(A2)   
with $`E>\mathrm{\Phi }(r_0)`$ and $`\eta [0,\eta _{\mathrm{max}}]`$ (cf. van der Marel, Sigurdsson & Hernquist 1997).

Because of the quasi-separable nature of the DF, the probability function for the energies $`E`$ can be separated:

$$𝒫(E)=\frac{4\pi }{r_0\rho (r_0)}g(E)L_c(E)=\frac{r_0}{u_HV_c^2}\mathrm{exp}\left[\frac{E}{V_c^2}\right].$$
(A3)   
The normalized, cumulative probability distribution of the energy is

$$\widehat{𝒫}(<E)=1\mathrm{exp}\left[\frac{\mathrm{\Phi }(r_0)E}{V_c^2}\right].$$
(A4)   
Using the inversion of $`\widehat{𝒫}(<E)`$, energies are selected by drawing a random number $``$ and assigning an energy

$$E=\mathrm{\Phi }(r_0)V_c^2\mathrm{ln}(1)$$
(A5)   
Once the energy is known, the maximum value of $`\eta `$ corresponding to that energy can be calculated using equation (10). For our random number $``$ this yields

$$\eta _{\mathrm{max}}=\sqrt{2\mathrm{e}}(1)\sqrt{\mathrm{ln}(1)}$$
(A6)   
Note that $`\eta _{\mathrm{max}}`$ is independent of the radius $`r_0`$.

In order to draw $`\eta `$ from the probability distribution

$$𝒫(\eta )=\frac{\eta h(\eta )}{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}}$$
(A7)   
we choose the comparison probability distribution

$$𝒫_{\mathrm{comp}}(\eta )=\frac{\eta _{\mathrm{max}}}{\sqrt{\eta _{\mathrm{max}}^2\eta ^2}}$$
(A8)   
which has the property $`𝒫_{\mathrm{comp}}(\eta )𝒫(\eta )`$ for $`0\eta \eta _{\mathrm{max}}`$, as long as $`h(\eta )1`$. Trial values for $`\eta `$ are drawn from the probability distribution $`𝒫_{\mathrm{comp}}(\eta )`$ by inversion of its normalized, cumulative probability distribution,

$$\widehat{𝒫}_{\mathrm{comp}}(<\eta )=\frac{2}{\pi }\mathrm{arcsin}\left(\frac{\eta }{\eta _{\mathrm{max}}}\right),$$
(A9)   
and the rejection method (e.g., Press et al. 1992) is used to decide whether or not this trial value should be accepted.

The eccentricity of a random orbit in an anisotropic, singular isothermal sphere is thus obtained as follows:

1. Draw random numbers $`_1`$ and $`_2`$.
2. Calculate $`\eta _{\mathrm{max}}=\sqrt{2\mathrm{e}}(1_1)\sqrt{\mathrm{ln}(1_1)}`$. If $`_2>\eta _{\mathrm{max}}`$ return to step 1. This takes care of the fact that $`𝒫(E)`$ and $`𝒫(\eta )`$ are not independent: $`𝒫(\eta )`$ depends on $`\eta _{\mathrm{max}}`$ which is a function of energy.
3. Draw random numbers $`_3`$ and $`_4`$.
4. Calculate $`\eta _{\mathrm{try}}=\eta _{\mathrm{max}}\mathrm{sin}\left(\frac{\pi }{2}_3\right)`$. If $`_4>𝒫_{\mathrm{comp}}(\eta _{\mathrm{try}})/𝒫(\eta _{\mathrm{try}})`$ return to step 3
5. Accept $`\eta _{\mathrm{try}}`$ and compute the corresponding eccentricity by solving for the roots of equation (5).

When determining the distribution of orbital eccentricities for a power-law tracer population in a singular isothermal potential (see Section 2.2), one can use the same method as outlined above but with the equation in step 2 replaced by

$$\eta _{\mathrm{max}}=\sqrt{2\mathrm{e}}\sqrt{\beta }(1_1)^\beta \sqrt{\mathrm{ln}(1_1)},$$
(A10)   
where $`\beta =1/(\alpha 1)`$.