Dynamical Friction and Tidal Interactions

When a body of mass M moves with velocity v _M through an infinite homogeneous medium of stars with average stellar mass m, m ≪ M, it slows down due to the gravitational drag of the stars. The resulting dynamical friction force was first calculated by Chandrasekhar (Chandrasekhar 1943; Binney & Tremaine 2008) and is given by

$$\begin{eqnarray}&&M\displaystyle \frac{d{{\boldsymbol{v}}}_{M}}{{dt}}=-4\pi {G}^{2}{M}^{2}\displaystyle \frac{{{\boldsymbol{v}}}_{M}}{{v}_{{\rm{M}}}^{3}}\rho (\lt {v}_{M})\mathrm{ln}{\rm{\Lambda }},\end{eqnarray} \tag{ 1 }$$

where v_M = ∣ v _M ∣, ρ( < v_M ) is the density of stars with speeds less than v_M and $\mathrm{ln}{\rm{\Lambda }}$ is the so-called Coulomb logarithm. In fact, ${\rm{\Lambda }}={b}_{\max }/{b}_{\min }$, where b is the scattering impact parameter, ${b}_{\max }=D$ is the diameter of the smallest sphere that completely surrounds the whole system and ${b}_{\min }$ is the maximum of the size of the incident mass M and ${GM}/{v}_{M}^{2}$. For the background material and comprehensive treatment of Chandrasekhar's formula (1) and its limitations, see Binney & Tremaine (2008). Equation (1) can be interpreted in terms of gravitational wake, namely, the disturbance caused by the motion of M through the medium produces a density enhancement in the incident body's wake. This overdensity decelerates M via gravitational attraction; see Tremaine & Weinberg (1984) and the references cited therein.

Chandrasekhar's formula (1) is a consequence of two-body Newtonian scattering of M from each of the stars in the medium; indeed, it is an approximate result that ignores the attractive gravitational interaction of the stars. In the scattering of M from a star of mass m, the final (t = ∞ ) deflected momentum of M has a component along the initial (t = −∞ ) direction of its incident momentum M v _M that is invariant to linear order in the Newtonian constant G, but decreases to order G² and beyond in accordance with (Binney & Tremaine 2008)

$$\begin{eqnarray}&&M{\left(\delta {{\boldsymbol{v}}}_{M}\right)}_{| | }=-2\,m\,M\,{{\boldsymbol{v}}}_{0}\,\displaystyle \frac{{G}^{2}(m+M)}{{G}^{2}{\left(m+M\right)}^{2}+{b}^{2}{v}_{0}^{4}},\end{eqnarray} \tag{ 2 }$$

where v ₀ = v (t = −∞ ), v = v _M − v _m is the relative velocity, v₀ = ∣ v ₀∣ and b is, as before, the impact parameter. Here, $M{(\delta {{\boldsymbol{v}}}_{M})}_{| | }$ is the net change in the momentum of M along its initial direction of incidence as a consequence of gravitational scattering from m. The net loss of momentum of M along its initial direction of motion is the source of the dynamical friction force.

While numerous studies have verified that Chandrasekhar's formula (1) gives a rather accurate description in many different systems (e.g., van den Bosch et al. 1999; Boylan-Kolchin et al. 2008), the derivation of formula (1) is based on the assumption of an infinite homogeneous background medium and furthermore neglects the self-gravity of the background stars. In the meantime, the calculation of dynamical friction has been extended to inhomogeneous media. In spherically symmetric environments, for example, the phenomenon of "core stalling" has been found, where dynamical friction disappears within the central constant-density cores of massive halos (Read et al. 2006; Banik & van den Bosch 2021a). It has also been shown that a body orbiting outside a self-gravitating system experiences orbital decay and loses angular momentum (Lin & Tremaine 1983). Moreover, the influence of gravitational interaction of the background stars on dynamical friction has been considered in several analytical and numerical studies (see Tremaine & Weinberg 1984; Banik & van den Bosch 2021a, and the references cited therein).

In connection with tides, theoretical studies of dynamical friction in the presence of tidal interactions have generally used the impulse approximation for the sake of simplicity and have thus been restricted to the tidal shocking method (e.g., Spitzer 1958; Colpi & Pallavicini 1998; Gnedin et al. 1999; Binney & Tremaine 2008; Banik & van den Bosch 2021b). The impulse approximation was originally used by Fermi in quantum scattering calculations (Fermi 1936). In reference to tidal interaction, it signifies the change in the velocities of stars as a result of a brief tidal encounter of a galaxy. The time integral of tidal force is the impulse that results in the change of momentum. Within this approximation, the timescale of the encounter between the system and the external tidal source is much shorter than the dynamical timescale associated with the internal dynamics of the system. Naturally, this approximation is restrictive and cannot be applied to arbitrary astrophysical systems. Therefore, dynamical friction studies that deal with ordinary tidal interactions have been mostly based on N-body simulations (e.g., Muñoz et al. 2005; Renaud et al. 2011; Battaglia et al. 2015; Wang et al. 2017; Iorio et al. 2019). Indeed, various N-body simulations of Milky Way (MW) satellites have reported that the stellar components of these dwarf spheroidal (dSph) galaxies are not directly affected by the tidal field of the MW (e.g., Battaglia et al. 2015; Wang et al. 2017, 2019; Iorio et al. 2019); that is, the influence of tidal effects on the stellar kinematics must be rather small.

The main purpose of the present work is to extend Equation (2) in the presence of tidal interactions. Specifically, we study the two-body Newtonian scattering of M from a star of mass m within a stellar system that is under the tidal influence of a distant mass. To simplify matters, we expand the solution of the equation of relative motion of the binary in powers of the Newtonian gravitational constant and work to second order in G, since dynamical friction first appears at this order in accordance with Equation (2). In Section 2, we derive the equation of relative motion in the presence of the tides. To solve this equation, we present a perturbation scheme in powers of G in Section 3 and estimate the influence of tides on dynamical friction. The astrophysical implications of our results are briefly treated in Section 4. Finally, we discuss the limitations of our approach in Section 5.

Imagine a Newtonian two-body system with masses m₁ and m₂ within a stellar medium (of extent D) that is tidally perturbed by a distant galactic mass ${\mathbb{M}}$, ${\mathbb{M}}\gg {m}_{1}+{m}_{2}$. We focus attention on the two-body system and neglect the gravitational interaction of the two bodies with the stars of the background medium. The external source ${\mathbb{M}}$ generates a gravitational potential Φ at the location of the two-body system. The equations of motion of m₁ and m₂ are given in a local Cartesian coordinate system by

$$\begin{eqnarray}&&{m}_{1}\,\displaystyle \frac{{d}^{2}{{\boldsymbol{x}}}_{{m}_{1}}}{{{dt}}^{2}}=-\displaystyle \frac{{{Gm}}_{1}{m}_{2}{\boldsymbol{x}}}{{r}^{3}}-{m}_{1}{\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{1}\,\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{m}_{2}\,\displaystyle \frac{{d}^{2}{{\boldsymbol{x}}}_{{m}_{2}}}{{{dt}}^{2}}=\displaystyle \frac{{{Gm}}_{1}{m}_{2}{\boldsymbol{x}}}{{r}^{3}}-{m}_{2}{\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{2},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{x}}:= {{\boldsymbol{x}}}_{{m}_{1}}-{{\boldsymbol{x}}}_{{m}_{2}},\qquad r:= | {\boldsymbol{x}}| .\end{eqnarray} \tag{ 5 }$$

We assume that $r\leqslant D\ll {\mathbb{R}}$, where ${\mathbb{R}}$ is the distance of the binary system to ${\mathbb{M}}$. We expand the gravitational force of ${\mathbb{M}}$ on the two-body system about its center of mass

$$\begin{eqnarray}&&{\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{1}^{i}={\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{\mathrm{cm}}^{i}+{K}^{i}{}_{j}({{\boldsymbol{x}}}_{\mathrm{cm}})({x}_{{m}_{1}}^{j}-{x}_{\mathrm{cm}}^{j})+\cdots \end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{2}^{i}={\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{\mathrm{cm}}^{i}+{K}^{i}{}_{j}({{\boldsymbol{x}}}_{\mathrm{cm}})({x}_{{m}_{2}}^{j}-{x}_{\mathrm{cm}}^{j})+\cdots ,\end{eqnarray} \tag{ 7 }$$

where K_ij = ∇_i ∇_j Φ is the symmetric and traceless tidal matrix evaluated at the center of mass of the binary system

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{\mathrm{cm}}=\displaystyle \frac{{m}_{1}{{\boldsymbol{x}}}_{{m}_{1}}+{m}_{2}{{\boldsymbol{x}}}_{{m}_{2}}}{{m}_{1}+{m}_{2}}.\end{eqnarray} \tag{ 8 }$$

We assume that ${(D/{\mathbb{R}})}^{2}$ is negligibly small compared to unity; therefore, we henceforth neglect terms of second and higher orders in Equations (6) and (7). Then, the equations of motion imply

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{{\boldsymbol{x}}}_{\mathrm{cm}}}{{{dt}}^{2}}=-{\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{\mathrm{cm}}\,\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{x}^{i}}{{{dt}}^{2}}=-\displaystyle \frac{G({m}_{1}+{m}_{2}){x}^{i}}{{r}^{3}}-{K}^{i}{}_{j}({{\boldsymbol{x}}}_{\mathrm{cm}}){x}^{j}.\end{eqnarray} \tag{ 10 }$$

Appendix A contains the natural extension of the two-body approach developed here to an N-body system. However, it is not known how to derive the dynamical friction force within this scheme; therefore, in conformity with Chandrasekhar's approach, we hereafter ignore the gravitational interaction between the background stars.

To interpret Equations (9) and (10), we note that our dynamical system has a slow motion given by Equation (9) and a fast motion given by Equation (10). The slow motion involves the Keplerian revolution of the whole binary system about ${\mathbb{M}}$. The fast motion involves the internal scattering orbit of the tidally perturbed gravitational two-body system. We assume that during the fast motion, the slow motion can be considered to be essentially uniform. That is, the deviation of the Keplerian orbit of the center of mass of the binary from a straight line with constant velocity can be neglected during the fast motion of the binary. This approximation scheme can be naturally extended to include the background stellar medium in order to facilitate the calculation of dynamical friction. Moreover, in Equation (10), we neglect the temporal dependence of the tidal matrix Kⁱ _j ( x _cm) in our approach. In general, the temporal variation of the tidal field can lead to the possibility of resonance with the internal pulsation of the system; in fact, such a resonant orbital coupling has been investigated in Kuhn & Miller (1989) and Fleck & Kuhn (2003).

In the inertial reference frame that moves with the constant velocity of the center of mass, mass m₁ with state $({{\boldsymbol{x}}}_{{m}_{1}},{{\boldsymbol{v}}}_{{m}_{1}})$ interacts gravitationally with mass m₂ with state $({{\boldsymbol{x}}}_{{m}_{2}},{{\boldsymbol{v}}}_{{m}_{2}})$ and the relative motion is given by Equation (10). The solution of this equation and the scattering process will be described in detail in the next section. However, it is simple to conclude from the uniform motion of the center of mass and the definition of relative motion that the net change in the velocities of m₁ and m₂ are given by $\delta {{\boldsymbol{v}}}_{{m}_{1}}={m}_{2}\delta {\boldsymbol{v}}/({m}_{1}+{m}_{2})$ and $\delta {{\boldsymbol{v}}}_{{m}_{2}}=-{m}_{1}\delta {\boldsymbol{v}}/({m}_{1}+{m}_{2})$, where ${\boldsymbol{v}}:= \dot{{\boldsymbol{x}}}$ is the relative velocity.

It remains to evaluate the external force acting on the center of mass in Equation (9) and the tidal matrix in Equation (10). Let ${\mathbb{R}}\,{\boldsymbol{n}}$ be the vector that connects ${\mathbb{M}}$ to the center of mass of the binary system and n be the corresponding unit vector. Then,

$$\begin{eqnarray}&&{\left({\rm{\nabla }}{\rm{\Phi }}\right)}_{\mathrm{cm}}=\displaystyle \frac{G\,{\mathbb{M}}}{{{\mathbb{R}}}^{2}}\,{\boldsymbol{n}}\,\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}&&{K}^{i}{}_{j}({{\boldsymbol{x}}}_{\mathrm{cm}})=\displaystyle \frac{G\,{\mathbb{M}}}{{{\mathbb{R}}}^{3}}({\delta }_{j}^{i}-3{n}^{i}{n}_{j}).\end{eqnarray} \tag{ 12 }$$

In the local Cartesian coordinate system x = (x, y, z) established in the inertial frame that moves with a constant velocity of the center of mass of the binary, we can choose the unit vector n to point in the z direction with no loss in generality. The equation for relative motion of the binary (10) now takes the form

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\boldsymbol{x}}}{{{dt}}^{2}}=-G\left(\displaystyle \frac{M^{\prime} }{{r}^{3}}+\lambda \right){\boldsymbol{x}}+3G\lambda ({\boldsymbol{x}}\cdot {\boldsymbol{n}}){\boldsymbol{n}},\end{eqnarray} \tag{ 13 }$$

