© The Authors 2023

1. Introduction

Dynamical friction (DF) is an important physical phenomenon that has several consequences in stellar dynamics (and in plasma physics). It can be qualitatively thought of as the slowing-down of a test particle of mass M, moving at v_T in a background of field particles of mass m, mean number density n, and velocity distribution f(v_F), due to the cumulative effect of their long-range gravitational (or Coulomb) interactions.

An analytical estimation of the DF was evaluated for the first time for stellar systems by Chandrasekhar (1943), who found that M must experience a slowing down along its initial direction of propagation as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}=-4\pi G^{2}nm(M+m)log\mathrm{\Lambda }\frac{\mathrm{\Xi }(v_{\mathrm{T}})}{v{3}_{\mathrm{T}}}{\mathit{v}}_{\mathrm{T}}.\end{array}$$(1)

In this equation, G is the gravitational constant, logΛ is the Coulomb logarithm (analogously with the same quantity in plasma physics, Spitzer 1965), of the ratio Λ of the maximum and minimum impact parameters b_max and b_min and

$$\begin{array}{c}\hfill \mathrm{\Xi }(v_{\mathrm{T}})\equiv 4\pi {\int }_0^{v_{\mathrm{T}}}f(v_{\mathrm{F}})v_{\mathrm{F}}^2\mathrm{d}{v}_{\mathrm{F}}\end{array}$$(2)

is the fractional velocity volume function.

The process of DF is crucial for the evolution of collisional systems (see e.g. Binney & Tremaine 1987), from the large scales of galaxies clusters (Ostriker & Tremaine 1975; Richstone 1976; Gunn & Tinsley 1976; Adhikari et al. 2016) to the smaller scales of dense stellar systems, for its consequences on the motion of supermassive black holes (SMBHs) in galactic cores (Antonini & Merritt 2012; Tremmel et al. 2018; Di Cintio et al. 2020; Ricarte et al. 2021; Chen et al. 2022), globular clusters (GCs) orbiting their host galaxies (Weinberg 1989; Colpi & Pallavicini 1998; Bertin et al. 2003; Arena et al. 2006; Arena & Bertin 2007), or exotic stellar objects, for example blue straggler stars (BSS, Ferraro et al. 1995) in GCs (see Paresce et al. 1992, Ferraro et al. 1995, 2001, 2009; Procter Sills et al. 1995; Ransom et al. 2005; Pooley & Hut 2006; Alessandrini et al. 2014, 2016; Miocchi et al. 2015; Pasquato et al. 2018; Pasquato & Di Cintio 2020).

Since the pioneering work of Chandrasekhar, the DF formalism, initially conceived for a point-like particle in an infinitely extended background of scatterers, has been extended to the case of finite-sized objects (see e.g. Mulder 1983; Zelnikov & Kuskov 2016) sinking in the host stellar system, flattened or spherical models with self-gravity (see e.g. Kalnajs 1971; Tremaine & Weinberg 1984), or spheroids with anisotropic velocity distribution (Binney 1977).

More recently, Ciotti & Binney (2004) and Nipoti et al. (2008) derived an expression for the DF formula in the case of modified Newtonian dynamics (MOND, Milgrom 1983; Bekenstein & Milgrom 1984), and performed N-body simulations of sinking satellites in MOND, while Ciotti (2010) considered the effect of a mass spectrum for the field particles and Silva et al. (2016) that of a non-thermal power-law-like velocity distribution. The main conclusion of these works is that, in general, using the original formulation of the DF for idealized infinite systems (see Eq. (1)), leads to a substantial underestimation (even by a factor of 10) of its effectiveness when applied to more realistic models.

Prompted by the detection of gravitational waves (GWs) from binary compact objects announced by the LIGO and VIRGO collaborations Abbott et al. (2016a,b), a renewed interest in relativistic stellar dynamics (Shapiro & Teukolsky 1985; Hamers et al. 2014) has recently emerged, in particular with respect to the formation and migration processes of single and binary black holes (BHs) in GCs (Samsing & D’Orazio 2018; Rodriguez et al. 2018; Antonini et al. 2019; Torniamenti et al. 2022), supermassive BHs in galactic cores (Merritt 2015; Fang & Huang 2020; Kelley et al. 2021; Liu & Lai 2022), or runaway objects from dense star clusters (Ryu et al. 2017; Farias et al. 2020; Bhat et al. 2022).

From the theoretical point of view, if on the one hand the distribution function-based approach to relativistic stellar dynamics has been widely explored since Fackerell (1968) made a first attempt (see e.g. Katz et al. 1975; Ellis et al. 1983; Israel & Kandrup 1984; Kandrup 1994, 1984; Kandrup & Morrison 1993; Chavanis 2020a,b and references therein), on the other hand much less has been done in the context of relativistic collisional systems in the original Chandrasekhar picture. Lee (1969) formally extended Eq. (1) accounting for the first post-Newtonian corrections to the gravitational force, though did not evaluate it for an explicit choice of the velocity distribution f(v_F) or include the effects of relativistic velocities. Nevertheless, he concluded that DF in a dense stellar system, where strong deflections induced by small impact parameters may happen, is always enhanced when considering corrections to the Newtonian force of order 1/c², where c is the speed of light.

Syer (1994) derived an alternative expression for the DF in the case of a relativistic test particle crossing an isotropic medium of lighter particles with m ≪ M in the limit of small scattering angles. More recently, Barausse (2007), Katz et al. (2019), Traykova et al. (2021), Vicente & Cardoso (2022), and Correia (2022) explored the onset of DF in relativistic fluids where the test particle induces a wake, and Cashen et al. (2017) included the precession effect induced by the contribution of the gravitomagnetic effect (see e.g. Costa & Natário 2014 and references therein).

In this work we explore this matter further. Our aim is to formulate a general treatment of DF that can be applied in different regimes dominated by collisional effects, involving very large test particle masses and/or relativistic background systems.

The paper is structured as follows. In Sect. 2 we introduce the astrophysical systems where relativistic collisional processes are relevant. In Sect. 3 we revise the classical formalism of DF based on the hyperbolic two-body problem to introduce the fundamental quantities used in the following. In Sect. 4 we derive the (special) relativistic DF expression, complete with a relativistic velocity distribution. In Sect. 5 we compute the post-Newtonian 1/c² corrections to the classical DF expression, and we show an alternative derivation of the relativistic post-Newtonian DF expression using the Darwin Lagrangian. Finally, in Sect. 6 we summarize our findings and outline the possible applications of this work.

2. Relativistic stellar dynamical systems

When discussing relativistic effects in stellar dynamics (in particular with respect to DF) three set-ups can be considered: (i) the test mass M moves at high speed v so that its relativistic (Lorentz) factor

$$\begin{array}{c}\hfill {\gamma }_v=\frac{1}{\sqrt{1-\frac{v^2}{c^2}}}\end{array}$$(3)

becomes significantly larger than unity; (ii) the test mass M is so large that some degree of post-Newtonian approximation should be used when resolving the close approaches with the background stars; (iii) the velocity distribution of the system f(v) has non-negligible relativistic tails.

An important example of the first case is represented by hypervelocity stars (HVSs). These could be produced during multiple strong interactions between the central galactic SMBH and an infalling GC, which can accelerate many stars belonging to the GC to high velocities, and can even eject them in jets from the inner galactic region. This usually occurs as a result of three- or four-body interactions with the SMBH (or SMBH-binary). Alternatively, dynamical kicks due to supernova explosions or close encounters of a hard massive binary star and a single massive star could accelerate stars to extreme velocities. In addition, some galactic HVSs are thought to be the result of a merger event with another nearby galaxy with a high-velocity relative to the present galactic environment (see Fragione & Capuzzo-Dolcetta 2016).

Typically, HVSs observed in the galactic halo have velocities that can reach ∼1.2 × 10³ km s⁻¹; the typical stellar velocity in the Galaxy is roughly 100 km s⁻¹. Therefore, the size of the relativistic corrections, quantified in γ − 1, on a HVS are of the order of 8 × 10⁻⁶, while for an average star it is 5.6 × 10⁻⁸.

The second case is exemplified by a massive black hole moving through a dense stellar system, for example a nuclear star cluster (NSCs) or a core-collapsed GC. BHs can be accelerated to large velocities, for instance during the final stage of SMBHs coalescence by the large recoil due to anisotropic emission of GWs, with v_recoil ≈ 10³ km s⁻¹ for the coalesced object, that may also be displaced from the minimum of the host galaxy potential-well, or even ejected Lena et al. (2014), Kim et al. (2017).

