© The Authors 2023

1. Introduction

A perturber moving through a host – an environment of collisionless particles – is subject to dynamical friction. The perturber distorts the host’s density distribution, and this distortion alters the velocity of the perturber. Dynamical friction was first described by Chandrasekhar (1943) for an infinite, homogeneous, and isotropic environment and later advanced by Tremaine & Weinberg (1984) for a spherically symmetrical host. Since then, several other generalisations have been developed, for example, multi-component systems and ellipsoidal velocity distribution (Bonetti et al. 2021; Moreno et al. 2022). Heyvaerts et al. (2017) evaluated dynamical friction via diffusion coefficients and showed equivalence of many representations. Banik & van den Bosch (2022) reformulated the dynamical friction modelling in an orbital-based framework with the aim of describing the central regions of galaxies. The response of a galaxy to a massive perturber has been studied with simulations (see e.g. Vasiliev et al. 2021; Vasiliev 2023) as well as semi-analytically (Bonetti et al. 2021). Most dynamical friction studies emphasise the loss of integrals of motion of the perturber or its influence on the host. The aim of this paper drifts from these goals, as we aim to infer the host properties based on the distorting forces induced by the wake of the dynamical friction towards the perturber.

Dynamical friction depends on the velocity distribution of a host’s particles, and one can use this dependency to infer the host’s dark matter (DM) kinematics. Since the counterpart particle of DM has not yet been found, an enormous effort is being made to determine its nature from astronomical observations. Dark matter inference studies on large and galactic scales have limitations due to the precision of observations (Fischer et al. 2022) and/or the modelling methods (Roper et al. 2023; Downing & Oman 2023). On smaller scales, a particularly powerful discriminator of the DM nature is the mass function of the DM substructures (Benito et al. 2020; Lovell et al. 2021; Nadler et al. 2021). Dynamical friction offers an alternative way to study the small-scale effects induced by hypothetical microphysical properties of DM (Hui et al. 2017; Hartman et al. 2021). Moreover, its applications are much broader than DM inference, as dynamical friction is the underlying phenomenon that influences a wide variety of astrophysical systems, from stellar streams and globular clusters to the orbits of satellite galaxies and galaxy merging times (Kipper et al. 2019, 2020; Pflamm-Altenburg et al. 2023; Correa Magnus & Vasiliev 2022).

In essence, dynamical friction is a process that makes a heavier object slow down in comparison to its environment. As the slowing down effect depends on the host’s environment, dynamical friction can also be used as a tool to study both host and perturber properties. To put this tool into action, we distinguished two observational regimes.

In the first regime, the inference can be based on a long-term influence on the perturber orbit. Weak and persistent dynamical friction influences an object over an extended period of time, thus altering the radial location of the object in a galaxy or cluster. A prominent example of a long-term effect is core stalling. The globular clusters in the Fornax dwarf are an example of this, as they have not yet fallen to the centre (see Meadows et al. 2019; Bose et al. 2019). Another example of long-term dynamical friction manifestation is its influence on the relative position of a globular cluster compared to a stream centre (Kipper et al. 2019).

The inference can be based on the tidal field by the wake. Dynamical friction decreases the velocity of the perturber and additionally causes tidal forces to affect the perturber. Some parts of the perturber slow down faster than others, causing distortions in shape. This was first studied by Mulder (1983) and in the context of inferring DM properties by Garavito-Camargo et al. (2019) and Kipper et al. (2020). Kipper et al. (2020) studied an isolated galaxy moving in a void environment and showed that dynamical friction due to intergalactic DM can produce lopsided gas and stellar distributions. In this paper, our attention is focused on the tidal aspect of the dynamical friction. We develop a model with the emphasis to use as few simplifying assumptions as possible when using the method.

Simplifying assumptions and their validity are central for most dynamical methods. Solving the equations of motion that include dynamical friction requires using simplifying assumptions that might influence the accuracy of the results. For example, the classical methods of Chandrasekhar and Tremaine & Weinberg (1984) assume that the underlying medium is uniform, spherical, or axisymmetric. In many galaxies, this is not the case. For example, the Milky Way disc has perturbations and non-equilibrium features (Antoja et al. 2018). Another common assumption is to use nearly circular orbits for the perturber. However, most globular clusters do not have circular orbits (Vasiliev & Baumgardt 2021). Lastly, several studies assume that the density distribution of a galaxy is time independent and thus discard secular evolution of the host (Chandrasekhar 1943; Tremaine & Weinberg 1984; Banik & van den Bosch 2021; Desjacques et al. 2022), but some evolution has been shown to occur (Gunawardhana et al. 2013; Kipper et al. 2021a; Hashemizadeh et al. 2022). These kinds of assumptions lead to biases in dynamical friction calculations. Nonetheless, relaxing these assumptions poses problems for analytical calculations. Analytical description of the model is close to mandatory in case one wishes to infer something about specific objects using likelihood, although there are a few exceptions, such as Vasiliev et al. (2021), Reza et al. (2022). Moreover, numerical simulations can be used for calculating dynamical friction, but they are time-consuming, and since slight deviations from initial conditions may give quite different results (Genel et al. 2019), they can only provide a statistical understanding (Amorisco 2017).

The tidal aspect of dynamical friction is not well studied. The pioneering paper by Mulder (1983) solved it in the case of an isotropic velocity distribution and concluded that tidal effects are mainly from regions near the perturber. Another paper discussing tidal aspects of dynamical friction is Kipper et al. (2020), where we showed that isolated galaxies can form lopsidedness without mergers just by moving through large-scale structures.

In this paper, we propose a semi-analytic method to calculate dynamical friction that relaxes many of the commonly adopted assumptions. This enables a more accurate estimation of the tides of dynamical friction and the sensitivity needed to investigate the nature of DM. The method can be used with real Gaia data to study the Milky Way. Testing of the method indicates that the process of the dynamical friction is subject to fluctuations of forces from background stars. This aspect has not obtained sufficient attention yet.

The paper is organised as follows. Section 2 presents our method for calculating dynamical friction and its extension for tidal field evaluations. Section 3 summarises the simulated tools that are used in Sect. 4 to validate the method and evaluate its limitations. Then, Sect. 5 assesses the dependence of the dynamical friction estimates on sampling noise. We conclude with a discussion and a summary in Sects. 6 and 7.

2. Method

From this point and onwards, we refer to the tides from the wake of dynamical friction as ‘tidal dynamical friction’. Preserving the complexity of phase-space density with maximal accuracy is important for realistic tidal dynamical friction estimates. For example, Petts et al. (2016) and Leaman & van de Ven (2022) showed that deviations from the simple Maxwellian velocity distribution significantly affect dynamical friction estimates. Hence, we must be able to accurately account for all the complexities. We describe the host as a set of particles, each with its own position and velocity, which allows for an accurate characterisation of the phase-space density. However, this particle-based description is challenging to acquire observationally, which may necessitate the use of smooth distribution functions. In this work, we develop a framework for tidal dynamical friction calculations suitable with both phase-space and particle-based descriptions of the system.

Our aim is to obtain the tidal dynamical friction estimation, which is mostly determined by regions in the vicinity of the perturber (Mulder 1983). We kept this feature in mind when deriving the methodology.

2.1. Basic framework

The setup for estimating dynamical friction consists of an environment – the host – and a massive perturber moving through the host. The acceleration field caused by dynamical friction is the difference between the acceleration fields generated by the host’s density distribution when the perturber is present or absent. Our definition of dynamical friction is slightly wider than what is commonly used, as it includes evaluation points other than the perturber itself. We denote the position of a particle in the absence of the perturber as x₀, and in the presence of the perturber, we denote it as x₁. The total number of particles in the system is N_tot. We describe the accelerations at an arbitrary point x in the host as a₀ in the absence of the perturber and as a₁when the perturber exists:

