© The Authors 2024

1. Introduction

Modified Newtonian dynamics (hereafter MOND, Milgrom 1983) is an alternative theory of (classical) gravity that was originally introduced to solve the missing mass problem on astrophysical and cosmological scales without resorting to dark matter theory (hereafter DM). In a Lagrangian formulation (see Bekenstein & Milgrom 1984), MOND can be amended to a substitution of the Poisson equation relating the gravitational potential, Φ, and mass density, ρ, with the non-linear field equation:

$$\begin{array}{c}\hfill \mathbf{\nabla }·\left[\mu \left(\frac{||\mathbf{\nabla }\mathrm{\Phi }||}{a_0}\right)\mathbf{\nabla }\mathrm{\Phi }\right]=4\pi G\rho ,\end{array}$$(1)

where the acceleration, scaling as a₀ ≈ 1.2 × 10⁻²ms⁻², is a new universal constant; whereas the interpolating function, μ, (which is, in principle, unknown) is monotonic with the asymptotic behaviour, as follows:

$$\begin{array}{c}\hfill \mu (x){\to }^{x\gg 1}1,\mu (x){\sim }^{x\ll 1}x.\end{array}$$(2)

A widely adopted form used for μ is:

$$\begin{array}{c}\hfill \mu (x)=\frac{x}{1+x}.\end{array}$$(3)

According to the equations above, for ∇Φ ≫ a₀, we can recover the Newtonian limit and Eq. (1) reduces to the Poisson equation. Vice versa, for ∇Φ ≪ a₀ the system is in the so-called deep-MOND limit (hereafter, dMOND) and Eq. (1) becomes the p-Laplace equation:

$$\begin{array}{c}\hfill \mathbf{\nabla }·\left(||\mathbf{\nabla }\mathrm{\Phi }||\mathbf{\nabla }\mathrm{\Phi }\right)=4\pi a_{0}G\rho .\end{array}$$(4)

For a given mass density, ρ, the right-hand-side of the Poisson equation and Eq. (1) are the same. We can therefore eliminate ρ, thereby obtaining the relation:

$$\begin{array}{c}\hfill \mu \left(\frac{||{\mathit{g}}_M||}{a_0}\right){\mathit{g}}_M={\mathit{g}}_N+\mathit{S},\end{array}$$(5)

between the MOND and Newtonian force fields g_M and g_N, where S ≡ ∇ × h(ρ) is a density-dependent solenoidal field that zeros-out for systems in spherical, cylindrical, or planar symmetry, while it is typically non-zero for more general density profiles.

MOND has been rather successful in reproducing the kinematics of galaxies without the need for DM to be included (Bugg 2015; Zhu et al. 2023). In fact, it has proven to be a valid alternative to the ΛCDM paradigm where the latter is challenged, across a broad range of gravitational problems, such as the structural properties of dwarf galaxies (Milgrom 1994; Sánchez-Salcedo & Hernandez 2007; McGaugh & Milgrom 2013; Asencio et al. 2022); radial acceleration relation of galaxies (RAR, Lelli et al. 2017; Mistele et al. 2024b,a); profiles of low-surface-brightness galaxies (Sanders 2021); dynamics of wide binary stars (Hernandez et al. 2012, 2023; Chae 2024); evolution of globular clusters of ultra-diffuse galaxies (Bílek et al. 2021); and the orbital velocity of interacting galaxy pairs (Scarpa et al. 2022).

MOND has been extensively investigated numerically (Brada & Milgrom 1999; Knebe & Gibson 2004; Tiret & Combes 2007; Nipoti et al. 2007a,b,c; Sanders 2008; Wu et al. 2009; Malekjani et al. 2009; Roshan et al. 2021; Kroupa et al. 2022; Banik et al. 2022; Nagesh et al. 2023; Re & Di Cintio 2023; Sollima et al. 2023), for collisionless mean field processes; for example: violent relaxation and phase-mixing in monolithic collapse, galaxy merging, and the vertical dynamics of disk galaxies. Unfortunately, due to the non-linearity of the theory (cf. Eq. 1) and, hence, the absence of the superposition principle, there is much less information available on the collisional processes (see e.g. Bílek et al. 2021). For example, Ciotti & Binney (2004) used an approach based on the fluctuations of uniform field to estimate the MOND two-body relaxation (t_2b, Chandrasekhar 1941a,b) and dynamical friction (t_DF, Chandrasekhar 1943a,b,c) time scales in the dMOND limit. The outcome of this study (see also the numerical simulations in Nipoti et al. 2008) shows that a test mass, M, crossing a purely baryonic system would feel a stronger dynamical friction (hereafter, DF) force due to the encounters with the system’s stars in MOND than in the parent

Newtonian system without DM; however, it would undergo an only slightly more efficient DF than in the equivalent Newtonian system (ENS). The is comprised of the baryonic plus DM system constructed such that the potential is the same as the purely baryonic MOND model. The Newtonian and MOND DF time scales $t_{\mathrm{DF}}^{\mathrm{N}}$ and $t_{\mathrm{DF}}^{\mathrm{M}}$, are related as:

$$\begin{array}{c}\hfill t_{\mathrm{DF}}^{\mathrm{M}}=\frac{\sqrt{2}}{1+\mathcal{R}}{t}_{\mathrm{DF}}^{\mathrm{N}},\end{array}$$(6)

where ℛ is the ratio of the amount of DM to baryons in the ENS.

In MOND, the collisional direct N-body simulations are intrinsically unfeasible; thus, results such as those given by Eq. (6) cannot be exhaustively explored via numerical experiments. The reason for this is the absence of a MONDian expression for the force exerted by a point-like particle. To address this issue, Milgrom (1986) proposed the approximated expression for the force exchanged between two masses m₁ and m₂ placed at distance, r,

$$\begin{array}{c}\hfill F_{1,2}\approx \frac{m_1{m}_2}{\sqrt{m_1+m_2}}\frac{\sqrt{G{a}_0}}{r}.\end{array}$$(7)

The latter is however valid only in the far field limit, thus making it unusable in a direct N-body code, as it lacks a regime bridging to the Newtonian 1/r² limit.

More recently, Milgrom (2010) formulated a quasi-linear MOND theory (hereafter QuMOND), where the governing field equation is:

$$\begin{array}{c}\hfill \mathbf{\nabla }·\left[\nu \left(\frac{||\mathbf{\nabla }{\mathrm{\Phi }}^N||}{a_0}\right)\mathbf{\nabla }{\mathrm{\Phi }}^N\right]=\mathrm{\Delta }{\mathrm{\Phi }}^M,\end{array}$$(8)

formally identical to Eq. (1) where now, Φ^N is the Newtonian potential generated by the density ρ through the usual Poisson equation:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Phi }}^N=-\mathbf{\nabla }·{\mathit{g}}^N=4\pi G\rho ,\end{array}$$(9)

and the new interpolating function ν(y) is related to μ(x) via:

$$\begin{array}{c}\hfill \{\begin{array}{c}x=y\nu (y)\hfill \\ y=x\mu (x)\hfill \end{array}\Rightarrow \mu (x)=\frac{y}{x}=\frac{1}{\nu (y)}.\end{array}$$(10)

The asymptotic behaviours for ν in the Newtonian and in the dMOND limits are

$$\begin{array}{c}\hfill \nu (y){\to }^{y\gg 1}1,\nu (y){\sim }^{y\ll 1}{y}^{-1/2}.\end{array}$$(11)

We note that Eq. (8) can be rearranged as

$$\begin{array}{c}\hfill \stackrel{\sim }{\rho }=-\frac{1}{4\pi G}\mathrm{\nabla }·[\nu (g^N/a_0){\mathit{g}}^N];{\mathit{g}}^N=-\mathrm{\nabla }{\mathrm{\Phi }}^N,\end{array}$$(12)

that is, the QuMOND potential satisfies the Poisson equation for the auxiliary density $\stackrel{\sim }{\rho }$. In practice, in QuMOND, we have to evaluate the Newtonian potential Φ^N generated by ρ, evaluate its gradient and then via a non-linear algebraic step obtain the auxiliary MOND density $\stackrel{\sim }{\rho }$ via Eq. (12) that acts as source for the potential Φ^M.

At variance with the usual MOND formulation, in its quasi-linear formulation the interpolating function ν has a stronger effect on the form of the gravitational potential than its counterpart μ in Eq. (1). For the usual choice (3), we obtain:

$$\begin{array}{c}\hfill \nu (y)=\frac{1}{2}+\frac{1}{2}\sqrt{1+\frac{4}{y}}.\end{array}$$(13)

We note that QuMOND contains a sort of superposition principle for the auxiliary densities $\stackrel{\sim }{\rho }$. It is therefore tempting to try to derive an expression for the DF force in this context, making use of purely kinetic arguments as in the original derivation. In this paper, we first discuss a simple expression of the QuMOND DF summing up the contributions of the effective densities for point-particles of mass, m, interacting with a test mass, with some simple choices for ν. Moreover, we extend the mean field formulation of DF, originally derived in Newtonian gravity by Kandrup (1983) to the case of QuMOND and discuss the results of simple numerical experiments.

2. Dynamical friction: Classical approach

2.1. Newtonian case

Before tackling the QuMOND DF formula, we start by briefly recalling the ingredients of the derivation of its Newtonian expression (for a more detailed review, see e.g. Binney & Tremaine 2008 or Ciotti 2021). We consider a test mass, M, travelling at v_M through an infinite background medium of particles of a mass, m, and (by now constant) number density, n, with velocity distribution function, f(v). In each encounter, M and a background particle exchange a force expressed as F ≈ GMm/b², acting for an interval of time, Δt ≈ b/v_M, where b is the impact parameter of the dynamical collision. The associated variation of momentum is Δp ≈ FΔt ≈ GMm/v_Mb. The rate of collisions of M is then

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{enc}}=\frac{2\pi nb\mathrm{d}b\mathrm{d}x}{\mathrm{d}t}=2\pi n{v}_{M}b\mathrm{d}b.\end{array}$$(14)