Nuclear star clusters have typical mass densities in the range 3 × 10⁴ − 2 × 10⁵ M_⊙ pc⁻³. Assuming an average stellar mass of about 0.5 M_⊙ implies that their average interstar distance r_int is of the order of 2 × 10⁻² pc, while the Schwarzschild radius r_s = 2GM_BH/c² for a 10⁵ − 10⁶ M_⊙ massive BH is roughly 10⁻⁷ − 10⁻⁶ pc. Therefore, a close encounter with a star at about 10⁻²r_int happens at ∼10²r_s, where relativistic corrections should be taken into account.

An example of the third case is the distribution of dark matter in proximity of a massive black hole. In general, in the Lambda-cold dark matter (ΛCDM) paradigm, dark matter is supposed to be in a non-relativistic and collisionless regime (see e.g. Binney & Tremaine 2008). However, in the presence of the deep gravitational potential of a BH, (dark) matter particles could reach relativistic velocities (Zamir 1993). Quantifying the relativistic corrections to DF in a dark matter cusp with a central BH is therefore worth exploring in the context of seeding mechanisms on primordial BHs (Kavanagh et al. 2020; Cole et al. 2023).

3. Dynamical friction: Classical case

Before tackling the problem of its relativistic generalization, it is instructive to revisit the classical treatment of the DF and retrive Eq. (1). As usual, we consider each single encounter between the test particle and the field particles as a hyperbolic two-body problem in the frame centred on the field particle. Let (x_T, v_T) and (x_F, v_F) be the positions and velocities of M and m, respectively, and let

$$\begin{array}{c}\hfill \mathit{r}={\mathit{x}}_{\mathrm{T}}-{\mathit{x}}_{\mathrm{F}};\mathit{V}=\dot{\mathit{r}}={\mathit{v}}_{\mathrm{T}}-{\mathit{v}}_{\mathrm{F}}\end{array}$$(4)

be their relative position and velocity. We recall that the equation of motion for a fictitious particle of reduced mass μ = mM/(M + m) moving in the Keplerian potential of the fixed body of mass M + m, is

$$\begin{array}{c}\hfill \frac{\mathit{mM}}{m+M}\ddot{\mathit{r}}=-\frac{\mathit{GMm}}{r^2}{\widehat{\mathit{e}}}_{\mathit{r}}.\end{array}$$(5)

The energy conservation along the orbit of μ for a given encounter with impact parameter b, implies that the relative velocity vector V is deflected by an angle π − 2ψ, in the orbital plane defined by

$$\begin{array}{c}\hfill cos\psi =\frac{1}{\sqrt{1+\frac{b^2{V}^4}{G^2{(M+m)}^2}}}.\end{array}$$(6)

Using finite differences, we can always express the relative velocity change as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathit{V}=\mathrm{\Delta }{\mathit{v}}_{\mathrm{F}}-\mathrm{\Delta }{\mathit{v}}_{\mathrm{T}},\end{array}$$(7)

where Δv_F and Δv_T are the velocity variations of m and M during the encounter. Since the velocity of the centre of mass is constant (by definition) during the encounter, we have that

$$\begin{array}{c}\hfill m\mathrm{\Delta }{\mathit{v}}_{\mathrm{F}}+M\mathrm{\Delta }{\mathit{v}}_{\mathrm{T}}=0.\end{array}$$(8)

Eliminating Δv_F in the two equations above yields

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{v}}_{\mathrm{T}}=-\left(\frac{m}{m+M}\right)\mathrm{\Delta }\mathit{V}.\end{array}$$(9)

We must now evaluate ΔV in order to find Δv_T. The conserved angular momentum per unit mass of the reduced particle is L = bV. Let us now label with θ_defl the deflection angle. The relation between the radius and azimuthal angle of a particle on a Keplerian orbit becomes

$$\begin{array}{c}\hfill \frac{1}{r}=Ccos(\psi -{\psi }_0)+\frac{G(M+m)}{b^2{V}^2},\end{array}$$(10)

where the constant C and the phase angle ψ₀ = ψ(t = 0) are determined by the initial conditions. Deriving (10) with respect to the time we obtain

$$\begin{array}{c}\hfill \frac{\mathrm{d}r}{\mathrm{d}t}=C{r}^2\dot{\psi }sin(\psi -{\psi }_0)=CbVsin(\psi -{\psi }_0),\end{array}$$(11)

where the last term arises from $L=r^2\dot{\psi }$. If we impose that ψ = 0 when t → −∞, we obtain from (11)

$$\begin{array}{c}\hfill -V=Cbsin(-{\psi }_0).\end{array}$$(12)

Evaluating Eq. (10) we then have

$$\begin{array}{c}\hfill 0=Ccos({\psi }_0)+\frac{G(M+m)}{b^2{V}^2},\end{array}$$(13)

and eliminating C from the equations above we obtain

$$\begin{array}{c}\hfill tan{\psi }_0=-\frac{b{V}^2}{G(M+m)}.\end{array}$$(14)

From Eqs. (10) and (11) we find that the point of closest approach is reached when ψ = ψ₀ and, since the orbit is symmetrical about this point, the deflection angle is θ_defl = 2ψ₀ − π. Thanks to the conservation of energy, after the encounter, the modulus of the relative velocity, equals the modulus of the initial relative velocity, and therefore the components of ΔV (perpendicular and parallel to the initial relative velocity vector V: ΔV_∥ and ΔV_⊥)

$$\begin{array}{c}\hfill ||\mathrm{\Delta }{\mathit{V}}_{\perp }||=\frac{2b{V}^3}{G(M+m)}[1+\frac{b^2{V}^4}{G^2{(M+m)}^2}{]}^{-1};\end{array}$$(15)

and

$$\begin{array}{c}\hfill ||\mathrm{\Delta }{\mathit{V}}_{\Vert }||=2V[1+\frac{b^2{V}^4}{G^2{(M+m)}^2}{]}^{-1}.\end{array}$$(16)

In a homogeneous background of particles equal masses m, all Δv_T⊥ sum to zero by symmetry, (using the ‘Jeans swindle’; see e.g. Binney & Tremaine 1987), while the parallel velocity changes add up, and thus the mass M will experience a deceleration (see Eq. (9)) as a result of the DF. Therefore, it is sufficient to evaluate Δv_T∥ as

$$\begin{array}{c}\hfill ||\mathrm{\Delta }{\mathit{v}}_{\mathrm{T}\Vert }||=\frac{2mV}{M+m}[1+\frac{b^2{V}^4}{G^2{(M+m)}^2}{]}^{-1}.\end{array}$$(17)

In a system defined by the phase-space distribution function F = nf(v_F), where n is a constant number density and f(v_F) is the velocity distribution, the rate at which the mass M encounters stars with impact parameter between b and b + db, and velocities between v_F and v_F + dv_F, is

$$\begin{array}{c}\hfill n_{\mathrm{enc}}=2\pi b\mathrm{d}bVnf({\mathit{v}}_{\mathrm{F}})d^3{\mathit{v}}_{\mathrm{F}},\end{array}$$(18)

where d³v_F is the velocity-space element. The total change in velocity undergone by M is found by adding all the contributions of ||Δv_T∥|| due to particles with impact parameters from b_min to a b_max, and then summing over all velocities of stars as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}{|}_{{\mathit{v}}_{\mathrm{F}}}& =\mathit{V}nf({\mathit{v}}_{\mathrm{F}})d^3{\mathit{v}}_{\mathrm{F}}{\int }_{b_{\min }}^{b_{\max }}||\mathrm{\Delta }{\mathit{v}}_{\mathrm{T}\Vert }||2\pi b\mathrm{d}b\hfill \\ \hfill & =\mathit{V}nf({\mathit{v}}_{\mathrm{F}})d^3{\mathit{v}}_{\mathrm{F}}{\int }_{b_{\min }}^{b_{\max }}\frac{2mV}{M+m}[1+\frac{b^2{V}^4}{G^2{(M+m)}^2}{]}^{-1}2\pi b\mathrm{d}b.\hfill \end{array}$$(19)

We first perform the integral over b:

$$\begin{array}{cc}& {\int }_{b_{\min }}^{b_{\max }}\frac{2mV}{M+m}[1+\frac{b^2{V}^4}{G^2{(M+m)}^2}{]}^{-1}2\pi b\mathrm{d}b\hfill \\ \hfill & =\frac{2\pi G^2(M+m)}{V^3}log\left[\frac{1+\frac{b_{\max }^2{V}^4}{G^2{(M+m)}^2}}{1+\frac{b_{\min }^2{V}^4}{G^2{(M+m)}^2}}\right].\hfill \end{array}$$(20)

We note that choosing the minimum and maximum impact parameters is a rather delicate step. When using the impulsive approximation (i.e. μ||ΔV_⊥|| = 2GMm/bV; see e.g. Ciotti 2010), the ‘ultraviolet divergence’ occurs when performing the integral over b in (19) and setting b_min = 0. This divergence is actually artificial as it disappears when the full solution of the hyperbolic two-body problem is taken into account as in this discussion. Conversely, the ‘infrared divergence’ appearing for b → ∞ cannot be eliminated in an infinite system, and therefore an upper cutoff must be used with a suitable choice of b_max.

Not surprisingly, this point is still a source of debate (see e.g. the discussion in Van Albada & Szomoru 2020 and references therein). The problem lies in what should dominate between a few strong encounters with nearby stars (see Chandrasekhar 1941, 1942, 1943, see also Kandrup 1983) or many weak encounters with distant stars (see Spitzer 1987; Binney & Tremaine 1987). In the first interpretation b_max should be of the order of the average interparticle distance, while in the second, it should be of the order of the size of the system. Both views are in principle plausible, the former as it is more intuitive to think that the largest contribution must be due to the nearest stars and the latter because there is no screening in gravitational systems, at variance with (quasi-)neutral plasmas where charges of opposite signs are present. In this work we follow the Spitzer approach.

We note that, under most conditions of practical interest in astrophysics, the quantity Λ² = b²V⁴/G²(M + m)² is typically much greater than unity. For this reason, from now on we replace log1 + Λ² with 2logΛ, recovering the widely used definition of the Coulomb logarithm as $log[b_{\max }^2{V}^4/G^2{(M+m)}^2]$¹.

In this approximation, the right-hand side of Eq. (20) becomes

$$\begin{array}{c}\hfill \frac{2\pi G^2(M+m)}{V^3}log[\frac{1+\frac{b_{\max }^2{V}^4}{G^2{(M+m)}^2}}{1+\frac{b_{\min }^2{V}^4}{G^2{(M+m)}^2}}]\approx \frac{4\pi G^2(M+m)}{V^3}log\mathrm{\Lambda }.\end{array}$$(21)

Combining the expression above with Eq. (19) yields

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}=-4\pi G^{2}nm(M+m)log\mathrm{\Lambda }\int f({\mathit{v}}_{\mathrm{F}})\frac{({\mathit{v}}_{\mathrm{T}}-{\mathit{v}}_{\mathrm{F}})}{||{\mathit{v}}_{\mathrm{T}}-{\mathit{v}}_{\mathrm{F}}{||}^3}{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}},\end{array}$$(22)

where we have replaced V with its definition (see. Eq. (4)) and assumed logΛ as the velocity averaged Coulomb logarithm (see e.g. the discussion in Ciotti 2021). The velocity integral in Eq. (22) is often referred to as the first Rosenbluth potential (see e.g. Rosenbluth et al. 1957).

Remarkably, the problem of computing the acceleration dv_T/dt integrating over all field star velocities is formally equivalent to that of evaluating the gravitational field at v_T generated by the mass density: ρ(v_F) = 4πlogΛGm(M + m)f(v_F). Assuming an isotropic (spherically symmetric) velocity distribution, in virtue of Newton’s second theorem (see e.g. Chandrasekhar 1995) we have that only the stars such that v_F < v_T contribute to the slowing down of M, and hence

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}=-16{\pi }^2{G}^{2}nm(M+m)log\mathrm{\Lambda }\frac{{\mathit{v}}_{\mathrm{T}}}{v{3}_{\mathrm{T}}}{\int }_0^{v_{\mathrm{T}}}f(v_{\mathrm{F}})v{2}_{\mathrm{F}}\mathrm{d}{v}_{\mathrm{F}}.\end{array}$$(23)

In the special case where f(v_F) is a Maxwellian with dispersion σ,

$$\begin{array}{c}\hfill f(v_{\mathrm{F}})=\frac{v_{\mathrm{F}}^2}{{(2\pi {\sigma }^2)}^{3/2}}exp(-\frac{v_{\mathrm{F}}^2}{2{\sigma }^2}),\end{array}$$(24)

evaluating the velocity volume function integral in Eq. (23) yields

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}=-4\pi G^2(M+m)\rho log\mathrm{\Lambda }[\mathrm{Erf}\left(\frac{v_{\mathrm{T}}}{\sqrt{2}\sigma }\right)-\frac{2v_{\mathrm{T}}{e}^{-\frac{v_{\mathrm{T}}^2}{2{\sigma }^2}}}{\sqrt{2\pi }\sigma }]\frac{{\mathit{v}}_{\mathrm{T}}}{v_{\mathrm{T}}^3},\end{array}$$(25)

where we have condensed nm in ρ (the mean mass density of the field particles), and Erf(x) is the standard error function defined (see Arfken et al. 2012) as

$$\begin{array}{c}\hfill \mathrm{Erf}(x)=\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{-t^2}\mathrm{d}t.\end{array}$$(26)

4. Dynamical friction: Special relativistic generalization

As mentioned above, the Chandrasekhar DF formula was extended to relativistic velocities by Syer (1994), but only in the weak scattering limit (i.e. b ≫ r_s ≈ GM/c² and small deflection angle θ_defl). We will now derive a more general expression for a generic θ_defl, therefore also accounting for the case of strong scattering, (i.e. θ_defl → π/2).

In this derivation we keep the classical 1/r² Newtonian force and replace velocity composition with its relativistic counterpart. For this purpose we define an inertial frame S′, in which the test star of mass M is stationary at the beginning of the encounter (i.e. the field star m is at infinity). In this frame, θ_defl = π − 2ψ (with ψ given by Eq. (6)) is the deflection angle according to the classical unbound two-body problem. This assumption is motivated by the fact that in the astrophysically relevant case where M > m, the relativistic scattering angle in an elastic collision has the same majorant as the classical case (see Landau & Lifshitz 1976, Chap. 2) given by sin θ_defl, max = m/M (see also Thornton & Marion 2004). In this approximation, however, we do not assume a small angle limit (i.e. weak scattering) as no assumption has been made on the relative angles between v_T and v_F.

In the (special) relativistic encounter the relative velocity V becomes (see e.g. Landau & Lifshitz 1976)

$$\begin{array}{c}\hfill V^2=\frac{||{\mathit{v}}_{\mathrm{T}}-{\mathit{v}}_{\mathrm{F}}{||}^2-\frac{1}{c^2}||{\mathit{v}}_{\mathrm{T}}\wedge {\mathit{v}}_{\mathrm{F}}{||}^2}{(1-\frac{{\mathit{v}}_{\mathrm{T}}·{\mathit{v}}_{\mathrm{F}}}{c^2}{)}^2},\end{array}$$(27)

for arbitrary choices of v_T and v_F. We note that, in the relativistic case, the velocity v_A|B of a body A with respect to another B is not equal to −v_B|A of B with respect to A. This loss of symmetry is related to the Thomas (1927) precession², and the fact that two subsequent Lorentz transformations rotate the coordinate system (see Weinberg 1972). This rotation, however, conveniently has no effect on the magnitude of a vector, and hence the modulus of the relative velocity is symmetrical.

We let v_F, μ = γ_vF(c, v_F) and v_T, μ = γ_vT(c, v_T) be the four-velocities of the field and test stars, respectively, and where ${\gamma }_{v_{\text{F};\text{T}}}={(1-v_{\text{F};\text{T}}^2/c^2)}^{-1/2}$ are the Lorentz factors of M and m. We consider a single encounter in the ‘laboratory frame’ S. This process can be expressed as a product of a Lorentz boost Γ in S′, a rotation ℛ(θ_defl) in the three-space and, finally, an inverse boost Γ⁻¹, reverting back to S.