$$\begin{array}{cc}& {\mathbf{a}}_0(\mathbf{x},t)={\displaystyle \sum _{i=1,i\ne j}^{N_{\mathrm{tot}}}}{K}_0\frac{{\mathbf{x}}_{0,i}-\mathbf{x}}{|{\mathbf{x}}_{0,i}-\mathbf{x}|}\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\mathbf{a}}_1(\mathbf{x},t)={\displaystyle \sum _{i=1,i\ne j}^{N_{\mathrm{tot}}}}{K}_1\frac{{\mathbf{x}}_{1,i}-\mathbf{x}}{|{\mathbf{x}}_{1,i}-\mathbf{x}|},\hfill \end{array}$$(2)

respectively. The condition i ≠ j is only used when the acceleration is evaluated for particle j. The term K denotes a kernel describing a particle’s acceleration, which can take different forms depending on the system we are interested in. For example, in simulation comparison, the kernel should include softening. For point sources with mass m, the kernel is Newton’s gravitational acceleration:

$$\begin{array}{c}\hfill K=\frac{\mathit{Gm}}{r^2},\end{array}$$(3)

where r in K₀ corresponds to |x_0, i − x| for the unperturbed state, while r in K₁ corresponds to |x_1, i − x| for the perturbed state.

These equations hold at any fixed time. Dynamical friction manifests as the difference between these acceleration fields:

$$\begin{array}{c}\hfill {\mathbf{a}}_{\mathrm{DF}}\equiv {\mathbf{a}}_1-{\mathbf{a}}_0.\end{array}$$(4)

We call the a_DF total ‘one’ when summing in Eqs. (1) and (2) is done over all the points and refer to the total as ‘local’ when only points in the vicinity of the perturber are included. In the first case, a₀ and a₁ correspond to the acceleration field of the host $-\nabla {\Phi }_{\text{h}}^0(\mathbf{x},t)={\mathbf{a}}_0$ or $-\nabla {\Phi }_{\text{h}}^1(\mathbf{x},t)={\mathbf{a}}_1$ in the perturbed state. The difficulty regarding dynamical friction calculations is that, generally, only one of these acceleration fields is known, either Eqs. (1) or (2). In order to get one from the other, we devised a method that is based on initially integrating the particles and perturber back in time until the perturber has a negligible influence on the particles. Then, we remove or add the perturber’s potential to the host and re-integrate the orbits back to the present. A similar method was developed by Correa Magnus & Vasiliev (2022) to evaluate the mass of the Milky Way. Although from its theoretical aspects, this is a simple procedure, from the practical aspects, such as selection function, uncertainties, and the accuracy of gravitational potential, this method poses several questions during its implementation.

We define the time t = 0 as the moment at which we wish to calculate the tidal dynamical friction and assume that we know the unperturbed distribution of particles at that time, that is, their positions x₀(t = 0) and velocities v₀(t = 0). The equations of motion describe the movement of particles, where the host potential fully determines the acceleration of these particles:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{v}}_0}{\mathrm{d}t}=-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{h}}^0({\mathbf{x}}_0,t).\end{array}$$(5)

The host potential depends on position and time. The explicit time dependence can include contributions from a rotating bar or other external processes (i.e. non-stationarity). However, in the presence of the perturber, described with potential Φ_p and coordinate x_p, the particle accelerates according to

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{v}}_1}{\mathrm{d}t}=-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{h}}^1({\mathbf{x}}_1,t)-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{p}}({\mathbf{x}}_1|{\mathbf{x}}_{\mathrm{p}}(t),t).\end{array}$$(6)

Unfortunately, regarding the use of Eq. (6) to calculate the orbit of a disturbed particle, we do not yet know the potential Φ¹. Thus, we calculate it iteratively and replace Φ¹ with a guessed potential Φ^1*, namely,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{v}}_1}{\mathrm{d}t}=-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{h}}^{1\ast }({\mathbf{x}}_1,t)-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{p}}({\mathbf{x}}_1|{\mathbf{x}}_{\mathrm{p}}(t),t).\end{array}$$(7)

For the first iteration we take

$$\begin{array}{c}\hfill {\mathrm{\Phi }}^{1\ast }={\mathrm{\Phi }}^0.\end{array}$$(8)

Taking this approach means we approximate the potentials without taking into account internal deformations of the galaxies, that is, for the first iteration, the self-gravity of the wake is not accounted for. This regime is satisfied when two galaxies gravitate towards each other, and the smaller has undergone only one pericentre passage (Vasiliev et al. 2022; Correa Magnus & Vasiliev 2022). For subsequent iterations, we use the following approximation to Φ^1*

$$\begin{array}{c}\hfill \mathrm{\nabla }{\mathrm{\Phi }}^{1\ast }=\mathrm{\nabla }{\mathrm{\Phi }}^0+{\mathbf{a}}_{\mathrm{DF}\ast },\end{array}$$(9)

where a_DF* is evaluated¹ using Eq. (2), but the coordinates we use are ${\mathbf{x}}_1^{*}$ instead of x₁. Upon several iterations, Φ^1* approaches Φ¹. The perturber’s trajectory x_p(t) should also be calculated according to

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbf{v}}_{\mathrm{p}}}{\mathrm{d}t}=-\mathrm{\nabla }{\mathrm{\Phi }}^{1\ast }.\end{array}$$(10)

We assume that the initial conditions of the perturber (coordinates x_p and velocity v_p) are known.

Solving Eqs. (5) or (6) numerically requires knowing the initial conditions or the conditions at time t = 0. In the case where the perturber is observed, the initial conditions are x_1, i(t = 0) = x_obs, i. If we want to calculate the dynamical friction of a hypothetical perturber located at x_p, then the initial conditions are x_0, i(t = 0) = x_obs, i. From now on, we assume the latter case, but we can derive the equations for the former case in a similar manner.

In order to convert the position of a particle from the unperturbed case (x₀) to the perturbed one (x₁), we make the assumption that there is a time t_far ≤ 0 where the particle on its orbit is so far from the perturber that the perturber’s contribution to the particle acceleration is negligible. Mathematically, this is equivalent to the following condition:

$$\begin{array}{c}\hfill |\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{h}}({\mathbf{x}}_0(t_{\mathrm{far}}))|\gg |\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{p}}({\mathbf{x}}_0(t_{\mathrm{far}})|.\end{array}$$(11)

The practical limitations of this condition to exist are explained in Sect. 2.3. At this point, we may integrate the orbit for a particle i at the acceleration field a₀(x, t) with the initial condition x_0, i(t = 0) = x_obs, i until t_far and derive x_0, i(t_far). At t = t_far the particle has no information on whether the perturber exists or not. Therefore, at t = t_far,

$$\begin{array}{c}\hfill {\mathbf{x}}_1(t_{\mathrm{far}})={\mathbf{x}}_0(t_{\mathrm{far}}),\end{array}$$(12)

and we may take x_1 * ,i(t_far) as the initial condition to integrate the orbit x_1 * ,i until t = 0. Again, the asterisk denotes the potential approximation Φ^1* in use. With this procedure, we are able to calculate x_1*(t = 0) from x₀(t = 0) for all the particles. We then enter the newly acquired positions x_1* into Eq. (2), which are the basis to calculate (tidal) dynamical friction with Eq. (4). This method is illustrated in Fig. 1.

Fig. 1. Illustration of the dynamical friction framework. The magenta points describe the points of the host in the vicinity of a perturber (depicted as a blue dot in the enlarged sub-plot). Different coloured lines show a small sub-sample of integrated orbits, both backwards and forwards in time. The dashed circle shows the distance from the perturber where the tfar is reached. The acceleration field difference caused by unperturbed (x0) and perturbed (x1) points gives the dynamical friction. In the enlarged image, we connect the corresponding points with a grey line.

The equations above correspond to the dynamical friction for a particle-based description. In the case of a continuous description of the system, the equations change slightly:

$$\begin{array}{cc}& {\mathbf{a}}_0(\mathbf{x})={\displaystyle \int }{K}_0\frac{{\mathbf{x}}_0-\mathbf{x}}{|{\mathbf{x}}_0-\mathbf{x}|}f({\mathbf{x}}_0,{\mathbf{v}}_0)\mathrm{d}{\mathbf{x}}_0\mathrm{d}{\mathbf{v}}_0,\hfill \end{array}$$(13)

$$\begin{array}{cc}& {\mathbf{a}}_1(\mathbf{x})={\displaystyle \int }{K}_1\frac{{\mathbf{x}}_1-\mathbf{x}}{|{\mathbf{x}}_1-\mathbf{x}|}f({\mathbf{x}}_1,{\mathbf{v}}_1)\mathrm{d}{\mathbf{x}}_1\mathrm{d}{\mathbf{v}}_1.\hfill \end{array}$$(14)

We denote f as the phase-space distribution normalised so that integration over velocities is the ordinary spatial density: ∫fdv = ρ.

2.2. The tidal aspect of dynamical friction

The last section describes the dynamical friction in a broad context. In this section, we shift the attention to tidal aspects that originate locally. As local regions are near the perturber and they are able to provide high-leverage contributions, noise, and fluctuations, we can estimate how much the sampling noise propagates to dynamical friction.

To gain insights into the forming of dynamical friction and its fluctuations, we return to Eqs. (1), (2), and (4). Rewriting Eq. (4) in a simplified, one-dimensional form for a single particle contribution, we obtain

$$\begin{array}{c}\hfill \delta a=\frac{\mathit{Gm}}{r^2}-\frac{\mathit{Gm}}{{(r-\delta x)}^2}=Gm\left[\frac{\delta x^2-2r\delta x}{r^2(r^2-2r\delta x+\delta x^2)}\right].\end{array}$$(15)

We use r to denote the distance of the unperturbed position to the evaluated point, and δx = x₁ − x₀ is the shift between perturbed and unperturbed positions (i.e. the source of dynamical friction). In case the force-inducing pair is far from the evaluated point (δx ≪ r), the above relation becomes

$$\begin{array}{c}\hfill \delta a=\frac{-2Gm\delta x}{r^3}.\end{array}$$(16)

That is, the dynamical friction asymptotically has a tidal-like (r⁻³) behaviour. If the force-inducing pair is acting as the source of the tidal force (i.e. different distances (δy) are affected differently by this pair), then the contribution analogously becomes

$$\begin{array}{c}\hfill \delta a_{\mathrm{tide}}=\frac{-2Gm\delta x}{r^3}-\frac{-2Gm\delta x}{{(r-\delta y)}^3}\propto r^{-4}.\end{array}$$(17)

We use δa_tide in order to denote the contribution of the pair to the tidal dynamical friction.

The dynamical friction does not result from a single pair of points but from all points in the host. We divided the points among shells with a mean distance to the perturber denoted as r. Combining the shell contributions into the total signal is done repeatedly (for Newton acceleration, dynamical friction, and tidal dynamical friction); hence, we denote the total signal as F and its variance as σ², the kernel (gravitational force ∝r⁻², dynamical friction ∝r⁻³, or tidal force from dynamical friction ∝r⁻⁴) as w, and the contribution from each shell i as W_i.

Each shell’s contribution W_i is proportional to the number of points in the shell and the kernel ∝N_ptw. We considered a shell with thickness Δr. The number of particles in the shell is proportional to the density (n) and the volume of the shell: W ∝ nr²Δrw. The corresponding uncertainty for each shell originates from sampling, as the Poisson distribution gives contributions only from the density part, while the kernel is constant $\mathrm{\Delta }{W}_i\propto \sqrt{N_{\mathrm{pt}}}w\propto \sqrt{r^2\mathrm{\Delta }r}w=rw\sqrt{\mathrm{\Delta }r}$ and variance $\mathrm{\Delta }{W}_i^2\propto w^2{r}^2\mathrm{\Delta }r$. After finding each shell contribution, they need to be combined. Summing or integrating the contributions over all shells at the thin shell limit (Δr → dr) to either signal (W) or its variance (ΔW) gives F ∝ ∑W_i = ∫r²wdr or σ² ∝ ∑(ΔW_i)² = ∫r²w²dr. In order to find the signal-to-noise ratio, we divided the two $\sqrt{{\sigma }^2}/F$.

Next, we compared the convergence of signal and noise in cases of Newton gravity, dynamical friction, and the tidal dynamical friction in the context of signal, noise, and signal-to-noise ratio (see Table 1). As can be seen in the case of dynamical friction, the forces accumulate proportionally with lnr, precisely as in case of Chandrasekhar dynamical friction. But a new result is the formation of the signal-to-noise ratio. In the case of dynamical friction, the ratio accumulates faster, indicating that the inner parts contribute more to noise accumulation than in the case of Newton’s law of gravity. For the tidal dynamical friction, one can see that the total contribution behaves as r⁻¹, causing inner parts to be far more relevant in the signal accumulation than the outer parts. This was reported in Mulder (1983), and we reproduce it in this work. Analogously, one can also see that the noise originates from the vicinity of the perturber.

The previous line of reasoning relies on the assumption of a homogeneous environment, but this is not generally the case, and this analysis only provides an understanding of what to expect. For the general uncertainty estimations from sampling noise, we used bootstrapping (Efron 1979).

The tidal dynamical friction provides the force distorting the perturber. How exactly the perturber is distorted is shown by the shapes of the orbits bound to the perturber. Denoting δx_tidal as the distance from the perturber where we want to evaluate a particle orbit, the tidal dynamical friction term

$$\begin{array}{c}\hfill \delta {\mathbf{a}}_{\mathrm{DF}}=\frac{\mathrm{d}{\mathbf{a}}_{\mathrm{DF}}}{\mathrm{d}\mathbf{x}}{|}_{\mathbf{x}={\mathbf{x}}_{\mathrm{p}}}\delta {\mathbf{x}}_{\mathrm{tidal}}\end{array}$$(18)

must be included in the acceleration field. We take the derivative at the location of the perturber. The tidal field from the dynamical friction affects the orbits in or around the perturber, and they are evaluated based on acceleration

$$\begin{array}{c}\hfill \dot{\mathbf{v}}-{\dot{\mathbf{v}}}_{\mathrm{p}}=-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{p}}+\delta {\mathbf{a}}_{\mathrm{DF}}-\mathrm{\nabla }({\mathrm{\Phi }}_{\mathrm{h}}^0-{\mathrm{\Phi }}_{\mathrm{h}}^0{|}_{{\mathbf{x}}_{\mathrm{p}}}).\end{array}$$(19)

The first term describes the perturber’s potential, the second describes the tidal dynamical friction, and the third is the host’s tidal field.

To summarise, the current method is aimed at inferring the host properties using the shape distortions the wake induces to the perturber. The shape distortions originate from the dynamical friction differences at points in the vicinity of the perturber. The mean dynamical friction averaged over all points near the perturber does not contribute to the distortions but causes the slowing down of the perturber. In this way, we can split the dynamical friction into two parts:

$$\begin{array}{c}\hfill {\mathbf{a}}_{\mathrm{DF}}={\overline{\mathbf{a}}}_{\mathrm{DF}}+{\mathbf{a}}_{\mathrm{DF}}^{\star },\end{array}$$(20)

where ${\overline{\mathbf{a}}}_{\mathrm{DF}}\equiv ⟨{\mathbf{a}}_{\mathrm{DF}}⟩$ describes the average pull of the perturber and ${\mathbf{a}}_{\text{DF}}^{\star }$ describes the differences from the average. The mean value ${\overline{\mathbf{a}}}_{\mathrm{DF}}$ originates from quite different regions in the host due to its ∝lnr dependence, but the tidal field caused by the wake (∇a_DF) originates mainly from the vicinity of the perturber since contributions from different distances are accumulated with weights ∝r⁻¹. Therefore, in order to recover the total dynamical friction, we need to recover the distortions all around the galaxy, while tidal field recovery is local. In the current approach, we are interested in tidal fields, hence ${\mathbf{a}}_{\text{DF}}^{\star }$, and we may therefore evaluate dynamical friction only within a finite region.

