Dynamical friction and scratches of orbiting satellite galaxies on host systems

Ogiya, Go; Burkert, Andreas

doi:10.1093/mnras/stw091

Abstract

We study the dynamical response of extended systems, hosts, to smaller systems, satellites, orbiting around the hosts using extremely high-resolution N-body simulations with up to one billion particles. This situation corresponds to minor mergers which are ubiquitous in the scenario of hierarchical structure formation in the universe. According to Chandrasekhar, satellites create density wakes along the orbit and the wakes cause a deceleration force on satellites, i.e. dynamical friction. This study proposes an analytical model to predict the dynamical response of hosts as reflected in their density distribution and finds not only traditional wakes but also mirror images of over- and underdensities centred on the host. Our controlled N-body simulations with high resolutions verify the predictions of the analytical model. We apply our analytical model to the expected dynamical response of nearby interacting galaxy pairs, the Milky Way–Large Magellanic Cloud system and the M31–M33 system.

galaxies: dwarf, galaxies: evolution, galaxies: formation, galaxies: kinematics and dynamics, Local Group, Magellanic Clouds

1 INTRODUCTION

Chandrasekhar (1943) first discussed a fundamental physical process called dynamical friction for collisionless systems. According to his calculation, a massive object moving through a sea of particles produces density enhancements, ‘wakes’ behind itself due to its gravitational force and the gravitationally induced wakes generate a decelerating force on the moving massive object. As a consequence, the massive object loses its orbital energy and angular momentum. Ostriker (1999) later on studied the process of dynamical friction for gaseous systems.

Dynamical friction arises in various astronomical phenomena (Binney & Tremaine 2008). Considering galaxy mergers, satellite galaxies orbiting around larger ones lose orbital energy and angular momentum by the effects of dynamical friction and eventually fall into the centre of their hosts (e.g. White 1976; Colpi, Mayer & Governato 1999; Taylor & Babul 2001). A prominent example is the orbit decay of the Large Magellanic Cloud (LMC) that has been studied in great details (e.g. Tremaine 1976; Murai & Fujimoto 1980; Gardiner, Sawa & Fujimoto 1994; Hashimoto, Funato & Makino 2003). Dynamical friction also plays important roles for the formation and evolution of black hole binaries (Begelman, Blandford & Rees 1980; Yu 2002; Makino & Funato 2004; Tanikawa & Umemura 2011; Fiacconi et al. 2013). Studies based on numerical simulations reveal the validity but also problems of Chandrasekhar's theory (Fujii, Funato & Makino 2006; Read et al. 2006; Inoue 2009).

Important insight has been gained, applying the arguments of dynamical friction to flattened systems, such as galactic stellar discs. Binney (1977) proposed a correction of Chandrasekhar's formulation of dynamical friction for systems with anisotropic velocity distributions. Peñarrubia, Just & Kroupa (2004) confirmed the reliability of Binney's treatment numerically. Dynamical friction generates the heating source to explain the thickness of stellar discs (e.g. Velazquez & White 1999; Mori & Rich 2008). Peñarrubia, Kroupa & Boily (2002) and Read et al. (2008) found that orbiting satellite galaxies around a larger galaxy with disc components are dragged into the disc plane by their dynamical friction. This may explain the distribution of satellite galaxies around the Milky Way (MW; Peñarrubia et al. 2002) and predicts the existence of stable disc structures of dark matter which significantly boost capturing dark matter particles (Read et al. 2008, 2009).

Gravitational wakes were investigated in details analytically by early work (e.g. Mulder 1983; Weinberg 1989). However, only a few studies have confirmed gravitational wakes directly by using numerical simulations in spite of the importance to get a better understanding of dynamical friction. Weinberg & Katz (2002, 2007) showed over- and underdensities induced by bars around the centre of galaxies. Antonini & Merritt (2012) illustrated the density response induced by black holes in galactic nuclei. For gaseous systems, Kim & Kim (2009) and Kim (2010) analysed density wakes by hydrodynamic simulations.

Previous numerical studies were limited by particle resolution. The motivation of this study is to investigate the features of scratches induced by orbiting smaller systems (satellites) on larger ones (hosts) in great details adopting one billion particles. The structure of this paper is as follows. Section 2 provides an analytical description to predict the response of the hosts as reflected in their density distribution. We find not only traditional gravitational wakes, but also a mirror image of over- and underdensities around the centre of hosts. In Section 3, we perform high-resolution N-body simulations to test the analytical prediction. The results of our simulations well match to the predictions of the analytical description. We use the analytical model for the galaxy pairs, the MW–LMC system and the M31–M33 system as applications in Section 4. Section 5 summarizes and discusses the results.

2 ANALYTICAL MODEL