Defining p_Fμ = mγ_vF(c, v_F) as the four-momentum of m before the encounter, we have that ${{\mathit{p}}_{\mathrm{F}}\prime }^{\mu }={\mathrm{\Lambda }}_{\nu }^{\mu }{\mathit{p}}_{\mathrm{F}}^{\nu }$, where Λ = Γ⁻¹ℛ(θ_defl)Γ, and thus we formally obtain

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{p}}_{\mathrm{F}}^{\mu }=({\mathrm{\Gamma }}^{-1}\mathcal{R}({\theta }_{\mathrm{defl}})\mathrm{\Gamma }-1){\mathit{p}}_{\mathrm{F}}^{\mu }.\end{array}$$(28)

Since the motion is planar, as we are still dealing with a classical two-body problem, we can simplify the notation involving four vectors by using three vectors instead, where only two of the space dimensions are maintained, one parallel to v_T and one perpendicular. However, as argued before, by reasons of symmetry, only the parallel component of V contributes to the DF. Denoting with ϕ the angle between v_T and v_F, we can now write

$$\begin{array}{c}\hfill {\mathit{p}}_{\mathrm{F}}^{\mu }=m{\gamma }_{v_{\mathrm{F}}}\left(\begin{array}{c}c\\ v_{\mathrm{F}}cos\varphi \\ v_{\mathrm{F}}sin\varphi \end{array}\right).\end{array}$$(29)

With this choice of Γ, ℛ and Γ⁻¹ read

$$\begin{array}{cc}\hfill \mathrm{\Gamma }& =\left(\begin{array}{ccc}{\gamma }_{v_{\mathrm{T}}}& -\frac{v_{\mathrm{T}}}{c}{\gamma }_{v_{\mathrm{T}}}& 0\\ -\frac{v_{\mathrm{T}}}{c}{\gamma }_{v_{\mathrm{T}}}& {\gamma }_{v_{\mathrm{T}}}& 0\\ 0& 0& 1\end{array}\right),\mathcal{R}({\theta }_{\mathrm{defl}})=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos{\theta }_{\mathrm{defl}}& -sin{\theta }_{\mathrm{defl}}\\ 0& sin{\theta }_{\mathrm{defl}}& cos{\theta }_{\mathrm{defl}}\end{array}\right),\hfill \\ \hfill {\mathrm{\Gamma }}^{-1}& =\left(\begin{array}{ccc}{\gamma }_{v_{\mathrm{T}}}& \frac{v_{\mathrm{T}}}{c}{\gamma }_{v_{\mathrm{T}}}& 0\\ \frac{v_{\mathrm{T}}}{c}{\gamma }_{v_{\mathrm{T}}}& {\gamma }_{v_{\mathrm{T}}}& 0\\ 0& 0& 1\end{array}\right).\hfill \end{array}$$(30)

After the encounter the four-momentum of the field star ${\mathit{p}}_{\mathrm{F}}^{\mu }$ changes by $\mathrm{\Delta }{\mathit{p}}_{\mathrm{F}}^{\mu }$ defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{p}}_{\mathrm{F}}^{\mu }=m{\gamma }_{v_{\mathrm{F}}}\left(\begin{array}{c}A\\ B\\ C\end{array}\right),\end{array}$$(31)

where, respectively,

$$\begin{array}{cc}\hfill A& =c{\gamma }_{v_{\mathrm{T}}}^2-\frac{v_{\mathrm{T}}{v}_{\mathrm{F}}}{c}{\gamma }_{v_{\mathrm{T}}}^{2}cos\varphi +\frac{v_{\mathrm{T}}}{c}{\gamma }_{v_{\mathrm{T}}}[{\gamma }_{v_{\mathrm{T}}}(v_{\mathrm{F}}cos\varphi -v_{\mathrm{T}})\hfill \\ \hfill & \times cos{\theta }_{\mathrm{defl}}-v_{\mathrm{F}}sin\varphi sin{\theta }_{\mathrm{defl}}]-c,\hfill \\ \hfill B& =v_{\mathrm{T}}{\gamma }_{v_{\mathrm{T}}}^2(1-\frac{v_{\mathrm{T}}{v}_{\mathrm{F}}}{c^2}cos\varphi )+{\gamma }_{v_{\mathrm{T}}}[{\gamma }_{v_{\mathrm{T}}}(v_{\mathrm{F}}cos\varphi -v_{\mathrm{T}})\hfill \\ \hfill & \times cos{\theta }_{\mathrm{defl}}-v_{\mathrm{F}}sin\varphi sin{\theta }_{\mathrm{defl}}]-v_{\mathrm{F}}cos\varphi ,\hfill \\ \hfill C& ={\gamma }_{v_{\mathrm{T}}}(v_{\mathrm{F}}cos\varphi -v_{\mathrm{T}})sin{\theta }_{\mathrm{defl}}+v_{\mathrm{F}}sin\varphi cos{\theta }_{\mathrm{defl}}-v_{\mathrm{F}}sin\varphi ,\hfill \end{array}$$(32)

and where

$$\begin{array}{c}\hfill {\theta }_{\mathrm{defl}}=\pi -2{cos}^{-1}\left(\frac{1}{\sqrt{1+\frac{b^2{V}^4}{G^2{(M+m)}^2}}}\right)\end{array}$$(33)

is the deflection angle. We now have to multiply $\mathrm{\Delta }{\mathit{p}}_{\mathrm{F}}^{\mu }$ for the differential number of encounters $\text{d}{n}_{\text{enc}}^{\prime }=2\pi n\mathcal{V}{d}^3{\mathit{v}}_{\text{F}}b\text{d}b\text{d}t$ in the laboratory frame. This quantity, at variance with that given by Eq. (18), is a Lorentz-invariant quantity³, where

$$\begin{array}{c}\hfill \mathcal{V}\equiv V(1-{\mathit{v}}_{\mathrm{T}}·{\mathit{v}}_{\mathrm{F}}/c^2)=V(1-v_{\mathrm{T}}{v}_{\mathrm{F}}cos\varphi /c^2).\end{array}$$(34)

In integral form, the momentum variation of the field particle m is now given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{p}}_{\mathrm{F}}^{\mu }}{\mathrm{d}t}=2\pi n\int \int b\mathrm{\Delta }{\mathit{p}}_{\mathrm{F}}^{\mu }\mathcal{V}f({\mathit{v}}_{\mathrm{F}})\mathrm{d}b{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}}.\end{array}$$(35)

The momentum of the test particle M, $p_{\mathrm{T}}^{\mu }=M{\gamma }_{v_{\mathrm{T}}}(c,{\mathit{v}}_{\mathrm{T}})$, will undergo the opposite change

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{p}}_{\mathrm{T}}^{\mu }}{\mathrm{d}t}=-\frac{\mathrm{d}{\mathit{p}}_{\mathrm{F}}^{\mu }}{\mathrm{d}t}.\end{array}$$(36)

To evaluate Eq. (35) we need first to substitute the expression for the cosine and sine of the deflection angle Eq. (33) in Eqs. (30)–(32) that read

$$\begin{array}{c}\hfill cos{\theta }_{\mathrm{defl}}=\frac{\frac{b^2{V}^4}{G^2{(M+m)}^2}-1}{\frac{b^2{V}^4}{G^2{(M+m)}^2}+1};sin{\theta }_{\mathrm{defl}}=\frac{\frac{2b{V}^2}{G(M+m)}}{1+\frac{b^2{V}^4}{G^2{(M+m)}^2}}.\end{array}$$(37)

We note that for the derivation of the DF formula, only the parallel component of v_T contributes, so it is sufficient to evaluate the following expression for the parallel component p_F∥:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{p}_{\mathrm{F}\Vert }}{\mathrm{d}t}& =2\pi nm\int \int b{\gamma }_{v_{\mathrm{F}}}\mathcal{V}f({\mathit{v}}_{\mathrm{F}})\{v_{\mathrm{T}}{\gamma }_{v_{\mathrm{T}}}^2(1-\frac{v_{\mathrm{T}}{v}_{\mathrm{F}}}{c^2}cos\varphi )\hfill \\ \hfill & +{\gamma }_{v_{\mathrm{T}}}[{\gamma }_{v_{\mathrm{T}}}(v_{\mathrm{F}}cos\varphi -v_{\mathrm{T}})\frac{\frac{b^2{V}^4}{G^2{(M+m)}^2}-1}{\frac{b^2{V}^4}{G^2{(M+m)}^2}+1}-v_{\mathrm{F}}sin\varphi \hfill \\ \hfill & \times \frac{\frac{2b{V}^2}{G(M+m)}}{\frac{b^2{V}^4}{G^2{(M+m)}^2}+1}]-v_{\mathrm{F}}cos\varphi \}\mathrm{d}b{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}}.\hfill \end{array}$$(38)