2.3. Determining the extent of orbit integration

To use our dynamical friction framework, one needs to use the condition (11) to decide the time span of orbit integration or t_far. In this section, we examine the basis to find t_far in more detail.

We assume that in any practical application, we do not know the gravitational potential fully and have to rely on some approximation for it. The main ones in the literature are axial or spherical symmetry and equilibrium, but we are also concerned with more fine-grained mismatches. By neglecting asymmetries and the time evolution of the acceleration field in our dynamical friction implementation, we biased the estimation via orbit integration of the perturbed positions x₁. The mismatch between x_{1, recovered} and x₁ scales with the integration time, and an upper limit is $x_{\text{mismatch}}=a_{\sigma }{t}_{\text{far}}^2/2$ in case of a single larger fluctuation or a random walk for many smaller mismatches. We denote the fluctuation levels of the acceleration field as a_σ. We should minimise x_mismatch. That is, we needed to find the minimum value of t_far (or the maximum value of a_far) that satisfies condition (11). Further, we required that the shift in the position of a star due to the perturber (the dynamical friction) be greater than x_mismatch; otherwise, the change in perturbed and unperturbed positions is an artefact of the modelling or implementation assumptions that enter the acceleration description. As $x_1-x_0\gtrsim a_{\text{far}}{t}_{\text{far}}^2/2$, a_far > a_σ, the condition (11) can be rewritten as

$$\begin{array}{c}\hfill |\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{h}}({\mathbf{x}}_0(t_{\mathrm{far}}))|\gg |\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{p}}({\mathbf{x}}_0(t_{\mathrm{far}})|>|\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{h},\mathrm{fluctuations}}|.\end{array}$$(21)

In this way, we minimised mismatches in the positions of stars due to irregularities in the acceleration field of the host galaxy not included in our smooth description that enters the orbit integration in Eqs. (5) and (6). We write the above condition as

$$\begin{array}{c}\hfill a_{\mathrm{far}}=N{a}_{\sigma },\end{array}$$(22)

where N is the free parameter of the modelling, and it depends on the noise level of the gravitational potential of the host (a_σ).

2.4. An alternative orbit integration scheme

To evaluate the tidal dynamical friction with the current method, the condition (11) or (21) must be fulfilled. For stars that are near the host’s centre or have some specific co-movement with the perturber, the condition can never be fulfilled or is always fulfilled. Hence, we needed an alternative for these cases.

To overcome this issue, we modified how these conditions can be fulfilled in a similar manner to Tremaine & Weinberg (1984). These authors introduced the perturber adiabatically to the system by including the term exp(ηt) in the density of the perturber with the η value being very small in order to mimic the buildup of the system. That way, the perturber starts with almost zero initial mass when t ≪ 0 and develops gradually to its final mass when time equals zero. We applied a similar recipe to the acceleration field of the perturber by replacing it with

$$\begin{array}{c}\hfill \mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{pert}}\to \tau (t)\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{pert}}.\end{array}$$(23)

The τ(t) describes multiplicative correction to the acceleration (or, equivalently, density). In this case, the conditions (11) or (21) are fulfilled using the adiabatic introduction, not from physical separation. This fulfilment can be achieved by taking t_far as a fixed time when the perturber had a certain, sufficiently small mass or by retaining individual values of t_far for each particle. The latter possibility does not solve the issue when the conditions are always fulfilled nor when many of the very minor contributions to the acceleration field are included.

This modification poses some changes to the calculation’s algorithm, as it is necessary to specify the form of τ. The selection of the form can have a different basis, depending on the condition, namely, either to complete the backward integration until the perturber is not in the system or, if it is not possible, until some specific time in the past (e.g. until particles that are currently unbound to the perturber were still in its Roche limit).

3. Simulations

To test our dynamical friction framework in realistic conditions, such as the presence of an evolving halo, we built two galaxy N-body simulations using the GADGET-4 tree code (Springel et al. 2021). In one simulation, a galaxy with no perturber acquires the x₀ coordinates, and in the other, a galaxy with a perturber obtains perturbed x₁ coordinates. The initial positions of all particles in both galaxies are identical. We estimated the dynamical friction by comparing the acceleration field of the host particles between the unperturbed and perturbed simulations; see Eq. (4). We denote this version of the dynamical friction values as ’actual’ onward.

3.1. General considerations and initial conditions

We selected a simple model for the simulations so that most effects are due to dynamical friction and not from simulation or model specifics. For the host galaxy, we selected an isothermal sphere model as the observed dynamical mass of galaxies that can be described by isothermal sphere-like distributions (Shajib et al. 2022). The density decreases for the isothermal sphere with the radius as ρ ∝ r⁻², and all the locations have the same velocity dispersion σ. The relationship between the enclosed mass within radius r and the velocity dispersion is M(< r) = 2σ²r/G. The total mass of the host is 10¹¹ M_⊙, distributed in the form of an isothermal sphere within 50 kpc (the cut-off radius) with 4 × 10⁶ particles. Hence, the mass of a particle is 2.5 × 10⁴ M_⊙.

These simulations were aimed at studying the dynamical friction behaviour in various regions, including sparsely populated outskirts. To cope with these regions and verify the method with minimally noisy data (see Sect. 2.2 for a discussion on small-scale fluctuations), we picked a sizeable softening length, b_soft = 1 kpc, to suppress small-scale fluctuations. The softening is described in Springel et al. (2021). The parameters controlling the simulation’s accuracy were integration accuracy (0.0012), maximum time step (10 Myr), and error tolerance to the acceleration (0.005). To see the practical accuracy of the test, we re-ran the simulation backwards, from the end to the beginning, and calculated the shifts and velocity differences between the initial and integrated cases. For 8 Gyr, the shift was ≈1.5 kpc and the velocity difference was ≈1.5 km s⁻¹. As most of our orbit integration is about 0.3 Gyr, we consider this accuracy acceptable.

Initially, we let the simulation reach equilibrium by running it for 2.5 Gyr. The model retained its spherical shape². Then we reset the clock to t = 0 and introduced the perturber described by the softening potential in Springel et al. (2021; i.e. the perturber has a point-like density distribution). As mentioned, in both simulated galaxies, the initial conditions for the particles were identical; thus, the only differences between these simulations concern the perturber with the mass M_pert = 10⁹ M_⊙. The initial conditions for the perturber were x = (0, 60, 0) kpc and v = ( − 40, −20, 0) km s⁻¹. Figure 2 shows the orbit of the perturber.

Fig. 2. Orbit for the perturber with mass Mpert = 109 M⊙. The orange cross shows the point corresponding to the snapshot where the acceleration fields are more thoroughly studied.

3.2. Approximating the acceleration field

As input for our dynamical friction calculations, we used the acceleration field of the simulated host galaxies. For a given simulation snapshot, we approximated the host’s acceleration field by following the steps outlined in the next paragraphs.

We first set the centre of the coordinates (0, 0, 0) to the system’s highest density peak to account for the reflex motion of the host. We then divided the galaxy into radial bins. The mean radial acceleration in each bin was calculated by averaging over all particle accelerations in the bin. The accelerations originate directly from the GADGET-4 output. Finally, we fit a second-order polynomial to the distribution of radial accelerations in each bin separately. The standard deviation of the distribution of radial accelerations minus the fitted polynomial defines the noise in the radial acceleration field in that particular bin as quantified by a_σ.

Through these steps, we derived a smooth, spherical acceleration distribution for the host’s potential. Figure 3 shows the spherical acceleration field for the unperturbed galaxy as a function of distance to the host’s centre. The strength of the estimated noise a_σ is approximately two to three orders of magnitude smaller than the average radial acceleration field in the case of the unperturbed galaxy.

