11institutetext:  Universitäts-Sternwarte, Fakultät für Physik, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, D-81679 München, Germany
11email<EMAIL_ADDRESS>22institutetext: Excellence Cluster ORIGINS,
Boltzmannstrasse 2, D-85748 Garching, Germany 33institutetext: Institute for
Theoretical Physics, Goethe University, 60438 Frankfurt am Main Germany
33email<EMAIL_ADDRESS>
# Dynamical friction from self-interacting dark matter

Moritz S. Fischer and Laura Sagunski Dynamical friction from self-interacting
dark matterDynamical friction from self-interacting dark matterDynamical
friction from self-interacting dark matterDynamical friction from self-
interacting dark matterDynamical friction from self-interacting dark
matterDynamical friction from self-interacting dark matter

(Received XX Month, 20XX / Accepted XX Month, 20XX)

###### Abstract

Context. Merging compact objects such as binary black holes provide a
promising probe for the physics of dark matter (DM). The gravitational waves
emitted during inspiral potentially allow to detect DM spikes around black
holes. This is because the dynamical friction force experienced by the
inspiraling black hole alters the orbital period and thus the gravitational
wave signal.

Aims. The dynamical friction arising from DM can potentially differ from the
collisionless case when DM is subject to self-interactions. This paper aims to
understand how self-interactions impact dynamical friction.

Methods. To study the dynamical friction force, we use idealized $N$-body
simulations, where we include self-interacting dark matter.

Results. We find that the dynamical friction force for inspiraling black holes
would be typically enhanced by DM self-interactions compared to a
collisionless medium (ignoring differences in the DM density). At lower
velocities below the sound speed, we find that the dynamical friction force
can be reduced by the presence of self-interactions.

Conclusions. DM self-interactions have a significant effect on the dynamical
friction for black hole mergers. Assuming the Chandrasekhar formula may
underpredict the deceleration due to dynamical friction.

###### Key Words.:

methods: numerical –- dark matter – black hole physics – gravitational waves

##  1 Introduction

Since the first direct detection of gravitational waves (GWs) from merging
stellar-mass black holes (BHs) and neutron stars with the LIGO detectors
(Abbott et al., 2016, 2017a, 2017b, 2019), 100 years after the formulation of
Einstein’s theory of General Relativity (Einstein, 1916) and about 40 years of
the first indirect evidence (Hulse & Taylor, 1975; Taylor et al., 1976;
Weisberg et al., 1981), GW astronomy has become an important and very
promising field to probe fundamental physics. The GWs emitted by binary
mergers do not only allow for testing general relativity at unprecedented
levels, but also to probe other physical processes. This in particular
includes the physics of neutron stars. Binary mergers of those objects could
allow probing the equation of state of ultrahigh-density matter (e.g. Baiotti,
2022). In the future, GW detectors will become sensitive to sub-Hertz GW
frequencies, thanks to space-based GW observatories such as LISA (Amaro-Seoane
et al., 2017) and Taiji (Hu & Wu, 2017). This will allow us to observe the
evolution of intermediate- and extreme-mass ratio inspirals (I/EMRIs), in
which a small stellar-mass compact object inspirals into a much larger
intermediate/supermassive black hole, over long times and provide a sensitive
probe of BH environments.

BHs might be “dressed” with DM, i.e., surrounded by a dense DM spike (e.g.
Gondolo & Silk, 1999; Bertone et al., 2005; Bertone, 2024). The first robust
evidence for an observation of a DM spike around a BH is claimed to be found
by Chan & Lee (2024), see also the work by Alachkar et al. (2023). They
studied the supermassive black hole binary system OJ 287 with outburst
observations covering a time span of 104 years. The presence of a DM spike can
impact the orbit of an inspiralling BH and make the GW signal sensitive to DM
physics. This could in principle be detectable, e.g., for M87* (Daghigh &
Kunstatter, 2024). The inspiralling BH is decelerated by the dynamical
friction force, which depends on the density and the velocity dispersion of
the DM spike (Chandrasekhar, 1943). In consequence, the orbital decay of an
inspiralling BH allows probing properties of the DM spike and thus the
particle nature of DM.

Potentially, such DM spikes could lead to a faster decay of the orbit at small
separations ($\approx 1$ pc). This would provide an explanation for the final
parsec problem, i.e., how binary BHs can lose enough energy and angular
momentum to reach separations small enough that the energy loss due to the
emission of GWs is sufficient to make them merge within a sufficiently small
timescale (e.g. Milosavljević & Merritt, 2003). In this context, Alonso-
Álvarez et al. (2024) claim that collisionless DM is not sufficient to make
the BH orbit decay fast enough because the DM spike would be destroyed by the
friction energy injected from the BHs. Instead, to solve the final parsec
problem, DM self-interactions would be needed to recreate the spike on a short
enough timescale.

Self-interacting dark matter (SIDM) is a class of particle physics models
assuming that DM has self-interactions beyond gravity. Initially, SIDM was
proposed to solve problems on galactic scales, namely reducing the abundance
of Milky Way satellites and the density in the centre of DM halos (Spergel &
Steinhardt, 2000). DM self-interactions have been studied in systems covering
many orders of magnitude in mass, ranging from dwarf galaxies to galaxy
clusters (e.g. Vogelsberger et al., 2014; Fitts et al., 2019; Robertson et
al., 2019; Zeng et al., 2022; Ragagnin et al., 2024; Sabarish et al., 2024).
Observations of those systems have been used to constrain the strength of the
self-interactions (e.g. Sagunski et al., 2021; Despali et al., 2022; Kong et
al., 2024; Yang et al., 2024; Zhang et al., 2024). For a review on SIDM, we
refer the reader to Tulin & Yu (2018); Adhikari et al. (2022).

Potentially, the range of systems used to constrain SIDM can be extended by
merging binary BHs. Based on pulsar timing array observations (Agazie et al.,
2023; Reardon et al., 2023; Antoniadis et al., 2023), it was claimed that a
self-interaction cross-section of order
$0.5$–$5.0\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ would be needed to solve the
final parsec problem (Alonso-Álvarez et al., 2024). We note that this claim
has been made based on the assumption of the density spike being in the long-
mean-free-path regime (see Section 2.3). Interestingly, those binary systems
could be a very sensitive probe for DM scattering, as we discuss in this
paper.

When DM is self-interacting, the density profile of the DM spike may change
compared to collisionless matter. For SIDM, the slope of the density profile
of the DM spike could be smaller than in the case of collisionless DM (Shapiro
& Paschalidis, 2014). The lower density would allow relaxing constraints from
DM annihilation (Alvarez & Yu, 2021). Importantly, this also affects the
dynamical friction, as it depends on the density of the DM spike, and thus
also the GW signal (e.g., Macedo et al., 2013; Yue & Han, 2018; Becker et al.,
2022). However, given the properties of the density spike, one may still need
to calibrate the strength of the dynamical friction with simulations. Kavanagh
et al. (2020, 2024) run $N$-body simulations of collisionless DM to measure
the Coulomb logarithm for a BH orbiting in a DM spike.

Several studies of dynamical friction for various DM models have been
conducted. This includes ultralight DM (Wang & Easther, 2022; Boey et al.,
2024), also with self-interactions (Glennon et al., 2024), self-interacting
scalar DM (Boudon et al., 2022, 2023a, 2023b), super-fluid DM (Berezhiani et
al., 2023), fuzzy DM (Lancaster et al., 2020). In this paper, we investigate
the strength of the dynamical friction for a medium consisting of SIDM.

Often the binary BHs have been described in the Newtonian limit, however,
depending on their evolution stage this may not be valid. Studies beyond
Newtonian gravity have been undertaken by various authors (e.g., Traykova et
al., 2023; Montalvo et al., 2024; Mukherjee et al., 2024). For this study, we
focus on the Newtonian regime, i.e., when the separation of the two BHs is
sufficiently large. This allows us to build upon earlier work and simplifies
the numerical description compared to the relativistic case. Extending our
study to the relativistic case is left for future work.

The rest of this paper is structured as follows: Section 2 deals with the
analytic description of dynamical friction and discusses it in the context of
a BH moving through a DM spike. In Section 3, we describe the numerical setup
that we use to study dynamical friction in the context of DM self-
interactions. This includes the simulation code and the initial conditions and
parameters for the simulations. It follows the analysis of our simulations and
the results (Section 4). In Section 5, we discuss the limitations of this work
as well as directions for further research. Finally, we summarise and conclude
in Section 6. Additional information is provided in the appendix.

##  2 Analytic Estimates

In this section, we first describe the properties of DM spikes around BHs
relevant to our study. Secondly, we review analytic estimates of the dynamical
friction force for collisionless and gaseous media. Last, we discuss how to
describe the strength of self-interactions and their impact on the dynamical
friction arising from DM spikes acting on inspiralling BHs.

###  2.1 DM spike

Figure 1: The density and one-dimensional velocity dispersion are shown for a
DM density spike as a function of radius according to Eq. (2). We assume
$r_{\mathrm{sp}}=1\mu\mathrm{pc}$, $\rho_{\mathrm{sp}}=5.45\times
10^{-3}\,\mathrm{M_{\odot}}\,\mu\mathrm{pc}^{-3}$, $\alpha_{\mathrm{sp}}=7/3$,
and $M_{\mathrm{BH}}=10^{3}\,\mathrm{M}_{\odot}$. In addition, we indicate the
Schwarzschild radius $R_{s}=2\mathrm{G}M/c^{2}$. For computing the velocity
dispersion, we follow Kremer (2022). Here, we assume an isotropic velocity
distribution. In contrast to Kremer (2022), we solve their eq. 27 including
higher-order terms instead of eq. 28. Moreover, we show the results for
$m/(k_{\mathrm{B}}T)\approx 6.53\times
10^{3}\,\mathrm{s}^{2}\,\mathrm{km}^{-2}$, with $m$ being the rest mass of a
physical DM particle, $k_{\mathrm{B}}$ the Boltzmann constant and $T$ the
absolute temperature of the system.

In the vicinity of a BH the DM density can potentially reach high values due
to the gravitational influence of the BH. Such a DM spike is commonly
parameterised by a power law,

$\rho_{\mathrm{DM}}(r)=\rho_{6}\left(\frac{r}{r_{6}}\right)^{-\alpha_{\mathrm{sp}}}\,.$
(1)

The radius $r_{6}=10^{-6}\,\mathrm{pc}$ is a reference radius at which the
density is given by $\rho_{6}$. Note, with this choice, we follow Coogan et
al. (2022).

When taking relativistic effects into account (Sadeghian et al., 2013) the
spike profile becomes,

$\rho_{\mathrm{DM}}(r)=\begin{cases}\,0&\textnormal{for}\quad r\leq
2\,R_{s}\\\ \,\rho_{\rm
sp}\left(1-\frac{2R_{s}}{r}\right)^{3}\left(\frac{r}{r_{\rm
sp}}\right)^{-\alpha_{\mathrm{sp}}}&\textnormal{for}\quad 2\,R_{s}<r\leq
r_{\rm sp}\\\ \,\frac{\rho_{s}\,r_{s}}{r}&\textnormal{for}\quad r_{\rm
sp}<r\ll r_{s}\\\ \end{cases}\,.$ (2)

For the relativistic corrections on the properties of the DM spike, we refer
the reader to Eda et al. (2015); Tang et al. (2021); Capozziello et al.
(2023); John et al. (2023). In Fig. 1 we show the corresponding spike density
and velocity dispersion as a function of radius. We note that the density
profile of the DM around the BH could be even more complicated when
considering a realistic formation scenario of the BH (Bertone et al., 2024).

As we will see in the next subsection the Mach number of the secondary BH
moving through the density spike plays a crucial role in determining the
dynamical friction force acting on it. The Mach number,
$\mathcal{M}=v_{\mathrm{BH}}/c_{\mathrm{s}}$, is given by the velocity,
$v_{\mathrm{BH}}$, of the BH and the sound speed, $c_{\mathrm{s}}$. The local
speed of sound depends on the adiabatic index $\gamma$, and is given by

$c_{\mathrm{s}}=\sqrt{\gamma}\,\sigma_{\mathrm{DM}}\,.$ (3)

To estimate $c_{\mathrm{s}}$, we assume for simplicity that the local velocity
distribution obeys a Maxwell-Boltzmann distribution with a one-dimensional
velocity dispersion $\sigma_{\mathrm{DM}}$. To estimate the Mach number,
$\mathcal{M}$, we assume that the DM spike at sufficiently large radii is well
described by Eq. 1.

