11institutetext: CNR-ISC, via Madonna del Piano 17 50022 Sesto Fiorentino,
Italy 22institutetext: INAF-Osservatorio Astronomico di Arcetri, Largo Enrico
Fermi 5 50125 Firenze Italy 33institutetext: Dipartimento di Fisica e
Astronomia, Universitá di Firenze & INFN-Sezione di Firenze, via Sansone 1
50022 Sesto Fiorentino, Italy
33email<EMAIL_ADDRESS>44institutetext: Dipartimento di
Fisica ”Giuseppe Occhialini”, Universitá di Milano Bicocca, Piazza della
Scienza 3 20126, Milano, Italy 55institutetext: INFN-Sezione di Milano Via
Celoria 15 20133, Milano, Italy
55email<EMAIL_ADDRESS>66institutetext: Dipartimento di Scienze
Fisiche, Informatiche e Matematiche, Universitá di Modena e Reggio Emilia, Via
Campi 213/A, I-41125 Modena, Italy 77institutetext: CNR-NANO, via Campi 213/A
I-41125, Modena, Italy
77email<EMAIL_ADDRESS>
# Dynamical friction in the quasi-linear formulation of MOND

Pierfrancesco Di Cintio 112233 Federico Re 4455 Caterina Chiari 6677

(Received ??; accepted ??)

###### Abstract

Aims. We explore the dynamical friction on a test mass in gravitational
systems in the Quasi linear formulation of Modified Newtonian Dynamics
(QuMOND).

Methods. Exploiting the quasi linearity of QuMOND we derive a simple
expression for the dynamical friction in akin to its Newtonian counterpart in
the standard Chandrasekhar derivation. Moreover, adopting a mean field
approach based on the Liouville equation we obtain a more rigorous (though in
integral form) dynamical friction formula that can be evaluated numerically
for a given choice of the QuMOND interpolation function.

Results. Consistently with previous work, we observe that dynamical friction
is stronger in MOND with respect to a baryon only Newtonian system with the
same mass distribution. This amounts to a correction of the Coulomb
logarithmic factor via extra terms proportional to the MOND radius of the
system. Morover, with the aid of simple numerical experiments we confirm our
theoretical predictions and those of previous work on MOND.

###### Key Words.:

Galaxies: kinematics and dynamics - Stars: kinematics and dynamics -
Gravitation - Methods: analytical

##  1 Introduction

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

$\@vec{\nabla}\cdot\left[\mu\left(\frac{||\@vec{\nabla}\Phi||}{a_{0}}\right)\@vec{\nabla}\Phi\right]=4\pi
G\rho,$ (1)

where the acceleration scale $a_{0}\approx 1.2\times 10^{-2}{\rm ms}^{-2}$ is
a new universal constant, and the in principle unknown interpolating function
$\mu$ is monotonic with the asymptotic behavior

$\mu(x)\rightarrow^{x\gg 1}1,\quad\mu(x)\sim^{x\ll 1}x.$ (2)

A widely adopted form of $\mu$ is

$\mu(x)=\frac{x}{1+x}.$ (3)

According to the equations above, for $\nabla\Phi\gg a_{0}$ one recovers the
Newtonian limit and Eq. (1) reduces to the Poisson equation. Vice versa, for
$\nabla\Phi\ll a_{0}$ the system is in the so-called deep-MOND limit
(hereafter dMOND) and Eq. (1) becomes the p-Laplace equation

$\@vec{\nabla}\cdot\left({||\@vec{\nabla}\Phi||}\@vec{\nabla}\Phi\right)=4\pi
a_{0}G\rho.$ (4)

For a given mass density $\rho$ the right-hand-side of the Poisson equation
and Eq. (1) are the same. One can therefore eliminate $\rho$ obtaining the
relation

$\mu\left(\frac{||\mathbf{g}_{M}||}{a_{0}}\right)\mathbf{g}_{M}=\mathbf{g}_{N}+\mathbf{S}$
(5)