Fig. 3. Comparison of different acceleration components in the simulations. The red line depicts the spherical acceleration field of the unperturbed host’s galaxy at t ≈ 0.8 Gyr (orange cross in Fig. 2). The origin of r is the host’s centre. The blue line shows the estimated noise from the host (see text for details), and the green points show the acceleration from various wake points, as in Fig. 5. For comparison, the perturber causes an acceleration of 43 km2 s−2 kpc−1 at 10 kpc from the system’s centre.

4. Validation

The validation was conducted in two parts. First, we validated the method using the simulations described in the previous section. In this more complex scenario, we quantified the ability of the proposed method to recover the tidal field generated by dynamical friction. Second, we tested our results against Chandrasekhar dynamical friction in Appendix A.

4.1. Tidal field origin

For practical purposes, we divided the dynamical friction evaluation into three distinct parts based on the distance to the perturber. In this way, the Coulomb logarithm can be split into

$$\begin{array}{cc}& ln\mathrm{\Lambda }=ln\frac{b_{\max }}{b_{\min }}=ln\frac{b_{\max }}{R_{\mathrm{dat}}}\frac{R_{\mathrm{dat}}}{b_{\mathrm{soft}}}\frac{b_{\mathrm{soft}}}{b_{\min }}\hfill \\ \hfill & =ln{\mathrm{\Lambda }}_{\mathrm{glob}}+ln{\mathrm{\Lambda }}_{\mathrm{loc}}+ln{\mathrm{\Lambda }}_{\mathrm{unres}},\hfill \end{array}$$(24)

with b_max and b_min as the size of the system and the minimum impact parameter, respectively, b_soft as the softening length, and R_dat as the radius of a sphere centred on the perturber containing the particles used to calculate the dynamical friction. In practical calculations, we added a 5 kpc edge to R_dat when selecting the sample in order to account for particles leaking in or out of R_dat, but for the acceleration calculations, we strictly used the data within it. We define R_dat as the local region around the perturber:

$$\begin{array}{c}\hfill |\mathbf{x}-{\mathbf{x}}_{\mathrm{p}}(t=0)|\le R_{\mathrm{dat}}.\end{array}$$(25)

Figure 4 shows how the different splits capture different amounts of dynamical friction. The fraction of 100% is with R_dat = ∞ capturing all of dynamical friction. The acceleration field is affected by the global perturbations induced by the perturber; therefore, R_dat must be large in order to obtain an unbiased estimate of the acceleration field. The larger R_dat, the less computationally efficient the orbit rewind method will be, where R_dat = ∞ means effectively reconducting the entire simulation. For the current case with a 10 kpc sized region and 50 kpc sized galaxy, the average fraction of stars to evaluate is 0.8%. Nonetheless, the tidal field is less prone to global perturbations, and it converges to its total value for smaller R_dat than for the acceleration field. For this reason, we assessed the accuracy of the method by recovering the tidal field, and we fixed R_dat = 10 kpc, which is sufficient for recovering most of the tides of dynamical friction. Discarding the global part of the dynamical friction estimation leaves the doubt as to whether a sufficient maximum R_dat value is used. From the practical perspective, the question can be answered ‘on the fly’ at every evaluation. Having a set of points x₀ and x₁, the R_dat dependence within some maximum value can be constructed. If the dependence reaches some plateau (as in Fig. 4), then the chosen R_dat is sufficient.

Fig. 4. Recovery fraction of total dynamical friction or its tides in case of 0.3 kpc or 1 kpc softening measured by dimensionless acceleration ratios. The simulation was run with 1 kpc softening; hence, the contributions from regions dampened by softening would be unresolvable for the simulation. To give an estimate of how much is unresolved, we re-evaluated the dynamical friction using a smaller softening length of b = 0.3 kpc. Hence, the difference of the lines between different softening values (0.3 and 1.0 kpc) characterise the contribution of the unresolved part of the dynamical friction. The band around the lines is formed from the estimates from all the simulation snapshots available.

To quantify the recovery of the tidal field at the perturber’s position due to dynamical friction, we compared the fields constructed using the perturbed and unperturbed states. The unperturbed states were obtained directly from the simulation, while the perturbed positions are from the orbit-rewinding of the unperturbed points or from the perturbed simulation.

4.2. Reconstruction properties

The recovery of tidal dynamical friction assumes some choices in the implementation. We used a particle-based description with Eqs. (1) and (2), and we assumed the kernel to be a Plummer sphere (see Eq. (A.1)) with b = 1 kpc to account for softening. The orbit integration was performed using the Runge-Kutta-Felberg method with an adaptive step size while keeping the velocity accuracy below 10⁻⁴ km s⁻¹. The N parameter that controls the extent of orbit integration from Eq. (22) is N = 10 in our fiducial case. The step size for the numerical derivative (see Eq. (18)) is Δ = 2 kpc. We note that it should be larger than the softening radius and similar to the radii around the perturber at which the lopsidedness is to be evaluated. The previous values were varied in order to observe the sensitivity of the recovery. The total number of iterations of the method (see Sect. 2.1) was one, and we did not have to approximate the local dynamical friction estimate with Eq. (4) but instead used the sum directly.

The acceleration field due to dynamical friction at 0.8 Gyr (orange cross in Fig. 2) is shown in Fig. 5. As the perturber orbits in the x-y plane, we only took interest in the x and y components, namely, ∂a_DF, x/∂x, ∂a_DF, x/∂y, ∂a_DF, y/∂x, and ∂a_DF, y/∂y. The comparison of the actual and recovered derivatives as a function of time are shown in Fig. 6.

Fig. 5. Acceleration field from the dynamical friction at 0.8 Gyr from the beginning of the simulation. The point in the orbit is shown in Fig. 2. In the current example, the wake is about 10 kpc behind the perturber and behaving as a sink for aDF.

Fig. 6. Recovery of the tides from dynamical friction for the fiducial case (see Table 2 for numeric characterisation of the recovery). Each panel shows one component of the tides. The ratio of residuals and uncertainty are shown in Fig. 7. The uncertainty (green corridor) was evaluated from the summing in Eqs. (2) and (1) in the recovery case (green line). We note that the uncertainty is not normally distributed. The blue line shows the actual value.

For numerical estimates of the recovery, we used the following estimator

$$\begin{array}{c}\hfill \xi =(\mathrm{\nabla }{\mathbf{a}}_{\mathrm{DF}}^{\mathrm{actual}}-\mathrm{\nabla }{\mathbf{a}}_{\mathrm{DF}}^{\mathrm{recovered}})/{\sigma }_{\mathrm{\nabla }{\mathrm{a}}_{\mathrm{DF}}}.\end{array}$$(26)

Estimates were calculated at each simulation snapshot for each of the derivatives considered. In the equation, σ_∇aDF denotes the uncertainty of $\mathrm{\nabla }{\mathbf{a}}_{\mathrm{DF}}^{\mathrm{recovered}}$ and was evaluated using the bootstrap method (Efron 1979). We did not use the relative error, as the values oscillate around zero, and divisions by small values do not describe the actual recovery. The evolution of ξ values is illustrated in Fig. 7.

Fig. 7. Evolution of the recovery parameter defined with Eq. (26). In other words, this figure shows the residuals of Fig. 6 divided by the uncertainty. This figure complements Table 2.

Figure 7 shows ξ values as a function of time for the considered derivatives. In addition, Table 2 summarises the value of ξ averaged over time and its components. We note that in the case of perfect recovery, the mean of ξ is zero, and the standard deviation of the values is minimal. Standard normal distribution has a standard deviation of one, but we emphasise that due to the existence of high leverage contributors in the sum in Eqs. (1) and (2), the normal assumption is not valid.

We concluded that the current method is able to recover the tides from the wake of dynamical friction. The first stages (before first pericentre passage) have a better quality of recovery, which is expected, as the system remains in a simpler state.

4.3. Verification of the adiabatic introduction

In Sect. 2.4, we pointed out that for some stars, the condition for the extent of the orbit integration cannot be achieved and introduced an alternative to it. In this section, we test this alternative: the adiabatic introduction.

We used the sigmoid function for the τ