Following the classical derivation discussed in Sect. 3, we perform first the integral over the impact parameter b. As we are considering an isotropic f(v_F), we assume ϕ such that all odd terms involving sin ϕ zero-out. We then apply Eq. (36) obtaining the deceleration on the momentum of particle M, and then we obtain the DF formula for the velocity v_T dividing by γ_vTM as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}=-4\pi G^2\frac{{(M+m)}^2}{M}\rho {\gamma }_{v_{\mathrm{T}}}log\mathrm{\Lambda }\int \frac{{\gamma }_{v_{\mathrm{F}}}\mathcal{V}f({\mathit{v}}_{\mathrm{F}})({\mathit{v}}_{\mathrm{T}}-{\mathit{v}}_{\mathrm{F}})}{V^4}{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}},\\ \hfill \end{array}$$(39)

where again we made use of the limit log(1 + Λ²)→logΛ for large values of Λ and substituted ρ for mn.

As we are accounting for relativistic velocities and transformations, when performing the velocity integral in Eq. (39), with V given by Eq. (27), we need to use a velocity distribution function in covariant form. For example, when considering a thermalized relativistic gas, a natural choice is the Maxwell-Jüttner distribution (see Jüttner 1911; see also Fackerell 1968):

$$\begin{array}{c}\hfill f(v_{\mathrm{F}})=\frac{{\gamma }_{v_{\mathrm{F}}}^5{v}_{\mathrm{F}}^2}{c^3\mathrm{\Theta }{\mathcal{K}}_2({\mathrm{\Theta }}^{-1})}exp(-\frac{{\gamma }_{v_{\mathrm{F}}}}{\mathrm{\Theta }}).\end{array}$$(40)

In the expression above Θ = σ²/3c² and 𝒦₂(x) is the modified Bessel function of the second kind (see e.g. Arfken et al. 2012), often somewhat improperly dubbed the Neumann function.

In the limit of c → +∞, Eq. (39) becomes Eq. (22) with (M + m)²/M in lieu of (M + m). This discrepancy with the classical case occurs because in the relativistic treatment we used the total relativistic momentum conservation of Eq. (36), whereas classically Eqs. (7)–(9) are used. The application of the latter would in principle commute with the limit c → +∞ (and the term (M + m) at the denominator in Eq. (9) would elide with another identical factor in Eq. (39), while Eq. (36) loses meaning in such a limit, as the time component would diverge.

We note that Kalnajs (1972) also encounters a similar disagreement between the Chandrasekhar expression and the particle limit of his fluid formalism (i.e. the frictional force on the test mass is proportional to M²m rather than Mm(M + m)). However, the reason for this difference is ascribed to the fluid picture where the force on M is due to a continuous mass density distribution (see e.g. Kandrup 1983). We also note that, in the fully kinetic approach on DF based on the Fluctuation-Dissipation theorem pioneered by Bekenstein & Maoz (1992), among others, one recovers the original expression formulated by Chandrasekhar with the term Mm(M + m).

Having established this, we estimated the relativistic DF for three cases of astrophysical interest: a 10 M_⊙ BH in a star cluster with velocity dispersion σ ∼ 10 km s⁻¹ and core density of 10² M_⊙ pc⁻³; a 10⁷ M_⊙ BH in a galactic core with σ ∼ 500 km s⁻¹ and density of 4 × 10⁶ M_⊙ pc⁻³; and, finally, stellar mass BH in a relativistic dark matter cusp around a SMBH with σ ∼ 2 × 10⁴ km s⁻¹ and mean density of 7 × 10³ g cm⁻¹. We always assumed an isotropic relativistic velocity distribution as given in Eq. (40) when solving (numerically) the integral in Eq. (39). In Fig. 1 (red solid lines) we compare it with the classical expression (Eq. (1), black dashed lines). In agreement with Syer (1994) we observe for all systems, even if v < σ, that the relativistic DF is augmented with respect to its classical counterpart, up to a factor 1.1, which is in general much smaller than the 16/3γ_vT found by Syer. For v ≳ σ the relativistic DF force is slightly smaller than the classical force, while at large v_T increases again due to the diverging pre-factor γ_vT in the limit v_T → c.

Fig. 1. Relativistic (Eq. (39) red solid lines, and Eq. (43) blue dot-dashed lines) and classical (dashed lines) dynamical friction force as a function of the test particle velocity vT for three different cases with isotropic relativistic distribution. Shown are a 10 M⊙ BH in a star cluster with velocity dispersion σ ∼ 10 km s−1 and core density of 102 M⊙ pc−3 (left panel); a 107 M⊙ BH in a galactic core with σ ∼ 500 km s−1 and density of 4 × 106 M⊙ pc−3 (mid panel); and a stellar mass BH in a relativistic dark matter cusp around a SMBH with σ ∼ 2 × 104 km s−1 and mean density of 7 × 10t g cm−1 (right panel). In all cases the relativstic and classical curves are normalized to the value of the main peak of the classical.

We note that, for relativistic power-law velocity distributions, Eq. (39) would become for a given v_T considerably larger than (1), due to the contribution of large v tails. We speculate that this could be relevant in the context of plasma physics where (multiple) power-law velocity distributions are often encountered (see e.g. Sandquist et al. 2006). The derivation of Eq. (39) can be carried out in a similar fashion for a charge q_T deflected in a plasma, as the impact parameter and velocity integrals are the same, the only difference being the dimensional factor containing the masses and the gravitational constant. We also note that Eq. (39) was obtained under the assumption that the classical angle ψ₀ given in Eq. (14) holds even in the case of relativistic velocities. To be more rigorous, one should use instead its relativistic generalization given by

$$\begin{array}{c}\hfill tan{\stackrel{\sim }{\psi }}_0=-\frac{L\sqrt{2E}}{G(M+m)},\end{array}$$(41)

where L = γ_VVb and E = c²(γ_V − 1) are the norm of the spatial part of the specific relativistic angular momentum and the specific relativistic kinetic energy, respectively, and γ_V is the Lorentz factor of the relative velocity V.

With this choice, tanψ₀ and $tan{\stackrel{\sim }{\psi }}_0$ differ by the multiplicative factor ${\gamma }_V\sqrt{2({\gamma }_V-1)}c/V$. The latter increases significantly the value of tanψ₀ only for V ≫ 0.5c. In this limit, expanding the square roots arising from γ_V, the terms containing b²V⁴/G²(M + m)² in Eqs. (33)–(37) in the derivation above are augmented by a factor (1 + V²/2c²)², so that Eq. (38) becomes

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{p}_{\mathrm{F}\Vert }}{\mathrm{d}t}& =2\pi nm\int \int b{\gamma }_{v_{\mathrm{F}}}\mathcal{V}f({\mathit{v}}_{\mathrm{F}})\{v_{\mathrm{T}}{\gamma }_{v_{\mathrm{T}}}^2(1-\frac{v_{\mathrm{T}}{v}_{\mathrm{F}}}{c^2}cos\varphi )\hfill \\ \hfill & +{\gamma }_{v_{\mathrm{T}}}[{\gamma }_{v_{\mathrm{T}}}(v_{\mathrm{F}}cos\varphi -v_{\mathrm{T}})\frac{\frac{b^2{V}^4}{G^2{(M+m)}^2}-(1+\frac{V^2}{2c^2}{)}^{-2}}{\frac{b^2{V}^4}{G^2{(M+m)}^2}+(1+\frac{V^2}{2c^2}{)}^{-2}}-v_{\mathrm{F}}sin\varphi \hfill \\ \hfill & \times \frac{\frac{2b{V}^2}{G(M+m)}(1+\frac{V^2}{2c^2})}{\frac{b^2{V}^4}{G^2{(M+m)}^2}(1+\frac{V^2}{2c^2}{)}^2+1}]-v_{\mathrm{F}}cos\varphi \}\mathrm{d}b{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}}.\hfill \end{array}$$(42)

The integration over the impact parameter b can be carried out first in the same way as above, yielding a velocity averaged Coulomb logarithm now multiplied by the pre-factor (1 + V²/2c²)⁻², so that Eq. (39) becomes

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}=-4\pi G^2\frac{{(M+m)}^2}{M}\rho {\gamma }_{v_{\mathrm{T}}}log\mathrm{\Lambda }\int \frac{{\gamma }_{v_{\mathrm{F}}}\mathcal{V}f({\mathit{v}}_{\mathrm{F}})({\mathit{v}}_{\mathrm{T}}-{\mathit{v}}_{\mathrm{F}})}{(1+\frac{V^2}{c^2}{)}^2{V}^4}{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}}.\end{array}$$(43)