The formula of dynamical friction proposed by Chandrasekhar (1943) can be written as

\begin{equation} \frac{\mathrm{d} {\boldsymbol v}_{\rm sat}}{\mathrm{d}t} = -4 \pi G^2 M_{\rm sat} \ln {\Lambda } \rho _{\rm host} \frac{{\boldsymbol v}_{\rm sat}}{v_{\rm sat}^3}, \end{equation}

(1)

where M_sat and ρ_host are the mass of the satellite and the density of the host, and |${\boldsymbol v}_{\rm sat}$| represents the velocity of the satellite in the rest frame of the host. G and ln Λ are the gravitational constant and the Coulomb logarithm, respectively. Here, the host systems are assumed to have an isotropic velocity field in the initial equilibrium state. We suppose that the mass of a particle which belongs to the host is negligible compared with M_sat and assume that all particles in the host contribute to the force of dynamical friction for simplicity although only host particles which satisfy |$v < |{\boldsymbol v}_{\rm sat}|$| cause dynamical friction in the original Chandrasekhar's theory. These assumptions do not influence our conclusion as shown in Section 3. In Chandrasekhar's derivation of dynamical friction, an infinite homogeneous particle distribution is assumed. However, galaxies and dark haloes have finite sizes and gradients in their density structures. Thus, the Coulomb logarithm may have a dependence on the satellite position within the host (e.g. Hashimoto et al. 2003; Just & Peñarrubia 2005). One may modify equation (1) as follows:

\begin{equation} \frac{\mathrm{d} {\boldsymbol v}_{\rm sat}}{\mathrm{d}t} = -4 \pi G^2 M_{\rm sat} \ln {\Lambda ({\boldsymbol r}_{\rm sat})} \rho _{\rm host}({\boldsymbol r}_{\rm sat}) \frac{{\boldsymbol v}_{\rm sat}}{v_{\rm sat}^3}, \end{equation}

(2)

where |${\boldsymbol r}_{\rm sat}$| is the position of the satellite in the frame whose origin is the centre of the host. Hereafter, we study the evolution in this rest frame.

Orbiting satellites are perturbers for hosts. Because of the gravitational force of the satellite, induced density perturbations will arise in the host. We label physical quantities of the host in the equilibrium state as ‘0’ and induced ones as ‘1’, respectively, i.e.

\begin{eqnarray} \rho _{\rm host} = \rho _0 + \rho _1, \end{eqnarray}

(3)

\begin{eqnarray} \Phi _{\rm host} = \Phi _0 + \Phi _1, \end{eqnarray}

(4)

where Φ is the gravitational potential. Poisson's equation connects the gravitational potential with the density field. For the induced quantities, we get

\begin{equation} {\bf \nabla }^2 \Phi _1 = -{\bf \nabla } \cdot {\boldsymbol g}_1 = 4 \pi G \rho _1, \end{equation}

(5)

where |${\boldsymbol g}_1$| is the specific gravitational force caused by the induced density perturbations. We can regard the deceleration force of the dynamical friction process as |${\boldsymbol g}_1$|⁠. Substituting equation (2) into equation (5), the induced density field is

\begin{eqnarray} \rho _1({\boldsymbol r}, {\boldsymbol v}_{\rm sat}) = G M_{\rm sat} [ \rho _0({\boldsymbol r}) {\bf \nabla } \ln {\Lambda ({\boldsymbol r})} + \ln {\Lambda ({\boldsymbol r})} {\bf \nabla } \rho _0({\boldsymbol r})] \cdot \frac{{\boldsymbol v}_{\rm sat}}{v_{\rm sat}^3}. \end{eqnarray}

(6)

Here, we assume that the absolute values of induced quantities are much smaller than those in the equilibrium state.

We need to provide the density distribution of the host in the equilibrium state, |$\rho _0({\boldsymbol r})$| and the Coulomb logarithm, |$\ln {\Lambda ({\boldsymbol r})}$| to use equation (6). The density distribution of the host galaxy may be expressed well by the following double power-law formula,

\begin{equation} \rho _0(r) = \frac{\rho _{\rm s}}{r^{\alpha } [1 + (r / r_{\rm s})]^{\beta }}, \end{equation}

(7)

where ρ_s and r_s are the scale density and length, respectively. The model of α = 1, β = 2 is known as the Navarro–Frenk–White (NFW) profile which well matches the density structure of cold dark matter haloes obtained in dissipationless cosmological simulations (Navarro, Frenk & White 1997). The Hernquist profile (Hernquist 1990) which reproduces the de Vaucouleurs law (de Vaucouleurs 1948) and is frequently used to express density structures of elliptical galaxies or bulges, corresponds to the model of α = 1, β = 3 (see however e.g. Burkert 1995). The derivative of equation (7) is given by