where we have introduced

$$\begin{eqnarray}&&{m}_{1}+{m}_{2}:= M^{\prime} ,\qquad \lambda := \displaystyle \frac{{\mathbb{M}}}{{{\mathbb{R}}}^{3}}.\end{eqnarray} \tag{ 14 }$$

Inspection of the equation of relative motion of the binary (13) reveals that the external tidal force becomes comparable to the internal Newtonian inverse-square force at r ≈ R₀, where R₀ is the tidal radius given by

$$\begin{eqnarray}&&\displaystyle \frac{M^{\prime} }{{R}_{0}^{3}}=\displaystyle \frac{{\mathbb{M}}}{{{\mathbb{R}}}^{3}},\qquad {R}_{0}={\left(\displaystyle \frac{M^{\prime} }{\lambda }\right)}^{1/3}.\end{eqnarray} \tag{ 15 }$$

Indeed, for r > R₀ the tidal force is dominant, while for r < R₀ the Newtonian attractive force between m₁ and m₂ is dominant. The main purpose of the present paper is to determine how the tidal interaction modifies the dynamical friction force that is a consequence of gravitational drag as m₁ scatters from m₂. Therefore, we must limit our considerations to regions within a sphere of radius R₀. In such regions, the tidal force is smaller than the internal Newtonian force of attraction and can be treated as a linear perturbation on the binary. This approach is consistent with the limiting situation in which the exterior tidal force disappears; that is, for λ → 0, we find R₀ → ∞ and our treatment of dynamical friction reduces to the standard method (Chandrasekhar 1943; Tremaine & Weinberg 1984; Binney & Tremaine 2008). We treat the influence of tides as a first-order perturbation; therefore, we expect that the external tidal interaction will have a generally small effect on the internal dynamical friction force. The main purpose of this work is to determine the extent to which dynamical friction is modified by external tides.

To calculate the dynamical friction force in the presence of tides, we seek a perturbative solution of the equation of relative motion (13) within the domain bounded by the tidal radius R₀. Let us first observe that Equation (13) is invariant under a constant temporal translation as well as a constant rotation in the (x, y) plane. We use these symmetries below.

3.1. Initial Motion in the (x, y) Plane

Let us first assume for the sake of simplicity that the unperturbed uniform relative motion occurs in the (x, y) plane. We can therefore rotate the Cartesian coordinate system about the z-axis such that the unperturbed relative motion occurs along the x-axis from (−x₀, y₀, z₀) to (x₀, y₀, z₀). Here, x₀ = v₀ t₀, y₀, and z₀ are constants that characterize the unperturbed motion. That is,

$$\begin{eqnarray}&&x(t)={v}_{0}\,t,\,y={y}_{0},\,z={z}_{0},\,t=-{t}_{0}\to t={t}_{0}\,\end{eqnarray} \tag{ 16 }$$

such that

$$\begin{eqnarray}&&{r}_{0}={\left({v}_{0}^{2}\,{t}_{0}^{2}+{y}_{0}^{2}+{z}_{0}^{2}\right)}^{1/2}\lt {R}_{0}.\end{eqnarray} \tag{ 17 }$$

Let us note here that we have used the invariance under time translation to set x(t = 0) = 0. Henceforth, we will express the unperturbed motion in the background (x, y, z) coordinate system in the form

$$\begin{eqnarray}&&{\boldsymbol{\xi }}(t)={\boldsymbol{b}}+{{\boldsymbol{v}}}_{0}\,t,\,{\boldsymbol{b}}=(0,{y}_{0},{z}_{0}),\,{{\boldsymbol{v}}}_{0}=({v}_{0},0,0),\end{eqnarray} \tag{ 18 }$$

where b · v ₀ = 0 and ∣ b ∣ = b is the impact parameter.

Let us recall here that dynamical friction is a phenomenon that first appears at order G² in the standard treatment (Chandrasekhar 1943; Binney & Tremaine 2008). We therefore solve Equation (13) perturbatively to second order in the gravitational constant G and to linear order in the tidal parameter λ. Specifically, we assume a solution of the form

$$\begin{eqnarray}&&{\boldsymbol{x}}(t)={\boldsymbol{\xi }}(t)+G\,{{\boldsymbol{x}}}_{1}(t)+{G}^{2}\,{{\boldsymbol{x}}}_{2}(t)+\cdots ,\end{eqnarray} \tag{ 19 }$$

where we have treated G as the expansion parameter. Similarly, we write

$$\begin{eqnarray}&&{\boldsymbol{v}}=\displaystyle \frac{d{\boldsymbol{x}}}{{dt}}={{\boldsymbol{v}}}_{0}+G\,\displaystyle \frac{d{{\boldsymbol{x}}}_{1}}{{dt}}+{G}^{2}\,\displaystyle \frac{d{{\boldsymbol{x}}}_{2}}{{dt}}+\cdots .\end{eqnarray} \tag{ 20 }$$

Henceforth, we drop terms of order G³ and higher. To proceed, it is useful to introduce the quantities

$$\begin{eqnarray}&&U(t):= {\left({\boldsymbol{\xi }}\cdot {\boldsymbol{\xi }}\right)}^{1/2}={\left({v}_{0}^{2}\,{t}^{2}+{b}^{2}\right)}^{1/2}\,\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}&&W(t):= \displaystyle \frac{{\boldsymbol{\xi }}\cdot {{\boldsymbol{x}}}_{1}}{{U}^{2}},\end{eqnarray} \tag{ 22 }$$

in terms of which we can write

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{3}}=\displaystyle \frac{1}{{U}^{3}}(1-3{GW})+O({G}^{2}).\end{eqnarray} \tag{ 23 }$$

Plugging these expressions in Equation (13), we find to first order in G

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{{\boldsymbol{x}}}_{1}}{{{dt}}^{2}}=-\left(\lambda +\displaystyle \frac{M^{\prime} }{{U}^{3}}\right){\boldsymbol{\xi }}(t)+3\lambda ({\boldsymbol{\xi }}\cdot {\boldsymbol{n}}){\boldsymbol{n}}\end{eqnarray} \tag{ 24 }$$

and to second order in G

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{{\boldsymbol{x}}}_{2}}{{{dt}}^{2}}=-\left(\lambda +\displaystyle \frac{M^{\prime} }{{U}^{3}}\right){{\boldsymbol{x}}}_{1}+3\lambda ({{\boldsymbol{x}}}_{1}\cdot {\boldsymbol{n}}){\boldsymbol{n}}+3\displaystyle \frac{M^{\prime} \,W}{{U}^{3}}{\boldsymbol{\xi }}.\end{eqnarray} \tag{ 25 }$$

We must now solve these equations with the boundary conditions that initially at t = −t₀, we have x ₁ = 0, x ₂ = 0, ${\dot{{\boldsymbol{x}}}}_{1}=0$ and ${\dot{{\boldsymbol{x}}}}_{2}=0$. Here, ${\dot{x}}_{1}:= {{dx}}_{1}/{dt}$, etc.

Using standard integrals, given for convenience in Appendix B, we find

$$\begin{eqnarray}&&{\dot{x}}_{1}=\displaystyle \frac{M^{\prime} }{{v}_{0}}({U}^{-1}-{U}_{0}^{-1})-\displaystyle \frac{1}{2}\,\displaystyle \frac{\lambda }{{v}_{0}}({u}^{2}-{u}_{0}^{2}),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\dot{y}}_{1}=-\displaystyle \frac{M^{\prime} \,{y}_{0}}{{v}_{0}{b}^{2}}(u\,{U}^{-1}+{u}_{0}\,{U}_{0}^{-1})-\displaystyle \frac{\lambda }{{v}_{0}}\,{y}_{0}(u+{u}_{0}),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\dot{z}}_{1}=-\displaystyle \frac{M^{\prime} \,{z}_{0}}{{v}_{0}{b}^{2}}(u\,{U}^{-1}+{u}_{0}\,{U}_{0}^{-1})+2\displaystyle \frac{\lambda }{{v}_{0}}\,{z}_{0}(u+{u}_{0}),\end{eqnarray} \tag{ 28 }$$

where, for the sake of convenience, we have introduced

$$\begin{eqnarray}\begin{array}{rcl}u & := & {v}_{0}\,t,\quad {u}_{0}:= {x}_{0}={v}_{0}\,{t}_{0},\\ {U}_{0} & := & {r}_{0}={\left({u}_{0}^{2}+{b}^{2}\right)}^{1/2}.\end{array}\end{eqnarray} \tag{ 29 }$$

Let us note from Equation (26) that ${\dot{x}}_{1}({t}_{0})=0$ and hence $\delta {\dot{x}}_{1}=0$, which is the expected result at linear order in G from the standard treatment of dynamical friction. In general, the situation is different in the presence of tides, as will be demonstrated in the last part of this section; however, the tidal contribution does vanish in certain special cases that include initial motions along the x and z axes.

Integration of Equations (26)–(28) with the appropriate initial conditions results in

$$\begin{eqnarray}&&{x}_{1}=\displaystyle \frac{M^{\prime} }{{v}_{0}^{2}}\,{{\mathbb{A}}}_{1}-\displaystyle \frac{1}{6}\,\displaystyle \frac{\lambda }{{v}_{0}^{2}}\,{{\mathbb{B}}}_{1},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{y}_{1}=-\displaystyle \frac{M^{\prime} \,{y}_{0}}{{v}_{0}^{2}\,{b}^{2}}\,{{\mathbb{A}}}_{2}-\displaystyle \frac{1}{2}\,\displaystyle \frac{\lambda }{{v}_{0}^{2}}\,{y}_{0}\,{{\mathbb{B}}}_{2},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{z}_{1}=-\displaystyle \frac{M^{\prime} \,{z}_{0}}{{v}_{0}^{2}\,{b}^{2}}\,{{\mathbb{A}}}_{2}+\displaystyle \frac{\lambda }{{v}_{0}^{2}}\,{z}_{0}\,{{\mathbb{B}}}_{2}.\end{eqnarray} \tag{ 32 }$$