between the MOND and Newtonian force fields $\mathbf{g}_{M}$ and
$\mathbf{g}_{N}$, where $\mathbf{S}\equiv\nabla\times\mathbf{h}(\rho)$ 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 of DM (Bugg 2015; Zhu et al. 2023), and has proven to be a
valid alternative to the $\Lambda$CDM paradigm in broad range of gravitational
phenomena (Milgrom 1994; Sánchez-Salcedo & Hernandez 2007; McGaugh & Milgrom
2013; Sanders 2021; Bílek et al. 2021; Asencio et al. 2022; 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; Kroupa et al. 2022; Banik et al.
2022; Nagesh et al. 2022, 2023; Re & Di Cintio 2023; Sollima et al. 2023), for
collisionless mean field processes, such for example violent relaxation and
phase-mixing in monolithic collapse, galaxy merging, vertical dynamics of disk
galaxies. Unfortunately, due to the non-linearity of the theory (cfr. Eq. 1
and hence the absence of the superposition principle, much less is known for
what concerns the collisional processes (see e.g. Bílek et al. 2021). For
example, Ciotti & Binney (2004) using an approach based on the fluctuations of
uniform field estimated 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) is 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, while would undergo an only slightly more
efficient DF than in the equivalent Newtonian system (hereafter ENS, i.e. 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^{\rm
N}_{\rm DF}$ and $t^{\rm M}_{\rm DF}$, are related as

$t^{\rm M}_{\rm DF}=\frac{\sqrt{2}}{1+\mathcal{R}}t^{\rm N}_{\rm DF},$ (6)

where $\mathcal{R}$ is the ratio of the amount of DM to baryons in the ENS.
Since in MOND collisional direct $N-$body simulations are intrinsically
unfeasible, results such as Eq. (6) can not be exhaustively explored via
numerical experiments; the reason being the absence of a MONDian expression
for the force exerted by a point-like particle. Milgrom (1986) proposed the
approximated expression for the force exchanged between two masses $m_{1}$ and
$m_{2}$ placed at distance $r$

$F_{1,2}\approx\frac{m_{1}m_{2}}{\sqrt{m_{1}+m_{2}}}\frac{\sqrt{Ga_{0}}}{r},$
(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^{2}$ limit.
More recently, Milgrom (2010) formulated a quasi-linear MOND theory (hereafter
QuMOND) where the governing field equation is

$\@vec{\nabla}\cdot\left[\nu\left(\frac{||\@vec{\nabla}\Phi^{N}||}{a_{0}}\right)\@vec{\nabla}\Phi^{N}\right]=\Delta\Phi^{M},$
(8)

formally identical to Equation (1) where now, $\Phi^{N}$ is the Newtonian
potential generated by the density $\rho$ through the usual Poisson equation

$\Delta\Phi^{N}=-\@vec{\nabla}\cdot\@vec{g}^{N}=4\pi G\rho.$ (9)

and the new interpolating function $\nu(y)$ is related to $\mu(x)$ by

$\cases{x}=y\nu(y)\\\
y=x\mu(x){}\Rightarrow\mu(x)=\frac{y}{x}=\frac{1}{\nu(y)}.$ (10)

The asymptotic behaviours for $\nu$ in the Newtonian and in the dMOND limits
are

$\nu(y)\rightarrow^{y\gg 1}1,\quad\nu(y)\sim^{y\ll 1}y^{-1/2}.$ (11)

We note that, Equation (8) can be rearranged as

$\tilde{\rho}=-\frac{1}{4\pi
G}\nabla\cdot[\nu(g^{N}/a_{0})\mathbf{g}^{N}];\quad\mathbf{g}^{N}=-\nabla\Phi^{N},$
(12)

that is, the QuMOND potential satisfies the Poisson equation for the auxiliary
density $\tilde{\rho}$. In practice, in QuMOND one has to evaluate the
Newtonian potential $\Phi^{N}$ generated by $\rho$, evaluate its gradient and
then via a non-linear algebraic step obtain the auxiliary MOND density
$\tilde{\rho}$ via Eq. (12) that acts as source for the potential $\Phi^{M}$.
At variance with the usual MOND formulation, in its quasi-linear formulation
the interpolating function $\nu$ has a stronger effect on the form of the
gravitational potential than its counterpart $\mu$ in Eq. (1). For the usual
choice (3), one obtains

$\displaystyle\nu(y)=\frac{1}{2}+\frac{1}{2}\sqrt{1+\frac{4}{y}}.$ (13)

We note that, QuMOND contains a sort of superposition principle for the
auxiliary densities $\tilde{\rho}$. It is therefore tantalizing 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 the contributions of the
effective densities for point-particles of mass $m$ interacting with a test
mass, with some simple choices of $\nu$. 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, let us briefly recall the ingredients
of the derivation of its Newtonian expression (for a more detailed review see
e.g. Binney & Tremaine 2008 or Ciotti 2021). Let us consider a test mass $M$
travelling at $v_{M}$ through an infinite background medium of particles with
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 $F\approx{GMm}/{b^{2}}$, acting for an interval of time $\Delta
t\approx{b}/{v_{M}}$, where $b$ is the impact parameter of the dynamical
collision. The associated variation of momentum is $\Delta p\approx F\Delta
t\approx{GMm}/{v_{M}b}$. The rate of collisions of $M$ is then

$\dot{N}_{\rm enc}=\frac{2\pi nbdbdx}{dt}=2\pi nv_{M}bdb.$ (14)

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

$D_{v}=\int\Delta p^{2}\frac{dN_{\rm enc}}{dt}=2\pi
G^{2}M^{2}\frac{m^{2}n}{v_{M}}\int\frac{db}{b}.$ (15)

In the Equation above, the integral in the impact parameter $b$ diverges for
both $b\rightarrow 0$ and $b\rightarrow\infty$ (infrared and ultraviolet
divergence, respectively). It is however possible to impose suitable cut-offs
$b_{\rm min}$ and $b_{\rm max}$. This factor, after integration takes the form
of the Coulomb logarithm (in analogy with the same quantity in plasma physics)
$\ln\Lambda$, where $\Lambda=b_{\rm max}/b_{\rm min}$.
Let us now assume that $f(v)$, defined such that

$\int f(\@vec{x},\@vec{v})d^{3}\@vec{x}d^{3}\@vec{v}=N$ (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

$\displaystyle\frac{d\mathbf{v}_{M}}{dt}$
$\displaystyle=-16\pi^{2}G^{2}\rho(M+m)\ln\Lambda\frac{\@vec{v}_{M}}{v_{M}^{3}}\int_{0}^{v_{M}}v^{2}f(v)dv.$
(17)

In the assumption of velocity isotropy, only particles moving slower than $M$
contribute to its slowing down. Moreover for $M\gg m$, in Eq. (17) the factor
$m(M+m)n$ becomes $M\rho$, where $\rho$ is the background mass density.

###  2.2 QuMOND case

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

$F=\nu\left(\frac{Gm}{a_{0}b^{2}}\right)\frac{GMm}{b^{2}}.$ (18)

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

$r_{\rm M}=\sqrt{\frac{Gm}{a_{0}}}$ (19)

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

$\Delta p=\nu(r_{\rm M}^{2}/b^{2})GMm/v_{M}b.$ (20)

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

$I_{b}=\int_{b_{\rm min}}^{b_{\rm max}}\nu(r_{\rm M}^{2}/b^{2})^{2}db/b.$ (21)

Once substituting in the DF formula, one obtains

$\frac{d\mathbf{v}_{M}}{dt}=-16\pi^{2}I_{b}G^{2}\rho(M+m)\frac{\@vec{v}_{M}}{v_{M}^{3}}\int_{0}^{v_{M}}v^{2}f(v)dv.$
(22)

From Equation (22) above we notice immediately that the QuMOND dynamical DF is
always larger than its Newtonian expression for the same mass distribution.
This happens because $\nu^{2}>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 Equation (22) for a simple choice of $\nu$. Knowing that $\nu$
has the asymptotic trends given in Eq. (11) we start with a simplified form

$\nu(y)=1+1/\sqrt{y},$ (23)

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

$\displaystyle I_{b}$ $\displaystyle=\ln\frac{b_{\rm max}}{b_{\rm
min}}+2\frac{b_{\rm max}-b_{\rm min}}{r_{\rm M}}+\frac{b_{\rm max}^{2}-b_{\rm
min}^{2}}{2r_{\rm M}^{2}}=$ (24) $\displaystyle\cong\ln\Lambda+2\frac{b_{\rm
max}}{r_{\rm M}}+\frac{b_{\rm max}^{2}}{2r_{\rm M}^{2}};$ (25)

where we assumed $b_{\rm min}\ll r_{\rm M}$ in the last equality. We note
that, the maximum impact parameter $b_{\rm max}$ in the context of Newtonian
gravity is usually assumed to be some (arbitrary) scale distance111Such cutoff
length in neutral plasma physics is unequivocally constrained by the Debye
length, see Spitzer (1965) of the system at hand. In the original derivation
of Chandrasekhar (1943a), where the system is infinite, the definition of
$b_{\rm max}$ remains somewhat unclear. Some authors (see e.g. Van Albada &
Szomoru 2020 and references therein) typically take $b_{\rm max}$ of the order
of the average inter-particle distance. Hereby, we follow the same assumption
noting that for a star of a solar mass the MOND radius $r_{\rm M}$ is roughly
$10^{3}$ AU, 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

$\displaystyle\left(\frac{d\mathbf{v}_{M}}{dt}\right)\approx$
$\displaystyle-16\pi^{2}\left(\ln\Lambda+\frac{b_{\rm max}^{2}}{2r_{\rm
M}^{2}}+2\frac{b_{\rm max}}{r_{\rm M}}\right)\times$ (26) $\displaystyle\times
G^{2}\rho(M+m)\frac{\@vec{v}_{M}}{v_{M}^{3}}\int_{0}^{v_{M}}v^{2}f(v)dv.$ (27)

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
$\nu$. Following this approach, the calculations as shown in Appendix A leads
to the same expression of Eq. (26). We note that Equation (26) differs from
what one would obtain using the simple derivation of $t_{2b}$ in Sect. 2 of
Ciotti & Binney 2004 that is proportional to the crossing time $t_{\rm
cross}=1/\sqrt{G\rho}$ to evaluate the dynamical friction time scale with the
Spitzer (1987) relation $t_{\rm DF}=2t_{2b}m/(m+M)$. This is because their
expression of $t_{2b}$ makes use of the force (7), assuming implicitly the
dMOND regime.
Unlike the usual MOND interpolating function $\mu(x)$, the form of the QuMOND
function $\nu(y)$ affects the behaviour of the MONDian force field. Adopting
the commonly used Eq. (13) (see Kroupa et al. 2022), the DF expression becomes

$\displaystyle\frac{d\mathbf{v}_{M}}{dt}\approx$
$\displaystyle-16\pi^{2}\left(\ln\Lambda+\frac{b_{\rm max}^{2}}{2r_{\rm
M}^{2}}+\frac{b_{\rm max}}{r_{\rm M}}-\frac{1}{2}\ln\frac{b_{\rm max}}{r_{\rm
M}}-\frac{2+\ln 2}{4}\right)\times$ (28) $\displaystyle\times
G^{2}\rho(M+m)\frac{\@vec{v}_{M}}{v_{M}^{3}}\int_{0}^{v_{M}}v^{2}f(v)dv.$ (29)

The explicit derivation is given in Appendix B below.
Noteworthy, we observe in both cases that the MOND correction is proportional
to a function of the only $b_{\rm max}/r_{\rm M}$

$\displaystyle\left(\frac{d\mathbf{v}_{M}}{dt}\right)_{\rm QuMOND}$
$\displaystyle-\left(\frac{d\mathbf{v}_{M}}{dt}\right)_{\rm
Newt}\propto\frac{b_{\rm max}^{2}}{2r_{\rm M}^{2}}+2\frac{b_{\rm max}}{r_{\rm
M}}\quad,$ (30) $\displaystyle\frac{b_{\rm max}^{2}}{2r_{\rm
M}^{2}}+\frac{b_{\rm max}}{r_{\rm M}}-\frac{1}{2}\ln\frac{b_{\rm max}}{r_{\rm
M}}-\frac{1}{2}.$ (31)

The same term $b_{\rm max}^{2}/2r_{\rm M}^{2}$ always dominates the MOND
correction of DF, regardless of the exact form of $\nu$, while the smaller
terms are affected by the particular choice. Recently some authors (see e.g.
Lelli et al. 2017) proposed an expression for $\nu$ 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.
(22) particularly cumbersome.

Figure 1: Ratio of the QuMOND correction factor $X$ to the standard Coulomb
logarithm as function of the stellar density. The dependence on the local
velocity dispersion $\sigma$ is color coded.

As an example, in Figure 1 we show for $M\approx 2m$ the ratio of the
corrective QuMOND term $X=b_{\rm max}^{2}/2r_{M}^{2}+b_{\rm
max}/r_{M}-1/2\ln(b_{\rm max}/r_{M})-(2+\ln 2)/4$ of Eq. (28) with the
standard Coulomb logarithm as function of the stellar density, where we have
assumed $b_{\rm max}=(\rho/M_{\odot})^{-1/3}$ and $b_{\rm
min}=G(M+m)/\sigma^{2}$, where $\sigma$ (colour-coded in Figure) is the local
velocity dispersion (see e.g. Binney & Tremaine 2008). It is evident that for
an isolated baryon only system the DF in QuMOND is augmented (up to a factor
$\approx 50$) at low density, comparable to the typical intra-galactic density
of $10^{-1}$ stars per pc3, 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 $\sigma$.

##  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.
Let us consider the usual system of $N$ equal mass $m$ (field) particles, with
coordinates $\mathbf{r}_{i}$ and momenta $\mathbf{p}_{i}=m\mathbf{v}_{i}$
described by the (time-dependent) phase-space distribution function222Note
that, $\mathcal{F}$ is not the one particle distribution function $f(r,v)$
entering the collisionless Boltzmann equation, but rather the phase-space
distribution of a $6N$ degrees of freedom Hamiltonian system entering the
Liouville equation. $\mathcal{F}(\mathbf{r},\mathbf{p};t)$; and a test
particle $M$ with coordinate $\mathbf{r}_{0}=\mathbf{R}$ and momentum
$\mathbf{p}_{0}=\mathbf{P}=M\mathbf{v}_{M}$ that perturbs the initial phase-
space distribution $\mathcal{F}_{0}$. In 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 assumes an external potential $\Phi$ confining
the system; see Appendix C below. We note that, MOND systems with self-
consistent gravitational field $g_{\rm in}$ embedded in an external
gravitational field $g_{\rm ext}$ are prone to the so-called external field
effect (hereafter EFE). The latter implies that for $g_{\rm in}<a_{0}<g_{\rm
ext}$ the system is purely Newtonian, while for $g_{\rm in}<g_{\rm ext}<a_{0}$
the system behaves as a Newtonian model with rescaled gravitational constant
$G^{\prime}=Ga_{0}/g_{\rm ext}$.
Under the assumptions given above, as usual, the average number density of
particles $n=N/{\rm Vol}$ should be taken constant as $N,{\rm
Vol}\rightarrow\infty$. The dynamical friction force $F_{0}^{\rm fr}$ can be
therefore obtained simply as $\langle F_{0}^{\rm 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 $\hat{\mathbf{r}},\hat{\mathbf{p}}$. The linear evolution
operator for $\mathcal{F}$ takes the form
$\dot{\mathcal{F}}(\hat{\mathbf{r}}_{i},\hat{\mathbf{p}}_{i};t)=-i\mathcal{L}[\mathcal{F}]$.
See Eqs. (67) and (69) 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 local333MOND
is local in the sense that the density at point $\xi$ equals a differential
operator acting on the potential in $\xi$. Vice versa, for the case of
$1/r^{\alpha}$ forces with $\alpha\neq 2$, obeying the superposition principle
this is not true (Stein 1970). In practice, for such general long-range
interactions, density and potential are related only via integral relation
over the whole domain occupied by the system (Di Cintio & Ciotti 2011; Di
Cintio et al. 2013). 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 one cannot substitute the potentials in a form

$W_{i}=\sum_{j\neq i}W_{ij}(|\mathbf{r}_{i}-\mathbf{r}_{j}|),$ (32)

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

$\displaystyle\mathbf{F}_{0}^{\rm fr}=-\beta\int
d\Sigma\mathcal{F}_{0}\left(\sum_{j=1}^{N}\mathbf{F}_{i}^{\rm
tot}\right)\int_{0}^{t}d\tau G_{\mathcal{L}}(\tau\rightarrow
t)\left[\sum_{i=1}^{N}\mathbf{v}_{i}\cdot\mathbf{F}_{i}^{\rm tot}\right],$
(33)

where $G_{\mathcal{L}}$ is the Greenian of the operator $\mathcal{L}$ and
$d\Sigma=d^{3N}\mathbf{r}d^{3N}\mathbf{p}$ 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. (33). Assuming that the system under consideration has a
finite memory (as implied by its ergodicity), we can replace the time
integration in $\tau$ on the finite interval $[0;t]$ with an integration
extended over the semi-infinite interval $[0;\infty)$. Moreover, making use of
the standard linear trajectory approximation (see e.g. Ter Haar 1977, see also
Syer 1994), allows to simplify $G_{\mathcal{L}}(\tau\rightarrow t)[Q]\cong
Q(t-\tau)$. 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 a non-linear theory.
Under the assumptions listed above, the DF force formula becomes

$\displaystyle M\frac{d\mathbf{v}_{M}}{dt}=\mathbf{F}_{0}^{\rm fr}\cong$
$\displaystyle-\frac{3(GM)^{2}\rho}{\langle v^{2}\rangle}\int
v_{\alpha}f(\mathbf{v})d^{3}\mathbf{v}\int_{0}^{\infty}d\tau$ (34)
$\displaystyle\times\int d^{3}\mathbf{s}\nu\left(\frac{r_{\rm
M}^{2}}{|\mathbf{s}-\tilde{\mathbf{v}}\tau|^{2}}\right)\frac{s^{\alpha}-\tilde{v}^{\alpha}\tau}{|\mathbf{s}-\tilde{\mathbf{v}}_{1}\tau|^{3}}\nu\left(\frac{r_{\rm
M}^{2}}{s^{2}}\right)\frac{\mathbf{s}}{s^{3}},$ (35)

Figure 2: Evolution of the galactocentric radius $r_{g}$ 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.

where $\tilde{\mathbf{v}}=\mathbf{v}-\mathbf{V}$. At this stage, the problem
is reduced to solving a three dimensional integral for a given choice of the
QuMOND interpolation function $\nu$.
Equation (34) in the Newtonian limit (i.e. $\nu\rightarrow 1$) can be solved
by an integration by parts applying

$\nabla\cdot\frac{\mathbf{s}}{s^{3}}=-4\pi\delta^{3}(\mathbf{s})$ (36)

and

$\mathbf{v}\cdot\frac{\mathbf{s}-\tilde{\mathbf{v}}\tau}{|\mathbf{s}-\tilde{\mathbf{v}}_{1}\tau|^{3}}=\nabla\cdot\frac{\mathbf{v}}{|\mathbf{s}-\tilde{\mathbf{v}}_{1}\tau|},$
(37)

so that one has

$\int\left(\mathbf{v}\cdot\frac{\mathbf{s}-\tilde{\mathbf{v}}\tau}{|\mathbf{s}-\tilde{\mathbf{v}}_{1}\tau|^{3}}\right)\frac{\mathbf{s}}{s^{3}}d^{3}\mathbf{s}=4\pi\frac{\mathbf{v}}{\tilde{v}\tau}.$
(38)

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

$\int_{0}^{\infty}\frac{d\tau}{\tau}=\int_{0}^{\infty}\frac{db}{b}$ (39)

returns the usual Coulomb logarithm (see e.g. Kalnajs 1972), while the
integral in $d^{3}\mathbf{v}$ is rewritten as in the Newtonian case, therefore
yielding Eq. (17).
The QuMOND case is, even for the simple form of $\nu$ given in Eq. (23), is
plagued by its intrinsic complexity. Analogously to the naive formula (22) we
notice the presence of two factors $\nu$. In this mean field formalism the
arguments for the two functions $\nu$ in (34) are however slightly different.
From such difference $-\tau\tilde{\mathbf{v}}$ one recovers the logarithmic
term in $\tau$, as well as any additional QuMOND terms. At variance with the
Newtonian case, which requires some cutoff only for the integral over $d\tau$,
in the QuMOND DF formula (34) the integral over $d^{3}\mathbf{s}$ is also
diverging. Therefore needing a cutoff at $b_{\rm max}$. This yields terms of
the form $b_{\rm max}/r_{\rm M}$ and $b_{\rm max}^{2}/r_{\rm M}^{2}$,
substantially in agreement with the simplified approach performed in Sect.
2.2.

##  4 Numerical experiments

To clarify the different behaviour of DF in the QuMOND and Newtonian gravity
model of a given star system, it is instructive perform a simple numerical
experiment in a set-up such that in Newtonian gravity with DM the friction on
the test object $M$ would be negligible while the gravitational field of the
parent system is in the MONDian regime $g\ll a_{0}$ over a broad interval of
radii.
We consider a $M=10^{2}M_{\odot}$ star cluster moving through a dwarf galaxy
with stellar mass $M_{*}=2\times 10^{5}M_{\odot}$ and scale radius $r_{s}=0.8$
kpc (parameters compatible to those of the Draco ultrafaint dwarf galaxy,
Mashchenko et al. 2006), with stellar distribution given by

$\rho_{*}(r)=\frac{3}{4\pi}\frac{M_{i}r_{s}}{(r+r_{s})^{4}},$ (40)

corresponding to a Dehnen (1993) $\gamma-$model for $\gamma=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 (40) plus
a DM halo with $\rho_{\rm DM}=\tilde{\rho}-\rho_{*}$ (cfr. Eq. 12). By
contrast, the DF force experienced by $M$ is different in the two paradigms
-Newtonian and MOND- being Eq. (17) with $\rho_{DM}+\rho_{*}$ in lieu of
$\rho$ in the first case444We recall that, at equipartition in a multimass
system (e.g. stars and elementary particle sized DM), when the test mass $M$
is much larger than the mean mass $\langle m\rangle$ the prefactor of Eq. (17)
is dominated by $M$ times the total mass density, see e.g. Ciotti (2021). and
Eq. (28) in the second.
Using the numerical approach discussed in Appendix D, see also Pasquato & Di
Cintio (2020); Di Cintio et al. (2020); Di Cintio & Casetti (2022), we
integrated different orbits under the effect of the same gravitational field
and the two different DF expressions, in 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||\Phi||}$). In Figure 2 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 pergalactic radius falling
from $\approx 1.63$ kpc down to $\approx 1.2$ kpc at 12 Gyr. This is
remarkably in agreement with the analytical estimates of Ciotti & Binney
(2004), that predicted a MOND inspiral time for an object considerably more
massive than a star to the central regions of a Draco-like system of less than
a Hubble time, much shorter than its Newtonian analog and argued that this
could be the reason behind the absence of star clusters in similar dwarf
galaxies.

##  5 Conclusions

We have 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 a typical stellar mass $r_{\rm M}$. Such
additional term becomes relevant when the maximum impact parameter
(proportional to the average interparticle distance, Spitzer 1987) $b_{\rm
max}\gg r_{\rm M}$, cfr. Eqs. (26,28). In practice, in a dense star cluster or
in galactic nuclei where the mean stellar density is such that the typical
inter-stellar distance is always smaller than $r_{\rm M}$, the DF is always
purely Newtonian. Vice versa, in low density systems, such as ultrafaint Dwarf
galaxies (Simon 2019) the $b_{\rm max}/r_{\rm 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 $\nu$. However, the strength of the enhancing term
does not vary significantly for such different forms of $\nu$.
In additions, we have extended the mean field treatment of DF, pioneered by
Kandrup (1983) in the Newtonian case to QuMOND. In this framework, alternative
to the fluctuation-based approach by Ciotti & Binney (2004), we recovered an
integral expression that contains the usual dependence on $\nu$.
Unfortunately, the explicit evaluation of the DF force, becomes rapidly
cumbersome (at variance with the simpler Newtonian case) even in the dMOND
limit, one is therefore forced to integrate Eq. (34) numerically. Some
specific cases will be presented and discussed elsewhere. Working out in
QuMOND the two-body relaxation time (and therefore the DF coefficient of Eq. 6
via the relation $t_{DF}=2t_{2b}m/M$) with the fluctuation approach in Fourier
space of Ciotti & Binney (2004), though in principle possible, as 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
$\tilde{\rho}$.
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). We stress the fact however that, in MOND one can not 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 Alessandrini et al. (2014); Arca-Sedda &
Capuzzo-Dolcetta (2014); Di Cintio & Casetti (2022). and therefore our simple
numerical estimates could possibly overestimate the MONDian DF. Using the
quasi-linear formulation of MOND could in principle allow to incorporate the
contribution of DF (and possibly density fluctuations on a scale smaller than
the particle-resolution) adding it before evaluating $\tilde{\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_{\rm int}<g_{\rm
ext}<a_{0}$, the dynamics essentially Newtonian with a rescaled $G$. In those
cases, the DF force would be given by Eq. (17) augmented by the multiplicative
factor $(a_{0}/g_{\rm ext})^{2}$, and thus a full mean field QuMOND treatment
including the DF expression discussed in Sect. 3 is valid only for a isolated
system.

###### Acknowledgements.

We express gratitude to Jan Pflamm-Altenburg, Pavel Kroupa and Luca Ciotti for
the useful discussions at an early stage of this work.

## References

  * Alessandrini et al. (2014) Alessandrini, E., Lanzoni, B., Miocchi, P., Ciotti, L., & Ferraro, F. R. 2014, ApJ, 795, 169
  * Arca-Sedda & Capuzzo-Dolcetta (2014) Arca-Sedda, M. & Capuzzo-Dolcetta, R. 2014, ApJ, 785, 51
  * Asencio et al. (2022) Asencio, E., Banik, I., Mieske, S., et al. 2022, MNRAS, 515, 2981
  * Banik et al. (2022) Banik, I., Thies, I., Truelove, R., et al. 2022, Monthly Notices of the Royal Astronomical Society, 513, 129–158
  * Bekenstein & Milgrom (1984) Bekenstein, J. & Milgrom, M. 1984, ApJ, 286, 7
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Brada & Milgrom (1999) Brada, R. & Milgrom, M. 1999, ApJ, 519, 590
  * Bugg (2015) Bugg, D. 2015, Canadian Journal of Physics, 93, 119–125
  * Bílek et al. (2021) Bílek, M., Zhao, H., Famaey, B., et al. 2021, Astronomy & Astrophysics, 653, A170
  * Chandrasekhar (1941a) Chandrasekhar, S. 1941a, ApJ, 93, 285
  * Chandrasekhar (1941b) Chandrasekhar, S. 1941b, ApJ, 93, 323
  * Chandrasekhar (1943a) Chandrasekhar, S. 1943a, ApJ, 97, 255
  * Chandrasekhar (1943b) Chandrasekhar, S. 1943b, ApJ, 97, 263
  * Chandrasekhar (1943c) Chandrasekhar, S. 1943c, ApJ, 98, 54
  * Chiari & Di Cintio (2023) Chiari, C. & Di Cintio, P. 2023, A&A, 677, A140
  * Ciotti (2021) Ciotti, L. 2021, Introduction to Stellar Dynamics
  * Ciotti & Binney (2004) Ciotti, L. & Binney, J. 2004, MNRAS, 351, 285
  * Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250
  * Di Cintio (2023) Di Cintio, P. 2023, Disruption of star clusters in Newtonian and MOND gravity
  * Di Cintio & Casetti (2022) Di Cintio, P. & Casetti, L. 2022, in Multi-Scale (Time and Mass) Dynamics of Space Objects, ed. A. Celletti, C. Galeş, C. Beaugé, & A. Lemaître, Vol. 364, 152–157
  * Di Cintio & Ciotti (2011) Di Cintio, P. & Ciotti, L. 2011, International Journal of Bifurcation and Chaos, 21, 2279
  * Di Cintio et al. (2013) Di Cintio, P., Ciotti, L., & Nipoti, C. 2013, MNRAS, 431, 3177
  * Di Cintio et al. (2020) Di Cintio, P., Ciotti, L., & Nipoti, C. 2020, in Proceedings of the IAU Symposia, Vol. 351, Star Clusters: From the Milky Way to the Early Universe, ed. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, 93–96
  * Gilbert (1968) Gilbert, I. H. 1968, ApJ, 152, 1043
  * Kalnajs (1972) Kalnajs, A. J. 1972, in Astrophysics and Space Science Library, Vol. 31, IAU Colloq. 10: Gravitational N-Body Problem, ed. M. Lecar, 13
  * Kandrup (1980) Kandrup, H. E. 1980, Phys. Rep, 63, 1
  * Kandrup (1983) Kandrup, H. E. 1983, Ap&SS, 97, 435
  * Kandrup et al. (2000) Kandrup, H. E., Pogorelov, I. V., & Sideris, I. V. 2000, MNRAS, 311, 719
  * Knebe & Gibson (2004) Knebe, A. & Gibson, B. K. 2004, MNRAS, 347, 1055
  * Kroupa et al. (2022) Kroupa, P., Jerabkova, T., Thies, I., et al. 2022, MNRAS, 517, 3613
  * Lelli et al. (2017) Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152
  * Malekjani et al. (2009) Malekjani, M., Rahvar, S., & Haghi, H. 2009, ApJ, 694, 1220
  * Mannella (2004) Mannella, R. 2004, Phys. Rev. E, 69, 041107
  * Mashchenko et al. (2006) Mashchenko, S., Sills, A., & Couchman, H. M. 2006, ApJ, 640, 252
  * McGaugh & Milgrom (2013) McGaugh, S. & Milgrom, M. 2013, ApJ, 766, 22
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
  * Milgrom (1986) Milgrom, M. 1986, ApJ, 302, 617
  * Milgrom (1994) Milgrom, M. 1994, ApJ, 429, 540
  * Milgrom (2010) Milgrom, M. 2010, MNRAS, 403, 886
  * Nagesh et al. (2022) Nagesh, S. T., Kroupa, P., Banik, I., et al. 2022, Monthly Notices of the Royal Astronomical Society, 519, 5128–5148
  * Nagesh et al. (2023) Nagesh, S. T., Kroupa, P., Banik, I., et al. 2023, MNRAS, 519, 5128
  * Nipoti et al. (2008) Nipoti, C., Ciotti, L., Binney, J., & Londrillo, P. 2008, MNRAS, 386, 2194
  * Nipoti et al. (2007a) Nipoti, C., Londrillo, P., & Ciotti, L. 2007a, ApJ, 660, 256
  * Nipoti et al. (2007b) Nipoti, C., Londrillo, P., & Ciotti, L. 2007b, MNRAS, 381, L104
  * Nipoti et al. (2007c) Nipoti, C., Londrillo, P., Zhao, H., & Ciotti, L. 2007c, MNRAS, 379, 597
  * Pasquato & Di Cintio (2020) Pasquato, M. & Di Cintio, P. 2020, A&A, 640, A79
  * Re & Di Cintio (2023) Re, F. & Di Cintio, P. 2023, A&A, 678, A110
  * Sánchez-Salcedo & Hernandez (2007) Sánchez-Salcedo, F. J. & Hernandez, X. 2007, ApJ, 667, 878
  * Sanders (2008) Sanders, R. H. 2008, MNRAS, 386, 1588
  * Sanders (2021) Sanders, R. H. 2021, Monthly Notices of the Royal Astronomical Society, 507, 803–808
  * Scarpa et al. (2022) Scarpa, R., Falomo, R., & Treves, A. 2022, MNRAS, 512, 544
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375
  * Sollima et al. (2023) Sollima, A., Nipoti, C., Calura, F., Pascale, R., & Baumgardt, H. 2023, MNRAS, 524, 5291
  * Spitzer (1965) Spitzer, L. 1965, Physics of fully ionized gases (Interscience Publishers)
  * Spitzer (1987) Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton University Press)
  * Stein (1970) Stein, E. M. 1970, Singular Integrals and Differentiability Properties of Functions (PMS-30) (Princeton University Press)
  * Syer (1994) Syer, D. 1994, MNRAS, 270, 205
  * Ter Haar (1977) Ter Haar, D. 1977, Lectures on selected topics in statistical mechanics.
  * Tiret & Combes (2007) Tiret, O. & Combes, F. 2007, A&A, 464, 517
  * Van Albada & Szomoru (2020) Van Albada, T. S. & Szomoru, A. 2020, Proceedings of the International Astronomical Union, 351, 532
  * Wu et al. (2009) Wu, X., Zhao, H., Wang, Y., Llinares, C., & Knebe, A. 2009, MNRAS, 396, 109
  * Zhu et al. (2023) Zhu, Y., Ma, H.-X., Dong, X.-B., et al. 2023, MNRAS, 519, 4479

##  Appendix A Effective Dark Matter halos

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

$\tilde{\rho}=-\frac{1}{4\pi G}\@vec{\nabla}\cdot\@vec{g}=-\frac{1}{4\pi
G}\@vec{\nabla}\cdot\left[\nu(|\<EMAIL_ADDRESS>(41)

For a point particle, the density becomes
$\rho(\@vec{r})=m\delta^{3}(\@vec{r})$, so that
$\@vec{g}_{N}(\@vec{r})=-{Gm}/{r^{3}}\@vec{r}$, and

$\tilde{\rho}=-\frac{1}{4\pi
G}\@vec{\nabla}\cdot\left[-\frac{Gm}{r^{3}}\@vec{r}\nu\left(\frac{Gm}{a_{0}r^{2}}\right)\right]=\frac{m}{4\pi}\@vec{\nabla}\cdot\left[\frac{\@vec{r}}{r^{3}}\nu\left(\frac{r_{\rm
M}^{2}}{r^{2}}\right)\right]$ (42)

Substituting the simple interpolation (23), we then find

$\tilde{\rho}=\frac{m}{4\pi}\@vec{\nabla}\cdot\left(\frac{\@vec{r}}{r^{3}}\right)+\frac{m}{4\pi
r_{\rm
M}}\@vec{\nabla}\cdot\left(\frac{\@vec{r}}{r^{2}}\right)=m\delta^{3}(\@vec{r})+\frac{m}{4\pi
r_{\rm M}r^{2}},$ (43)

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

$m_{DM}(r)=\int\frac{m}{4\pi r_{\rm M}r^{2}}4\pi r^{2}dr=m\frac{r}{r_{\rm
M}}.$ (44)

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

$F=\frac{GMm(b)}{b^{2}}=\frac{GMm}{b^{2}}+\frac{GMm}{r_{\rm M}b}.$ (45)

The associated momentum variation is $\Delta
p=Fb/v_{M}={GMm}/{v_{M}}\left(b^{-1}+r^{-1}_{\rm M}\right)$ and then the
diffusion coefficient becomes

$\displaystyle D_{v}$ $\displaystyle=\int\Delta p^{2}2\pi nv_{M}ddb=2\pi
G^{2}M^{2}\frac{m^{2}n}{v_{M}}\int\left(\frac{1}{b}+\frac{1}{r_{\rm
M}}\right)^{2}bdb.$ (46)

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

$\displaystyle\frac{d\mathbf{v}_{M}}{dt}$
$\displaystyle=-2\rho(M+m)D_{v}\frac{N(|\@vec{v}|<v_{M})}{v_{M}^{2}}\@vec{v}_{M}=$
(47) $\displaystyle=-16\pi^{2}\left(\ln\Lambda+2\frac{b_{\rm max}}{r_{\rm
M}}+\frac{b_{\rm max}^{2}}{2r_{\rm
M}^{2}}\right)G^{2}\rho(M+m)\frac{\@vec{v}_{M}}{v_{M}^{3}}\int_{0}^{v_{M}}v^{2}f(v)dv.$
(48)

##  Appendix B Explicit choice of QuMOND interpolating function

If the usual choice (13) of the QuMOND interpolating function $\nu$ is used in
Eq. (22), we have

$\nu(r_{\rm M}^{2}/b^{2})=\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{b^{2}}{r_{\rm
M}}}\Rightarrow\nu(b)^{2}=\frac{1}{2}+\frac{1}{2}\sqrt{1+4\frac{b^{2}}{r_{\rm
M}^{2}}}+\frac{b^{2}}{r_{\rm M}^{2}}.$ (50)

In the Equation above, unfortunately, we can not take the usual approximation
for small or large impact parameters $b$, since the integral runs from $b_{\rm
min}<r_{\rm M}$ to $b_{\rm max}>r_{\rm M}$. Vary large values of $b_{\rm 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\Lambda$.
This must also hold true _a fortiori_ in MOND, as in the far-field the forces
decay as $1/r$ rather than $1/r^{2}$. Therefore, we can not expand $\nu$ to
low order terms in $b/r_{\rm M}$.
We must then solve the exact integral $\int\nu(b)^{2}db/b$, where the square
root term yields

$\displaystyle\int\sqrt{1+4b^{2}/r_{\rm M}^{2}}\frac{db}{b}$
$\displaystyle=\int\frac{\sqrt{1+x^{2}}}{x}dx=$ (51)
$\displaystyle=\sqrt{1+x^{2}}-\text{arctanh}(\sqrt{1+x^{2}})+{\rm const}.=$
(52)
$\displaystyle=\sqrt{1+x^{2}}-\frac{1}{2}\ln\frac{\sqrt{x^{2}+1}+1}{\sqrt{x^{2}+1}-1}+\left({\rm
const}-\frac{\pi}{2}i\right),$ (53)

where $x=2b/r_{\rm M}$. Evaluating Eq. (51) for $x=2b_{\rm min}/r_{\rm
M}\rightarrow 0$ returns approximately

$\displaystyle(1+x^{2}/2)-\frac{1}{2}\ln\frac{1+1}{(1+x^{2}/2)-1}=\ln\frac{b_{\rm
min}}{r_{\rm M}}+1,$ (55)

while for $x=2b_{\rm max}/r_{\rm M}\rightarrow\infty$ becomes asymptotically

$\displaystyle\left(x+\frac{1}{2x}\right)-\frac{1}{2}\ln\frac{x+1}{x-1}=2\frac{b_{\rm
max}}{r_{\rm M}}.$ (56)

The integral (51) is

$\int\sqrt{1+4b^{2}/r_{\rm M}^{2}}\frac{db}{b}\sim 2\frac{b_{\rm max}}{r_{\rm
M}}-\ln\frac{b_{\rm min}}{r_{\rm M}}-1.$ (57)

The integral over the impact parameters is then evaluated as

$\displaystyle I_{b}$
$\displaystyle=\frac{1}{2}\int\frac{db}{b}+\frac{1}{2}\int\sqrt{1+4b^{2}/r_{\rm
M}^{2}}\frac{db}{b}+\frac{1}{r_{\rm M}^{2}}\int bdb=$ (58)
$\displaystyle\sim\frac{1}{2}\ln\frac{b_{\rm max}}{b_{\rm min}}+\frac{b_{\rm
max}}{r_{\rm M}}-\frac{1}{2}\ln\frac{b_{\rm min}}{r_{\rm
M}}-\frac{1}{2}+\frac{b_{\rm max}^{2}}{2r_{\rm M}^{2}}=$ (59)
$\displaystyle=\ln\Lambda+\frac{b_{\rm max}^{2}}{2r_{\rm M}^{2}}+\frac{b_{\rm
max}}{r_{\rm M}}-\frac{1}{2}\ln\frac{b_{\rm max}}{r_{\rm M}}-\frac{1}{2}.$
(60)

##  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 $\mathcal{F}(\mathbf{r},\mathbf{p};t)$ the average of a
given phase-space observable $Q$, taken with respect to said distribution
function is

$\langle
Q\rangle_{\mathcal{F}}(t):=\int\mathcal{F}(\mathbf{r},\mathbf{p};t)Q(\mathbf{r},\mathbf{p})d\Sigma.$
(61)

where $d\Sigma=d^{3N}\mathbf{r}d^{3N}\mathbf{p}$ is the differential phase-
space element.
The initial configuration $\mathcal{F}_{0}$ 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 $\mathcal{F}_{0}$ is subtracted

$F_{0}^{fr}=\langle F_{0}^{tot}\rangle_{\mathcal{F}}-\langle
F_{0}^{tot}\rangle_{\mathcal{F}_{0}},\quad
F_{0}^{tot}=-\sum_{i=1}^{N}F_{i}^{tot}.$ (62)

The field particles are assumed statistically uncorrelated in their initial
state
$\mathcal{F}_{0}(\mathbf{r},\mathbf{p})=\prod_{i=1}^{N}f(\mathbf{r}_{i},\mathbf{p}_{i})$.
We also assume their velocity distribution to be a Maxwellian

$f(\mathbf{r}_{i},\mathbf{p}_{i})=\mathcal{N}f_{1}(\mathbf{r}_{i})e^{-\beta\mathbf{p}_{i}^{2}/2m_{i}},$
(63)

with temperature $\beta$ and normalization $\mathcal{N}$. 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_{1}(\mathbf{r}_{i})=e^{-\beta m_{i}\Phi_{i}(\mathbf{r}_{i})}$ with an
efficient potential $\Phi_{i}$.
For our purposes, we will $f_{1}$ as the limit of an infinite, homogeneous
distribution of field particles. Such homogeneous distribution on a finite
volume $Vol\subseteq\mathbb{R}^{3}$ is given for
$f_{1}(\mathbf{r_{i}})=\frac{1}{|Vol|}\chi_{Vol}(\mathbf{r}_{i})$, i.e. when
$\Phi_{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|\rightarrow\infty$. 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

$H(\mathbf{r},\mathbf{p};t)=\sum_{i=1}^{N}\frac{\mathbf{p}_{i}^{2}}{2m_{i}}+\sum_{i=1}^{N}W_{i}(\mathbf{r}).$
(64)

As discussed in Sect. (3) one can not make the usual substitution
$W_{i}=\sum_{j\neq i}W_{ij}(|\mathbf{r}_{i}-\mathbf{r}_{j}|)$, with
$W_{ij}(r)=-Gm_{i}m_{j}/r$, as in MOND the Superposition Principle does not
hold. The potentials should be rather expressed in using the QuMOND formalism
(8), so that

$-\nabla
W_{i}=m_{i}\nu\left(\frac{|\mathbf{g}_{Ni}|}{a_{0}}\right)\mathbf{g}_{Ni},\quad
s.t.\quad-\sum_{j\neq i}\nabla W_{ij}=m_{j}\mathbf{g}_{Ni}.$ (65)

As a consequence, it is useless to define the forces between couples of
particles $\mathbf{F}(j\rightarrow
i)=-\partial_{\mathbf{r}_{i}}W_{ij}(|\mathbf{r}_{i}-\mathbf{r}_{j}|)=-\mathbf{F}(i\rightarrow
j)$, given that the total QuMOND force on the $i-$th particle is not simply
$\mathbf{F}_{i}^{tot}\neq\sum_{j\neq 0}\mathbf{F}(j\rightarrow i)$. The force
$\mathbf{F}_{i}^{tot}$ can be calculated instead, only as

$\mathbf{F}_{i}^{tot}=-\partial_{\mathbf{r}_{i}}W_{i}(\mathbf{r}),$ (66)

where $W_{i}(\mathbf{r})$ is given by the nonlinear QuMOND formula (65). The
evolution Equation is therefore expressed as

$\dot{\mathcal{F}}=-\sum_{i=1}^{N}\frac{\mathbf{p}_{i}}{m_{i}}\cdot\frac{\partial\mathcal{F}}{\partial\mathbf{r}_{i}}-\sum_{i=1}^{N}\mathbf{F}_{i}^{tot}\cdot\frac{\partial\mathcal{F}}{\partial\mathbf{p}_{i}}:=-iL[\mathcal{F}].$
(67)

In the equation above the sums arise from those in the Hamiltonian (64). In
practice, the meaning of Eq. (67) 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

$\hat{\mathbf{r}}_{i}:=\mathbf{r}_{i}-\mathbf{R}(t)+\mathbf{R}(0),\quad\hat{\mathbf{p}}_{i}:=\mathbf{p}_{i}-\mathbf{P}(t).$
(68)

The evolution operator for $\mathcal{F}$ hence takes the form

$\dot{\mathcal{F}}(\hat{\mathbf{r}}_{i},\hat{\mathbf{p}}_{i};t)=-iL[\mathcal{F}]+\sum_{i=1}^{N}\dot{\mathbf{P}}(t)\cdot\frac{\partial\mathcal{F}}{\partial\mathbf{p}_{i}}:=-i\mathcal{L}[\mathcal{F}],$
(69)

where $\dot{\mathbf{P}}$ acts as an apparent force on the $i$-th particle.
Except for this, in this frame the forces do not change
$\hat{F}_{i}^{tot}\equiv F_{i}^{tot}$, since even in MOND they depend only on
the respective particles positions. Moreover, it is useful to define the
perturbation $\tilde{\mathcal{F}}(t):=\mathcal{F}(t)-\mathcal{F}_{0}$ of the
distribution function so that

$\langle Q\rangle_{\mathcal{F}}=\langle Q\rangle_{\mathcal{F}_{0}}+\langle
Q\rangle_{\tilde{\mathcal{F}}}$ (70)

for any quantity $Q$. This will be used for evaluating Eq. (62), so that
$\mathbf{F}_{0}^{fr}=\langle\mathbf{F}_{0}^{tot}\rangle_{\tilde{\mathcal{F}}}$.
The distribution function $\tilde{\mathcal{F}}$ evolves as

$\displaystyle\dot{\tilde{\mathcal{F}}}=-i\mathcal{L}[\tilde{\mathcal{F}}]+\beta\mathcal{F}_{0}(\mathbf{r},\mathbf{p})\sum_{i=1}^{N}\left[\frac{\hat{\mathbf{p}}_{i}}{m_{i}}+\mathbf{V}(t)\right]\cdot\left[F_{i}^{tot}+m_{i}\frac{\partial\hat{\Phi}_{i}}{\partial\hat{\mathbf{r}}_{i}}(t)\right],$
(71)

where again $\hat{\mathbf{p}}_{i}/m_{i}+\mathbf{V}(t)=\mathbf{p}_{i}/m_{i}$,
and the term $\partial\hat{\Phi}_{i}/\partial\hat{\mathbf{r}}_{i}$ vanishes in
the infinite, homogeneous limit. In principle, keeping the contribution of
$\hat{\Phi}$ would allow one to investigate in a self-consistent fashion the
MOND EFE, see Sect. 3. At this point, Equation (71) can be rewritten in its
simpler form

$\dot{\tilde{\mathcal{F}}}=-i\mathcal{L}[\tilde{\mathcal{F}}]+\beta\mathcal{F}_{0}\sum_{i=1}^{N}\mathbf{p}_{i}\cdot\mathbf{g}_{i}^{MOND}.$
(72)

Its solution can thus be formally expressed as

$\tilde{\mathcal{F}}(t)=\beta\mathcal{F}_{0}\int_{0}^{t}d\tau
G_{\mathcal{L}}(\tau\rightarrow
t)\left[\sum_{i=1}^{N}\mathbf{p}_{i}\cdot\mathbf{g}_{i}\right]$ (73)

with the Greenian $G_{\mathcal{L}}$ of the operator $\mathcal{L}$. After all
the substitution we have a formal, general expression for the dynamical
friction on the test particle in the form

$\displaystyle\mathbf{F}_{0}^{fr}$ $\displaystyle=\int
d\Sigma\mathbf{F}_{0}^{tot}\tilde{\mathcal{F}}=\int
d\Sigma\mathbf{F}_{0}^{tot}\beta\mathcal{F}_{0}\int_{0}^{t}d\tau
G_{\mathcal{L}}(\tau\rightarrow
t)\left[\sum_{i=1}^{N}\mathbf{v}_{i}\cdot\mathbf{F}^{*}_{i}\right]=$ (74)
$\displaystyle=-\beta\int
d\Sigma\mathcal{F}_{0}\left(\sum_{j=1}^{N}\mathbf{F}_{i}^{tot}\right)\int_{0}^{t}d\tau
G_{\mathcal{L}}(\tau\rightarrow
t)\left[\sum_{i=1}^{N}\mathbf{v}_{i}\cdot\mathbf{F}^{*}_{i}\right],$ (75)

where
$\mathbf{F}^{*}_{i}=\mathbf{F}_{i}^{tot}+m_{i}\partial\hat{\Phi}_{i}/\partial\hat{\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 $\sum_{i=0}^{N}\mathbf{F}_{i}^{tot}=0$.
In order to perform the integrals appearing in Eqs. (74), 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 $\tau$ on $[0;t]$ with an integration of $[0;\infty)$.
Moreover, we make use of the standard linear trajectory approximation, which
allows to simplify $G_{\mathcal{L}}(\tau\rightarrow t)[Q]\cong Q(t-\tau)$.
Under these approximations, the DF force becomes

$\mathbf{F}_{0}^{fr}\cong-\beta\int
d\Sigma\mathcal{F}_{0}\left(\sum_{j=1}^{N}\mathbf{F}_{j}^{tot}(t)\right)\int_{0}^{\infty}d\tau\left(\sum_{i=1}^{N}\mathbf{v}_{i}(t-\tau)\cdot\mathbf{F}_{i}^{tot}(t-\tau)\right).$
(76)

This can be further simplified since $\langle F_{i}F_{j}\rangle=0$ for the
off-diagonal combinations $i\neq j$, since the field particles are by
definition statistically uncorrelated. We are then left with

$\displaystyle\mathbf{F}_{0}^{fr}$ $\displaystyle\cong-\beta\int
d\Sigma\mathcal{F}_{0}\sum_{i=1}^{N}\mathbf{F}_{i}^{tot}(t)\int_{0}^{\infty}d\tau\mathbf{v}_{i}(t-\tau)\cdot\mathbf{F}_{i}^{tot}(t-\tau)=$
(77)
$\displaystyle=-\beta\int_{0}^{\infty}d\tau\left\langle\sum_{i=1}^{N}\mathbf{F}_{i}^{tot}(t)\left(\mathbf{v}_{i}(t-\tau)\cdot\mathbf{F}_{i}^{tot}(t-\tau)\right)\right\rangle_{\mathcal{F}_{0}}=$
(78) $\displaystyle=-\beta\int_{0}^{\infty}d\tau
N\left\langle\mathbf{F}_{1}^{tot}(t)\left(\mathbf{v}_{1}(t-\tau)\cdot\mathbf{F}_{1}^{tot}(t-\tau)\right)\right\rangle_{\mathcal{F}_{0}}.$
(79)

Note that, the last step is justified only if the $N$ field particles have all
the same mass $m$, so that they are statistically equivalent in the initial
phase-space distribution. Averaging $\mathcal{F}_{0}$ for the $i=1$ particle,
practically means that only the $f(\mathbf{r}_{1},\mathbf{p}_{1})$
distribution function is relevant, hence

$\displaystyle\left\langle\mathbf{F}_{1}^{tot}\left(\mathbf{v}_{1}\cdot\mathbf{F}_{1}^{tot}\right)\right\rangle_{\mathcal{F}_{0}}$
$\displaystyle=\int
m\mathbf{g}_{1}^{M}\left(m\mathbf{v}_{1}\cdot\mathbf{g}_{1}^{M}\right)\prod_{i=1}^{N}f(\mathbf{r}_{i},\mathbf{p}_{i})d^{3N}\mathbf{r}d^{3N}\mathbf{p}=$
(80)
$\displaystyle=m^{2}\int\mathbf{g}_{1}^{M}\left(\mathbf{v}_{1}\cdot\mathbf{g}_{1}^{M}\right)f(\mathbf{r}_{1},\mathbf{p}_{1})d^{3}\mathbf{r}_{1}d^{3}\mathbf{p}_{1}=$
(81) $\displaystyle=m^{2}\int d^{3}\mathbf{r}_{1}\int
d^{3}\mathbf{v}_{1}\frac{n}{N}f(\mathbf{v}_{1})\mathbf{g}_{1}^{M}\left(\mathbf{v}_{1}\cdot\mathbf{g}_{1}^{M}\right)=$
(82) $\displaystyle=m^{2}\frac{n}{N}\left[\int
v_{\alpha}f(\mathbf{v})d^{3}\mathbf{v}\right]\left[\int g_{1}^{\alpha
M}\mathbf{g}_{1}^{M}d^{3}\mathbf{r}_{1}\right],$ (83)

where $f(\mathbf{v}_{1})=\mathcal{N}e^{-\beta mv_{1}^{2}}$ is again the
Maxwellian velocity distribution. The mean square speed is thus given by
$\langle v^{2}\rangle=3/m\beta$ that can be substituted to obtain

$\displaystyle\mathbf{F}_{0}^{fr}\cong$ $\displaystyle-\frac{3n}{m\langle
v^{2}\rangle}m^{2}\int_{0}^{\infty}d\tau\left[\int
v_{\alpha}f(\mathbf{v})d^{3}\mathbf{v}\right]\left[\int g_{1}^{\alpha
M}(t-\tau)\mathbf{g}_{1}^{M}(t)d^{3}\mathbf{r}_{1}\right]=$ (84)
$\displaystyle\cong$ $\displaystyle-\frac{3\rho G^{2}M^{2}}{\langle
v^{2}\rangle}\int_{0}^{\infty}d\tau\int
v_{\alpha}f(\mathbf{v})d^{3}\mathbf{v}$ (85) $\displaystyle\times\int
d^{3}\mathbf{s}\nu\left(\frac{r_{\rm
M}^{2}}{|\mathbf{s}-\tilde{\mathbf{v}}\tau|^{2}}\right)\frac{s^{\alpha}-\tilde{v}^{\alpha}\tau}{|\mathbf{s}-\tilde{\mathbf{v}}_{1}\tau|^{3}}\nu\left(\frac{r_{\rm
M}^{2}}{s^{2}}\right)\frac{\mathbf{s}}{s^{3}},$ (86)

since
$\mathbf{g}_{1}^{M}(\mathbf{s},t)=GM\nu\left(GM/a_{0}s^{2}\right){\mathbf{s}}/{s^{3}}$
and $\mathbf{s}(t-\tau)\cong\mathbf{s}-\tilde{\mathbf{v}}_{1}\tau$, being
$\tilde{\mathbf{v}}=\mathbf{v}-\mathbf{V}$.

##  Appendix D Numerical methods

The dynamics of a particle confined by a static potential $\Phi$ under the
combined effect of force fluctuations and friction is given by the Langevin
type equations

$\frac{{\rm d}^{2}\mathbf{r}}{{\rm
d}t^{2}}=-\nabla\Phi_{\mathbf{r}}-\nu_{\mathbf{r},\mathbf{v}}\mathbf{v}+\delta
F;\quad\mathbf{v}=\frac{{\rm d}\mathbf{r}}{{\rm d}t}.$ (87)

In the most general case the DF coefficient $\nu$ depends explicitly on the
position and velocity through the phase-space distribution of the embedding
system (also sourcing the potential $\Phi$) and $\delta F$ is a fluctuating
force per unit mass with distribution and amplitude connected to $\nu$ via a
fluctuation-dissipation relation. For the systems of interest (Newtonian or
MOND) we solve Eq. (87) using the Mannella (2004) scheme

$\displaystyle x^{\prime}=x(t+\Delta t/2)=x(t)+\frac{\Delta t}{2}v(t)$
$\displaystyle v(t+\Delta t)=c_{2}\left[c_{1}v(t)+\Delta
t\nabla\Phi(x^{\prime})+d_{1}\tilde{F}(x^{\prime})\right]$ $\displaystyle
x(t+\Delta t)=x^{\prime}+\frac{\Delta t}{2}v(t+\Delta t),$ (88)

written here for simplicity for a 1D system and fixed time step $\Delta t$,
where the coefficients $c_{1},c_{2}$ and $d_{1}$ are given by

$c_{1}=1-\frac{\eta\Delta t}{2};\quad c_{2}=\frac{1}{1+\eta\Delta t/2};\quad
d_{1}=\sqrt{2\zeta\eta\Delta t}.$ (89)

In the equations above $\tilde{F}$ is sampled from a norm 1 Gaussian and
$\zeta$ in the case of a delta correlated noise is fixed by the standard
deviation of the distribution of $F$ as

$\langle F(x,t)F(x,t^{\prime})\rangle=2\eta\zeta\delta(t-t^{\prime}).$ (90)

In practice, if such distribution is unknown one takes, following Kandrup et
al. (2000) one could assume $\zeta=\mathcal{E}$ where $\mathcal{E}$ is the
instantaneous relative (positive) energy per unit mass of the test particle
along the orbit. In the simulations discussed in this paper the contribution
of the fluctuating term is artificially set to 0 as it would be negligible for
the type of gravitational systems under consideration. Note that, if
$\nu=\zeta=0$ Eqs. (D) simply become the standard Verlet second order scheme
in the drift-kick-drift form.