A&A 405, 73-88 (2003)
DOI: 10.1051/0004-6361:20030596

Slow evolution of elliptical galaxies induced by dynamical friction

I. Capture of a system of satellites

G. Bertin - 1,2 T. Liseikina - 2,3 - F. Pegoraro 3,[*]


1 - Università degli Studi di Milano, Dipartimento di Fisica, via Celoria 16, 20133 Milano, Italy
2 - Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy
3 - Università degli Studi di Pisa, Dipartimento di Fisica, via Buonarroti 2, 56100 Pisa, Italy

Received 6 September 2002 / Accepted 10 April 2003

Abstract
The main goal of this paper is to set up a numerical simulation for the study of the slow evolution of the density and of the pressure tensor profiles of an otherwise collisionless stellar system, as a result of the interactions with a minority component of heavier "particles". The effects that we would like to study are those attributed to slow collisional relaxation and generically called "dynamical friction"; in real cases, or in numerical simulations, the processes involved are complex, so that the relaxation associated with the granularity in phase space is generally mixed with and masked by evolution resulting from lack of equilibrium or from a variety of instabilities and collective processes. We start by revisiting the problem of the sinking of a satellite inside an initially isotropic, non-rotating, spherical galaxy, which we follow by means of N-body simulations using about one million particles. We then consider a quasi-spherical problem, in which the satellite is fragmented into a set of many smaller masses with a spherically symmetric initial density distribution. In a wide set of experiments, designed in order to bring out effects genuinely associated with dynamical friction, we are able to demonstrate the slow evolution of the density profile and the development of a tangentially biased pressure in the underlying stellar system, while we briefly address the issue of the circularization of orbits induced by dynamical friction on the population of fragments. The results of the simulations presented here and others planned for future investigations allow us to study the basic mechanisms of slow relaxation in stellar systems and thus may be of general interest for a variety of problems, especially in the cosmological context. Here our experiments are conceived with the specific goal to clarify some mechanisms that may play a role in the evolution of an elliptical galaxy as a result of the interaction between the stars and a significant population of globular clusters or of the merging of a large number of small satellites.

Key words: galaxies: evolution - galaxies: interactions - galaxies: dynamics

1 Introduction

If we consider a globular cluster of mass $m_{\rm gc} \sim 10^5~ M_{\odot}$ orbiting inside a spherical galaxy, given the mass ratio 105 with respect to the mass of a typical star, we find that it should suffer dynamical friction on a time-scale $T_{\rm fr} \sim 10^{9}$ yr, through scattering of the stars. (The term dynamical friction refers to the slow relaxation process associated with discrete two-body encounters investigated in the pioneering studies of Chandrasekhar 1943. In these classical analyses, estimate is given of the cumulative effect of discrete encounters on a test particle passing through a homogeneous system of field particles. This effect, associated with the granularity of the system in phase space, is ignored in pure Vlasov descriptions of stellar systems. The term dynamical friction is often extended to describe loosely the physical origin of orbital decays of satellites of finite mass also in real inhomogeneous systems and simulations, which fall obviously outside the classical description of Chandrasekhar.) If the cluster starts from a circular orbit located at ri, in a time $\sim$ $T_{\rm fr}$ it will reach a lower orbit close to $r_{\rm f}$, while a fraction of the energy initially associated with such a cluster $E_i = m_{\rm gc}\Phi(r_i)+ (1/2) m_{\rm gc}r_i ({\rm d}\Phi/{\rm d}r)_{r_i}$ is released into the distribution function of the hosting galaxy. If we refer to the realistic possibility of the presence of several thousands of such globular clusters we may thus argue that on a fraction of the Hubble time the underlying stellar system may well absorb a non-negligible fraction of its binding energy and thus should slowly evolve (see also Lotz et al. 2001).

Under the suggestive scenario just outlined, the issue we would like to address in this paper is the following conceptually challenging problem. In a two-component system, where the more massive component (the "galaxy") is collisionless and the minority component (the globular cluster system or a large number of mini-satellites caught in minor mergers) is slowly dragged in toward the galaxy center by processes of dynamical friction with the stars of the galaxy, how does the distribution function of the galaxy evolve during the process? A resolution of this challenging problem requires a clarification of the basic mechanism of dynamical friction under inhomogeneous conditions and a discussion of the collective effects involved in the reaction of the underlying distribution function. Note that while the infall of a single cluster (or satellite) is a problem with a preferred direction (given by the angular momentum of the satellite), the capture of a large number of randomly distributed clusters (or mini-satellites) can preserve spherical symmetry, if present to begin with.

It might be argued that the slow evolution in the quasi-isotropic many-cluster case described above could be modeled in terms of an adiabatic evolution of a distribution function within spherical symmetry. For the different problem of the adiabatic growth of a central massive black hole, such a procedure has indeed been developed (e.g., see Cipollina & Bertin 1994 and references therein). Here it is not clear how a semi-analytic procedure to calculate such slow evolution, taking advantage of the conservation of the relevant actions, can be formulated. In addition, given the relation between dynamical friction and resonant interactions (e.g., see Weinberg 1989), in spite of its slowness the process is likely to be inherently non-adiabatic. Therefore, it seems that the only viable way to an answer is that of approaching the problem by means of suitable N-body simulations.

In the long run, results from this project could also shed light on another theoretically interesting scenario, which so far has found major applications only in the context of the internal dynamics of globular clusters. This is a picture, pioneered by Bonnor (1956) in his analysis of gas spheres, by which stellar systems may be subject to a process of gravitational catastrophe (Lynden-Bell & Wood 1968), when their central concentration exceeds a certain threshold value. The catastrophe consists in an instability process, where evolution leads to cores that become more and more concentrated. The source of this instability is thermodynamical, in the sense that it results from a counter-intuitive property of self-gravitating systems, which can be described in terms of negative specific heat. Recently, this general scenario has been demonstrated to hold for an astrophysically interesting family of partially relaxed stellar systems (Bertin & Trenti 2003).

The study of relaxation in N-body simulations and by means of N-body simulations is a vast field of research (e.g., see Hohl 1973 and Huang et al. 1993). The study of the sinking of a single satellite inside a galaxy, as a result of dynamical friction, was undertaken systematically about two decades ago, initially through simplified simulations, within a restricted three-body problem and with the center of mass of the primary galaxy being kept fixed, and then by means of more realistic, self-consistent simulations (Bontekoe & van Albada 1987; Zaritsky & White 1988; Bontekoe 1988). The results showed that the general scaling predicted by the classical formulae derived in the "homogeneous limit" (or "local limit") (Chandrasekhar 1943) was basically applicable, although the sinking time scale turned out to be significantly longer (by a factor of 2-3 for the regimes studied by the numerical experiments) if the self-consistent response of the galaxy was properly incorporated. A convincing interpretation of the sinking process can be made, at least in the limit where the mass of the sinking satellite is small, by investigating the linear response in the galaxy distribution function through stellar dynamical techniques (Fridman & Polyachenko 1984; Palmer & Papaloizou 1985; Weinberg 1986, 1989; Bertin et al. 1994, hereafter BPRV94; Palmer 1994). Such response consists of an "in-phase" contribution, which has little relevance to the mechanism responsible for dynamical friction, and an "out-of-phase" contribution, associated with the "resonant" interaction, which produces a wake able to set a significant torque on the satellite. These techniques demonstrate that the collective wake in the galaxy, largely dominated by the dipolar (l = 1) response, can trail significantly away from the orbiting object (when the satellite is outside the galaxy, the response is concentrated around a location opposite to the object with respect to the global center of mass), so that the resulting torque on the satellite that is responsible for dynamical friction can be significantly reduced with respect to other more naive calculations. Discrepancies between earlier simulations and the predictions of the stellar dynamic semi-analytic theory were ascribed to non-linearity effects (in the sense that the satellite-to-galaxy mass ratio $M_{\rm s}/M_{\rm G} \approx 0.1$ considered in the simulations was judged to be excessive; see Hernquist & Weinberg 1989), but the general picture appears to have been clarified. It is noted that the differences in behavior with respect to the naive expectations of the traditional formulae for dynamical friction become substantial when the size of the satellite is finite (Weinberg 1989 refers to the length ratio $R_{\rm s}/R_{\rm G}$; the multipoles associated with the perturbation with l above a certain threshold are bound to be unimportant), while the differences should become negligible for very compact satellites (in which case all multipoles contribute, with a divergence that is removed by the Coulomb logarithm).

