# Acoustic wake in an isothermal profile: dynamical friction and gravitational
wave emission

Gali Eytan<EMAIL_ADDRESS>Physics Department, Technion –
Israel Institute of Technology, Haifa 3200003, Israel  Vincent Desjacques
<EMAIL_ADDRESS>Physics Department, Technion – Israel Institute
of Technology, Haifa 3200003, Israel  Robin Buehler
<EMAIL_ADDRESS>Physics Department, Technion – Israel
Institute of Technology, Haifa 3200003, Israel

###### Abstract

We consider the motion of a circularly-moving perturber in a self-gravitating,
collisional system with spherically symmetric density profile. We concentrate
on the singular isothermal sphere which, despite its pathological features,
admits a simple polarization function in linear response theory. This allows
us to solve for the acoustic wake trailing the perturber and the resulting
dynamical friction, in the limit where the self-gravity of the response can be
ignored. In steady-state and for subsonic velocities $v_{p}<c_{s}$, the
dynamical friction torque $F_{\varphi}\propto v_{p}^{3}$ is suppressed for
perturbers orbiting in an isothermal sphere relative to the infinite,
homogeneous medium expectation $F_{\varphi}\propto v_{p}$. For highly
supersonic motions, both expectations agree and are consistent with a local
approximation to the gravitational torque. At fixed resolution (a given
Coulomb logarithm), the response of the system is maximal for Mach numbers
near the constant circular velocity of the singular isothermal profile. This
resonance maximizes the gravitational wave (GW) emission produced by the
trailing acoustic wake. For an inspiral around a massive black hole of mass
$10^{6}M_{\odot}$ located at the center of a (truncated) isothermal sphere,
this GW signal could be comparable to the vacuum GW emission of the black hole
binary at sub-nanohertz frequencies when the small black hole enters the Bondi
sphere of the massive one. The exact magnitude of this effect depends on
departures from hydrostatic equilibrium and on the viscosity present in any
realistic astrophysical fluid, which are not included in our simplified
description.

##  1 Introduction

Astrophysical bodies from the size of black holes to galaxies and dark matter
halos merge by transferring orbital energy and angular momentum to their
gravitational interaction with the astrophysical environment. This physical
process is known as Dynamical Friction (DF) since Chandrasekhar (1943), and is
a key ingredient in cosmic structure formation (Tremaine et al., 1975;
Paczynski, 1976; Murai & Fujimoto, 1980; Binney & Tremaine, 1987; Kauffmann et
al., 1993; Ibata & Lewis, 1998; Somerville & Primack, 1999; van den Bosch et
al., 1999; Ostriker, 1999; Cole et al., 2000; Goldreich et al., 2004; Croton
et al., 2006; Boylan-Kolchin et al., 2008; Mo et al., 2010). In a seminal
work, Chandrasekhar (1943) derived an expression for the gravitational drag or
DF force produced by a point-like perturber moving in linear motion in a
infinite, homogeneous and collisionless medium. His calculation has been
extended to different boundary conditions and backgrounds to gain
understanding on the physics of DF. The latter includes stellar distributions
(Lynden-Bell & Kalnajs, 1972; Tremaine & Weinberg, 1984; Palmer & Papaloizou,
1985; Weinberg, 1986; Bekenstein & Maoz, 1992; Nelson & Tremaine, 1999; Banik
& van den Bosch, 2021; Ginat et al., 2023; Kipper et al., 2023), gaseous media
(Dokuchaev, 1964; Ruderman & Spiegel, 1971; Rephaeli & Salpeter, 1980; Just &
Kegel, 1990; Ostriker, 1999; Sánchez-Salcedo & Brandenburg, 2001; Kim & Kim,
2007; Lee & Stahler, 2011; Vicente et al., 2019; Sánchez-Salcedo, 2019;
Desjacques et al., 2022) and, more recently, exotic dark matter backgrounds
(Hui et al., 2017; Bar-Or et al., 2019; Berezhiani et al., 2019; Annulli et
al., 2020; Chavanis, 2021a, b; Traykova et al., 2021; Hartman et al., 2021;
Bar et al., 2021; Vicente & Cardoso, 2022; Buehler & Desjacques, 2023; Foote
et al., 2023; Tomaselli et al., 2023; Traykova et al., 2023; Boudon et al.,
2024). In all cases, the DF force is given by the gravitational deceleration
produced by a nonlocal or global density wake trailing the perturber, as was
already pointed out by Kalnajs (1971, 1972).

In this work, we calculate the global density wake produced by a circularly-
moving perturber in an ideal fluid with a $r^{-2}$ isothermal profile. The
case of a general, realistic fluid with finite mean free path
$\lambda_{\text{mfp}}$ explored in, e.g., Katz et al. (2019), is beyond the
scope of this paper. Nonetheless, our study complements Tremaine & Weinberg
(1984); Weinberg (1986); Kaur & Stone (2022), who focused on collisionless
spherical systems, and provides useful insights on dynamical friction for
self-gravitating spherical systems with vanishing $\lambda_{\text{mfp}}$. We
assume that the self-gravity of the density wake can be neglected and use
linear response theory to compute the polarization function (Kalnajs, 1971;
Tremaine & Weinberg, 1984; Binney & Tremaine, 1987; Weinberg, 1989; Nelson &
Tremaine, 1999). Following Tremaine & Weinberg (1984); Desjacques et al.
(2022) for the collisionless and collisional limit, respectively, we expand
the perturbations in multipoles to solve for the trailing acoustic wake and
the friction coefficient. We concentrate on the singular isothermal profile,
which admits a simple polarization function at the expense of its well-known
pathological features (Chandrasekhar, 1939). Our approach can, in principle,
be applied to any spherically symmetric profile in hydrostatic equilibrium.

The paper is organized as follows. In §2, we summarize the assumptions behind
our approach and spell out the derivation of the polarization function, the
acoustic wake density and the dynamical friction. We validate our approach
with the infinite, homogeneous case for which semi-numerical simulations are
available. In §3, we apply our method to a singular isothermal sphere, and
derive the trailing acoustic wake and the dynamical friction in the steady-
state regime. In §4, we discuss the gravitational wave emission emanating from
a trailing acoustic wake induced by a compact binary embedded in a (truncated)
isothermal sphere. We conclude in §5.

##  2 Theory

Consider a spherically symmetric, inviscid gaseous medium in hydrostatic
equilibrium of density $\rho_{0}(r)$ and (adiabatic) sound speed $c_{0}(r)$.
The equation of state for the pressure $P_{0}(r)$ is the polytrope
$P_{0}=K\rho_{0}^{\gamma}$ where $K$ is a constant and $\gamma$ is the
adiabatic index. The potential of this self-gravitating medium is
$\phi_{0}(r)$.

We are interested in the li near response of this gaseous system to an
external, pointlike perturbation of mass $M$ moving on a fixed circular orbit
of radius $a$ and angular velocity $\Omega>0$, in the $x-y$ plane.

Let $\textbf{\em r}_{p}=\textbf{\em r}_{p}(t)$ and ${\textbf{\em
r}_{p}}^{\prime}=\textbf{\em r}_{p}(t^{\prime})$ be the position of the
perturber at ”present-day” time $t$ and ”retarded” time $t^{\prime}$,
respectively, relative to the center of the density profile. We have
$\textbf{\em r}_{p}(t)=a(\cos\Omega t\,\mathrm{\hat{\bf{x}}}+\sin\Omega
t\,\mathrm{\hat{\bf{y}}})$. The Mach number is $\mathcal{M}\equiv
v_{p}/c_{s}$, where $v_{p}=\Omega a$ is the perturber’s circular velocity. The
point mass has a density $\rho_{p}(\textbf{\em
r},t)=Mh(t)\delta^{D}(\textbf{\em r}-\textbf{\em r}_{p}(t))$, where the time-
dependent function $h(t)$ is unity (zero) when the perturber is turned on
(off).

###  2.1 Methodology

####  2.1.1 Sound wave equation

Adding the perturber to the system leads to linearized continuity and momentum
conservation equations involving perturbations in the gas density $\delta\rho$
and pressure $\delta P=c_{0}^{2}\delta\rho$ as well as the perturber’s
gravitational potential $\phi_{p}$,

$\displaystyle\partial_{t}\delta\rho$
$\displaystyle=-\boldsymbol{\nabla}\cdot\big{(}\rho_{0}\delta\textbf{\em
v}\big{)}$ (1) $\displaystyle\partial_{t}\big{(}\rho_{0}\delta\textbf{\em
v}\big{)}$ $\displaystyle=-\boldsymbol{\nabla}\delta
P-\rho_{0}\boldsymbol{\nabla}\delta\phi-\delta\rho\boldsymbol{\nabla}\phi_{0}-\rho_{0}\boldsymbol{\nabla}\phi_{p}\;.$
(2)

These equations are supplemented by the perturbed Poisson equations

$\nabla^{2}\delta\phi=4\pi
G\delta\rho\quad\mbox{and}\quad\nabla^{2}\delta\phi_{P}=4\pi
G\delta\rho_{P}\;,$ (3)

where $\delta\phi$ is the perturbation to the potential of the gaseous sphere.
Using the hydrostatic equilibrium condition
$(1/\rho_{0})\boldsymbol{\nabla}P_{0}+\boldsymbol{\nabla}\phi_{0}=0$, they can
be combined to yield a perturbation equation of the form

$\ddot{\delta}\rho-\nabla^{2}\delta
P-\boldsymbol{\nabla}\phi_{0}\cdot\boldsymbol{\nabla}\delta\rho-4\pi
G\rho_{0}\delta\rho\\\
=\boldsymbol{\nabla}\big{(}\rho_{0}\boldsymbol{\nabla}\phi_{p}\big{)}+\boldsymbol{\nabla}\big{(}\rho_{0}\boldsymbol{\nabla}\delta\phi\big{)}\;.$
(4)

The right-hand side involves the total perturbation in the gravitational
potential, which is the sum of the perturber’s potential $\phi_{P}$ and the
perturbation $\delta\phi$ induced by the trailing density wake. For infinite,
homogeneous systems, the right-hand side is simply the total density
perturbation $\rho_{P}+\delta\rho$. For realistic, finite and inhomogeneous
systems, there is an additional ”non-local” term
$\boldsymbol{\nabla}\rho_{0}\cdot\boldsymbol{\nabla}(\phi_{p}+\delta\phi)$
which significantly complicates the analysis. Note also that the contribution
of the unperturbed gravitational field $\boldsymbol{\nabla}\phi_{0}$ is
neglected in the treatment of the infinite, homogeneous case following the
”Jeans swindle”. Eq. (4) describes a driven harmonic oscillator, but we shall
refer to it as the ”sound wave” equation for simplicity.

For the polytropes considered here, the sound wave equation becomes

$\ddot{\delta}\rho-\nabla^{2}(c_{0}^{2}\delta\rho)-4\pi
G\rho_{0}\delta\rho-\boldsymbol{\nabla}\phi_{0}\cdot\boldsymbol{\nabla}\delta\rho=S(\rho_{0},\phi_{p},\delta\phi)$
(5)