Furthermore, we assume that the gravitational potential is dominated by the
mass, $M_{\mathrm{BH}}$ of the primary (more massive) BH, and that the
contribution of the DM spike and of the less massive inspiralling BH are
negligible. Under these assumptions, the one-dimensional velocity dispersion
of the DM is

$\sigma_{\mathrm{DM}}=\sqrt{\frac{\mathrm{G}\,M_{\mathrm{BH}}}{1+\alpha_{\mathrm{sp}}}\frac{1}{r}}\,.$
(4)

Assuming that the BH is moving on a circular orbit
($v_{\mathrm{circ}}=\sqrt{\mathrm{G}\,M(<r)/r}$), it follows that the Mach
number is

$\mathcal{M}=\sqrt{\frac{1+\alpha_{\mathrm{sp}}}{\gamma}}\,.$ (5)

The adiabatic growth model for a DM spike assuming collisionless DM suggests,
for an NFW halo, a spike index of $\alpha_{\mathrm{sp}}=7/3$ (Gondolo & Silk,
1999; Fields et al., 2014; Eda et al., 2015). For SIDM, Shapiro & Paschalidis
(2014) derived the DM spike to be $\alpha_{\mathrm{sp}}=7/4$. We note that
this is based on the assumption that the self-interaction cross-section has
the velocity dependence $\sigma\propto v^{-4}$. Moreover, Shapiro &
Paschalidis (2014) assumed that the spike resides in the long-mean-free-path
regime, in Section 2.3, we explain that this may often not be the case. This
implies that the motion of the inspiralling BH is supersonic when being on a
circular orbit (see Eq. 5) given a spike index of $\alpha_{\mathrm{sp}}=7/3$
or $\alpha_{\mathrm{sp}}=7/4$. Even for a smaller value of
$\alpha_{\mathrm{sp}}$, it is not feasible to reach the subsonic regime on a
circular orbit within the range typically considered for the spike index
($1<\alpha_{\mathrm{sp}}<3$). However, a lower value for
$\alpha_{\mathrm{sp}}$ might be possible depending on the formation scenario
(Ullio et al., 2001).

###  2.2 Dynamical friction

The dynamical friction force for a perturber moving through a collisionless
medium of constant density was derived by Chandrasekhar (1943). The strongly
collisional case with an object moving through a gaseous medium was studied
analytically by several authors (e.g., Rephaeli & Salpeter, 1980; Just et al.,
1986; Just & Kegel, 1990). Ostriker (1999) derived an analytic description for
this regime. Following her, the dynamical friction force can be expressed as

$F_{\mathrm{df}}=-4\pi\,\rho_{\mathrm{DM}}\left(\frac{\mathrm{G}\,m_{\mathrm{BH}}}{v_{\mathrm{BH}}}\right)^{2}\,I\,.$
(6)

Here, $m_{\mathrm{BH}}$, is the mass of the perturber, in our case a BH. Its
velocity is given by $v_{\mathrm{BH}}$ and the density of the perturbed
medium, in our case DM, is denoted by $\rho_{\mathrm{DM}}$. The factor $I$
depends besides other things such as the Mach number, $\mathcal{M}$, on the
collisional nature of the perturbed medium and can be specified for the
different cases. In the collisionless case, $I$, is given by

$I_{\mathrm{collisionless}}=\ln{\left(\frac{r_{\mathrm{max}}}{r_{\mathrm{min}}}\right)}\left[\mathrm{erf}(X)-\frac{2X}{\sqrt{\pi}}\,e^{-X^{2}}\right]\\\
\textnormal{with}\quad
X\equiv\frac{v_{\mathrm{BH}}}{\sigma_{\mathrm{DM}}\sqrt{2}}\,.$ (7)

The effective size of the BH is indicated by $r_{\mathrm{min}}$, and
$r_{\mathrm{max}}$ denotes the effective size of the surrounding medium.
Where, $r_{\mathrm{min}}$, can be understood as the $90\degree$ deflection
radius and $r_{\mathrm{max}}$ as the maximum impact parameter (Binney &
Tremaine, 2008). We note, that while Chandrasekhar’s formula given by Eq. 7
approximates the dynamical friction well, corrections to the Coulomb
logarithm, $\ln(r_{\mathrm{max}}/r_{\mathrm{min}})$, can be necessary (Just &
Peñarrubia, 2005; Just et al., 2010), e.g., for velocity distributions
deviating from a Maxwell–Boltzmann distribution. For corrections arising from
fast moving DM particles on the Chandrasekhar formula, we refer the reader to
the work by Dosopoulou (2023) on the dynamical friction in DM spikes.

For gaseous media, one distinguishes the subsonic case where the BH is moving
with a velocity $v$ that is smaller than the speed of sound $c_{\mathrm{s}}$
and the supersonic case where it is faster than $c_{\mathrm{s}}$. In the first
one, the Mach number, $\mathcal{M}=v/c_{\mathrm{s}}$ is smaller than unity and
$I$ is given as

$I_{\mathrm{subsonic}}=\frac{1}{2}\ln{\left(\frac{1+\mathcal{M}}{1-\mathcal{M}}\right)}-\mathcal{M}\,.$
(8)

Note, this equation is based on the assumption that
$(c_{s}-v_{\mathrm{BH}})\,t$ exceeds the effective size, $r_{\mathrm{min}}$,
of the perturber and that $(c_{s}+v_{\mathrm{BH}})\,t$ is smaller than the
effective size, $r_{\mathrm{max}}$, of the surrounding medium. The time, $t$,
indicates for how long the BH is perturbing the DM distribution. In the
supersonic case with $\mathcal{M}>1$, we have

$I_{\mathrm{supersonic}}=\frac{1}{2}\,\ln{\left(1-\frac{1}{\mathcal{M}^{2}}\right)}+\ln{\left(\frac{v_{\mathrm{BH}}\,t}{r_{\mathrm{min}}}\right)}\,.$
(9)

This implicitly assumes $(v_{\mathrm{BH}}-c_{s})\,t>r_{\mathrm{min}}$ and
$(v_{\mathrm{BH}}+c_{s})\,t<r_{\mathrm{max}}$.

In Fig. 2, we plot the drag force computed according to Ostriker (1999) for
our test setup used in Section 4. We can see that the dynamical friction for
velocities slightly above the sound speed is becoming fairly strong in the
gaseous case and sharply drops when the velocity falls below the speed of
sound.

In the previous subsection, we have seen that the BH would typically move with
a supersonic velocity. In consequence, the inspiralling BH will experience
enhanced dynamical friction if the DM is not collisionless. However, if the BH
is on an eccentric orbit, it might partially, i.e., when it is close to its
apocentre, move with a subsonic velocity. We also note, that dynamical
friction can contribute to a circularization of the orbit (Becker et al.,
2022; Karydas et al., 2024).

In general, the dynamical friction acting on inspiralling BHs may not only
arise from the DM spike but also from the accretion disk. The latter can also
be described with the Ostriker formula (e.g., Szölgyén et al., 2022; Becker &
Sagunski, 2023).

Figure 2: The dynamical friction force is shown as a function of the Mach
number, i.e. the velocity in units of the sound speed. The dashed lines
indicated the dynamical friction for a collisionless medium according to the
Chandrasekhar formula (Chandrasekhar, 1943) and the solid lines are for a
gaseous medium as given by Ostriker (1999). Note, for the plot, we assume
$r_{\mathrm{max}}=\mathcal{M}\,c_{s}\,t$ and
$\ln(r_{\mathrm{max}}/r_{\mathrm{min}})=\ln(\mathcal{M})+C$. The line colours
correspond to different values of $C\in\\{4.0,6.0,8.0,10.0,12.0\\}$. In
addition, the red dotted lines indicate the velocity of a secondary BH moving
on a circular orbit for different adiabatic indices.

###  2.3 The strength of the self-interactions

A remaining question is how close a given cross-section would be to the
collisionless or gaseous regime. This means we need to understand how much the
trajectory of the DM particles is impacted by gravity compared to the
scattering. A common measure for this in the context of SIDM halos is the
Knudsen number (e.g. Koda & Shapiro, 2011). It is given by

$\mathrm{Kn}=\sqrt{\frac{4\pi\mathrm{G}}{\rho\sigma_{\mathrm{DM}}^{2}}}\left(\frac{\sigma_{\mathrm{V}}}{m}\right)^{-1}\,.$
(10)

This is comparing the ratio of the mean free path to the Jeans length. The
latter is a measure for how susceptible the DM medium is to gravitational
collapse. We note that the Knudsen number, Kn, is independent of the mass of
the perturber, i.e., the BH mass. However, the strength of the dynamical
friction depends on it. For Eq. 10 we used the viscosity cross-section per
mass, $\sigma_{\mathrm{V}}/m$, for defining the Knudsen number. It is given by

$\sigma_{\mathrm{V}}=3\pi\int_{-1}^{1}\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}\sin^{2}\theta\,\mathrm{d}\cos\theta\,.$
(11)

Here, we follow the normalisation of Yang & Yu (2022) implying that viscosity
and total cross-section have equal values for isotropic scattering
($\left.\sigma_{\mathrm{V}}\right|_{\rm
iso}=\left.\sigma_{\mathrm{tot}}\right|_{\rm iso}$). The viscosity cross-
section has the advantage that it often provides a good matching between
different angular dependencies (e.g., Yang & Yu, 2022; Colquhoun et al., 2021)
and is thus a good characteristic for the strength of the scattering. If the
Knudsen number is large, i.e. $\mathrm{Kn}\gg 1$, a system is considered to be
in the long-mean-free-path regime, i.e., gravity dominates over the
scattering. In contrast, if $\mathrm{Kn}\ll 1$, the system would be in the
short-mean-free-path regime and scattering dominates over gravity. So we
expect that the system is close to the collisionless case when the Knudsen
number is larger, and for small numbers falls into the gaseous regime.

For the DM spike, the Knudsen number becomes smaller when getting closer to
the primary BH as density and velocity dispersion increase. This implies that
the behaviour of the dynamical friction could transition from the
collisionless regime to the gaseous one while the secondary BH is
inspiralling. However, it depends on the strength of the DM self-interactions.

Last, we note that also in the case of a primary BH hosting a DM spike, the
heat conduction due to self-interactions being in the long-mean-free-path or
short-mean-free-path regime is independent of the mass of the primary BH. In
other words, the Knudsen number does not depend on the mass of the primary BH.
We note that this is in contradiction to the assumption by Shapiro &
Paschalidis (2014); Shapiro (2018). In consequence, a DM spike given a cross-
section within the typically studied range of $\approx$ 0.1–100 cm2g-1 (e.g.
Tulin & Yu, 2018; Adhikari et al., 2022) would be to a large extent rather in
the short-mean-free-path regime than in the long-mean-free-path regime. This
is in line with the work by Pollack et al. (2015), where they studied a DM
model with a fraction of the matter being ultra-strongly self-interacting.
Their heat conductivity does only depend on the density of the ultra-strongly
interacting matter component.

##  3 Numerical Setup

In this section, we first describe our simulation code with its implementation
of DM self-interactions. Next, it follows the setup for our simulations to
understand the role of self-interactions for dynamical friction.

###  3.1 Simulation Code

We use the cosmological $N$-body code OpenGadget3 (see Groth et al., 2023, and
the references therein). The SIDM module of the code has been described by
Fischer et al. (2021a, b, 2022, 2024b). Note, it can simulate both large and
small-angle scatterings. To the first one, we also refer to as rare scattering
since a single scattering event on average leads to a larger energy and
momentum transfer between the particles compared to small-angle scattering,
which must be much more frequent to have a similar effect. For the first one,
the angular dependence is modelled explicitly by drawing random angles with
respect to the differential cross-section, for the second one an effective
description is used based on a drag-like behaviour (see also Kahlhoefer et
al., 2014). The interactions between the numerical particles are modelled on a
pair-wise basis. We compute the scattering probability for the rarely self-
interacting dark matter (rSIDM) and the drag force for the frequently self-
interacting dark matter (fSIDM) using a spline kernel (Monaghan & Lattanzio,
1985). Its size is determined with the next neighbour approach using
$N_{\mathrm{ngb}}=64$ neighbours. This choice is known to provide accurate
results (Fischer et al., 2021a). We employ a separate time step criterion for
SIDM to keep the interaction probability per interaction sufficiently low. The
SIDM module is designed to explicitly conserve linear momentum and energy,
even when executed in parallel.