Here, we have introduced, for the sake of convenience, the functions

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{A}}}_{1}(u) & := & \mathrm{ln}\left(\displaystyle \frac{U+u}{{U}_{0}-{u}_{0}}\right)-\displaystyle \frac{u+{u}_{0}}{{U}_{0}},\\ {{\mathbb{A}}}_{2}(u) & := & U-{U}_{0}+\displaystyle \frac{{u}_{0}}{{U}_{0}}(u+{u}_{0})\end{array}\end{eqnarray} \tag{ 33 }$$

and

$$\begin{eqnarray}&&{{\mathbb{B}}}_{1}(u):= {u}^{3}-3{u}_{0}^{2}\,u-2{u}_{0}^{3},\qquad {{\mathbb{B}}}_{2}(u):= {\left(u+{u}_{0}\right)}^{2}.\end{eqnarray} \tag{ 34 }$$

To go further, we can integrate the second-order perturbation equations for (x₂, y₂, z₂). On the other hand, for the purpose of calculating dynamical friction, we only need the equation for x₂; that is, we must integrate the first component of Equation (25) from t = − t₀ to t = t₀. More precisely, we need to calculate

$$\begin{eqnarray}\begin{array}{rcl}\delta {v}_{x} & = & {v}_{x}({t}_{0})-{v}_{x}(-{t}_{0})={v}_{x}({t}_{0})-{v}_{0}\\ & = & G\,\displaystyle \frac{{{dx}}_{1}}{{dt}}({t}_{0})+{G}^{2}\,\displaystyle \frac{{{dx}}_{2}}{{dt}}({t}_{0})={G}^{2}\,\displaystyle \frac{{{dx}}_{2}}{{dt}}({t}_{0}),\end{array}\end{eqnarray} \tag{ 35 }$$

since we already know that ${\dot{x}}_{1}({t}_{0})=0$, i.e., the term linear in G vanishes in this case. We treat the tidal parameter λ to first order in our perturbation analysis; therefore, after some algebra we find

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{x}_{2}}{{{dt}}^{2}}=\displaystyle \frac{M{{\prime} }^{2}}{{v}_{0}^{2}{U}_{0}{U}^{5}}\,{ \mathcal A }+\displaystyle \frac{M^{\prime} \,\lambda }{{v}_{0}^{2}}\left(\displaystyle \frac{{{ \mathcal B }}_{1}}{{U}^{5}}+{{ \mathcal B }}_{2}\right)+O({\lambda }^{2}).\end{eqnarray} \tag{ 36 }$$

Here,

$$\begin{eqnarray}&&{ \mathcal A }={U}_{0}[(2{u}^{2}-{b}^{2}){{\mathbb{A}}}_{1}-3u\,{{\mathbb{A}}}_{2}],\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{{ \mathcal B }}_{1}=-\displaystyle \frac{1}{6}(2{u}^{2}-{b}^{2}){{\mathbb{B}}}_{1}+\displaystyle \frac{3}{2}(2{z}_{0}^{2}-{y}_{0}^{2})u\,{{\mathbb{B}}}_{2},\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{{ \mathcal B }}_{2}=-{{\mathbb{A}}}_{1}.\end{eqnarray} \tag{ 39 }$$

Next, we integrate Equation (36) from t = − t₀ to t = t₀ with the boundary condition that ${\dot{x}}_{2}(-{t}_{0})=0$. The result is

$$\begin{eqnarray}&&{\dot{x}}_{2}({t}_{0})=\displaystyle \frac{2M{{\prime} }^{2}{u}_{0}}{{v}_{0}^{3}{U}_{0}^{3}}\,{ \mathcal S }+\displaystyle \frac{2M^{\prime} \,\lambda \,{u}_{0}}{{v}_{0}^{3}}{ \mathcal T },\end{eqnarray} \tag{ 40 }$$

where

$$\begin{eqnarray}&&{ \mathcal S }={\mathbb{C}}-\displaystyle \frac{{u}_{0}^{3}}{{b}^{2}{U}_{0}},\qquad { \mathcal T }={\mathbb{C}}+\displaystyle \frac{{u}_{0}^{3}}{3{b}^{2}{U}_{0}^{3}}(5{z}_{0}^{2}-4{y}_{0}^{2}).\end{eqnarray} \tag{ 41 }$$

Here, ${\mathbb{C}}(b,{u}_{0})\leqslant 0$ is given by

$$\begin{eqnarray}&&{\mathbb{C}}(b,{u}_{0}):= \mathrm{ln}\left(\displaystyle \frac{{U}_{0}-{u}_{0}}{b}\right)+\displaystyle \frac{{u}_{0}}{{U}_{0}}.\end{eqnarray} \tag{ 42 }$$

We can now derive the modification of Equation (2) in the presence of tides to second order in G. Suppose that m₁ = M, m₂ = m and $M+m=M^{\prime}$. Furthermore, our perturbative approach implies

$$\begin{eqnarray}&&{\left(\delta {{\boldsymbol{v}}}_{M}\right)}_{| | }=(m/M^{\prime} ){\left(\delta {\boldsymbol{v}}\right)}_{| | }=(m/M^{\prime} )\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{{v}_{0}}{G}^{2}{\dot{x}}_{2}({t}_{0}).\end{eqnarray} \tag{ 43 }$$

Using Equation (40), we find

$$\begin{eqnarray}&&{\left(\delta {{\boldsymbol{v}}}_{M}\right)}_{| | }=2{G}^{2}m\,\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{{v}_{0}^{4}}\left[\displaystyle \frac{M^{\prime} \,{u}_{0}}{{U}_{0}^{3}}\,{ \mathcal S }+\lambda \,{u}_{0}{ \mathcal T }\right],\end{eqnarray} \tag{ 44 }$$

which is the analog of Equation (2) in the presence of tides. We note here the remarkable fact that ${ \mathcal S }$ is negative; more specifically, ${ \mathcal S }$ as a function of u₀/b starts out from zero at u₀/b = 0 with vanishing slope and then monotonically decreases to −∞ as u₀/b → ∞. Moreover, let us note that in the presence of tides t₀ > 0 is finite and is limited by Equation (17); however, in the absence of tides λ = 0 and we can let t₀ → ∞ as in the standard treatment (Binney & Tremaine 2008). Therefore, in the absence of tides, t₀ → ∞ , u₀ = v₀ t₀ → ∞ and ${ \mathcal S }\to -{u}_{0}^{2}/{b}^{2}$. Then, Equation (44) reduces to

$$\begin{eqnarray}&&{\left(\delta {{\boldsymbol{v}}}_{M}\right)}_{| | }=-\displaystyle \frac{2{G}^{2}m(m+M){{\boldsymbol{v}}}_{0}}{{b}^{2}\,{v}_{0}^{4}},\end{eqnarray} \tag{ 45 }$$

in agreement with Equation (2) to second order in G.

A comment is in order here regarding the fact that t₀ = u₀/v₀ is essentially a free parameter in Equation (44). Indeed, in the application of Chandrasekhar's formula (1) to actual astrophysical systems, the scattering of M from each star of mass m in the system does involve a different finite t₀ that we simply ignore. Instead, in each such case, we employ the result obtained from the ideal case with t₀ = ∞ . This approximation does not ordinarily encounter any obstacles. However, the presence of a finite tidal radius R₀ in the case under consideration implies that t₀ is a free parameter subject to the restriction contained in the inequality ${v}_{0}^{2}\,{t}_{0}^{2}+{b}^{2}\lt {R}_{0}^{2}$.

To find the analog of Chandrasekhar's formula (1) in the present case, we need to extend Equation (44) for an arbitrary star m to the background stars by taking the flux of the stars into account and integrating over the appropriate state space of the background stars. Let f( x _m , v _m , t) be the state-space number density of background stars. In the idealized infinite homogeneous background medium employed in the original derivation of Equation (1), the distribution function f reduces to f₀, which is independent of temporal and spatial variables and is an isotropic function of v _m ; however, in the presence of tides, the situation is quite different. Starting from an initial distribution function, it is possible to take the time dependence of tides into account and then basically integrate the tidal force (per unit mass of a star) over a certain interval of time and in this way come up with the change in the velocity of the star due to the corresponding tidal impulse. In principle, the evolution of the distribution function due to the presence of tides can thereby be determined (Spitzer 1958; Gnedin et al. 1999; Binney & Tremaine 2008; Banik & van den Bosch 2021b). The determination of the appropriate distribution function in the case under consideration requires a separate investigation and is beyond the scope of the present work.

To get a rough estimate of the influence of tides on dynamical friction, we henceforth ignore the difference between f( x _m , v _m , t) and f₀(v_m ). The net rate of change of momentum M v _M per unit time is then given by

$$\begin{eqnarray}\begin{array}{rcl}M\displaystyle \frac{d{{\boldsymbol{v}}}_{M}}{{dt}} & \approx & 2{G}^{2}{mM}\displaystyle \int \left[\displaystyle \frac{M^{\prime} {u}_{0}}{{U}_{0}^{3}}{ \mathcal S }+\lambda {u}_{0}{ \mathcal T }\right]{{dy}}_{0}\,{{dz}}_{0}\\ & & \times \,\displaystyle \int \displaystyle \frac{{{\boldsymbol{v}}}_{M}-{{\boldsymbol{v}}}_{m}}{| {{\boldsymbol{v}}}_{M}-{{\boldsymbol{v}}}_{m}{| }^{3}}{f}_{0}({v}_{m}){d}^{3}{v}_{m},\end{array}\end{eqnarray} \tag{ 46 }$$

where v ₀ has been replaced by v _M − v _m , as in the standard procedure (Binney & Tremaine 2008). For an isotropic function f₀ = f₀(v_m ), Newton's shell theorem implies (Binney & Tremaine 2008)

$$\begin{eqnarray}&&\int \displaystyle \frac{{{\boldsymbol{v}}}_{M}-{{\boldsymbol{v}}}_{m}}{| {{\boldsymbol{v}}}_{M}-{{\boldsymbol{v}}}_{m}{| }^{3}}{f}_{0}\,{d}^{3}{v}_{m}=4\pi \displaystyle \frac{{{\boldsymbol{v}}}_{M}}{{v}_{M}^{3}}{\int }_{0}^{{v}_{M}}{f}_{0}\,{v}_{m}^{2}\,{{dv}}_{m}.\end{eqnarray} \tag{ 47 }$$

Using

$$\begin{eqnarray}&&\rho (\lt {v}_{M})=4\pi m{\int }_{0}^{{v}_{M}}{f}_{0}({v}_{m}){v}_{m}^{2}\,{{dv}}_{m},\end{eqnarray} \tag{ 48 }$$

we finally arrive at

$$\begin{eqnarray}\begin{array}{rcl}M\displaystyle \frac{d{{\boldsymbol{v}}}_{M}}{{dt}} & \approx & 2{G}^{2}M\displaystyle \frac{{{\boldsymbol{v}}}_{M}}{{v}_{M}^{3}}\,\rho (\lt {v}_{M})\\ & & \times \,\displaystyle \int \left[\displaystyle \frac{M^{\prime} \,{u}_{0}}{{U}_{0}^{3}}\,{ \mathcal S }+\lambda \,{u}_{0}{ \mathcal T }\right]{{dy}}_{0}\,{{dz}}_{0}.\end{array}\end{eqnarray} \tag{ 49 }$$

Let us proceed with the integration over dy₀ ∧ dz₀. Introducing the azimuthal angle φ, we have

$$\begin{eqnarray}&&{y}_{0}=b\,\cos \varphi ,\qquad {z}_{0}=b\,\sin \varphi .\end{eqnarray} \tag{ 50 }$$

Then, dy₀ ∧ dz₀ = b db ∧ d φ and one can show that