From this, the velocity diffusion coefficient along the trajectory of M can be written as:

$$\begin{array}{c}\hfill D_v={\displaystyle \int \mathrm{\Delta }}{p}^2\frac{\mathrm{d}{N}_{\mathrm{enc}}}{\mathrm{d}t}=2\pi G^2{M}^2\frac{m^{2}n}{v_M}{\displaystyle \int \frac{\mathrm{d}b}{b}}.\end{array}$$(15)

In the equation above, the integral in the impact parameter b diverges for both b → 0 and b → ∞ (infrared and ultraviolet divergence, respectively). It is however possible to impose suitable cut-offs b_min and b_max. This factor, after integration takes the form of the Coulomb logarithm (in analogy with the same quantity in plasma physics) lnΛ, where Λ = b_max/b_min.

Let us now assume that f(v), defined such that

$$\begin{array}{c}\hfill {\displaystyle \int f}(\mathit{x},\mathit{v})d^3\mathit{x}{d}^3\mathit{v}=N\end{array}$$(16)

is the total number of particles, is position independent and isotropic for all velocities. For such a system, the DF force per unit mass acting on M is given as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}& =-16{\pi }^2{G}^{2}n(M+m)mln\mathrm{\Lambda }\frac{\mathrm{\Psi }(v_M)}{v_M^3}{\mathit{v}}_M,\hfill \end{array}$$(17)

where

$$\begin{array}{c}\hfill \mathrm{\Psi }(v_M)={\displaystyle {\int }_0^{v_M}}{v}^{2}f(v)\mathrm{d}v\end{array}$$(18)

is the so-called fractional velocity volume function (Ciotti 2021). Therefore, in the assumption of velocity isotropy, only particles moving slower than M contribute to its slowing down. Moreover for M ≫ m, in Eq. (17) the factor m(M + m)n becomes Mρ, where ρ = nm is the background mass density. In Newtonian systems with DM, two species of particles exert DF on the test mass M. Following Ciotti (2010), Eq. (17) can be rewritten for a model with a mass spectrum as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}& =-16{\pi }^2{G}^{2}n(M+⟨m⟩)⟨m⟩ln\overline{\mathrm{\Lambda }}\frac{\stackrel{\sim }{\mathrm{\Psi }}(v_M)}{v_M^3}{\mathit{v}}_M,,\hfill \end{array}$$(19)

where $ln\overline{\mathrm{\Lambda }}$ is a velocity averaged Coulomb logarithm, ⟨m⟩ is the mean mass of the background particles, and $\stackrel{\sim }{\mathrm{\Psi }}({\mathit{v}}_M)$ is the total fractional velocity volume function¹ of the system. For the case of interest of a binary mass spectrum (i.e. stars and DM particles), we have

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Psi }}(v_M)=\frac{\mathcal{M}{H}_1+H_2}{\mathcal{M}+1},\end{array}$$(20)

where ℳ = M/⟨m⟩ and

$$\begin{array}{cc}\hfill H_1(v_M)& =\frac{{\mathrm{\Psi }}_{\ast }(v_M)+xy{\mathrm{\Psi }}_{\mathrm{DM}}(v_M)}{1+xy};\hfill \\ \hfill H_2(v_M)& =(1+x)\frac{{\mathrm{\Psi }}_{\ast }(v_M)+x{y}^2{\mathrm{\Psi }}_{\mathrm{DM}}(v_M)}{{(1+xy)}^2}.\hfill \end{array}$$(21)

In the expressions above, x = n_DM/n_* and y = m_DM/m_* are the ratios of the number densities and masses of stars and DM particles, respectively and the Ψ_DM and Ψ_* their fractional velocity volume functions. Stars and DM are not at energy equipartition, as the thermalization time scale would exceed the age of the universe of several order of magnitude due to the small mass of DM particles. Therefore, both species have the same velocity distribution dictated by their collective potential, hence Ψ_DM = Ψ_* = Ψ_DM. As the DM is constituted by unknown elementary particles, with typical masses of the order of a few GeV, we have m_DM ≪ m_* and n_DM ≫ n_*. The average mass ⟨m⟩ is expressed as:

$$\begin{array}{c}\hfill \frac{1+xy}{x}{m}_{\ast }\ll m_{\ast },\end{array}$$(22)

while

$$\begin{array}{c}\hfill H_1(v_M)=\mathrm{\Psi }(v_M);H_2(v_M)=\frac{x}{{(1+xy)}^2}\mathrm{\Psi }(v_M),\end{array}$$(23)

so that the total fractional velocity volume functions becomes:

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Psi }}(v_M)=H_1+\frac{⟨m⟩}{M}{H}_2=\mathrm{\Psi }(v_M)+\frac{m_{\ast }}{M}\frac{1}{1+xy}\mathrm{\Psi }(v_M).\end{array}$$(24)

For systems where the test mass M ≫ m_* the latter further simplifies and $\stackrel{\sim }{\mathrm{\Psi }}({\mathit{v}}_M)=\mathrm{\Psi }({\mathit{v}}_M)$ and the DF formula can be simplified to:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}& =-16{\pi }^2{G}^2{\rho }_{\mathrm{Tot}}Mln\overline{\mathrm{\Lambda }}\frac{\mathrm{\Psi }(v_M)}{v_M^3}{\mathit{v}}_M.\hfill \end{array}$$(25)

In practice, the presence of a DM distribution has enhances the DF of a factor, ρ_Tot/ρ_*, with respect to the baryonic system with the same phase-space distribution. We note that Sánchez-Salcedo et al. (2006) also derived a simplified two-component DF formula for the cases of circular orbits in DM plus baryons King (1962) cores, concluding that DF is reduced by the presence of DM of a factor (1 + ℛ)⁻¹ with respect to the case where the (same) mass is in one single component, following the light distribution.

2.2. QuMOND case

We first present a naive general expression of the MONDian DF, making use of the fact that in QuMOND, we have:

$$\begin{array}{c}\hfill F=\nu \left(\frac{\mathit{Gm}}{a_0{b}^2}\right)\frac{\mathit{GMm}}{b^2}.\end{array}$$(26)

For each particle of the background system, we defined the so-called MOND radius as

$$\begin{array}{c}\hfill r_{\mathrm{M}}=\sqrt{\frac{\mathit{Gm}}{a_0}}.\end{array}$$(27)

With the same arguments of Sect. 2.1, we now find that the momentum variation of M is:

$$\begin{array}{c}\hfill \mathrm{\Delta }p=\nu (r_{\mathrm{M}}^2/b^2)GMm/v_{M}b.\end{array}$$(28)

The form of the collision rate factor remains unchanged, as it is a geometric quantity. The diffusion coefficient D_v retains the same definition as in Eq. (15), but with a different integral over the impact parameter:

$$\begin{array}{c}\hfill I_b={\displaystyle {\int }_{b_{\min }}^{b_{\max }}}\nu {(r_{\mathrm{M}}^2/b^2)}^{2}db/b.\end{array}$$(29)

After substituting the expression in the DF formula, we obtain:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}=-16{\pi }^2{I}_b{G}^2\rho (M+m)\frac{\mathrm{\Psi }(v_M)}{v_M^3}{\mathit{v}}_M.\end{array}$$(30)