where
$S(\rho_{0},\phi_{p},\delta\phi)=\boldsymbol{\nabla}\big{(}\rho_{0}\boldsymbol{\nabla}(\phi_{p}+\delta\phi)\big{)}$
is the source term. We now apply Green’s method and look for the Green’s
function $G(\textbf{\em r},\textbf{\em r}^{\prime},t,t^{\prime})$ which solves

$\left(\partial_{t}^{2}-\nabla^{2}c_{0}^{2}-4\pi
G\rho_{0}-\boldsymbol{\nabla}\phi_{0}\cdot\boldsymbol{\nabla}\right)G=\delta^{D}(t-t^{\prime})\delta^{D}(\textbf{\em
r}-\textbf{\em r}^{\prime})$ (6)

Here, $\boldsymbol{\nabla}\phi_{0}\cdot\boldsymbol{\nabla}$ involves only the
radial part owing to the spherical symmetry of the unperturbed profile, and it
is understood that $\nabla^{2}$ acts on the product $c_{0}^{2}G$. Physically,
$G$ is the polarization function (usually denoted by $P$ as in Binney &
Tremaine, 1987) since it relates the wake overdensity $\delta\rho$ to the
total perturbation.

The polarization function $G$ depends on zeroth order quantities and,
therefore, inherits the stationarity (invariance under time-translations) and
isotropy (invariance under rotations) of the hydrostatic equilibrium.
Therefore, $G$ must be of the form

$G(\textbf{\em r},\textbf{\em r}^{\prime},t,t^{\prime})=\sum_{\ell
m}\int_{\omega}g_{\omega\ell}(r,r^{\prime})e^{-i\omega(t-t^{\prime})}Y_{\ell}^{m}(\mathrm{\hat{\bf{r}}})Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}^{\prime})\;.$
(7)

Here, $g_{\omega\ell}(r,r^{\prime})$ are the radial components of the
polarization function, $Y_{\ell}^{m}(\mathrm{\hat{\bf{r}}})$ are spherical
harmonics, and the shorthand for the integral over the (circular) frequency of
the sound waves $\omega$ is
$\int_{\omega}\equiv\frac{1}{2\pi}\int_{-\infty}^{+\infty}\\!d\omega$.
Substituting this ansatz into Eq. (6) and using the representations of the
Dirac distributions $\delta^{D}(t-t^{\prime})$ and $\delta^{D}(\textbf{\em
r}-\textbf{\em r}^{\prime})$ in terms of plane waves and spherical harmonics,

$\displaystyle\delta^{D}(t-t^{\prime})$
$\displaystyle=\int_{\omega}e^{-i\omega(t-t^{\prime})}$ (8)
$\displaystyle\delta^{D}(\textbf{\em r}-\textbf{\em r}^{\prime})$
$\displaystyle=\frac{1}{r^{2}}\delta^{D}(r-r^{\prime})\sum_{\ell
m}Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}})Y_{\ell}^{m}(\mathrm{\hat{\bf{r}}}^{\prime})\;,$
(9)

the radial Green’s equation can be recast into

$-r\partial_{r}^{2}\Big{(}r\,c_{0}^{2}\,g_{\omega\ell}\Big{)}-r^{2}\,\partial_{r}\phi_{0}\,\partial_{r}g_{\omega\ell}-\Big{[}\Big{(}\omega^{2}+c_{0}^{2}k_{J}^{2}\Big{)}r^{2}\\\
-\ell\big{(}\ell+1\big{)}c_{0}^{2}\Big{]}g_{\omega\ell}=\delta^{D}(r-r^{\prime})$
(10)

where $k_{J}^{2}(r)=4\pi G\rho_{0}(r)/c_{0}^{2}(r)$ is a (position dependent)
Jeans wavenumber.

####  2.1.2 Acoustic wake and DF force

Once the polarization function is known, the (over)density
$\delta\rho(\textbf{\em r},t)$ of the acoustic wake is determined by the
(integro-differential) equation

$\delta\rho(\textbf{\em
r},t)=\int\\!\\!d^{3}r^{\prime}\,\int\\!\\!dt^{\prime}\,G(\textbf{\em
r},\textbf{\em
r}^{\prime},t-t^{\prime})\,S(\rho_{0}^{\prime},\phi_{P}^{\prime},\delta\phi^{\prime})$
(11)

where the primes signify that the source term is evaluated at the coordinates
$(\textbf{\em r}^{\prime},t^{\prime})$. When the self-gravity of the density
wake can be neglected, we have
$S\approx\boldsymbol{\nabla}(\rho_{0}\boldsymbol{\nabla}\phi_{P})$ and the
polarization function is equal to the response function of the system:

$\delta\rho(\textbf{\em
r},t)\approx\int\\!\\!d^{3}r^{\prime}\,\int\\!\\!dt^{\prime}\,G(\textbf{\em
r},\textbf{\em
r}^{\prime},t-t^{\prime})\,\boldsymbol{\nabla}\big{(}\rho_{0}^{\prime}\boldsymbol{\nabla}\phi_{P}^{\prime}\big{)}\;.$
(12)

For the pointlike perturber considered here, the solution to Eq. (12) in terms
of the density multipoles $\delta\rho_{\ell}^{m}(r,t)$ reads

$\delta\rho_{\ell}^{m}(r,t)=4\pi
GM\int_{\omega}\int\\!\\!dt^{\prime}\,e^{-i\omega(t-t^{\prime})}\,h(t^{\prime})\\\
\times
Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}_{p}^{\prime})A_{\omega\ell}(r,r_{p}^{\prime})$
(13)

where the function $A_{\omega\ell}(r,r_{p}^{\prime})$ is given by

$\displaystyle A_{\omega\ell}(r,r_{p}^{\prime})$
$\displaystyle=\rho_{0}(r_{p}^{\prime})\,g_{\omega\ell}(r,r_{p}^{\prime})$
(14)
$\displaystyle\quad-\frac{1}{2\ell+1}\int\\!\\!dss^{2}\,g_{\omega\ell}(r,s)\partial_{s}\rho_{0}(s)\partial_{s}\\!\\!\left(\frac{s_{<}^{\ell}}{s_{>}^{\ell+1}}\right)$

Here, $r_{p}^{\prime}=|\textbf{\em r}_{p}^{\prime}|$ and
$\mathrm{\hat{\bf{r}}}_{p}^{\prime}$ is the unit vector along $\textbf{\em
r}_{p}^{\prime}$. Furthermore, $s_{<}={\rm min}(s,r_{p}^{\prime})$ and
$s_{>}={\rm max}(s,r_{p}^{\prime})$. The acoustic wake pattern contributes a
gravitational potential

$\delta\phi(\textbf{\em
r},t)=-G\int\\!\\!d^{3}r^{\prime}\frac{\delta\rho(\textbf{\em
r}^{\prime},t)}{|\textbf{\em r}-\textbf{\em r}^{\prime}|}\;,$ (15)

the gradient of which gives the DF force experienced by the perturber. Using
the multipole expansion of $|\textbf{\em r}-\textbf{\em r}^{\prime}|^{-1}$,
the multipoles $\delta\phi_{\ell}^{m}(r,t)$ of this self-gravity perturbation
can be recast into the form

$\displaystyle\delta\phi_{\ell}^{m}(r,t)$ $\displaystyle=-\frac{4\pi
G}{2\ell+1}\int\\!\\!ds\,s^{2}\,\frac{s_{<}^{\ell}}{s_{>}^{\ell+1}}\,\delta\rho_{\ell}^{m}(s,t)$
(16)

where $s_{<}={\rm min}(s,r)$ and $s_{>}={\rm max}(s,r)$.

The perturbation $\delta\phi$ generates a gravitational drag force (DF) on the
perturber given by

${\bf
F}_{\text{DF}}(t)=F_{r}(t)\,\mathrm{\hat{\bf{r}}}(t)+F_{\varphi}(t)\,\hat{\boldsymbol{\varphi}}(t)$
(17)

where $\mathrm{\hat{\bf{r}}}(t)=\textbf{\em r}_{p}(t)/|\textbf{\em r}_{p}(t)|$
and $\hat{\boldsymbol{\varphi}}(t)$ are unit vectors in the radial and
tangential direction. By definition, the azimuthal unit vector
$\hat{\boldsymbol{\varphi}}(t)$ is oriented counterclockwise.

The tangential component $F_{\varphi}(t)$ of the DF force at time $t$ can be
obtained from the torque ${\boldsymbol{\tau}}(t)$ exerted by the acoustic wake
on the perturber:

$\displaystyle{\boldsymbol{\tau}}(t)$ $\displaystyle=-M\textbf{\em
r}_{p}\times\boldsymbol{\nabla}\delta\phi(\textbf{\em r}_{p},t)$ (18)
$\displaystyle=-iML_{z}\delta\phi(\textbf{\em
r}_{p},t)\,\mathrm{\hat{\bf{z}}}$ $\displaystyle\equiv
aF_{\varphi}(t)\,\mathrm{\hat{\bf{z}}}\;,$ (19)

where $L_{z}\equiv-i\partial/\partial\varphi$ is the $z$-component of the
angular momentum operator ${\bf L}=-i\,\textbf{\em
r}\times\boldsymbol{\nabla}$. Therefore,

$F_{\varphi}(t)=-i\frac{M}{a}L_{z}\delta\phi(\textbf{\em r}_{p},t)$ (20)

which, in essence, is the Lynden-Bell & Kalnajs (1972) formula for the DF
torque. We shall also consider the radial component given by

$F_{r}(t)=-M\partial_{r}\delta\phi(\textbf{\em r}_{p},t)\;,$ (21)

although it is dynamically irrelevant for circular orbits. These expressions
are completely general insofar it applies to any spherically symmetric medium
and perturber motion.

###  2.2 Consistency check: the homogeneous medium

The acoustic wake pattern and DF force can be computed easily in the case of
steady-state circular motion in an infinite homogeneous background of density
$\bar{\rho}$. This case has been explored in e.g. Kim & Kim (2007); Kim et al.
(2008); Desjacques et al. (2022) (see Buehler et al., 2024; O’Neill et al.,
2024, for eccentric motion). Ignoring the self-gravity of the medium (which
leads to a frequency shift), the radial Green’s equation simplifies to

$-\big{(}r^{2}g_{\omega\ell}^{\prime}\big{)}^{\prime}-\left[\left(\frac{\omega}{c_{s}}\right)^{2}r^{2}-l\big{(}l+1)\right]g_{\omega\ell}\\\
=\frac{\delta^{D}(r-r^{\prime})}{c_{s}^{2}}$ (22)

which is of the Sturm-Liouville type. For purely outgoing wave boundary
conditions at infinity and a finite response in the limit $r\to 0$, the radial
polarization function takes the form

$g_{\omega\ell}(r,r^{\prime})=\frac{i}{c_{s}^{2}}\bigg{(}\frac{\omega}{c_{s}}\bigg{)}\,j_{\ell}\\!\bigg{(}\frac{\omega}{c_{s}}r_{<}\bigg{)}\,h_{\ell}^{(1)}\\!\bigg{(}\frac{\omega}{c_{s}}r_{>}\bigg{)}$
(23)