Several aspects of the general problem of dynamical friction have been revisited recently. The issue of the contrast between "local" and "non-local" effects has been re-discussed (Prugniel & Combes 1992; Maoz 1993; Cora et al. 1997), with the confirmation that the description of the friction process in terms of the "homogeneous limit" (Chandrasekhar 1943) should be adequate in the limit of low-mass, relatively compact satellites. A very recent investigation (Cora et al. 2001) also suggests that the process depends very little on the regularity or the stochasticity of the stellar orbits in the galaxy, in contrast with some earlier conjectures (Pfenniger 1986). Other interesting problems are the problem of circularization (of the satellite orbit), the relation of the friction effect to the past history of the dynamical event under consideration, and the issue of a direct calculation of the relevant dissipative force from the response of the stellar system to the perturbing satellite (e.g., see Séguin & Dupraz 1994,1996; see also Colpi 1998; Colpi & Pallavicini 1998; Colpi et al. 1999; Nelson & Tremaine 1999).

Note that so far investigations of processes of dynamical friction have mostly focused on the issue of the fate of the sinking objects (and mostly for the case of one object; but see Tremaine et al. 1975), while the problem of the reaction of the distribution function of the underlying stellar system (especially in the quasi-spherical limit of a system of many heavy objects outlined above) remains practically unexplored.

Our paper is organized as follows. In Sect. 2 we summarize the theoretical framework that allows us to understand the mechanisms of dynamical friction within the context of a single sinking satellite. In Sect. 3 we describe the model adopted. For the present exploratory investigation, the choice of the model (both for the galaxy and for the satellite or the fragments) is mostly dictated by convenience (for the study of dynamical mechanisms and for an easier comparison with previous analyses) rather than by realistic requirements. In Sect. 4 we describe the diagnostic tools that will be used in the simulations and some preliminary tests. In Sect. 5 we proceed to illustrate the results of the simulations for the case of one sinking satellite and for the quasi-spherical case for which the sinking of many smaller mini-satellites is considered. In Sect. 6 we briefly draw our conclusions and set out the plans for future investigations.  

   
2 Theoretical framework for the problem of a single sinking satellite

Here we briefly recall the main issues involved in the problem of a single sinking satellite of "diameter" $2 R_{\rm s}$, initially located on a circular orbit at distance $r_{\rm s} \approx R_{\rm G}$ from the center of mass of the entire system. Here $R_{\rm G}$ denotes the "radius" of the galaxy. Let $\Omega_{\rm s} > 0$ be the angular velocity of the satellite along its circular orbit. The case of a very compact satellite is likely to be within the reach of the classical formulae for dynamical friction (Chandrasekhar 1943). In the case where the satellite is sufficiently extended, for the description of its density distribution and the associated gravitational potential, only the first harmonics in a multipole expansion are important, up to $l_{\max} \sim \pi r_{\rm s}/(2 R_{\rm s})$ (see Weinberg 1989).

A "friction" torque acting on the satellite is expected to result from the action of suitable resonances with star orbits. To understand this, we may refer to the action of the satellite as a perturbation over a basic state that is non-rotating, spherically symmetric, and characterized by isotropic distribution function f0 (E), where $E = v^2/2 + \Phi_0 (r)$ is the star energy, per unit mass, in the absence of the perturbation, and we may then follow the analysis developed earlier (BPRV94).

Analytical results can be derived in the limit when the mass of the satellite is vanishingly small (with respect to the mass of the galaxy, so that the resulting perturbation in the galaxy distribution function can be treated by a linear theory; see Sect. 2.2). In practice, for real systems or for numerical simulations, some of the evolution effects are going to be associated with the finite mass of the satellite. Once we move to consider the possibility of a quasi-equilibrium configuration for a composite system (galaxy plus satellite or a population of mini-satellites), we then have to worry about genuine instabilities that may occur, independent of the granularity present in the system. All of this will contribute to evolution although, in principle, it is not related to the original gedanken experiment at the basis of the classical definition of dynamical friction.

Another subtle aspect of the problem is related to geometry. Broadly speaking, we may think that the difference between the case of straight orbits (considered in the classical description of dynamical friction) and that of orbits in the inhomogeneous potential of the galaxy is taken into account in the following discussion of the resonant response. Yet, the cumulative effects on a single star, because of multiple encounters with the satellite (for example, for stars on quasi-circular orbits), grow in time differently from the classical picture, even before dealing with the collective behavior of the entire set of stars that characterizes the resonant response.

   
2.1 Single-star resonances

The frequency $\omega$ associated with the mth component of a multipolar term of order l in the expansion of the potential $\Phi_{\rm s}$ of the satellite is readily related to the angular frequency of the satellite $\omega = - m \Omega_{\rm s}$. Therefore, if we refer to Eq. (80) of BPRV94, we find that the appropriate resonance condition is $ - m \Omega_{\rm s} - p \Omega_r(E,J) + s \Omega_{\theta}(E,J) = 0$; here s, m, p are integers ( $- l \le s,m \le +l$; $- \infty < p < \infty$) and $\Omega_r(E,J) > 0$, $\Omega_{\theta}(E,J) > 0$ are the two frequencies (dependent on the star specific energy E and angular momentum J) that characterize the star orbit and depend on the properties of the basic state. As discussed in BPRV94, only the components with s = l, l-2, .. -l contribute to the density perturbation.

For an extended satellite, the p = 0 component, which corresponds to an average along the star radial orbit excursion, is expected to dominate so that the resonance condition should reduce to $\Omega_{\theta}(E,J) = m \Omega_{\rm s}/s$. Note that in the course of time, because of dynamical friction, the angular frequency of the satellite is going to change. Clearly, the resonance conditions that we have just stated are applicable provided that $\dot{\Omega}_{\rm s}/\Omega_{\rm s} \ll \Omega_{\theta}(E,J)$, which we can rewrite as $\dot{\Omega}_{\rm s} \ll \Omega_{\rm s}^2/l$.

   
2.2 "Shielding" and global resonances