\begin{equation} {\bf \nabla } \rho _0 (r) = -\rho _0(r) \left [ \frac{\alpha }{r} + \frac{\beta }{r + r_{\rm s}} \right ] \frac{{\boldsymbol r}}{r}. \end{equation}

(8)

The Coulomb logarithm is defined as ln Λ ≡ ln (b_max/b_min) where b_max and b_min are the maximum and minimum impact parameters, respectively. The maximum impact parameter, b_max should depend on the vector pointing from the satellite to given points, |${\boldsymbol r}$|⁠, |${\boldsymbol d} = {\boldsymbol r} - {\boldsymbol r}_{\rm sat}$| since b_max determines the region affected by the gravitational force of the satellite. For simplicity, we suppose that b_min is a constant,

\begin{equation} b_{\rm min} = A l, \end{equation}

(9)

where A and l mean a constant and the size of the satellite, respectively (see also Tremaine 1976; Hashimoto et al. 2003).

We define the maximum impact parameter, b_max by the following procedure. Let us consider a plane which is perpendicular to the velocity vector of the satellite, |${\boldsymbol v}_{\rm sat}$|⁠. The vector, |${\boldsymbol v}_{\rm sat} \Delta t$| measures the distance between the satellite and plane and determines the plane on which a point is given by |${\boldsymbol r}_{\rm sat} - {\boldsymbol v}_{\rm sat} \Delta t$|⁠. Given points, |${\boldsymbol r}$| on the plane satisfy the condition,

\begin{equation} ({\boldsymbol r} - {\boldsymbol r}_{\rm sat} + {\boldsymbol v}_{\rm sat} \Delta t) \cdot {\boldsymbol v}_{\rm sat} = 0. \end{equation}

(10)

From equation (10), Δt is derived by

\begin{equation} \Delta t = -\frac{{\boldsymbol v}_{\rm sat} \cdot ({\boldsymbol r} - {\boldsymbol r}_{\rm sat})}{v_{\rm sat}^2} = -\frac{{\rm d} \cos {\phi }}{v_{\rm sat}}, \end{equation}

(11)

where ϕ is the angle between the vectors, |${\boldsymbol d}$| and |${\boldsymbol v}_{\rm sat}$|⁠. We define b_max as the length of a free-fall motion in the time interval, Δt, i.e. b_max = (1/2)aΔt², behind the satellite at which |${\boldsymbol v}_{\rm sat} \cdot {\boldsymbol d} < 0$| is satisfied. The absolute value of the gravitational acceleration in the perpendicular direction to |${\boldsymbol v}_{\rm sat}$|⁠, a is

\begin{equation} a = \frac{G M_{\rm sat}}{[b_{\rm max}^2 + (v_{\rm sat} \Delta t)^2]^{3/2}} b_{\rm max}. \end{equation}

(12)

For a simple evaluation of the gravitational acceleration, we assume that the satellite is located at |${\boldsymbol r}_{\rm sat}$| from t = T − Δt to t = T. Solving equation (12) for b_max, we obtain

\begin{equation} b_{\rm max} = \sqrt{\biggl [ \frac{G^2 M_{\rm sat}^2 \Delta t^4}{4} \biggr ]^{1/3} - (v_{\rm sat} \Delta t)^2} \end{equation}

(13)

for points at which |${\boldsymbol v}_{\rm sat} \cdot {\boldsymbol d} < 0$| is satisfied. Also,

\begin{equation} b_{\rm max} = B b_{\rm min} \end{equation}

(14)

for points at which |${\boldsymbol v}_{\rm sat} \cdot {\boldsymbol d} \ge 0$| is satisfied. Here, B is a constant. The derivative is

\begin{equation} {\bf \nabla } \ln {\Lambda ({\boldsymbol r})} = \frac{1}{b_{\rm max}^2} \biggl [ \biggl (\frac{2\,G^2 M_{\rm sat}^2 \cos ^4{\phi }}{27 v_{\rm sat}^4} d \biggr )^{1/3} - d \cos ^2{\phi } \biggr ] \frac{{\boldsymbol d}}{d} \end{equation}

(15)

for points at which |${\boldsymbol v}_{\rm sat} \cdot {\boldsymbol d} < 0$| and

\begin{equation} {\bf \nabla } \ln {\Lambda ({\boldsymbol r})} = 0 \end{equation}

(16)

for points at which |${\boldsymbol v}_{\rm sat} \cdot {\boldsymbol d} \ge 0$| is satisfied, respectively. When b_max < Bb_min in equation (13), we set b_max = Bb_min and |${\bf \nabla } \ln {\Lambda ({\boldsymbol r})} = 0$|⁠. Hence, a constant, B defines the minimum value of the Coulomb logarithm.