$$\begin{array}{c}\hfill \tau =\frac{1}{1+exp[(t-t_0)/\mathrm{\Delta }t]},\end{array}$$(27)

which has two plateaus and a smooth transition region between them. In our implementation, the first plateau corresponds to the time when the perturber was not in the system, and the second plateau corresponds to when asymptotic current perturber mass was achieved. This mimics the buildup and provides a well-defined orbit integration end condition and a constant perturber mass region in order to get the highest accuracy of the part of the wake closest to the perturber, which presumably has the highest impact on the tidal dynamical friction and forms last. In Eq. (27), the t₀ = 0.35 Gyr denotes the mid-point of introducing the perturber, and Δt = 50 Myr is the buildup speed.

The two flavours of the recovery method of the (tidal) dynamical friction are similar. In Fig. 8, we present the dynamical friction values along the path traced by the perturber. The dynamical friction vector is projected to the same trajectory line. The recovered value does not match with the actual one, having an about 30% mismatch and indicating that the global effects influence recovery. However, the slopes of the dynamical friction match well, indicating the tidal field that has a local origin (see Sect. 2.2) is well recovered.

Fig. 8. Dynamical friction as a function of the distance to the perturber along the movement vector of the perturber. The blue line denotes the actual value (see Sect. 4.2), and the golden line (denoted as AdI) shows the value recovered using the adiabatic introduction. The main aspect of this result is that the total values of the dynamical friction do not match as expected from the absence of the global wake. However, the slopes match, indicating that the tides from the dynamical friction are robust in the current case.

Overall, we conclude that the adiabatic introduction recovers the (tidal) dynamical friction as well as the Eq. (21)-based version, but it has both advantages and disadvantages. It has a wider application range and can be used without knowing the noise levels of the galaxy a_σ (in such a case, t_far is shared between all particles). Regarding disadvantages, the integration of orbits is no longer minimal, causing an accumulation of x_mismatch and possibly biasing the results further. The longer orbit integrations also increase calculation time (by about an order of magnitude for the current example).

5. Sampling noise

As shown in Sect. 2.2, in dynamical friction estimates, one also needs to be aware of the sampling noise. As the previous section showed, we were able to recover the local wake of dynamical friction adequately; hence, the current method is able to study how the local sampling noise affects dynamical friction. In this section, we study noise only from a physical basis without any attempt to compensate it with softening, that is, the kernel in Eqs. (1) and (2) is the kernel of the point source.

If we have a region where we want to evaluate the dynamical friction as accurately as possible, a problem arises when we do not observe all the stars, as the dynamical friction evaluation then has some uncertainty and inaccuracy. Even if the overall particle distribution is statistically identical, realisations of this distribution vary, thus varying the estimates of the dynamical friction.

From the simulated galaxy described in Sect. 3, we created different realisations of a sphere with radius R_dat = 2 kpc centred on the perturber at t ≈ 0.8 Gyr (orange cross in Fig. 2) and studied the influence of the sampling noise. This region contains N₀ = 1943 particles, corresponding to a mass density of 0.0014 M_⊙ pc⁻³ (for comparison, the mass density in the solar neighbourhood is 0.097 M_⊙ pc⁻³McKee et al. 2015, or it is the same value as the galactic DM density 24 kpc, McMillan 2017), and has a total mass of M_tot = 5 × 10⁷ M_⊙. The statistical analogues share these parameters and are calculated by sampling from the phase-space distribution function of the simulated particles in that region. Specifically, by placing the perturber at the centre of the reference frame and using spherical coordinates, we first sampled the radial number density and then the angles from the perturber. At the obtained position, we sampled from the constructed radial and tangential velocity distribution functions. Hence, we generated statistical analogues of the region of interest. For each analogue i, we calculated the local acceleration field due to the dynamical friction at the perturber’s position (i.e. $a_{\text{DF}}^i$) to where we added the global contribution from the simulation data. To assess the dependence of the dynamical friction estimate on the number density of the particles, we combined realisations in the following manner. We combined N_stack random realisations and estimated the acceleration field as

$$\begin{array}{c}\hfill a_{\mathrm{stack}}=\frac{1}{N_{\mathrm{stack}}}{\displaystyle \sum _i}{a}_{\mathrm{DF}}^i.\end{array}$$(28)

In this way, the effective particle mass for the stacked set of particles is

$$\begin{array}{c}\hfill m_{\mathrm{particle},\mathrm{eff}.}=\frac{M_{\mathrm{tot}}}{N_0{N}_{\mathrm{stack}}}.\end{array}$$(29)

As we were summing in a_stack, we were able to evaluate its uncertainty by stacking different analogues, similar to bootstrapping. If the m_particle reached the physical average of the stellar initial mass function (IMF), we would reach the noise levels the physical dynamical friction can have.

Figure 9 shows the mean of this distribution of $a_{\text{DF}}^i$ values with the corresponding uncertainty as a function of the particle mass in a fixed physical density. From the figure, we observed that the sampling noise increases with a decreasing number density (larger mass). We note that for some setups, the dynamical friction is physically indeterminable. At a smaller particle mass, the sampling noise introduces a 10% uncertainty at the 1σ level. In case of solar mass particles, the uncertainty is about 23%.

Fig. 9. Dependence on the physical fluctuations in the dynamical friction due to the discreteness the of stars. The aDF describes the estimate of the dynamical friction by combining different analogues of the region. The dark- and light-blue bands represent the 1σ and 2σ uncertainty of the aDF. See Sect. 5 for further details.

We note that the described stacking procedure has a physical limitation: The mass of a single particle cannot fall below the average mass of a star of a brown dwarf. Therefore, Eq. (29) poses a strict limit on how many realisations we can combine or stack. A similar restriction should be posed to numerical integration over phase-space density:

$$\begin{array}{c}\hfill f(\mathbf{x},\mathbf{v})\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{v}\gtrsim ⟨m_{\mathrm{particle}}⟩.\end{array}$$(30)

The ⟨m_particle⟩ denotes the physical mass of an average star (or DM particle in case of halo dynamical friction inference). Suppressing the noise further is not physically justified. Analogously, the spatial element of adaptive integration cannot be smaller than the element containing the mass of the smallest stars in order to keep the integration accuracy compatible with physics.

6. Discussion

This paper has two aims: to introduce a method to calculate local dynamical friction that is usable for perturber distortion evaluation and to show that dynamical friction has intrinsic fluctuations from the sampling noise. As with all dynamical friction methods, the current one has its limitations but also regions where it excels. Method limitations originate mainly from using only the positions of particles in the vicinity of the perturber (within R_dat). Having this limitation, we are only able to include contributors to dynamical friction that are inside of the perturber and not the effects of a global wake. To summarise, the limitations and advantages are following.

The method can capture slightly more than half of the total drag force arising from dynamical friction from a finite R_dat (see Sect. 4.1 and Fig. 4), but at the same time, the method provides input to recover the tidal distortions of the perturber due to the wake of dynamical friction. It assumes the unresolved softening part produces no contributions. Upon a significant suppression of noise from the softening, we recover a larger fraction of the total (tidal) dynamical friction.

The method does not hold any assumption about the host’s symmetry, stationarity, or equilibrium. Hence, the dynamical friction effects can be inferred in very complex environments, and thus the practical applications are limited only by our ability to describe the system fully, including observational limitations. This method can cope with systems without simplifications. For example, if the velocity distribution is not Gaussian but is described as Gaussian, dynamical friction misevaluation alters infall times (and hence dynamical friction) by a factor of 50%−300% (Just et al. 2011), or the absence of the long-term tidal effects alter the friction by about 10% (Roshan & Mashhoon 2022). The tidal effects caused by dynamical friction are not a well-studied topic; hence, we are unable to compare the tidal contributions directly, and we are satisfied with the total dynamical friction. It pushes the limit to evaluate the tidal dynamical friction to cases where complexity is very high, for example, where a pair of perturbers, such as the Large and Small Magellanic Cloud, produce wakes and tides.