If we consider the linearized Vlasov-Poisson equation for our problem, we see that Eq. (54) of BPRV94 is modified into $f_1 = \mathcal{L}(\Phi_{\rm s} + \Phi_1)$, where $\mathcal{L}$ is a linear operator that basically describes one integration along the unperturbed characteristics. The quantity f1 can be called the linear response of the distribution function. This response has a "driven" part, associated with  $\Phi_{\rm s}$, and a "self-gravitating" part, related to the potential $\Phi_1$, that must be consistent with the response itself, following the Poisson equation $\nabla^2 \Phi_1 = 4 \pi G \int f_1 {\rm d}^3 v$.

In plasma physics, the latter contribution is often called the "shielding" associated with the collective behavior of the plasma. Such shielding is often ignored, although it should be emphasized that this step is improper, because the omitted term is a linear contribution. A relatively simple analysis can then be made, by neglecting the self-gravity of the response, showing the connection between single-star resonances and dynamical friction. Note that this is the only ground where a comparison with Chandrasekhar's (1943) work can be made.

The general case, which includes the possibility of resonance with discrete global modes of the system, is currently beyond the reach of analytical investigations (see also Vesperini & Weinberg 2000).

   
3 Model

As a model of the primary stellar system (the "galaxy"), following Bontekoe & van Albada (1987) we consider an isotropic polytropic model of index n = 3. The distribution function is thus of the form $f_{\rm G}(\vec{r},\vec{v}) = A(E_0-E)^{n-3/2}$ for E<E0, with $f_{\rm G}(\vec{r},\vec{v}) = 0$ for E>E0, where $E = \Phi_G(r) + v^2/2$ is the single-star energy, $E_0 = \Phi_G(R) = - GM/R $ is the gravitational potential at the boundary r=R, and M is the total mass of the system. Polytropes can be described in terms of the Lane-Emden function $\theta(\xi),$ which satisfies the differential equation:

 \begin{displaymath} \frac{1}{\xi^2} \frac{\rm d}{{\rm d} \xi} (\xi^2 \frac{{\rm d} \theta}{{\rm d} \xi}) = - \theta^n, \end{displaymath} (1)

where n is the polytropic index, with the boundary conditions $\theta = 1$, ${\rm d} \theta / {\rm d} \xi = 0$ at $\xi = 0$. From the function $\theta(\xi)$ one can find the potential $\Phi_{\rm G}(r)$and the associated density $\rho_{\rm G}(r)$.

The secondary object (the "satellite") is described by a rigid Plummer density distribution (corresponding to a polytrope of index n = 5), characterized by $\theta = (1 + 0.5 \xi^2)^{-1/2}$. The potential and the mass distribution of the satellite are then defined by two scales, $R_{\rm s}$ and mass $M_{\rm s}$, so that $\Phi_{\rm s} (r) = - G M_{\rm s} (R_{\rm s}^2 + r ^2 )^{-1/2}$, with cumulative mass $M_{\rm s}(r) = M_{\rm s} r^3(R_{\rm s}^2 + r^2 ) ^{-3/2}$. This distribution extends to infinity, but 90% of the mass is contained within  $3.7 R_{\rm s}.$

As a result of the interaction with the satellite, the distribution function of the galaxy is a time-evolving f that will be associated with a mean field $- \nabla \Phi$, to be derived from the Poisson equation $\Delta\Phi(\vec{r},t) = 4 \pi G \int f {\rm d}^3 v = 4 \pi G \rho(\vec{r},t)$. We solve the Poisson equation on a spherical polar mesh of $24\times 12\times 24$ cells in the $r,\theta$, and $\varphi$ directions. The radial mesh is logarithmic. The angular mesh is fixed, with 12 equal divisions in $\cos\theta$ and 24 in $\varphi.$ On this mesh the density, the potential, and the force field are expanded in associated Legendre functions $P_l^m(\cos\theta)$ up to P66. At each time-step the center of the spherical mesh used to calculate the potential and the force field in the galaxy is repositioned so as to coincide with the center of mass of the galaxy. The equations of the motion (for the particles representative of the galaxy and, separately, for the satellite) are integrated in Cartesian coordinates with the so called leapfrog scheme. This is basically the N-body code described by van Albada & van Gorkom (1977) and then improved and used by van Albada (1982) and by Bontekoe & van Albada (1987). Therefore, we do not provide here a full description of the method employed, since it is available in the papers quoted.

   
3.1 Sinking of a single satellite

A first set of numerical experiments addresses the sinking of a single satellite, initially located on a circular orbit at $\theta = \pi/2$.

The galaxy is considered as a collisionless stellar system, so that its representative particles (i = 1,..., N) interact with each other through the mean gravitational potential $\Phi$ and directly with the satellite:

 \begin{displaymath} \frac{{\rm d}^2 \vec{r}_i}{{\rm d} t^2} = - \nabla\Phi(\vec{... ...s})}{[R_{\rm s}^2 + (\vec{r}_i-\vec{r}_{\rm s})^2]^{3/2}}\cdot \end{displaymath} (2)

Note that $\Phi$ is updated separately by the Poisson-solver scheme, following the original method described by van Albada & van Gorkom (1977) (this step effectively closes the system of equations), which guarantees a pure time-reversible evolution.

In turn, the satellite is taken to interact with all the representative particles of the galaxy directly via two-body forces:

 \begin{displaymath} \frac{{\rm d}^2 \vec{r}_{\rm s}}{{\rm d}t^2} = \sum_{i=1}^{i... ...{\rm s}^2 + (\vec{r}_i - \vec{r}_{\rm s})^2\right]^{3/2}}\cdot \end{displaymath} (3)

In the equations above, $\vec{r}_i, m$ are the position vectors and the mass of the representative particles of the galaxy, while $\vec{r}_{\rm s}$ is the position vector of the satellite. All position vectors are measured with respect to an inertial frame of reference with its origin at the center of mass of the combined (galaxy + satellite) system.

   
3.2 Sinking of many fragments

Another set of experiments studies an initial configuration where the satellite is fragmented into several ($N_{\rm f}$) pieces ("mini-satellites" or "fragments"). For simplicity, we take the fragments to be characterized by the same Plummer density distribution with mass $M_{\rm f}$ (such that $M_{\rm s} = N_{\rm f} M_{\rm f}$) and the same lengthscale $R_{\rm f} = R_{\rm s}$ (a smaller lengthscale $R_{\rm f}$ would be more realistic for the astronomical system we have in mind, but would be less easy to simulate).

The mutual interactions between mini-satellites can be kept "on'' by considering

 \begin{displaymath} \frac{{\rm d}^2 \vec{r}_i}{{\rm d}t^2} = - \nabla\Phi(\vec{r... ...)} {[R_{\rm f}^2 + (\vec{r}_i - \vec{r}_{{\rm f}j})^2]^{3/2}}, \end{displaymath} (4)

for the representative particles of the galaxy (labeled by i = 1, ...., N), and

 \begin{displaymath} \frac{{\rm d}^2 \vec{r}_{{\rm f}j}}{{\rm d}t^2} = \vec{a}_j^{(1)} + \vec{a}_j^{(2)} \end{displaymath} (5)

for the mini-satellites (labeled by $j = 1,..., N_{\rm f}$), where

 \begin{displaymath} \vec{a}_j^{(1)} = \sum_{i=1}^{i=N} \frac{G m (\vec{r}_i - \v... ...R_{\rm f}^2 + (\vec{r}_i - \vec{r}_{{\rm f}j})^2\right]^{3/2}} \end{displaymath} (6)