To model the gaseous regime when self-interactions are very strong, we use
smoothed-particle hydrodynamics (SPH, Gingold & Monaghan, 1977). Compared to
gadget-2 (Springel, 2005), an updated SPH formulation (Beck et al., 2016) with
an updated viscosity scheme (Dolag et al., 2005; Beck et al., 2016) is
employed. We choose a Wendland $C^{6}$ kernel (Dehnen & Aly, 2012) with its
size being determined by the, $N_{\mathrm{ngb}}=230$, next neighbours.

###  3.2 Simulation Setup

To test the dynamical friction, we simulate a BH travelling through a constant
density consisting of DM. The mass density is set to $\rho=4.815\times
10^{-4}\,\mathrm{M_{\odot}}\,{\mu\mathrm{pc}}^{-3}$. The velocities of the DM
particles follow a Maxwell-Boltzmann distribution, and their one-dimensional
velocity dispersion is chosen to be
$\sigma_{\mathrm{1D}}=625.5\,\mathrm{km}\,\mathrm{s}^{-1}$. The BH is
represented by a single particle with a mass of $m_{\mathrm{BH}}=7\times
10^{3}\,\mathrm{M_{\odot}}$. Initially, it is moving with a velocity of
$v_{\mathrm{ini}}=1037\,\mathrm{km}\,\mathrm{s}^{-1}$. These values are
selected because they lead to a strong dynamical friction force, and we expect
to find a significant impact in our simulations. Furthermore, we chose a cubic
simulation volume with a side length of
$l_{\mathrm{box}}=829.5\,\mu\mathrm{pc}$ for most of our simulations but also
use a box length of $l_{\mathrm{box}}=663.6\,\mu\mathrm{pc}$. We vary the box
length to study its impact on the dynamical friction force. In addition,
periodic boundary conditions are employed. The DM is resolved by
$N_{\mathrm{DM}}=5\times 10^{6}$ particles, which leads to a numerical
particle mass of $m_{\mathrm{DM}}=5.5\times 10^{-2}\,\mathrm{M_{\odot}}$. For
testing the accuracy of our simulations, we also perform runs with
$N_{\mathrm{DM}}=1\times 10^{6}$ DM particles, i.e., a five times higher
particle mass. We set the Plummer-equivalent gravitational softening length to
$\epsilon=24\,\mu\mathrm{pc}$. The gravitational forces are computed between
all particles, i.e., we also take the self-gravity of the perturbed medium
into account. In principle, the self-gravity can be relevant for the resulting
dynamical friction (see Just & Kegel, 1990). However, we do not separately
study its contribution.

Note, that the $N$-body representation of the constant density is subject to
fluctuation due to the sampling noise. The fluctuations can collapse under
gravitational influence. We have checked that, for our choice of the box
length and resolution, this does not occur within the time we simulate.

We run simulations with a very anisotropic cross-section by employing the
numerical scheme for frequent self-interactions developed by Fischer et al.
(2021a). Here the strength of the self-interactions is specified in terms of
the momentum transfer cross-section,

$\sigma_{\mathrm{T}}=2\pi\int_{-1}^{1}\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega_{\mathrm{cms}}}\left(1-|\cos\theta_{\mathrm{cms}}|\right)\mathrm{d}\cos\theta_{\mathrm{cms}}\,.$
(12)

We simulate velocity-independent elastic scattering with cross-sections of
$\sigma_{\mathrm{V}}/m\in\\{3\times 10^{-9},3\times 10^{-8},3\times
10^{-7},9\times 10^{-7},3\times
10^{-6}\\}\,(\mathrm{cm}^{2}\,\mathrm{g}^{-1})$.111To make it easier to
compare cross-sections of different angular dependencies, we do not specify
the cross-section in terms of the momentum transfer cross-section but use the
viscosity cross-section (see Eq. 11). In addition, we simulate velocity-
independent isotropic cross-sections with a total cross-section of
$\sigma_{\mathrm{V}}/m\in\\{2\times 10^{-8},3\times 10^{-8},2\times
10^{-7},3\times 10^{-7}\\}\,(\mathrm{cm}^{2}\,\mathrm{g}^{-1})$ to investigate
the role of the angular dependence. We summarise all the cross-sections used
in Tab. 1.

In addition, we also investigate velocity-dependent self-interactions. Here,
we use the implementation by Fischer et al. (2024b). We parameterise the
velocity dependence of the viscosity cross-section as

$\frac{\sigma_{\mathrm{V}}}{m}=\frac{\sigma_{0}}{m}\left(1+\left(\frac{v}{w}\right)^{2}\right)^{-2}\,.$
(13)

For the velocity parameter, $w$, we choose
$w\in\\{300,461.1,683.2,1483.5\\}\,(\mathrm{km}\,\mathrm{s}^{-1})$. Note, the
mean relative velocity of the DM particles in the ICs is $\langle
v\rangle=1411.6\,\mathrm{km}\,\mathrm{s}^{-1}$. For all simulations with a
velocity-dependent cross-section, we choose the value of $\sigma_{0}$ to
result in an effective cross-section (see Yang & Yu, 2022) of
$\sigma_{\mathrm{eff}}/m=3\times 10^{-8}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$. A
summary of all cross-sections is given in Tab. 1.

Note, that all the cross-sections mentioned above are many orders of magnitude
smaller than the tightest constraints on SIDM in the current literature (e.g.
Harvey et al., 2019; Andrade et al., 2021; Sagunski et al., 2021; Eckert, D.
et al., 2022; Gopika & Desai, 2023). Cross-sections as large as the ones
typically discussed in the astrophysical context imply a mean free path
substantially shorter than the kernel size, and may no longer be faithfully
simulated with the standard SIDM schemes (see our Appendix A as well as
Fischer et al., 2024a). In addition, the simulations become more costly with
increasing cross-section as the required time-step decreases. To model the
regime of such a large cross-section, we use SPH. For the ICs, we set the
initial velocity of the SPH particles to zero and choose their internal energy
to correspond to the same velocity dispersion as in the SIDM setup.

Type | $\sigma_{\mathrm{V}}/m$ | $w$
---|---|---
| [cm2g-1] | [km s-1]
Collisionless | 0.0 | –
Frequent | $3\times 10^{-9}$ | –
Rare | $2\times 10^{-8}$ | –
Frequent | $3\times 10^{-8}$ | –
Rare | $3\times 10^{-8}$ | –
Rare | $2\times 10^{-7}$ | –
Frequent | $3\times 10^{-7}$ | –
Rare | $3\times 10^{-7}$ | –
Frequent | $9\times 10^{-7}$ | –
Frequent | $3\times 10^{-6}$ | –
Frequent | $6\times 10^{-5}$ | 300
Rare | $6\times 10^{-5}$ | 300
Frequent | $1.2\times 10^{-5}$ | 461.1
Rare | $1.2\times 10^{-5}$ | 461.1
Frequent | $3\times 10^{-6}$ | 683.2
Rare | $3\times 10^{-6}$ | 683.2
Frequent | $3\times 10^{-7}$ | 1483.5
Rare | $3\times 10^{-7}$ | 1483.5
Fluid | – | –
Table 1: The table gives the different cross-sections we use for the
simulations of our test setup. The first column specifies the type of
scattering, where frequent refers to a very anisotropic cross-section and rare
to an isotropic one. The second row gives the viscosity cross-section (Eq.
11). In the case of the velocity-dependent cross-sections, the second row
corresponds to $\sigma_{0}/m$ of Eq. (13). The last column gives the parameter
$w$ of Eq. (13) for the velocity-dependent cross-sections.

##  4 Test Problems

In this section, we present the results from our simulations following the
setup described in section 3. We begin with the collisionless case and
investigate the impact of variations in the numerical setup. Subsequently, we
study the role of the strength of the self-interactions and its impact on the
density wake. Furthermore, we investigate the matching of different angular
and velocity dependencies onto each other in the case of dynamical friction.

###  4.1 The collisionless case

Figure 3: The velocity and position of the BH in our test setup are shown as a
function of time. All simulations are for collisionless DM. The solid lines
are from the setup with a box length of $l=829.5\mu\mathrm{pc}$. In contrast,
the dashed lines correspond to a setup with a smaller box size of
$l=663.6\mu\mathrm{pc}$. The high-resolution results correspond to
$N_{\mathrm{DM}}=5\times 10^{6}$ DM particles, and the low-resolution ones to
$N_{\mathrm{DM}}=10^{6}$. Note that in the smaller box setup, the low
resolution implies a higher mass resolution than in the large box setup.

In our test setup as described in Section 3, the BH is travelling through an
initially constant density. The BH perturbs the DM and creates a density wake,
leading effectively to a force decelerating the BH. This dynamical friction
force depends on the velocity of the BH relative to the velocity dispersion of
the perturbed DM density (see Fig. 2). In our setup, we first expect the
dynamical friction to become stronger while the BH is slowing down. After the
BH velocity has decreased enough, it enters the phase where the dynamical
friction becomes weaker with decreasing velocity.

In Fig. 3, we show the position and velocity of the BH in our test setup for
collisionless matter. As expected, the velocity of the BH decreases over the
course of the simulation. The results for two different box sizes are shown.
We find that the dynamical friction is reduced for a smaller box size. This is
not surprising, as $r_{\mathrm{max}}$ in the Coulomb logarithm of the
Chandrasekhar formula, Eq. (7), depends for our setup effectively on the box
size because the box size limits the maximal possible impact parameter. The
results agree qualitatively with the Chandrasekhar formula, after an initial
phase during which the density wake forms, i.e., Eq. (7) provides a good
description when the Coulomb logarithm is taken as a free parameter. For the
large box, we show the results using two different resolutions. In both cases,
the trajectory of the BH is very small, such that we can safely assume that
the resolution for our study is large enough.

###  4.2 Dependence on the strength of the self-interactions

Figure 4: Similar to Fig. 3, we show the position and velocity of the BH
particle as a function of time. The results are shown for several interaction
strengths of the DM, indicated by the opacity of the lines. All simulations
are for the large box with the high resolution of $N_{\mathrm{DM}}=5\times
10^{6}$ particles. For the simulation of gaseous matter, we used the same
number of SPH particles.

Given that the dynamical friction force differs between the collisionless and
the gaseous case as described in Section 2, we expect that self-interactions
can have an impact on the strength of the dynamical friction force. To study
this in detail, we run simulations with SIDM. In Fig. 4, we show the results
for fSIDM with different values for the momentum transfer cross-section. We
also include the collisionless and gaseous case. It is visible that the self-
interactions can enhance or reduce the dynamical friction, depending on the
strength of the cross-section. For the gaseous regime, we expect the dynamical
friction to be enhanced above the sound speed and reduced below. Accordingly,
we find a larger velocity of the BH at the end of the simulation, when the BH
velocity has fallen well below the sound speed of the DM.

Figure 5: We show the surface density (top row) and the velocity dispersion
(bottom row) for our test simulations at a time of $t=0.49\,\mathrm{yr}$. The
panels display the different types of matter we investigate, with increasing
collisionality from the left to the right. The density and one dimensional
velocity dispersion are indicated following a logarithmic colour scale
spanning a factor of two. For all panels of the same quantity, the same colour
scale is used. In addition, we indicate the position of the BH particle with a
red dot. In the online ancillary files, we provide the time evolution of the
surface density as a video.

A detailed investigation of the dynamic friction force will follow. But
before, we take a look at the DM density perturbed by the BH. For the same
simulations as shown in Fig. 4, we show the surface density of the DM in the
top row of Fig. 5. Here, we can see that the shape and density of the wake
depend on the strength of the self-interactions. We find that for
$\sigma_{\mathrm{V}}/m\in\\{3\times 10^{-8},3\times
10^{-7}\\}\,(\mathrm{cm}^{2}\,\mathrm{g}^{-1})$ the densities in the wake are
higher and at the same time the deceleration is stronger. For these small
cross-sections, the scattering leads to more low-velocity particles in the
density wake compared to CDM effectively enhancing the deceleration of the BH.

In contrast, for the cases of even stronger self-interactions, we find that
the densities in the wake are decreasing and the deceleration of the BH is
weaker. But also note that the shape of the density wake changes for an
increasing self-interaction strength, impacting the dynamical friction force
as well. The reduced dynamical friction for the strong cross-sections is in
line with the expectation of a weaker deceleration in the gaseous case for the
subsonic regime.

In addition to the surface density, we show the velocity dispersion in the
bottom row of Fig. 5. To compute the velocity dispersion, we first subtracted
the bulk motion. In the gaseous case, we computed the velocity dispersion
directly from the temperature of the SPH particles. Similar to the density
plots, it is visible how the wake changes when self-interactions are present.
For example, the velocity dispersion in the vicinity of the perturber
increases with cross-section and the shape of the wake changes.