From Eq. (30) above, we notice immediately that the QuMOND DF is always larger than its Newtonian equivalent for the same baryonic mass distribution. This happens because ν² > 1, thus returning a larger factor within the integral over the impact parameters. In addition, this is also in agreement with the fact that in MOND theories the far field behaviour of a given mass distribution falls off less rapidly than in Newtonian gravity, as implied by the logarithmic trend of the MOND potential.

We now evaluate Eq. (30) for a simple choice of ν. Knowing that ν has the asymptotic trends given in Eq. (11) we start with a simplified form:

$$\begin{array}{c}\hfill \nu (y)=1+1/\sqrt{y},\end{array}$$(31)

which is essentially the sum of the Newtonian and the deep MOND régimes. In terms of μ, Eq. (31) corresponds to $\mu (x)=1+(1-\sqrt{1+4x})/2x$. For such a choice, the correction factor is $\nu (r_{\text{M}}^2/b^2)=1+b/r_{\text{M}}$, and the integral on the impact parameters becomes:

$$\begin{array}{cc}\hfill I_b& =ln\frac{b_{\max }}{b_{\min }}+2\frac{b_{\max }-b_{\min }}{r_{\mathrm{M}}}+\frac{b_{\max }^2-b_{\min }^2}{2r_{\mathrm{M}}^2}=\hfill \\ \hfill & \cong ln\mathrm{\Lambda }+2\frac{b_{\max }}{r_{\mathrm{M}}}+\frac{b_{\max }^2}{2r_{\mathrm{M}}^2};\hfill \end{array}$$(32)

where we assumed b_min ≪ r_M in the last equality. We note that the maximum impact parameter b_max in the context of Newtonian gravity is usually assumed to be some (arbitrary) scale distance² of the system at hand. In the original derivation of Chandrasekhar (1943a), where the system is infinite, the definition of b_max remains somewhat unclear. Some authors (see e.g. Van Albada & Szomoru 2020 and references therein) typically take b_max of the order of the average inter-particle distance. Here, we follow the same assumption, noting that for a star of a solar mass, the MOND radius r_M is roughly 10³ AU; this is much smaller than the typical scale radius of a galaxy, but comparable to the mean inter-particle distance among stars. This established, the QuMOND DF then becomes

$$\begin{array}{cc}\hfill \left(\frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}\right)\approx & -16{\pi }^2(ln\mathrm{\Lambda }+\frac{b_{\max }^2}{2r_{\mathrm{M}}^2}+2\frac{b_{\max }}{r_{\mathrm{M}}})\times \hfill \\ \hfill & \times G^2\rho (M+m)\frac{\mathrm{\Psi }(v_M)}{v_M^3}{\mathit{v}}_M.\hfill \end{array}$$(33)

In practice, at first order in QuMOND, each particle behaves as a point-like source exerting the usual Newtonian 1/r potential plus an infinitely extended ‘phantom DM’ halo whose contribution depends on the specific form of ν. Following this approach, the calculations as shown in Appendix A leads to the same expression of Eq. (33). We note that Eq. (33) differs from what we would obtain using the simple derivation of t_2b in Sect. 2 of Ciotti & Binney 2004; namely, this derivation is proportional to the crossing time $t_{\mathrm{cross}}=1/\sqrt{G\rho }$ to evaluate the dynamical friction time scale with the Spitzer (1987) relation t_DF = 2t_2bm/(m + M). This is because their expression of t_2b makes use of the force (7), implicitly assuming the dMOND regime.

Unlike the usual MOND interpolating function μ(x), the form of the QuMOND function ν(y) affects the behaviour of the MONDian force field. Adopting the commonly used Eq. (13) (see Kroupa et al. 2022), the DF expression becomes:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}\approx & -16{\pi }^2(ln\mathrm{\Lambda }+\frac{b_{\max }^2}{2r_{\mathrm{M}}^2}+\frac{b_{\max }}{r_{\mathrm{M}}}-\frac{1}{2}ln\frac{b_{\max }}{r_{\mathrm{M}}}-\frac{2+ln2}{4})\times \hfill \\ \hfill & \times G^2\rho (M+m)\frac{\mathrm{\Psi }(v_M)}{v_M^3}{\mathit{v}}_M.\hfill \end{array}$$(34)

The explicit derivation is given in Appendix B.

Noteworthy, we observe in both cases that the MOND correction with respect to the Newtonian baryon only case is proportional to a function of the only b_max/r_M

$$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}\right)}_{\mathrm{QuMOND}}& -{\left(\frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}\right)}_{\mathrm{Newt}}\propto \frac{b_{\max }^2}{2r_{\mathrm{M}}^2}+2\frac{b_{\max }}{r_{\mathrm{M}}},\hfill \\ \hfill & \frac{b_{\max }^2}{2r_{\mathrm{M}}^2}+\frac{b_{\max }}{r_{\mathrm{M}}}-\frac{1}{2}ln\frac{b_{\max }}{r_{\mathrm{M}}}-\frac{1}{2}.\hfill \end{array}$$(35)

The same term $b_{\text{max}}^2/2r_{\text{M}}^2$ always dominates the MOND correction of DF, regardless of the exact form of ν, while the smaller terms are affected by the particular choice. Recently, a number of authors (see e.g. Lelli et al. 2017) proposed an expression for ν supported by observational data proportional to the form of radial acceleration relation. However, its functional form is rather complex thus making the explicit evaluation of Eq. (30) particularly cumbersome.

As an example, in the top panel of Figure 1 we show for a test mass M = 10³ M_⊙ the ratio η_QuMOND/η_N of the QuMOND and Newtonian dynamical friction coefficient as a function of the stellar density, ρ_*, for different values of the velocity dispersion σ (see e.g. Binney & Tremaine 2008) indicated in the figure with the colour map. In the corrective QuMOND term, $X=b_{\text{max}}^2/2r_M^2+b_{\text{max}}/r_M-1/2\ln (b_{\text{max}}/r_M)-(2+\ln 2)/4$ of Eq. (34) and in the standard Coulomb logarithm, we assumed b_max = (ρ/M_⊙)^−1/3 and b_min = G(M + m)/σ². It is evident that for an isolated baryon-only system (corresponding to a dark to luminous matter ratio, ℛ = 0) the DF in QuMOND is augmented (up to a factor ≈10²) at low density, comparable to the typical intra-galactic density of 10⁻¹ stars per pc³, with respect to the Newtonian case. Vice-versa (and as expected) in denser systems/regions the QuMOND correction is negligible, in particular for large values of the velocity dispersion, σ. Of course, a fair comparison between Newtonian and MONDian DF has to be made considering the presence of DM. In the bottom panel of Figure 1 we show η_QuMOND/η_N for the same values of ρ_* and σ and ℛ = 10 (a typical DM to baryonic matter ratio). In this case, the Newtonian DF is given by Eq. (25). Remarkably, in this case the QuMOND DF is substantially lower than the Newtonian one including the dark mass contribution (substantially confirming the results of Sánchez-Salcedo et al. 2006 for a simplified circular orbit model and of Nipoti et al. 2008; Roshan et al. 2021 for the dynamics of bars in spiral galaxies), for almost all astrophysically relevant combinations of stellar density and velocity dispersion. In the low ρ_* regime (say below ≈8 M_⊙pc⁻³), the simplified QuMOND DF expression dominates over its Newtonian counterpart. This leads us to speculate that DM-dominated low ρ_* systems, such as ultra-faint dwarf galaxies, could be used as a testing ground for MOND theories versus ΛCDM using the different behaviour of the DF (see e.g. Bílek et al. 2021; Oehm & Kroupa 2024). We stress the fact that in finite systems, both in Newtonian gravity and MOND, the orbital structure (i.e. the velocity anisotropy profile of the background particle distribution) has a relevant effect on the magnitude of DF. For example, a test mass, M, moving on a highly radial orbit in a radially anisotropic spherical system suffers less DF than in a systems with isotropic velocity distribution. In this case, in MOND, where only stars contribute to DF, the drag is further reduced when M reaches radii where no mass is present. A similar argument has been applied (among others) by Angus & Diaferio (2009) to explain in MOND the survival of the globular cluster systems of the Fornax dwarf spheroidal.

Fig. 1. Ratio of the QuMOND and Newtonian DF coefficients as function of the stellar density, for ℛ = 0 (top) and 10 (bottom), for the case of a M = 103 M⊙ test mass. The dependence on the local velocity dispersion σ is colour coded.

3. Mean field approach

A more rigorous evaluation of the DF expression can be carried out using the mean field formalism developed by Kandrup (1983) (see also Gilbert 1968 and Kandrup 1980) in the context of Newtonian gravity that can be applied to any kind of long-range force, not necessarily obeying the superposition principle.