for $\omega\neq 0$. Here, $r_{<}={\rm min}(r,r^{\prime})$ and $r_{>}={\rm
max}(r,r^{\prime})$, whereas $j_{\ell}(x)$ and $h_{\ell}^{(1)}(x)$ are
spherical Bessel and Hankel functions, respectively. When $\omega=0$, the
solution is a powerlaw. We find

$g_{0\ell}(r,r^{\prime})=\frac{c_{s}^{-2}}{(2\ell+1)}\frac{r_{<}^{\ell}}{r_{>}^{\ell+1}}$
(24)

upon requiring that $g_{0\ell}$ decays to zero at infinity.

Figure 1: Tangential (left panel) and radial (right panel) friction
coefficient $I_{\varphi}(\mathcal{M})$ and $I_{r}(\mathcal{M})$ as a function
of the Mach number $\mathcal{M}$ in the steady-state regime for a circularly-
moving perturber in an infinite homogeneous background. The dashed curves
represent the contribution from the inner and outer regions (assuming
$\ell_{\text{max}}=12$), whereas the (cyan) dotted and (green) dot-dashed
curves are the contribution from the $\ell=0$ and $\ell\leq 1$ multipoles,
respectively. The data is from Kim & Kim (2007) (see text for details).

On setting $h(t^{\prime})\equiv~{}1$ and $\rho(\textbf{\em
r}_{p}^{\prime})\equiv\bar{\rho}$ in Eq. (13), the integral over $t^{\prime}$
simplifies to

$\int\\!\\!dt^{\prime}\,e^{-i\omega(t-t^{\prime})}e^{-im\Omega
t^{\prime}}=2\pi\,\delta^{D}(\omega-m\Omega)\,e^{-im\Omega t}\;,$

where the second complex exponential arises from

$Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}_{p}^{\prime})=Y_{\ell}^{m*}\\!\left(\frac{\pi}{2},0\right)e^{-im\Omega
t^{\prime}}\;.$ (25)

Note that $Y_{\ell}^{0}(\pi/2,0)=0$ for $\ell$ even, which reflects the planar
symmetry of the response, i.e. $\delta\rho$ is invariant under a reflexion
along the $z$ axis. Inserting this result into the expression of $\delta\rho$
and taking advantage of the fact that $r_{p}=r_{p}^{\prime}=a$ for circular
motion, the multipole expansion of the acoustic wake density leads to

$\delta\rho_{\ell}^{m}(r,t)=4\pi
GM\,\bar{\rho}\,g_{m\Omega,\ell}(r,a)\,Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}_{p})$
(26)

Observe the appearance of all the higher harmonics $m\Omega$ of the
fundamental (orbital) frequency. Furthermore,
$g_{m\Omega,\ell}(r,a)^{*}=g_{-m\Omega,\ell}(r,a)$, which ensures that
$\delta\rho$ is real.

The monopole ($\ell=0$) contribution to $\delta\rho$ is independent of the
Mach number. It reflects the accretion onto the perturber’s orbit and, in
steady-state, leads to an overdensity

$\delta\rho_{\ell=0}(\textbf{\em
r},t)\equiv\delta\rho_{0}^{0}(r,t)\,Y_{0}^{0}(\mathrm{\hat{\bf{r}}})=\frac{GM}{c_{s}^{2}\,r_{>}}\bar{\rho}$
(27)

where $r_{>}={\rm max}(r,a)$. Furthermore, while it does not contribute to the
tangential DF force $F_{\varphi}(t)$, it yields a force in the radial
direction which, as we will see shortly, is the dominant contribution to
$F_{r}(t)$ in the limit $\mathcal{M}\to\infty$.

The multipole coefficients of the perturbation $\delta\phi$ to the medium’s
gravitational potential follows from Eq. (15) and read

$\delta\phi_{\ell}^{m}(r,t)=-\big{(}4\pi
G\big{)}^{2}M\bar{\rho}\,\big{(}2\ell+1\big{)}^{-1}\,Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}_{p})\\\
\times\int_{0}^{\infty}\\!\\!ds\,s^{2}\frac{s_{<}^{\ell}}{s_{>}^{\ell+1}}\,g_{m\Omega,\ell}(s,a)\;,$
(28)

where $s_{<}={\rm min}(s,r)$ and $s_{>}={\rm max}(s,r)$.

Starting from the relation Eq. (20) and using the property
$L_{z}Y_{\ell}^{m}=mY_{\ell}^{m}$, and the relations

$\displaystyle\int_{0}^{1}\\!\\!du\,u^{\ell+2}\,j_{\ell}(xu)$
$\displaystyle=\frac{1}{x}\,j_{\ell+1}(x)$ (29)
$\displaystyle\int_{1}^{\infty}\\!\\!du\,u^{-\ell+1}\,h_{\ell}^{(1)}\\!(xu)$
$\displaystyle=\frac{1}{x}\,h_{\ell-1}^{(1)}\\!(x)$

which emerge after splitting the integral over $s$ into the intervals $[0,a]$
and $[a,\infty[$, the tangential DF force can be recast into

$F_{\varphi}(t)=-4\pi\left(\frac{GM}{\Omega
a}\right)^{2}\bar{\rho}\,I_{\varphi}(\mathcal{M})$ (30)

where

$\displaystyle I_{\varphi}(\mathcal{M})$
$\displaystyle=\mathcal{M}^{2}\sum_{\ell=1}^{\infty}\sum_{m\neq
0}m\,\frac{(\ell-m)!}{(\ell+m)!}\,\big{\lvert}P_{\ell}^{m}(0)\big{\lvert}^{2}$
(31)
$\displaystyle\quad\times\Big{[}j_{\ell+1}(m\mathcal{M})h_{\ell}^{(1)}\\!(m\mathcal{M})+j_{\ell}(m\mathcal{M})h_{\ell-1}^{(1)}\\!(m\mathcal{M})\Big{]}$

is the tangential friction coefficient. The first (second) term in the square
brackets of Eq. (31) is the contribution from the acoustic wake inside
(outside) the orbital radius. The monopole does not contribute to
$I_{\varphi}(\mathcal{M})$ due to Newton’s shell (or Birkhoff’s) theorem. The
dipole ($\ell=1$) term dominates the friction for $\mathcal{M}\ll 1$. In
addition, for $\ell\geq 1$, all the contributions with $m=0$ vanish because
$F_{\varphi}$ cannot be sourced by a density perturbation with azimuthal
symmetry. Finally, in the limit $\mathcal{M}\to 0$, the tangential friction
scales as $I_{\varphi}(\mathcal{M})\sim\mathcal{M}^{3}/3$, which implies that
the tangential drag is proportional to the perturber’s velocity for $v_{p}\ll
c_{s}$.

The friction coefficient $I_{\varphi}(\mathcal{M})$ is displayed in the left
panel of Fig. 1 as the solid (red) curve. It is compared to the tangential DF
extracted from a numerical schemes implementing the Liénard-Wiechert potential
(Kim & Kim, 2007) as indicated by the triangles. The agreement is excellent as
was already demonstrated in Desjacques et al. (2022). For supersonic motion,
$I_{\varphi}(\mathcal{M})$ is affected by a short distance Coulomb divergence,
which is regularized by truncating the series expansion Eq. (31) at a maximum
multipole $\ell_{\text{max}}$. Choosing $\ell_{\text{max}}=12$ yields a
reasonable fit to the simulation data.

Similarly, the radial DF force can be expressed as

$F_{r}(t)=-4\pi\left(\frac{GM}{\Omega
a}\right)^{2}\bar{\rho}\,I_{r}(\mathcal{M})$ (32)

where

$\displaystyle I_{r}(\mathcal{M})$
$\displaystyle=\frac{\mathcal{M}^{2}}{3}+i\mathcal{M}^{2}\sum_{\ell=1}^{\infty}\sum_{m}\frac{(\ell-m)!}{(\ell+m)!}\,\big{\lvert}P_{\ell}^{m}(0)\big{\lvert}^{2}$
$\displaystyle\quad\times\Big{[}(\ell+1)\,j_{\ell+1}(m\mathcal{M})\,h_{\ell}^{(1)}\\!(m\mathcal{M})$
$\displaystyle\qquad-\ell\,j_{\ell}(m\mathcal{M})\,h_{\ell-1}^{(1)}\\!(m\mathcal{M})\Big{]}$
(33)

is the ”radial friction coefficient”. The term $\mathcal{M}^{2}/3$ arises from
the monopole Eq. (27) of the density perturbation $\delta\rho(\textbf{\em
r},t)$ which gives rise to a radial force

$F_{r}(t)\supset-\frac{GM}{a^{2}}\left(\frac{4\pi}{3}\frac{GM}{c_{s}^{2}a}\bar{\rho}a^{3}\right)=-4\pi\left(\frac{GM}{\Omega
a}\right)^{2}\bar{\rho}\,\frac{\mathcal{M}^{2}}{3}$ (34)

as follows immediately from Newton’s shell theorem. Furthermore, when $m=0$,
the terms in the square brackets of Eq. (2.2) should be replaced by

$\displaystyle(\ell+1)\,j_{\ell+1}\,h_{\ell}^{(1)}$
$\displaystyle\to\left(-\frac{i}{4}\right)\frac{\big{(}\ell+1\big{)}}{\left(\ell+\frac{3}{2}\right)\left(\ell+\frac{1}{2}\right)}$
(35) $\displaystyle\ell\,j_{\ell}\,h_{\ell-1}^{(1)}$
$\displaystyle\to\left(-\frac{i}{4}\right)\frac{\ell}{\left(\ell+\frac{1}{2}\right)\left(\ell-\frac{1}{2}\right)}$

for practical evaluation.

For a static perturber ($\mathcal{M}=0$), the frequencies are
$\omega=m\Omega=0$ so that the density multipoles are given by

$\delta\rho_{\ell}^{m}(r,t)=4\pi
GM\bar{\rho}\,g_{0\ell}(r,a)\,Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}_{p})\;,$
(36)

which yields an overdensity

$\delta\rho(\textbf{\em
r})=\frac{GM}{c_{s}^{2}}\bar{\rho}\,\frac{1}{\big{|}\textbf{\em r}-\textbf{\em
r}_{p}\big{|}}\qquad(\mathcal{M}=0)$ (37)

that exceeds unity within the Bondi radius $2GM/c_{s}^{2}$ of the perturber
(where linear response theory becomes invalid). This density profile
originates from the accretion onto the perturber, which is present regardless
of the perturber’s motion. Since it is isotropic around the constant perturber
position $\textbf{\em r}_{p}$, DF vanishes in the limit $\mathcal{M}\to 0$.
Therefore, the $\mathcal{M}^{2}$ term in the multipole expansion of
$I_{r}(\mathcal{M})$ must also vanish in the same limit. This translates into
a non-trivial constraint among the different multipoles. Furthermore, a series
expansion around $\mathcal{M}=0$ reveals that
$I_{r}(\mathcal{M})\propto\mathcal{M}^{4}$ for $\mathcal{M}\ll 1$. The absence
of a $\mathcal{M}^{3}$ term originates from the fact that, in the limit
$v_{p}\ll c_{s}$, the perturber can be thought of as moving on a (very short)
straight line trajectory on which there is no perpendicular DF force (see also
the discussion in Desjacques et al., 2022).

