Issue 
A&A
Volume 680, December 2023



Article Number  A91  
Number of page(s)  13  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202347235  
Published online  15 December 2023 
Back to the present: A general treatment for the tidal field from the wake of dynamical friction
^{1}
Tartu Observatory, University of Tartu, Observatooriumi 1, Tõravere 61602, Estonia
email: rain.kipper@ut.ee
^{2}
Department of Physics and the Asher Space Science Institute, Technion – Israel Institute of Technology, Haifa 3200003, Israel
^{3}
Estonian Academy of Sciences, Kohtu 6, 10130 Tallinn, Estonia
Received:
20
June
2023
Accepted:
10
October
2023
Context. Dynamical friction can be a valuable tool for inferring dark matter properties that are difficult to constrain by other methods. Most applications of dynamical friction calculations are concerned with the longterm angular momentum loss and orbital decay of the perturber within its host. This, however, assumes knowledge of the unknown initial conditions of the system.
Aims. We advance an alternative methodology to infer the host properties from the perturber’s shape distortions induced by the tides of the wake of dynamical friction, which we refer to as the tidal dynamical friction.
Methods. As the shape distortions rely on the tidal field that has a predominantly local origin, we present a strategy to find the local wake by integrating the stellar orbits back in time along with the perturber, then removing the perturber’s potential and reintegrating them back to the present. This provides perturbed and unperturbed coordinates and hence a change in coordinates, density, and acceleration fields, which yields the backreaction experienced by the perturber.
Results. The method successfully recovers the tidal field of the wake based on a comparison with Nbody simulations. We show that similar to the tidal field itself, the noise and randomness of the dynamical friction force due to the finite number of stars is also dominated by regions close to the perturber. Stars near the perturber influence it more but are smaller in number, causing a high variance in the acceleration field. These fluctuations are intrinsic to dynamical friction. We show that a stellar density of 0.0014 M_{⊙} kpc^{−3} yields an inherent variance of 10% to the dynamical friction.
Conclusions. The current method extends the family of dynamical friction methods that allow for the inference of host properties from tidal forces of the wake. It can be applied to specific galaxies, such as Magellanic Clouds, with Gaia data.
Key words: galaxies: kinematics and dynamics / dark matter / methods: miscellaneous
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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, multicomponent 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 orbitalbased 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 semianalytically (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 smallscale 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; PflammAltenburg 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 longterm 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 longterm 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 longterm 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 GaravitoCamargo 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 nonequilibrium 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 timeconsuming, 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 largescale structures.
In this paper, we propose a semianalytic 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 phasespace 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 phasespace density. However, this particlebased 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 phasespace and particlebased 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_{0}, and in the presence of the perturber, we denote it as x_{1}. 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_{0} in the absence of the perturber and as a_{1}when the perturber exists:
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:
where r in K_{0} corresponds to x_{0, i} − x for the unperturbed state, while r in K_{1} 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:
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_{0} and a_{1} correspond to the acceleration field of the host or 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 reintegrate 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_{0}(t = 0) and velocities v_{0}(t = 0). The equations of motion describe the movement of particles, where the host potential fully determines the acceleration of these particles:
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. nonstationarity). However, in the presence of the perturber, described with potential Φ_{p} and coordinate x_{p}, the particle accelerates according to
Unfortunately, regarding the use of Eq. (6) to calculate the orbit of a disturbed particle, we do not yet know the potential Φ^{1}. Thus, we calculate it iteratively and replace Φ^{1} with a guessed potential Φ^{1*}, namely,
For the first iteration we take
Taking this approach means we approximate the potentials without taking into account internal deformations of the galaxies, that is, for the first iteration, the selfgravity 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*}
where a_{DF*} is evaluated^{1} using Eq. (2), but the coordinates we use are instead of x_{1}. Upon several iterations, Φ^{1*} approaches Φ^{1}. The perturber’s trajectory x_{p}(t) should also be calculated according to
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_{0}) to the perturbed one (x_{1}), 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:
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_{0}(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},
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_{0}(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 subplot). Different coloured lines show a small subsample of integrated orbits, both backwards and forwards in time. The dashed circle shows the distance from the perturber where the t_{far} is reached. The acceleration field difference caused by unperturbed (x_{0}) and perturbed (x_{1}) 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 particlebased description. In the case of a continuous description of the system, the equations change slightly:
We denote f as the phasespace 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 highleverage 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, onedimensional form for a single particle contribution, we obtain
We use r to denote the distance of the unperturbed position to the evaluated point, and δx = x_{1} − x_{0} is the shift between perturbed and unperturbed positions (i.e. the source of dynamical friction). In case the forceinducing pair is far from the evaluated point (δx ≪ r), the above relation becomes
That is, the dynamical friction asymptotically has a tidallike (r^{−3}) behaviour. If the forceinducing 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
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 σ^{2}, the kernel (gravitational force ∝r^{−2}, dynamical friction ∝r^{−3}, or tidal force from dynamical friction ∝r^{−4}) 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_{pt}w. 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^{2}Δ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 and variance . 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^{2}wdr or σ^{2} ∝ ∑(ΔW_{i})^{2} = ∫r^{2}w^{2}dr. In order to find the signaltonoise ratio, we divided the two .
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 signaltonoise 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 signaltonoise 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^{−1}, 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.
Contribution from all shells to different quantities.
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
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
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:
where describes the average pull of the perturber and describes the differences from the average. The mean value 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^{−1}. 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 , 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 finegrained 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_{1}. The mismatch between x_{1, recovered} and x_{1} scales with the integration time, and an upper limit is 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 , a_{far} > a_{σ}, the condition (11) can be rewritten as
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
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 comovement 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
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 Nbody simulations using the GADGET4 tree code (Springel et al. 2021). In one simulation, a galaxy with no perturber acquires the x_{0} coordinates, and in the other, a galaxy with a perturber obtains perturbed x_{1} 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 spherelike distributions (Shajib et al. 2022). The density decreases for the isothermal sphere with the radius as ρ ∝ r^{−2}, 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σ^{2}r/G. The total mass of the host is 10^{11} M_{⊙}, distributed in the form of an isothermal sphere within 50 kpc (the cutoff radius) with 4 × 10^{6} particles. Hence, the mass of a particle is 2.5 × 10^{4} 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 smallscale fluctuations), we picked a sizeable softening length, b_{soft} = 1 kpc, to suppress smallscale 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 reran 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^{−1}. 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^{2}. 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 pointlike 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^{9} M_{⊙}. The initial conditions for the perturber were x = (0, 60, 0) kpc and v = ( − 40, −20, 0) km s^{−1}. Figure 2 shows the orbit of the perturber.
Fig. 2. Orbit for the perturber with mass M_{pert} = 10^{9} 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 GADGET4 output. Finally, we fit a secondorder 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 km^{2} 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
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:
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_{0} and x_{1}, 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 reevaluated 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 orbitrewinding 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 particlebased 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 RungeKuttaFelberg method with an adaptive step size while keeping the velocity accuracy below 10^{−4} km s^{−1}. 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 xy 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 a_{DF}. 
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
Estimates were calculated at each simulation snapshot for each of the derivatives considered. In the equation, σ_{∇aDF} denotes the uncertainty of 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.
Recovery of tides of dynamical friction.
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 τ
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 welldefined 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} = 0.35 Gyr denotes the midpoint 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_{0} = 1943 particles, corresponding to a mass density of 0.0014 M_{⊙} pc^{−3} (for comparison, the mass density in the solar neighbourhood is 0.097 M_{⊙} pc^{−3}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^{7} M_{⊙}. The statistical analogues share these parameters and are calculated by sampling from the phasespace 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. ) 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
In this way, the effective particle mass for the stacked set of particles is
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 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 a_{DF} describes the estimate of the dynamical friction by combining different analogues of the region. The dark and lightblue bands represent the 1σ and 2σ uncertainty of the a_{DF}. 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 phasespace density:
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 longterm tidal effects alter the friction by about 10% (Roshan & Mashhoon 2022). The tidal effects caused by dynamical friction are not a wellstudied 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 locationdependent 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 timedependent 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 wakeinduced 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 Rampressure stripping: Galaxies moving through the intergalactic medium are affected by rampressure 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 NbodySPH simulations focusing on the rampressure scenario Mastropietro et al. (2009) on the Large Magellanic Cloud have tried to explain observations of HI regions and recent starformation 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: PflammAltenburg 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 WayOC 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 twofold, 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 (comoving 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 JPAS 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 semianalytical 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 nonperturbed state corresponds to the absence of the perturber^{3}. 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 baseassumption 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^{−3} 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^{−3} 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
 Amorisco, N. C. 2017, MNRAS, 464, 2882 [Google Scholar]
 Antoja, T., Helmi, A., RomeroGómez, M., et al. 2018, Nature, 561, 360 [Google Scholar]
 Banik, U., & van den Bosch, F. C. 2021, ApJ, 912, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Banik, U., & van den Bosch, F. C. 2022, ApJ, 926, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, G., & Starkenburg, E. 2012, A&A, 539, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benitez, N., Dupke, R., Moles, M., et al. 2014, ArXiv eprints [arXiv:1403.5237] [Google Scholar]
 Benito, M., Criado, J. C., Hütsi, G., Raidal, M., & Veermäe, H. 2020, Phys. Rev. D, 101, 103023 [NASA ADS] [CrossRef] [Google Scholar]
 Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAM Rev., 59, 65 [Google Scholar]
 Bonetti, M., Bortolas, E., Lupi, A., & Dotti, M. 2021, MNRAS, 502, 3554 [NASA ADS] [CrossRef] [Google Scholar]
 Bose, S., Frenk, C. S., Jenkins, A., et al. 2019, MNRAS, 486, 4790 [NASA ADS] [CrossRef] [Google Scholar]
 Boselli, A., Fossati, M., & Sun, M. 2022, A&ARv, 30, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
 Correa Magnus, L., & Vasiliev, E. 2022, MNRAS, 511, 2610 [NASA ADS] [CrossRef] [Google Scholar]
 Cullinane, L. R., Mackey, A. D., Costa, G. S. D., Koposov, S. E., & Erkal, D. 2022, MNRAS, 518, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Danisch, S., & Krumbiegel, J. 2021, J. Open Source Software, 6, 3349 [NASA ADS] [CrossRef] [Google Scholar]
 Desjacques, V., Nusser, A., & Bühler, R. 2022, ApJ, 928, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Downing, E. R., & Oman, K. A. 2023, MNRAS, 522, 3318 [NASA ADS] [CrossRef] [Google Scholar]
 Efron, B. 1979, Ann. Stat., 7, 1 [Google Scholar]
 Erkal, D., Deason, A. J., Belokurov, V., et al. 2021, MNRAS, 506, 2677 [NASA ADS] [CrossRef] [Google Scholar]
 Fischer, M. S., Brüggen, M., SchmidtHoberg, K., et al. 2022, MNRAS, 516, 1923 [NASA ADS] [CrossRef] [Google Scholar]
 GaravitoCamargo, N., Besla, G., Laporte, C. F. P., et al. 2019, ApJ, 884, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Genel, S., Bryan, G. L., Springel, V., et al. 2019, ApJ, 871, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Gunawardhana, M. L. P., Hopkins, A. M., BlandHawthorn, J., et al. 2013, MNRAS, 433, 2764 [NASA ADS] [CrossRef] [Google Scholar]
 Hartman, S. T. H., Winther, H. A., & Mota, D. F. 2021, A&A, 647, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hashemizadeh, A., Driver, S. P., Davies, L. J. M., et al. 2022, MNRAS, 515, 1175 [NASA ADS] [CrossRef] [Google Scholar]
 Heyvaerts, J., Fouvry, J.B., Chavanis, P.H., & Pichon, C. 2017, MNRAS, 469, 4193 [NASA ADS] [CrossRef] [Google Scholar]
 Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541 [NASA ADS] [CrossRef] [Google Scholar]
 Ihaka, R., & Gentleman, R. 1996, J. Comput. Graphical Stat., 5, 299 [Google Scholar]
 Just, A., Khan, F. M., Berczik, P., Ernst, A., & Spurzem, R. 2011, MNRAS, 411, 653 [Google Scholar]
 Kipper, R., Tenjes, P., Hütsi, G., Tuvikene, T., & Tempel, E. 2019, MNRAS, 486, 5924 [NASA ADS] [CrossRef] [Google Scholar]
 Kipper, R., Benito, M., Tenjes, P., Tempel, E., & de Propris, R. 2020, MNRAS, 498, 1080 [NASA ADS] [CrossRef] [Google Scholar]
 Kipper, R., Tamm, A., Tempel, E., de Propris, R., & Ganeshaiah Veena, P. 2021a, A&A, 647, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kipper, R., Tenjes, P., Tempel, E., & de Propris, R. 2021b, MNRAS, 506, 5559 [NASA ADS] [CrossRef] [Google Scholar]
 Leaman, R., & van de Ven, G. 2022, MNRAS, 516, 4691 [Google Scholar]
 Lovell, M. R., Cautun, M., Frenk, C. S., Hellwing, W. A., & Newton, O. 2021, MNRAS, 507, 4826 [NASA ADS] [CrossRef] [Google Scholar]
 Ma, L., Hopkins, P. F., Ma, X., et al. 2021, MNRAS, 508, 1973 [NASA ADS] [CrossRef] [Google Scholar]
 Mastropietro, C., Burkert, A., & Moore, B. 2009, MNRAS, 399, 2004 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., Parravano, A., & Hollenbach, D. J. 2015, ApJ, 814, 13 [Google Scholar]
 McMillan, P. J. 2017, MNRAS, 465, 76 [Google Scholar]
 Meadows, N., Navarro, J. F., SantosSantos, I., BenítezLlambay, A., & Frenk, C. 2019, MNRAS, 491, 3336 [Google Scholar]
 Moreno, E., FernándezTrincado, J. G., PérezVillegas, A., ChavesVelasquez, L., & Schuster, W. J. 2022, MNRAS, 510, 5945 [NASA ADS] [CrossRef] [Google Scholar]
 Mulder, W. A. 1983, A&A, 117, 9 [NASA ADS] [Google Scholar]
 Nadler, E., DrlicaWagner, A., Bechtol, K., et al. 2021, Phys. Rev. Lett., 126, 091101 [NASA ADS] [CrossRef] [Google Scholar]
 Petts, J. A., Read, J. I., & Gualandris, A. 2016, MNRAS, 463, 858 [NASA ADS] [CrossRef] [Google Scholar]
 PflammAltenburg, J., Kroupa, P., Thies, I., et al. 2023, A&A, 671, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reza, M., Zhang, Y., Nord, B., et al. 2022, Proc. 39th International Conference on Machine Learning (ICML 2022), July 22nd, Baltimore, MD, online at https://ml4astro.github.io/icml2022, 20 [Google Scholar]
 Roper, F. A., Oman, K. A., Frenk, C. S., et al. 2023, MNRAS, 521, 1316 [NASA ADS] [CrossRef] [Google Scholar]
 Roshan, M., & Mashhoon, B. 2022, ApJ, 926, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Shajib, A. J., Vernardos, G., Collett, T. E., et al. 2022, ArXiv eprints [arXiv:2210.10790] [Google Scholar]
 Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729 [Google Scholar]
 Vasiliev, E. 2023, MNRAS, 527, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978 [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279 [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E., Belokurov, V., & Evans, N. W. 2022, ApJ, 926, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Wan, Z., Oliver, W. H., Baumgardt, H., et al. 2021, MNRAS, 502, 4513 [NASA ADS] [CrossRef] [Google Scholar]
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^{5} M_{⊙} is moving with a velocity v = 30 km s^{−1} relative to the solar neighbourhood in an environment with velocity dispersion σ = 20 km s^{−1} and density ρ = 0.1 M_{⊙}pc^{−3}. Thus, the perturber is almost comoving with the rotation of the disc. We characterised the perturber with a Plummer profile
with c = 0.04 kpc.
The Chandrasekhar dynamical friction equation assumes that a pointmassed 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
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 phasespace description of the system. The orbit integration was evaluated using a RungeKuttaFehlberg 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^{2} s^{−2} kpc^{−1}; for comparison, the acceleration towards the Galactic centre in the solar neighbourhood is ≈6400 km^{2} s^{−2} kpc^{−1} (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 subplot). Different coloured lines show a small subsample of integrated orbits, both backwards and forwards in time. The dashed circle shows the distance from the perturber where the t_{far} is reached. The acceleration field difference caused by unperturbed (x_{0}) and perturbed (x_{1}) 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 M_{pert} = 10^{9} 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 km^{2} 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 reevaluated 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 a_{DF}. 

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 a_{DF} describes the estimate of the dynamical friction by combining different analogues of the region. The dark and lightblue bands represent the 1σ and 2σ uncertainty of the a_{DF}. 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.