We consider the usual system of N equal mass m (field) particles, with coordinates r_i and momenta, p_i = mv_i, described by the (time-dependent) phase-space distribution function³, ℱ(r, p; t), a test particle M, with a coordinate, r₀ = R, and momentum, p₀ = P = Mv_M, which perturbs the initial phase-space distribution ℱ₀. By virtue of the third law of dynamics, such perturbation corresponds to the force decelerating the test mass, M. We assume the field particles to be statistically uncorrelated in their initial state with Maxwellian velocity distribution. Moreover, we also take the limit of infinite, homogeneous distribution of field particles. In his original work, Kandrup (1983) initially assumed an external potential, Φ, confining the system (see Appendix C below). We note that MOND systems with self-consistent gravitational field, g_in, embedded in an external gravitational field g_ext are prone to the so-called external field effect (hereafter EFE). The latter implies that for g_in < a₀ < g_ext the system is purely Newtonian, while for g_in < g_ext < a₀ the system behaves as a Newtonian model with rescaled gravitational constant, G′=Ga₀/g_ext.

Under the assumptions given above, as usual, the average number density of particles, n = N/Vol, should be taken constant as N, Vol → ∞. The dynamical friction force $F_0^{\text{fr}}$ can be therefore obtained simply as ${\langle F_0^{\text{tot}}\rangle }_{\mathcal{F}}$.

In the same fashion of the special relativistic extensions of DF (see e.g. Syer 1994 and Chiari & Di Cintio 2023), we shift to the frame of reference of M such that for the new coordinates, $\widehat{\mathit{r}},\widehat{\mathit{p}}$, the linear evolution operator for ℱ takes the form $\dot{\mathcal{F}}({\widehat{\mathit{r}}}_i,{\widehat{\mathit{p}}}_i;t)=-i\mathcal{L}[\mathcal{F}]$. See Eqs. (C.7) and (C.9) in the appendix below for its explicit expression.

The linear formulation of the evolution equation is allowed since MOND (and therefore QuMOND), although non-linear, is still a Lagrangian and local⁴ theory. At this stage, it is not necessary to consider explicitly the Bekenstein & Milgrom (1984) (or Milgrom 2010) Langrangian, since it is enough to consider the particle’s energies.

The significant difference with respect to the Newtonian case explored in Kandrup (1983) is that in MOND, we cannot substitute the potentials in a form:

$$\begin{array}{c}\hfill W_i={\displaystyle \sum _{j\ne i}}{W}_{\mathit{ij}}(|{\mathit{r}}_i-{\mathit{r}}_j|),\end{array}$$(36)

due to the non-linearity of MOND implied by the absence of a Superposition principle. Such potentials however can be rather calculated in the QuMOND formalism (8).

A general expression for the DF force on the test particle can be written formally as

$$\begin{array}{c}\hfill {\mathit{F}}_0^{\mathrm{fr}}=-\beta {\displaystyle \int \mathrm{d}\mathrm{\Sigma }}{\mathcal{F}}_0\left({\displaystyle \sum _{j=1}^N}{\mathit{F}}_i^{\mathrm{tot}}\right){\displaystyle {\int }_0^t}\mathrm{d}\tau G_{\mathcal{L}}(\tau \to t)[{\displaystyle \sum _{i=1}^N}{\mathit{v}}_i·{\mathit{F}}_i^{\mathrm{tot}}],\end{array}$$(37)

where G_ℒ is the Greenian of the operator ℒ and dΣ = d^3Nrd^3Np is the differential element in phase space. The details of the derivation are discussed in Appendix C.

We now sketch the main simplifying hypotheses needed to perform practically the integrals in Eq. (37). Assuming that the system under consideration has a finite memory (as implied by its ergodicity), we can replace the time integration in τ on the finite interval [0; t] with an integration extended over the semi-infinite interval [0; ∞). Moreover, making use of the standard linear trajectory approximation (see e.g. Ter Haar 1977, see also Syer 1994), allows us to simplify G_ℒ(τ → t)[Q]≅Q(t − τ). The latter is the most delicate step of the present calculation, since in this approximation, the effects of the different field particles decouple that is, in principle, not valid in non-linear theory.

Under the assumptions listed above, the DF formula becomes:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}=\frac{{\mathit{F}}_0^{\mathrm{fr}}}{M}\cong & -\frac{3G^2{M\rho |}_{v<V_M}}{⟨v^2⟩}{\displaystyle {\int }_0^{\infty }}{\displaystyle \int {\int }_{\mathrm{\Omega }}}f(\mathit{v})\frac{\mathit{v}·(\mathit{s}-\stackrel{\sim }{\mathit{v}}\tau )}{|\mathit{s}-\stackrel{\sim }{\mathit{v}}{\tau |}^3}\hfill \\ \hfill & \times \nu \left(\frac{r_{\mathrm{M}}^2}{s^2}\right)\nu \left(\frac{r_{\mathrm{M}}^2}{|\mathit{s}-\stackrel{\sim }{\mathit{v}}{\tau |}^2}\right)\frac{\mathit{s}}{s^3}{d}^3\mathit{v}{d}^3\mathit{s}\mathrm{d}\tau ,\hfill \end{array}$$(38)

where $\stackrel{\sim }{\mathit{v}}=\mathit{v}-{\mathit{v}}_M$, Ω is the phase space volume occupied by the system and only the field particles slower than V_M induce the drag, so that we only for account their density, ρ|_v < VM. At this stage, the problem is reduced to solving a three dimensional integral for a given choice of the QuMOND interpolation function, ν.

Eq. (38) in the Newtonian limit (i.e. ν → 1) can be solved by an integration, incrementally applying:

$$\begin{array}{c}\hfill \mathrm{\nabla }·\frac{\mathit{s}}{s^3}=-4\pi {\delta }^3(\mathit{s})\end{array}$$(39)

and

$$\begin{array}{c}\hfill \mathit{v}·\frac{\mathit{s}-\stackrel{\sim }{\mathit{v}}\tau }{|\mathit{s}-\stackrel{\sim }{\mathit{v}}{\tau |}^3}=\mathrm{\nabla }·\frac{\mathit{v}}{|\mathit{s}-\stackrel{\sim }{\mathit{v}}\tau |},\end{array}$$(40)

thus, we have:

$$\begin{array}{c}\hfill {\displaystyle \int (\mathit{v}·\frac{\mathit{s}-\stackrel{\sim }{\mathit{v}}\tau }{|\mathit{s}-\stackrel{\sim }{\mathit{v}}{\tau |}^3})}\frac{\mathit{s}}{s^3}{d}^3\mathit{s}=4\pi \frac{\mathit{v}}{\stackrel{\sim }{v}\tau }.\end{array}$$(41)

Here, the term $\stackrel{\sim }{v}\tau =b$ is the standard impact parameter (see again Binney & Tremaine 2008), thus, the integral

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}\tau }{\tau }={\displaystyle {\int }_0^{\infty }}\frac{\mathit{db}}{b}\end{array}$$(42)

returns the usual Coulomb logarithm (see e.g. Kalnajs 1972), while the integral in d³v is rewritten as in the Newtonian case, therefore yielding Eq. (17).

The QuMOND case is, even for the simple form of ν given in Eq. (31), plagued by its intrinsic complexity. Analogously to the naive formula in Eq. (30) we notice the presence of two factors of ν. In this mean field formalism, the arguments for the two functions ν in (38) are still slightly different. From such difference $-\tau \stackrel{\sim }{\mathit{v}}$, we recovers the logarithmic term in τ, as well as any additional QuMOND term. At variance with the Newtonian case, which requires some cut-off value only for the integral over dτ, in the QuMOND DF formula of Eq. (38), the integral over d³s is also diverging; therefore, it requires a cut-off at b_max. This yields terms of the form b_max/r_M and $b_{\text{max}}^2/r_{\text{M}}^2$, which are vastly in agreement with the simplified approach described in Sect. 2.2.

In the deep MOND regime, where ν(y)→y^−1/2, we can simplify Eq. (38) as:

$$\begin{array}{cc}& \frac{\mathrm{d}{\mathit{v}}_M}{\mathrm{d}t}=-\frac{3G^{2}M{\rho }_{\ast }{|}_{v<V_M}}{⟨v^2⟩r_{\mathrm{M}}^2}{\displaystyle {\int }_0^{\infty }}{\displaystyle \int {\int }_{\mathrm{\Omega }}}f(\mathit{v})\frac{\mathit{v}·(\mathit{s}-\stackrel{\sim }{\mathit{v}}\tau )}{|\mathit{s}-\stackrel{\sim }{\mathit{v}}{\tau |}^2}\frac{\mathit{s}}{s^2}{d}^3\mathit{v}{d}^3\mathit{s}\mathrm{d}\tau =\hfill \\ \hfill & =-\frac{12{\pi }^2{G}^{2}M{\rho }_{\ast }{b}_{\mathit{max}}^2}{⟨v^2⟩r_M^2}\mathrm{\Psi }(V_M){\displaystyle {\int }_{v<V_M}}\frac{f(v)}{\stackrel{\sim }{v}}[\frac{{\pi }^2}{8}\mathit{v}-(1+\frac{3}{16}{\pi }^2)\frac{\mathit{v}·\stackrel{\sim }{\mathit{v}}}{{\stackrel{\sim }{v}}^2}\stackrel{\sim }{\mathit{v}}]d^3\mathbf{v}.\hfill \\ \hfill \end{array}$$(43)

We recall that Sánchez-Salcedo et al. (2006) (see also the numerical simulations of Bílek et al. 2021), used the circular orbit approximation and the estimates from Ciotti & Binney (2004) to formulate a local dMOND DF expression for a spherical system where the radial friction on M is expressed as:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{v}_c}{\mathrm{d}t}=-\frac{16{\pi }^{2}log\mathrm{\Lambda }{G}^{2}M{\rho }_{\ast }(r)}{v_c^2}\left(\frac{a_0^2}{\sqrt{2}{a}_{\mathrm{typ}}^2}\right)\mathrm{\Psi }(v_c).\end{array}$$(44)

In the expression above, $a_{\text{typ}}\approx G{M}_{*}/r_s^2$ is the typical Newtonian acceleration of the system with scale radius r_s, and the circular velocity profile is given (see Milgrom 1983) as a function of the stellar radial mass profile M_*(r), as follows:

$$\begin{array}{c}\hfill v_c(r)={\left[G{a}_0{M}_{\ast }(r)\right]}^{1/4}.\end{array}$$(45)

Taking the same assumptions applied in Sánchez-Salcedo et al. 2006, we can simplify Eq. (43) as v and s are generally likely to be perpendicular, while v and $\stackrel{\sim }{\mathit{v}}$ are almost parallel. This allows us to then simplify the expression of the term $\mathit{v}·(\mathit{s}-\stackrel{\sim }{\mathit{v}}\tau )/{|\mathit{s}-\stackrel{\sim }{\mathit{v}}\tau |}^2$ in scalar form and perform the integral in dτ before those in velocity, dv, and in space, ds, having imposed spherical symmetry. Standard (albeit tedious) algebra then yields

$$\begin{array}{c}\hfill \frac{\mathrm{d}{v}_c}{\mathrm{d}t}=-\frac{12{\pi }^2{G}^{2}M{\rho }_{\ast }(r)b_{\max }^{2}log2}{⟨v^2⟩r_M^2}\mathrm{\Psi }(v_c),\end{array}$$(46)

where the dMOND mean square velocity is given (Milgrom 1995) as

$$\begin{array}{c}\hfill {\sigma }^2=⟨v^2⟩=\sqrt{4G{a}_0{M}_{\ast }/81},\end{array}$$(47)

for a system with total baryonic mass M_* and b_max (typically of the order of the local mean inter-particle distance) is fixed to a scale radius r_s.

In Fig. 2, we show the DF acceleration in the limit of circular orbits on a mass of M = 2 × 10⁵ M_⊙, as function of the radial distance, r, as given in Eqs. (44) and (46), sinking in a cored spherical model (cfr. also Eq. 49 in the next section) with M_* = 2 × 10⁷ M_⊙ and a core radius of r_s = 0.3 kpc, corresponding to the model of a Fornax-like dwarf galaxy considered by Sánchez-Salcedo et al. 2006. With such a choice of parameters, the constant dMOND velocity dispersion given by Eq. (47) is ≈11km/s. We note that at large radii (e.g. above the core radius indicated in Figure by the vertical dashed line), the two expressions are somewhat comparable; however, in the inner part of the system, the DF is considerably lower in our approximated formulation. This is because for r → 0, the l.h.s. of Eq. (46) is mostly dominated by the contribution of Ψ(v_c), that for the specific (but questionable) choice of a Maxwellian velocity distribution is expressed as:

$$\begin{array}{c}\hfill \mathrm{\Psi }(v_c)={(4\pi )}^{-1}[\mathrm{erf}\left(\frac{v_c}{\sqrt{2}\sigma }\right)-\sqrt{\frac{2}{\pi }}\frac{v_c}{\sigma }exp(-\frac{v_c^2}{2{\sigma }^2})],\end{array}$$(48)

Fig. 2. Radial profiles of the DF force (per unit mass) in the circular orbit approximation as given by Eqs. (44,46); purple and cyan lines, respectively.

then falls rapidly to zero, while in Eq. (44), it is partially compensated by the factor $1/v_c^2$.

Interestingly, the QuMOND N-body simulations of Bílek et al. (2021) offer a hint that Eq. (44) overestimates the DF in the inner region of cored systems; whereas at large radii, the semianalytical estimate of Sánchez-Salcedo et al. (2006) is well matched by the simulations. The discrepancy is interpreted as an effect of the so-called dynamical buoyancy or core stalling (Banik & van den Bosch 2021) in a nearly harmonic core, preventing the further sinking of the test particle, M. We note that in MOND (and in particular in the dMOND limit assumed here), an almost constant density core does not exert a harmonic potential, thereby potentially undermining the core-stalling argument in MOND simulations. This leads to conjecture that the intrinsically less effective dMOND DF in the N-body simulation is indeed responsible for the discrepancy with the prediction of Eq. (44).

4. Numerical experiments

To clarify the different behaviour of DF in the QuMOND and Newtonian gravity model of a given stellar system, it is useful to perform a simple numerical experiment in a set-up where Newtonian gravity with DM the friction on the test object, M, would be negligible, while the gravitational field of the parent system would be in the MONDian regime, g ≪ a₀, over a broad interval of radii.

We consider a M = 10² M_⊙ star cluster moving through a dwarf galaxy with stellar mass of M_* = 2 × 10⁵ M_⊙ and scale radius of r_s = 0.8 kpc (parameters compatible to those of the Draco ultra-faint dwarf galaxy, Mashchenko et al. 2006), with a stellar distribution given by

$$\begin{array}{c}\hfill {\rho }_{\ast }(r)=\frac{3}{4\pi }\frac{M_{\ast }{r}_s}{{(r+r_s)}^4},\end{array}$$(49)

corresponding to a Dehnen (1993)γ-model for γ = 0. The associated QuMOND gravitational potential given by Eq. (8) is (by construction) identical to the potential exerted by a Newtonian system with stellar density, as in Eq. (49), plus a DM halo with ${\rho }_{\mathrm{DM}}=\stackrel{\sim }{\rho }-{\rho }_{\ast }$ (cfr. Eq. 12). By contrast, the DF force experienced by M is different in the two paradigms: Newtonian and MOND, with Eq. (25) displaying the first case⁵ and Eq. (34) the second.

Using the numerical approach discussed in Appendix D, (see also Pasquato & Di Cintio (2020), Di Cintio et al. (2020, 2022)), we integrated different orbits under the effect of the same gravitational field and the two different DF expressions, under the assumption of local Maxwellian approximation (i.e. at all radii the velocity distribution is an isotropic Maxwell-Boltzmann with velocity dispersion approximated by $\sigma =\sqrt{2||\mathrm{\Phi }||}$). In Figure 3, we show the evolution of the galactocentric distance, r_g, of M over 12 Gyr for an orbit of initial ellipticity, e = 0.49. The green and blue curves mark the Newtonian and QuMOND simulations, respectively, while the black curve is the unperturbed orbit in the static potential of the model. We observe that while the Newtonian DF only minimally alters the orbit of M at around 9 Gyr, in QuMOND, the test particle suffers a rather strong orbital decay, with the pergalactic radius falling from ≈1.63 kpc down to ≈1.2 kpc at 12 Gyr. This is remarkably in agreement with the analytical estimates of Ciotti & Binney (2004), who predicted a MOND inspiral time of less than a Hubble time for an object that is considerably more massive than a star to the central regions of a Draco-like system. This time span is much shorter than its Newtonian analog and argued that this could be the reason behind the absence of star clusters in similar dwarf galaxies.

Fig. 3. Evolution of the galactocentric radius, rg, for a star cluster orbiting through a dwarf galaxy in Newtonian gravity (green curve) and QuMOND (blue curve). The thin black line marks the case without dynamical friction.

5. Conclusions