$$\begin{eqnarray}\begin{array}{rcl}{\displaystyle \int }_{{b}_{\min }}^{{b}_{\max }}\displaystyle \frac{{u}_{0}\,{ \mathcal S }}{{U}_{0}^{3}}\,b\,{db}\,d\varphi & = & {\mathbb{S}}({b}_{\max },{u}_{0}{| }_{{b}_{\max }})\\ & & -\,{\mathbb{S}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }}),\end{array}\end{eqnarray} \tag{ 51 }$$

where

$$\begin{eqnarray}&&{\mathbb{S}}(b,{u}_{0}):= -\displaystyle \frac{2\pi {u}_{0}}{{U}_{0}}\,{\mathbb{C}}(b,{u}_{0}).\end{eqnarray} \tag{ 52 }$$

Let us observe that ${\mathbb{S}}(b,{u}_{0})$ vanishes for u₀ = 0, while for u₀ > 0, ${\mathbb{S}}(b,{u}_{0})$ is a positive function of b/u₀ that monotonically decreases from ∞ at b/u₀ = 0 and tends to zero as b/u₀ → ∞ .

In the absence of tides, λ = 0, u₀ → ∞ and ${ \mathcal S }\to -{u}_{0}^{2}/{b}^{2};$ then, as expected, we recover $-\mathrm{ln}({b}_{\max }/{b}_{\min })$ from Equation (51). However, in the presence of tides ${u}_{0}^{2}+{b}^{2}={r}_{0}^{2}\lt {R}_{0}^{2}$ according to Equation (17); therefore, in Equation (51), ${\mathbb{S}}({b}_{\max },{u}_{0}{| }_{{b}_{\max }})$ vanishes, since ${b}_{\max }={r}_{0}$ and hence ${u}_{0}{| }_{{b}_{\max }}=0$. To interpret this result physically, let us first note that u₀ = v₀ t₀ vanishes when either v₀ = 0 or t₀ = 0. It is clear from Equation (2) that if v₀ = 0, there is no scattering and dynamical friction force vanishes. On the other hand, if t₀ vanishes, there is no time for interaction and hence no change in momentum is possible. Thus, for the nontidal part of the integral in Equation (49), we find

$$\begin{eqnarray}&&\int \displaystyle \frac{M^{\prime} \,{u}_{0}}{{U}_{0}^{3}}\,{ \mathcal S }\,{{dy}}_{0}\,{{dz}}_{0}=-(m+M){\mathbb{S}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }}),\end{eqnarray} \tag{ 53 }$$

where

$$\begin{eqnarray}&&{b}_{\min }^{2}+{u}_{0}{| }_{{b}_{\min }}^{2}={r}_{0}^{2}.\end{eqnarray} \tag{ 54 }$$

In our case, the nontidal part is thus − (m + M) times the analog of the Coulomb logarithm. As in the standard treatment of dynamical friction (Binney & Tremaine 2008), we assume ${b}_{\min }$ is the maximum of the size of the incident mass M and ${GM}/{v}_{M}^{2}$, while ${r}_{0}=\min (D,{R}_{0})$, where D is the diameter of the galactic system. In the absence of tides, the tidal radius R₀ goes to infinity and, consequently, we have r₀ = D.

Let us next consider the tidal part of Equation (49), namely, the integral of $\lambda \,{u}_{0}\,{ \mathcal T }$ that turns out to be

$$\begin{eqnarray}\begin{array}{rcl}\lambda \,{\displaystyle \int }_{{b}_{\min }}^{{b}_{\max }}{u}_{0}\,{ \mathcal T }b\,{db}\,d\varphi & = & \pi \,\lambda [{\mathbb{T}}({b}_{\max },{u}_{0}{| }_{{b}_{\max }})\\ & & -\,{\mathbb{T}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})],\end{array}\end{eqnarray} \tag{ 55 }$$

where ${\mathbb{T}}$ is proportional to u₀ and is given by

$$\begin{eqnarray}&&{\mathbb{T}}(b,{u}_{0}):= {b}^{2}\,{u}_{0}\left({\mathbb{C}}+\displaystyle \frac{2}{3}\displaystyle \frac{{u}_{0}^{3}}{{b}^{2}{U}_{0}}\right).\end{eqnarray} \tag{ 56 }$$

As illustrated in Figure 1, the function $3{\mathbb{T}}(b,{u}_{0})/(2{u}_{0}^{3})$ versus b/u₀ starts from unity with zero slope at b/u₀ = 0 and monotonically decreases to zero as b/u₀ → ∞. As noted before, ${u}_{0}{| }_{{b}_{\max }}=0$ implies that ${\mathbb{T}}({b}_{\max },{u}_{0}{| }_{{b}_{\max }})=0;$ therefore,

$$\begin{eqnarray}&&\lambda \,{\int }_{{b}_{\min }}^{{b}_{\max }}{u}_{0}\,{ \mathcal T }b\,{db}\,d\varphi =-\pi \,\lambda \,{\mathbb{T}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }}).\end{eqnarray} \tag{ 57 }$$

Zoom In Zoom Out Reset image size

Finally, to determine the relative strength of the tidal part compared to the nontidal part in the calculation of dynamical friction force in the case of motion along the x-axis, we must calculate Ψ_x ,

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{x}=\displaystyle \frac{\pi \,\lambda }{M^{\prime} }\displaystyle \frac{{\mathbb{T}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})}{{\mathbb{S}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})},\qquad {R}_{0}={\left(\displaystyle \frac{M^{\prime} }{\lambda }\right)}^{1/3},\end{eqnarray} \tag{ 58 }$$

once we ignore the influence of tides on the distribution function of the stars. In Section 4, we estimate Ψ_x for the case of the Fornax dSph galaxy tidally perturbed by the Galaxy (MW).

3.2. Initial Motion Along the z-axis

Let us assume that the unperturbed motion is given by

$$\begin{eqnarray}\begin{array}{rcl}x & = & {\bar{x}}_{0},\,y={\bar{y}}_{0},\,z={v}_{0}\,t,{\bar{x}}_{0}^{2}+{\bar{y}}_{0}^{2}={b}^{2},\\ t & = & -{t}_{0}\to t={t}_{0}\end{array}\end{eqnarray} \tag{ 59 }$$

such that ${r}_{0}={U}_{0}={({u}_{0}^{2}+{b}^{2})}^{1/2}\lt {R}_{0}$ as before. That is,

$$\begin{eqnarray}&&{\boldsymbol{\xi }}(t)={\boldsymbol{b}}+{{\boldsymbol{v}}}_{0}\,t,\,{\boldsymbol{b}}=({\bar{x}}_{0},{\bar{y}}_{0},0),\,{{\boldsymbol{v}}}_{0}=(0,0,{v}_{0}).\end{eqnarray} \tag{ 60 }$$

Let us note here that we have used the invariance under time translation to set z(t = 0) = 0.

Following the same procedures as in the previous subsection, we can obtain x ₁ and x ₂ for the perturbed motion. In the case of the linear deviation x ₁, an inspection of Equation (13) reveals that we can find the results in this case from Equations (30)–(32) for motion along the x direction by simply letting λ → − 2 λ. Indeed,

$$\begin{eqnarray}&&{x}_{1}=-\displaystyle \frac{M^{\prime} \,{\bar{x}}_{0}}{{v}_{0}^{2}\,{b}^{2}}\,{{\mathbb{A}}}_{2}-\displaystyle \frac{1}{2}\,\displaystyle \frac{\lambda }{{v}_{0}^{2}}\,{\bar{x}}_{0}\,{{\mathbb{B}}}_{2},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&{y}_{1}=-\displaystyle \frac{M^{\prime} \,{\bar{y}}_{0}}{{v}_{0}^{2}\,{b}^{2}}\,{{\mathbb{A}}}_{2}-\displaystyle \frac{1}{2}\,\displaystyle \frac{\lambda }{{v}_{0}^{2}}\,{\bar{y}}_{0}\,{{\mathbb{B}}}_{2},\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&{z}_{1}=\displaystyle \frac{M^{\prime} }{{v}_{0}^{2}}\,{{\mathbb{A}}}_{1}+\displaystyle \frac{1}{3}\,\displaystyle \frac{\lambda }{{v}_{0}^{2}}\,{{\mathbb{B}}}_{1}.\end{eqnarray} \tag{ 63 }$$

We are interested in the third component of Equation (25) in the present case. After some straightforward calculations, we find

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{z}_{2}}{{{dt}}^{2}}=\displaystyle \frac{M{{\prime} }^{2}}{{v}_{0}^{2}{U}_{0}{U}^{5}}\,\bar{{ \mathcal A }}+\displaystyle \frac{M^{\prime} \,\lambda }{{v}_{0}^{2}}\left(\displaystyle \frac{{\bar{{ \mathcal B }}}_{1}}{{U}^{5}}+{\bar{{ \mathcal B }}}_{2}\right)+O({\lambda }^{2}).\end{eqnarray} \tag{ 64 }$$

Here, the nontidal part is the same as in the case of motion along the x direction by the spherical symmetry of the Newtonian gravitational potential of a point mass; that is,

$$\begin{eqnarray}&&\bar{{ \mathcal A }}={ \mathcal A },\end{eqnarray} \tag{ 65 }$$

while changes occur in the tidal part. Indeed,

$$\begin{eqnarray}&&{\bar{{ \mathcal B }}}_{1}=\displaystyle \frac{1}{3}(2{u}^{2}-{b}^{2}){{\mathbb{B}}}_{1}-\displaystyle \frac{3}{2}\,{b}^{2}u\,{{\mathbb{B}}}_{2},\quad {\bar{{ \mathcal B }}}_{2}=-2\,{{ \mathcal B }}_{2}.\end{eqnarray} \tag{ 66 }$$

It follows that in this case

$$\begin{eqnarray}&&{\left(\delta {{\boldsymbol{v}}}_{M}\right)}_{| | }=2{G}^{2}m\,\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{{v}_{0}^{4}}\left[\displaystyle \frac{M^{\prime} \,{u}_{0}}{{U}_{0}^{3}}\,{ \mathcal S }+\lambda \,{u}_{0}\bar{{ \mathcal T }}\right],\end{eqnarray} \tag{ 67 }$$

where

$$\begin{eqnarray}&&\bar{{ \mathcal T }}=-2\,{\mathbb{C}}-\displaystyle \frac{{u}_{0}^{3}}{3{U}_{0}^{3}}.\end{eqnarray} \tag{ 68 }$$

Following essentially the same arguments as before leads to the next step which involves

$$\begin{eqnarray}\begin{array}{rcl}M\displaystyle \frac{d{{\boldsymbol{v}}}_{M}}{{dt}} & \approx & 2{G}^{2}M\displaystyle \frac{{{\boldsymbol{v}}}_{M}}{{v}_{M}^{3}}\,\rho (\lt {v}_{M})\\ & & \times \,\displaystyle \int \left[\displaystyle \frac{M^{\prime} \,{u}_{0}}{{U}_{0}^{3}}\,{ \mathcal S }+\lambda \,{u}_{0}\bar{{ \mathcal T }}\right]d{\bar{x}}_{0}\,d{\bar{y}}_{0}.\end{array}\end{eqnarray} \tag{ 69 }$$

With ${\bar{x}}_{0}=b\,\cos \bar{\varphi }$, ${\bar{y}}_{0}=b\,\sin \bar{\varphi }$, and $d{\bar{x}}_{0}\wedge d{\bar{y}}_{0}=b\,{db}\wedge d\bar{\varphi }$, we find

$$\begin{eqnarray}\begin{array}{rcl}\lambda \,{\displaystyle \int }_{{b}_{\min }}^{{b}_{\max }}{u}_{0}\,\bar{{ \mathcal T }}b\,{db}\,d\bar{\varphi } & = & \pi \,\lambda [\bar{{\mathbb{T}}}({b}_{\max },{u}_{0}{| }_{{b}_{\max }})\\ & & -\,\bar{{\mathbb{T}}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})].\end{array}\end{eqnarray} \tag{ 70 }$$