The friction coefficient $I_{r}(\mathcal{M})$ is shown in the right panel of
Fig. 1. In the limit $\mathcal{M}\gg 1$, the $m=0$ contribution to
$I_{r}(\mathcal{M})$ dominate because the perturber moves so fast that it
essentially behaves like a one-dimensional, infinitely thin ring of matter
with density $M/2\pi a$, whence the absence of any azimuthal dependence. The
dominant contribution originates from the monopole. In the limit
$\mathcal{M}\to 0$, the vanishing of $I_{r}(\mathcal{M})$ is not perfect
because multipoles up to $\ell_{\text{max}}=12$ solely have been included in
order to match the resolution of the numerical data.

We have numerically checked that $I_{\varphi}(\mathcal{M})$ and
$I_{r}(\mathcal{M})$ precisely match the imaginary and real part,
$\Im\big{(}I(\mathcal{M})\big{)}$ and $\Re\big{(}I(\mathcal{M})\big{)}$, of
the (complex) friction coefficient $I(\mathcal{M})$ calculated by Desjacques
et al. (2022). However, our result provides explicit expressions for the
contributions from the $r<a$ and $r>a$ region: they are given by the first and
the second term in the square brackets of Eqs. (31) and (2.2).

##  3 The singular isothermal sphere

We now turn to the singular isothermal sphere (SIS) with $\gamma=1$ and
$K=c_{s}^{2}$, for which the unperturbed density, potential and pressure are

$\displaystyle\rho_{0}(r)$ $\displaystyle=\frac{c_{s}^{2}}{2\pi G}\,r^{-2}$
(38) $\displaystyle\phi_{0}(r)$ $\displaystyle=2c_{s}^{2}\,{\rm
ln}(r)+\mbox{const}$

and $P_{0}(r)=c_{s}^{2}\rho_{0}(r)$. This yields a Jeans wavenumber
$k_{J}^{2}=2/r^{2}$. The Green’s equation is more complex in this case but, as
we shall see below, the methodology is the same as in the infinite,
homogeneous case.

###  3.1 Polarization function

For the SIS profile, the radial Green’s equation becomes

$-r^{2}g_{\omega\ell}^{\prime\prime}-4rg_{\omega\ell}^{\prime}-\bigg{[}\left(\frac{\omega}{c_{s}}\right)^{2}r^{2}-\ell(\ell+1)+2\bigg{]}\,g_{\omega\ell}\\\
=\frac{\delta^{D}(r-r^{\prime})}{c_{s}^{2}}$ (39)

The constancy of the sound speed implies $g_{\omega\ell}(\lambda r,\lambda
r^{\prime})=\lambda^{-1}g_{(\lambda\omega)\ell}(r,r^{\prime})$ for any
$\lambda>0$. The solution to the homogeneous equation (that is, for $r\neq
r^{\prime}$) reads

$g_{\omega\ell}=\frac{A}{r}\,j_{\ell}\\!\left(\frac{\omega}{c_{s}}r\right)+\frac{B}{r}\,y_{\ell}\\!\left(\frac{\omega}{c_{s}}r\right)\;.$
(40)

To construct the radial polarization function, we impose the same boundary
conditions as in the infinite homogeneous case. Since $y_{\ell}$ diverges at
the origin, $B$ must be to zero for $r<r^{\prime}$. Furthermore, the outgoing
wave boundary condition at infinity implies that the solution is proportional
to $h_{\ell}^{(1)}$ for $r>r^{\prime}$. Hence, we have

$g_{\omega\ell}(r,r^{\prime})=\begin{cases}A(r^{\prime})\,r^{-1}\,j_{\ell}\\!\left(\frac{\omega}{c_{s}}r\right)&(r<r^{\prime})\\\
B(r^{\prime})\,r^{-1}\,h_{\ell}^{(1)}\\!\left(\frac{\omega}{c_{s}}r\right)&(r>r^{\prime})\end{cases}\;.$
(41)

The continuity and jump conditions at $r=r^{\prime}$ determine $A(r^{\prime})$
and $B(r^{\prime})$ such that, for $\omega\neq 0$, the radial polarization
function is given by

$g_{\omega\ell}(r,r^{\prime})=\frac{i}{c_{s}^{2}}\,\left(\frac{\omega}{c_{s}}\right)\frac{r^{\prime}}{r}j_{\ell}\\!\bigg{(}\frac{\omega}{c_{s}}r_{<}\bigg{)}\,h_{\ell}^{(1)}\\!\bigg{(}\frac{\omega}{c_{s}}r_{>}\bigg{)}$
(42)

Except for the asymmetric term $r^{\prime}/r$ arising from the radial
dependence of the unperturbed density profile, it is identical to the
polarization function Eq. (23) of infinite, homogeneous systems. Therefore, we
have

$g_{0\ell}(r,r^{\prime})=\frac{c_{s}^{-2}}{(2\ell+1)}\frac{r^{\prime}}{r}\frac{r_{<}^{\ell}}{r_{>}^{\ell+1}}$
(43)

in the special case $\omega=0$.

Figure 2: Slice through the density profile of the acoustic wake for a static
perturber ($\mathcal{M}=0$). The unperturbed density is normalized to
$\rho_{0}(a)=1$. A value of $GM/c_{s}^{2}a=0.1$ is assumed. Dotted (blue)
lines indicate negative values (of $\delta\rho$).

###  3.2 Acoustic wake

In the steady-state limit, the multipole coefficients of the wake
(over)density are given by

$\delta\rho_{\ell}^{m}(r,t)=4\pi
GM\,A_{m\Omega,\ell}(r,a)\,Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}_{p})$ (44)

where $A_{m\Omega,\ell}(r,a)$ is the auxiliary function, Eq. (14) with
$\omega=m\Omega$ and $r_{p}^{\prime}\equiv a$. Introducing

$\displaystyle\mathcal{J}_{\ell}^{n}(x;a,b)$
$\displaystyle=\int_{a}^{b}\\!\\!du\,u^{n}\,j_{\ell}(xu)$ (45)
$\displaystyle\mathcal{H}_{\ell}^{n}(x;a,b)$
$\displaystyle=\int_{a}^{b}\\!\\!du\,u^{n}\,h_{\ell}^{(1)}\\!(xu)$

for shorthand convenience, $A_{m\Omega,\ell}(r,a)$ is given by

$\displaystyle A_{m\Omega,\ell}(r,a)$
$\displaystyle=\rho_{0}(a)\frac{im\mathcal{M}}{c_{s}^{2}r}\bigg{\\{}j_{\ell}\\!\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}\,h_{\ell}^{(1)}\\!(m\mathcal{M})$
(46)
$\displaystyle\quad+\frac{2}{2\ell+1}\bigg{[}\ell\,h_{\ell}^{(1)}\\!\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}\,\mathcal{J}_{\ell}^{\ell-1}\\!\Big{(}m\mathcal{M};0,\frac{r}{a}\Big{)}$
$\displaystyle\quad+\ell\,j_{\ell}\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}\,\mathcal{H}_{\ell}^{\ell-1}\\!\Big{(}m\mathcal{M};\frac{r}{a},1\Big{)}$
$\displaystyle\quad-\big{(}\ell+1\big{)}j_{l}\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}\,\mathcal{H}_{\ell}^{-\ell-2}\\!\big{(}m\mathcal{M};1,\infty\big{)}\bigg{]}\bigg{\\}}$

for $r<a$, and

$\displaystyle A_{m\Omega,\ell}(r,a)$
$\displaystyle=\rho_{0}(a)\frac{im\mathcal{M}}{c_{s}^{2}r}\bigg{\\{}j_{\ell}(m\mathcal{M})\,h_{\ell}^{(1)}\\!\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}$
(47)
$\displaystyle\quad+\frac{2}{2\ell+1}\bigg{[}\ell\,h_{\ell}^{(1)}\\!\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}\,\mathcal{J}_{\ell}^{\ell-1}\\!\big{(}m\mathcal{M};0,1\big{)}$
$\displaystyle\quad-\big{(}\ell+1\big{)}h_{\ell}^{(1)}\\!\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}\,\mathcal{J}_{\ell}^{-\ell-2}\\!\Big{(}m\mathcal{M};1,\frac{r}{a}\Big{)}$
$\displaystyle\quad-\big{(}\ell+1\big{)}j_{l}\Big{(}m\mathcal{M}\frac{r}{a}\Big{)}\,\mathcal{H}_{\ell}^{-\ell-2}\\!\Big{(}m\mathcal{M};\frac{r}{a},\infty\Big{)}\bigg{]}\bigg{\\}}$

for $r>a$. Eqs .(44) – (47) are the equivalent of the density response
computed in Weinberg (1986); Kaur & Stone (2022) for self-gravitating
collisionless systems. For the boundary conditions assumed here (steady-
state), the acoustic wake pattern is purely resonant and depends only on
harmonics of the fundamental frequency $\Omega$.

At large separations $r\gg a$, the outgoing wave (encoded by the Hankel
functions) behaves like $h_{\ell}^{(1)}(x)\sim(-i)^{\ell+1}e^{ix}/x\sim a/r$,
which implies that the fractional density perturbation $\delta\rho(\textbf{\em
r},t)/\rho_{0}(r)$ asymptotes to a constant at large radii. Using Eq. (38), we
find

$\delta\rho(\textbf{\em r},t)\sim\epsilon\,\rho_{0}(r)\quad\mbox{for}\quad
r\gg a\;,$ (48)

where

$\epsilon\equiv\frac{GM}{c_{s}^{2}a}\;.$ (49)

Notice also that $M/M(<a)=\epsilon/2$, where $M(<a)$ is the mass enclosed in
the profile within a radius $a$. Therefore, we must have $\epsilon\ll 1$ to
ensure the validity of linear response theory, at least across a reasonable
range of radii. In practice, the medium is not an ideal fluid and sound waves
are damped by viscosity before they escape to infinity.

Figure 3: Steady-state acoustic wake overdensity $\delta\rho/\rho_{0}$ in the
orbital plane of a perturber moving circularly with a Mach number of
$\mathcal{M}=0.2$ (left panels) and $\mathcal{M}=2$ (right panels). The left
panels display the dipole contribution, while the right panels show the
acoustic wake for $\ell_{\text{max}}=8$. The mass $M$ of the perturber, the
sound speed and the orbital radius $a$ are chosen such that $\epsilon=10^{-4}$
and $10^{-4}$ for $\mathcal{M}=0.2$ and 2, respectively. The perturber’s
position is represented as a white symbol, and its circular orbit as a white
circle.

The monopole of the density wake follows from $A_{0,0}(r,t)$, which turns out
to vanish within the orbital radius. Namely,

$A_{0,0}(r,a)=\frac{1}{2\pi
G}\frac{1}{r^{2}}\left(\frac{1}{r}-\frac{1}{a}\right)\Theta\big{(}r-a\big{)}$
(50)