We investigated the dynamical friction in the Quasi linear formulation of MOND. Using a simple dimensional analysis, we find that the expression for the MONDian DF is augmented with respect to its Newtonian counterpart of a factor proportional to the MOND radius of the star, r_M. This additional term becomes relevant when the maximum impact parameter (proportional to the average interparticle distance, Spitzer 1987) is b_max ≫ r_M, as seen by Eqs. (33 and 34). In practice, in a dense star cluster or in galactic nuclei with a mean stellar density and where the typical inter-stellar distance is always smaller than r_M, the DF is always purely Newtonian. Vice versa, in low-density systems, such as ultra-faint dwarf galaxies (Simon 2019) the b_max/r_M correction is always of the order of the Coulomb logarithm, thus enhancing the (stellar) DF for a system (where, in principle, it should be negligible). In addition, we find that the explicit form of the QuMOND correction is dependent on the specific choice of the interpolating function ν. However, the strength of the enhancing term does not vary significantly for such different forms of ν. We also compared our results with a Newtonian plus DM expression obtained from the Ciotti (2010) formalism for DF in presence of a mass spectrum, confirming that in almost all astrophysically relevant regimes, the Newtonian DF (in presence of dark matter) dominates over its MOND counterpart.

Moreover, we have extended the mean field treatment of DF, pioneered by Kandrup (1983) in the Newtonian case to QuMOND. In this framework, alternatively to the fluctuation-based approach by Ciotti & Binney (2004), we recovered an integral expression that contains the usual dependence on ν. Unfortunately, the explicit evaluation of the DF force becomes rapidly cumbersome (at variance with the simpler Newtonian case) even in the dMOND limit; thus, we are forced to integrate Eq. (38) numerically (specific cases of this will be discussed in a future work). Working out in QuMOND the two-body relaxation time (and therefore the DF coefficient of Eq. 6 via the relation t_DF = 2t_2bm/M) with the fluctuation approach in Fourier space of Ciotti & Binney (2004) may, in principle, be possible; however, it would imply a double application of the classical Poisson equation and a single non-linear algebraic step is hindered by the implicit relation between the Fourier transform of the QuMOND potential and the associated auxiliary density $\stackrel{\sim }{\rho }$. Using the mean field approach, we worked out an approximate formula, valid in the deep MOND regime for circular orbits analogous to the one presented in Sánchez-Salcedo et al. (2006). For the case of a globular cluster orbiting in a isolated Fornax-like dwarf galaxy, the two expressions are comparable, at least for radii larger than the core radius of the dwarf, below which our expression underestimates the DF. This result appears to be in agreement with the N-body simulation of Bílek et al. (2021).

Interestingly, simple numerical integration of a test particle in a QuMOND and its Newtonian equivalent system evidence that dynamical friction acts in a different way in the two paradigms, being considerably stronger in the first if the model is in a deep MOND regime (such as for DM-dominated dwarf galaxies). However, we stress the fact that in MOND, we cannot simply add a semi-analytic DF force to a mean field potential obtained solving numerically Eq. (1) for a given density (either imposed extrapolated from particles position) in the same fashion as in Alessandrini et al. (2014), Arca-Sedda & Capuzzo-Dolcetta (2014), Di Cintio et al. (2022). Therefore, our simple numerical estimates could potentially overestimate the MONDian DF. Using the quasi-linear formulation of MOND could, in principle, allow us to incorporate the contribution of DF (and possibly density fluctuations on a scale smaller than the particle-resolution) by adding it before evaluating $\stackrel{\sim }{\rho }$ to solve Eq. (8) (see Di Cintio 2023). In conclusion, it is also important to recall that as a consequence of the EFE, in systems with g_int < g_ext < a₀, the dynamics is essentially Newtonian with a rescaled G. In those cases, the DF force would be given by Eq. (17) augmented by the multiplicative factor (a₀/g_ext)². Therefore, a full mean field QuMOND treatment including the DF expression (discussed in Sect. 3) is only valid for an isolated system.

Acknowledgments

We express gratitude to Jan Pflamm-Altenburg, Pavel Kroupa and Luca Ciotti for the useful discussions at an early stage of this work. We also thank the anonymous Referee for his/her useful remarks that helped improving the presentation of our results.

References

Appendix A: Effective dark matter halos

Combining Eqs. (9) and (8), we evaluate the equivalent matter distribution from Eq. (12) as

$$\begin{array}{c}\hfill \stackrel{\sim }{\rho }=-\frac{1}{4\pi G}\mathbf{\nabla }·\mathit{g}=-\frac{1}{4\pi G}\mathbf{\nabla }·\left[\nu (|{\mathit{g}}_N|/a_0){\mathit{g}}_N\right].\end{array}$$(A.1)

For a point particle, the density becomes ρ(r) = mδ³(r), so that g_N(r) = − Gm/r³r, and

$$\begin{array}{c}\hfill \stackrel{\sim }{\rho }=-\frac{1}{4\pi G}\mathbf{\nabla }·[-\frac{\mathit{Gm}}{r^3}\mathit{r}\nu \left(\frac{\mathit{Gm}}{a_0{r}^2}\right)]=\frac{m}{4\pi }\mathbf{\nabla }·\left[\frac{\mathit{r}}{r^3}\nu \left(\frac{r_{\mathrm{M}}^2}{r^2}\right)\right].\end{array}$$(A.2)

Substituting the simple interpolation (31), we then find

$$\begin{array}{c}\hfill \stackrel{\sim }{\rho }=\frac{m}{4\pi }\mathbf{\nabla }·\left(\frac{\mathit{r}}{r^3}\right)+\frac{m}{4\pi r_{\mathrm{M}}}\mathbf{\nabla }·\left(\frac{\mathit{r}}{r^2}\right)=m{\delta }^3(\mathit{r})+\frac{m}{4\pi r_{\mathrm{M}}{r}^2},\end{array}$$(A.3)

where the first term here is the real point particle, and the second one represents its phantom dark matter halo ${\rho }_{\mathit{pDM}}=\stackrel{\sim }{\rho }-\rho$. Its cumulative mass function is

$$\begin{array}{c}\hfill m_{\mathit{DM}}(r)={\displaystyle \int \frac{m}{4\pi r_{\mathrm{M}}{r}^2}}4\pi r^{2}dr=m\frac{r}{r_{\mathrm{M}}}.\end{array}$$(A.4)

Following the approach of Chandrasekhar, the impulsive force generated by the particle and its effective DM halo is then

$$\begin{array}{c}\hfill F=\frac{GMm(b)}{b^2}=\frac{\mathit{GMm}}{b^2}+\frac{\mathit{GMm}}{r_{\mathrm{M}}b}.\end{array}$$(A.5)

The associated momentum variation is $\mathrm{\Delta }p=Fb/v_M=\mathit{GMm}/v_M(b^{-1}+r_{\mathrm{M}}^{-1})$ and then the diffusion coefficient becomes

$$\begin{array}{cc}\hfill D_v& ={\displaystyle \int \mathrm{\Delta }}{p}^{2}2\pi n{v}_{M}ddb=2\pi G^2{M}^2\frac{m^{2}n}{v_M}{\displaystyle \int {(\frac{1}{b}+\frac{1}{r_{\mathrm{M}}})}^2}bdb.\hfill \end{array}$$(A.6)

Keeping in mind that b_m ≪ r_M < b_max, the integral can be evaluated as $\approx \ln \Lambda +2b_{\text{max}}/r_{\text{M}}+b_{\text{max}}^2/2r_{\text{M}}^2$. The QuMOND DF force is then

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{v}}_M}{\mathit{dt}}=-2\rho (M+m)D_v\frac{N(|\mathit{v}|<v_M)}{v_M^2}{\mathit{v}}_M=\\ \hfill =-16{\pi }^2(ln\mathrm{\Lambda }+2\frac{b_{\max }}{r_{\mathrm{M}}}+\frac{b_{\max }^2}{2r_{\mathrm{M}}^2})G^2\rho (M+m)\frac{{\mathit{v}}_M}{v_M^3}{\displaystyle {\int }_0^{v_M}}{v}^{2}f(v)dv.\end{array}$$(A.7)

Appendix B: Explicit choice of the QuMOND interpolating function

If the usual choice (13) of the QuMOND interpolating function ν is used in Eq. (30), we have

$$\begin{array}{c}\hfill \nu (r_{\mathrm{M}}^2/b^2)=\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{b^2}{r_{\mathrm{M}}}}\Rightarrow \nu {(b)}^2=\frac{1}{2}+\frac{1}{2}\sqrt{1+4\frac{b^2}{r_{\mathrm{M}}^2}}+\frac{b^2}{r_{\mathrm{M}}^2}.\end{array}$$(B.1)