Equation (43) differs significantly from the simplified expression Eq. (39) only in the limit of large velocity dispersion (and for high velocities), as shown by the blue dot-dashed lines in the right panel of Fig. 1.

5. Post-Newtonian approximation

So far, we have derived a formal generalization of the dynamical friction formula in the limit of high test particle velocities or relativistic velocity distributions assuming classical Newtonian forces. We now carry out an alternative derivation involving strong gravitational scattering in the post-Newtonian regime. Introduced by Einstein (1915; see also Weinberg 1972; Blanchet 2010 and references therein) to study the precession of the perihelion of Mercury, the post-Newtonian approximation consists of an expansion in orders of the parameter v/c, such that at the zeroth order it reduces to Newtonian gravity, while at higher orders (nPN) the acceleration on the mass m due to the mass M is augmented by corrections of order (v/c)²ⁿ ∼ (GM/rc²)ⁿ.

5.1. Non-relativistic velocities

The cores of dense star clusters are often dominated by massive objects due to dynamical mass segregation; it is therefore interesting to evaluate the DF on a test particle in such an environment where strong scattering by larger masses are likely to happen, even though the velocity distribution might not be relativistic. To do so, we begin with a naive derivation of dv_T∥/dt in impulsive approximation keeping the Galilean transformations of velocities, but using the 1PN acceleration

$$\begin{array}{cc}\hfill {\mathit{a}}_{1\mathrm{PN}}& =-\frac{G(M+m)}{r^2}\frac{\mathit{r}}{r}+\frac{G(M+m)}{c^2{r}^2}\{[(4+2\eta )\frac{G(M+m)}{r}+\hfill \\ \hfill & -(1+3\eta )V^2+\frac{3}{2}\eta {\dot{\mathit{r}}}^2]\frac{\mathit{r}}{r}+(4-2\eta )\dot{\mathit{r}}\mathit{V}\}\hfill \end{array}$$(44)

in the frame centred on the field particle m (see Mora & Will 2004; Cashen et al. 2017). In the equation above, the parameters μ, V, and r have the same meaning as in Sect. 2, and η = μ/(M + m). In impulsive approximation we can express the velocity change of the test particles in a discrete time interval Δt = 2b/V as ||ΔV_T|| ∼ aΔt = 2ab/V (see e.g. Ciotti 2021), so that its parallel component becomes

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{v}}_{\mathrm{T}\Vert }\sim -\frac{\mu }{M}\frac{||\mathrm{\Delta }{V}_{\perp }{||}^2}{2V^2}\mathit{V}.\end{array}$$(45)

Defining $\mathit{r}=r\widehat{\mathit{r}}\sim b\widehat{\mathit{r}}$ and $\mathit{V}=V\widehat{\mathit{V}}$, we obtain the perpendicular relative velocity change as

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\mathit{V}}_{\perp }& \sim {\mathit{a}}_{1\mathrm{PN}}\frac{2b}{V}\sim -\frac{2G(M+m)}{\mathit{Vb}}\widehat{\mathit{r}}+\frac{2G(M+m)}{c^{2}Vb}\hfill \\ \hfill & \times \{[(4+2\eta )\frac{G(M+m)}{b}+-(1+3\eta )V^2+\frac{3}{2}\eta V^2]\widehat{\mathit{r}}\hfill \\ \hfill & +(4-2\eta )V^2\widehat{\mathit{V}}\},\hfill \end{array}$$(46)

and its square as

$$\begin{array}{cc}& ||\mathrm{\Delta }{\mathit{V}}_{\perp }{||}^2\sim \frac{4G^2{(M+m)}^2}{V^2{b}^2}-\frac{8G^2{(M+m)}^2}{c^2{V}^2{b}^2}\hfill \\ \hfill & \times [(4+2\eta )\frac{G(M+m)}{b}+3V^2-\frac{7}{2}\eta V^2],\hfill \end{array}$$(47)

where we assume that $\widehat{\mathit{r}}·\widehat{\mathit{V}}\sim 1$ and drop the terms proportional to 1/c⁴. The parallel velocity change of the test particle becomes

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\mathit{v}}_{\mathrm{T}\Vert }& \sim -\frac{\mu }{M}\frac{||\mathrm{\Delta }{\mathit{V}}_{\perp }{||}^2}{2V^2}\mathit{V}\hfill \\ \hfill & \sim [-\frac{2m{G}^2(M+m)}{b^2{V}^4}+\frac{4G^{3}m{(M+m)}^2(4+2\eta )}{c^2{b}^3{V}^4}\hfill \\ \hfill & +\frac{4G^{2}m(M+m)(3-\frac{7}{2}\eta )}{c^2{b}^2{V}^2}]\mathit{V},\hfill \end{array}$$(48)

so that, assuming again as in Sect. 2 that the number of encounters is given by Eq. (18), the finite differences velocity change of particle M is expressed as

$$\begin{array}{cc}\hfill \frac{\mathrm{\Delta }{\mathit{v}}_{\mathrm{T}}}{\mathrm{\Delta }t}& =2\pi b\mathrm{d}bVnf({\mathit{v}}_{\mathrm{F}})d^3{\mathit{v}}_{\mathrm{F}}[-\frac{2m{G}^2(M+m)}{b^2{V}^4}\mathit{V}\hfill \\ \hfill & +\frac{4G^{3}m{(M+m)}^2(4+2\eta )}{c^2{b}^3{V}^4}\mathit{V}+\frac{4G^{2}m(M+m)(3-\frac{7}{2}\eta )}{c^2{b}^2{V}^2}\mathit{V}].\hfill \end{array}$$(49)

With the standard integration over the impact parameter b, we easily obtain

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}& =-4\pi mn{G}^2(M+m)log\mathrm{\Lambda }\int \frac{f({\mathit{v}}_{\mathrm{F}})\mathit{V}}{V^3}{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}}\hfill \\ \hfill & +\frac{16\pi mn{G}^3{(M+m)}^2(2+\eta )}{{\stackrel{\sim }{b}}_{\min }{c}^2}\int \frac{f({\mathit{v}}_{\mathrm{F}})\mathit{V}}{V^3}{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}}\hfill \\ \hfill & +\frac{8\pi mn{G}^2(M+m)(3-\frac{7}{2}\eta )log\mathrm{\Lambda }}{c^2}\int \frac{f({\mathit{v}}_{\mathrm{F}})\mathit{V}}{V}{\mathrm{d}}^3{\mathit{v}}_{\mathrm{F}},\hfill \end{array}$$(50)

where the term (b_max − b_min)/b_minb_max, arising from the integration in b, is substituted in the second addendum with ${\stackrel{\sim }{b}}_{\min }^{-1}$. The latter is a bona fide minimum impact parameter defined by ${\stackrel{\sim }{b}}_{\min }=2G(M+m)/{\sigma }^2$, where σ is the velocity dispersion of f.

In practice, a naive 1PN extension of the classical case, independently of the specific choice of f(v_F) augments the Chandrasekhar expression of the DF (first line of Eq. (50)) of two additional terms. The first has a different dependence on the impact parameter but the same integral on V, while the second contains the classical Coulomb logarithm but a different integral on V. For the special case of an isotropic Maxwellian distribution with velocity dispersion σ, the 1PN DF expression can be easily integrated in the limit of v_T < σ and rewritten as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathit{v}}_{\mathrm{T}}}{\mathrm{d}t}=-4\pi \rho G^2(M+m)\frac{log\mathrm{\Lambda }}{{\sigma }^3}\zeta {\mathit{v}}_{\mathrm{T}},\end{array}$$(51)

where the corrective factor ζ, plotted in Fig. 2 as a function of σ/c is defined by

$$\begin{array}{c}\hfill \zeta =1-[\frac{8+4\eta }{log\mathrm{\Lambda }}+6-14\eta ]\frac{{\sigma }^2}{c^2}.\end{array}$$(52)

Fig. 2. Correcting factor ζ in the low-velocity limit as a function of σ.