We showed that the noise in dynamical friction originates predominantly from the vicinity of the perturber, Sect. 2.2 and Fig. 4, arising from statistical fluctuations of the density field due to a finite number of stars. The current modelling enables us to quantify this effect.

The current method relies heavily on an orbit integration between the past and present time, and one can ask if it would not be simpler to only make two simulations in cases where one wishes to evaluate dynamical friction. For studies that are aimed at examining the behaviour of dynamical friction in various setups, it is a valid approach. But if one wishes to infer the physical aspects of some specific host or perturber, this does not suffice. The main reason is that the initial conditions corresponding to the observations are not known.

The current method is applicable to several studies, with the highest impact potentially on problems that involve the evaluation of tidal dynamical friction in a time- or location-dependent setup.

• Lopsided galaxies: Kipper et al. (2020) have demonstrated that in some cases, a galaxy moving through a medium of uniform density can produce a wake behind it (consisting of unbound particles) that both slow down the galaxy and cause a tidal field, producing lopsidedness in the galaxy. Due to its ability to describe complex situations, the current method allows for diversification of the studies of environments where the lopsidedness can form.

• Simulation recipes: Upon calibrating (using machine learning methods or otherwise) the dependence of precise local dynamical friction on host environmental parameters and incorporating such information into simulations, the calibration will lead to an improvement of the resolution of processes involving the sinking of super massive black holes, or their seeds (Ma et al. 2021), to the centres of galaxies.

• Heating of globular cluster outskirts: Changes in gravitational potential, either through fluctuations in dynamical friction or the change of the location of the wake with respect to the globular cluster, cause heating in stellar systems. The wake is one example of a time-dependent potential. Upon resolving the time dependency of the gravitational potential in the vicinity of globular clusters, the heating can be estimated. The current method allows for the wake-induced potential changes to be evaluated. For example, Wan et al. (2021) showed that processes included in their simulation cannot explain the velocities in the globular cluster outskirts, but the current method can provide some additional possibilities for the explanation.

• Magellanic Clouds: A lot of interest has emerged in studies of the Large Magellanic Cloud due to the possibility that its mass is significantly larger than previous studies have indicated (Erkal et al. 2021). Also, Cullinane et al. (2022) reported the perturbed state of the Magellanic Cloud outskirts. If the dynamical friction is capable of causing the Magellanic Clouds to sink into the Milky Way, it is reasonable to speculate that the wake can alter its outskirts. The current methodology is able to cope with multiple perturbers and evaluate the contribution of wakes to the shape distortions of the perturbers. A similar methodology has been used by Correa Magnus & Vasiliev (2022) to estimate the mass of our Galaxy.

• Tail of Ram-pressure stripping: Galaxies moving through the intergalactic medium are affected by ram-pressure stripping. A jellyfish galaxy tail should coincide with the wake of dynamical friction. The tidal wake caused by the dynamical friction could lead the streams and tails of ejected material to an overdensity that follows the acceleration field Boselli et al. (2022), Triantafyllaki & Longobardi (in prep.). Previous works using N-body-SPH simulations focusing on the ram-pressure scenario Mastropietro et al. (2009) on the Large Magellanic Cloud have tried to explain observations of HI regions and recent star-formation events. Using this method, it is possible to test whether this star formation can be localised in a region where the streams are tugged by the acceleration from the wake.

• DM inference: Once it is possible to evaluate the dynamical friction from the wake’s tidal field of lopsided galaxies, the formation of lopsidedness can be studied in the context of different DM scenarios. Various DM models show differences in the wake (Hui et al. 2017) of the DM; hence, the current method has a potential to distinguish among them based on the effects they cause on the lopsided galaxy.

• Open cluster (OC) stream asymmetry: Pflamm-Altenburg et al. (2023) investigated the asymmetric evaporation of stars in the leading and trailing tidal arms from four open star clusters. With this method, we plan to measure if this asymmetry could be a result of a difference in the dynamical friction effects in the tidal arms. The tidal radius is at most a few hundred parsecs from the clusters, depending on the location and the velocity of the OC. The wake can be at similar distances but not necessarily aligned with the galactic tidal field and should have an overlap with the trailing arm of the stream. In addition, the location of the Lagrangian points becomes altered compared to the isolated Milky Way-OC system, affecting the formation and shape of the stream.

To use the proposed method in real galaxies, one of the concerns is the selection of stars. Namely, this pertains to how to separate the stars belonging to the perturber (or the ones that have just left the perturber via tidal stripping) from the host stars. The strategies for coping with this problem are two-fold, depending on if the individual stars are resolved or not. If the stars are resolved, such as in the case of the Milky Way halo, the split can be done by using either kinematics (co-moving host stars typically have larger velocities with respect to the perturber than the escape velocity from the perturber) or detailed chemical composition (metallicities of the perturber stars are, as a rule, in a narrow range). For example, Battaglia & Starkenburg (2012) found the recipe to separate the interlopers using specific spectral lines. If the single stars of a galaxy are not resolved, the only possible approach is to compare the stellar population properties of the host and the perturber. For example, the J-PAS survey Benitez et al. (2014) provides an excellent set of filters to obtain the basis for the population separation.

7. Summary and conclusions

We developed a flexible semi-analytical method to calculate tidal dynamical friction of a perturber in a complex matter distribution. The current method reformulates dynamical friction as the difference between the acceleration fields of the stars with and without the perturber. Hence, we need two acceleration fields, or particle configurations, to infer dynamical friction. One configuration comes from observations, and the second has to be computed. In Sect. 4 we validated the method.

The central assumption of the proposed method is that for each particle, there was a time in the past (t_far < 0) when the influence of the perturber was insignificant compared to the host. The t_far is obtained by integrating the particle’s and perturber’s orbits towards the past and continuously checking the validity of the condition (21). Once achieved, the particle coordinates at t_far form the initial conditions for calculating the two configurations. The perturbed state corresponds to integrating back to the present in the presence of the perturber potential, and the non-perturbed state corresponds to the absence of the perturber³. After repeating this procedure for all the particles, we obtained the basis for calculating the acceleration field and hence the dynamical friction.

The main conclusions of this study are as follows:

• We designed a method that allows for accurate local dynamical friction estimates in a system with various complexity levels where the phase space density is describable either analytically or as a set of particles.

• We successfully recovered the Chandrashekhar dynamical friction in a homogeneous and continuous medium (as shown in Fig. A.1).

• A drastically different base-assumption set (see Sect. 2.1) allowed us to calculate local dynamical friction in a broad range of possible environments and configurations (see Sect. 6), such as infalling systems.

• The method is limited to evaluating the local aspects of dynamical friction that are the sources of tides from the wake of dynamical friction. These tides lead to shape distortions of the perturber.

• Our assumption set and dynamical friction formulation suggest that the points contributing to dynamical friction have an r⁻³ distance dependence from the perturber, which makes dynamical friction more prone to sampling noise than total acceleration (see Sect. 2.2). The current method is able to distinguish the issues arising from the sampling noise.

• The uncertainty from sampling depends on the number density of particles. In Sect. 5, we showed the sampling noise associated with dynamical friction for various number densities. The fluctuations from the sampling noise do not reach insignificant levels with an increasing number density, showing that the fluctuations are inherent in the nature of the dynamical friction of point sources. The level of noise for our setup with density 0.0014 M_⊙ pc⁻³ has an intrinsic scatter of 10% for dynamical friction (Sect. 5).

Although there are many possible applications of the presented method for computing dynamical friction, this paper covers only one: the intrinsic noise of dynamical friction. We showed that the uncertainty (fluctuations) in dynamical friction calculations is related to the number density of particles. In reality, the number density of stars is not sufficient to firmly calculate dynamical friction. For example, when assuming that a stellar halo consists of stars with the mass of an of an average brown dwarf, there is always a ∼10% uncertainty in dynamical friction (see Fig. 9).

Acknowledgments