where $\Theta(x)$ is the Heaviside step function. As a result, the angle-
average wake density profile is

$\delta\rho_{\ell=0}(r,t)=\frac{GM}{c_{s}^{2}}\rho_{0}(r)\left(\frac{1}{r}-\frac{1}{a}\right)\Theta\big{(}r-a\big{)}$
(51)

regardless of the value of $\mathcal{M}$. This shows that the outer region of
the profile is depleted as gas is accreted onto the perturber. By contrast,
the mass $M(<a)$ enclosed in the inner region is unchanged. Note that, for
$\ell\geq 1$, we have

$A_{0,\ell}(r,a)=\frac{c_{s}^{2}}{2\ell+1}\,\rho_{0}(r)\,\frac{r_{<}^{\ell}}{r_{>}^{\ell+1}}\;.$
(52)

Therefore, in the limit where the wake self-gravity can be neglected, the
acoustic wake overdensity is

$\displaystyle\delta\rho(\textbf{\em r})$ $\displaystyle=4\pi GM\sum_{\ell
m}A_{0,\ell}(r,t)\,Y_{\ell}^{m*}(\mathrm{\hat{\bf{r}}}_{p})\,Y_{\ell}^{m}(\mathrm{\hat{\bf{r}}})$
(53)
$\displaystyle=\frac{GM}{c_{s}^{2}}\,\rho_{0}(r)\left(\frac{1}{\big{|}\textbf{\em
r}-\textbf{\em r}_{p}\big{|}}-\frac{1}{a}\right)\quad(\mathcal{M}=0)$

for a static perturber. This overdensity is shown in Fig. 2, together with the
unperturbed and perturbed density profiles $\rho_{0}$ and
$\rho=\rho_{0}+\delta\rho$, along the line of symmetry going through the
center of the SIS and the perturber. A value of $GM/c_{s}^{2}a=0.1$ is assumed
for illustration.

The acoustic wave patterns shown in Fig 3 illustrate the response of the self-
gravitating system for supersonic and subsonic motion. Results are shown for
$\mathcal{M}=0.2$ (left panels) and 2 (right panels). For $\mathcal{M}>1$, the
acoustic wake is confined to a Mach cone, and its density diverges near the
perturber’s position. By contrast, a perturber moving at subsonic speed
imprints a thinner, slightly overdense spiral tail analogous to the infinite,
homogeneous result. The fractional density $\delta\rho/\rho_{0}$ converge to
zero toward the center of the profile. However, one should bear in mind that
$\delta\rho\sim r^{-1}$ for $r\ll a$ due to the divergence of the SIS model.
The first terms in the right-hand side of Eqs. (46) and (47) dominate the
acoustic wake density, Eq. (44). They correspond to the ”local source”
$\rho_{0}\boldsymbol{\nabla}^{2}\delta\phi_{P}$ in the quantity
$S(\rho_{0},\phi_{p},\delta\phi)$ (see Eq. 5). By contrast, the remaining
terms proportional to $2/(2\ell+1)$ originate from the ”nonlocal source”
$\boldsymbol{\nabla}\rho_{0}\cdot\boldsymbol{\nabla}\delta\phi_{P}$, and
contribute to the total wake density at the $\sim 10$% level.

Figure 4: Tangential (left panel) and radial (right panel) friction
coefficient $I_{\varphi}(\mathcal{M})$ and $I_{r}(\mathcal{M})$ as a function
of the Mach number $\mathcal{M}$ in the steady-state regime for a circularly-
moving perturber in a SIS profile. The dashed curves indicate the contribution
from the inner and outer regions (assuming $\ell_{\text{max}}=12$). The
(green) dot-dashed curve is the dipole contribution. For comparison, the
infinite, homogeneous result is shown as the (red) dotted curve.

###  3.3 Friction coefficients

The perturbation $\delta\phi$ to the gravitational potential of the isothermal
sphere induced by a circularly-moving perturber at radius $a$ follows from Eq.
(16). First, we compute the multipole coefficients $\delta\phi_{\ell}^{m}$
assuming $r<a$. Next, we obtain the components $F_{\varphi}$ and $F_{r}$ of
the DF force upon applying Eqs. (20) and (21), and subsequently taking the
limit $r\to a$.

The DF force can eventually be expressed as

${\bf F}_{\text{DF}}(t)=-4\pi\left(\frac{GM}{\Omega
a}\right)^{2}\rho_{0}(a)\,\Big{(}I_{r}\,\mathrm{\hat{\bf{r}}}(t)+I_{\varphi}\,\hat{\boldsymbol{\varphi}}(t)\Big{)}\;.$
(54)

In the steady-state regime considered here, the tangential friction
coefficient is given by

$\displaystyle I_{\varphi}(\mathcal{M})$
$\displaystyle=\mathcal{M}^{2}\sum_{\ell=1}^{\infty}\sum_{m\neq
0}m\,\frac{(\ell-m)!}{(\ell+m)!}\,\big{\lvert}P_{\ell}^{m}(0)\big{\lvert}^{2}$
(55)
$\displaystyle\qquad\times\Big{[}\mathcal{I}_{\ell}^{-}(m\mathcal{M})+\mathcal{I}_{\ell}^{+}\\!(m\mathcal{M})\Big{]}\;.$

The functions $\mathcal{I}_{\ell}^{\pm}\\!(x)$ defined in Appendix §A are
related to the contribution from the inner ($-$) and outer ($+$) region.
Likewise, the radial friction coefficient is

$\displaystyle I_{r}(\mathcal{M})$
$\displaystyle=i\mathcal{M}^{2}\sum_{\ell=1}^{\infty}\sum_{m}\frac{(\ell-m)!}{(\ell+m)!}\,\big{\lvert}P_{\ell}^{m}(0)\big{\lvert}^{2}$
(56)
$\displaystyle\qquad\times\Big{[}(\ell+1)\,\mathcal{I}_{\ell}^{-}\\!(m\mathcal{M})-\ell\,\mathcal{I}_{\ell}^{+}\\!(m\mathcal{M})\Big{]}\;.$

Unlike the infinite, homogeneous case, the monopole contribution vanishes
because there is no excess mass within the orbital radius. For practical
evaluations, the terms in the square brackets of Eq. (2.2) should be replaced
by

$\displaystyle(\ell+1)\,\mathcal{I}_{\ell}^{-}$
$\displaystyle\to-\frac{i\,(\ell+1)}{\big{(}2\ell+1\big{)}^{2}}$ (57)
$\displaystyle\ell\,\mathcal{I}_{\ell}^{+}$
$\displaystyle\to-\frac{i\,\ell}{\big{(}2\ell+1\big{)}^{2}}$

when $m=0$. Appendix §A outlines how the limiting behavior of
$\mathcal{I}^{\pm}(x)$ can be derived. Note also that
$I_{\varphi}(\mathcal{M})$ ($I_{r}(\mathcal{M})$) generally is an odd (even)
function of $\mathcal{M}$ since a change in the direction of the circular
motion leaves $F_{r}$ unchanged while it flips the sign of $F_{\varphi}$.
Therefore, a series expansion of $I_{\varphi}(\mathcal{M})$
($I_{r}(\mathcal{M})$) around $\mathcal{M}=0$ only contains odd (even) powers
of $\mathcal{M}$.

Fig. 4 displays the friction coefficients $I_{\varphi}(\mathcal{M})$ (left
panel) and $I_{r}(\mathcal{M})$ (right panel) as a function of the Mach
number. The dashed curves indicate the contributions from the inner ($r<a$)
and outer ($r>a$) region, while the dashed-dotted curve shows the contribution
from the dipole. For comparison, the infinite, homogeneous results are also
shown as the dotted curves. In the SIS case,
$I_{r}(\mathcal{M})/\mathcal{M}^{2}$ does not vanish in the limit
$\mathcal{M}\to 0$ due to the radial asymmetry of the acoustic wake around the
perturber (see figure 2). For $\mathcal{M}\ll 1$, we have

$\mathcal{I}_{r}(\mathcal{M})\approx\mathcal{M}^{2}\sum_{\ell=1}^{\ell_{\text{max}}}\sum_{m}\frac{\big{(}\ell-m\big{)}!}{\big{(}\ell+m\big{)}!}\frac{\big{\lvert}P_{\ell}^{m}(0)\big{\lvert}^{2}}{(2\ell+1)^{2}}\;,$
(58)

which returns $I_{r}(\mathcal{M})\approx 0.2145\mathcal{M}^{2}$ for the value
$\ell_{\text{max}}=12$ assumed in the figure. For $\mathcal{M}\gg 1$, the
radial friction $I_{r}(\mathcal{M})$ decays to zero since, unlike the infinite
homogeneous case, there is no excess mass within the orbital radius.
Interestingly, deviations from the infinite, homogeneous case also occur for
$I_{\varphi}(\mathcal{M})$, which encodes the drag relevant to the dynamic of
circular orbits. The striking feature is the suppression of
$I_{\varphi}(\mathcal{M})$ in the limit $\mathcal{M}\ll 1$ relative to the
infinite, homogeneous expectation. Namely, the tangential friction scales as
$I_{\varphi}(\mathcal{M})\sim\mathcal{M}^{5}/5$, that is, $F_{\varphi}\propto
v_{p}^{3}$ for $v_{p}\ll c_{s}$. This, however, would not affect the orbital
decay timescale of e.g. a supermassive black hole in an $r^{-2}$ isothermal
profile relative to the homogeneous case, because the circular orbit velocity
$v_{c}=\sqrt{2}c_{s}$ corresponds to a Mach number $\mathcal{M}=\sqrt{2}$ for
which the coefficients $I_{\varphi}(\mathcal{M})$ are close to each other.
Finally, note that the feature at $\mathcal{M}\approx 4.6033$ visible in Figs
1 and 4 corresponds to a critical Mach number at which the number of zeros of
the polarization function changes (see Kim & Kim, 2007, for a related
discussion).

##  4 Gravitational wave emission

An interesting aspect of our findings is the $r^{-2}$ decay of the wake
density multipoles $\delta\rho_{\ell}^{m}$ in the outer region. In a realistic
astrophysical system, such a scaling would, of course, hold only over a
limited range of scales. Notwithstanding, interesting effects could arise even
for a truncated isothermal profile, as we will illustrate now.

###  4.1 Truncated isothermal profile

Consider a circular, binary black hole (BBH) with mass $M_{2}\ll M_{1}$, such
that the massive companion $M_{1}$ is located at the center of the circular
orbit of radius $a$ drawn by the perturber $M_{2}$. This BBH is embedded in an
$1/r^{2}$ isothermal profile extending from $r_{\text{min}}$ out to a finite
radial distance $r_{\text{max}}$. The mass enclosed in a sphere of radius $a$
thus is $M(<a)=M_{1}+(2c_{s}^{2}/G)a=M_{1}(1+2/\epsilon_{1}(a))$. Here and
henceforth, $\epsilon_{i}(a)=GM_{i}/c_{s}^{2}a$ designates the ratio of (half)
the Bondi radius of $M_{i}$ to the orbital radius. For a black hole of mass
$M_{1}=10^{6}\ {\rm M_{\odot}}$ and a sound speed $c_{s}=100\ {\rm
km\,s^{-1}}$ for example, $\epsilon_{1}=1$ for $a\simeq 0.43\ {\rm pc}$. At
this orbital radius however, such a system would emit at frequencies much
smaller than the nanohertz (nHz) band proved by pulsar timing arrays (PTAs).
In general, $\epsilon_{1}\gg 1$ is required to get GW emission in the nHz
frequency band, but one should relax the assumption of hydrostatic equilibrium
to explore this interesting regime. Such an extension is beyond the scope of
this paper, and we leave it for future work. Likewise, the validity of linear
response theory requires $\epsilon_{2}(a)\ll 1$, which will not be satisfied
for small values of $a$.

