© The Authors 2025

1. Introduction

Dynamical friction (hereafter DF, Chandrasekhar 1943a,b,c) is the process that drives a massive particle of mass m_t, for example a compact object (black hole, neutron star) or a satellite dwarf galaxy, to sink into the host stellar system via multiple dynamical collisions with its constituting field particles of (mean) mass m. In its original formulation, evaluated for an infinitely extended system of constant number and mass density, n, and ρ = mn with isotropic velocity distribution f(v), DF amounts to a velocity dependent drag term of the form

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-16{\pi }^2{G}^2\rho (m_t+m)ln\mathrm{\Lambda }\frac{\mathrm{\Xi }(v)}{v^3}\mathbf{v}.\end{array}$$(1)

In the equation above, lnΛ is the (velocity averaged) Coulomb logarithm of the maximum to minimum impact parameters, Λ = b_max/b_min, in the hyperbolic scattering with 1/r² forces, v = ||v||, and

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

is the fractional velocity volume function (e.g. see Ciotti 2021) that counts the portion of field particles moving slower than m_t and contributing to DF.

Since its inception, Equation (1) has been applied with success to a wide plethora of astrophysical system ranging from massive stars in globular clusters (GCs) to satellite galaxies in dark matter (DM) halos. However, numerical experiments have so far revealed that a naive correction to the Chandrasekhar expression, whereby essentially one substitutes the local values of ρ(r) and f(v) in Eq. (1), might lead to either an overestimation or underestimation of DF. For example, in a series of papers also accounting for the feedback of m_t on the host system by Bertin et al. (2003), Arena et al. (2006), Arena & Bertin (2007), it has been shown that the classical idealized DF formulation works rather well in models with broad density cores, while it fails (i.e. predicts shorter sinking timescales) when the density profile is highly concentrated at small radii. In addition, in the latter case the feedback on the density distribution is less and less important the higher its initial concentration. Conversely, Read et al. (2006) and Goerdt et al. (2006, 2010) observed that massive objects moving on initially circular orbits in cored models first suffer a stronger DF, with respect to the Chandrasekhar estimate – often referred to as super-Chandrasekhar – to then ‘stall’ at a radius, r_s, which appears to depend on the specific details of the galaxy model and on m_t. This apparent suppression of DF¹ in stellar (or DM) cores, frequently dubbed core-stalling, has also been reported among others by Zelnikov & Kuskov (2016) for satellites sinking in flat galactic cores, by Kaur & Sridhar (2018), Dutta Chowdhury et al. (2019), and Cole et al. (2012), Bílek et al. (2024) in simulations of GCs orbiting dwarf galaxies, and by Genina et al. (2024) in simulations of supermassive black holes (SMBHs) in nuclear star clusters.

Banik & van den Bosch (2021, 2022) (see also Dattathri et al. 2025) interpret the core-stalling (at least in very weakly collisional systems) as an effect of the Lynden-Bell & Kalnajs (1972) torque generalized such that its original secular and adiabatic assumptions are relaxed to account for the contribution of the non-resonant orbits to the slowing down of the massive test particle. In practice once the test particle reaches the central core, the structure of the relative orbits of the near-resonant field particles, m, with respect to m_t (e.g. see Inoue 2011) changes significantly. Since the phase space distribution in nearly harmonic cores is characterized by rather gentle gradients, the net torque from said orbits is statistically almost always positive, thereby countering the friction (an effect labeled as dynamical buoyancy). The core-stalling is therefore happening once the torque induced buoyancy balances DF. Remarkably, another resonance-based argument (see Tremaine & Weinberg 1984), which was sketched out by Read et al. (2006) but later ruled out by Petts et al. (2015, 2016), interpreted the so-called super-Chandrasekhar regime observed in some N-body simulations as a consequence of using a local Maxwellian approximation to evaluate Ξ(v), instead of integrating over the phase-space distribution, f, which supports the specific density profile of the host model (see Appendix A, see also Kazantzidis et al. 2004).

For the tightly related problem of a heavy particle (BH) initially sitting at rest at the centre of a flat-cored system then migrating outwards due to the combined action of force fluctuations and torques², N-body simulations (e.g. Chatterjee et al. 2002a,b, 2003; Brockamp et al. 2011; Di Cintio et al. 2023) revealed a strong resolution dependence of the wander radius, r_wan (Bahcall & Wolf 1976), loosely defined as the typical displacement of m_t from the centre of the potential well, in the form

$$\begin{array}{c}\hfill r_{\mathrm{wan}}\approx r_c\sqrt{\frac{m}{m_t}}=r_c\sqrt{\frac{M}{N{m}_t}},\end{array}$$(3)

where r_c is the typical scale length of the model, M its total mass, and N the number of stars (or simulation particles in a N-body experiment). Another radial scale considered in BH dynamics that is usually compared with r_wan is the so-called influence radius (i.e. the radius below which the potential, Φ_t = −Gm_t/r, exerted by m_t dominates over the contribution of the host stellar system, e.g. see Peebles 1972), defined as

$$\begin{array}{c}\hfill r_{\inf }=\frac{G{m}_t}{{\sigma }^2}\approx r_c\frac{m_t}{M}=r_c\frac{m_t}{\mathit{Nm}},\end{array}$$(4)

where σ is the velocity dispersion. This established, it remains unclear whether, at least in some regimes of N and m_t, the core stalling is a discreteness effect linked to the rapid and incoherent force fluctuations in dense cores, or a resolution-independent process related to the instantaneous fraction of orbits with frequencies nearly resonant with the test particle m_t. In this work we explore this matter further by means of variable resolution N-body simulations and single particle propagation in smooth potentials with tunable DF and mean force fluctuations.

2. Setting the stage

2.1. Galaxy models

In most of the simulations discussed in this paper, we employ the frequently used γ models (Dehnen 1993; Tremaine et al. 1994) with a density distribution defined by

$$\begin{array}{c}\hfill \rho (r)=\frac{3-\gamma }{4\pi }\frac{M{r}_c}{r^{\gamma }{(r+r_c)}^{4-\gamma }};0\le \gamma <3,\end{array}$$(5)