Figure 6: We show the dynamical friction force as a function of the Mach
number for our simulations. This plot is analogue to Fig. 2, which gives the
analytic description of the dynamical friction force. Different cross-sections
are shown, indicated by the Knudsen number. Note, $\mathrm{Kn}=\infty$
corresponds to a collisionless medium and $\mathrm{Kn}=0.0$ to a gaseous
medium.

For a more detailed analysis of the dynamical friction force and its
dependence on the cross-section, we compute the force and show it as a
function of the velocity (Fig. 6). We note, the strongest SIDM cross-section
shown is $\sigma_{\mathrm{V}}/m=3\times
10^{-6}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$. Given our resolution, we are not
able to faithfully simulate a cross-section as large as
$\sigma_{\mathrm{V}}/m=3\times 10^{-5}\mathrm{cm}^{2}\mathrm{g}^{-1}$ because
the mean free path becomes much smaller than the size of the kernel used for
the self-interactions (see Appendix A and Fischer et al., 2024a).

We specify the strength of the cross-section in terms of the Knudsen number
(Eq. 10). The collisionless case is indicated by $\mathrm{Kn}=\infty$ and the
gaseous case by $\mathrm{Kn}=0.0$. In Fig. 6 we can see that the dynamical
friction is initially increasing (starting at the highest velocity) and
decreases for later times of the simulation. At the beginning of the
simulation, where the BH moves with a supersonic velocity, the dynamical
friction is the smallest in the collisionless case. At later stages, it
depends on the cross-section whether the self-interactions enhance or reduce
the deceleration of the BH. In Fig.6, we plot the dynamical friction as a
function of the velocity. This shows that the decelerating force reaches its
maximum roughly at the sound speed, as expected. Interestingly, the self-
interacting run with an intermediate cross-section ($\mathrm{Kn}=1.35\times
10^{-1}$) shows the strongest dynamical friction, stronger than in the
collisionless or gaseous regime. In principle, the figure is built analogously
to Fig. 2, but a direct comparison is limited by the fact that our simulation
setup may not fulfill the assumptions of the Ostriker formula for the
dynamical friction in the gaseous case (see Section 2.2). This is for example
illustrated by Fig. 3, where we have seen that the exact results depend on the
choice of the setup, in particular on the box size.

Overall, we find that scatterings of the DM particles can alter the strength
of the dynamical friction in a non-trivial way. Importantly, we have shown
that the Knudsen number (Eq. 10) provides a good measure for the medium being
collisionless or collisional and its impact on the dynamical friction arising
from the self-interactions.

###  4.3 Properties of the wake

Figure 7: The density of the wake, i.e., the mean density about the point with
the highest density, is shown as a function of the velocity of the perturbing
BH. The simulation results of the different cross-sections are indicated by
their Knudsen number. Figure 8: The one-dimensional velocity dispersion of the
wake, i.e., the velocity dispersion in the vicinity of the point with the
highest density, is shown as a function of the velocity of the perturber. We
note that the velocity dispersion is computed after the bulk motion has been
subtracted. The simulation results for different cross-sections are shown,
indicated by the Knudsen number. Figure 9: The mean velocity of the wake in
the direction of the motion of the BH, i.e., the mean velocity in the vicinity
of the point with the highest density, is shown. The results for different
cross-sections are shown as a function of the velocity of the perturber. They
are indicated by their corresponding Knudsen number. Figure 10: The DM mass
gravitational bound to the perturbing BH is shown as a function of the
velocity of the perturber. We indicate the different cross-sections by their
Knudsen number.

To provide more insight into the increased dynamical friction for mildly
collisional media, we compute further quantities such as the density of the
wake, its velocity dispersion, and bulk motion as well as the DM mass
gravitational bound to the BH.

In Fig. 7, we show the mean density within a sphere of radius
$r=43.09\,\mu\mathrm{pc}$ about the point of the maximum DM density in the
wake as a function of the velocity of the perturber. The simulations showing a
strong dynamical friction force also show a larger density compared to the
ones with relatively weak dynamical friction. This is not surprising, as a
more massive wake can exert a stronger gravitational pull on the BH and thus
faster decelerate it. The presence of self-interactions allows for the growth
of a larger wake density, as the scattering produces more low-velocity
particles compared to the collisionless case. However, for sufficiently large
self-interaction when entering the short-mean-free-path regime, the wake
density starts to decrease with increasing cross-section.

Furthermore, we show the velocity dispersion in the density wake in Fig. 8. We
use the same volume as for the density computation. The velocity dispersion is
computed after subtracting the mean velocity of the particles. The velocity
dispersion in the wake increases compared to the initial conditions and is
larger for stronger DM self-interactions. Interestingly, the cross-section
with $\mathrm{Kn}=1.35\times 10^{-1}$, which gives the strongest dynamical
friction, shows a substantially lower velocity dispersion compared to the
larger cross-sections.

The bulk motion of the density wake is displayed in Fig. 9. The stronger the
self-interactions the more efficiently the BH drags the DM along its
trajectory, except for the largest cross-section for which the bulk velocity
is again smaller. This is in line with the strength of the dynamical friction
force found in Fig. 6, i.e., when the DM is dragged along by the BH, it can
more efficiently decelerate the BH.

Last, we compute the DM mass that is gravitationally bound to the BH. To do
so, we compare the relative velocity of a DM particle and the BH to the escape
velocity,

$v_{\mathrm{esc}}=\sqrt{\frac{2\,\mathrm{G}\,m_{\mathrm{BH}}}{r}}\,,$ (14)

of the BH. Here, the perturbing BH has a mass of $m_{\mathrm{BH}}$ and the
gravitational constant is given by G. The distance between the DM particle and
the BH is denoted by $r$. This may underestimate the bound mass to some extent
as the DM bound to the BH increases the effective mass and thus may lead to
further particles being gravitationally bound.

Our estimate for the gravitational bound mass is shown in Fig. 10. It becomes
visible that the bound mass depends on the cross-section. In particular,
strong self-interactions can effectively increase the bound mass. This can be
understood by the effect that the scattering has on the velocity distribution.
We note that the velocities in the wake do not follow a Maxwell–Boltzmann
distribution, but the scattering evolves the distribution towards it. This
implies that the self-interactions create more low-velocity particles compared
to the collisionless case. In consequence, more DM particles are
gravitationally bound to the BH. This effectively enhances the dynamical
friction acting on the perturber given the cross-section is small enough such
that the behaviour known from the gaseous regime (Ostriker, 1999) does not
become dominant. Although the dynamical friction decreases with increasing
cross-section for strong self-interactions, the gaseous regime nevertheless
yields a stronger dynamical friction force in the supersonic regime compared
to the collisionless case.

###  4.4 Angular dependence

The angular dependence of the self-interaction cross-section can in principle
lead to qualitative differences in the evolution of astrophysical systems.
This is for example well known for merging DM halos (e.g., Fischer et al.,
2023; Sabarish et al., 2024) or the abundance of satellites (Fischer et al.,
2022). However, for systems that are rather relaxed, the angular dependence
plays a minor role only and the effect of self-interactions with different
angular dependencies can be mapped onto each other using the viscosity cross-
section $\sigma_{\mathrm{V}}$ (Eq. 11, see also Yang & Yu, 2022).

To test whether this mapping also works for our setup to study dynamical
friction, we run simulations with a very anisotropic cross-section (using the
fSIDM framework) and an isotropic cross-section. The cross-sections are chosen
to have the same value in terms of $\sigma_{\mathrm{V}}/m$ or
$\sigma_{\mathrm{T}}/m$ to test the matching with the momentum transfer cross-
section (see Eq. 12) as well.

In Fig. 11, we compare the velocity of the BH for the different angular
dependencies. We compute the ratio of the change in BH velocities for
simulations with a matched isotropic cross-section to the fiducial fSIDM
simulations. The two panels in Fig. 11 display the ratio of the velocities as
a function of time for a weaker fSIDM cross-section of
$\sigma_{\mathrm{V}}/m=3\times 10^{-8}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ (top
panel) and a stronger one of $\sigma_{\mathrm{V}}/m=3\times
10^{-7}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ (bottom panel). For a perfect
matching, we would expect the ratio to be unity. As we can see, the viscosity
cross-section (Eq. 11) provides a better match than the transfer cross-section
(Eq. 12). This is the case for the weaker and stronger cross-section that we
have simulated to test this.

Figure 11: We compute the ratio between the velocity change, $\Delta
v=v(0)-v(t)$, due to dynamical friction of an fSIDM cross-section and the
transfer (blue) / viscosity (orange) matched isotropic cross-sections. The
upper panel is for an fSIDM viscosity cross-section of
$\sigma_{\mathrm{V}}/m=3\times 10^{-8}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ and
the lower panel is for $\sigma_{\mathrm{V}}/m=3\times
10^{-7}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$.

###  4.5 Velocity dependence

The effective cross-section $\sigma_{\mathrm{eff}}$ (introduced by Yang & Yu,
2022) allows to map cross-sections with different velocity dependencies onto
each other. It integrates the cross-section with its velocity dependence over
the assumed Maxwell-Boltzmann distribution of the scattering velocities,
applying a weighting factor of $v^{5}$. The effective cross-section is given
by

$\sigma_{\mathrm{eff}}=\frac{\langle
v^{5}\sigma_{\mathrm{V}}(v)\rangle}{\langle v^{5}\rangle}\,.$ (15)

Note that this is very similar to the matching procedure proposed by Yang et
al. (2023); see also the work by Outmezguine et al. (2023). It has been shown
that this works fairly well for isolated halos and, in turn, velocity-
dependent cross-sections can be approximately described with constant cross-
sections. However, it is known that this does not work for systems that
involve multiple velocity scales, such as merging galaxy clusters (Sabarish et
al., 2024) or halos containing substructure (Fischer et al., 2024b).

Here, we investigate how well the matching works for our setup to study
dynamical friction. In Fig. 12, we show simulations involving different
velocity dependencies specified in terms of $w$, but they are all matched
using Eq. (15). We have simulated very anisotropic cross-sections (top panel)
and isotropic ones (bottom panel). For times up to $t\approx
0.5\,\mathrm{yr}$, the difference between the runs is fairly small and
increases at later stages when the BH has slowed down. We find that the
matching of the velocity dependence using the effective cross-section does not
work as well as for the angular dependence using the viscosity cross-section.
However, in our case where the velocity of the BH is falling below half its
initial value over the course of the simulation, the deviation between the
runs hardly exceeds 2.5%. This might be enough for a rough estimate, and
separate simulations for velocity-dependent cases may only be needed when a
higher accuracy is required.

Figure 12: The ratio of the velocity change, $\Delta v=v(0)-v(t)$, due to
dynamical friction acting on the BH for a velocity-dependent and a velocity-
independent (constant) cross-section is shown. All cross-sections were matched
using the effective cross-section (Eq. 15), i.e., all cross-sections
correspond to $\sigma_{\mathrm{eff}}/m=3\times
10^{-8}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$. The upper panel displays runs with
frequent self-interactions, and the lower panel is for isotropic scattering.

##  5 Discussion

In this section, we discuss the limitations of our study, mainly the
simulation setup. Moreover, we elaborate on the implications of our results
and further directions to provide a better understanding of the dynamical
friction in merging binary BHs.

Currently, there does not exist a theory for the non-equilibrium dynamics of
SIDM to describe dynamical friction. This is why we have to use simulations to
understand the impact of a mildly collisional medium on the dynamical friction
force. Ideally, one would develop a more general formalism that contains both
the Chandrasekhar and the Ostriker formula in the limits of a collisionless
and a gaseous medium but also describes the intermediate regime as a function
of the self-interaction cross-section. But even then one may still depend on
simulations to calibrate such a formalism and make it applicable to problems
such as merging binary BHs. Moreover, we want to note that the Ostriker
formula assumes the trajectory of the perturber to be a straight line. Kim &
Kim (2007) improved on this and provided a generalisation to circular orbits.

For binary BHs, we find that the motion of the secondary BH through the DM
spike would typically happen at Mach numbers larger than unity, i.e., in the
supersonic regime. Moreover, the Mach number does only very slightly depend on
the radius, as the mass of the DM spike in the relevant radial range is
typically much smaller than the mass of the primary BH. Also, the spike index
does not play a significant role in the motion of the secondary BH being
supersonic. But more importantly, a non-zero eccentricity could lead to a
subsonic phase about the apocentre of the orbit of the secondary BH. This
implies that the supersonic regime where the dynamical friction is enhanced by
self-interactions plays an important role and would lead to a faster decay of
the orbit compared to a collisionless medium with the same density.