Combined with equations (6) and (15), the first term in equation (6) generates density perturbations around the satellite. Since equation (13) defines b_max on the back side of the satellite motion, the induced density arises only behind the satellite. Hence, the first term represents Chandrasekhar's original gravitational wake. As indicated by combining equations (6) and (8), the sign of the density fluctuations caused by the second term of equation (6) depends on the angle between the vectors of the satellite velocity and the position. For spherical systems, the density fluctuations are symmetric with respect to the centre of the host.

3 SIMULATIONS

3.1 Set up

We perform controlled collisionless N-body simulations with a parallelized code optimized for graphic processing unit (GPU) clusters. The numerical code employs the tree algorithm proposed by Barnes & Hut (1986) and the second-order Runge–Kutta scheme in time integration. Along the lines of Nakasato (2012), CPU cores construct tree structures of particles and GPU cards compute the gravitational force among particles through traversing the tree structures (Ogiya et al. 2013).

We simulate mergers between two systems, the host and the satellite. In order to generate N-body systems which follow the NFW density profile, i.e. α = 1, β = 2 in equation (7), in the equilibrium states, we use the method proposed by (Kazantzidis, Zentner & Kravtsov 2006, see also Eddington 1916). The phase–space distribution function is assumed to depend only on energy and the systems have an isotropic velocity dispersion initially. The distribution of particles is truncated at the virial radius, R_v that is related to the virial mass, M_v

\begin{equation} M_{\rm v} = \frac{4 \pi }{3} \Delta \rho _{\rm crit} (1+z)^3 R_{\rm v}^3, \end{equation}

(17)

where ρ_crit and z are the critical density of the universe and redshift, respectively. We adopt a conventional value of the overdensity, Δ = 200 in this study. The concentration parameter, c is defined by c = R_v/r_s. The host and satellite have c = 10 and 15, respectively.

The initial separation between the centres of the host and satellite is the virial radius of the host, R_v,host. The orbit has circularity η = 0.5 initially. This value is very typical, as suggested by cosmological N-body simulations (Khochfar & Burkert 2006). The initial velocity of the satellite is V_ini = ηV_c(R_v,host) = η(GM_v,host/R_v,host)^1/2. Here, V_c(r) is the circular velocity measured at r. In the coordinate system, centred on the host, the initial position and velocity vectors of the satellite are |${\boldsymbol X} = (R_{\rm v,host}, 0, 0)$| and |${\boldsymbol V} = (0, V_{\rm ini}, 0)$|⁠, respectively.

The position and velocity of the host and satellite are determined by the following procedure. We calculate the bound mass of each system with the method proposed by Funato, Makino & Ebisuzaki (1993). Particles initially belong either to the host or to the satellite system. At each snapshot, we compute the gravitational potential of each particle by particles which are still bound to the system at the previous snapshot. The bulk velocity of each system is determined by an iterative procedure (Fujii et al. 2006). When the binding energy of a particle is positive, the particle is regarded as an escaper. We define the centre of mass and bulk velocity of bound particles as the position and velocity of each system at given time.

We construct the host systems which have the virial mass, M_v,host with N_host particles. For satellite systems with a virial mass, M_v,sat, N_sat = (M_v,sat/M_v,host)N_host particles are employed. Hence, all particles have equal masses in each simulation and the total number of particles is N_tot = (1 + M_v,sat/M_v,host)N_host. Table 1 summarizes the simulations. Because we use sufficient numbers of particles, artificial effects such as two-body relaxation are negligible. In the simulations, the softening length is |$\epsilon = 4 R_{\rm v, host} / \sqrt{N_{\rm host}}$| and the tolerance parameter of the tree algorithm is θ = 0.6 (Power et al. 2003).

Table 1.

Open in new tab

Summary of simulation runs. Description of columns: (1) name of simulation runs. (2) mass ratio between the host and satellite, M ≡ M_v,host/M_v,sat. (3) Total number of particles, N_tot = N_host + N_sat.

Run	M	N_tot
(1)	(2)	(3)
A	100.0	105 906 176
B	50.0	106 954 752
C	20.0	110 100 480
D	10.0	115 343 360
E	10.0	1038 090 240

Run	M	N_tot
(1)	(2)	(3)
A	100.0	105 906 176
B	50.0	106 954 752
C	20.0	110 100 480
D	10.0	115 343 360
E	10.0	1038 090 240

Table 1.

Open in new tab

Summary of simulation runs. Description of columns: (1) name of simulation runs. (2) mass ratio between the host and satellite, M ≡ M_v,host/M_v,sat. (3) Total number of particles, N_tot = N_host + N_sat.