Interestingly, the net effect of a simple 1PN correction is to reduce the DF drag force with respect to its classical counterpart given by Eq. (1). This has the relevant consequence that a highly massive object, for which the GR corrections to its gravitational field cannot be neglected (e.g. a SMBH moving at non-relativistic speed in a star system), undergos a less effective gravitational drag than what would be estimated using the Chandrasekhar formula. This implies in that specific case an even longer in-spiral timescale for the BHs in galactic cores. The reason for a lower DF effect in a simple 1PN approximation is ascribed to the fact that the relativistic precession due to the velocity dependent force term in a two-body scatter always acts in the opposite direction with respect to the deflection caused by the 1/r² force term.

The behaviour of the 1PN DF cannot be extrapolated for high (relativistic) velocities since for such high values of v the velocity composition should be made in terms of Lorentz transforms, as in Sect. 4.

5.2. Relativistic velocities

We now extend the DF formula to the case where the particle velocities are relativistic and the density is such that the force is to be evaluated at 1PN order during close encounters. Following Lee (1969), we consider a mass M moving at v_T in a uniform medium of field stars, of constant number density n, with isotropic velocity distribution f(v_F), in an inertial frame S. As usual, in a time interval Δt, its velocity changes by an amount ∑Δv_T after n_enc encounters.

As the effect of General Relativity is accounted for only during the deflection of the test particle M, its velocity is transformed according to the Lorentz transformations of Special Relativity.

As in the classical case, the isotropy of f(v_F) allows us to write

$$\begin{array}{c}\hfill \sum \mathrm{\Delta }{\mathit{v}}_{\mathrm{T}}\wedge {\mathit{v}}_{\mathrm{T}}=0,\end{array}$$(53)

so it is sufficient to perform

$$\begin{array}{c}\hfill \sum \mathrm{\Delta }{v}_{\mathrm{T}\Vert }=\frac{\sum \mathrm{\Delta }{\mathit{v}}_{\mathrm{T}}·{\mathit{v}}_{\mathrm{T}}}{v_{\mathrm{T}}}.\end{array}$$(54)

Let us consider a single encounter of a test star with a field star: the velocity v_T becomes ${\mathit{v}}_{\mathrm{T}}^f$, so we have

$$\begin{array}{c}\hfill \mathrm{\Delta }{v}_{\mathrm{T}\Vert }=\frac{({\mathit{v}}_{\mathrm{T}}^f-{\mathit{v}}_{\mathrm{T}})·{\mathit{v}}_{\mathrm{T}}}{v_{\mathrm{T}}}.\end{array}$$(55)

It is now useful to introduce a second reference frame, S′, in which the test star is, initially, at rest. Let ${\mathit{w}}_{\mathrm{T}}^f$ be the velocity of the test star after the encounter, in the frame S′. With our assumptions (the motion of all the stars is approximately described by straight line) we can write, without losing generality, that

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{T}}^f=\frac{({\mathit{w}}_{\mathrm{T}}^f+{\mathit{v}}_{\mathrm{T}})}{1+\frac{{\mathit{w}}_{\mathrm{T}}^f·{\mathit{v}}_{\mathrm{T}}}{c^2}},\end{array}$$(56)

and therefore

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{T}}^f·{\mathit{v}}_{\mathrm{T}}=\frac{({\mathit{w}}_{\mathrm{T}}^f+{\mathit{v}}_{\mathrm{T}})·{\mathit{v}}_{\mathrm{T}}}{1+\frac{{\mathit{w}}_{\mathrm{T}}^f·{\mathit{v}}_{\mathrm{T}}}{c^2}}.\end{array}$$(57)

Noting that $\mathrm{\Delta }{\mathit{w}}_{\mathrm{T}}={\mathit{w}}_{\mathrm{T}}^f$, (i.e. the test star is initially at rest in S′) and keeping only terms of order 1/c², after some algebra Eq. (55) becomes

$$\begin{array}{cc}\hfill \mathrm{\Delta }{v}_{\mathrm{T}\Vert }& =\frac{1}{v_{\mathrm{T}}}[\frac{({\mathit{w}}_{\mathrm{T}}^f+{\mathit{v}}_{\mathrm{T}})·{\mathit{v}}_{\mathrm{T}}}{1+\frac{{\mathit{w}}_{\mathrm{T}}^f·{\mathit{v}}_{\mathrm{T}}}{c^2}}-v_{\mathrm{T}}^2]\underset{{\mathit{w}}_{\mathrm{T}}^f=\mathrm{\Delta }{\mathit{w}}_{\mathrm{T}}}{\approx }\hfill \\ \hfill & \approx \frac{1}{v_{\mathrm{T}}}[\mathrm{\Delta }{\mathit{w}}_{\mathrm{T}}·{\mathit{v}}_{\mathrm{T}}{\gamma }_{v_T}^{-1}-\frac{{(\mathrm{\Delta }{\mathit{w}}_{\mathrm{T}}·{\mathit{v}}_{\mathrm{T}})}^2}{c^2}].\hfill \end{array}$$(58)

We must now sum Eq. (58) over all possible values of Δw_T. We recall that v_F and w_F are the initial velocities of field stars in the frames S and S′, related by the following Lorentz transformation (see e.g. Landau & Lifshitz 1976)

$$\begin{array}{c}\hfill {\mathit{w}}_{\mathrm{F}}=\frac{1}{1-\frac{{\mathit{v}}_{\mathrm{F}}·{\mathit{v}}_{\tau }}{c^2}}[\frac{{\mathit{v}}_{\mathrm{F}}}{{\gamma }_{\mathrm{v}\tau }}-{\mathit{v}}_{\tau }+\frac{1}{c^2}\frac{{\gamma }_{\mathrm{v}\tau }({\mathit{v}}_{\mathrm{F}}·{\mathit{v}}_{\tau }){\mathit{v}}_{\tau }}{{\gamma }_{\mathrm{v}\tau }+1}],\end{array}$$(59)

where v_τ is the translational velocity of S′ relative to S; ${\gamma }_{\mathrm{v}\tau }={(1-{\mathit{v}}_{\tau }^2/c^2)}^{-1/2}$; and b is, as usual, the impact parameter of the encounter, defined in S′. We also have that w_F ⋅ b = w_F b cos φ, while Δw_T, decomposed into its components parallel and perpendicular to w_F, becomes

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{w}}_{\mathrm{T}}=\mathrm{\Delta }{w}_{\mathrm{T}\Vert }\frac{{\mathit{w}}_{\mathrm{F}}}{w_{\mathrm{F}}}+\mathrm{\Delta }{w}_{\mathrm{T}\perp }\frac{\mathit{b}}{b}.\end{array}$$(60)

The velocity distribution of field particles in S′ is (see e.g. Landau & Lifshitz 1976)

$$\begin{array}{c}\hfill f^{\prime }({\mathit{w}}_{\mathrm{F}})={\gamma }_{v_{\mathrm{T}}}f({\mathit{v}}_{\mathrm{F}})(\frac{{\gamma }_{w_{\mathrm{F}}}}{{\gamma }_{v_{\mathrm{F}}}}{)}^4,\end{array}$$(61)

where v_F = v_F(w_F, v_T). Let us denote with

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}^{\prime }={\gamma }_{v_T}^{-1}\mathrm{\Delta }t\end{array}$$(62)

the transformed time interval during which the encounters with the field particles are summed.

It is important to note that $f\prime ({\mathit{w}}_{\mathrm{F}})$ in the frame S′ is only homogeneous, but no longer isotropic. As a consequence of this, the impact parameter b and the deflection angle ϕ are randomly distributed for a given value of w_F. Integrating Eq. (58) over all values of Δw_T, we obtain the formal result

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{v}_{\mathrm{T}\Vert }}{\mathrm{d}t}& =\int n{f}^{\prime }({\mathit{w}}_{\mathrm{F}}){\gamma }_{v_T}^{-1}{w}_{\mathrm{F}}{d}^3{\mathit{w}}_{\mathrm{F}}{\int }_{b_{\min }}^{b_{\max }}2\pi b\mathrm{d}b\hfill \\ \hfill & \times {\int }_0^{2\pi }\frac{1}{2\pi }\frac{1}{v_{\mathrm{T}}}[\mathrm{\Delta }{\mathit{w}}_{\mathrm{T}}·{\mathit{v}}_{\mathrm{T}}{\gamma }_{v_T}^{-2}-\frac{{(\mathrm{\Delta }{\mathit{w}}_{\mathrm{T}}·{\mathit{v}}_{\mathrm{T}})}^2}{c^2}]\mathrm{d}\phi .\hfill \end{array}$$(63)