Unfortunately, here we cannot take the usual approximation for small or large impact parameters b, since the integral runs from b_min < r_M to b_max > r_M. Vary large values of b_max (say of the order of the system size), return an important contribution already in the Newtonian case, in the form of large values of logΛ. This must also hold true a fortiori in MOND, as in the far-field the forces decay as 1/r rather than 1/r². Therefore, we cannot expand ν to low-order terms in b/r_M.

We must then solve the exact integral ∫ν(b)²db/b, where the square root term yields

$$\begin{array}{cc}\hfill {\displaystyle \int \sqrt{1+4b^2/r_{\mathrm{M}}^2}}\frac{\mathit{db}}{b}& ={\displaystyle \int \frac{\sqrt{1+x^2}}{x}}dx=\hfill \\ \hfill & =\sqrt{1+x^2}-\text{arctanh}(\sqrt{1+x^2})+\mathrm{const}.=\hfill \\ \hfill & =\sqrt{1+x^2}-\frac{1}{2}ln\frac{\sqrt{x^2+1}+1}{\sqrt{x^2+1}-1}+(\mathrm{const}-\frac{\pi }{2}i),\hfill \\ \hfill \end{array}$$(B.2)

where x = 2b/r_M. Evaluating Eq. (B.2) for x = 2b_min/r_M → 0 returns approximately

$$\begin{array}{c}\hfill (1+x^2/2)-\frac{1}{2}ln\frac{1+1}{(1+x^2/2)-1}=ln\frac{b_{\min }}{r_{\mathrm{M}}}+1,\end{array}$$(B.3)

while for x = 2b_max/r_M → ∞ becomes asymptotically

$$\begin{array}{c}\hfill (x+\frac{1}{2x})-\frac{1}{2}ln\frac{x+1}{x-1}=2\frac{b_{\max }}{r_{\mathrm{M}}}.\end{array}$$(B.4)

The integral (B.2) is

$$\begin{array}{c}\hfill {\displaystyle \int \sqrt{1+4b^2/r_{\mathrm{M}}^2}}\frac{\mathit{db}}{b}\sim 2\frac{b_{\max }}{r_{\mathrm{M}}}-ln\frac{b_{\min }}{r_{\mathrm{M}}}-1.\end{array}$$(B.5)

The integral over the impact parameters is then evaluated as

$$\begin{array}{cc}\hfill I_b& =\frac{1}{2}{\displaystyle \int \frac{\mathit{db}}{b}}+\frac{1}{2}{\displaystyle \int \sqrt{1+4b^2/r_{\mathrm{M}}^2}}\frac{\mathit{db}}{b}+\frac{1}{r_{\mathrm{M}}^2}{\displaystyle \int b}db=\hfill \\ \hfill & \sim \frac{1}{2}ln\frac{b_{\max }}{b_{\min }}+\frac{b_{\max }}{r_{\mathrm{M}}}-\frac{1}{2}ln\frac{b_{\min }}{r_{\mathrm{M}}}-\frac{1}{2}+\frac{b_{\max }^2}{2r_{\mathrm{M}}^2}=\hfill \\ \hfill & =ln\mathrm{\Lambda }+\frac{b_{\max }^2}{2r_{\mathrm{M}}^2}+\frac{b_{\max }}{r_{\mathrm{M}}}-\frac{1}{2}ln\frac{b_{\max }}{r_{\mathrm{M}}}-\frac{1}{2}.\hfill \end{array}$$(B.6)

Appendix C: Further details on the mean field approach

For the system of N + 1 particles considered in Section 3 described by the distribution function ℱ(r, p; t) the average of a given phase-space observable Q, taken with respect to said distribution function is

$$\begin{array}{c}\hfill {⟨Q⟩}_{\mathcal{F}}(t):={\displaystyle \int \mathcal{F}}(\mathbf{r},\mathbf{p};t)Q(\mathbf{r},\mathbf{p})\mathrm{d}\mathrm{\Sigma },\end{array}$$(C.1)

where dΣ = d^3Nrd^3Np is the differential phase-space element.

The initial configuration ℱ₀ is perturbed by the passage of the (N + 1)-th particle, (our usual test particle). In virtue of the Third Law of Dynamics, this corresponds to the force decelerating the test particle once the force due to the unperturbed configuration ℱ₀ is subtracted:

$$\begin{array}{c}\hfill F_0^{\mathit{fr}}={⟨F_0^{\mathit{tot}}⟩}_{\mathcal{F}}-{⟨F_0^{\mathit{tot}}⟩}_{{\mathcal{F}}_0},F_0^{\mathit{tot}}=-{\displaystyle \underset{N}{\overset{}{{\displaystyle \sum i=1}}}}{F}_i^{\mathit{tot}}.\end{array}$$(C.2)

The field particles are assumed statistically uncorrelated in their initial state ${\mathcal{F}}_0(\mathbf{r},\mathbf{p})={\displaystyle \underset{N}{\overset{}{{\displaystyle \prod i=1}}}}f({\mathbf{r}}_i,{\mathbf{p}}_i)$. We also assume their velocity distribution to be a Maxwellian:

$$\begin{array}{c}\hfill f({\mathbf{r}}_i,{\mathbf{p}}_i)=\mathcal{N}{f}_1({\mathbf{r}}_i)e^{-\beta {\mathbf{p}}_i^2/2m_i},\end{array}$$(C.3)

with temperature β and normalization 𝒩. The distribution of the spatial coordinates is chosen such that the field particles are limited within a certain volume; therefore, it can be also expressed as f₁(r_i) = e^{−βm_iΦ_i(r_i)}, with an efficient potential Φ_i.

For our purposes, we set f₁ as the limit of an infinite, homogeneous distribution of field particles. Such a homogeneous distribution on a finite volume Vol ⊆ ℝ³ is given for $f_1({\mathbf{r}}_{\mathbf{i}})=\frac{1}{|Vol|}{\chi }_{\mathit{Vol}}({\mathbf{r}}_i)$, namely, when Φ_i has the profile of a rigid box. The average number density of particles can be then defined as n = N/|Vol|. In the limit of an infinite homogeneous system, this parameters should be taken constant as N, |Vol|→∞. The unperturbed force ${\langle F_0^{tot}\rangle }_{{\mathcal{F}}_0}$ vanishes, so that the dynamical friction $F_0^{fr}$ can be obtained simply as ${\langle F_0^{tot}\rangle }_{\mathcal{F}}$.

The field particles evolves following the time-dependent Hamiltonian:

$$\begin{array}{c}\hfill H(\mathbf{r},\mathbf{p};t)={\displaystyle \sum _{i=1}^N}\frac{{\mathbf{p}}_{\mathbf{i}}^{\mathbf{2}}}{2m_i}+{\displaystyle \sum _{i=1}^N}{W}_i(\mathbf{r}).\end{array}$$(C.4)

As discussed in Sect. (3), we cannot make the usual substitution of W_i = ∑_j ≠ iW_ij(|r_i − r_j|) with W_ij(r) = − Gm_im_j/r. This is because, in MOND, the superposition principle does not hold. The potentials should be rather expressed in using the QuMOND formalism (8), so that:

$$\begin{array}{c}\hfill -\mathrm{\nabla }{W}_i=m_i\nu \left(\frac{|{\mathbf{g}}_{\mathit{Ni}}|}{a_0}\right){\mathbf{g}}_{\mathit{Ni}},s.t.-{\displaystyle \sum _{j\ne i}}\mathrm{\nabla }{W}_{\mathit{ij}}=m_j{\mathbf{g}}_{\mathit{Ni}}.\end{array}$$(C.5)

As a consequence, it is useless to define the forces between couples of particles F(j → i) = − ∂_riW_ij(|r_i − r_j|) = −F(i → j), given that the total QuMOND force on the i−th particle is not simply ${\mathbf{F}}_i^{tot}\ne {\displaystyle \sum _{j\ne 0}}\mathbf{F}(j\to i)$. The force ${\mathbf{F}}_i^{tot}$ can be calculated instead, only as

$$\begin{array}{c}\hfill {\mathbf{F}}_i^{\mathit{tot}}=-{\partial }_{{\mathbf{r}}_{\mathbf{i}}}{W}_i(\mathbf{r}),\end{array}$$(C.6)

where W_i(r) is given by the non-linear QuMOND formula (C.5). The evolution equation is therefore expressed as

$$\begin{array}{c}\hfill \dot{\mathcal{F}}=-{\displaystyle \sum _{i=1}^N}\frac{{\mathbf{p}}_{\mathbf{i}}}{m_i}·\frac{\partial \mathcal{F}}{\partial {\mathbf{r}}_i}-{\displaystyle \sum _{i=1}^N}{\mathbf{F}}_i^{\mathit{tot}}·\frac{\partial \mathcal{F}}{\partial {\mathbf{p}}_i}:=-iL[\mathcal{F}].\end{array}$$(C.7)