Here, $\bar{{\mathbb{T}}}$ turns out to be

$$\begin{eqnarray}&&\bar{{\mathbb{T}}}(b,{u}_{0})=-2\,{\mathbb{T}}(b,{u}_{0}).\end{eqnarray} \tag{ 71 }$$

The negative sign in Equation (71) means that the tidal stretching along the z direction is clearly an antifriction force that opposes the dynamical friction force.

Finally, the relative strength of the tidal part compared to the nontidal part for motion along the z direction is given by Ψ_z = − 2Ψ_x , which we will estimate for the tidal perturbation of the Fornax dSph galaxy by the MW in Section 4.

3.3. General Treatment

In the general case, the initial fixed relative velocity vector is arbitrary; however, the freedom in a constant rotation about the z-axis may be employed to render v ₀ in the (x, z) plane with no loss in generality. As before, we can write the unperturbed motion in the form

$$\begin{eqnarray}&&{\boldsymbol{\xi }}(t)={\boldsymbol{b}}+{{\boldsymbol{v}}}_{0}\,t,\qquad {\boldsymbol{b}}\cdot {{\boldsymbol{v}}}_{0}=0,\end{eqnarray} \tag{ 72 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{v}}}_{0} & = & {v}_{0}(\sin \theta ,0,\cos \theta ),\\ {\boldsymbol{b}} & = & b(-\cos \theta \sin \phi ,\cos \phi ,\sin \theta \sin \phi ).\end{array}\end{eqnarray} \tag{ 73 }$$

That is, ${\boldsymbol{b}}=b(\cos \phi \,{{\boldsymbol{e}}}_{1}+\sin \phi \,{{\boldsymbol{e}}}_{2})$, where e ₁ = (0, 1, 0) is the unit vector in the y direction and ${{\boldsymbol{e}}}_{2}=(-\cos \theta ,0,\sin \theta )$ is the unit vector in the (x, z) plane normal to v ₀. Here, θ: 0 → π is the fixed polar angle of the initial relative velocity vector, while the azimuthal angle ϕ: 0 → 2π indicates the range of the possible directions of the impact parameter. As before, the motion takes place from t = − t₀ → t = t₀ and the freedom in time translation has been employed such that the midpoint of the unperturbed motion is at t = 0. Let us note that we work within the region bounded by the tidal radius R₀, namely,

$$\begin{eqnarray}&&{r}_{0}={\left[{\boldsymbol{\xi }}(\pm {t}_{0})\cdot {\boldsymbol{\xi }}(\pm {t}_{0})\right]}^{1/2}={\left({v}_{0}^{2}\,{t}_{0}^{2}+{b}^{2}\right)}^{1/2}\lt {R}_{0}.\end{eqnarray} \tag{ 74 }$$

The general case considered here simply reduces to the two previous cases when (θ, ϕ) are chosen appropriately. Indeed, for θ = π/2 and ϕ = φ the unperturbed motion is along the x direction with ${y}_{0}=b\cos \varphi$ and ${z}_{0}=b\sin \varphi$, while for θ = 0 and $\phi =\bar{\varphi }-\pi /2$ we recover the unperturbed motion along the z-axis with ${\bar{x}}_{0}=b\cos \bar{\varphi }$ and ${\bar{y}}_{0}=b\sin \bar{\varphi }$.

We now use the general Equations (19)–(25) and the results of Appendix B to integrate the perturbation equations. We find

$$\begin{eqnarray}\begin{array}{rcl}{\dot{{\boldsymbol{x}}}}_{1} & = & -\displaystyle \frac{M^{\prime} }{{v}_{0}{b}^{2}}\left[\displaystyle \frac{u}{U}+\displaystyle \frac{{u}_{0}}{{U}_{0}}+\displaystyle \frac{\lambda \,{b}^{2}}{M^{\prime} }(u+{u}_{0})\right]{\boldsymbol{b}}\\ & & +\,\displaystyle \frac{M^{\prime} }{{v}_{0}^{2}}\left[\displaystyle \frac{1}{U}-\displaystyle \frac{1}{{U}_{0}}-\displaystyle \frac{\lambda }{2\,M^{\prime} }({u}^{2}-{u}_{0}^{2})\right]{{\boldsymbol{v}}}_{0}+3\,\displaystyle \frac{\lambda }{{v}_{0}}\\ & & \times \,\left[(u+{u}_{0})b\sin \theta \sin \phi +\displaystyle \frac{1}{2}({u}^{2}-{u}_{0}^{2})\cos \theta \right]\,{\boldsymbol{n}},\end{array}\end{eqnarray} \tag{ 75 }$$

where n = (0, 0, 1) is the unit vector along the z-axis defined in Equation (11). Let us note that $\delta {\dot{{\boldsymbol{x}}}}_{1}={\dot{{\boldsymbol{x}}}}_{1}({t}_{0})-{\dot{{\boldsymbol{x}}}}_{1}(-{t}_{0})={\dot{{\boldsymbol{x}}}}_{1}({t}_{0})$ can be simply calculated from Equation (75) and is such that

$$\begin{eqnarray}&&\delta {\dot{{\boldsymbol{x}}}}_{1}\cdot {{\boldsymbol{v}}}_{0}=6\,\lambda \,b\,{u}_{0}\sin \theta \cos \theta \sin \phi .\end{eqnarray} \tag{ 76 }$$

This quantity is in general nonzero in the presence of tides. It does vanish, however, for some special cases including θ = 0, i.e., for motion along the z-axis, and θ = π/2, i.e., for motion along the x-axis.

Integration of Equation (75) with the appropriate initial conditions gives

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{x}}}_{1} & = & \displaystyle \frac{M^{\prime} }{{v}_{0}^{3}}\,{{\mathbb{A}}}_{1}{{\boldsymbol{v}}}_{0}-\displaystyle \frac{M^{\prime} }{{v}_{0}^{2}{b}^{2}}{{\mathbb{A}}}_{2}\,{\boldsymbol{b}}\\ & & -\,\displaystyle \frac{\lambda }{6{v}_{0}^{3}}({{\mathbb{B}}}_{1}\,{{\boldsymbol{v}}}_{0}+3{v}_{0}{{\mathbb{B}}}_{2}\,{\boldsymbol{b}}-3{v}_{0}{{\mathbb{B}}}_{3}\,{\boldsymbol{n}}),\end{array}\end{eqnarray} \tag{ 77 }$$

where ${{\mathbb{B}}}_{3}$ is given by

$$\begin{eqnarray}&&{{\mathbb{B}}}_{3}:= \cos \theta \,{{\mathbb{B}}}_{1}+3b\sin \theta \sin \phi \,{{\mathbb{B}}}_{2}.\end{eqnarray} \tag{ 78 }$$

The net change in relative velocity is given by $\delta {\boldsymbol{v}}=G\,{\dot{{\boldsymbol{x}}}}_{1}({t}_{0})+{G}^{2}\,{\dot{{\boldsymbol{x}}}}_{2}({t}_{0})$. We are interested in

$$\begin{eqnarray}&&\delta {\boldsymbol{v}}\cdot {{\boldsymbol{v}}}_{0}=6G\lambda {{bu}}_{0}\sin \theta \cos \theta \sin \phi +{G}^{2}\,{\dot{{\boldsymbol{x}}}}_{2}({t}_{0})\cdot {{\boldsymbol{v}}}_{0},\end{eqnarray} \tag{ 79 }$$

where the linear term in G is given by Equation (76). To find the equation for x ₂, we must substitute Equation (77) into Equation (25). After some algebra, the result that we need is given by

$$\begin{eqnarray}&&{\ddot{{\boldsymbol{x}}}}_{2}\cdot {{\boldsymbol{v}}}_{0}=\displaystyle \frac{M{{\prime} }^{2}}{{v}_{0}{U}_{0}{U}^{5}}\,{ \mathcal A }-\displaystyle \frac{M^{\prime} \lambda }{6{v}_{0}{U}^{5}}\,{{\mathbb{B}}}_{4}+\displaystyle \frac{M^{\prime} \lambda }{2{v}_{0}b}\,{{\mathbb{B}}}_{5}+O({\lambda }^{2}).\end{eqnarray} \tag{ 80 }$$

Here, ${ \mathcal A }$ is given by Equation (37) and

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{B}}}_{4} & := & 6{U}^{5}\,{{\mathbb{A}}}_{1}+(2{u}^{2}-{b}^{2}){{\mathbb{B}}}_{1}+9{{ub}}^{2}\,{{\mathbb{B}}}_{2}\\ & & -\,9u(u\cos \theta +b\sin \theta \sin \phi ){{\mathbb{B}}}_{3},\end{array}\end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray}&&{{\mathbb{B}}}_{5}:= [6b\cos \theta {{\mathbb{A}}}_{1}-6\sin \theta \sin \phi \,{{\mathbb{A}}}_{2}-{{bU}}^{-3}{{\mathbb{B}}}_{3}]\cos \theta .\end{eqnarray} \tag{ 82 }$$

Next, we integrate Equation (80) from t = − t₀ to t = t₀ and keep only the linear terms in λ. The result is

$$\begin{eqnarray}&&\delta {\dot{{\boldsymbol{x}}}}_{2}\cdot {{\boldsymbol{v}}}_{0}=\displaystyle \frac{2M{{\prime} }^{2}{u}_{0}}{{v}_{0}^{2}{U}_{0}^{3}}\,{ \mathcal S }+\displaystyle \frac{2M^{\prime} \,\lambda \,{u}_{0}}{{v}_{0}^{2}}\,{ \mathcal T }^{\prime} ,\end{eqnarray} \tag{ 83 }$$

where ${ \mathcal S }$ is given by Equation (41) and ${ \mathcal T }^{\prime}$ can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal T }^{\prime} & = & {P}_{2}(\cos \theta )\bar{{ \mathcal T }}-\displaystyle \frac{3}{2}\displaystyle \frac{{u}_{0}^{3}}{{U}_{0}^{3}}{\sin }^{2}\theta \cos 2\phi \\ & & +\,\displaystyle \frac{3b}{2{u}_{0}}\left[\bar{{ \mathcal T }}-\displaystyle \frac{{u}_{0}^{3}(2{u}_{0}^{2}+3{b}^{2})}{{U}_{0}^{3}{b}^{2}}\right]\sin \theta \cos \theta \sin \phi .\end{array}\end{eqnarray} \tag{ 84 }$$

Here, P₂(x) is the Legendre polynomial of the second degree (Abramowitz & Stegun 1964)

$$\begin{eqnarray}&&{P}_{2}(x):= \displaystyle \frac{1}{2}(3{x}^{2}-1).\end{eqnarray} \tag{ 85 }$$

In Equation (84), ${ \mathcal T }^{\prime}$ reduces to $\bar{{ \mathcal T }}$ for motion along the z-axis (θ = 0) and to ${ \mathcal T }$ for motion along the x-axis (θ = π/2).

Let us recall that $\delta {{\boldsymbol{v}}}_{M}=(m/M^{\prime} )\delta {\boldsymbol{v}};$ therefore,

$$\begin{eqnarray}\begin{array}{rcl}{\left(\delta {{\boldsymbol{v}}}_{M}\right)}_{| | } & = & 6{Gm}\,\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{{v}_{0}^{2}}\displaystyle \frac{\lambda \,b\,{u}_{0}}{M^{\prime} }\sin \theta \cos \theta \sin \phi \\ & & +\,2{G}^{2}m\,\displaystyle \frac{{{\boldsymbol{v}}}_{0}}{{v}_{0}^{4}}\left[\displaystyle \frac{M^{\prime} {u}_{0}}{{U}_{0}^{3}}\,{ \mathcal S }+\lambda \,{u}_{0}{ \mathcal T }^{\prime} \right].\end{array}\end{eqnarray} \tag{ 86 }$$