###  4.2 Power radiated in GWs

The Quadrupole formula can be applied to compute the GW emission of the BBH,
had it been evolving in vacuum. In this approximation, the power radiated in
GW is

$P_{\text{gw}}=\frac{G}{5c^{5}}\big{\langle}\dddot{Q}_{ij}\dddot{Q}_{ij}\big{\rangle}\;,$
(59)

where the temporal average $\langle\dots\rangle$ is taken over several
periods, and the quadrupole moment is

$Q_{ij}=\int\\!d^{3}r\,\delta\rho(\textbf{\em
r},t)\left(x_{i}x_{j}-\frac{1}{3}r^{2}\delta_{ij}\right)\;.$ (60)

In vacuum, the circular binary would radiate GWs at a frequency $2\Omega$ with
a power (Peters, 1964)

$P_{\text{gw}}^{\text{Vac}}\approx\frac{32}{5c^{5}}\frac{G^{4}M_{2}^{2}M_{1}^{3}}{a^{5}}$
(61)

for the small mass ratio assumed here. The effect of the truncated isothermal
profile on the emission of GWs is twofold: (i) it perturbs the emission
arising from the two compact bodies through a non-vacuum metric; and (ii) it
provides an extra source of GWs through the time-varying quadrupole moment of
the acoustic wake. Let us assess the importance of (ii) in what follows.

Figure 5: Shape factor $\mathcal{F}(\mathcal{M},r_{+}/a)$ (Eq. (74)) as a
function of Mach number for three different values of $r_{+}/a$ (see text for
details) as quoted in the figure. The thick and thin curves represent the
predictions obtained for a truncated isothermal profile and for a uniform
sphere, respectively (see text for details). The vertical line indicate the
Mach number $\mathcal{M}=\sqrt{2}$ of a circularly-moving perturber in an
isothermal profile when the latter dominates the mass of the system.

The acoustic wake induced by the motion of the small companion $M_{2}$ (the
perturber) in the medium also produces GWs with a frequency
$\omega_{\text{gw}}^{2}\sim GM(<a)/a^{3}$ or, equivalently, a (reduced)
wavelength $\lambdabar_{\text{gw}}^{2}\sim a^{3}c^{2}/GM(<a)$. Let
$r_{-}<a<r_{+}<r_{\text{max}}$ be the radius of the region within which the
acoustic wake has not dissipated. GW emission by the wake will be coherent
only if $\lambdabar_{\text{gw}}/r_{+}\gg 1$. This condition is best satisfied
if $r_{+}\sim\alpha a$ scales with the orbit and $\alpha\lesssim 10$. For the
astrophysical system considered above for instance, taking $a=10^{-2}\ {\rm
pc}$ gives $\lambdabar_{\text{gw}}\sim 4.6\ {\rm pc}$ is still appreciably
larger than $r_{+}$ if $\alpha\sim$ a few. When
$\lambdabar_{\text{gw}}/r_{+}\gg 1$, GW emission by the wake is coherent and,
for the non-relativistic motions considered here, dominated by the quadrupole
moment of the mass distribution. The latter reads

$Q_{ij}=4\pi
GM_{2}\sum_{m}\mathcal{Q}_{ij}^{(m)}\,Y_{2}^{m*}(\mathrm{\hat{\bf{r}}}_{p})\,\int_{r_{-}}^{r_{+}}\\!\\!dr\,r^{4}\,A_{m\Omega,2}(r,a)$
(62)

where $-2\leq m\leq+2$. The integration domain of the radial integral will be
specified shortly. Only the $m=\pm 2$ multipoles contribute to GW emission.
For sake of completeness however, the matrices $\mathcal{Q}^{(m)}$ are given
by $\mathcal{Q}^{(0)}=(2/3)\sqrt{\pi/5}\,{\rm diag}\big{(}-1,-1,2\big{)}$,

$\displaystyle\mathcal{Q}^{(+1)}$
$\displaystyle=\sqrt{\frac{2\pi}{15}}\left(\begin{array}[]{ccc}0&0&-1\\\
0&0&i\\\ -1&i&0\end{array}\right)$ (66) $\displaystyle\mathcal{Q}^{(+2)}$
$\displaystyle=\sqrt{\frac{2\pi}{15}}\left(\begin{array}[]{ccc}1&-i&0\\\
-i&-1&0\\\ 0&0&0\end{array}\right)$ (70)

and $\mathcal{Q}^{(-m)}=\mathcal{Q}^{(+m)*}$. The temporal average becomes

$\displaystyle\langle\dddot{Q}_{ij}\dddot{Q}_{ij}\rangle$
$\displaystyle=\big{(}4\pi GM_{2}\big{)}^{2}\Omega^{6}\sum_{m=\pm
2}m^{6}\big{\lvert}Y_{2}^{m}(\mathrm{\hat{\bf{r}}}_{p})\big{\lvert}^{2}$ (71)
$\displaystyle\qquad\times\bigg{\lvert}\int_{r_{-}}^{r_{+}}\\!\\!dr\,r^{4}\,A_{m\Omega,2}(r,a)\bigg{\lvert}^{2}{\rm
tr}\big{(}\mathcal{Q}^{(m)}\mathcal{Q}^{(m)*}\big{)}\;,$

which follows from

$\Big{\langle}Y_{2}^{m*}(\mathrm{\hat{\bf{r}}}_{p})Y_{2}^{m^{\prime}}(\mathrm{\hat{\bf{r}}}_{p})\Big{\rangle}=\delta_{mm^{\prime}}\;.$
(72)

Taking advantage of the identity ${\rm tr}\big{(}\mathcal{Q}^{(\pm
2)}\mathcal{Q}^{(\pm 2)*}\big{)}=8\pi/15$, the GW power emitted by the
acoustic wake in the quadrupole approximation can be expressed as

$P_{\text{gw}}^{\text{Wake}}=\frac{8}{15c^{5}}\frac{G^{4}M_{2}^{2}M^{3}(<a)}{a^{5}}\,\mathcal{F}(\mathcal{M},r_{+}/a)$
(73)

where

$\mathcal{F}(\mathcal{M},r_{+}/a)=\sum_{m=\pm
2}m^{6}\frac{(2-m)!}{(2+m)!}\,\big{\lvert}P_{2}^{m}(0)\big{\lvert}^{2}\,\big{\lvert}\mathcal{A}_{m}\lvert^{2}$
(74)

is a ”shape” factor, which depends mainly on the Mach number $\mathcal{M}$ of
the perturber $M_{2}$ and the ratio $r_{+}/a$ for the $r^{-2}$ profile (and
the uniform medium)) considered below. There is, of course, no contribution
from the time-independent term $m=0$. Explicit expressions for the (complex)
quantity

$\mathcal{A}_{m}=\frac{2\pi
G}{a^{2}}\int^{r_{\text{max}}}\\!\\!dr\,r^{4}\,A_{m\Omega,2}(r,a)$ (75)

can be found in Appendix §A.

###  4.3 GW radiation from an acoustic wake

Fig. 5 displays $\mathcal{F}(\mathcal{M},r_{+}/a)$ for several values of
$r_{+}/a$ and a fixed ratio $r_{-}/a=0.5$. For sake of comparison, we also
show the prediction for an acoustic wake propagating in a uniform sphere of
radius $r_{\text{max}}$. In this case, we set the density assuming
$r_{+}=r_{\text{max}}$ and equating the sound-crossing time $r_{+}/c_{s}$ to
the free-fall time $(G\bar{\rho})^{-1}$. This gives
$\bar{\rho}=G^{-1}(c_{s}/r_{+})^{2}$. This way, the unperturbed density is
proportional to $c_{s}^{2}$ for both the truncated isothermal profile and the
uniform sphere.

We have checked that the results do not change noticeably when the inner
radius is set to $r_{-}=0$, or a region $[0.9a,1.1a]$ around the small
companion $M_{2}$ is excised. In other words, all the signals originates from
the outer region of the profile, which is unsurprising given the $r^{-2}$
dependence of $\delta\rho$. The vertical dotted line marks
$\mathcal{M}=\sqrt{2}$, which is the (constant) Mach number of a circularly-
moving perturber in a truncated isothermal profile so long as $M(<a)\gg
M_{1}$. For $\mathcal{M}\approx\sqrt{2}$, the motion of $M_{2}$ and the
acoustic wake are in resonance and maximize the shape factor
$\mathcal{F}(\mathcal{M},r_{+}/a)$. At these Mach numbers, the ratio

$\frac{P_{\text{gw}}^{\text{wake}}}{P_{\text{gw}}^{\text{Vac}}}=\frac{1}{12}\left(\frac{M(<a)}{M_{1}}\right)^{3}\,\mathcal{F}(\mathcal{M})$
(76)

can be as large as $\sim 10^{2}$ for $M(<a)\sim M_{1}$ and $r_{+}\sim 10a$,
regardless of whether the binary is embedded in an isothermal profile or a
uniform medium. At fixed $a$ and $M_{1}$, the limit $\mathcal{M}\to\infty$ is
equivalent to $c_{s}\to\infty$, that is, a vanishing profile density
$\rho_{0}$. In this regime, $P_{\text{gw}}^{\text{wake}}$ becomes negligible
as expected. The critical density $\rho_{*}$ above which
$P_{\text{gw}}^{\text{wake}}>P_{\text{gw}}^{\text{Vac}}$ decreases with
increasing $r_{+}$. For $r_{+}\lesssim 10a$ for instance, we have
$\mathcal{F}(\mathcal{M})\gtrsim 10$ for $\mathcal{M}\lesssim 10$. Therefore,
the condition $P_{\text{gw}}^{\text{wake}}>P_{\text{gw}}^{\text{Vac}}$ becomes
$c_{s}^{2}\gtrsim 10^{-2}\Omega^{2}a^{2}$, which implies
$\rho_{0}(a)>\rho_{*}\sim 10^{-2}M_{1}/a^{3}$. This shows that, in this case,
it is possible to get a significant effect even if the mass enclosed in the
profile within the orbital radius is much smaller than $M_{1}$.