To study the dynamical friction with self-interactions being present,
simulations are required. We used an idealized setup with a periodic box
hosting initially a constant density. We have to note that this setup comes
with drawbacks that limit the conclusion we can draw. In particular, it is
difficult to create a stable setup where the box size is small enough such
that the velocity dispersion of the DM counteracts the gravitational collapse.
This effectively limits the maximal box size that can be studied, which also
has an impact on the dynamical friction. This is because the box size is
related to the maximum impact parameter $r_{\mathrm{max}}$ (see Eq. 7).
Moreover, the setup does not always fulfill the assumptions of analytic
formulas discussed in section 2.2. Also, we cannot exclude that perturbations
arising from the supersonic phase still affect the trajectory of the BH in the
subsonic phase. Nevertheless, we were still able to make qualitative
statements and gain new insights into the dynamical friction. Especially for
this purpose, we tuned the setup to give a relatively larger impact on the BH
velocity and make it easier to obtain qualitative differences being much
larger than the numerical errors.

To obtain accurate quantitative predictions on the strength of the dynamical
friction in BH mergers, one would need to simulate a setup much closer to the
real system. Kavanagh et al. (2020, 2024) studied intermediate mass-ratio
inspirals resolving the DM spike around the primary BH and measured the
deceleration due to the dynamical friction of the secondary. In principle,
such a setup could be used for studies of SIDM as well. Potentially this does
not only allow us to calibrate the Chandrasekhar formula, but also the
Ostrikers formula for the gaseous regime and with this make it applicable for
merging BHs. Ideally, one would have a more general formalism for arbitrary
cross-sections that can be calibrated as well.

In our simulations, we find a regime of intermediate cross-sections where the
dynamical friction is stronger than in the collisionless case or the gaseous
case. Given the drawback our setup suffers from, it would be great to study
this in a setup closer to the physical systems as described above. Simulations
of BHs moving through a DM spike could potentially tell whether the
intermediate regime of enhanced dynamical friction is only due to our setup or
would also occur in real BH mergers.

The prediction of the exact strength of the dynamical friction force and its
impact on the GW signal is difficult. We find that not only does the DM
density matter but its collisionless or collisional nature, too. Moreover, the
latter plays a crucial role in the spike index and the DM density (Shapiro &
Paschalidis, 2014). But it must also be noted that the energy that is
transferred due to dynamical friction from the secondary BH into the DM spike
affects the DM spike as well. Alonso-Álvarez et al. (2024) argues that only
SIDM can effectively transport this heat outwards and replenish the DM spike.
In contrast to collisionless matter, this effectively allows a SIDM spike to
make the BH orbit decay and solve the final parsec problem. Overall, it
remains challenging to predict the exact GW signal, especially beyond the
Newtonian regime.

##  6 Conclusions

In this paper, we have reviewed analytic estimates for the dynamical friction
force in the collisionless and gaseous regime, and discussed the two regimes
in the context of binary BH mergers. Furthermore, we have performed idealized
simulations to study the dynamical friction as a function of the self-
interaction cross-section. In addition, we also studied how well cross-
sections with different angular and velocity dependencies can be matched onto
each other. Our main results can be summarized as follows:

  1. 1.

DM self-interactions can enhance but also lower the dynamical friction force,
depending on the velocity of the perturber relative to the local velocity
dispersion of the DM.

  2. 2.

For merging binary BHs, we expect the perturber to be in the supersonic regime
all the time when being in a circular orbit and partially when following an
eccentric orbit. In this regime, self-interactions enhance the dynamical
friction force.

  3. 3.

The Knudsen number, Kn, is a good proxy for how much the self-interactions
impact the dynamical friction, i.e., how close the system is to the
collisionless or gaseous regime.

  4. 4.

When the SIDM cross-section is large enough, $\mathrm{Kn}\lesssim 2$, we
expect the dynamical friction in the supersonic regime to be significantly
larger than in the collisionless case. Note, for typical binary BHs such a
cross-section is many orders of magnitude smaller than current upper bounds on
the self-interaction cross-section (Adhikari et al., 2022).

  5. 5.

This implies that observed BH mergers, when being consistent with
collisionless matter, can potentially imply extremely tight upper bounds on
the SIDM cross-section.

  6. 6.

We find that the viscosity cross-section provides a good matching procedure
between isotropic and very anisotropic cross-sections when computing the
effect of self-interactions on the dynamical friction.

  7. 7.

The matching of different velocity dependencies onto each other works less
well than the angular matching. However, depending on the application, the
deviation could still be within an acceptable range.

The results we obtained are solely based on simulations of a constant density
perturbed by a BH in a periodic box. This allowed us to gain insight into the
qualitative behaviour but does not allow for accurately quantifying how strong
the effect of self-interactions is on the dynamical friction force. In
consequence, the next step would be to simulate a SIDM spike with an orbiting
BH, similar to the work done by Kavanagh et al. (2020, 2024). Such a setup
would allow us to precisely determine the strength of the dynamical friction
force acting onto BHs orbiting in DM spikes when self-interactions are
present. This is the subject of a forthcoming paper (Sabarish et al., in
prep.).

###### Acknowledgements.

The authors thank the organizers of the Pollica 2023 SIDM Workshop, where part
of this work was done. In addition, the authors thank all participants of the
Darkium SIDM Journal Club for discussions. MSF and LS thank Niklas Becker,
Andrew Benson, Akaxia Cruz, Xiaolong Du, Frederick Groth, Charlie Mace, Annika
H. G. Peter, Sabarish Venkataramani, Daneng Yang, Shengqi Yang, Hai-Bo Yu,
Xingyu Zhang, and Zhichao Carton Zeng as well as the members of the CAST group
at the University Observatory Munich for helpful discussions. This work is
funded by the _Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation)_ under Germany’s Excellence Strategy – EXC-2094 “Origins” –
390783311. Software: NumPy (Harris et al., 2020), Matplotlib (Hunter, 2007).