Next, we need to sum the contribution of all of the stars of the background medium. As before, we neglect the influence of tides on the distribution function. In summing over the flux of stars, the spatial integration is over b db ∧ d ϕ, as before. When we integrate over the azimuthal angle ϕ, terms proportional to $\sin \phi$ and $\cos 2\phi$ in Equations (84) and (86) vanish. As in the special cases we have considered above, we find

$$\begin{eqnarray}\begin{array}{rcl}M\displaystyle \frac{d{{\boldsymbol{v}}}_{M}}{{dt}} & \approx & 2{G}^{2}M\displaystyle \frac{{{\boldsymbol{v}}}_{M}}{{v}_{M}^{3}}\,\rho (\lt {v}_{M})\left[M^{\prime} \,{\mathbb{S}}(b,{u}_{0})\right.\\ & & {\left.+\,\pi \lambda \,{P}_{2}(\cos \theta )\bar{{\mathbb{T}}}(b,{u}_{0})\right]}_{{b}_{\min }}^{{b}_{\max }},\end{array}\end{eqnarray} \tag{ 87 }$$

or, finally,

$$\begin{eqnarray}\begin{array}{rcl}M\displaystyle \frac{d{{\boldsymbol{v}}}_{M}}{{dt}} & \approx & -2{G}^{2}M\displaystyle \frac{{{\boldsymbol{v}}}_{M}}{{v}_{M}^{3}}\rho (\lt {v}_{M})\left[M^{\prime} {\mathbb{S}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})\right.\\ & & \left.+\,\pi \lambda \,{P}_{2}(\cos \theta )\bar{{\mathbb{T}}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})\right],\end{array}\end{eqnarray} \tag{ 88 }$$

which reduces to our previous results for θ = π/2 and θ = 0. It is interesting to note that the tidal contribution to dynamical friction in our calculation vanishes for ${P}_{2}(\cos \theta )=0$, which occurs for ${\cos }^{2}\theta =1/3$, i.e., when the polar angle of the initial velocity is θ = θ₀ or θ = π − θ₀, where θ₀ ≈ 54.7°.

In principle, the dynamical friction force can be used to distinguish between Newtonian gravitation with particle dark matter and the Newtonian limit of nonlocal gravity theory, where the nonlocal aspect of the gravitational interaction simulates dark matter. Consequently, within the context of nonlocal gravity, there is no particle dark matter halo around galaxies. On the other hand, in the standard dark matter picture, the main contribution to dynamical friction experienced by the galactic bars comes from the dark matter halo. Therefore, it is natural to expect different amounts of dynamical friction within the galactic systems modeled in the particle dark matter scenario or nonlocal gravity. Appendix C contains a brief discussion of how Equation (88) can be extended to the Newtonian regime of nonlocal gravity.

As an application of our approach, we consider the Fornax dSph galaxy, which is one the most luminous and widely studied dSph galaxies of the Local Group. It is generally assumed that the tidal forces of the MW have not had much impact on the stellar kinematics of Fornax; however, see Hammer et al. (2019, 2020). In any case, there is no direct observational evidence for a tidal disturbance of the stellar component of Fornax (Wang et al. 2017). However, tidal stripping can still significantly affect the more extended dark matter component of Fornax. It is necessary to mention that the dark matter content of Fornax inferred from its stellar kinematics is unexpectedly low, by a factor of ∼3, compared to corresponding galaxies in cosmological simulations (Read et al. 2019). Therefore, tidal stripping of its dark matter halo has been proposed as a solution to this discrepancy (Peñarrubia et al. 2008; Genina et al. 2020).

Another important feature of Fornax that is relevant to our discussion is the puzzling spatial distribution of its globular clusters (GCs). The dark matter content of Fornax is enough to induce strong dynamical friction in the motion of GCs and the timescale of orbital decay is shorter than the age of GCs. However, there are several GCs within Fornax that are far from its center.

We can provide certain estimates for the influence of MW tides on the dynamical friction experienced by the GCs. To do so, it is first necessary to check the validity of our general approximation scheme. Our method works if the dynamical timescale within Fornax is much shorter than its orbital timescale about MW. Fornax's distance from MW is ${\mathbb{R}}=147\pm 12\,$ kpc, its total luminosity is L = (1.4 ± 0.4) × 10⁷ L_⊙ and its half-light radius is r_half = 668 ± 34 pc (Walker et al. 2009, 2010). On the other hand, the stellar mass-to-light ratio of Fornax is ϒ_* ≃ 4.6 ϒ_⊙ (Peñarrubia et al. 2008). The baryonic matter density of Fornax is well fitted by a Plummer sphere (Walker et al. 2009, 2010), namely,

$$\begin{eqnarray}&&{\rho }_{b}(r)=\displaystyle \frac{3L{{\rm{\Upsilon }}}_{* }}{4\pi {r}_{\mathrm{half}}^{3}}{\left(1+\displaystyle \frac{{r}^{2}}{{r}_{\mathrm{half}}^{2}}\right)}^{-\tfrac{5}{2}}.\end{eqnarray} \tag{ 89 }$$

Finally, the dark matter content of Fornax can be described by universal two-parameter mass profiles presented in Walker et al. (2009, 2010) for dSphs. There are two different halo models for Fornax, a cuspy NFW halo and a cored one; for the present purpose, it does not matter which one is used. We employ the NFW dark matter halo:

$$\begin{eqnarray}&&{\rho }_{d}(r)=0.368\displaystyle \frac{\,{V}_{\max }^{2}}{G\,{r}_{h}^{2}}\left(\displaystyle \frac{{r}_{h}}{r}\right){\left(1+\displaystyle \frac{r}{{r}_{h}}\right)}^{-2},\end{eqnarray} \tag{ 90 }$$

where r_h is the halo lengthscale and ${V}_{\max }$ is the maximum circular velocity associated with this halo. In the case of Fornax, these parameters are given by ${V}_{\max }=18\pm 1\,\mathrm{km}\ {{\rm{s}}}^{-1}$ and r_h = 795 pc (Walker et al. 2009, 2010). Both densities (89) and (90) can be integrated to find the corresponding masses, i.e., M_b (r) and M_d (r), enclosed within radius r. Let us define the baryonic mean density as $\langle {\rho }_{b}\rangle =3{M}_{b}({r}_{\mathrm{half}})/(4\pi {r}_{\mathrm{half}}^{3})$. Using the above-mentioned quantities we find 〈ρ_b 〉 ≃ 1.8 × 10⁻² M_⊙ pc⁻³ and we can hence calculate the internal Fornax dynamical timescale, namely,

$$\begin{eqnarray}&&{t}_{\mathrm{dyn}}=\displaystyle \frac{1}{\sqrt{G\,\langle {\rho }_{b}\rangle }}\simeq 1.14\times {10}^{8}\,\mathrm{yr}.\end{eqnarray} \tag{ 91 }$$

This is consistent with the estimate reported in Hammer et al. (2018a, 2018b) using a different method. On the other hand, the Keplerian period of Fornax in orbit about the MW is given by

$$\begin{eqnarray}&&T=2\pi \sqrt{\displaystyle \frac{{{\mathbb{R}}}^{3}}{G\,{\mathbb{M}}}}\simeq 3.87\times {10}^{9}\,\mathrm{yr},\end{eqnarray} \tag{ 92 }$$

where we have used ${\mathbb{M}}\simeq 2\times {10}^{12}\,{M}_{\odot }$ for the mass of the MW. Consequently, we have T/t_dyn ≃ 33.8. Notice that by taking into account the dark matter density in Equation (91), the ratio T/t_dyn gets even larger. There are other possible methods to estimate T, such as using the Fornax orbital parameters from Gaia DR2 data (Gaia Collaboration et al. 2018, 2020a, 2020b). However, all these estimates indicate that T/t_dyn ≫ 1 and hence our approximation scheme is reasonably valid in this case.

We are now in a position to estimate the tidal influence of the MW on the dynamical friction within Fornax. More specifically, we compute

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{x}=-\displaystyle \frac{1}{2}\,{{\rm{\Psi }}}_{z}=\displaystyle \frac{\pi \,\lambda }{M^{\prime} }\displaystyle \frac{{\mathbb{T}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})}{{\mathbb{S}}({b}_{\min },{u}_{0}{| }_{{b}_{\min }})},\end{eqnarray} \tag{ 93 }$$

on the basis of Equations (58) and (71). To do so, we take the mean value of the King model core radii of Fornax GCs as representative of ${b}_{\min }\simeq 4.11\,$ pc (Boldrini et al. 2019). Similarly, we use the mean mass of the GCs for $M^{\prime} \simeq M\simeq 1.91\times {10}^{5}\,{M}_{\odot }$. Let us recall here that the tidal radius R₀ is independent of the size of the binary system and is purely determined by the net two-body mass and the characteristics of the external source. For all interactions between a GC and the stars of the background medium, we get essentially the same tidal radius R₀. That is,

$$\begin{eqnarray}&&{\left({u}_{0}^{2}+{b}^{2}\right)}^{1/2}\lt {R}_{0}={\left(\displaystyle \frac{M^{\prime} }{{\mathbb{M}}}\right)}^{1/3}\,{\mathbb{R}},\end{eqnarray} \tag{ 94 }$$

where for Fornax, R₀ = 0.672 kpc. On the other hand, the characteristic size of the Fornax galaxy is given by D = 2 r₂₀₀, where r₂₀₀ is the radius at which ρ_b + ρ_d = 200 ρ_crit. Here, ${\rho }_{\mathrm{crit}}=3{H}_{0}^{2}/(8\pi G)$ is the cosmic critical density and H₀ is the Hubble constant, which we assume is given by H₀ = 70 km s⁻¹ Mpc⁻¹. Using Equations (89) and (90), one can numerically determine r₂₀₀; hence, we find D ≈ 18 kpc. Therefore, D > R₀ and we find from Equation (93) that ${u}_{0}{| }_{{b}_{\min }}$ is at most ≈0.672 kpc. The determination of dynamical friction force in Fornax is not a simple task without implementing an appropriate distribution function. On the other hand, the ratios Ψ_x and Ψ_z are helpful in estimating the significance of tides in our problem. Given all the parameters discussed above, these ratios are Ψ_x ≈ 0.07 and Ψ_z ≈ −0.14. Therefore, for motion in the (x, y) plane the tidal effects enhance dynamical friction, while for motion along the z-axis tidal stretching diminishes dynamical friction. For motion along an arbitrary direction with polar angle θ, the corresponding ratio is given by ${{\rm{\Psi }}}_{z}{P}_{2}(\cos \theta )$. Thus far, our preliminary analysis for Fornax implies that tidal effects would affect dynamical friction by a factor of around 10%; however, an appropriate distribution function is required to reach firm conclusions. Nonetheless, this conclusion is consistent with the observations where no sign of tidal disturbances is found in the stellar component within Fornax (Wang et al. 2017, 2019). Furthermore, our result is consistent with the N-body simulations claiming that the internal kinematics of Fornax is only mildly influenced by the tidal effects of MW (Battaglia et al. 2015).

We have presented an approximation scheme based on expansion in powers of the gravitational constant G to study the Chandrasekhar dynamical friction force in the presence of tides. The tidal interaction is considered within the distant tide linear perturbation approach. To extend our results to a system of stars, we need to determine the distribution function of the stars when tides are present. However, in this initial analytic study of dynamical friction in the presence of tides, we neglect the influence of tides on the phase space distribution of stars. The self-gravity of the stellar system is neglected as well. Depending on circumstances, tidal forces can strengthen or weaken the dynamical friction force. In connection with possible astrophysical applications of our results, we consider the influence of the Milky Way on the Fornax dwarf galaxy. In this case, preliminary estimates suggest that the effect of Galactic tides on dynamical friction within the Fornax dSph galaxy could be around 10%.