is the acceleration produced by the galaxy and $\vec{a}_j^{(2)}$ denotes the gravitational acceleration exerted on the jth mini-satellite by the other fragments. The gravitational interaction between two extended bodies is complicated to describe. For the purposes of the present paper, where we are not interested in the phenomena associated with the details of such interaction, we simplify the simulations by referring to the following prescription:

 \begin{displaymath} \vec{a}_j^{(2)} = \sum_{k =1,k \ne j}^{k=N_{\rm f}}\frac{G M... ...(\vec{r}_{{\rm f}k} - \vec{r}_{{\rm f}j})^2\right]^{3/2}}\cdot \end{displaymath} (7)

All position vectors are measured with respect to an inertial frame with its origin at the center of mass of the entire system.

The mutual interaction between mini-satellites can be ignored and turned "off" simply by taking $\vec{a}_j^{(2)} = 0$. This procedure may be the more realistic procedure if we wish to simulate a system of globular clusters inside a galaxy. In fact, for such system, the few-body-problem effects associated with such interaction are expected to play a secondary role.

As to the initial distribution of the fragments we have considered three cases, to be described below. The first is the most natural generalization of the study of the sinking of a single satellite. A posteriori, we have found that in such a case the effects related to the lack of equilibrium in the initial conditions are too strong and confuse the picture of evolution. (Similar problems, related to the lack of equilibrium, were present also in the case of a single satellite, but those are not emphasized in this paper because that case has been well studied in the past and we prefer to focus here on the new quasi-spherical problem in the presence of many fragments.) Therefore we have devised quieter starts so as to better disentangle the effects associated with the granularity of the system from those associated with the lack of equilibrium.

Given the relatively small number of fragments involved, we anticipate that different realizations of the physical cases described below may be associated with different evolutions.

   
3.2.1 I: Fragments quasi-isotropically distributed on circular orbits on a thin shell
We first consider the case where all the fragments are distributed initially with random positions on a thin shell, located at radius $r_{\rm shell}(0)$. The fragments are given initial velocities appropriate for circular orbits, with random orientations, based on the gravitational field generated by the galaxy alone.

The purpose of simulations of this type is to derive information relevant for a more realistic situation where the fragments have a full three-dimensional space distribution, so that the interaction between the galaxy and the minority population is differential with radius. Such ideally desired simulations are not feasible with current codes. To move one step further in the desired direction we have then considered the case of two shells of fragments initially located at different radii.

Unfortunately, due to the finite mass of the set of fragments, these simulations are characterized by violent epicyclic oscillations that shake the system on the fast dynamical time scale and tend to persist long after the initially expected transient.

   
3.2.2 II: Fragments quasi-isotropically distributed on circular orbits on a thick shell with improved initial conditions
In order to set up a smoother quasi-equilibrium initial condition, we have decided to proceed as follows. We refer to a spherically symmetric density distribution for a shell of the form $\rho_{\rm shell} = \rho_0 \exp[-(r-r_{\rm shell}(0))^2/a^2]$ for $ \vert r-r_{\rm shell}(0)\vert \le 2 a$ (and vanishing for $\vert r-r_{\rm shell}(0)\vert> 2 a$). The normalization constant $\rho_0$ is taken in such a way that the total mass associated with the shell is $M_{\rm s}$. In view of the choice of model for the fragments, we have decided to take $2 a = R_{\rm s}$.

Such density distribution generates a potential $\Phi_{\rm shell}(r)$, which can be calculated numerically. We now consider a modified distribution function for the galaxy $f_{\rm G}(\vec{r},\vec{v})$ taken to be of the polytropic form as in Sect. 3, but with $E = \Phi_{\rm G}(r) + \Phi_{\rm shell}(r) +v^2/2$. We then have to integrate a Lane-Emden equation similar to Eq. (1) but modified by the presence of the external potential $\Phi_{\rm shell}(r)$. This yields a potential $\Phi_{\rm G}(r)$ different from that considered in Sect. 3, to be inserted into the definition of $f_{\rm G}(\vec{r},\vec{v})$.

Finally, we consider the clumpy realization of the shell density distribution by reducing it to a set of $N_{\rm f}$ fragments of mass $M_{\rm f}$. The jth fragment is initially distributed in r at a location  $r_{{\rm f}j}$ (with random angular coordinates) so that the set of fragments reproduces the assumed density distribution of the shell; its initial velocity is that appropriate for the circular velocity of a test particle in the combined potential $\Phi_{\rm G}(r) + \Phi_{\rm shell}(r)$.

This procedure ensures a quieter start while leaving the basic picture of the single shell case unchanged. It can be generalized to the case of two shells of fragments initially located at different radii.

   
3.2.3 III: Fragments constructed by clumping part of the galaxy distribution function

The cases considered above are aimed at describing the interaction between the galaxy and a minority component with phase-space properties distinct from those of the galaxy. As we have seen while addressing the need for a quieter start, such distinction is likely to be accompanied by effects that are due to the lack of equilibrium in the initial conditions, difficult to separate from the secular effects that are traditionally called dynamical friction. In this assessment, we can go even one step further and we may argue that some of these initial conditions (with the different phase-space properties considered) may actually be characterized by dynamical instabilities. If this were the case, it would be very hard or even impossible to distill the true effects of dynamical friction out of the simulations. In order to address this point, we will run parallel simulations of case II with a smooth, Vlasov realization of the shell density, which should allow us to single out the effects of possible Vlasov instabilities of such system (see Sect. 3.2.4).

In view of this discussion, we have decided to address a third set of simulations that is less relevant for the astrophysical case that has motivated this paper but is definitely more interesting for the study of the problem of dynamical friction.

Here we consider the initial equilibrium state to be basically that described by the galaxy model alone, without satellites or shells or fragments, as defined in the first paragraph of Sect. 3. The mini-satellites here are then introduced by clumping together part of the distribution function $f_{\rm G}$ in such a way that it corresponds to a diffuse "shell" in space. Formally, we separate out a piece of the distribution function $f_{\rm shell}$ (for simplicity, we keep the notation "shell" even if we are considering a situation different from that of case I and case II). Then the system made of $f_{\rm G} {-} f_{\rm shell}$ is simulated by the Vlasov code, while the contribution $f_{\rm shell}$ is reduced to clumps that are brought back to the form of the fragments and treated by means of direct interaction with the galaxy representative particles (as described in Sect. 3.2). Note that in this case the "clumps" or "fragments" are distributed in phase space exactly as the particles of the galaxy (thus minimizing the lack of equilibrium in the initial conditions and the risk of introducing undesired sources of instability); in particular, the clumps are initially characterized by an isotropic distribution in velocity space.

We should emphasize that this third class of simulations, while addressing the settling of heavy masses dragged inwards by dynamical friction, is conceptually distinct from the study of the mass stratification process in self-gravitating systems, which may be performed by means of direct or Fokker-Planck codes (e.g., see Spitzer 1987).

   
3.2.4 Purely collisionless reference models