Run	M	N_tot
(1)	(2)	(3)
A	100.0	105 906 176
B	50.0	106 954 752
C	20.0	110 100 480
D	10.0	115 343 360
E	10.0	1038 090 240

Run	M	N_tot
(1)	(2)	(3)
A	100.0	105 906 176
B	50.0	106 954 752
C	20.0	110 100 480
D	10.0	115 343 360
E	10.0	1038 090 240

We are free to scale the mass, length and time-scales since our simulations only take into account gravitational effects. For MW-sized haloes with M_v,host = 10¹² M_⊙ and cosmic redshifts, z = 0, the virial radius of the host halo is R_{v, host} ≈ 211kpc and the dynamical time of the host, T_d,host defined by

\begin{equation} T_{\rm d,host} \equiv \frac{R_{\rm v,host}}{V_{\rm c}(R_{\rm v,host})} = \sqrt{\frac{R_{\rm v,host}^3}{G M_{\rm v,host}}}, \end{equation}

(18)

is ≈1.45Gyr. Here, a Hubble constant of H₀ = 67.5 km s⁻¹ Mpc⁻¹ (Planck Collaboration XIII 2015) is adopted. Simulation data are output every 0.05T_d,host.

3.2 Results

Fig. 1 summarizes a typical case of galaxy merger. The satellite loses angular momentum due to dynamical friction and the orbit shrinks gradually (upper panel). In the meantime, the satellite is stripped by the tidal force of the host, especially when it approaches the pericentres. Eventually, it is completely destroyed at t = 5.65T_d,host in run D and 6.0T_d,host in run E, respectively (lower panel). Fig. 1 also tests the numerical convergence of the simulations. The orbital evolution is well converged. The evolution of the bound mass of the satellite for both resolutions deviates only at t > 4T_d,host, at which most of its mass has already been stripped away and significant density scratches do not arise (see Figs 2 and 3). Hence, the results of this paper do not depend on the number of particles, N. As shown in this figure, dynamical friction plays a key role during the merging process.

Figure 1.

Time evolution of the satellite orbit (upper panel) and bound mass of the satellite, M_b,sat (lower panel) derived from run D (black) and E (red). In the upper panel, solid and dashed lines represent the specific angular momentum, L of the satellite orbit relative to the host and the distance, D between the centres of the host and satellite. L and D are normalized by the initial values, L₀ and D₀. M_b,sat and time, t are scaled by the virial mass of the host, M_v,host and dynamical time, T_d,host, respectively. [A color version of this figure is available in the online version.]

Open in new tab Download slide

Figure 2.

Distribution of enhancement and reduction in the column density of the host system derived from run E. The colour bar represents enhancement and reduction in the column density, (Σ − Σ₀)/Σ₀, where Σ₀ is the initial column density at given position in the host frame. Spatial coordinates are scaled by the scalelength of the host system, r_s. Thick and thin black lines show the satellite orbit and motion of the minimum potential point in the centre-of-mass frame. Panels (a), (b), (c), (d), (e) and (f) demonstrate snapshots at T = 1.15, 1.25, 1.35, 1.5, 2.0 and 2.5T_d,host, respectively. [A color version of this figure is available in the online version.]

Open in new tab Download slide

Figure 3.

Same as Fig. 2. Panels (g), (h), (i), (j), (k) and (l) demonstrate snapshots at T = 3.0, 3.5, 4.0, 5.0, 6.0 and 7.0T_d,host, respectively. The thin black line which represents the motion of the minimum potential point in the centre-of-mass frame is not drawn for better visibility. [A color version of this figure is available in the online version.]

Open in new tab Download slide

Figs 2 and 3 illustrate that the orbiting satellite leaves significant scratches on the host. We find two kinds of scratches in its density distribution. The first one is the gravitational wake which has been found and discussed by many previous studies. The gravitational wake arises along the satellite orbit as imagined by Chandrasekhar (1943). It can be found more clearly in the early phase of the merger process (see panels a and b) since it mixes with another type of scratch described below in the later phase.

The second type of scratch caused by the gravitational field of the satellite is a pair of density enhancements and reductions around the centre of the host. Similar results have been obtained analytically by Weinberg (1989). This distribution of over- and underdensities is mirror image-like and the mirror plane is roughly perpendicular to the direction of the velocity vector of the satellite. The underdensity is located in the direction of the velocity vector of the satellite and the overdensity arises in the opposite direction. The directions are more visible in the early phase again due to the mixing of both, the wake and the central perturbation in the later phase.