M.R. is grateful to Elena Asencio and Indranil Banik for sharing their unpublished research on the tidal stability of Fornax cluster dwarf galaxies. The work of M.R. has been supported by the Ferdowsi University of Mashhad.

We consider a Newtonian astronomical system of N bodies with inertial masses M_α , α = 1, 2,...,N, in a background Cartesian coordinate system X = (X, Y, Z) as in Figure 2. The Newtonian gravitational force on M_α due to M_β is

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\alpha \beta }=\displaystyle \frac{G\,{M}_{\alpha }\,{M}_{\beta }({{\boldsymbol{X}}}_{\beta }-{{\boldsymbol{X}}}_{\alpha })}{| {{\boldsymbol{X}}}_{\beta }-{{\boldsymbol{X}}}_{\alpha }{| }^{3}}.\end{eqnarray} \tag{ A1 }$$

The N-body system is placed in the exterior gravitational field of an external source with potential Φ( X ) so that ∇²Φ = 0 at the N-body system. The equation of motion for M_α is

$$\begin{eqnarray}&&{M}_{\alpha }\,\displaystyle \frac{{d}^{2}{X}_{\alpha }^{i}}{{{dt}}^{2}}=-{M}_{\alpha }\,{\delta }^{{ij}}{{\rm{\nabla }}}_{j}{\rm{\Phi }}({{\boldsymbol{X}}}_{\alpha })+\sum _{\beta \ne \alpha }\,{F}_{\alpha \beta }^{i}.\end{eqnarray} \tag{ A2 }$$

We define the center of mass of the N-body system in the standard manner, namely,

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{{CM}}=\displaystyle \frac{{\sum }_{\alpha }\,{M}_{\alpha }\,{{\boldsymbol{X}}}_{\alpha }}{{\sum }_{\alpha }{M}_{\alpha }}.\end{eqnarray} \tag{ A3 }$$

Zoom In Zoom Out Reset image size

Let us now introduce the approximation that the linear size of the N-body system D is much smaller than ${\mathbb{R}}$, the distance to the source; that is $D\ll {\mathbb{R}}$. Therefore, we make a first-order tidal approximation

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{i}{\rm{\Phi }}({{\boldsymbol{X}}}_{\alpha })\approx {{\rm{\nabla }}}_{i}{\rm{\Phi }}({{\boldsymbol{X}}}_{{CM}})+{{\rm{\nabla }}}_{i}\,{{\rm{\nabla }}}_{j}\,{\rm{\Phi }}({{\boldsymbol{X}}}_{{CM}})({X}_{\alpha }^{j}-{X}_{{CM}}^{j}).\end{eqnarray} \tag{ A4 }$$

It is useful to define the symmetric and traceless tidal matrix K_ij ,

$$\begin{eqnarray}&&{K}_{{ij}}=\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {X}^{i}\,\partial {X}^{j}},\qquad {K}_{{ij}}={K}_{{ji}},\qquad \mathrm{tr}({K}_{{ij}})=0.\end{eqnarray} \tag{ A5 }$$

Moreover, we define

$$\begin{eqnarray}&&{x}_{\alpha }^{i}={X}_{\alpha }^{i}-{X}_{{CM}}^{i},\end{eqnarray} \tag{ A6 }$$

so that Equation (A2) now takes the form

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\alpha }\,\displaystyle \frac{{d}^{2}{x}_{\alpha }^{i}}{{{dt}}^{2}}+{M}_{\alpha }\,\displaystyle \frac{{d}^{2}{X}_{{CM}}^{i}}{{{dt}}^{2}} & = & -{M}_{\alpha }\,{\delta }^{{ij}}{{\rm{\nabla }}}_{j}{\rm{\Phi }}({{\boldsymbol{X}}}_{{CM}})\\ & & -{M}_{\alpha }\,{K}^{i}{}_{j}({{\boldsymbol{X}}}_{{CM}}){x}_{\alpha }^{j}+\displaystyle \sum _{\beta \ne \alpha }\,{F}_{\alpha \beta }^{i}.\end{array}\end{eqnarray} \tag{ A7 }$$

If we now sum this equation over α and use ∑_α M_α x _α = 0, we find to first order in the tides

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{{\boldsymbol{X}}}_{{CM}}}{{{dt}}^{2}}=-{\rm{\nabla }}{\rm{\Phi }}({{\boldsymbol{X}}}_{{CM}}),\end{eqnarray} \tag{ A8 }$$

which describes the motion of the center of mass of the system about the source as though the tides never existed. This is, of course, true at the linear order in tidal perturbation. We note that X _CM is a function of time t, since the center of mass of the system orbits about the distant mass ${\mathbb{M}}$. Combining Equations (A7) and (A8), we find the important result that

$$\begin{eqnarray}&&{M}_{\alpha }\,\displaystyle \frac{{d}^{2}{x}_{\alpha }^{i}}{{{dt}}^{2}}=-{M}_{\alpha }\,{K}^{i}{}_{j}({{\boldsymbol{X}}}_{{CM}}){x}_{\alpha }^{j}+\sum _{\beta \ne \alpha }\,{F}_{\alpha \beta }^{i},\end{eqnarray} \tag{ A9 }$$

where in ${F}_{\alpha \beta }^{i}$, we have X _β − X _α = x _β − x _α .

Let us now derive the energy equation for the internal motions of the N-body system. To this end, we introduce v _α = d x _α /dt in Equation (A9) in the standard manner and obtain

$$\begin{eqnarray}&&{M}_{\alpha }\,{{\boldsymbol{v}}}_{\alpha }\cdot \displaystyle \frac{d{{\boldsymbol{v}}}_{\alpha }}{{dt}}=-{M}_{\alpha }\,{K}_{{ij}}({{\boldsymbol{X}}}_{{CM}}){v}_{\alpha }^{i}\,{x}_{\alpha }^{j}+\sum _{\beta \ne \alpha }\,{{\boldsymbol{v}}}_{\alpha }\cdot {{\boldsymbol{F}}}_{\alpha \beta }.\end{eqnarray} \tag{ A10 }$$

The total internal energy ${ \mathcal E }$ is the sum of internal kinetic and potential energies,

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal E } & = & { \mathcal K }+{ \mathcal W },\qquad { \mathcal K }=\displaystyle \frac{1}{2}\displaystyle \sum _{\alpha }{M}_{\alpha }\,{v}_{\alpha }^{2},\\ { \mathcal W } & = & -\displaystyle \frac{1}{2}\,\displaystyle \sum _{\alpha ,\beta }^{^{\prime} }\,\displaystyle \frac{G\,{M}_{\alpha }\,{M}_{\beta }}{| {{\boldsymbol{x}}}_{\alpha }-{{\boldsymbol{x}}}_{\beta }| },\end{array}\end{eqnarray} \tag{ A11 }$$

where a prime over the summation sign indicates that α ≠ β in the sum. Then, Equation (A10) implies

$$\begin{eqnarray}&&\displaystyle \frac{d{ \mathcal E }}{{dt}}=-{K}_{{ij}}({{\boldsymbol{X}}}_{{CM}})\sum _{\alpha }{M}_{\alpha }\,{v}_{\alpha }^{i}\,{x}_{\alpha }^{j}=-\displaystyle \frac{1}{6}{K}_{{ij}}\displaystyle \frac{{{dQ}}^{{ij}}}{{dt}},\end{eqnarray} \tag{ A12 }$$

where Q^ij is the symmetric and traceless internal quadrupole tensor of the system given by

$$\begin{eqnarray}&&{Q}^{{ij}}=\sum _{\alpha }{M}_{\alpha }(3{x}_{\alpha }^{i}\,{x}_{\alpha }^{j}-{x}_{\alpha }^{2}\,{\delta }^{{ij}}).\end{eqnarray} \tag{ A13 }$$

Let us note that we can write

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left({ \mathcal E }+\displaystyle \frac{1}{6}{K}_{{ij}}\,{Q}^{{ij}}\right)=\displaystyle \frac{1}{6}\displaystyle \frac{{{dK}}_{{ij}}}{{dt}}\,{Q}^{{ij}}.\end{eqnarray} \tag{ A14 }$$

The temporal dependence of K_ij has to do with the motion of the center of mass of the system, which is typically very slow compared to motions within the N-body system. During the fast internal motions of the system, the slow motion of the center of mass may be considered to be approximately uniform; that is, it may be a reasonable approximation in some cases to neglect the temporal dependence of the slow motion. Then, Equation (A14) implies that we have a tidally disturbed N-body system in the first order of tidal perturbation that is approximately conservative; that is, the sum total of kinetic plus potential plus tidal energies remains constant in time.

Similarly, the internal angular momentum of the system can be defined via

$$\begin{eqnarray}&&{{ \mathcal L }}_{i}={\epsilon }_{{ijk}}\sum _{\alpha }{M}_{\alpha }\,{x}_{\alpha }^{j}\,\displaystyle \frac{{{dx}}_{\alpha }^{k}}{{dt}}\end{eqnarray} \tag{ A15 }$$

and equation of motion (A9) then implies

$$\begin{eqnarray}&&\displaystyle \frac{d{{ \mathcal L }}^{i}}{{dt}}=\displaystyle \frac{1}{3}{\epsilon }^{{ijk}}{K}_{j}{}^{l}\,{Q}_{{lk}}.\end{eqnarray} \tag{ A16 }$$

For astrophysical applications, this approach can be extended to include the virial theorem (Binney & Tremaine 2008) as follows. Let us define the quantities related to the moment of inertia of the system, namely,

$$\begin{eqnarray}\begin{array}{rcl}{\mathbb{I}} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{\alpha }{M}_{\alpha }\,{x}_{\alpha }^{2},\quad \displaystyle \frac{d\,{\mathbb{I}}}{{dt}}=\displaystyle \sum _{\alpha }{M}_{\alpha }\,{{\boldsymbol{x}}}_{\alpha }\cdot {{\boldsymbol{v}}}_{\alpha },\\ \displaystyle \frac{{d}^{2}{\mathbb{I}}}{{{dt}}^{2}} & = & 2{ \mathcal K }+\displaystyle \sum _{\alpha }{M}_{\alpha }\,{{\boldsymbol{x}}}_{\alpha }\cdot \displaystyle \frac{d{{\boldsymbol{v}}}_{\alpha }}{{dt}}.\end{array}\end{eqnarray} \tag{ A17 }$$

Then, from Equation (A9) we get

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\alpha }\,{{\boldsymbol{x}}}_{\alpha }\cdot \displaystyle \frac{d{{\boldsymbol{v}}}_{\alpha }}{{dt}} & = & -{K}_{{ij}}({{\boldsymbol{X}}}_{{CM}}(t)){M}_{\alpha }\,{x}_{\alpha }^{i}\,{x}_{\alpha }^{j}\\ & & +\,\displaystyle \sum _{\beta \ne \alpha }\,\displaystyle \frac{G\,{M}_{\alpha }\,{M}_{\beta }({{\boldsymbol{x}}}_{\beta }-{{\boldsymbol{x}}}_{\alpha })\cdot {{\boldsymbol{x}}}_{\alpha }}{| {{\boldsymbol{x}}}_{\beta }-{{\boldsymbol{x}}}_{\alpha }{| }^{3}}.\end{array}\end{eqnarray} \tag{ A18 }$$

We sum this equation over α and write the sum again with α and β exchanged. Adding the resulting equations we finally get in the standard manner