For the three cases just outlined, the simulations carried out following the description given in Sect. 3.2 will be accompanied by reference simulations of similar but purely collisionless models. These can be seen as cases for which the number of fragments formally diverges, while the total mass contained in the fragments remains constant. In the simulations the case of very large number of fragments cannot be treated in terms of direct interaction; instead, for these cases, the fragments are treated as simulation particles of the collisionless component. It is clear that for case III, in such reference collisionless models, the identification of $f_{\rm shell}$ is meaningless. In practice, to monitor the properties of our integration scheme, we have labeled the representative particles associated with $f_{\rm shell}$ by a different "color".

The purpose of these parallel reference runs is to check, as much as possible, that we are not confusing evolution effects associated with the initial conditions with the effects that are genuinely associated with the granularity of the system.

   
3.2.5 Initial set up

In order to derive the initial distribution of the particles in the simulation we start from the mass distribution M(r) which is known on a discrete set of points labelled by the index j, with Mj = M(rj).

First we determine the distance of a particle from the galaxy center. If all the particles of the galaxy have the same mass m = M/N, then for the ith particle we define $xx = m \cdot (i - 1/2),$ and then find the cell j such that Mj>xx>Mj-1.The distance from the center ri of the ith particle is then found by an interpolation between rj-1, rj, and rj+1. Next we choose two uniformly distributed random numbers: $q_{\varphi}\in[0,1]$ for the $\varphi$ coordinate of a particle, such that $\varphi=2\pi q_{\varphi}$, and $q_{\theta}\in[-1,1]$ for the $\theta$ coordinate, such that $\cos\theta = q_{\theta}.$ The Cartesian coordinates are then found from the spherical coordinates by standard formulae. The choice of the initial velocities is performed in a similar way, using a rejection method (e.g., see Press et al. 1992) to obtain the desired distribution function. For the potential $\Phi_{\rm G}(r)$ we proceed as for the mass distribution; the value of the potential at the radius of the particle location is found by interpolation between the values at rj-1, rj, and rj+1.

In case III, of a "subtracted" shell, the particles of the galaxy are divided in two "species": light particles, with mass $m= M/N =(M-M_{\rm s})/N_1$, and heavy particles with mass $M_{\rm f} = M_{\rm s}/N_{\rm f}$respectively, where $M_{\rm s}$ is the mass of the shell, $N_{\rm f}$ the number of heavy particles, and N1 the number of light particles. Initially, the heavy particles are all located inside a thin shell. Inside this shell there are no light particles. The ith light particle, for $i\in[1,i1=N\cdot M(r_{\rm shell})/M],$where $r_{\rm shell}$ is the position of the "shell", is located as in the non-subtracted case at the distance ri obtained by interpolation between the rj-1, rj, and rj+1, where jis such that the condition $M_{j}>m\cdot (i-1/2)>M_{j-1}$ is satisfied. Then, for the kth heavy particle, with $k\in[1,N_f]$, we adopt $xx = m \cdot (i1 - 1/2) + M_{\rm f} \cdot (k-1/2) $ and then find the cell j where Mj>xx>Mj-1. The distance from the center rk of the kth heavy particle is then found by the interpolation rule mentioned above. After that the remaining light particles are distributed according the previous rule. Finally, in the simulations, the light particles interact with one another through their mean field, while the heavy particles interact directly with the light particles and among themselves, as described at the beginning of Sect. 3.2.

   
3.2.6 Other possible initial conditions
Another possible initial condition is to turn on slowly the gravity associated with the massive objects or massive shells. In a collisionless simulation, it would lead to the "empirical construction" of equilibrium systems in which dynamical friction is absent and would thus be, for our purposes, essentially equivalent to the procedure described in Case II. Such construction of equilibrium models can also be approached analytically (see Cipollina & Bertin 1994 and references therein).

The case where the satellite or the fragments are treated as discrete objects would provide information on the friction force, different from the information provided in this paper, but such approach would have no special merits, in terms of physical justification and physical interpretation, with respect to the one adopted here.

   
3.3 The choice of the code

The choice of code that we have made (basically, the code used by Bontekoe & van Albada 1987) appears to be appropriate for the study of the effects of dynamical friction, which are dominated by low-l wakes (see Sect. 2). Still, one may ask the reason why we have resorted to that code instead of codes that are commonly utilized for a variety of studies in galactic dynamics or in the cosmological context. In particular, one popular class of codes, the tree-codes (e.g., see Barnes & Hut 1986), has found modern realizations (e.g., GADGET; see Springel et al. 2001) that are particularly appealing since they are flexible, widely used, and have gone through many tests. Tree-codes treat near particles in terms of direct interactions via a softened potential, because the use of the true Newtonian potential would introduce a generally undesired collisionality among the simulation particles. Tree codes are often used to study the evolution of collisionless systems; for this purpose, the softening length has to be chosen in an optimal way, depending on the physical characteristics of the system that is being investigated.

The use of a tree code in our project would then pose severe interpretation problems, because the non-physical effects associated with the introduction of the softening parameter would show up precisely in relation to the issues that we are trying to investigate.

   
4 Diagnostics and test cases

The results of our simulations can be studied in terms of the energy and angular momentum budget as a function of time and by means of figures that illustrate the spatial structure of the density and of the pressure distribution within the galaxy in the course of time. Another set of figures will provide plots of the azimuthally averaged density profile for the galaxy and of the location $r_{\rm s}$ of the satellite as a function of time; for the case of many mini-satellites, we will show a plot of the radius rx of the sphere that contains different fractions (xpercent) of the total mass in the form of mini-satellites. We first proceed by specifying units and definitions for the relevant quantities.

   
4.1 Units

Unless specified otherwise, the units adopted for mass, length, and time are the mass of the galaxy (M), the radius of the galaxy (R), and the natural crossing time ( $t_{\rm cr} = G M^{5/2}(2 K_{\rm gal})^{-3/2}$ evaluated at the beginning of the simulation; here $K_{\rm gal}$ is the total kinetic energy associated with the stars in the galaxy). For the galaxy model considered in this paper, the revolution period relative to a circular orbit at the periphery (at R) is $\approx $ $11 t_{\rm cr}$. To return to physical units, we may refer to the case where $M = 2 \times 10^{11}~M_{\odot}$, R = 10 kpc, so that $t_{\rm cr} = 0.18 \times 10^8$ yrs and the revolution period at the periphery is $\approx $ $ 2 \times 10^8$ yrs.

   
4.2 Conservation laws

In the course of the simulation the total energy and the total angular momentum should remain constant at their initial values  $E_{\rm tot}$ and $\vec{J}_{\rm tot}$.

In the case of a single satellite we write the energy conservation as

 
$\displaystyle E_{\rm tot}= E_{\rm gal} + K_{\rm sat} + E_{\rm int}= \sum_{i=1}^... ...frac{G M_{\rm s} m}{(R_{\rm s}^2 +(\vec{r}_i-\vec{r}_{\rm s})^2)^{1/2}}\right],$     (8)

where $E_{\rm gal} = K_{\rm gal} + W_{\rm gal}$ is the total (kinetic plus gravitational) energy of the galaxy in the absence of the satellite. From the point of view of the satellite, it is natural to introduce the energy $E_{\rm sat} = K_{\rm sat} + E_{\rm int}$, so that we may split the total energy into $E_{\rm tot}= E_{\rm gal} + E_{\rm sat}$. The angular momentum conservation is given simply by $\vec{J}_{\rm tot}= \vec{J}_{\rm gal} + \vec{J}_{\rm sat} =\sum_{i=1}^{i=N} m~\vec{r}_i \times \vec{v}_i + M_{\rm s} ~ \vec{r}_{\rm s}\times \vec{v}_{\rm s}$.