The quantity Δw_T(w_F, b) appearing in the equation above must be expressed as a function of w_F and b, and can be computed with the help of the two-body problem 1PN-Lagrangian. For this purpose, it is convenient to introduce a third reference frame S″, in which the centre of mass (c.o.m.) of the encounter is at rest at the origin of coordinates⁴, where we denote the particle velocities as u.

At 1PN order the equations of motion for the two-body encounter are usually derived from the Einstein-Infeld-Hoffmann Lagrangian (EIH, see Einstein et al. 1938, see also Eddington & Clark 1938)

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{EIH}}& =\frac{1}{2}m{\mathit{u}}_{\mathrm{F}}^2+\frac{1}{2}M{\mathit{u}}_{\mathrm{T}}^2+\frac{1}{8c^2}(m{\mathit{u}}_{\mathrm{F}}^4+M{\mathit{u}}_{\mathrm{T}}^4)+\frac{\mathit{GmM}}{r}\hfill \\ \hfill & \times \{1-\frac{1}{2c^2}[7({\mathit{u}}_{\mathrm{F}}·{\mathit{u}}_{\mathrm{T}})+({\mathit{u}}_{\mathrm{F}}·\mathit{n})({\mathit{u}}_{\mathrm{T}}·\mathit{n})-3{({\mathit{u}}_{\mathrm{F}}-{\mathit{u}}_{\mathrm{T}})}^2]\hfill \\ \hfill & -\frac{G(m+M)}{2r{c}^2}\}.\hfill \end{array}$$(64)

Remarkably, qualitatively similar (but slightly simpler) equations of motion can be derived using the gravitational analogue of the Darwin Lagrangian of electrodynamics (see e.g. Jackson 1975; see also Essén 2007, 2014), often referred to as the Fock Lagrangian (see Fock 1964; Kennedy 1972; Deruelle & Uzan 2018):

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{Darwin}}& =\frac{1}{2}m{\mathit{u}}_{\mathrm{F}}^2+\frac{1}{2}M{\mathit{u}}_{\mathrm{T}}^2+\frac{1}{8c^2}(m{\mathit{u}}_{\mathrm{F}}^4+M{\mathit{u}}_{\mathrm{T}}^4)+\frac{\mathit{GmM}}{r}\hfill \\ \hfill & \times [1-\frac{1}{2c^2}({\mathit{u}}_{\mathrm{F}}·{\mathit{u}}_{\mathrm{T}}+({\mathit{u}}_{\mathrm{F}}·\mathit{n})({\mathit{u}}_{\mathrm{T}}·\mathit{n}))].\hfill \end{array}$$(65)

In Eqs. (64) and (65), r = r_F − r_T is the (instantaneous) relative position vector and n = r/||r|| (see Fagundes et al. 1976; Damour & Deruelle 1985; Zürcher 2017, for details). In both cases the c.o.m. coordinates at the 1PN order are

$$\begin{array}{c}\hfill {\mathit{r}}_{\mathrm{cm}}=\frac{{\mathcal{E}}_{\mathrm{T}}{\mathit{r}}_{\mathrm{T}}+{\mathcal{E}}_{\mathrm{F}}{\mathit{r}}_{\mathrm{F}}}{{\mathcal{E}}_{\mathrm{T}}+{\mathcal{E}}_{\mathrm{F}}},\end{array}$$(66)

where

$$\begin{array}{cc}\hfill {\mathcal{E}}_{\mathrm{T}}& =M{c}^2+\frac{1}{2}M{\mathit{u}}_{\mathrm{T}}^2-\frac{1}{2}\frac{\mathit{GMm}}{||{\mathit{r}}_{\mathrm{T}}-{\mathit{r}}_{\mathrm{F}}||},\hfill \\ \hfill {\mathcal{E}}_{\mathrm{F}}& =m{c}^2+\frac{1}{2}m{\mathit{u}}_{\mathrm{F}}^2-\frac{1}{2}\frac{\mathit{GMm}}{||{\mathit{r}}_{\mathrm{T}}-{\mathit{r}}_{\mathrm{F}}||}.\hfill \end{array}$$(67)

Since the c.o.m is by definition at rest at the origin of S″ (i.e. r_cm = 0), in terms of the relative velocity u we can always write

$$\begin{array}{c}\hfill {\mathit{u}}_{\mathrm{T}}=\frac{m}{m+M}\mathit{u}+\mathcal{O}(1/c^2);{\mathit{u}}_{\mathrm{F}}=-\frac{M}{m+M}\mathit{u}+\mathcal{O}(1/c^2),\end{array}$$(68)

so that

$$\begin{array}{c}\hfill \frac{1}{2}(M{u}_{\mathrm{T}}^2+m{u}_{\mathrm{F}}^2)=\frac{1}{2}(\frac{\mathit{mM}}{M+m}{u}^2)+\mathcal{O}(1/c^4).\end{array}$$(69)

At variance with Lee (1969), we will assume hereafter the Lagrangian (65), that in terms of the relative velocity u, once conveniently transformed in polar coordinates, becomes

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{D}}& =\frac{1}{2}\mu [(\frac{\mathrm{d}r}{\mathrm{d}t}{)}^2+r^2(\frac{\mathrm{d}\vartheta }{\mathrm{d}t}{)}^2]+\mathrm{\Phi }+\frac{1}{8c^2}\mu (\frac{\mu }{{\mu }_3}{)}^3\hfill \\ \hfill & \times [(\frac{\mathrm{d}r}{\mathrm{d}t}{)}^2+r^2(\frac{\mathrm{d}\vartheta }{\mathrm{d}t}{)}^2{]}^2+\frac{\mathrm{\Phi }}{2c^2}\frac{\mu }{\mathcal{M}}[2(\frac{\mathrm{d}r}{\mathrm{d}t}{)}^2+r^2(\frac{\mathrm{d}\vartheta }{\mathrm{d}t}{)}^2].\hfill \end{array}$$(70)

In the equation above, μ is again the reduced mass, ℳ = M + m, μ₃ = [m³M³/(m³ + M³)]^1/3, and Φ = Gμℳ/r is the Newtonian gravitational potential energy. With this choice, the equations of motion of the associated particles are obtained in explicit form at the order 1/c², while conversely, in the EIH case the force of particle M on particle m depends on the force of particle m on particle M.

In analogy with the non-relativistic velocity case discussed above, we now obtain (see Appendix A for the explicit mathematical details) the DF formula in integral form as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{v}_{\mathrm{T}\Vert }}{\mathrm{d}t}& =\int n{f}^{\prime }({\mathit{w}}_{\mathrm{F}}){\gamma }_{v_{\mathrm{T}}}^{-1}{w}_{\mathrm{F}}{\mathrm{d}}^3{\mathit{w}}_{\mathrm{F}}{\int }_{b_{\min }}^{b_{\max }}2\pi b\mathrm{d}b\hfill \\ \hfill & {\int }_0^{2\pi }\frac{1}{2\pi }\frac{1}{v_{\mathrm{T}}}\{{\gamma }_{v_{\mathrm{T}}}^{-2}[\frac{{\mathit{v}}_{\mathrm{T}}·{\mathit{w}}_{\mathrm{F}}}{w_{\mathrm{F}}}\mathrm{\Delta }{w}_{\mathrm{T}\Vert }+\frac{{\mathit{v}}_{\mathrm{T}}·\mathit{b}}{b}\mathrm{\Delta }{w}_{\mathrm{T}\perp }]+\hfill \\ \hfill & -\frac{1}{c^2}[(\frac{{\mathit{v}}_{\mathrm{T}}·{\mathit{w}}_{\mathrm{F}}}{w_{\mathrm{F}}}{)}^2(\mathrm{\Delta }{w}_{\mathrm{T}\Vert }{)}^2-2(\frac{{\mathit{v}}_{\mathrm{T}}·{\mathit{w}}_{\mathrm{F}}}{w_{\mathrm{F}}})(\frac{{\mathit{v}}_{\mathrm{T}}·\mathit{b}}{b})\mathrm{\Delta }{w}_{\mathrm{T}\Vert }\mathrm{\Delta }{w}_{\mathrm{T}\perp }\hfill \\ \hfill & +(\frac{{\mathit{v}}_{\mathrm{T}}·\mathit{b}}{b}{)}^2(\mathrm{\Delta }{w}_{\mathrm{T}\perp }{)}^2]\}\mathrm{d}\phi .\hfill \end{array}$$(71)