## References

  * Abbott et al. (2016) Abbott, B., Abbott, R., Abbott, T., et al. 2016, Phys. Rev. Lett., 116
  * Abbott et al. (2019) Abbott, B., Abbott, R., Abbott, T., et al. 2019, Phys. Rev. X, 9
  * Abbott et al. (2017a) Abbott, B., Abbott, R., Abbott, T., et al. 2017a, Phys. Rev. Lett., 119
  * Abbott et al. (2017b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJ, 848, L12
  * Adhikari et al. (2022) Adhikari, S., Banerjee, A., Boddy, K. K., et al. 2022, Astrophysical Tests of Dark Matter Self-Interactions
  * Agazie et al. (2023) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023, ApJ, 951, L8
  * Alachkar et al. (2023) Alachkar, A., Ellis, J., & Fairbairn, M. 2023, Phys. Rev. D, 107
  * Alonso-Álvarez et al. (2024) Alonso-Álvarez, G., Cline, J. M., & Dewar, C. 2024, Self-interacting dark matter solves the final parsec problem of supermassive black hole mergers
  * Alvarez & Yu (2021) Alvarez, G. & Yu, H.-B. 2021, Phys. Rev. D, 104
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, Laser Interferometer Space Antenna
  * Andrade et al. (2021) Andrade, K. E., Fuson, J., Gad-Nasr, S., et al. 2021, MNRAS, 510, 54
  * Antoniadis et al. (2023) Antoniadis, J., Arumugam, P., Arumugam, S., et al. 2023, A&A, 678, A50
  * Baiotti (2022) Baiotti, L. 2022, Arab. J. Math., 11, 105
  * Beck et al. (2016) Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110
  * Becker & Sagunski (2023) Becker, N. & Sagunski, L. 2023, Phys. Rev. D, 107
  * Becker et al. (2022) Becker, N., Sagunski, L., Prinz, L., & Rastgoo, S. 2022, Phys. Rev. D, 105
  * Berezhiani et al. (2023) Berezhiani, L., Cintia, G., Luca, V. D., & Khoury, J. 2023, Dynamical friction in dark matter superfluids: The evolution of black hole binaries
  * Bertone (2024) Bertone, G. 2024, Dark matter, black holes, and gravitational waves
  * Bertone et al. (2024) Bertone, G., Wierda, A. R. A. C., Gaggero, D., et al. 2024, Dark Matter Mounds: towards a realistic description of dark matter overdensities around black holes
  * Bertone et al. (2005) Bertone, G., Zentner, A. R., & Silk, J. 2005, Phys. Rev. D, 72
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic dynamics, 2nd edn., Princeton Series in Astrophysics (Princeton University Press)
  * Boey et al. (2024) Boey, R., Wang, Y., Kendall, E., & Easther, R. 2024, Dynamical Friction and Black Holes in Ultralight Dark Matter Solitons
  * Boudon et al. (2022) Boudon, A., Brax, P., & Valageas, P. 2022, Phys. Rev. D, 106
  * Boudon et al. (2023a) Boudon, A., Brax, P., & Valageas, P. 2023a, Supersonic friction of a black hole traversing a self-interacting scalar dark matter cloud
  * Boudon et al. (2023b) Boudon, A., Brax, P., Valageas, P., & Wong, L. K. 2023b, Gravitational waves from binary black holes in a self-interacting scalar dark matter cloud
  * Capozziello et al. (2023) Capozziello, S., Zare, S., & Hassanabadi, H. 2023, Testing bumblebee gravity with global monopoles in a dark matter spike by EHT observations from M87 and Sgr A
  * Chan & Lee (2024) Chan, M. H. & Lee, C. M. 2024, The first robust evidence showing a dark matter density spike around the supermassive black hole in OJ 287
  * Chandrasekhar (1943) Chandrasekhar, S. 1943, ApJ, 97, 255
  * Colquhoun et al. (2021) Colquhoun, B., Heeba, S., Kahlhoefer, F., Sagunski, L., & Tulin, S. 2021, Phys. Rev. D, 103, 035006
  * Coogan et al. (2022) Coogan, A., Bertone, G., Gaggero, D., Kavanagh, B. J., & Nichols, D. A. 2022, Phys. Rev. D, 105
  * Daghigh & Kunstatter (2024) Daghigh, R. G. & Kunstatter, G. 2024, Who knows what dark matter lurks in the heart of M87: The shadow knows, and so does the ringdown
  * Dehnen & Aly (2012) Dehnen, W. & Aly, H. 2012, MNRAS, 425, 1068
  * Despali et al. (2022) Despali, G., Walls, L. G., Vegetti, S., et al. 2022, MNRAS, 516, 4543
  * Dolag et al. (2005) Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753
  * Dosopoulou (2023) Dosopoulou, F. 2023, Dynamical friction in dark matter spikes: corrections to Chandrasekhar’s formula
  * Eckert, D. et al. (2022) Eckert, D., Ettori, S., Robertson, A., et al. 2022, A&A, 666, A41
  * Eda et al. (2015) Eda, K., Itoh, Y., Kuroyanagi, S., & Silk, J. 2015, Phys. Rev. D, 91
  * Einstein (1916) Einstein, A. 1916, Annalen der Physik, 354, 769
  * Fields et al. (2014) Fields, B. D., Shapiro, S. L., & Shelton, J. 2014, Phys. Rev. Lett., 113
  * Fischer et al. (2021a) Fischer, M. S., Brüggen, M., Schmidt-Hoberg, K., et al. 2021a, MNRAS, 505, 851
  * Fischer et al. (2022) Fischer, M. S., Brüggen, M., Schmidt-Hoberg, K., et al. 2022, MNRAS, 516, 1923
  * Fischer et al. (2021b) Fischer, M. S., Brüggen, M., Schmidt-Hoberg, K., et al. 2021b, MNRAS, 510, 4080
  * Fischer et al. (2024a) Fischer, M. S., Dolag, K., & Yu, H.-B. 2024a, Numerical challenges for energy conservation in N-body simulations of collapsing self-interacting dark matter haloes
  * Fischer et al. (2023) Fischer, M. S., Durke, N.-H., Hollingshausen, K., et al. 2023, MNRAS, 523, 5915–5933
  * Fischer et al. (2024b) Fischer, M. S., Kasselmann, L., Brüggen, M., et al. 2024b, MNRAS, 529, 2327–2348
  * Fitts et al. (2019) Fitts, A., Boylan-Kolchin, M., Bozek, B., et al. 2019, MNRAS, 490, 962–977
  * Gingold & Monaghan (1977) Gingold, R. A. & Monaghan, J. J. 1977, MNRAS, 181, 375
  * Glennon et al. (2024) Glennon, N., Musoke, N., Nadler, E. O., Prescod-Weinstein, C., & Wechsler, R. H. 2024, Phys. Rev. D, 109
  * Gondolo & Silk (1999) Gondolo, P. & Silk, J. 1999, Phys. Rev. Lett., 83, 1719–1722
  * Gopika & Desai (2023) Gopika, K. & Desai, S. 2023, Physics of the Dark Universe, 42, 101291
  * Groth et al. (2023) Groth, F., Steinwandel, U. P., Valentini, M., & Dolag, K. 2023, MNRAS, 526, 616
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Harvey et al. (2019) Harvey, D., Robertson, A., Massey, R., & McCarthy, I. G. 2019, MNRAS, 488, 1572
  * Hu & Wu (2017) Hu, W.-R. & Wu, Y.-L. 2017, National Science Review, 4, 685
  * Hulse & Taylor (1975) Hulse, R. A. & Taylor, J. H. 1975, ApJ, 195, L51
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
  * John et al. (2023) John, I., Leane, R. K., & Linden, T. 2023, Dark Matter Scattering Constraints from Observations of Stars Surrounding Sgr A*
  * Just & Kegel (1990) Just, A. & Kegel, W. H. 1990, A&A, 232, 447
  * Just et al. (1986) Just, A., Kegel, W. H., & Deiss, B. M. 1986, A&A, 164, 337
  * Just et al. (2010) Just, A., Khan, F. M., Berczik, P., Ernst, A., & Spurzem, R. 2010, MNRAS, 411, 653–674
  * Just & Peñarrubia (2005) Just, A. & Peñarrubia, J. 2005, A&A, 431, 861–877
  * Kahlhoefer et al. (2014) Kahlhoefer, F., Schmidt-Hoberg, K., Frandsen, M. T., & Sarkar, S. 2014, MNRAS, 437, 2865
  * Karydas et al. (2024) Karydas, T. K., Kavanagh, B. J., & Bertone, G. 2024, Sharpening the dark matter signature in gravitational waveforms I: Accretion and eccentricity evolution
  * Kavanagh et al. (2024) Kavanagh, B. J., Karydas, T. K., Bertone, G., Cintio, P. D., & Pasquato, M. 2024, Sharpening the dark matter signature in gravitational waveforms II: Numerical simulations with the NbodyIMRI code
  * Kavanagh et al. (2020) Kavanagh, B. J., Nichols, D. A., Bertone, G., & Gaggero, D. 2020, Phys. Rev. D, 102
  * Kim & Kim (2007) Kim, H. & Kim, W. 2007, ApJ, 665, 432–444
  * Koda & Shapiro (2011) Koda, J. & Shapiro, P. R. 2011, MNRAS, 415, 1125
  * Kong et al. (2024) Kong, D., Yang, D., & Yu, H.-B. 2024, CDM and SIDM Interpretations of the Strong Gravitational Lensing Object JWST-ER1
  * Kremer (2022) Kremer, G. M. 2022, Universe, 8, 179
  * Lancaster et al. (2020) Lancaster, L., Giovanetti, C., Mocz, P., et al. 2020, J. Cosmology Astropart. Phys., 2020, 001–001
  * Macedo et al. (2013) Macedo, C. F. B., Pani, P., Cardoso, V., & Crispino, L. C. B. 2013, ApJ, 774, 48
  * Milosavljević & Merritt (2003) Milosavljević, M. & Merritt, D. 2003, ApJ, 596, 860–878
  * Monaghan & Lattanzio (1985) Monaghan, J. J. & Lattanzio, J. C. 1985, A&A, 149, 135
  * Montalvo et al. (2024) Montalvo, D., Smith-Orlik, A., Rastgoo, S., et al. 2024, Post-Newtonian effects in compact binaries with a dark matter spike: A Lagrangian approach
  * Mukherjee et al. (2024) Mukherjee, D., Holgado, A. M., Ogiya, G., & Trac, H. 2024, Examining the Effects of Dark Matter Spikes on Eccentric Intermediate Mass Ratio Inspirals Using $N$-body Simulations
  * Ostriker (1999) Ostriker, E. C. 1999, ApJ, 513, 252–258
  * Outmezguine et al. (2023) Outmezguine, N. J., Boddy, K. K., Gad-Nasr, S., Kaplinghat, M., & Sagunski, L. 2023, MNRAS, 523, 4786
  * Pollack et al. (2015) Pollack, J., Spergel, D. N., & Steinhardt, P. J. 2015, ApJ, 804, 131
  * Ragagnin et al. (2024) Ragagnin, A., Meneghetti, M., Calura, F., et al. 2024, Dianoga SIDM: galaxy cluster self-interacting dark matter simulations
  * Reardon et al. (2023) Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6
  * Rephaeli & Salpeter (1980) Rephaeli, Y. & Salpeter, E. E. 1980, ApJ, 240, 20
  * Robertson et al. (2019) Robertson, A., Harvey, D., Massey, R., et al. 2019, MNRAS, 488, 3646–3662
  * Sabarish et al. (in prep.) Sabarish, V. M., Brüggen, M., Schmidt-Hoberg, K., Fischer, M. S., & et al. in prep., arXiv e-prints
  * Sabarish et al. (2024) Sabarish, V. M., Brüggen, M., Schmidt-Hoberg, K., Fischer, M. S., & Kahlhoefer, F. 2024, MNRAS, 529, 2032
  * Sadeghian et al. (2013) Sadeghian, L., Ferrer, F., & Will, C. M. 2013, Phys. Rev. D, 88
  * Sagunski et al. (2021) Sagunski, L., Gad-Nasr, S., Colquhoun, B., Robertson, A., & Tulin, S. 2021, J. Cosmology Astropart. Phys., 2021, 024
  * Shapiro (2018) Shapiro, S. L. 2018, Phys. Rev. D, 98
  * Shapiro & Paschalidis (2014) Shapiro, S. L. & Paschalidis, V. 2014, Phys. Rev. D, 89
  * Spergel & Steinhardt (2000) Spergel, D. N. & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  * Szölgyén et al. (2022) Szölgyén, Á., MacLeod, M., & Loeb, A. 2022, MNRAS, 513, 5465–5473
  * Tang et al. (2021) Tang, M., Xu, Z., & Wang, J. 2021, Chinese Physics C, 45, 015110
  * Taylor et al. (1976) Taylor, J. H., Hulse, R. A., Fowler, L. A., Gullahorn, G. E., & Rankin, J. M. 1976, ApJ, 206, L53
  * Traykova et al. (2023) Traykova, D., Vicente, R., Clough, K., et al. 2023, Phys. Rev. D, 108
  * Tulin & Yu (2018) Tulin, S. & Yu, H.-B. 2018, Phys. Rep, 730, 1
  * Ullio et al. (2001) Ullio, P., Zhao, H., & Kamionkowski, M. 2001, Phys. Rev. D, 64
  * Vogelsberger et al. (2014) Vogelsberger, M., Zavala, J., Simpson, C., & Jenkins, A. 2014, MNRAS, 444, 3684
  * Wang & Easther (2022) Wang, Y. & Easther, R. 2022, Phys. Rev. D, 105
  * Weisberg et al. (1981) Weisberg, J. M., Taylor, J. H., & Fowler, L. A. 1981, Scientific American, 245, 74
  * Yang & Yu (2022) Yang, D. & Yu, H.-B. 2022, J. Cosmology Astropart. Phys., 2022, 077
  * Yang et al. (2023) Yang, S., Du, X., Zeng, Z. C., et al. 2023, ApJ, 946, 47
  * Yang et al. (2024) Yang, S., Jiang, F., Benson, A., et al. 2024, A quantitative comparison between velocity dependent SIDM cross sections constrained by the gravothermal and isothermal models
  * Yue & Han (2018) Yue, X.-J. & Han, W.-B. 2018, Phys. Rev. D, 97
  * Zeng et al. (2022) Zeng, Z. C., Peter, A. H. G., Du, X., et al. 2022, MNRAS, 513, 4845
  * Zhang et al. (2024) Zhang, X., Yu, H.-B., Yang, D., & An, H. 2024, Self-interacting dark matter interpretation of Crater II

##  Appendix A Unresolved mean free path

In this appendix, we show the results of two simulations with a larger cross-
section than the ones we have shown in the main text. These additional
simulations are for fSIDM and employ a cross-section of
$\sigma_{\mathrm{V}}/m=3\times 10^{-5}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$.
They differ in resolution, we use $N=1\times 10^{6}$ and $N=5\times 10^{6}$
particles. In Fig. 13, we can see, in contrast to the previous simulations,
that the results for the two resolutions differ significantly from each other.

With increasing cross-section, the mean free path is shrinking and may fall
below the kernel size used for SIDM (see Fischer et al. 2024a). The ratio
between the kernel size, $h$, and the mean free path $l$, is fairly large for
these simulations, $h/l\approx 172$ for the lower resolution run and
$h/l\approx 100$ for the higher resolution. Thus, it is not surprising that
the simulations deviate from each other. The value for $h/l$ at which the
simulations start to become inaccurate may in general depend on the problem at
hand. The values of our simulations here may not be indicative of other
setups.

Figure 13: The position and velocity of the BH in our test setup is shown for
a simulation with a fSIDM cross-section of $\sigma_{\mathrm{V}}/m=3\times
10^{-5}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$. The high resolution corresponds to
$N=5\times 10^{6}$ particles and the low resolution to $N=1\times 10^{6}$
particles. The larger box size of $l_{\mathrm{box}}=829.5\,\mu\mathrm{pc}$ was
used.


admits a unique smooth solution $w_0^m$. Then, it follows from the standard $H^2$-estimate to the above elliptic system and Sobolev's inequality that
$w_0^m\rightarrow 0$ in $H^2$ as $m\rightarrow \infty$.
Denote $u_0^m=v_0^m-w_0^m$, it holds that $u_0^m\in C^\infty$ satisfy
5d3 u_0^m→u_0  in  H^2  m→∞,     u_0^m·n=u_0^m×n=0  on  Ø.Hence, for any fixed  $\eta \in (0, \eta_0),$ there exists $m_1(\eta)\geq m_0(\eta)$
such that for any $m\geq m_1(\eta)$
5d4u_0^m_L^2≤(1+η)u_0_L^2.We set $u_0^{m,\eta}= \frac{1}{1+\eta}u_0^m$.

Similarly, we will construct  $\te_0^{m,\eta}$.
It's easy to find a sequence  $0\le f_0^m\in C^\infty$ such that $f_0^m\rightarrow \te_0$ in $ H^2$ as $m\rightarrow \infty$.
For any fixed $f_0^m$, considering the unique  smooth solution $h_0^m$ to the Poisson's equation

Δh_0^m=Δf_0^m,   in  Ø,

h_0^m·n=f_0^m·n,   on  Ø,
with $\overline{h_0^m}=0,$  it is easy to deduce that   $h_0^m\rightarrow 0$ in $H^2$ as $m\rightarrow \infty$.
Denote $\te_0^{m}= f_0^m-h_0^m$, it holds that $\te_0^m\in C^\infty$  satisfy
5d5 _0^m→_0  in  H^2  m→∞,     _0^m·n=0  on  Ø.Hence, for any  $\eta \in (0, \eta_0),$ there exists $m_2(\eta)\geq m_1(\eta)$
such that for any $m\geq m_2(\eta),$
5d6_0^m-_0_L^∞≤η/4.We finally set $\te_0^{m,\eta}= \frac{1}{1+\eta}(\te_0^m+\frac{\eta}{2})$.

Therefore, one deduces from (<ref>) and (<ref>)-(<ref>) that
$(\n_0^{m,\eta},u_0^{m,\eta}, \te_0^{m,\eta}, H_0^{m,\eta})$ satisfies
de3(_0^m,η,u_0^m,η, _0^m,η, H_0^m,η)∈C^∞,

η≤_0^m,η  <,    η/8≤_0^m,η ≤,   u_0^m,η_L^2 ≤M,

u_0^m,η·n  =u_0^m,η×n=_0^m,η·n=0   on Ø,
provided that $\eta \in (0, \eta_0)$, $m\geq m_2(\eta)$, and de03
lim_m→∞ lim_η→0
( _0^m,η - _0 _ W^2,q+ u_0^m,η-u_0_H^2+ _0^m,η- _0  _H^2)=0.
Moreover,  the initial norm $C_0^{m,\eta}$
for $(\n_0^{m,\eta},u_0^{m,\eta}, \te_0^{m,\eta}, H_0^{m,\eta}),$   the right hand side of (<ref>)
with $(\n_0,u_0,\te_0,H_0)$   replaced by
$(\n_0^{m,\eta},u_0^{m,\eta}, \te_0^{m,\eta}, H_0^{m,\eta}),$
satisfy lim_m→∞ lim_η→0C_0^m,η=C_0.Therefore, there exists  an  $\eta_1\in(0, \eta_0) $
such that, for any $\eta\in(0,\eta_1),$ we can find some $m_3(\eta)\geq m_2(\eta)$  such that   de1 C_0^m,η≤C_0+_0/3≤_0 , provided thatde70<η<η_1 ,   m≥m_3(η).

We assume that $m,\eta$ satisfy (<ref>).
Proposition <ref> together with (<ref>) and (<ref>) thus yields that
there exists a smooth solution  $(\n^{m,\eta},u^{m,\eta}, \te^{m,\eta}, H^{m,\eta}) $
of problem (<ref>)–(<ref>) with  initial data $(\n_0^{m,\eta},u_0^{m,\eta}, \te_0^{m,\eta}, H_0^{m,\eta})$
on $\Omega\times (0,T] $ for all $T>0. $
Moreover, one has (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>)   with $(\n,u,\te,H)$  being replaced by $(\n^{m,\eta},u^{m,\eta}, \te^{m,\eta}, H^{m,\eta}).$

Next, for the initial data $(\n_0^{m,\eta},u_0^{m,\eta}, \te_0^{m,\eta}, H_0^{m,\eta})$, the function $\tilde g$ in (<ref>)  is
co5g̃     ≜(_0^m,η)^-1/2(-μΔu_0^m,η-(μ+λ)÷u_0^m,η+R(_0^m,η^m,η_0)-H_0^m,η× H_0^m,η)

= (_0^m,η)^-1/2√(_0)g+μ(_0^m,η)^-1/2Δ(u_0-u_0^m,η)

+(μ+λ) (_0^m,η)^-1/2 ÷(u_0-u_0^m,η)+ R(_0^m,η)^-1/2 (_0^m,η_0^m,η-_0_0),

+(_0^m,η)^-1/2(H_0× H_0-H_0^m,η× H_0^m,η),where in the second equality we have used (<ref>).

Since $g \in L^2,$ one deduces from (<ref>),  (<ref>), and  (<ref>)  that for any $\eta\in(0,\eta_1),$ there exist some $m_3(\eta)\geq m_2(\eta)$ and a positive constant $C$ independent of $m$ and $\eta$ such that
g̃_L^2≤g_L^2+Cη^-1/2(m) + Cη^1/2,
with   $0\le\de(m) \rightarrow 0$ as
$m \rightarrow \infty.$ Hence,  for any  $\eta\in(0,\eta_1),$ there exists some $m_4(\eta)\geq m_3(\eta)$ such that for any $ m\geq m_4(\eta)$,
de9(m) <η.We thus obtain from (<ref>) and (<ref>) that
there exists some positive constant $C$ independent of $m$ and $\eta$ such that  de14  g̃ _L^2≤g _L^2+C,provided thatde10 0<η<η_1,    m≥m_4(η).

Now, we   assume that $m,$  $\eta$ satisfy (<ref>).
It thus follows from (<ref>)–(<ref>),  (<ref>), Proposition <ref>,
and Lemmas <ref>, <ref>–<ref> that for any $T>0,$
there exists some positive constant $C$ independent of $m$ and $\eta$ such that
(<ref>), (<ref>),   (<ref>),  (<ref>), (<ref>), (<ref>), (<ref>),  (<ref>),  (<ref>), (<ref>),
and  (<ref>)  hold for  $(\n^{m,\eta},u^{m,\eta}, \te^{m,\eta}, H^{m,\eta}) .$
Then passing  to the limit first $m\rightarrow \infty,$ then $\eta\rightarrow 0,$
together with standard arguments yields that there exists a solution $(\n,u,\te,H)$ of the problem (<ref>)–(<ref>)
on $\Omega\times (0,T]$ for all $T>0$, such that the solution  $(\n,u,\te,H)$
satisfies  (<ref>), (<ref>),    (<ref>),  (<ref>),  (<ref>), (<ref>), (<ref>),  (<ref>), (<ref>),  (<ref>),
and  the estimates of $A_i(T)\,(i=1,\cdots,4)$ in
(<ref>). Hence, $(\n,u,\te,H)$ satisfying  (<ref>) and (<ref>) can be proved in the same way as that in [19]. Moreover, one deduces from Proposition <ref> that the desired exponential decay property (<ref>).
Furthermore,  the proof of the uniqueness of $(\n,u,\te,H)$ is similar to that  of [7] and we omit the details.
The proof of Theorem <ref> is finished.

Proof of Theorem <ref>.
As is shown by [3], we sketch the proof for completeness. First, we show that, for $T>0$, the Lagrangian coordinates of the system are given by
c61  ∂/∂τX(τ; t,x) =u(X(τ; t,x),τ),     0≤τ≤T

X(t;t,x)=x,      0≤t≤T, x∈Ω̅.By (<ref>), the transformation (<ref>) is well-defined. In addition, by $\eqref{CMHD}_1$, we have
c62ρ(x,t)=ρ_0(X(0; t, x)) exp{-∫_0^t÷u(X(τ;t, x),τ)dτ}.
If there exists some point $x_0\in \Omega$ such that $\n_0(x_0)=0,$ then there is a point $x_0(t)\in \bar{\Omega}$ such that $X(0; t, x_0(t))=x_0$. Hence, by (<ref>), $\rho(x_0(t),t)\equiv 0$ for any $t\geq 0.$
By (<ref>), we get for $ \tilde{r}\in  (3,\infty)$ and $\theta=\frac{2(\tilde{r}-3)}{5\tilde{r}-6}$,

\begin{align*}
\displaystyle  1 \leq\|\rho-1\|_{C\left(\ol{\O }\right)} \le C
\|\rho-1\|_{L^2}^{\theta}\|\na \rho\|_{L^{\tilde{r}}}^{1-\theta},
\end{align*}

Combining this with (<ref>) gives (<ref>) and completes the proof of Theorem <ref>.

Proof of  Theorem   <ref>.  We adopt the idea in [24] to  prove  Theorem   <ref> in two steps.

Step 1. Construction  of approximate  solutions. Assume $(\n_0,u_0,\te_0,H_0)$ satisfying (<ref>) and (<ref>) is the initial data in Theorem <ref> and  $C_0$ satisfies (<ref>)   with  $\ve $  as in  (<ref>).   For $j_{m^{-1}}(x)$ being the standard mollifying kernel of width $m^{-1}$, we construct

\begin{align*}
\hat\n_0^{m,\eta} =\! (\n_0 1_\O)\ast j_{m^{-1}}1_\O\!+\! \eta,\ \  \hat u_0^{m,\eta}=\!\frac{u_0^m }{1\!+\!\eta},\ \  \hat\te_0^{m,\eta}=\! \frac{(\n_0\te_0 1_{\O_m})\ast j_{m^{-1}} \!+\! \eta}{(\n_01_{\O_m})\ast j_{m^{-1}}\!+ \!\eta},\ \  \hat H_0^{m,\eta}=\!\frac{H_0^m }{1\!+\!\eta},
\end{align*}

where  $\O_m=\{x\in\O| dist(x,\p\O)>2/m\}$ and $u_0^{m},H_0^{m}$ satisfies the boundary conditions (<ref>), (<ref>) and

(u_0^m,H_0^m)∈C^∞∩H^1  and  lim_m→∞(u_0^m-u_0,H_0^m-H_0)_H^1=0.

Then for any $\eta\in (0, \eta_0)$ with $\eta_0$  as in (<ref>), there exists $m(\eta)>1$ such that for $m \ge m(\eta)$, the approximate initial data
$(\hat\n_0^{m,\eta},\hat u_0^{m,\eta},\hat\te_0^{m,\eta},\hat{H}_0^{m,\eta})$ satisfies the boundary conditions (<ref>)-(<ref>) and

η≤_0^m,η  <,    η/+η≤_0^m,η ≤,   û_0^m,η_L^2 ≤M,  Ĥ_0^m,η_L^2 ≤M,

and for any $p\geq1$,
lim_η→0 lim_m→∞
( _0^m,η - _0 _ L^p+û_0^m,η-u_0_H^1+_0^m,η_0^m,η-_0_0  _L^2+Ĥ_0^m,η-u_0_H^1)=0
owing to (<ref>)  and (<ref>).

Now, we claim that  the initial norm $\hat C_0^{m,\eta}$
for $(\hat\n_0^{m,\eta},\hat u_0^{m,\eta},\hat\te_0^{m,\eta},\hat{H}_0^{m,\eta}),$  i.e., the right hand side of
(<ref>) with $(\n_0,u_0,\te_0,H_0)$  replaced by
$(\hat\n_0^{m,\eta},\hat u_0^{m,\eta},\hat\te_0^{m,\eta},\hat{H}_0^{m,\eta}),$   satisfies
uv9  lim_η→0
lim_m→∞Ĉ_0^m,η≤C_0,which leads to that there exists an $\hat\eta\in(0,\eta_0)$ such that, for any $\eta\in
(0,\hat\eta ),$ there exists some $\hat m (\eta)\geq m(\eta)$ such that
uv8Ĉ_0^m,η≤C_0+_0/2≤_0 , provided uv01
0<η<η̂,   m≥m̂ (η).Then if we assume (<ref>) holds, it directly follows from Proposition <ref>, (<ref>) and (<ref>) that there exists a classical solution  $(\hat\n^{m,\eta},\hat u^{m,\eta},\hat\te^{m,\eta},\hat H^{m,\eta})  $  of problem (<ref>)–(<ref>) with  initial data $(\hat\n_0^{m,\eta},\hat u_0^{m,\eta},\hat\te_0^{m,\eta},\hat{H}_0^{m,\eta})$   on $\O\times(0,T]$ for all $T>0$.  Furthermore, $(\hat\n^{m,\eta},\hat u^{m,\eta},\hat\te^{m,\eta},\hat H^{m,\eta}) $  satisfies (<ref>), (<ref>), (<ref>), (<ref>), (<ref>),  (<ref>), and (<ref>) respectively  for any $T>0$ and $t\geq1$ with $(\n,u,\te,H)$   replaced by $(\hat\n^{m,\eta},\hat u^{m,\eta},\hat\te^{m,\eta},\hat H^{m,\eta})$.

It remains to prove (<ref>). Indeed, we just need to infer
uv10lim_η→0lim_m→∞∫_0^m,ηΦ(_0^m,η)dx≤∫_0Φ(_0)dx ,since   the other terms in  (<ref>) can be proved in a  similar and even simpler way.
Note that

\begin{align}\notag
\begin{aligned}
\hat\n_0^{m,\eta}\Phi(\hat\te_0^{m,\eta})%=\hat\n_0^{m,\eta}\left(\hat\te_0^{m,\eta}- \log \hat\te_0^{m,\eta} -1\right)\\
&=\hat\n_0^{m,\eta}(\hat\te_0^{m,\eta}-1)^2 \int_0^1\frac{\al}{\al(\hat\te_0^{m,\eta}-1)+1}d\al\\
&= \frac{\hat\n_0^{m,\eta}(j_{m^{-1}}*(\n_0(\te_0-1)1_{\O_m}))^2}{j_{m^{-1}}* (\n_01_{\O_m} )+\eta}\\
&\quad\cdot\int_0^1\frac{\al}{\al j_{m^{-1}}*(\n_0(\te_0-1)1_{\O_m})+j_{m^{-1}}* (\n_01_{\O_m} )+\eta}d\al \\
&\in \left[0, \, \hat\n\eta^{-2}(j_{m^{-1}}*(\n_0(\te_0-1)1_{\O_m}))^2 \right],
\end{aligned}
\end{align}

which combined with  Lebesgue's dominated convergence theorem yields that

\begin{align}\notag
\begin{aligned}
&\lim_{m\rightarrow \infty}\int\hat\n_0^{m,\eta}\Phi(\hat\te_0^{m,\eta})dx\\
- \ln\frac{\n_0\te_0+\eta}{\n_0+\eta}-1 \right)dx\\
%&=\int\left(\n_0\te_0 -\n_0 -(\n_0+\eta)\log\frac{\n_0\te_0+\eta}{\n_0+\eta}\right)dx \\
&=\int\left(\n_0\te_0 -\n_0+(\n_0+\eta)\ln(\n_0+\eta)\right)dx\\
&\le \int\left(\n_0\te_0 -\n_0+(\n_0+\eta)\ln(\n_0+\eta)\right)dx-\int\left(\n_0\ln(\n_0\te_0)+\eta\ln \eta\right)dx\\
& \rightarrow \int \n_0\left(\te_0-  \ln\te_0 -1 \right)dx=\int \n_0\Phi(\te_0)dx,\quad \mbox{ as }\eta\rightarrow 0.
\end{aligned}
\end{align}