In the case of several mini-satellites, following the model outlined in Sect. 3, the energy takes the form

 
                               $\displaystyle E_{\rm tot}= E_{\rm gal}$ + $\displaystyle (K_{\rm fra} + W_{\rm fra}) + E_{\rm int}= \displaystyle\sum_{i=1}^{i=N} \frac{m}{2} \left[ v_i^2 +\Phi(\vec{r}_i) \right]$  
    $\displaystyle + \displaystyle\sum_{j=1}^{j=N_{\rm f}} \frac{M_{\rm f}}{2}\left[... ...{\rm f}}{(R_{\rm f}^2 +(\vec{r}_{{\rm f}j}-\vec{r}_{{\rm f}k})^2)^{1/2}}\right]$  
    $\displaystyle - \displaystyle\sum_{i=1}^{i=N} \displaystyle\sum_{j=1}^{j=N_{\rm f}} \frac{G m M_{\rm f}}{(R_{\rm f}^2 +(\vec{r}_i-\vec{r}_{{\rm f}j})^2)^{1/2}}~,$ (9)

while the angular momentum can be written as $\vec{J}_{\rm tot}= \vec{J}_{\rm gal} + \vec{J}_{\rm fra} = \sum_{i=1}^{i=N} m~ ... ...sum_{j=1}^{j=N_{\rm f}}M_{\rm f} ~ \vec{r}_{{\rm f}j} \times \vec{v}_{{\rm f}j}$. Note that the energy associated with the system of mini-satellites now includes a term $W_{\rm fra}$that describes the gravitational energy internal to the system of fragments. By analogy with the case of the single satellite, we may wish to refer to the energy $E_{\rm fra} = (K_{\rm fra} + W_{\rm fra}) + E_{\rm int}$ as the total energy associated with the fragments, so that $E_{\rm tot} = E_{\rm gal} + E_{\rm fra}$.

   
4.3 Density and pressure distributions

In order to describe the galaxy density perturbation, we refer to the density response defined as $\rho_1 (t,\vec{r}) = \rho(t,\vec{r}-\vec{r}_0(t))-\rho_0(\vec{r})$; here $\rho_0(\vec{r})$ is the adopted density distribution of the galaxy in the absence of satellites, $\vec{r_0}(t)$ is the position of the center of mass of the galaxy in an inertial frame of reference. To diagnose the structure of the response it will be useful to consider the characteristics of the first multipoles.

In the simulations the density on the spherical mesh is expanded in associated Legendre functions $P_l^m(\cos\theta)$ up to l=6

 
                           $\displaystyle \rho(r,\theta,\varphi)$ = $\displaystyle \displaystyle\sum_{l=0}^{l=6} {\bigg[}A_{l,0} (r)~ Y_l^0(\cos\theta)$  
    $\displaystyle + \displaystyle\sum_{m=1}^{m=l} \sqrt{2} A_{l,m}(r)~ {\rm Re} (Y_l^m(\cos\theta,\varphi))$  
    $\displaystyle + \displaystyle\sum_{m=1}^{m=l} \sqrt{2} B_{l,m}(r)~ {\rm Im} (Y_l^m(\cos\theta,\varphi)){\bigg]}~,$ (10)

where the coefficients of the expansion into real spherical harmonics are given by suitable integrations over the relevant solid angle

 \begin{displaymath} A_{l,m}(r) = \sqrt{2} \int \rho(r,\theta,\varphi) {\rm Re}(Y_l^m(\cos\theta,\varphi)) {\rm d}\Omega~, \end{displaymath} (11)


 \begin{displaymath} A_{l,0}(r) = \int\rho(r,\theta,\varphi) Y_l^0(\cos\theta) {\rm d}\Omega~, \end{displaymath} (12)


 \begin{displaymath} B_{l,m}(r) = \sqrt{2} \int \rho(r,\theta,\varphi) {\rm Im}(Y_l^m(\cos\theta,\varphi)) {\rm d}\Omega~. \end{displaymath} (13)

In particular, the dipole (l=1) component of the density response $\rho_1$ is given by
 
$\displaystyle \rho_{l=1}(r,\theta,\varphi)=A_{1,0} \sqrt{3/4\pi}\cos\theta - A_... ...qrt{3/4\pi}\sin\theta\cos\varphi - B_{1,1} \sqrt{3/4\pi}\sin\theta\sin\varphi~.$     (14)

The galaxy pressure tensor is defined in terms of the star distribution function in the standard way. In the case of a single satellite sinking into a galaxy along the equatorial plane we have $\langle v_r \rangle \simeq 0$ and $\langle v_{\theta} \rangle \simeq 0.$

   
4.4 Tests

We have performed many preliminary experiments aimed at testing the numerical code.

Table 1: Single satellite, $r_{\rm s} (0)= R$.

The numerical implementation of the polytropic equilibrium configuration of Sect. 3 with $N=5~ \times~ 10^3~ \div~ 5~ \times~ 10^5$ particles has been tested by generating an unperturbed galaxy that has remained quiescent for about 50 crossing times. The radii of the spheres containing 0.5%, 1%, 1.5%, ..., 5% of the total mass have been checked to remain constant in time to within 1%.

Then, we have tested that a satellite with $R_{\rm s} = 0.1~R$ and very small mass, $M_{\rm s} = 10^{-8}~M$, behaves as a light "test particle" without disturbing the main galaxy. This satellite, placed initially at the periphery of the galaxy $r_{\rm s} = R$ on an orbit with $v_{{\rm s}\varphi} = 0.536~R/t_{\rm cr}$ and $v_{\rm sr} = 0.954 \times 10^{-3}~R/t_{\rm cr}$, remains very close to its initial circular path, never exceeding an eccentricity of e=0.05. At the end of the run, at time $t \sim 50~t_{\rm cr}$ the orbital energy of this light satellite is conserved to within 1%.

   
4.4.1 Energy and momentum conservation

In addition to the specific tests mentioned above, we have checked that globally energy and angular momentum are conserved, at $t \sim 50~t_{\rm cr}$, to better than 1%.

   
4.4.2 Comparison with earlier work by Bontekoe & van Albada (1987)

The main results of the simulations of the sinking of a single satellite by Bontekoe & van Albada (1987) have been checked in detail, with the use of simulations with N up to $5 \times 10^5$(see Sect. 5.1 below). This adds much confidence, in relation both to the overall structure of the code used in this paper and to the modeling of dynamical friction used in these investigations (since results have been checked to be largely independent of N over a range significantly larger than that available earlier).

   
5 Results of the simulations

We have performed a number of runs, mostly based on $N = 5 \times 10^5$. Runs labeled by A refer to experiments with a single sinking satellite, initially located on a circular orbit at $r_{\rm s} = R$; various combinations of satellite mass $M_{\rm s}/M$ and satellite size  $R_{\rm s}/R$ have been considered. For these runs (some are listed in Table 1), we have recorded the time $t_{\rm fall}/t_{\rm cr}$ at which the satellite reaches the central regions; for some experiments (in particular, for those labeled A8 and A9) the satellite is unable to reach the center and we have thus studied the radius of minimum radial approach $r_{\rm min}$.