generating the gravitational potential

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=\frac{\mathit{GM}}{r_c}\times \{\begin{array}{c}{\displaystyle \frac{1}{(2-\gamma )}[{\left(\frac{r}{r+r_c}\right)}^{2-\gamma }-1]\mathrm{for}\gamma \ne 2.}\hfill \\ {\displaystyle ln\frac{r}{r+r_c}\mathrm{for}\gamma =2.}\hfill \end{array}\end{array}$$(6)

The parameter γ appearing in Eqs. (5) and (6) above controls the central logarithmic density slope, ranging from 0 (corresponding to a flat core) to 3 (corresponding to a singularity). Here we concentrate on the cored γ = 0 and mildly cuspy γ = 1 cases. When generating the simulations initial conditions, the radial coordinates, r_i, of the N field particles are obtained sampling ρ with the usual Monte Carlo technique applied to the radial mass profile,

$$\begin{array}{c}\hfill M(r)=M{\left(\frac{r}{r+r_c}\right)}^{3-\gamma },\end{array}$$(7)

while their angular positions, ϑ_i and φ_i, are assigned randomly from two independent uniform distributions. The modulus of the velocities, v_i, is assigned using the rejection method (for example, see Di Cintio et al. 2015 and references therein) to sample the isotropic equilibrium phase-space distribution,

$$\begin{array}{c}\hfill f(\mathcal{E})=\frac{(3-\gamma )M}{2{(2{\pi }^{2}GM{r}_c)}^{3/2}}{\displaystyle {\int }_0^{\mathcal{E}}}\frac{{(1-y)}^2}{y^{4-\gamma }\sqrt{\mathcal{E}-\mathrm{\Psi }}}\mathrm{d}\mathrm{\Psi },\end{array}$$(8)

where Ψ = −Φ and ℰ = Ψ − v²/2 are the relative potential and energy, respectively, and

$$\begin{array}{c}\hfill y\left(\mathrm{\Psi }\right)=\{\begin{array}{c}{\displaystyle {[1-(2-\gamma )\mathrm{\Psi }]}^{1/(2-\gamma )},\gamma \ne 2.}\hfill \\ {\displaystyle exp(-\mathrm{\Psi }),\gamma =2.}\hfill \end{array}\end{array}$$(9)

In addition, to compare with previous numerical work, we have also considered the simpler Plummer (1911) model defined by the density potential pair

$$\begin{array}{c}\hfill \rho (r)=\frac{3M}{4\pi r_c^3}{(1+\frac{r^2}{r_c^2})}^{-5/2};\mathrm{\Phi }(r)=-\frac{\mathit{GM}}{\sqrt{r^2+r_c^2}},\end{array}$$(10)

with the radial mass profile

$$\begin{array}{c}\hfill M(r)=M{\left(\frac{r^2}{r^2+r_c^2}\right)}^{3/2}.\end{array}$$(11)

The isotropic phase-space distribution function associated with the model (10) is

$$\begin{array}{c}\hfill f(\mathcal{E})=\frac{\sqrt{2}}{378{\pi }^{2}G{r}_c^2{\sigma }_0}{\left(\frac{\mathcal{E}}{{\sigma }_0^2}\right)}^{7/2},\end{array}$$(12)

where ${\sigma }_0=\sqrt{GM/6r_c}$ is the central velocity dispersion. We note that the radial gravitational accelerations of a_r = −GM(r)/r² in the Plummer model and the γ = 0 model are both nearly harmonic near their centres; however, for the same choice of M and r_c, the radial span over which a_r ∝ r is larger for the Plummer density profile.

2.2. N-body method and simulation set-up

We performed N-body simulations in which a test mass, m_t, sinked in an equilibrium system represented by N equal mass m = M/N field particles, with 10⁴ ≤ N ≤ 3 × 10⁶. The test particle mass values used in this work range from 2 × 10⁻⁴ to 10⁻² in units of the total mass of the host system, M, while the orbit of m_t is initialized with values of the eccentricity, e, ranging from 0 (i.e. circular orbit) to 0.9, corresponding to a relatively eccentric orbit. As usual e is defined as

$$\begin{array}{c}\hfill e=\frac{r_{\mathrm{apo}}-r_{\mathrm{peri}}}{r_{\mathrm{apo}}+r_{\mathrm{peri}}},\end{array}$$(13)

where r_apo and r_peri are the apocentre and pericentre of the unperturbed orbit in the smooth potential, Φ.

In our numerical experiments the field particles either interact with each other, with forces computed with the fast Dehnen (2000) tree scheme implemented in the FVFPS code (see Londrillo et al. 2003; Nipoti et al. 2006), or move as independent tracers in the smooth potential given by Eqs. (6), (10) in a similar fashion to White (1983) and Bontekoe & van Albada (1987). In this latter case, we aim to neglect the higher-order effects on DF induced by the formation of a self-gravitating wake (Mulder 1983; Kaur & Stone 2022) behind m_t, sometimes also referred to as a polarization cloud (Gilbert 1970; Kalnajs 1972; Heyvaerts et al. 2017). In both types of simulation, each field particle of mass m interacts with the test particle, m_t, with the gravitational force evaluated with a direct sum scheme. To prevent unwanted resolution-induced force singularities, the Newtonian 1/r² interaction is substituted with a quadratic spline below the softening length, ϵ. In the results discussed hereafter, the fixed choice was ϵ = 2 × 10⁻³r_c. We recall that in direct and tree codes the softening acts as a sort of spatial resolution. We recall also that the introduction of said softening effectively dictates a fixed minimum impact parameter in two body-encounters (see Huang et al. 1993; Chiron & Marcos 2019). All particles are propagated with a second-order leapfrog scheme with a fixed time step, Δt = t_c/150, where the reference timescale is $t_c=\sqrt{r_c^3/GM}$. We adopted units such that G = M = r_c = t_c = 1.

2.3. Stochastic equation of motion models

For reference we also performed simulations where m_t is propagated in the smooth potential under the additional effect of DF and force fluctuations modeling the ‘granularity’ of the system, in principle with an arbitrary resolution, without resulting to integrating the field particle orbits. In practice, following the approach of Sartorello et al. (2025) we assume that the dynamics of m_t is statistically equivalent to that induced by the stochastic differential Langevin equation (see Chandrasekhar 1943d; Kandrup 1980)

$$\begin{array}{c}\hfill \ddot{\mathbf{r}}=-\mathrm{\nabla }\mathrm{\Phi }(\mathbf{r})-\eta \mathbf{v}+\delta \mathbf{f}(\mathbf{r});\mathbf{v}=\dot{\mathbf{r}},\end{array}$$(14)

where ηv is essentially the DF term given in Eq. (1) and δf is a stochastic force (per unit mass) fluctuation. The formulation of the problem given by Eq. (14) is equivalent to the evolution equation for the particle energy fluctuations ⟨dE/dt⟩ derived by Bekenstein & Maoz (1992) and Maoz (1993) for homogeneous or inhomogeneous gravitational media using a fluctuation-dissipation approach.

In the case of an infinite and homogeneous system of equal mass particles interacting via 1/r² forces subjected to Poissonian density fluctuations, δf is distributed according (see Chandrasekhar & von Neumann 1942, 1943) to the Holtsmark (1919) distribution, originally evaluated for an infinite gas of charged particles. The Holtsmark distribution of the modulus of the force fluctuations, δf = ||δf||, reads in integral form as a function of the dimensionless auxiliary variable, s, as

$$\begin{array}{c}\hfill \mathcal{F}(\delta f)=\frac{2}{\pi \delta f}{\displaystyle {\int }_0^{\infty }}exp[-\alpha {(s/\delta f)}^{3/2}]ssin(s)\mathrm{d}s,\end{array}$$(15)

while its density dependent normalization factor, α, is defined as

$$\begin{array}{c}\hfill \alpha =\frac{4}{15}{(2\pi G)}^{3/2}{m}^{1/2}\rho .\end{array}$$(16)

Consistently with previous work, we accommodated the definition (15) in the γ models (5) and Plummer model (10) using the local values of ρ (e.g. see Petrovskaya 1986) and used the typical local Maxwellian approximation³ for the velocity distribution of the model in the Chandrasekhar formula so that

$$\begin{array}{c}\hfill \mathrm{\Xi }(v)={(4\pi )}^{-1}[\mathrm{Erf}\left(\frac{v}{\sqrt{2}\sigma }\right)-\sqrt{\frac{2}{\pi }}\frac{v}{\sigma }exp(-\frac{v^2}{2{\sigma }^2})].\end{array}$$(17)

For the γ models, the radial velocity dispersion profile is given in an integral form as

$$\begin{array}{c}\hfill {\sigma }^2(r)=GM{r}^{\gamma }{(r+r_c)}^{4-\gamma }{\displaystyle {\int }_r^{\infty }}\frac{r^{\prime 1-2\gamma }\mathrm{d}{r}^{\prime }}{{(r^{\prime }+r_c)}^{7-2\gamma }}·\end{array}$$(18)

In the two cases of interest explored here, γ = 0 and 1, we have

$$\begin{array}{c}\hfill {\sigma }^2(r)=\frac{\mathit{GM}}{30}\frac{r_c+6r}{{(r_c+r)}^2}\end{array}$$(19)

and (see also Hernquist 1990)

$$\begin{array}{cc}\hfill {\sigma }^2(r)& =\frac{\mathit{GM}}{12r_c}[\frac{12r{(r+r_c)}^3}{r_c^4}ln(\frac{r+r_c}{r})+\hfill \\ \hfill & -\frac{r}{r+r_c}(25+\frac{52r}{r_c}+\frac{42r^2}{r_c^2}+\frac{12r^3}{r_c^3})],\hfill \end{array}$$(20)

respectively, while for the Plummer sphere we have

$$\begin{array}{c}\hfill {\sigma }^2(r)=\frac{\mathit{GM}}{6\sqrt{r^2+r_c^2}}·\end{array}$$(21)

The Coulomb logarithm appearing in Eq. (1) was evaluated accounting for the exact integration over the impact parameter, b, in hyperbolic orbits (see e.g. Binney & Tremaine 2008) as

$$\begin{array}{c}\hfill log[1+(\frac{b_{\max }}{b_{\min }}{)}^2]=log\{1+[\frac{r_{\ast }{\sigma }^2}{G(m_t+m)}{]}^2\}\approx log\mathrm{\Lambda }.\end{array}$$(22)

The choice of minimum impact parameter, b_min, was made as usual to account for the angular momentum conservation in a φ = π/2 deflection at relative velocity of the order of the local velocity dispersion, σ, while b_max was set to r_*, a typical scale length of the order of the average inter-particle distance, ⟨r⟩=(8πmr_c³/3M)^1/3 ∝ r_c/N^1/3. We stress the fact that there is still not much of a consensus around the ‘correct’ choice of b_max (see Hamilton et al. 2018), with some authors assuming the scale length of the host, r_c, and others (e.g. Kandrup 1980, 1983) preferring the local average interparticle distance, $\overline{r}={(m/\rho )}^{1/3}$ (see the discussion in van Albada & Szomoru 2020). Another choice could be the distance at which the correlations in the fluctuations of stellar density, velocity, and gravitational potential induced by m_t vanish, in a similar fashion to the Debye-Hukel theory in neutral plasmas (see Marochnik 1968). Ad hoc direct N-body simulations by Marcos et al. (2017), in which b_min can be artificially fixed by the softening parameter, ϵ, seem to indicate that b_max ≈ r_c/3. Here we note that, due to their rather low highest resolution (i.e. N ≈ 1.7 × 10⁴), setting the maximum impact parameter to either r_c/3 or ⟨r⟩ inside the logarithm has little effect on the magnitude of the DF coefficient for resolutions below roughly ∼5 × 10⁶.

The stochastic equation of motion (14) was solved by computing the DF term at all times as a function of the local density and velocity distribution of the background model, and sampling the modulus of the force fluctuations from Eq. (15) with a standard Monte Carlo. The direction of the fluctuations was assumed to be isotropic (i.e. δf was given a random direction by sampling a uniform distribution of angles). We integrated Eq. (14) with the generalized leapfrog scheme implemented in Pasquato & Di Cintio (2020), Di Cintio et al. (2020), Sartorello et al. (2025) (for additional details and the original implementation see also Mannella 2004; Burrage et al. 2007). The benefits of applying this integration technique are essentially the fact that it converges on the standard second-order Verlet scheme for vanishing noise and friction, and its rather fast convergence on asymptotic distributions in the case of simple Brownian motion. Other higher-order and time-adaptive schemes, for example Runge-Kutta or Bulrisch-Stoer (see Burrage et al. 2007), require longer computational times, as is dictated by a smaller adaptive Δt, to obtain a comparable accuracy. It must be pointed out that the method introduced by Mannella (2004) in our implementation explicitly assumes a delta-correlated noise (i.e. subsequent encounters are independent); however, strictly speaking, it must be noted that the space-integrated force fluctuations are at zero-sum altogether, since δf is related to density fluctuations, δρ, through the usual Poisson equation. At constant total mass, all fluctuations must integrate out to zero. This implies that force fluctuations cannot be completely uncorrelated, making the delta-correlated noise a rough approximation that could be, in principle, sharpened by introducing a correlation time, τ_c (see e.g. Trani & Di Cintio 2025). Other authors, such as Pogorelov & Kandrup (1999) and Terzic & Kandrup (2003), used different forms of the correlation function of gravitational fluctuations, without finding significant qualitative differences with the delta-correlated noise in the case of strongly interacting systems. In some runs, in which the δf term was artificially set to zero, we instead used a modified mid-point leapfrog, typically adopted for the integration of orbits with velocity-dependent force terms (see for example Mikkola & Merritt 2006).

2.4. Recovering the distribution of force fluctuations

Using a locally tailored Holtsmark distribution in Eq. (14) might seem a rather crude approximation. In order to check the reliability of this assumption, inspired by previous work by Ahmad & Cohen (1973, 1974) and Del Popolo (1996a,b) we studied the numerically evaluated distribution of the force (per unit mass) fluctuations for some discrete N-body realization of the models discussed in the previous section.

We defined the force fluctuation at the position r_i of the ith particle in a discrete realization of the parent continuum system as

$$\begin{array}{c}\hfill \delta \mathbf{f}({\mathbf{r}}_i)=-\mathrm{\nabla }\mathrm{\Phi }({\mathbf{r}}_i)+Gm{\displaystyle \sum _{i\ne j=1}^N}\frac{{\mathbf{r}}_{\mathbf{i}}-{\mathbf{r}}_{\mathbf{j}}}{||{\mathbf{r}}_i-{\mathbf{r}}_j{||}^3},\end{array}$$(23)

where Φ is the smooth potential given by Eq. (5) or Eq. (10), with normalization such that M = Nm. We then partitioned the model in N_s logarithmically spaced radial shells and evaluated for each of them the probability density function, 𝒫, of the modulus of the force fluctuations, δf = ||δf||, by particle counts as well as its average value, ⟨δf⟩, and variance, $\overline{\delta f}=⟨\delta f^2⟩-{⟨\delta f⟩}^2$. In Fig. 1 we show ⟨δf⟩ (left) and $\overline{\delta f}$ (right) for a Plummer model with N = 10⁶ for N_s = 20 radial shells as a function of their radius, r. It appears that both the average and the variance of the force fluctuations grow for decreasing r, being almost constant inside the core (i.e. below r_c). The associated histograms of 𝒫(δf) are shown in Fig. 2 for 100 independent realizations (grey dots) and their average (blue dots), together with the Holtsmark distribution (red lines) computed evaluating numerically the integral in Eq. (15) with the choice of α where m = M/N and ρ = ρ(r) for the specific model at hand, in this case given by Eq. (10). Remarkably, the empirical force distribution⁴ appears to agree rather well with the Holtsmark distribution, in particular for the shells enclosing larger fractions of N for which the different realizations have a less pronounced scatter -the average number of field particles, $\overline{N}$, in each shell is indicated at the bottom of the panels, so that we can safely assume that the stochastic single particle simulations can be reasonably compared to honest direct N-body simulations, at least for what concerns the form of the fluctuation distribution.

Fig. 1. Radial profile of the average of the discrete force fluctuations (left panel) and its variance (right panel) evaluated for a single realization of a Plummer system with N = 106.

Fig. 2. Numerically evaluated distribution of force fluctuation, 𝒫(δf), in different shells of increasing radius (clockwise from top left, r = 3.1 × 10−2, 7.2 × 10−2, 0.38, 1.3, 6.9, and 54.8 in units of rc) for a Plummer model with N = 106 and 102 realizations (grey dots). The blue dots mark the PDF averaged over said ensemble, while the red curve is the corresponding the Holtsmark distribution (15) computed for the local value of the density, ρ(r), and average number of particles in the shell, $\overline{N}$.

3. Simulations and results

3.1. Preliminary tests

As a first set of numerical experiments, we compared the orbital decay of the test mass, m_t, in self-consistent N-body simulations and non-interacting tracer systems, for different resolutions and initial orbital eccentricities. We find that, in both cases, independently of the specific model or details of the orbit, the core-stalling can always be observed to some degree. For a given combination of initial orbital parameters and galaxy model, the same mass, m_t, systematically reaches a slightly smaller stalling radius, r_s, in the simulation, where it is slowed down by non-interacting tracers. In the simulations discussed hereafter, r_s is estimated by time-averaging the radial position of m_t over the last 10t_c of the simulation. Additionally we have also evaluated the PDF of the radial position over the same time span and extracted its mean and standard deviation, finding a rather good agreement with the time average. Overall, the orbital decay is qualitatively comparable across the two types of runs, first affecting only r_a. As an example, in Fig. 3 we show two test orbits with initial eccentricities of e = 0.001 and 0.3, for a particle of mass m_t = 10⁻³M propagated in γ = 0 and γ = 1 models (left and right panels, respectively). Both orbits plotted here stall at a r_s larger than the softening length, ϵ = 2 × 10⁻³r_c, used in the simulation, here indicated by the horizontal dashed line. Typically, the orbit begins to stall (i.e. the decrease in the pericentre becomes faster than that of the apocentre) once the enclosed mass fraction, M(r), becomes of the same order as m_t, as is indicated by the dotted black lines. The critical radius, r_mt, where this happens encloses a systematically smaller number of particles for decreasing N. Therefore for the low-resolution cases it may contain only order 10 particles. A further comparison between the two sets of numerical experiments is discussed in Appendix B.

Fig. 3. Radial decay of a massive tracer, mt = 10−3M, orbiting in a γ = 0 (left) and a γ = 1 (right) model in a self consistent N-body run (blue curves) and a simplified simulation where the same N = 106 field particles move independently in the static smooth potential (magenta curves). The horizontal dashed orange line marks the value of the softening parameter ϵ = 2 × 10−3rc, while the dotted black line indicates the radius, rmt, at which M(r) = mt, equal to ≈0.1 for the simulations with γ = 0 and ≈0.033 for those with γ = 1.

As a means to explore the feedback of m_t on the parent system, in Fig. 4 we plot the radial density profile, ρ(r) at t = 10³t_c (triangles), against the initial one (squares), corresponding to the simulations shown in Fig. 3. In the case of the cored system (γ = 0, top panel), ρ shows little to no evolution in both self-consistent (upward blue triangles) and non-interacting tracer (downward magenta triangles) runs. A small departure from the analytical profile around r ≈ 10⁻²r_c can be rather unambiguously interpreted as resolution effects as they decrease for higher resolutions and increase for lower ones (here N = 10⁶). Test particles, m_t, sinking in cuspy models appear to affect to a larger degree the density profile of the host system once reaching their inner regions, an effect that is related to the so-called BH scouring in elliptical galaxies hosting single or binary SMBHs (e.g. see Nasim et al. 2021 and references therein). In practice, strong scatters between stars on radial orbits and the massive BH deplete the stellar density profile, reducing the slope of a possibly initially steeper stellar cusp. In the simulations performed here (see bottom panel of Fig. 4, and the related inset), we observe a substantial depletion of the central cusp of the γ = 1 model in the simulation with non-interacting tracers, and an almost complete suppression of the cusp for the self-consistent case, where the background system becomes effectively ‘cored’ as in the N-body simulations reported by Goerdt et al. (2010).

Fig. 4. Initial density profile corresponding to a γ = 0 (top) and a γ = 1 (bottom) model (filled orange squares) and corresponding density profiles at t = 103tc in the self-consistent and tracers in smooth potential simulations (upward blue triangles and downward magenta triangles, respectively). The dashed black lines mark ρ(r) as given in Eq. (5). The points correspond to the same simulations of Fig. 3.

Prompted by the good agreement between the theoretical Holtsmark distribution for an infinite system with the force fluctuation distribution in spherical shells on a finite model (cfr. Fig. 2), it is natural to compare the orbital decay of a test mass in a N-body simulation with the trajectory of the same mass in the corresponding system modeled by Eq. (14). We performed several runs for different values of m_t/M and choices of the density profile, ρ, finding that, in general, the stochastic equation method acceptably recovers both the timescale over which the mean galactocentric distance of given test mass halves and, in the case of a cored model, its stalling radius, r_s. As an example, Figure 5 shows the evolution of r/r_c for a test particle, m_t = 3 × 10⁻³M, initially placed on a nearly circular orbit at r = 2r_c in a Plummer model with N = 10⁶, propagated in a N-body simulation (magenta curve) and in a stochastic simulation (cyan curve). For comparison, we also added the orbital decay evaluated adding the Chandrasekhar DF (1) to the smooth potential without fluctuations (dashed black line). In agreement with the literature discussed above, we observe that m_t stalls at around r_s ≈ 0.25r_c, with the N-body simulation always remaining in a sub-Chandrasekhar regime. Notably, at variance with the DF-only integration, where the particle falls to the centre of the Plummer system remaining on approximately circular orbits of decaying radius, in both the N-body and the stochastic simulations some finite change in eccentricity (more marked for the orbit obtained integrating Eq. (14)) can be appreciated.

Fig. 5. Decay of the galactocentric distance in a Plummer system for a particle with mt = 2 × 10−3M in a N-body simulation with N = 106 (magenta curve), and two semi-analytical integrations in the parent smooth potential with DF and fluctuations (solid cyan line) and with DF only (dashed black line).

3.2. N-body simulations

We investigated the resolution-induced partial suppression of DF by evolving the same initial condition for the test particle, m_t, in different realizations of γ models with an increasing number of particles, N, and fixed total mass, M. In the interval 10⁴ ≤ N ≤ 3 × 10⁶, for each explored value of m_t in the range between 2 × 10⁻⁴ and 10⁻² in units of M, the stalling radius, r_s, systematically settles at lower values for increasing N. We observe that, at fixed N and m_t and initial orbital parameters e and a, the stalling radius is always lower in models with a central cusp, in both self-consistent and simplified N-body simulations with tracers.

In Fig. 6 we show r_s/r_c as a function of N for cored and cuspy γ models for m_t = 10⁻³M initially placed on nearly circular orbits with the semi-major axis a ≈ 2r_c. We find that the stalling radius decreases with N, with a rather robust $1/\sqrt{N}$ trend, as is indicated by the dashed line. The square root tendency is indicative of a discreteness effect-related mechanism behind the core-stalling (e.g. see Gabrielli et al. 2014; El-Zant et al. 2019, and references therein). This becomes even more evident if the evolution of the galactocentric distance of the different realizations is rescaled by the estimated value of r_s times N^1/2. In Figure 7 we show the decay of r_s up to t = 10³t_c in a γ = 0 (top) and γ = 1 (bottom) model for the case of m_t = 10⁻³M in a system of non-interacting tracers and the associated rescaled curves for the last 50t_c (left panels). A good data collapse is obtained, as also the amplitude of the fluctuations of r_s fall on the same scale with a N^−1/2 trend (middle panels in Fig. 7). A similar picture was also recovered for the self-consistent simulations (not shown here).

Fig. 6. Core-stalling radius as function of N for a particle of mass mt = 10−3M in a cored model (γ = 0, circles) and a mildly cuspy model (γ = 1, triangles). Filled magenta symbols refer to the full self-consistent N-body simulations, and empty blue ones to the simulations with non-interacting tracers. To guide the eye, the dashed line marks the N−1/2 trend.

Fig. 7. Decay of the galactocentric distance of a tracer particle of mass mt = 10−3M in models with γ = 0 (top left panel) and 1 (bottom left) and the data collapse after a N1/2 scaling (mid panels) to the time-averaged scaling radius, rs, and the evolution of the corresponding mass enclosed with in a sphere of radius 5ϵ centred around mt. The horizontal dotted and dashed lines have the same meaning as in Fig. 3.

Aiming at comparing with the idealized Chandrasekhar DF, we have also integrated Eq. (14) for m_t in the parent smooth models (i.e. setting to 0 the fluctuating term, δf). The dashed black curves in the left panels of Fig. 7 mark the corresponding predicted orbital decay. Notably, in the γ = 0 cored case, the test particle might or might not experience the phase of super-Chandrasekhar DF reported by Read et al. (2006) and Zelnikov & Kuskov (2016), depending on the specific value of N in the simulation at hand. We observe that the onset of the super-Chandrasekhar regime happens at N ≲ 3 × 10⁵, suggesting that it might be connected with low-resolution effects artificially enhancing the collisionality. The same effect is also observable in self-consistent simulations using comparable values of ϵ. Interestingly, for all values of N explored here the test particle stalls at r_mt > r_s ≳ ϵ. We performed additional tests at fixed N and m_t tuning the softening length (and the associated optimal Δt) between 3 × 10⁻⁵r_c and 5 × 10⁻²r_c, finding that for increasing ϵ, the core stalling sets-in at later times and larger radii. According to the common wisdom, large values of the softening length render an N-body simulation ‘less collisional’ as it depletes the large δf tail of the two-body force distribution by smoothing the pair interactions. At the same time, a large ϵ also hinders the DF by reducing the effective Coulomb logarithm (Marcos et al. 2017; Ma et al. 2023).

Once m_t has stalled around r_s, the fraction of the host system’s mass that surrounds it (i.e. the mass fraction enclosed in a fixed sphere centred on the instantaneous position of m_t) remains almost constant, even in the simulations without self-gravity among the field particles. In the right panels of Fig. 7 we show as a function of time the mass fraction, M_d < 5ϵ/M, inside a bona fide radius, d = 5ϵ (roughly of the same order of r_inf for m_t = 10⁻³M). For the cored model (top right panel), modulo some resolution-dependent fluctuation, the value of M_d < 5ϵ is practically the same for all N, while in the case of a cuspy system (bottom right panel), M_d < 5ϵ at core stalling has the same value for N ≳ 10⁵. In other words, since at fixed m_t/M one has in our units approximately the same value of r_inf (cfr. Eq. (4)), this means that the test particle accretes the same amount of mass⁵ within its influence radius regardless of the specific resolution given by N. In practice, for increasing N a systematically larger number of orbits is trapped around m_t and essentially becomes a ‘dressed particle’, possibly contributing to a more significant erosion of the central density of cuspy systems via stronger field particle ejection (see Fig. 4).

As was expected, at fixed M and N, increasing values of m_t (marked by increasingly light shades of blue) result in lower stalling radii, as is shown in the top panel of Fig. 8 for γ = 0, N = 10⁶, and values of the test particle mass in the range of 2 × 10⁻⁴ ≤ m_t/M ≤ 10⁻². Consistently with the predictions of the Chandrasekhar DF formula (represented in the figure by the dashed lines), heavier particles also sink in a shorter time. The mass enclosed within 5ϵ around their position of m_t (bottom panel of Fig. 8) settles in general to slightly lower values for larger m_t (i.e. of the order of a 8% difference between the largest and lowest values of the test mass considered here). This might appear in contradiction with what is shown in Figure 7 (right panels) where at fixed m_t, M_d < 5ϵ is practically independent on N (and thus on m_t/m). However, one must notice that larger values of m_t at fixed m, associated with lower r_s, might scatter the field masses, m, more efficiently, depopulating the control sphere at d = 5ϵ.

Fig. 8. Decay of the galactocentric distance in a γ = 0 cored system for test particles with masses in the range of 2 × 10−4 ≤ mt/M ≤ 10−2 and N = 106 (top panel) indicated by increasingly light shades of blue. The dashed lines mark the radial decay predicted by the Chandrasekhar formula. Evolution of the mass enclosed with in a sphere of radius 5ϵ centred around mt (bottom panel).

In Figure 9 we show the relation between the stalling radius and the ratio m_t/m for two resolutions, N = 5 × 10⁵ (diamonds) and N = 10⁶ (squares). For both sets of points the trend is rather well described by $r_s=C{r}_c\sqrt{m/m_t}$, where C is a dimensionless constant of order unity. The stalling radius appears to have the same dependence on m/m_t as the wander radius⁶, r_wan, defined in Eq. (3). We interpret this as another indication of the origin of the core stalling as a balance between friction and noise in a central potential, as was originally envisaged for the derivation of r_wan with a Fokker-Planck approach (see again Bahcall & Wolf 1976). We note that other studies, such as Read et al. (2006) or Cole et al. (2012), report a general trend of r_s ∝ r_c/2, even in direct runs with a resolution as large as N = 4 × 10⁶. We conjecture that this apparent discrepancy with our numerical results could be related to their rather large softening length, ranging from ≈r_c/3 to 3r_c in some runs. We performed additional test runs at fixed N, m_t, and γ (not shown here), finding that the stalling radius, r_s, is systematically larger for increasing ϵ.

Fig. 9. Stalling radius in a γ = 0 model as a function of the mass ratio, mt/m, for the two system resolutions N = 5 × 105 (red diamonds) and N = 106 (blue squares). The dashed line marks the (mt/m)−1/2 trend.

3.3. Stochastic simulations

Aiming to explore larger values of N not accessible in direct N-body simulations, we have carried out additional numerical experiments evolving the same initial orbits of Fig. 6 using Equation (14) with resolutions up to N = 10¹³, roughly corresponding to the number of stars in a giant elliptical galaxy. Figure 10 presents r_s/r_c as function of N for the stochastic simulations with a γ = 0 (circles) and γ = 1 (triangles) parent galaxy model. Remarkably, in the interval of N where we have also performed fully self-consistent N-body and simpler simulations with non-interacting field particles, we recover the same N^−1/2 trend (dashed line) with similar values for r_s/r_c to that reported in Fig. 6. For the cuspy systems, such a trend extends up to N ≈ 10⁹, where for m_t = 10⁻³M the amplitude of the fluctuations, δf, stops balancing a significant fraction of the DF drag, then transitions to a shallower N^−1/6, indicated by the dotted brown line. The origin of this slope change remains unclear. We speculate that it could be a numerical artefact related to the choice of Δt in the solution of Eq. (14). In the cored model, the slope change happens at around N ≈ 10⁷, leading to a seemingly N-invariant r_s/r_c at fixed m_t for values of N larger than 10⁹. For larger values of m_t, r_s systematically decreases for fixed values of N, while the slope changes move at lower N. It must be pointed out that, for values of the test mass larger than ≈M/50, for which in large N models the effect of the force fluctuations would be essentially immaterial, the present linear approximation of the DF fails as m_t would typically induce significant perturbations of both the local density and velocity distribution, well before sinking into the core. Unfortunately, checking via N-body the nature of the slope change between N ≈ 10⁷ and 10⁹ is made harder by the fact that the associated value of r_s is around ϵ for N ≈ 10⁸, a radius that typically encloses for γ = 1 a rather small number (order unity) of simulation particles, possibly enhancing the amplitude of the discreteness effects while hardly being coherent with the perturbed smooth potential approximation modelled by Eq. (14).

Fig. 10. Core-stalling radius as function of N for a particle of mass mt = 10−3M in a cored model (γ = 0, cyan circles) and a mildly cuspy model (γ = 1, red triangles) as extrapolated from Langevin simulations. As in Fig. 6 the dashed lines mark the N1/2 best fit law in the range 104 ≤ N ≤ 3 × 106.

To illustrate the interplay between the contributions of the different force terms in Equation (14), in Figure 11 we plot for different N the radial dependence of the DF drag force acting on a mass, m_t = 10⁻³M, on a circular orbit with velocity

$$\begin{array}{c}\hfill v_c=\sqrt{\frac{GM{r}^{2-\gamma }}{{(r+r_c)}^{3-\gamma }}},\end{array}$$(24)

Fig. 11. Radial dependence of the DF force acting on a particle of mass mt = 10−3M on a circular orbit at vc (dashed lines), of the typical force fluctuation $\overline{\delta f}\approx {\alpha }^{2/3}$ (solid lines, increasingly light shades of red correspond to larger values of N), and radial dependence of the strength of the gravitational field, ||∇Φ|| (black heavy solid line). The left panel corresponds to the γ = 0 case and the right one to the γ = 1 case.

and the typical amplitude of the force fluctuation $\overline{\delta f}\approx {\alpha }^{2/3}$, together with the radial gravitational field ||∇Φ(r)||. In the γ = 0 cored model (left panel), the DF term drops significantly for r ≲ 0.1r_c as v_c tends to 0 while σ is nearly constant, resulting in a vanishing fractional velocity volume function (17). At the same time, for all values of N the typical size of the force fluctuations becomes almost independent on r as $\overline{\delta f}\propto {\alpha }^{2/3}$ and α ∝ ρ, with ρ nearly constant below r_c for a γ = 0 model. The picture is qualitatively different in the γ = 1 case (right panel), where the decrease in the DF force in the inner regions is way less pronounced and the amplitude of the fluctuating force increases as 1/r^2/3. To further clarify how the force fluctuations compete with DF, in Fig. 12 we plot for the same cases of Fig. 11 the difference in the amplitude of fluctuations and the DF force weighted over the local value of the gravitational field, ||∇Φ||. For increasingly large values of N the ratio (α^2/3 − ηv_c)/||∇Φ|| becomes positive (i.e. force fluctuations mostly balance DF) at smaller radii in units of r_c when γ = 1 (right panel). In the γ = 0 systems, for N ≳ 10⁹, the radius at which the change in the trend happens is almost N-independent (left panel), as is evident from Fig. 10, where r_s/r_c has little to no dependence on N in that resolution range. We noted that such intersections between the curves in Fig. 12 with the horizontal dashed line have the same trend with N as r_s/r_c in Fig. 10. Remarkably, in the γ = 1 case no apparent slope change at around N ≈ 10⁹ is observed, further indicating a possible numerical origin for that feature. We recall that (see e.g. Gardiner 1994 Chap. 4), the characteristic time, τ_cr, in which a system described by Eq. (14) relaxes is associated with the inverse of the friction coefficient, η (Spitzer 1969). The latter is essentially the first diffusion coefficient, D_f, in the parent Fokker-Planck picture (see Spitzer 1987). The other two coefficients, D_⊥ and D_∥, associated with perpendicular and parallel heating (i.e. noisy fluctuations of the force) along the trajectory of m_t, respectively, are directly proportional to m instead of to m_t as for D_f. This implies that the timescale on which they effectively operate is typically larger for the cases with smaller values of m_t/m, as the dependence on ρ and logΛ is the same as in D_f (see the integral expression in Genina et al. 2024). We stress the fact that in the simulation presented here (N-body and stochastic) the integration time is systematically larger or comparable to the inverse of the noise coefficients below r_c.

Fig. 12. Radial dependence of the difference between the DF term and typical force fluctuation, normalized by the local value of ||∇Φ||, for the same systems shown in Fig. 11. The colour-coding of the curves is the same as in the previous figure.

4. Outlook and conclusions

Using N-body simulations and stochastic ODE methods, we have investigated the DF acting on a heavy test mass in spherical systems characterized by a cored or cuspy density profile. Previous numerical studies of other authors have often reported a limited or in some cases substantial reduction of the gravitational drag in cored models, with respect to the estimate evaluated by the Chandrasekhar formula (1). Typically, such a discrepancy is interpreted as a consequence of the strong assumptions in the original formulation of DF in infinite and homogeneous systems. More recent semi-analytical work (Banik & van den Bosch 2021, 2022), on the other hand, indicates that DF in cores is naturally suppressed as its coefficient, dominated by the contribution of orbits in resonance with that of the test particle, vanishes towards the centre with the vanishing fraction of nearly resonant particles, in a resolution-independent fashion.

The main results of our work can be summarized as follows. The reduction of the DF force observed in N-body experiments appears to be resolution-dependent (i.e. related to the number of simulation particles when fixing the normalization G = M = r_c = 1) in the explored range, 10⁴ ≤ N ≤ 3 × 10⁶, in both cored and cuspy models. The radius at which the test mass, m_t, stalls without falling any further to the centre has a marked r_cN^−1/2 dependence indicative of an underlying Poissonian fluctuation-driven process. Single particle stochastic integrations with an analogous parameter reproduce the same trend, well beyond the N values accessible by the direct simulations (Fig. 10). We observe that approximating the distribution of force fluctuations in said experiments with a Holtsmark function, accommodated to the local values of density and mean particle mass, nicely recovers the behaviour of orbits obtained by N-body simulations. Moreover, the histogram of such force fluctuations, obtained from N-particle realizations of a given system, fits remarkably well the Holtsmark distribution over a large interval of radii. Comparing fully self-consistent runs with others in which the field particles do not interact with one another but propagate as independent tracers in the smooth potential of the model reveals that in the range of 3 × 10⁻⁴ ≲ m_t/M ≲ 10⁻² the contribution of a self-gravitating wake does not significantly enhance the DF. The stalling radius determined from N-body simulations has a trend that is reasonably compatible, $r_s\propto r_c/\sqrt{m_t/m}$. This is reminiscent of what was obtained for the wander radius of black holes in dense stellar systems. It must be pointed out that Merritt et al. (2007), in N-body experiments of BHs wandering in γ models, using a simulation set-up essentially analogous to ours, determined that said r_wan in cusps broken below the influence radius, r_inf, could scale with m_t with a much shallower power law. This arises as the BHs’ typical velocity at equilibrium was found to be ⟨v²⟩∝m_t^{−1/(γ − 1)}, while the rms of its radial position is given by $⟨r^2⟩\propto \sqrt{⟨v^2⟩}$ when a linear restoring force is present (i.e. inside a nearly flatted out cusp the potential could reasonably be taken as harmonic). Nevertheless, we interpret the similarity of the trend in r_s and (the expected one) in r_wan as an additional confirmation of the fact that core-stalling induced by a balance between friction and force fluctuations is the same phenomenon at the origin of a mass-dependent wander radius for a heavy particle starting at rest from the centre of the potential well of a star cluster or galactic nucleus.

From the stochastic simulations, we extrapolate that in cuspy models the interplay between friction and fluctuation-related noise balances at smaller radii for increasing N (up to a regime in which δf becomes virtually inexistent), while for cored models there is a critical size (N_* ≈ 3 × 10⁷ for m_t ≈ 10⁻³M) above which friction and noise counter each other seemingly independently of the resolution. We must point out that in the continuum limit (i.e. N → ∞; m → 0), in principle, both friction and noise should vanish, leaving m_t under the effect of Φ only. At variance with the analogous problems in plasma physics, where typically one deals with N of the order of the Avogadro constant, a stellar system with N ≈ 10¹² is still far from being considered a continuum, and is still subjected to a certain degree to a ‘discreteness effect’.

In conclusion, the results of the present work are not in contrast with the resonance-based and N-independent argument of Banik & van den Bosch (2022) that ascribes the core stalling to the vanishing friction coefficient in cores. Indeed, our findings reconcile the predicted N-independent cores stalling with the observed resolution dependent stalling radius, r_s, reported by other N-body studies. In other words, for system sizes of the order of N ∼ 10³ up to 10⁶, the radial interval over which both the friction and noise terms in Eq. (14) are non-negligible with respect to the local smooth force field is significantly larger, regardless of the fact that the model is cuspy or cored. In cored models, for a sufficiently small test mass, m_t, the further outward extent of this interval is larger than the critical radius at which the DF force on a nearly circular orbit begins to drop significantly. For increasing resolution (i.e. larger values of N at fixed normalization and m_t/M), in cored models the r-independent fluctuations are smaller and counter DF nearly the same radius below which the friction coefficient has dropped off about two decades over less than a decade in radii (cfr. Figs. 11 and 12). In the end, observing the onset of a resonance-driven core stalling in N-body experiments remains a challenging task, as serial direct numerical simulations are strongly limited in resolution. However, considering a test mass moving in a host system where field particles do not feel their mutual self gravity but orbit in a time-independent field is a problem that can be easily partitioned over different cores or graphic units, allowing for a N up to 10¹⁰. In this set-up, choosing a distribution function entirely constituted by circular orbits (that would likely be unstable when accounting for the self consistent forces) could serve as a computational test to probe a possible resonant-driven suppression of the DF.

Acknowledgments

We thank Uddipan Banik, Frank van den Bosch, and Carlo Nipoti for the important discussions at different stages of this project. The anonymous Referee is also thanked for his/her valuable comments that helped improving the presentation of our results. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). PFDC wishes to thank the hospitality of the “Laboratoire J.-A. Dieudonné” at the university of Côte d’Azur where part of this work was completed. BM acknowledges support by the grant Segal ANR-19-CE31-0017 of the French Agence Nationale de la Recherche.

References

Appendix A: Analytical dynamical friction coefficient for the Plummer sphere

The DF expressed in Eq. (1) is often evaluated assuming that the phase space distribution of the underlying model can be factorized as f(r, v) = ρ(r)f(v). In the commonly adopted local Maxwellian approximation, the position dependence of the velocity distribution f(v) is restored by assuming a Maxwellian distribution with velocity dispersion proportional to r via the Jeans equations for the associated density potential pair (ρ, Φ). In principle, in a spherically symmetric model with isotropic phase space distribution, the integral in Eq. (2) could be instead carried out over f(ℰ) = f(r, v). However, most often f has a rather complicated form even for simple density profiles such as the γ-models. For the Plummer distribution, where f(ℰ) is given by Eq. (12) the integration becomes rather easy (e.g. see Miocchi et al. 2015) and one has

$$\begin{array}{c}\hfill {\mathrm{\Xi }}^{\ast }(r,v)=4\pi 2^{3/2}k\mathrm{\Psi }{\displaystyle {\int }_0^{v{(2\mathrm{\Psi })}^{-1/2}}}{w}^2{(1-w^2)}^{7/2}\mathrm{d}w=\\ \hfill =4\pi 2^{3/2}k\mathrm{\Psi }{\displaystyle {\int }_u^{\infty }}\frac{y^8}{{(y^2+1)}^6}\mathrm{d}w.\end{array}$$(A.1)

In the equation above $k=24\sqrt{2}{r}_c^2/7{\pi }^3{G}^5{M}^4$, $w=v/\sqrt{2\mathrm{\Psi }}$ and $y=\sqrt{w^{-2}-1}$. Setting $u=\sqrt{2\mathrm{\Psi }/v^2-1}$ one has

$$\begin{array}{c}\hfill {\mathrm{\Xi }}^{\ast }(r,v)=\alpha {\mathrm{\Psi }}^5[\frac{1}{2}{tan}^{-1}(u)+(\frac{u^9}{2}-\frac{79}{21}{u}^7-\frac{64}{15}{u}^5-\frac{7}{3}{u}^3-\frac{u}{2})\\ \hfill \times {(1+u^2)}^{-5}-\frac{\pi }{4}],\end{array}$$(A.2)

where $\alpha =-7\pi \sqrt{2}k/16$.

In Fig. A.1 we show the radial profile of the DF force acting on a test mass m_t = 10⁻³M moving on a circular orbit at the local circular velocity $v_c=\sqrt{GMr{(r^2+r_c^2)}^{-3/2}}$, evaluated in the local Maxwellian approximation (dashed lines) and integrating over the distribution function (solid line), so that the term ρΞ(v) in Equation (1) is substituted by Ξ^*(r, v) defined above. At r ≫ r_c the two expressions do not give appreciably different drag force, while below r ≈ 5 × 10⁻¹r_c, the local Maxwellian approximation overestimates the DF of a factor about 1.4, we conclude that using such approximate form for the friction force does not alter significantly its effect, in particular since the region where the error is greater corresponds to a radial range where its contribution itself becomes under-dominant.

Fig. A.1. Radial profile of the DF force on a mass mt = 10−3M on a circular orbit in a Plummer model.

Appendix B: Further comparison between self-consistent and tracer simulations

As the asymptotic stalling radius is roughly the same in self-consistent simulation and simplified smooth potential and tracer ones, we expect that the (relative) energy of m_t relaxes to comparable values once the its radial coordinate has settled at ≈r_s. In Fig. B.1 we show the evolution of ℰ for the same runs of Fig. 3. In this case, where the resolution is N = 10⁶ we observe a remarkably good agreement of the final values of the relative energy for both values of γ with relatively small fluctuations. Smaller values of N (not shown here) also yield rather similar values, though bearing larger fluctuations around the final value for N ≲ 10⁵. By contrast, the norm of the angular momentum J shown for the same simulations in Fig. B.2 shows wilder fluctuations (even at N as large as 3 × 10⁶), of about the 10% of its asymptotic value. We interpret this as a result of the more chaotic behaviour of orbits in the inner region of both models that reflects in violent changes in orbital eccentricity and inclination while the energy remains nearly constant.

Fig. B.1. Evolution of the relative orbital energy per unit mass ℰ for a test mass mt = 10−3 in a γ = 0 (left) and γ = 1 (right) model. In both cases N = 106 and the blue and magenta curves refer to the tracer in smooth potential and self-consistent runs, respectively.

Fig. B.2. Evolution of the norm of the orbital angular momentum J for the same systems of Fig. B.1.

All Figures

	Fig. 1. Radial profile of the average of the discrete force fluctuations (left panel) and its variance (right panel) evaluated for a single realization of a Plummer system with N = 106.
In the text

Fig. 2. Numerically evaluated distribution of force fluctuation, 𝒫(δf), in different shells of increasing radius (clockwise from top left, r = 3.1 × 10−2, 7.2 × 10−2, 0.38, 1.3, 6.9, and 54.8 in units of rc) for a Plummer model with N = 106 and 102 realizations (grey dots). The blue dots mark the PDF averaged over said ensemble, while the red curve is the corresponding the Holtsmark distribution (15) computed for the local value of the density, ρ(r), and average number of particles in the shell, $\overline{N}$.

In the text

Fig. 3. Radial decay of a massive tracer, mt = 10−3M, orbiting in a γ = 0 (left) and a γ = 1 (right) model in a self consistent N-body run (blue curves) and a simplified simulation where the same N = 106 field particles move independently in the static smooth potential (magenta curves). The horizontal dashed orange line marks the value of the softening parameter ϵ = 2 × 10−3rc, while the dotted black line indicates the radius, rmt, at which M(r) = mt, equal to ≈0.1 for the simulations with γ = 0 and ≈0.033 for those with γ = 1.

In the text

Fig. 4. Initial density profile corresponding to a γ = 0 (top) and a γ = 1 (bottom) model (filled orange squares) and corresponding density profiles at t = 103tc in the self-consistent and tracers in smooth potential simulations (upward blue triangles and downward magenta triangles, respectively). The dashed black lines mark ρ(r) as given in Eq. (5). The points correspond to the same simulations of Fig. 3.

In the text

	Fig. 5. Decay of the galactocentric distance in a Plummer system for a particle with mt = 2 × 10−3M in a N-body simulation with N = 106 (magenta curve), and two semi-analytical integrations in the parent smooth potential with DF and fluctuations (solid cyan line) and with DF only (dashed black line).
In the text

	Fig. 6. Core-stalling radius as function of N for a particle of mass mt = 10−3M in a cored model (γ = 0, circles) and a mildly cuspy model (γ = 1, triangles). Filled magenta symbols refer to the full self-consistent N-body simulations, and empty blue ones to the simulations with non-interacting tracers. To guide the eye, the dashed line marks the N−1/2 trend.
In the text

Fig. 7. Decay of the galactocentric distance of a tracer particle of mass mt = 10−3M in models with γ = 0 (top left panel) and 1 (bottom left) and the data collapse after a N1/2 scaling (mid panels) to the time-averaged scaling radius, rs, and the evolution of the corresponding mass enclosed with in a sphere of radius 5ϵ centred around mt. The horizontal dotted and dashed lines have the same meaning as in Fig. 3.

In the text

	Fig. 8. Decay of the galactocentric distance in a γ = 0 cored system for test particles with masses in the range of 2 × 10−4 ≤ mt/M ≤ 10−2 and N = 106 (top panel) indicated by increasingly light shades of blue. The dashed lines mark the radial decay predicted by the Chandrasekhar formula. Evolution of the mass enclosed with in a sphere of radius 5ϵ centred around mt (bottom panel).
In the text

	Fig. 9. Stalling radius in a γ = 0 model as a function of the mass ratio, mt/m, for the two system resolutions N = 5 × 105 (red diamonds) and N = 106 (blue squares). The dashed line marks the (mt/m)−1/2 trend.
In the text

	Fig. 10. Core-stalling radius as function of N for a particle of mass mt = 10−3M in a cored model (γ = 0, cyan circles) and a mildly cuspy model (γ = 1, red triangles) as extrapolated from Langevin simulations. As in Fig. 6 the dashed lines mark the N1/2 best fit law in the range 104 ≤ N ≤ 3 × 106.
In the text

Fig. 11. Radial dependence of the DF force acting on a particle of mass mt = 10−3M on a circular orbit at vc (dashed lines), of the typical force fluctuation $\overline{\delta f}\approx {\alpha }^{2/3}$ (solid lines, increasingly light shades of red correspond to larger values of N), and radial dependence of the strength of the gravitational field, ||∇Φ|| (black heavy solid line). The left panel corresponds to the γ = 0 case and the right one to the γ = 1 case.

In the text

	Fig. 12. Radial dependence of the difference between the DF term and typical force fluctuation, normalized by the local value of ||∇Φ||, for the same systems shown in Fig. 11. The colour-coding of the curves is the same as in the previous figure.
In the text

	Fig. A.1. Radial profile of the DF force on a mass mt = 10−3M on a circular orbit in a Plummer model.
In the text

	Fig. B.1. Evolution of the relative orbital energy per unit mass ℰ for a test mass mt = 10−3 in a γ = 0 (left) and γ = 1 (right) model. In both cases N = 106 and the blue and magenta curves refer to the tracer in smooth potential and self-consistent runs, respectively.
In the text

	Fig. B.2. Evolution of the norm of the orbital angular momentum J for the same systems of Fig. B.1.
In the text