This central dipole scratch is caused by the motion of the location of the minimum of the potential, i.e. highest density point in the initial state. The tidal force of the satellite halo perturbs the position of the highest density point and it well matches the position of the highest overdensity. Hence, the motion of the minimum potential point well traces the overdensity around the centre of the host in the centre-of-mass frame (see thin black line in Fig. 2). Assuming that the satellite is a less massive particle in a two-body problem, the thin line looks like an orbit of the corresponding more massive particle. Because of mass conservation, the overdensity on one side causes a corresponding underdensity on the opposite side of the minimum potential point with respect to the centre of mass in the host. Since there is a single point where the potential has its minimum in the host, the central over- and underdensities have a dipole structure. The effect does not affect the bulk structure of the host and it retains the initial spherical shape on the whole.

Figs 2 and 3 also show that the amplitude of over- and underdensities decreases with time. This is because of the decreasing satellite mass as a result of tidal stripping (see equation 6). After tidal disruption of the satellite, little scratches remain for some time (panels k and l).

We study the relation between the amplitude of the scratches and the satellite mass. Fig. 4 represents the maximum value of enhancement in column density as a function of the initial satellite mass, M_v,sat and shows that the maximum amplitude is proportional to M_v,sat. The maximum enhancements are obtained when the satellite approaches the first pericentre (t ∼ 1.5T_d,host) in all simulation runs. At that time, the gravitational wake merges with the central overdensity of the host system. In the analytical model, the satellite is regarded as a point mass and the amplitude in the density perturbation is proportional to the satellite mass at given points. As a consequence, the amplitude in the column density perturbation should also be proportional to the satellite mass. The results of our simulations validate the assumption in the analytical model. The comparison of run D and E (see Fig. 4) shows that the numerical simulations are well converged.

Figure 4.

Maximum values of enhancements in the column density. Red, green and blue lines show the results for XY, YZ and ZX projection planes, respectively. The dashed line represents the scaling, ∝ M_v,sat which is expected from equation (6). Crosses correspond to runs A, B, C, D and the open circle is run E. [A color version of this figure is available in the online version.]

Open in new tab Download slide

3.3 Comparison with analytical predictions

In order to test the validity of the analytical model described in Section 2 and to understand the simulation results, we compare our analytical model predictions with results of the simulation. A spherical system which follows an NFW density profile is assumed as the initial unperturbed state. This corresponds to assuming that the centre of the system is the centre of mass. Fig. 5 demonstrates the predicted enhancement and reduction in the distribution of column density of the host system. Comparing the simulation results, panel (a) in Fig. 2, with the top panel in Fig. 5, the analytical prediction well reproduces the results of the simulation, not only the distribution of enhancement and reduction but also the amplitude. Middle and bottom panels show the contributions from the first and second terms in equation (6). The results clearly indicate that the first and second terms make the gravitational wake and the mirror image of the over- and underdensities, respectively. Because most of previous studies have assumed homogeneous background density, the effects of the second term, the mirror image of the over- and underdensities have not been found and discussed. However, the feature should arise in many astrophysical systems such as galaxies and galaxy clusters since their density distributions have gradients.

Figure 5.

Distribution of enhancement and reduction in the column density derived by the analytical model. The position, velocity and mass of the satellite are taken from the snapshot at t = 1.15T_d,host in run E. The top panel demonstrates the total enhancement and reduction. Middle and bottom panels show the contribution from the first and second terms in equation (6), respectively. Constant numbers, A = 3.0 and B = 1.5 are adopted and the size of the satellite is set to be l = r_s,sat = R_v,sat/15. [A color version of this figure is available in the online version.]

Open in new tab Download slide

As described above, the predictions well match the simulation results and the analytical model provides a clear understanding. However, there still remain small deviations between the simulation results and the analytical predictions. The direction of the mirror plane is slightly different. This is mainly due to two effects. First of all, the stripped mass is not considered in the analytical model. The stripped mass is distributed along the satellite orbit in the simulation. On the other hand, the analytical model treats the satellite as a point mass and does not consider the stripped mass. This effect should become more important in the last phase of the merger process since more satellite mass has been stripped. The second effect is changes in the density distribution of the host. In the analytical model, we assume the initial density distribution of the host as the background field (physical quantities labelled ‘0’ in Section 2). The density distribution of the host however changes with time. Hence, it is sensible to avoid applying the analytical model to systems in violently changing dynamical states.

4 APPLICATION

In this section, we apply the analytical model to nearby interacting galaxy pairs, the MW–LMC system and the M31–M33 system. Information about the position, velocity and mass of the satellites, the LMC and M33 and the background density field of the hosts, the MW and M31 are needed. We assume an NFW model with a virial mass, M_v = 1.26 × 10¹² M_⊙ and concentration parameter c = 10 for the unperturbed density field of the MW and M31 (van der Marel et al. 2012). Actually, the observed density field of the hosts has been already perturbed, i.e. ρ_obs = ρ_host = ρ₀ + ρ₁. We assume here that |ρ₁| is much smaller than ρ₀. The distance between the Solar system and the Galactic Center is assumed to be R_sol = 8.5 kpc (Kerr & Lynden-Bell 1986). The parameters in equations (9) and (14) are the same as those used to plot Fig. 5, A = 3.0 and B = 1.5.