Table 2: Simulations with shells of fragments.


  \begin{figure} \par\epsfig{file=fig1.eps,height=5.9cm,angle=0}\end{figure} Figure 1: Numerical measurement of the Coulomb logarithm (left) and of the satellite energy loss by dynamical friction (right) for run A1 (characterized by $M_{\rm s}/M = 0.1$, $R_{\rm s}/R = 0.1$, $t_{\rm fall}/t_{\rm cr} = 48.125$). The quantities are plotted against a Lagrangian radial coordinate in the same format as in the paper by Bontekoe & van Albada (1987).
Open with DEXTER

Runs labeled by B are listed in Table 2 and refer to the case of the slow evolution of the system in the presence of many fragments, initially located in a single shell, following the procedure of case I described in Sect. 3.2.1. Runs of type BS are based on the smoother initialization of case II outlined in Sect. 3.2.2. Runs of type BT refer to the procedure of case III of Sect. 3.2.3. Here the fourth column records the initial radius of the shell and the fifth column lists the time at the end of the experiment; all runs are made with the mutual interaction, among fragments, on, but we have also performed a few runs for which the interaction is turned off. Runs of type C involve two shells of fragments.

   
5.1 The case of one sinking satellite

We have studied the energy transfer between the galaxy and the sinking satellite and the overall energy balance, using the total energy conservation as a diagnostic of the code accuracy. In particular we have considered the quantity $\Delta E_{\rm sat} = K_{\rm sat} + E_{\rm int} - K_{\rm sat}(t = 0) - E_{\rm int}(t = 0)$. For example, for run A1 (characterized by $M_{\rm s}/M = 0.1$, $R_{\rm s}/R = 0.1$, $t_{\rm fall}/t_{\rm cr} = 48.125$) the energy lost by the satellite is initially very small (on the order of 0.02, in the units following the conventions described in Sect. 4.1, at time $t \approx 38.5~t_{\rm cr}$), and then rapidly rises (to 0.075 at time $t \approx 49.5~t_{\rm cr}$); during the run the total energy is conserved to better than 0.002. Similarly, we have considered the angular momentum balance. Here the relatively large amount of angular momentum already associated with the galaxy right after the beginning of the simulation ( $J_{\rm gal} \approx 0.0054$ at $t \approx 5.5~t_{\rm cr}$) is not yet due to the scattering of star orbits by the sinking satellite, but is mostly related to the overall orbital motion of the galaxy in the inertial frame of reference of the center of mass; eventually, the satellite is dragged in, down to the center, and loses completely its angular momentum, which is acquired by the galaxy.


  \begin{figure} { \epsfig{file=fig2a.eps,height=5.55cm,angle=0}\hspace*{1mm} \eps... ...{1mm} \epsfig{file=fig2f.eps,height=5.55cm,angle=0}\vspace*{1mm} } \end{figure} Figure 2: Six frames describing the equatorial density response $\rho _1 = \rho - \rho _0$ of the galaxy to the infalling satellite for run A1. The density wakes are portrayed at time $t/t_{\rm cr} = {\bf a)}$ 38.5, b) 39.875, c) 41.25, d) 42.625, e) 44, and f) 53.625. The white circle denotes the location of the satellite.
Open with DEXTER