$$\begin{eqnarray}&&\sum _{\alpha }{M}_{\alpha }\,{{\boldsymbol{x}}}_{\alpha }\cdot \displaystyle \frac{d{{\boldsymbol{v}}}_{\alpha }}{{dt}}={ \mathcal W }-\displaystyle \frac{1}{3}\,{K}_{{ij}}({{\boldsymbol{X}}}_{{CM}}(t)){Q}^{{ij}},\end{eqnarray} \tag{ A19 }$$

which is equal to $\ddot{{\mathbb{I}}}-2{ \mathcal K }$ by Equation (A17). Finally, assuming the "fast" system relaxes over timescales short compared to the timescale of the "slow" center-of-mass motion, we can average our result over time. Assuming that the average of $\ddot{{\mathbb{I}}}$ over time vanishes (Landau & Lifshitz 1988), we finally have the result

$$\begin{eqnarray}&&2\,\langle { \mathcal K }\rangle +\langle { \mathcal W }\rangle =\displaystyle \frac{1}{3}\,{K}_{{ij}}\,\langle {Q}^{{ij}}\rangle .\end{eqnarray} \tag{ A20 }$$

The integration of equations of motion in Section 3 is simplified using the indefinite integrals given below. Let us assume $F={({x}^{2}+{b}^{2})}^{1/2};$ then,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \int \displaystyle \frac{{dx}}{F} & = & \mathrm{ln}(F+x),\\ \displaystyle \int \mathrm{ln}(F+x){dx} & = & x\,\mathrm{ln}(F+x)-F,\end{array}\end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray}&&\int \displaystyle \frac{{dx}}{{F}^{3}}=\displaystyle \frac{x}{{b}^{2}\,F},\ \int \displaystyle \frac{{dx}}{{F}^{5}}=\displaystyle \frac{1}{{b}^{4}{F}^{3}}\,x\left(\displaystyle \frac{2}{3}{x}^{2}+{b}^{2}\right),\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&\int \displaystyle \frac{{xdx}}{F}=F,\ \int \displaystyle \frac{{xdx}}{{F}^{3}}=-\displaystyle \frac{1}{F},\qquad \int \displaystyle \frac{{x}^{2}{dx}}{{F}^{5}}=\displaystyle \frac{{x}^{3}}{3{b}^{2}{F}^{3}},\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&\int (\mathrm{ln}x){xdx}=\displaystyle \frac{1}{2}{x}^{2}\mathrm{ln}x-\displaystyle \frac{1}{4}{x}^{2},\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}&&\int \mathrm{ln}(F-x)b\,{db}=\displaystyle \frac{1}{2}{b}^{2}\mathrm{ln}(F-x)-\displaystyle \frac{1}{2}x\,F-\displaystyle \frac{1}{4}{b}^{2},\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&\int \displaystyle \frac{(\mathrm{ln}b)b\,{db}}{{F}^{3}}=-\displaystyle \frac{\mathrm{ln}b}{F}+\displaystyle \frac{1}{2x}\,\mathrm{ln}\left(\displaystyle \frac{F-x}{F+x}\right),\end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray}&&\int \displaystyle \frac{\mathrm{ln}(F+x)}{{F}^{3}}b\,{db}=-\displaystyle \frac{\mathrm{ln}(F+x)}{F}+\displaystyle \frac{1}{x}[\mathrm{ln}F-\mathrm{ln}(F+x)].\end{eqnarray} \tag{ B7 }$$

We also need to evaluate the definite integral

$$\begin{eqnarray}&&I={\int }_{-{x}_{0}}^{{x}_{0}}\displaystyle \frac{\mathrm{ln}(F+x)}{{F}^{3}}\,{dx}={\int }_{-{x}_{0}}^{{x}_{0}}\displaystyle \frac{\mathrm{ln}(F-x)}{{F}^{3}}\,{dx}.\end{eqnarray} \tag{ B8 }$$

From $\mathrm{ln}(F+x)+\mathrm{ln}(F-x)=2\mathrm{ln}b$, we conclude

$$\begin{eqnarray}&&I=\mathrm{ln}b\,{\int }_{-{x}_{0}}^{{x}_{0}}\displaystyle \frac{{dx}}{{F}^{3}}=2\displaystyle \frac{\mathrm{ln}b}{{b}^{2}}\displaystyle \frac{{x}_{0}}{{F}_{0}},\end{eqnarray} \tag{ B9 }$$

where F₀ = F(x₀, b). Similarly,

$$\begin{eqnarray}&&{\int }_{-{x}_{0}}^{{x}_{0}}\displaystyle \frac{\mathrm{ln}(F+x)}{{F}^{5}}\,{dx}=2\displaystyle \frac{\mathrm{ln}b}{{b}^{4}}\displaystyle \frac{{x}_{0}}{{F}_{0}^{3}}\left(\displaystyle \frac{2}{3}{x}_{0}^{2}+{b}^{2}\right).\end{eqnarray} \tag{ B10 }$$

The corresponding definite integrals can be easily evaluated via integration by parts using Equation (B2).

In a previous paper (Roshan & Mashhoon 2021), we extended Chandrasekhar's formula for dynamical friction to the Newtonian regime of nonlocal gravity theory and briefly studied its implications for barred spiral galaxies. Nonlocal gravity (NLG) is a classical nonlocal generalization of Einstein's theory of gravitation that takes the past history of the gravitational field into account. The gravitational field in NLG is local but satisfies partial integro-differential field equations. Moreover, NLG has been constructed in close analogy with the nonlocal electrodynamics of media. It turns out that such a classical nonlocal aspect of the gravitational interaction simulates dark matter. A detailed description of NLG is contained in Mashhoon (2017). It seems worthwhile to indicate briefly how the formal results of the present work can carry over to NLG. To this end, we must first describe the analog of the Newtonian inverse-square law in NLG.

In the Newtonian limit, NLG in effect involves the standard nonrelativistic gravitational force

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{NLG}}({\boldsymbol{x}})=-m{\rm{\nabla }}{{\rm{\Phi }}}_{\mathrm{NLG}}({\boldsymbol{x}})\end{eqnarray} \tag{ C1 }$$

on a test particle of inertial mass m in the gravitational potential Φ_NLG( x ), which satisfies the nonlocal Poisson equation that can be expressed in the form

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\nabla }}}^{2}{{\rm{\Phi }}}_{\mathrm{NLG}} & = & 4\pi G(\rho +{\rho }_{D}),\\ {\rho }_{D}({\boldsymbol{x}}) & = & \displaystyle \int q({\boldsymbol{x}}-{\boldsymbol{y}})\rho ({\boldsymbol{y}}){d}^{3}y.\end{array}\end{eqnarray} \tag{ C2 }$$

Here, ρ is the density of matter and ρ_D is the density of effective dark matter in NLG. In this theory, what appears as dark matter in astrophysics and cosmology is in reality the nonlocal aspect of gravity itself and its density is given by the convolution of a certain reciprocal kernel q with the density of matter ρ. The field equations of nonlocal gravity in the Newtonian regime of the theory reduce to the nonlocal Poisson equation (C2) provided the functions involved be smooth and satisfy certain reasonable mathematical properties (Mashhoon 2017). Moreover, the reciprocal kernel q must be determined on the basis of observational data.

To simplify matters, let us assume that q( x − y ) is spherically symmetric; then, for q(r), r = ∣ x − y ∣, two possible forms have been discussed in detail. These are (Mashhoon 2017)

$$\begin{eqnarray}\begin{array}{rcl}{q}_{1}(r) & = & \displaystyle \frac{1}{4\pi {\lambda }_{0}}\,\displaystyle \frac{1+{\mu }_{0}\,({a}_{0}+r)}{r({a}_{0}+r)}\,{e}^{-{\mu }_{0}\,r},\\ {q}_{2}(r) & = & \displaystyle \frac{r}{{a}_{0}+r}\,{q}_{1}(r),\end{array}\end{eqnarray} \tag{ C3 }$$

which contain three constant parameters. The basic NLG lengthscale is determined by λ₀, since nonlocality disappears when λ₀ tends to infinity; moreover, a₀ moderates the short distance behavior of the kernel and μ₀ is the "Yukawa" parameter. It proves useful to define a dimensionless parameter α₀ = 2/(λ₀ μ₀). Solar system data provide a lower bound for a₀, namely, a₀ > 10¹⁴ cm (Chicone & Mashhoon 2016). With a₀ = 0, q₁ = q₂ and the rotation curves of nearby spiral galaxies can be used to find μ₀ and λ₀. Indeed, observational data regarding nearby spiral galaxies and clusters of galaxies are consistent with (Rahvar & Mashhoon 2014)

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{0} & = & 10.94\pm 2.56,\quad {\mu }_{0}=0.059\pm 0.028\,{\mathrm{kpc}}^{-1},\\ {\lambda }_{0} & = & \displaystyle \frac{2}{{\alpha }_{0}\,{\mu }_{0}}\approx 3\pm 2\,\mathrm{kpc}.\end{array}\end{eqnarray} \tag{ C4 }$$

It is straightforward to calculate the force of gravity according to NLG on a point mass m due to another point mass $m^{\prime}$ at position r . The result is a modification of Newton's law of universal gravitational attraction given by

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{NLG}}=\displaystyle \frac{{Gmm}^{\prime} {\boldsymbol{r}}}{{r}^{3}}[1+{\rm{\Delta }}(r)].\end{eqnarray} \tag{ C5 }$$

A detailed physical interpretation of this result is contained in Roshan & Mashhoon (2021). In Equation (C5), the net contribution of the effective dark matter is contained in Δ(r) ≥ 0. The function Δ(r) starts from zero at r = 0, increases monotonically with increasing r and approaches α₀ w asymptotically as r → ∞. Here, w is a positive constant such that w = 1 for a₀ = 0 and for a₀ > 0, w depends upon whether the reciprocal kernel is chosen to be q₁ or q₂; however, for reasonable values of a₀ less than a few parsecs, w is very close to unity. Henceforth, we ignore the deviation of w from unity. That is, NLG in the Newtonian regime is such that the magnitude of the force of gravity asymptotically approaches Newtonian gravity except that now the effective constant of gravitation is about an order of magnitude larger than the standard Newtonian constant of gravitation G.

With these introductory remarks about F _NLG, we are now in a position to indicate how our Newtonian approach in this paper is affected by the presence of Δ(r). In Equations (3)–(4) and (10)–(12), the internal forces get multiplied by [1 + Δ(r)], while the external forces get multiplied by $[1+{\rm{\Delta }}({\mathbb{R}})]$. In the distant tide approximation, if ${\mathbb{R}}\gg {\mu }_{0}^{-1}$, then ${\rm{\Delta }}({\mathbb{R}})\approx {\alpha }_{0}$. In this case, the tidal radius is smaller than the Newtonian R₀ given by Equation (15). Extending these considerations to dynamical friction, it has been shown in Roshan & Mashhoon (2021) that the nontidal "Chandrasekhar" dynamical friction force gets multiplied by ${(1+\hat{\eta })}^{2}$, where $\hat{\eta }$, $0\lt \hat{\eta }\lt {\alpha }_{0}$, is a constant that depends on the parameters of the stellar system and has been discussed in detail in Section 4.1 of Roshan & Mashhoon (2021). In the present case, we find that in Equation (88) the nontidal dynamical friction force gets multiplied by ${(1+\hat{\eta })}^{2}$ as before, while the tidal part of dynamical friction force gets multiplied by $(1+\hat{\eta })(1+{\alpha }_{0})$. Finally, the matter density in Equation (88) for dynamical friction would then simply refer to the baryonic density in this case, as there is no actual dark matter in NLG.