We take into account only a dark matter halo in the analysis for simplicity. Baryon components of host galaxies, such as bulges, discs and stellar haloes also react to the gravitational force of satellite galaxies and density scratches arise in them. Ongoing observations, e.g. Gaia and Subaru Hyper Suprime-Cam., may find not only density fluctuations, but also fluctuations in the velocity field caused by the induced density fields. Combining observational data with our analytical model might be interesting in order to constrain the orbits and masses of the satellite galaxies.

4.1 MW–LMC

It is useful to adopt a Cartesian coordinate system (X, Y, Z), the so-called Galactocentric rest frame (e.g. Gardiner et al. 1994). In this coordinate system, the origin corresponds to the Galactic Center and the X-, Y- and Z-axes point in the direction from the Solar system to the Galactic Center, in the direction of the Galactic Rotation of the Solar system and towards the Galactic North Pole, respectively. The position of the Solar system is given by |${\boldsymbol R}_{\rm sol} = (-R_{\rm sol}, 0, 0)$|⁠. van der Marel et al. (2002) provide the position of the LMC,

\begin{equation} {\boldsymbol r}_{\rm LMC} = (-0.78, -41.55, -26.95) \ {\rm kpc}, \end{equation}

(19)

and the relative velocity of the LMC with respect to the Galactic Center is obtained by Kallivayalil et al. (2013),

\begin{equation} {\boldsymbol v}_{\rm LMC} = (-57 \pm 13, -226 \pm 15, 221 \pm 19) \ {\rm km \ s^{-1}}. \end{equation}

(20)

The total dynamical mass of the LMC is uncertain by a factor of 10. The enclosed mass within 8.7 kpc from the centre of the LMC is (1.7 ± 0.7) × 10¹⁰ M_⊙ (van der Marel & Kallivayalil 2014). The total mass should be greater than this value. Determining the total mass of the LMC by using the abundance matching technique (Guo et al. 2010), the upper mass limit of the LMC is 2.5 × 10¹¹ M_⊙ (Kallivayalil et al. 2013). This is consistent with the estimation by Peñarrubia et al. (2016). We assume that the mass and size of the LMC are M_LMC = 10¹¹ M_⊙ and l = 8.7 kpc, respectively. As shown in Fig. 4 , the amplitude of the density enhancement can be scaled by ∝ M_LMC.

Fig. 6 demonstrates the predicted enhancement and reduction in the column density distribution of the MW. XY and ZX planes in the Galactocentric coordinate are good to find clear density scratches of the LMC. When one sees the south-side sky, the column density in the direction of the Galactic Rotation of the Solar system (plus Y) is expected to be systematically greater than that in the opposite direction (upper panel). Also, the column density on the side of the Galactic North Pole (plus Z) should be systematically lower than that on the opposite side when one looks into the opposite direction of the Galactic Rotation of the Solar system (lower panel).

Figure 6.

Distribution of enhancement and reduction in the column density of the MW induced by the Large Magellanic Cloud. The colour bar represents the enhancement and reduction in the column density, (Σ − Σ₀)/Σ₀, where Σ₀ is the column density of the background field at given position in the Galactocentric coordinate. Top, middle and bottom panels show the results for XY, YZ and ZX planes, respectively. In each panel, the column density is derived by integration of the density field in the range of Z = [−200, 0] kpc, X = [−200, −R_sol] kpc and Y = [−200, 0] kpc. [A color version of this figure is available in the online version.]

Open in new tab Download slide

4.2 M31–M33

The position and velocity vectors of M31 and M33 are obtained by van der Marel et al. (2012). We adopt a coordinate system, (X_M31, Y_M31, Z_M31) in which the origin is the centre of M31 and the Z_M31-axis points in the direction from the Solar system to the centre of M31. X_M31- and Y_M31-axes are perpendicular to the Z_M31-axis and a right-handed system is constructed. In the M31 rest frame, the position and velocity of M33 are

\begin{eqnarray} {\boldsymbol r}_{\rm M33} &=& (140.5, 146.1, -2.3) \ {\rm kpc}, \nonumber \\ {\boldsymbol v}_{\rm M33} &=& (-147.4, -72.2, 117.9) \ {\rm km \ s^{-1}}. \end{eqnarray}

(21)