In a real fluid, acoustic wakes will be dissipated by viscosity and thermal
conductivity (see e.g. Landau & Lifshitz, 1975; Narayan & Medvedev, 2001;
Fabian et al., 2003; Sijacki & Springel, 2006). Our results are valid in the
limit of small viscosity, i.e. $\gamma c_{s}/\omega\ll 1$ where
$\gamma\sim\omega^{2}/c_{s}^{3}\nu$ is the fluid absorption coefficient and
$\omega$ is the circular frequency of the sound waves (Landau & Lifshitz,
1975). In this case, the amplitude of the sound waves is damped according to
$e^{-\gamma x}$. Since the kinematic viscosity $\nu$ approximately scales as
$\nu\sim c_{s}\lambda_{\text{mfp}}$, where $\lambda_{\text{mfp}}$ is the mean
free path of the fluid particles, the condition $\gamma c_{s}/\omega\ll 1$ is
equivalent to $\lambda_{\text{mfp}}\ll c_{s}/\Omega$ or, equivalently,
$\lambda_{\text{mfp}}\ll a$. Furthermore, the condition
$\lambda_{\text{mfp}}\lesssim GM_{2}/c_{s}^{2}$ is required for the ideal
fluid approximation (see e.g. Katz et al., 2019).

For sound waves launched in the intracluster medium (ICM), the typical mean
free path of particles due to collisions, $\lambda_{\text{mfp}}\sim 0.1-10\
{\rm Kpc}$ (Schekochihin & Cowley, 2006; Andrade-Santos et al., 2013; Sanders
et al., 2013), are much larger than the orbital sizes of interest. For such
systems, one should start from the Chandrasekhar formula appropriate to
collisionless systems. In the (warm, diffuse) interstellar medium (ISM) of a
galaxy, the typical mean free path is of order $10^{15}-10^{16}{\ {\rm cm}}$,
i.e. $\lambda_{\text{mfp}}\sim 10^{-2}-10^{-3}\ {\rm pc}$ for densities $n\sim
1\ {\rm cm^{-3}}$ (McKee & Ostriker, 2007). For such densities however, the
mass in a gaseous sphere of radius $1\ {\rm pc}$ is $\sim 0.1\ {\rm
M_{\odot}}$, which is too small to influence the dynamical evolution of
compact binaries under consideration. This could drastically change if the
compact binary is embedded in a dense gaseous cloud, a gas-rich galaxy at high
redshift or, alternatively, in a dark matter core or spike.

##  5 Conclusions

We have applied linear response theory to solve for the acoustic wake induced
by a pointlike perturber in circular motion around the center of an isothermal
$r^{-2}$ density profile (for simplicity). We consider the limit of an ideal
fluid. Our work thus complements the analysis of Tremaine & Weinberg (1984);
Weinberg (1986), who focused on collisionless systems.

Our results are valid in the limit where the self-gravity of the acoustic wake
can be neglected. Taking into account the wake’s self-gravity leads to a
(Fredholm integral) equation of the form $\delta\rho=\delta\tilde{\rho}+{\rm
K}\,\delta\rho$, where $\delta\rho$ and $\delta\tilde{\rho}$ are the solution
with and without self-gravity, respectively, whereas ${\rm K}$ is a nonlocal
operator (due to Poisson’s equation $\boldsymbol{\nabla}^{2}\delta\phi=4\pi
G\delta\rho$). Depending on the properties of ${\rm K}$, this equation could
be solved iteratively (i.e. $\delta\rho=\sum_{n\geq 0}{\rm
K}^{n}\,\delta\tilde{\rho}$) in analogy with the matrix method of Kalnajs
(1971); Weinberg (1989) for collisionless systems. We have not explored this
avenue here.

To calculate the acoustic wake pattern, we have expanded the polarization
function in multipoles. This allows us to isolate the contribution of each
multipole to the dynamical friction experienced by the perturber. The outer
region $r>a$ ($a$ is the orbital radius) dominates the (tangential) drag force
$F_{\varphi}$ for subsonic motion. For Mach numbers $\mathcal{M}\ll 1$, the
dipole of the acoustic wake yields the largest contribution to the drag. In
this limit, the (tangential) friction coefficient scales as
$I_{\varphi}(\mathcal{M})\sim\mathcal{M}^{5}$ for the SIS profile, unlike the
$\mathcal{M}^{3}$ scaling in the infinite, homogeneous case. For
$\mathcal{M}\gg 1$ however, the two backgrounds yield similar results, in
agreement with the intuition that the local approximation should always be
valid deep in the supersonic regime. By contrast, predictions differ sharply
for the radial component $F_{r}$ of dynamical friction (which is dynamically
irrelevant for circular orbits). All our results assume steady-state.
Comparison with numerical simulations, which requires different boundary
conditions, is deferred to future work.

The circular velocity $\Omega a$ of the perturber is determined by the mass
$M(<a)$ enclosed within the orbit. In the absence of a massive companion at
the center, the $r^{-2}$ dependence of $\rho_{0}(r)$ leads to a flat rotation
curve and, thereby, a constant Mach number $\mathcal{M}\sim\sqrt{2}$ so long
as an adiabatic sequence of quasi-circular orbit is a reasonable approximation
to the inspiral of the perturber. At this Mach number, the acoustic wake
pattern is in resonance with the orbital motion and maximizes the drag.
Therefore, the GW emission produced by the quadrupole of the wake density
pattern should also be maximized around $\mathcal{M}\sim 1$. To estimate this
effect, we have considered an extreme mass ratio inspiral around a massive
black hole residing at the center of a (truncated) isothermal profile. Our
calculation suggests that GW emission from the acoustic wake could be large
relative to the vacuum predictions when $\mathcal{M}\sim$ a few (i.e. in the
early phase of an inspiral). These findings must be interpreted with caution
as our simplified setting does not take into account neither the viscosity of
real astrophysical fluids, nor the radial infall onto the massive companion.
Nevertheless, it would be interesting to assess how much of this signal
remains for real astrophysical systems, and the extent to which it could leave
a signal in the nanohertz frequency range probed by pulsar timing arrays.

It is a pleasure to thank Shmuel Bialy, Frans van Die, Alexei Bobrick, Adi
Nusser and Ivan Rappoport for discussions. G.E. would like to thank thank Anna
and Alfred Gray, and Mr. Vincent Meyer for their scholarships. V.D. and R.B.
acknowledge funding from the Israel Science Foundation (grant no. 2562/20).

##  Appendix A Useful relations