We are grateful to referee for useful recommendations and an idea which has immensely improved the paper. The present study was supported by the ETAG projects PRG1006, PSG700, and by the European Structural Funds grant for the Centre of Excellence “The Dark Side of the Universe” (TK133). We applied in this study R statistical environment (Ihaka & Gentleman 1996), C++ and Julia (Bezanson et al. 2017; Danisch & Krumbiegel 2021).

References

Appendix A: Chandrasekhar dynamical friction

We validated our method against the classical Chandrashekhar formulation under the assumptions of an infinite, homogeneous, and isotropic environment. We picked the Chandrasekhar dynamical friction for comparison, as it is an intrinsically local method without any global modes, which our method is capable of reproducing asymptotically. Reproducing Chandrasekhar results allowed us to validate the method in the absence of global modes.

We adopted a setup that mimics a globular cluster in the solar neighbourhood: a perturber with mass M = 2 × 10⁵ M_⊙ is moving with a velocity v = 30 km s⁻¹ relative to the solar neighbourhood in an environment with velocity dispersion σ = 20 km s⁻¹ and density ρ = 0.1 M_⊙pc⁻³. Thus, the perturber is almost co-moving with the rotation of the disc. We characterised the perturber with a Plummer profile

$$\begin{array}{c}\hfill \mathrm{\Phi }=\frac{-GM}{\sqrt{r^2+c^2}},\end{array}$$(A.1)

with c = 0.04 kpc.

The Chandrasekhar dynamical friction equation assumes that a point-massed perturber with mass M and velocity v moves in an infinite, uniform, and isotropic collisionless background. Under these assumptions, the acceleration of the perturber can be calculated as

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\frac{4\pi G^2\rho Mln\mathrm{\Lambda }}{v^3}[\text{erf}\left(\frac{v}{\sqrt{2{\sigma }^2}}\right)-\sqrt{\frac{2}{\pi }}\frac{v}{\sigma }{e}^{-0.5v^2/{\sigma }^2}]\mathbf{v},\end{array}$$(A.2)

where the only unknown parameter of the Chandrasekhar equation is the Coulomb logarithm lnΛ. The slight difference arising from our use of the Plummer profile instead of a point mass is negligible.

We calculated the dynamical friction as outlined in Sect. 2 using a continuous phase-space description of the system. The orbit integration was evaluated using a Runge-Kutta-Fehlberg adaptive step approach. To evaluate Equations (13) and (14) numerically, the integrals had to converge, but they are divergent, just as in the Chandrasekhar case. Therefore, we integrated up to a distance r_∞ = |x(t_far)−x_p(t_far)| = 4 kpc from the perturber, as at the distance r_∞ from the perturber, we consider conditions (11) and (21) to be fulfilled. At that distance, the acceleration from the perturber is about 0.05 km² s⁻² kpc⁻¹; for comparison, the acceleration towards the Galactic centre in the solar neighbourhood is ≈6400 km² s⁻² kpc⁻¹ (Kipper et al. 2021b).

The Chandrasekhar dynamical friction contains an undetermined constant that characterises the extent of the regions in which the dynamical friction is evaluated. As the constant is unknown, the direct values are not comparable; hence, we relied on the recovery of the dependency variables M, v, and σ for validation. Figure A.1 shows a comparison between the estimates of the dynamical friction from the current method and from the Chandrashekhar method when varying the perturber’s mass M (left panel), its velocity (middle panel), and the velocity dispersion of the particles in the background medium (right panel). In our calculations, we assumed the Coulomb logarithm is lnΛ = 0.65ln(r_∞/c) in order to match the normalisation of the Chandrasekhar method to ours. This figure shows, as indicated by the vertical bars, the range of acceleration values within 50 pc around the perturber. From the figure, we observed that the local dynamical friction values are consistent for both methods, and the trends remain robust. In conclusion, the current method successfully recovers the Chandrashekhar dynamical friction in the case of an infinite, homogeneous, and isotropic medium; hence, the local wake is reproduced.

Fig. A.1. Dynamical friction in an infinite, homogeneous, and isotropic medium calculated from the Chandrashekhar method (grey solid line) and our method (blue circles). The vertical bars in all the panels depict the dynamical friction changes in the vicinity of 50 pc from the perturber. We note that they do not describe error or uncertainty. The dynamical friction dependency on the mass (left panel), velocity (central panel), and velocity dispersion (right panel) is recovered.

All Tables

All Figures

Fig. 1. Illustration of the dynamical friction framework. The magenta points describe the points of the host in the vicinity of a perturber (depicted as a blue dot in the enlarged sub-plot). Different coloured lines show a small sub-sample of integrated orbits, both backwards and forwards in time. The dashed circle shows the distance from the perturber where the tfar is reached. The acceleration field difference caused by unperturbed (x0) and perturbed (x1) points gives the dynamical friction. In the enlarged image, we connect the corresponding points with a grey line.

In the text

	Fig. 2. Orbit for the perturber with mass Mpert = 109 M⊙. The orange cross shows the point corresponding to the snapshot where the acceleration fields are more thoroughly studied.
In the text

Fig. 3. Comparison of different acceleration components in the simulations. The red line depicts the spherical acceleration field of the unperturbed host’s galaxy at t ≈ 0.8 Gyr (orange cross in Fig. 2). The origin of r is the host’s centre. The blue line shows the estimated noise from the host (see text for details), and the green points show the acceleration from various wake points, as in Fig. 5. For comparison, the perturber causes an acceleration of 43 km2 s−2 kpc−1 at 10 kpc from the system’s centre.

In the text

Fig. 4. Recovery fraction of total dynamical friction or its tides in case of 0.3 kpc or 1 kpc softening measured by dimensionless acceleration ratios. The simulation was run with 1 kpc softening; hence, the contributions from regions dampened by softening would be unresolvable for the simulation. To give an estimate of how much is unresolved, we re-evaluated the dynamical friction using a smaller softening length of b = 0.3 kpc. Hence, the difference of the lines between different softening values (0.3 and 1.0 kpc) characterise the contribution of the unresolved part of the dynamical friction. The band around the lines is formed from the estimates from all the simulation snapshots available.

In the text

	Fig. 5. Acceleration field from the dynamical friction at 0.8 Gyr from the beginning of the simulation. The point in the orbit is shown in Fig. 2. In the current example, the wake is about 10 kpc behind the perturber and behaving as a sink for aDF.
In the text

Fig. 6. Recovery of the tides from dynamical friction for the fiducial case (see Table 2 for numeric characterisation of the recovery). Each panel shows one component of the tides. The ratio of residuals and uncertainty are shown in Fig. 7. The uncertainty (green corridor) was evaluated from the summing in Eqs. (2) and (1) in the recovery case (green line). We note that the uncertainty is not normally distributed. The blue line shows the actual value.

In the text

	Fig. 7. Evolution of the recovery parameter defined with Eq. (26). In other words, this figure shows the residuals of Fig. 6 divided by the uncertainty. This figure complements Table 2.
In the text

Fig. 8. Dynamical friction as a function of the distance to the perturber along the movement vector of the perturber. The blue line denotes the actual value (see Sect. 4.2), and the golden line (denoted as AdI) shows the value recovered using the adiabatic introduction. The main aspect of this result is that the total values of the dynamical friction do not match as expected from the absence of the global wake. However, the slopes match, indicating that the tides from the dynamical friction are robust in the current case.

In the text

	Fig. 9. Dependence on the physical fluctuations in the dynamical friction due to the discreteness the of stars. The aDF describes the estimate of the dynamical friction by combining different analogues of the region. The dark- and light-blue bands represent the 1σ and 2σ uncertainty of the aDF. See Sect. 5 for further details.
In the text

Fig. A.1. Dynamical friction in an infinite, homogeneous, and isotropic medium calculated from the Chandrashekhar method (grey solid line) and our method (blue circles). The vertical bars in all the panels depict the dynamical friction changes in the vicinity of 50 pc from the perturber. We note that they do not describe error or uncertainty. The dynamical friction dependency on the mass (left panel), velocity (central panel), and velocity dispersion (right panel) is recovered.

In the text