The sums here arise from those in the Hamiltonian (C.4). In practice, the meaning of Eq. (C.7) is just that the position of each field particle evolves according to the velocity at each instant, as well as the momentum follows the total force.

Following Kandrup (1983), it is useful to change the phase-space coordinates by taking the test particle’s reference:

$$\begin{array}{c}\hfill {\widehat{\mathbf{r}}}_i:={\mathbf{r}}_i-\mathbf{R}(t)+\mathbf{R}(0),{\widehat{\mathbf{p}}}_i:={\mathbf{p}}_i-\mathbf{P}(t).\end{array}$$(C.8)

The evolution operator for ℱ hence takes the form

$$\begin{array}{c}\hfill \dot{\mathcal{F}}({\widehat{\mathbf{r}}}_i,{\widehat{\mathbf{p}}}_i;t)=-iL[\mathcal{F}]+{\displaystyle \sum _{i=1}^N}\dot{\mathbf{P}}(t)·\frac{\partial \mathcal{F}}{\partial {\mathbf{p}}_i}:=-i\mathcal{L}[\mathcal{F}],\end{array}$$(C.9)

where Ṗ acts as an apparent force on the i-th particle. Except for this, in this frame the forces do not change ${\widehat{F}}_i^{\mathit{tot}}\equiv F_i^{\mathit{tot}}$, since even in MOND they depend only on the respective particles positions. Moreover, it is useful to define the perturbation $\stackrel{\sim }{\mathcal{F}}(t):=\mathcal{F}(t)-{\mathcal{F}}_0$ of the distribution function so that

$$\begin{array}{c}\hfill {⟨Q⟩}_{\mathcal{F}}={⟨Q⟩}_{{\mathcal{F}}_0}+{⟨Q⟩}_{\stackrel{\sim }{\mathcal{F}}}\end{array}$$(C.10)

for any quantity Q. This will be used for evaluating Eq. (C.2), so that ${\mathbf{F}}_0^{\mathit{fr}}={⟨{\mathbf{F}}_0^{\mathit{tot}}⟩}_{\stackrel{\sim }{\mathcal{F}}}$. The distribution function $\stackrel{\sim }{\mathcal{F}}$ evolves as

$$\begin{array}{c}\hfill \dot{\stackrel{\sim }{\mathcal{F}}}=-i\mathcal{L}[\stackrel{\sim }{\mathcal{F}}]+\beta {\mathcal{F}}_0(\mathbf{r},\mathbf{p}){\displaystyle \sum _{i=1}^N}[\frac{{\widehat{\mathbf{p}}}_i}{m_i}+\mathbf{V}(t)]·[F_i^{\mathit{tot}}+m_i\frac{\partial {\widehat{\mathrm{\Phi }}}_i}{\partial {\widehat{\mathbf{r}}}_i}(t)],\end{array}$$(C.11)

where again ${\widehat{\mathbf{p}}}_i/m_i+\mathbf{V}(t)={\mathbf{p}}_i/m_i$, and the term $\partial {\widehat{\mathrm{\Phi }}}_i/\partial {\widehat{\mathbf{r}}}_i$ vanishes in the infinite, homogeneous limit. In principle, keeping the contribution of $\widehat{\mathrm{\Phi }}$ would allow us to investigate in a self-consistent fashion the MOND EFE, see Sect. 3. At this point, Eq. (C.11) can be rewritten in its simpler form:

$$\begin{array}{c}\hfill \dot{\stackrel{\sim }{\mathcal{F}}}=-i\mathcal{L}[\stackrel{\sim }{\mathcal{F}}]+\beta {\mathcal{F}}_0{\displaystyle \sum _{i=1}^N}{\mathbf{p}}_i·{\mathbf{g}}_i^{\mathit{MOND}}.\end{array}$$(C.12)

Its solution can thus be formally expressed as

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathcal{F}}(t)=\beta {\mathcal{F}}_0{\displaystyle {\int }_0^t}d\tau G_{\mathcal{L}}(\tau \to t)[{\displaystyle \sum _{i=1}^N}{\mathbf{p}}_i·{\mathbf{g}}_i]\end{array}$$(C.13)

with the Greenian G_ℒ of the operator ℒ. After all the substitution we have a formal, general expression for the dynamical friction on the test particle in the form

$$\begin{array}{cc}\hfill {\mathbf{F}}_0^{\mathit{fr}}& ={\displaystyle \int d}\mathrm{\Sigma }{\mathbf{F}}_0^{\mathit{tot}}\stackrel{\sim }{\mathcal{F}}={\displaystyle \int d}\mathrm{\Sigma }{\mathbf{F}}_0^{\mathit{tot}}\beta {\mathcal{F}}_0{\displaystyle {\int }_0^t}d\tau G_{\mathcal{L}}(\tau \to t)[{\displaystyle \sum _{i=1}^N}{\mathbf{v}}_i·{\mathbf{F}}_i^{\ast }]=\hfill \\ \hfill & =-\beta {\displaystyle \int d}\mathrm{\Sigma }{\mathcal{F}}_0\left({\displaystyle \sum _{j=1}^N}{\mathbf{F}}_i^{\mathit{tot}}\right){\displaystyle {\int }_0^t}d\tau G_{\mathcal{L}}(\tau \to t)[{\displaystyle \sum _{i=1}^N}{\mathbf{v}}_i·{\mathbf{F}}_i^{\ast }],\hfill \end{array}$$(C.14)

where ${\mathbf{F}}_i^{\ast }={\mathbf{F}}_i^{\mathit{tot}}+m_i\partial {\widehat{\mathrm{\Phi }}}_i/\partial {\widehat{\mathbf{r}}}_i$ becomes ${\mathbf{F}}_i^{tot}$ in the infinite homogeneous system limit. Here we have exploited the fact that the full set of N + 1 particles is a closed, classic dynamical system, so that ${\displaystyle \sum _{i=0}^N}{\mathbf{F}}_i^{\mathit{tot}}=0$.

To perform the integrals appearing in Eqs. (C.14), we now make some simplifying hypotheses. First of all, we perform the infinite homogeneous limit. Second, if the system has only a finite memory, we can replace the integration of τ on [0; t] with an integration of [0; ∞). Moreover, we make use of the standard linear trajectory approximation, which allows us to simplify G_ℒ(τ → t)[Q]≅Q(t − τ). Under these approximations, the DF force becomes

$$\begin{array}{c}\hfill {\mathbf{F}}_0^{\mathit{fr}}\cong -\beta {\displaystyle \int d}\mathrm{\Sigma }{\mathcal{F}}_0\left({\displaystyle \sum _{j=1}^N}{\mathbf{F}}_j^{\mathit{tot}}(t)\right){\displaystyle {\int }_0^{\infty }}d\tau ({\displaystyle \sum _{i=1}^N}{\mathbf{v}}_i(t-\tau )·{\mathbf{F}}_i^{\mathit{tot}}(t-\tau )).\end{array}$$(C.15)

This can be further simplified since ⟨F_iF_j⟩ = 0 for the off-diagonal combinations i ≠ j, since the field particles are by definition statistically uncorrelated. We are then left with

$$\begin{array}{cc}\hfill {\mathbf{F}}_0^{\mathit{fr}}& \cong -\beta {\displaystyle \int d}\mathrm{\Sigma }{\mathcal{F}}_0{\displaystyle \sum _{i=1}^N}{\mathbf{F}}_i^{\mathit{tot}}(t){\displaystyle {\int }_0^{\infty }}d\tau {\mathbf{v}}_i(t-\tau )·{\mathbf{F}}_i^{\mathit{tot}}(t-\tau )=\hfill \\ \hfill & =-\beta {\displaystyle {\int }_0^{\infty }}d\tau {\langle {\displaystyle \sum _{i=1}^N}{\mathbf{F}}_i^{\mathit{tot}}(t)({\mathbf{v}}_i(t-\tau )·{\mathbf{F}}_i^{\mathit{tot}}(t-\tau ))\rangle }_{{\mathcal{F}}_0}=\hfill \\ \hfill & =-\beta {\displaystyle {\int }_0^{\infty }}d\tau N{\langle {\mathbf{F}}_1^{\mathit{tot}}(t)({\mathbf{v}}_1(t-\tau )·{\mathbf{F}}_1^{\mathit{tot}}(t-\tau ))\rangle }_{{\mathcal{F}}_0}.\hfill \end{array}$$(C.16)