During the process, following the steps indicated by Bontekoe & van Albada (1987), we can "measure" the Coulomb logarithm. The measurement of $\ln{\lambda}$ consists of the following steps: (i) Determine the velocity $\hat{v}_{\rm s}$ of the satellite with respect to the local streaming motion of the stars (the streaming motion is determined as a mean velocity of the stars in the "vicinity'' of the satellite); (ii) Measure the energy $E_{\rm sat}$ of the satellite at many instants, then determine ${\rm d}E_{\rm sat}/{\rm d}t;$ (iii) Determine the density $\rho_G$ of the galaxy at the position of the satellite by linear interpolation of densities between the 8 nearest mesh corners; (iv) Calculate $F(\hat{v}_{\rm s})$, the fraction of particles with velocity smaller than $\hat{v}_{\rm s}$. Finally the Coulomb logarithm is estimated from the relation $ 4\pi G^2 \ln\lambda = - v_{\rm s} ({\rm d}E_{\rm sat}/{\rm d} t)/(M_{\rm s}^2\rho_{\rm G} F(\hat{v}_{\rm s}))$.


  \begin{figure} { \epsfig{file=fig3a.eps,height=5.1cm,angle=0}\epsfig{file=fig3b.... ...height=5.1cm,angle=0}\epsfig{file=fig3f.eps,height=5.1cm,angle=0} } \end{figure} Figure 3: For each frame of Fig. 2 we illustrate in detail the monopole (l=0) and the dipole (l=1) contributions to the density wake in the galaxy (the various curves represent the coefficients A10, A11, A00, and B11 defined in the text; see Sect. 4.3).
Open with DEXTER

The main results that we have obtained by examining the experiments for the case of a single sinking satellite are the following:

   
5.2 The case of many fragments

   
5.2.1 Case I
For various runs of type B and C we have checked that the total energy conservation generally holds, over a period of $\approx $ $100~t_{\rm cr}$, to better than one part over 500. Yet, it is difficult to estimate the impact of dynamical friction on the galaxy evolution, because secular effects are masked by effects related to the initial lack of equilibrium. From this class of numerical experiments we have obtained the following results:

In conclusion, simulations of case I turn out to be of little direct use for the main objectives of our project. However, they are extremely useful and instructive, because they demonstrate that, in order to properly single out slow relaxation processes, one has to be very careful about the choice of initial conditions and, in particular, because they provide concrete examples of some misleading effects that can occur as a result of the use of improper models.


  \begin{figure} { \epsfig{file=fig6.eps,width=8.8cm,angle=0} } \end{figure} Figure 6: The curves illustrate the orbital decay process for the satellite for three different runs (A6, A8, and A9, labeled as 1, 2, and 3 respectively) characterized by different spatial size of the infalling satellite ( $R_{\rm s} = 0.1,~0.05,~0.025~R$; $M_{\rm s}/M = 0.05$).
Open with DEXTER

   
5.2.2 Case II

The use of the smoother initialization described as case II in Sect. 3.2.2 leads to runs where the secular effects associated with the granularity of the fragment population are more convincingly identified. The orbital decay of the fragments is illustrated in Fig. 8, where we also show for comparison a parallel run where the shell is treated as collisionless. The interaction which generates dynamical friction of the galaxy on the fragments is responsible for the slow, but systematic, development of a tangential bias in the pressure tensor, as illustrated in Fig. 9 and for an overall change in the density profile of the galaxy (see Fig. 10).


  \begin{figure} { \epsfig{file=fig7a.eps,height=6cm,angle=0}\psfig{file=fig7b.eps,height=6cm,angle=0} } \end{figure} Figure 7: Left: The density distribution change induced in the galaxy by the infall of the satellite (run A1). The thick upper curve (1) represents the initial distribution, while the thin lower curve (2) is associated with the final state. Right: Orbital decay for run A1 (cf. Fig. 2) as in Fig. 6.
Open with DEXTER


  \begin{figure} \mbox{ \epsfig{file=fig13_a.eps,height=5.3cm,angle=0}\hspace*{5mm... ...mm}\epsfig{file=fig13_d.eps,height=5.2cm,angle=0}\hspace*{5.5mm} } \end{figure} Figure 8: The fall of a shell of fragments towards the galaxy center, for the case of $N_{\rm f } = 20$ (left frame on the top, run BS1) and for $N_{\rm f} = 100$ (right frame on the top, run BS2). For comparison, the left frame at the bottom describes the time evolution of a collisionless shell (a run with $N_{\rm f} = 5 \times 10^4$). Displayed lines represent 5, 15, 25, 35, 45, 55, 65, 75, 85, and 95% of the shell mass. The right frame at the bottom displays, for run BS2, the correlation between single particle energy at time $t = 14~ t_{\rm cr}$ and energy at time $t = 118 ~t_{\rm cr}$, showing the effects of scattering associated with dynamical friction; units for energy are those indicated in Sect. 4.1.
Open with DEXTER

The check against the collisionless run is worth a special digression. In order to diagnose the collisionality present in our system better, we have performed the following test. We have studied the correlation between single particle energy (for the particles simulating the galaxy) at two different times. If the simulation is truly collisionless, the single particle energy $E_{\rm p} = m v^2/2 + m \Phi(r)$, where $\Phi(r)$ is the mean field potential, should be strictly conserved. Indeed, by means of a plot similar to the lower right frame of Fig. 8, we have checked that for the run illustrated by the lower left frame of Fig. 8 more than 90% of the particles have an energy correlation such that $\vert[E_{\rm p}(t=118~t_{\rm cr})-E_{\rm p}(t=14~t_{\rm cr})]/E_{\rm p}(t=14~t_{\rm cr})\vert< 0.06$. In turn, for run BS2, otherwise similar, but characterized by the presence of 100 fragments, a similar level of correlation for single particle energy (at 6%) is preserved within the same time interval only by 50% of the particles; the lower right frame of Fig. 8 thus illustrates the "damage" associated with dynamical friction.

   
5.2.3 Case III
The even smoother procedure of case III leads to slower orbital decays, as illustrated in the lower frames of Fig. 11. For this type of simulations, the initial distribution of fragment orbits replicates the isotropic star distribution of the underlying galaxy. Because of dynamical friction, the fragment orbits are subject to a slow process of circularization, as illustrated in the upper frames of Fig. 11.


  \begin{figure} \mbox{ \epsfig{file=fig14_a.eps,height=5.7cm,angle=0}\epsfig{file=fig14_b.eps,height=5.7cm,angle=0} } \end{figure} Figure 9: The development of pressure anisotropy in the galaxy as a result of the interaction with a shell of fragments dragged in towards the galaxy center, for the case of $N_{\rm f } = 20$ (left frame, run BS1) and for $N_{\rm f} = 100$ (right frame, run BS2). The dashed curves represent the evolving value of KT/2, where $K_{\rm T}$ is the total kinetic energy associated with the star motions in the tangential directions; the solid curve represents the evolving value of $K_{\rm r}$, the total kinetic energy associated with the star motions in the radial direction.
Open with DEXTER


  \begin{figure} { \epsfig{file=fig15.eps,height=7cm,angle=0} }\end{figure} Figure 10: The change in density profile for the galaxy as a result of the interaction with a shell of fragments dragged in towards the galaxy center, for the case of $N_{\rm f } = 20$ (run BS1). The thicker curve (1) represents the initial density profile, the thinner (2) gives the final profile.
Open with DEXTER

For each of the fragments that are involved in the simulation described as case II, we can repeat the analysis that has led to the measurement of the Coulomb logarithm from the study of the decay of the orbit of a single satellite. For example, for run BS2 we have 100 independent such determinations. In Fig. 12 we record, for run BS2 a measurement of an average value of the Coulomb logarithm. We note that the fact that the value of the Coulomb logarithm thus determined remains of the same order of magnitude as that obtained in Fig. 1 proves that the classical expression for the effects of dynamical friction captures the correct scaling with respect to the mass of the object subject to friction.

   
6 Discussion and conclusions

In this paper we have set up a numerical laboratory for the study of dynamical friction in inhomogeneous quasi-spherical galaxies and made several experiments that have elucidated the physical mechanisms that participate in the process.

The general character of the simulations presented here is distinctly different from that adopted in previous studies in this general research area, because we have addressed a physical scenario in which the slow collapse induced by dynamical friction approximately preserves the spherical symmetry present initially. This configuration is inspired by the astrophysical problem of studying the evolution of a spherical galaxy in the presence of a substantial globular cluster system or as a result of the capture of a large set of mini-satellites dragged in along random directions.

From the theoretical point of view, the framework adopted here has an important advantage with respect to the traditional case of the study of the decay of the orbit of a single satellite, because it allows us to run smoother simulations that are basically free from other effects unrelated to dynamical friction, such as those associated with lack of equilibrium in the initial configuration.

In the past, most of the attention has been focused on the effect of the host galaxy on the object subject to dynamical friction. Here, we have been able not only to test the adequacy of the classical formulae of dynamical friction (derived for ideal homogeneous models) in the inhomogeneous galaxy environment, but also to measure the evolution of the density and pressure tensor profiles in the galaxy induced by the presence of friction processes on a minority population of heavier particles. We have made quantitative tests that convince us that we are detecting genuine effects associated with secular evolution.

   
6.1 Plans for future investigations

Starting from the encouraging quantitative results obtained in this paper, there are now two major issues that we would like to address in investigations planned for the near future, which appear to be approachable by the use of the numerical laboratory developed here:


  \begin{figure} {\hspace*{7mm}\epsfig{file=fig16_a.eps,height=5cm,angle=0}\epsfig... ...e=0}\epsfig{file=fig16_d.eps,height=55mm,width=77mm,angle=0} }\par\end{figure} Figure 11: The effects of dynamical friction on a shell of fragments created out of the initial distribution function (case III of Sect. 3). The lower frames show the general decay of orbits of the fragments for the case $N_{\rm f } = 20$ (run BT3) and $N_{\rm f} = 100$ (run BT4); displayed lines represent 5, 15, 25, 35, 45, 55, 65, 75, 85, and 95% of the shell mass. For the two runs, the upper frames show the corresponding effect of circularization of orbits; here the histograms record the number of fragments (crosses represent initial conditions, open circles the conditions at the end of the run) with orbit aspect ratio above $R_{\rm min}/R_{\rm max}$, where $R_{\rm min}$ and $R_{\rm max}$ are consecutive minimum and maximum radii attained by a fragment along its orbit.
Open with DEXTER


  \begin{figure} {\epsfig{file=fig17a.eps,height=5.2cm,angle=0}\epsfig{file=fig17b.eps,height=5.2cm,angle=0} } \end{figure} Figure 12: The energies associated with the system of minisatellites (left frame) and measurement of the Coulomb logarithm $<\ln(\lambda)>~ = \sum_{j=1}^{j=N_{\rm f}} \ln(\lambda)/N_{\rm f}$ (right frame) (run BS2).
Open with DEXTER

Acknowledgements
Part of this work was supported by the Italian MIUR. We thank Tjeerd van Albada for providing us with a copy of his code and for a number of useful conversations.

References


Copyright ESO 2003