It thus gives (<ref>).

Step 2. Compactness results.
With the approximate solutions $(\hat\n^{m,\eta},\hat u^{m,\eta},\hat\te^{m,\eta},\hat H^{m,\eta}) $ obtained in the previous step at hand,  we can derive the global existence of weak solutions by passing to the limit  first   $m\rightarrow \infty,$ then $\eta\rightarrow 0 .$   Since the two steps are similar, we will only   sketch the  arguments for $m\rightarrow \infty.$  For any fixed  $\eta\in (0,\hat\eta)$, we simply denote  $(\hat\n^{m,\eta},\hat u^{m,\eta},\hat\te^{m,\eta},\hat H^{m,\eta}) $
by $(\n^m,u^m,\te^m,H^m).$ Then  the combination of Aubin-Lions Lemma with (<ref>), (<ref>),  (<ref>),  (<ref>),  and Lemma <ref> yields that there exists some appropriate subsequence $ m_j \rightarrow \infty$ of $m\rightarrow \infty$ such that, for any $0<\tau<T<\infty $, $p\in[1,\infty)$, and $\tilde p\in[1,6)$,
u^m_j⇀u, H^m_j⇀H      L^∞(0,T; H^1),
^m_j⇀,  H^m_j⇀H      L^2(0,T;H^1),
^m_j→     C([0,T];L^p)∩C([0,T];H^-1),vu5
^m_ju^m_j→u,  ^m_j^m_j
→     C([0,T];L^2 )∩C([0,T];H^-1),
vu26 ^m_j |u^m_j|^2→|u|^2       C([0,T];L^3 )∩C([0,T];H^-1),ghh
F^m_j→F, ω^m_j→ω    C([τ,T];H^1)∩C([τ,T];L^p̃),
lllu^m_j→u   C([τ,T];W^1,6)∩C(Ø×[τ,T]),and
^m_j→,  H^m_j→H    C([τ,T];H^2)∩C([τ,T];W^1,p̃).