The total mass of M33 is uncertain by a factor of 10 similar to the LMC mass. Corbelli (2003) found that dark halo mass out to 17 kpc from the centre of M33 is ∼5 × 10¹⁰ M_⊙. Seigar (2011) obtained a virial mass of the dark halo surrounding M33 of (2.2 ± 0.1) × 10¹¹ M_⊙ from the H i rotation curve. We assume that the mass and size of M33 are M_M33 = 10¹¹ M_⊙ and l = 17 kpc, respectively.

The analytical model predicts a mirror image of the density enhancement and reduction around the centre of M31 as shown in Fig. 7. The result appears to violate the basic assumption, ρ₀ ≫ |ρ₁|, but the amplitude of density enhancement can be scaled by ∝ M_M33.

Figure 7.

Distribution of enhancement and reduction in the column density of M31 induced by M33. The colour bar represents the enhancement and reduction in the column density, (Σ − Σ₀)/Σ₀, where Σ₀ is the column density of the background field at given position in the M31 rest frame. The column density is derived by integration in the range of Z_M31 = [−200, 200] kpc. [A color version of this figure is available in the online version.]

Open in new tab Download slide

5 SUMMARY AND DISCUSSION

We have investigated the dynamical response of extended host systems to the gravitational force of orbiting satellite systems, ‘scratches’. The scratches are classified into two types: the first one is the gravitational wake along the orbit of satellites as discussed in Chandrasekhar (1943). The second type is a mirror image of the over- and underdensities which become more evident in the centre of the hosts. The mirror plane is perpendicular to the direction of the satellite velocity. We derive features analytically from Chandrasekhar's formula of dynamical friction. Our N-body simulations validate the analytical predictions well.

The scratches may be found in nearby interacting galaxies by observations. The dynamical mass including a dark halo of the satellite galaxies, the LMC and M33 is still uncertain by a factor of 10 (e.g. Corbelli 2003; Seigar 2011; Kallivayalil et al. 2013; van der Marel & Kallivayalil 2014). As indicated by equation (6) and shown in Fig. 4, the amplitude of the induced density is proportional to the satellite mass. Combining the analytical model with observations, new constraints for the satellite masses may be provided.

The form of the Coulomb logarithm is important in order to determine the features and amplitudes of scratches. In this paper, we adopt a simple formula to provide the Coulomb logarithm as a function of position. A caveat is the constant minimum impact parameter, b_min in equation (9) and the parameter, B in equation (14). We determine them by fitting analytical predictions to the simulation result but they may vary from system to system. Actually, b_min may depend on the local density since it should have similar values as the typical distance between the satellite and nearby particles. More systematic studies can help to improve the form of the Coulomb logarithm.

In a later step, it also would be important to consider more realistic configurations of host systems with many satellite systems orbiting around them. The analytical arguments in this paper might also help to understand the dynamical phenomena in these more complexed systems.

We are grateful to the anonymous referee for providing many helpful comments and suggestions. We thank Alessandro Ballone, Manuel Behrendt, Masashi Chiba, Jorge Cuadra, Guinevere Kauffmann, Takanobu Kirihara, Lucio Mayer, Yohei Miki, Masao Mori, Daisuke Nagai, Simon White and Kohji Yoshikawa for fruitful discussions. Numerical simulations were performed with HA-PACS at the Center for Computational Sciences at University of Tsukuba. This work was supported by Grant-in-Aid for JSPS Fellows (25-1455 GO) and the DFG cluster of excellence ‘Origin and Structure of the Universe’ (www.universe-cluster.de).

†

Max–Planck Fellow.

REFERENCES

Antonini

F.

Merritt

D.

2012

ApJ

745

83

Month:	Total Views:
November 2016	1
February 2017	1
March 2017	1
May 2017	3
June 2017	6
July 2017	2
August 2017	4
September 2017	1
October 2017	1
November 2017	3
December 2017	1
January 2018	2
February 2018	1
March 2018	1
April 2018	2
May 2018	5
August 2018	1
September 2018	7
October 2018	2
November 2018	2
January 2019	3
February 2019	4
April 2019	4
May 2019	8
June 2019	1
July 2019	10
August 2019	5
September 2019	15
October 2019	14
November 2019	5
December 2019	11
January 2020	10
February 2020	9
March 2020	10
April 2020	5
May 2020	3
June 2020	4
July 2020	3
August 2020	9
September 2020	12

Article Contents

Dynamical friction and scratches of orbiting satellite galaxies on host systems Free

Abstract

1 INTRODUCTION

2 ANALYTICAL MODEL

3 SIMULATIONS

3.1 Set up

3.2 Results

3.3 Comparison with analytical predictions

4 APPLICATION

4.1 MW–LMC

4.2 M31–M33

5 SUMMARY AND DISCUSSION

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Dynamical friction and scratches of orbiting satellite galaxies on host systems