To derive the DF force ${\bf F}_{\text{DF}}(t)$ acting on a circularly-moving
perturber in the steady-state regime, we calculate the perturbation
$\delta\phi(\textbf{\em r},t)$ to the gravitational potential assuming $r<a$
and split accordingly the radial integral Eq. (16) into the intervals $[0,r]$,
$[r,a]$ and $[a,\infty[$. The components $F_{\varphi}$ and $F_{r}$ of the DF
force are obtained from Eqs. (20) and (21) upon taking the limit $r\to a$. In
the SIS case, this last step leads to the introduction of two (complex)
functions $\mathcal{I}_{\ell}^{\pm}(x)$ defined as

$\displaystyle\frac{1}{x}\,{\cal I}_{\ell}^{-}\\!(x)$
$\displaystyle=h_{\ell}^{(1)}\\!(x)\,\mathcal{J}_{\ell}^{\ell+1}(x;0,1)+\frac{2}{2\ell+1}\bigg{[}\ell\int_{0}^{1}\\!\\!du\,u^{\ell+1}\,h_{\ell}^{(1)}\\!(xu)\,\mathcal{J}_{\ell}^{\ell-1}(x;0,u)$
(A1)
$\displaystyle\qquad+\ell\int_{0}^{1}\\!\\!du\,u^{\ell+1}\,j_{\ell}(xu)\,\mathcal{H}_{\ell}^{\ell-1}(x;u,1)-\big{(}\ell+1\big{)}\,\mathcal{J}_{\ell}^{\ell+1}(x;0,1)\,\mathcal{H}_{\ell}^{-\ell-2}(x;1,\infty)\bigg{]}$
$\displaystyle\frac{1}{x}\,{\cal I}_{\ell}^{+}\\!(x)$
$\displaystyle=j_{\ell}(x)\,\mathcal{H}_{\ell}^{-\ell}(x;1,\infty)+\frac{2}{2\ell+1}\bigg{[}\ell\,\mathcal{J}_{\ell}^{\ell-1}(x;0,1)\,\mathcal{H}_{\ell}^{-\ell}(x;1,\infty)$
(A2)
$\displaystyle\qquad-\big{(}\ell+1\big{)}\int_{1}^{\infty}\\!\\!du\,u^{-\ell}\,h_{\ell}^{(1)}\\!(xu)\,\mathcal{J}_{\ell}^{-\ell-2}(x;1,u)-\big{(}\ell+1\big{)}\int_{1}^{\infty}\\!\\!du\,u^{-\ell}\,j_{\ell}(xu)\,\mathcal{H}_{\ell}^{-\ell-2}(x;u,\infty)\bigg{]}\;,$

which encode the contribution from the inner (-) and outer (+) regions. Here,
$\mathcal{J}_{\ell}^{n}(x;a,b)$ and $\mathcal{H}_{\ell}^{n}(x;a,b)$ are
defined in Eq. (45).

To evaluate $I_{\ell}^{\pm}(x)$ for $x=0$ (which arises whenever $m=0$ or
$\mathcal{M}=0$), we take the limit $x\to 0$ in the above expressions. This
step is relatively straightforward except for the last term in the right-hand
side of Eq. (A2). To handle it, we note that $u^{-\ell-2}h_{\ell}^{(1)}(u)\sim
u^{-2\ell-2}$ for $u\ll 1$ and write

$x\int_{1}^{\infty}\\!\\!du\,u^{-\ell}\,j_{\ell}(xu)\,\mathcal{H}_{\ell}^{-\ell-2}(x;u,\infty)=x^{2\ell+1}\int_{x}^{\infty}\\!\\!ds\,s^{-\ell}j_{\ell}(s)\,s^{-2\ell-2}\bigg{(}s^{2\ell+2}\int_{s}^{\infty}\\!\\!dt\,t^{-\ell-2}\,h_{\ell}^{(1}\\!(t)\bigg{)}\;.$
(A3)

As $x$ approaches zero, the integral in the brackets converges to a constant
plus subleading corrections in $s$. The constant is given by

$\lim_{s\to
0}\,s^{2\ell+1}\int_{s}^{\infty}\\!\\!dt\,t^{-\ell-2}h_{\ell}^{(1)}\\!(t)=-\frac{i(2\ell-1)!!}{2(\ell+1)}\;.$
(A4)

As a result, we have

$\displaystyle\lim_{x\to
0}\,x\int_{1}^{\infty}\\!\\!du\,u^{-\ell}\,j_{\ell}(xu)\,\mathcal{H}_{\ell}^{-\ell-2}(x;u,\infty)$
$\displaystyle=\lim_{x\to
0}\,x^{2\ell+1}\int_{x}^{\infty}\\!\\!ds\,s^{-\ell}j_{\ell}(s)\,s^{-2\ell-2}\bigg{(}-\frac{i(2\ell-1)!!}{2(\ell+1)}\bigg{)}$
(A5) $\displaystyle=-\frac{i}{2(\ell+1)(2\ell+1)^{2}}\;.$

The computation of the other terms in the limit $x\to 0$ proceeds analogously.
Adding up all the contributions leads to Eq. (57).

For a truncated $r^{-2}$ isothermal profile, the (complex) function
$\mathcal{A}(x)$ which appears in the shape factor $\mathcal{F}(\mathcal{M})$,
Eq. (74), reads

$\displaystyle\mathcal{A}(x)$
$\displaystyle=ix\left(\int_{r_{-}/a}^{1}\\!\\!du\,u^{3}\,j_{2}(xu)\,h_{2}^{(1)}(x)+\int_{1}^{r_{+}/a}\\!\\!du\,u^{3}\,j_{2}(x)\,h_{2}^{(1)}(xu)\right)$
(A6)
$\displaystyle\quad+\frac{2ix}{5}\int_{r_{-}/a}^{1}\\!\\!du\,u^{3}\bigg{[}2h_{2}^{(1)}(xu)\,\mathcal{J}_{2}^{1}(x;r_{-}/a,u)+2j_{2}(xu)\,\mathcal{H}_{2}^{1}(x;u,1)-3j_{2}(xu)\,\mathcal{H}_{2}^{-4}(x;1,r_{+}/a)\bigg{]}$
$\displaystyle\quad+\frac{2ix}{5}\int_{1}^{r_{+}/a}\\!\\!du\,u^{3}\bigg{[}2h_{2}^{(1)}(xu)\,\mathcal{J}_{2}^{1}(x;r_{-}/a,1)-3h_{2}^{(1)}(xu)\,\mathcal{J}_{2}^{-4}(x;1,u)-3j_{2}(xu)\,\mathcal{H}_{2}^{-4}(x;u,r_{+}/a)\bigg{]}$

where $r_{-}$ and $r_{+}$ define the radial domain over which the acoustic
wake contribute to the GW emission. For a uniform sphere of density
$\bar{\rho}$, this becomes

$\mathcal{A}(x)=\left(\frac{2\pi
G\bar{\rho}a^{2}}{c_{s}^{2}}\right)ix\bigg{[}h_{2}^{(1)}(x)\,\mathcal{J}_{2}^{4}(x;r_{-},1)+j_{2}(x)\,\mathcal{H}_{2}^{4}(x;1,r_{+}/a)\bigg{]}\;.$
(A7)

For the uniform density $\bar{\rho}=G^{-1}(c_{s}/r_{+})^{2}$ assumed in §4,
the prefactor of the last expression becomes $2\pi(a/r_{+})^{2}$.

## References

  * Andrade-Santos et al. (2013) Andrade-Santos, F., Nulsen, P. E. J., Kraft, R. P., et al. 2013, Astrophys. J., 766, 107
  * Annulli et al. (2020) Annulli, L., Cardoso, V., & Vicente, R. 2020, Phys. Rev. D, 102, 063022
  * Banik & van den Bosch (2021) Banik, U., & van den Bosch, F. C. 2021, Astrophys. J., 912, 43
  * Bar et al. (2021) Bar, N., Blas, D., Blum, K., & Kim, H. 2021, Phys. Rev. D, 104, 043021
  * Bar-Or et al. (2019) Bar-Or, B., Fouvry, J.-B., & Tremaine, S. 2019, Astrophys. J., 871, 28
  * Bekenstein & Maoz (1992) Bekenstein, J. D., & Maoz, E. 1992, Astrophys. J., 390, 79
  * Berezhiani et al. (2019) Berezhiani, L., Elder, B., & Khoury, J. 2019, JCAP , 2019, 074
  * Binney & Tremaine (1987) Binney, J., & Tremaine, S. 1987, Galactic dynamics
  * Boudon et al. (2024) Boudon, A., Brax, P., Valageas, P., & Wong, L. K. 2024, Phys. Rev. D, 109, 043504
  * Boylan-Kolchin et al. (2008) Boylan-Kolchin, M., Ma, C.-P., & Quataert, E. 2008, Mon. Not. R. Astron. Soc., 383, 93
  * Buehler & Desjacques (2023) Buehler, R., & Desjacques, V. 2023, Phys. Rev. D, 107, 023516
  * Buehler et al. (2024) Buehler, R., Kolyada, R., & Desjacques, V. 2024, Mon. Not. R. Astron. Soc., 530, 1151
  * Chandrasekhar (1939) Chandrasekhar, S. 1939, An introduction to the study of stellar structure
  * Chandrasekhar (1943) —. 1943, Astrophys. J., 97, 255
  * Chavanis (2021a) Chavanis, P.-H. 2021a, European Physical Journal Plus, 136, 703
  * Chavanis (2021b) —. 2021b, European Physical Journal Plus, 136, 703
  * Cole et al. (2000) Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, Mon. Not. R. Astron. Soc., 319, 168
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, Mon. Not. R. Astron. Soc., 365, 11
  * Desjacques et al. (2022) Desjacques, V., Nusser, A., & Bühler, R. 2022, Astrophys. J., 928, 64
  * Dokuchaev (1964) Dokuchaev, V. P. 1964, Soviet Astron., 8, 23
  * Fabian et al. (2003) Fabian, A. C., Sanders, J. S., Allen, S. W., et al. 2003, Mon. Not. R. Astron. Soc., 344, L43
  * Foote et al. (2023) Foote, H. R., Besla, G., Mocz, P., et al. 2023, Astrophys. J., 954, 163
  * Ginat et al. (2023) Ginat, Y. B., Panamarev, T., Kocsis, B., & Perets, H. B. 2023, Mon. Not. R. Astron. Soc., 525, 4202
  * Goldreich et al. (2004) Goldreich, P., Lithwick, Y., & Sari, R. 2004, Ann. Rev. Astron. Astrophys., 42, 549
  * Hartman et al. (2021) Hartman, S. T. H., Winther, H. A., & Mota, D. F. 2021, Astron. Astrophys., 647, A70
  * Hui et al. (2017) Hui, L., Ostriker, J. P., Tremaine, S., & et al. 2017, Phys. Rev., D95, 043541
  * Ibata & Lewis (1998) Ibata, R. A., & Lewis, G. F. 1998, Astrophys. J., 500, 575
  * Just & Kegel (1990) Just, A., & Kegel, W. H. 1990, Astron. Astrophys., 232, 447
  * Kalnajs (1971) Kalnajs, A. J. 1971, Astrophys. J., 166, 275
  * Kalnajs (1972) —. 1972, Astrophys. J., 175, 63
  * Katz et al. (2019) Katz, A., Kurkela, A., & Soloviev, A. 2019, JCAP, 08, 017
  * Kauffmann et al. (1993) Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, Mon. Not. R. Astron. Soc., 264, 201
  * Kaur & Stone (2022) Kaur, K., & Stone, N. C. 2022, Mon. Not. R. Astron. Soc., 515, 407
  * Kim & Kim (2007) Kim, H., & Kim, W.-T. 2007, Astrophys. J., 665, 432
  * Kim et al. (2008) Kim, H., Kim, W.-T., & Sánchez-Salcedo, F. J. 2008, Astrophys. J. Lett., 679, L33
  * Kipper et al. (2023) Kipper, R., Tenjes, P., Benito, M., et al. 2023, Astron. Astrophys., 680, A91
  * Landau & Lifshitz (1975) Landau, L. D., & Lifshitz, E. M. 1975, The classical theory of fields (Elsevier Science & Technology)
  * Lee & Stahler (2011) Lee, A. T., & Stahler, S. W. 2011, Mon. Not. R. Astron. Soc., 416, 3177
  * Lynden-Bell & Kalnajs (1972) Lynden-Bell, D., & Kalnajs, A. J. 1972, Mon. Not. R. Astron. Soc., 157, 1
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, Ann. Rev. Astron. Astrophys., 45, 565
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution
  * Murai & Fujimoto (1980) Murai, T., & Fujimoto, M. 1980, PASJ, 32, 581
  * Narayan & Medvedev (2001) Narayan, R., & Medvedev, M. V. 2001, Astrophys. J. Lett., 562, L129
  * Nelson & Tremaine (1999) Nelson, R. W., & Tremaine, S. 1999, Mon. Not. R. Astron. Soc., 306, 1
  * O’Neill et al. (2024) O’Neill, D., D’Orazio, D. J., Samsing, J., & Pessah, M. E. 2024, arXiv e-prints, arXiv:2401.16166
  * Ostriker (1999) Ostriker, E. C. 1999, Astrophys. J., 513, 252
  * Paczynski (1976) Paczynski, B. 1976, in Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, Vol. 73, 75
  * Palmer & Papaloizou (1985) Palmer, P. L., & Papaloizou, J. 1985, Mon. Not. R. Astron. Soc., 215, 691
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
  * Rephaeli & Salpeter (1980) Rephaeli, Y., & Salpeter, E. E. 1980, Astrophys. J., 240, 20
  * Ruderman & Spiegel (1971) Ruderman, M. A., & Spiegel, E. A. 1971, Astrophys. J., 165, 1
  * Sánchez-Salcedo (2019) Sánchez-Salcedo, F. J. 2019, Astrophys. J., 885, 152
  * Sánchez-Salcedo & Brandenburg (2001) Sánchez-Salcedo, F. J., & Brandenburg, A. 2001, Mon. Not. R. Astron. Soc., 322, 67
  * Sanders et al. (2013) Sanders, J. S., Fabian, A. C., Churazov, E., et al. 2013, Science, 341, 1365
  * Schekochihin & Cowley (2006) Schekochihin, A. A., & Cowley, S. C. 2006, Physics of Plasmas, 13, 056501
  * Sijacki & Springel (2006) Sijacki, D., & Springel, V. 2006, Mon. Not. R. Astron. Soc., 371, 1025
  * Somerville & Primack (1999) Somerville, R. S., & Primack, J. R. 1999, Mon. Not. R. Astron. Soc., 310, 1087
  * Tomaselli et al. (2023) Tomaselli, G. M., Spieksma, T. F. M., & Bertone, G. 2023, JCAP , 2023, 070
  * Traykova et al. (2021) Traykova, D., Clough, K., Helfer, T., et al. 2021, Phys. Rev. D, 104, 103014
  * Traykova et al. (2023) Traykova, D., Vicente, R., Clough, K., et al. 2023, arXiv e-prints, arXiv:2305.10492
  * Tremaine & Weinberg (1984) Tremaine, S., & Weinberg, M. D. 1984, Mon. Not. R. Astron. Soc., 209, 729
  * Tremaine et al. (1975) Tremaine, S. D., Ostriker, J. P., & Spitzer, L., J. 1975, Astrophys. J., 196, 407
  * van den Bosch et al. (1999) van den Bosch, F. C., Lewis, G. F., Lake, G., & Stadel, J. 1999, Astrophys. J., 515, 50
  * Vicente & Cardoso (2022) Vicente, R., & Cardoso, V. 2022, arXiv e-prints, arXiv:2201.08854
  * Vicente et al. (2019) Vicente, R., Cardoso, V., & Zilhão, M. 2019, Mon. Not. R. Astron. Soc., 489, 5424
  * Weinberg (1986) Weinberg, M. D. 1986, Astrophys. J., 300, 93
  * Weinberg (1989) —. 1989, Mon. Not. R. Astron. Soc., 239, 549