Now we consider the approximate solutions $(\n^{m_j},u^{m_j},\te^{m_j},H^{m_j})$ in the weak forms, i.e. (<ref>)–(<ref>), then  take appropriate limits. Standard arguments as well as (<ref>) and (<ref>)–(<ref>) thus conclude that the limit $(\n,u,\te,H)$ is a weak solution  of   (<ref>),   (<ref>)-(<ref>)  in the sense of Definition <ref> and satisfies (<ref>)–(<ref>) and the exponential decay property (<ref>).
Moreover, we obtain the estimates (<ref>)–(<ref>) with the aid of (<ref>), (<ref>), (<ref>),    and (<ref>)–(<ref>). Finally, (<ref>) shall be obtained by adopting the similar way as in [19].
The proof of Theorem <ref> is finished.

§ ACKNOWLEDGEMENTS
The research  is partially supported by the National Natural Science Foundation of China (No. 11901025).

[1]
S. Agmon, A. Douglis, and L. Nirenberg.
Estimates near the boundary for solutions of elliptic partial
differential equations satisfying general boundary conditions II.
Comm. Pure Appl. Math., 17:35–92, 1964.

[2]
J. Aramaki.
$L^p$ theory for the div-curl system.
Int. J. Math. Anal., 8(6):259–271, 2014.

[3]
G. Cai and J. Li.
Existence and exponential growth of global classical solutions to the
compressible Navier-Stokes equations with slip boundary conditions in 3D
bounded domains.
Indiana Univ. Math. J., in press.

[4]
G. Q. Chen and D. Wang.
Global solutions of nonlinear Magnetohydrodynamics with large
initial data.
J. Differential Equations, 182:344–376, 2002.

[5]
Y. Chen, B. Huang, and X. Shi.
Global strong and weak solutions to the initial-boundary-value
problem of two-dimensional compressible MHD system with large initial data
and vacuum.
SIAM J. Math. Anal., 54(3):3817–3847, 2022.

[6]
Y. Chen, B. Huang, and X. Shi.
Global strong solutions to the compressible magnetohydrodynamic
equations with slip boundary conditions in 3D bounded domains.

[7]
Y. Cho and H. Kim.
Existence results for viscous polytropic fluids with vacuum.
J. Differ. Equations, 228(2):377–411, 2006.

[8]
C. Dou, S. Jiang, and Q. Ju.
Global existence and the low mach number limit for the compressible
magnetohydrodynamic equations in a bounded domain with perfectly conducting
Z. Angew. Math. Phys., 64(6):1661–1678, 2013.

[9]
B. Ducomet and E. Feireisl.
The equations of magnetohydrodynamics: On the interaction between
matter and radiation in the evolution of gaseous stars.
Comm. Math. Phys., 266:595–629, 2006.

[10]
J. Fan and W. Yu.
Strong solution to the compressible magnetohydrodynamic equations
with vacuum.
Nonlinear Anal. Real World Appl., pages 392–409, 2009.

[11]
E. Feireisl.
Dynamics of viscous compressible fluids.
Oxford Science Publication, Oxford, 2004.

[12]
E. Feireisl and Y. Li.
On global-in-time weak solutions to the magnetohydrodynamic system of
compressible inviscid fluids.
Nonlinearity, 33(1):139–155, 2020.

[13]
G. P. Galdi.
An introduction to the mathematical theory of the
Navier-Stokes equations. Steady-state problems.
New York, NY: Springer, 2nd ed. edition, 2011.

[14]
G. Hong, X. Hou, H. Peng, and C. Zhu.
Global existence for a class of large solutions to three-dimensional
compressible magnetohydrodynamic equations with vacuum.
Siam J Math. Anal., 49(4):2409–2441, 2017.

[15]
X. Hou, M. Jiang, and H. Peng.
Global strong solution to 3d full compressible magnetohydrodynamic
flows with vacuum at infinity.
Z. Angew. Math. Phys., 73(1):18, 2022.
Id/No 13.

[16]
X. Hu and D. Wang.
Global solutions to the three-dimensional full compressible
magnetohydrodynamic flows.
Commun. Math. Phys., 283:255–284, 2008.

[17]
X. Hu and D. Wang.
Global existence and large-time behavior of solutions to the
three-dimensional equations of compressible magnetohydrodynamic flows.
Arch. Ration. Mech. Anal., 197:203–238, 2010.

[18]
B. Huang, X. Shi, and Y. Sun.
Large-time behavior of magnetohydrodynamics with
temperature-dependent heat-conductivity.
J. Math. Fluid Mech., 23(3):23, 2021.

[19]
X. Huang and J. Li.
Global classical and weak solutions to the three-dimensional full
compressible Navier-Stokes system with vacuum and large oscillations.
Arch. Ration. Mech. Anal., 227(3):995–1059, 2018.

[20]
X. Huang, J. Li, and Z. Xin.
Global well-posedness of classical solutions with large oscillations
and vacuum to the three-dimensional isentropic compressible Navier-Stokes
Comm. Pure Appl. Math., 65(4):549–585, 2012.

[21]
S. Kawashima.
Smooth global solutions for two-dimensional equations of
electro-magneto-fluid dynamics.
Japan J. Appl. Math., 1:207–222, 1984.

[22]
S. Kawashima and M. Okada.
Smooth global solutions for the one-dimensional equations in
Proc. Japan Acad. Ser. A Math. Sci., 53(9):384–387, 1982.

[23]
H. Li, X. Xu, and J. Zhang.
Global classical solutions to 3D compressible magnetohydrodynamic
equations with large oscillations and vacuum.
SIAM J. Math. Anal., 45:1356–1387, 2013.

[24]
J. Li, L. Boqiang, and W. Xue.
Global existence of classical solutions to full compressible
Navier-Stokes system with large oscillations and vacuum in 3D bounded

[25]
J. Li and Z. Xin.
Some uniform estimates and blowup behavior of global strong solutions
to the Stokes approximation equations for two-dimensional compressible
J. Differ. Equations, 221(2):275–308, 2006.

[26]
J. Li and Z. Xin.
Global Existence of Regular Solutions with Large Oscillations
and Vacuum.
In: Giga Y., Novotny A. (eds) Handbook of Mathematical Analysis in
Mechanics of Viscous Fluids. Springer, 2016.

[27]
Y. Li and Y. Sun.
Global weak solutions to a two-dimensional compressible MHD
equations of viscous non-resistive fluids.
J. Differ. Equations, 267(6):3827–3851, 2019.

[28]
Y. Li and Y. Sun.
On global-in-time weak solutions to a two-dimensional full
compressible nonresistive MHD system.
SIAM Journal on Mathematical Analysis, 53(4):4142–4177, 2021.

[29]
H. Liu, T. Luo, and H. Zhong.
Global solutions to an initial boundary problem for the compressible
3-D MHD equations with Navier-slip and perfectly conducting boundary
conditions in exterior domains.

[30]
S. Liu, H. Yu, and J. Zhang.
Global weak solutions of 3D compressible MHD with discontinuous
initial data and vacuum.
J. Differential Equations, 254:229–255, 2013.

[31]
Y. Liu and X. Zhong.
Global well-posedness to three-dimensional full compressible
magnetohydrodynamic equations with vacuum.
Z. Angew. Math. Phys., 71(6):1–25, 2020.

[32]
Y. Liu and X. Zhong.
Global existence and decay estimates of strong solutions for
compressible non-isentropic magnetohydrodynamic flows with vacuum.

[33]
Y. Liu and X. Zhong.
Global strong solution for 3D compressible heat-conducting
magnetohydrodynamic equations revisited.

[34]
L. Lu and B. Huang.
On local strong solutions to the Cauchy problem of the
two-dimensional full compressible magnetohydrodynamic equations with vacuum
and zero heat conduction.
Nonlinear Anal., Real World Appl., 31:409–430, 2016.

[35]
B. Lv and B. Huang.
On strong solutions to the Cauchy problem of the two-dimensional
compressible magnetohydrodynamic equations with vacuum.
Nonlinearity, 28(2):509–530, 2015.

[36]
B. Lv, X. Shi, and X. Xu.
Global existence and large-time asymptotic behavior of strong
solutions to the compressible magnetohydrodynamic equations with vacuum.
Indiana Univ. Math. J., 65(3):925–975, 2016.

[37]
C. Navier.
Mémoire sur les lois du mouvement des fluides.
Mém. Acad. Re. Sci. Inst. Fr., 6:389–416, 1823.

[38]
L. Nirenberg.
On elliptic partial differential equations.
Annali della Scuola Normale Superiore di Pisa, 13(2):115–162,

[39]
A. Suen and D. Hoff.
Global low-energy weak solutions of the equations of
three-dimensional compressible magnetohydrodynamics.
Arch. Rational Mech. Anal., 205:27–58, 2012.

[40]
Z. Tan and Y. Wang.
Global well-posedness of an initial-boundary value problem for
viscous non-resistive MHD systems.
SIAM J. Math. Anal., 50(1):1432–1470, 2018.

[41]
T. Tang and H. Gao.
Strong solutions to 3D compressible magnetohydrodynamic equations
with Navier-slip condition.
Math. Methods Appl. Sci., 39(10):2768–2782, 2016.

[42]
W. von Wahl.
Estimating $\nabla u$ by ${\rm div}\, u$ and ${\rm curl}\, u$.
Math. Methods Appl. Sci., 15(2):123–143, 1992.

[43]
R. W.
Real and complex analysis, Third Edition.
McGraw-Hill Book Company, New York, 1987.

[44]
D. Wang.
Large solutions to the initial-boundary value problem for planar
SIAM J. Math. Anal., 63(4):1424–1441, 2003.

[45]
J. Wu and Y. Wu.
Global small solutions to the compressible 2D magnetohydrodynamic
system without magnetic diffusion.
Adv. Math., 310:759–888, 2017.

[46]
S. Xi and X. Hao.
Existence for the compressible magnetohydrodynamic equations with
J. Math. Anal. Appl., 453:410–433, 2017.

[47]
S. Zhu.
On classical solutions of the compressible magnetohydrodynamic
equations with vacuum.
SIAM J. Math. Anal., 47(4):2722–2753, 2015.