FLUID-RIGID BODY INTERACTION IN A COMPRESSIBLE ELECTRICALLY CONDUCTING FLUID
JAN SCHERZ1,2,3

  * 1

Department of Mathematical Analysis, Faculty of Mathematics and Physics,
Charles University in Prague, Sokolovská 83, Prague 8, 18675, Czech Republic

  * 2

Mathematical Institute, Academy of Sciences, Žitná 25, Prague 1, 11567, Czech
Republic

  * 3

Institute of Mathematics, University of Würzburg, Emil-Fischer-Str. 40, 97074
Würzburg, Germany

Abstract

We consider a system of multiple insulating rigid bodies moving inside of an
electrically conducting compressible fluid. In this system we take into
account the interaction of the fluid with the bodies as well as with the
electromagnetic fields trespassing both the fluid and the solids. The main
result of this article yields the existence of weak solutions to the system.
While the mechanical part of the problem can be dealt with via a classical
penalization method, the electromagnetic part requires an approximation by
means of a hybrid discrete-continuous in time system: The discrete part of the
approximation enables us to handle the solution-dependent test functions in
our variational formulation of the induction equation, whereas the continuous
part makes sure that the non-negativity of the density and subsequently a
meaningful energy inequality is preserved in the approximate system.

##  1 Introduction

The goal of this article is the proof of the existence of weak solutions to a
system of PDEs modelling the motion of multiple non-conducting rigid objects
inside of an electrically conducting compressible fluid taking the interplay
with the electromagnetic fields in these materials into account. The problem
can be regarded as belonging to both the research areas of
magnetohydrodynamics (MHD) and fluid-structure interaction (FSI). Indeed, MHD
(c.f. [5, 8, 35]), on the one hand, describes a coupling between the Maxwell
system and the Navier-Stokes equations which models the influence of the
electromagnetic fields on the electrically conducting fluid and vice versa.
FSI, on the other hand, models the interaction between the fluid and the solid
bodies through a coupling between the Navier-Stokes system and the balance of
mass and momentum of the rigid bodies.

The interplay between the electromagnetic fields and the solids occurs
implicitly via the fluid due to the non-conductivity of the solid material.
This paper represents the extension of [1], where the corresponding problem
was studied and solved for incompressible fluids, to the compressible case.
Potential applications can be found in the area of biomechanics. A specific
example is constituted by capsule endoscopy, a medical procedure during which
small cameradevices are sent through the body with the aim of detecting
diseases. Such endocapsules of a microscopic scale can be navigated through
the (electrically conducting) blood by controlling their movement via the
application of electromagnetic forces, c.f. [6, Section 4.4], [4]. This
technique can also be applied in the problem of drug delivery, in which
microrobots are constructed to deliver drugs directly to the targeted area of
the body without affecting healthy tissue. However, in view of the usage of
electrically conducting microrobots for these purposes an extension of the
current model to electrically conducting rigid bodies would be needed.

In order to classify our result within the wide range of works in MHD and FSI,
we give a brief summary of the related literature. The existence of weak
solutions to the MHD problem for a compressible fluid but without any rigid
bodies involved is for example shown in [41] and [2]. The fluids considered in
the latter one of these articles are electrically as well as thermally
conducting and moreover the existence of strong solutions is addressed
therein. The existence of weak solutions to the MHD problem in the
incompressible case is proved in [25].

On the FSI side of the problem we first mention [23, 42] for a general
introduction to the fluid-rigid body interaction problem. Early existence
results for weak solutions to this problem were obtained mainly in the
incompressible case and include the articles [7, 9, 29, 32], wherein the proof
of local in time existence, i.e. existence up to a contact between several
bodies or a body and the domain boundary, is achieved in the case of two and
three spatial dimensions. The corresponding problem in the compressible case
was considered in [10]. A proof of the global in time existence of weak
solutions to a model describing the interaction between multiple rigid bodies
and a compressible fluid is given in [17]. Corresponding results are also
available for incompressible fluids, c.f. [40] for the two-dimensional and
[18] for the three-dimensional case. The article [40] in particular touches
upon the question whether contacts between the bodies with each other or the
domain boundary are possible and shows that such collisions can only occur if
the relative acceleration and velocity between the colliding objects vanish.
Moreover, the problem has been studied in the case of the Navier-slip instead
of the classical no-slip boundary condition, c.f. [36], and also the question
about the existence of strong solutions has been investigated in both the
compressible, c.f. [3, 30, 31, 39], and the incompressible case, c.f. [24, 44,
45].

The articles [27] (for the case of two spatial dimensions) and [28] (for the
$3$D case) can be considered as a first step towards the coupling between the
MHD and the FSI problem. The authors thereof studied a model of an
incompressible electrically conducting fluid flowing around a fixed non-
conducting solid region. This model served as the basis for the
electromagnetic part of the model in [1], wherein Benešová, Nečasová,
Schlömerkemper and the author of this article showed the (local in time)
existence of weak solutions to the problem of one movable insulating rigid
object travelling through an electrically conducting incompressible fluid. The
main goal of the present article is to extend the latter result to the case of
a compressible fluid. More precisely, we are able to prove the global in time
existence of weak solutions to the interaction problem of finitely many
insulating rigid bodies, an electrically conducting compressible fluid and the
electromagnetic fields trespassing these materials, see Section 2.3.

As in the incompressible case in [1], the main difficulty in the proof of the
existence result is caused by the test functions in the weak formulation of
the induction equation, c.f. (22), (33) below, which are chosen such that they
depend on the solid region. While such test functions do not generate any
problems in the case of an immovable solid region (see e.g. the proofs of [27,
Theorem 2.1] and [28, Theorem 2.3]), difficulties arise in our scenario, where
the solid domain depends on the overall solution to the system, which causes
our problem to be highly coupled. Following the proof in the incompressible
situation, we thus make use of a time discretization, which allows us to deal
with this problem by decoupling it: At each discrete time we first calculate
the domain of the solid bodies, which suggests a suitable definition for the
test functions in the induction equation at that specific time. Only after
this we solve the induction equation itself, which can then be achieved via
standard methods. In the compressible situation, however, this procedure turns
out to cause more problems. In particular, the author could not find a
suitable way to discretize the Navier-Stokes system while preserving the non-
negativity of the density. This is essential to obtain the uniform bounds from
the energy inequality required to pass to the limit in the approximate system.
We handle this problem by choosing a hybrid approximation system, in which the
induction equation is discretized in time via the Rothe method ([38, Section
8]), whereas the mechanical part of the system is studied as a time-dependent
problem on the small intervals between the discrete time points. The non-
negativity of the density can then be derived by classical arguments and, by
choosing the coupling terms in a suitable way, the discrete electromagnetic
part and the continuous mechanical part of the system can be combined to an
energy inequality with all desired features. To a smaller extent, we already
used such a hybrid approximation in the incompressible situation [1], where,
however, the time dependency was restricted to the transport equation for the
characteristic function of the solid region. The expansion of this idea to the
whole mechanical part of the system, in order to deal with the problems
outlined above, is what constitutes the main novelty in the proof of our main
result.

The problem of the solution dependent test functions also appears in the weak
formulation of the momentum equation, c.f. (21), (32) below. In this
situation, however, we have the penalization method used e.g. in [17] and [40]
readily available which allows us to evade the problem. More precisely, in
this penalization method a sequence of approximate solutions to some fluid-
only problems with classical test functions is constructed. Passing to the
limit in this approximation one then returns to a fluid-rigid body interaction
system by letting the viscosity of the fluid rise to infinity in the later
solid region.

The paper is structured as follows: We begin by summarizing the notation
needed for our model in Section 1.1 and subsequently present the model,
divided into a mechanical and an electromagnetic part, in Sections 1.2 and
1.3. After introducing some additional notation in Section 2.1, we present the
variational formulation of the above model in Section 2.2 as well as our main
result in Section 2.3. In Section 3 we explain the main ideas for the proof of
this result, which is based on an approximation of the original system. In
Section 4 we solve this approximate problem and finally, in the remaining
Sections 5–9, we pass to the limit in the approximation, proving the existence
of a weak solution to the original system.

###  1.1 Model

The model we consider describes several non-conducting rigid bodies travelling
through an electrically conducting compressible fluid as well as the involved
electromagnetic fields. This model is a combination of (i) the mechanical
fluid-rigid body interaction model used in [17] and (ii) the Maxwell system in
the model used in [27, 28] for the description of an electrically conducting
fluid surrounding an immovable, non-conducting solid region. It is further the
extension of the corresponding model for the incompressible case in [1] to the
compressible situation. Let $T>0$ and let $\Omega\subset\mathbb{R}^{3}$ be a
bounded domain. Inside of $\Omega$ we consider $m\in\mathbb{N}$ insulating
rigid bodies, the positions of which at time $t\in[0,T]$ are described through
subsets $S^{i}(t)\subset\Omega$, $i=1,...,m$. The complement of the solid
domain,

$\displaystyle F(t):=\Omega\setminus\overline{S}(t),\quad\quad\text{with}\quad
S(t):=\bigcup_{i=1}^{m}S^{i}(t),$

contains an electrically conducting viscous nonhomogeneous compressible fluid.
We denote by $Q$ the time-space domain $Q:=(0,T)\times\Omega$, which we split
into a solid part $Q^{s}$ and a fluid part $Q^{f}$,

$\displaystyle Q^{s}=Q^{s}\left(S\right):=\left\\{(t,x)\in Q:\ x\in
S(t)\right\\},\quad\quad Q^{f}=Q^{f}(S):=\left\\{(t,x)\in Q:\ x\in
F(t)\right\\}.$ (1)

For any function defined on $Q$ we mark its restriction to $Q^{f}$ or $Q^{s}$
by the superscript $f$ or $s$, respectively.

Figure 1: A domain filled with an electrically conducting compressible fluid
and three insulating rigid bodies.

The interaction between the fluid, the solids and the electromagnetic fields
in the domain is characterized through the mass density
$\rho:Q\rightarrow\mathbb{R}$, the velocity field
$u:Q\rightarrow\mathbb{R}^{3}$, the magnetic induction
$B:Q\rightarrow\mathbb{R}^{3}$, the electric field
$E:Q\rightarrow\mathbb{R}^{3}$ and the electric current
$j:Q\rightarrow\mathbb{R}^{3}$. As indicated above, our overall model, which
determines these functions, can be divided into two subsystems: The mechanical
subsystem for the description of $\rho$ and $u$ and the electromagnetic
subsystem for the description of $B$, $E$ and $j$.

###  1.2 Mechanical subsystem

The mechanical quantities evolve according to the compressible Navier-Stokes
equations in the fluid domain and the balance of linear and angular momentum
of the rigid bodies in the solid region, respectively,

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)=0\quad\quad$
$\displaystyle\text{in }Q^{f},$ (2) $\displaystyle\partial_{t}(\rho
u)+\text{div}(\rho u\otimes u)+\nabla p=\text{div}\mathbb{T}+\rho
g+\frac{1}{\mu}\text{curl}B\times B\quad\quad$ $\displaystyle\text{in }Q^{f},$
(3) $\displaystyle m^{i}\frac{d}{dt}V^{i}(t)=\frac{d}{dt}\int_{S^{i}(t)}\rho
u\ dx=\int_{\partial S^{i}(t)}\left[\mathbb{T}-p\text{Id}\right]\cdot\text{n}\
d\sigma+\int_{S^{i}(t)}\rho g\ dx,$ $\displaystyle\quad t\in[0,T],\
i=1,...,m,$ (4)
$\displaystyle\frac{d}{dt}\left(\mathbb{J}^{i}(t)w^{i}(t)\right)=\frac{d}{dt}\int_{S^{i}(t)}\rho\left(x-h^{i}\right)\times
u\ dx\ $ $\displaystyle=\int_{\partial
S^{i}(t)}(x-h^{i})\times\left[\mathbb{T}-p\text{Id}\right]\text{n}\
d\sigma+\int_{S^{i}(t)}\rho\left(x-h^{i}\right)\times g\ dx,$
$\displaystyle\quad t\in[0,T],\ i=1,...,m,$ (5)

combined with the boundary and interface conditions

$\displaystyle u(t)$ $\displaystyle=0\quad\text{on }\partial\Omega,\quad\quad$
$\displaystyle\ \ \ \ \ \ \ \ \ \ u^{f}(t)-u^{s}(t)=0\quad\text{on }\partial
S(t).$ (6)

The identity (2) is known as the continuity equation. In the momentum equation
(3) we see the pressure $p$, for which we assume the isentropic constitutive
relation

$\displaystyle p=p\left(\rho^{f}\right)=a\left(\rho^{f}\right)^{\gamma},\quad
a>0,\ \gamma>\frac{3}{2},$

and the stress tensor

$\displaystyle\mathbb{T}=\mathbb{T}(u):=2\nu\mathbb{D}(u)+\lambda\text{Id}\
\operatorname{div}u,\quad\quad\mathbb{D}(u):=\frac{1}{2}\nabla
u+\frac{1}{2}(\nabla u)^{T}$

with the viscosity coefficients $\nu,\lambda\in\mathbb{R}$ which satisfy
$\nu>0,\ \lambda+\nu\geq 0$. Moreover, we have two forcing terms: The external
force $g$ and the Lorentz force $\frac{1}{\mu}\text{curl}B\times B$ with the
magnetic permeability $\mu>0$. We assume that

$\displaystyle\mu\ \text{is constant in the whole domain}\ Q.$

Since, in general, the magnetic permeability takes different values in
conducting and insulating materials, this assumption is physically
restrictive. However, it is necessary for the transition conditions on the
magnetic induction $B$, c.f. Section 1.3 below. The Lorentz force constitutes
the connection of the mechanical to the electromagnetic subsystem. In the
balance of mass (4) and momentum (5) of the rigid bodies instead the Lorentz
force does not appear, which is in accordance with the assumption of the
bodies being non-conducting. The notation used in these relations includes the
total mass $m^{i}$ of the $i$-th body, its center of mass $h^{i}$ and the
associated inertia tensor $\mathbb{J}^{i}$,

$\displaystyle m^{i}:=\int_{S^{i}(t)}\rho(t,x)\ dx,\quad
h^{i}(t):=\frac{1}{m^{i}}\int_{S^{i}(t)}\rho(t,x)x\ dx,\quad\quad$
$\displaystyle t\in[0,T],$ $\displaystyle\mathbb{J}^{i}(t)a\cdot
b:=\int_{S^{i}(t)}\rho(t,x)\left[a\times\left(x-h^{i}(t)\right)\right]\cdot\left[b\times\left(x-h^{i}(t)\right)\right]\
dx,\quad\quad$ $\displaystyle t\in[0,T],\quad a,b\in\mathbb{R}^{3}.$

The equations (4) and (5) determine the translational and rotational
velocities $V^{i}$ and $w^{i}$ of the $i$-th rigid body, respectively,
allowing us to express its total velocity as a rigid velocity field

$\displaystyle
u(t,x)=u^{s^{i}}(t,x):=V^{i}(t)+w^{i}(t)\times\left(x-h^{i}(t)\right)\quad\quad\text{for
}t\in[0,T],\ x\in S^{i}(t).$

The coupling between the fluid and the solids is incorporated into the surface
integrals in (4) and (5) and into the no-slip interface condition in (6).
Indeed, on the one hand, the presence of the Cauchy-stress
$\mathbb{T}-p\text{Id}$ of the fluid in (4) and (5) shows how the velocity and
the pressure of the fluid affect the motion of the bodies. On the other hand,
considering the Navier-Stokes system (2), (3), (6) by itself, one can regard
the interface condition in (6) as a part of the boundary conditions, which
describes the impact of the velocity of each solid on the fluid velocity. The
first identity in (6) is the standard no-slip boundary condition.

###  1.3 Electromagnetic subsystem

The electromagnetic quantities are determined by the following version of the
Maxwell system,

$\displaystyle\frac{1}{\mu}\text{curl}B$
$\displaystyle=\left\\{\begin{array}[]{ll}j+J&\text{in }Q^{f}\\\ 0&\text{in
}Q^{s}\end{array}\right.,$ (9) $\displaystyle\partial_{t}B+\text{curl}E$
$\displaystyle=0\quad\quad\quad\quad\quad\text{in }Q^{f}\ \text{and }Q^{s},$
(10) $\displaystyle\text{div}E$
$\displaystyle=0\quad\quad\quad\quad\quad\text{in }Q^{s},$ (11)
$\displaystyle\text{div}B$ $\displaystyle=0\quad\quad\quad\quad\quad\text{in
}Q^{f}\ \text{and }Q^{s},$ (12)

together with Ohm’s law

$\displaystyle j=\sigma(E+u\times B)\quad\text{in
}Q,\quad\sigma=\left\\{\begin{array}[]{ll}\sigma^{f}>0&\text{in }Q^{f}\\\
\sigma^{s}=0&\text{in }Q^{s}\end{array}\right.$ (15)

and the boundary and interface conditions

$\displaystyle B(t)\cdot\text{n}$ $\displaystyle=0\quad\text{on
}\partial\Omega,\quad\quad$ $\displaystyle\ \ \ \ \ \ \ \ \
B^{f}(t)-B^{s}(t)=0\quad\text{on }\partial S(t),$ (16) $\displaystyle
E(t)\times\text{n}$ $\displaystyle=0\quad\text{on }\partial\Omega,\quad\quad$
$\displaystyle\left(E^{f}(t)-E^{s}(t)\right)\times\text{n}=0\quad\text{on
}\partial S(t).$ (17)

In this system we have Ampère’s law (9), the Maxwell-Faraday equation (10),
Gauss’s law (11) and Gauss’s law for magnetism (12). In comparison to the
general situation, these equations have undergone two kinds of reductions,
c.f. [27], [28]: First, in the solid region, the system is adjusted to the
assumption of the rigid bodies being insulating and second, in the fluid
domain, the system is reduced according to the magnetohydrodynamic
approximation, see [8, 13]. A justification for the latter simplification,
which is inherent to magnetohydrodynamics, from a physical point of view is
for example given in [33, 34]. The mathematical use of the magnetohydrodynamic
approximation consists of the fact that it enables us to summarize the whole
electromagnetic problem in $Q^{f}$ into a problem for the magnetic induction
$B$, c.f. the induction equation (33) in our definition of variational
solutions below. Once $B$ has been determined, the remaining unknowns $E$ and
$j$ are given explicitly through the relations (9) and (15). The quantity $J$
in (9) represents, as in [27], [28], a supplementary external force. Ohm’s law
(15) is what links the electromagnetic subsystem to the mechanical subsystem
(2)–(6) by describing the influence of the (fluid) velocity $u=u^{f}$ on the
electromagnetic quantities in $Q^{f}$. Since the electrical conductivity
$\sigma$ satisfies $\sigma=\sigma^{s}=0$ in $Q^{s}$, it also shows that the
electromagnetic fields are not affected by the solid velocities $u^{s^{i}}$.
In the boundary and interface conditions (16), (17) the assumption of the
magnetic permeability $\mu>0$ being constant in the whole domain $Q$ (c.f.
Section 1.2) comes into play. Indeed, while the boundary condition for $B$ in
(16) as well as the conditions (17) on $E$ are standard, the continuity of $B$
across the interface in (16) is not. However, it is standard to assume
continuity of the normal component of $B$ and of the tangential component of
$\frac{1}{\mu}B$ and hence, when $\mu$ is constant across $\partial S(t)$, we
infer also the second relation in (16). The reason why we cannot allow $B$ to
have a jump across the interface is because otherwise it could not be a
Sobolev function in $\Omega$, c.f. (28) below.
Both the mechanical subsystem (2)–(6) and the electromagnetic subsystem
(9)–(17) can be studied in their own right, if one considers $B$ and $u$ as a
prescribed external forcing term in Ohm’s law (15) and the momentum equation
(3), respectively. In this article, however, we study the combined system
(2)–(17), in which both $B$ and $u$ are regarded as unknowns, coupled via
Ohm’s law and the presence of the Lorentz force
$\frac{1}{\mu}\operatorname{curl}B\times B$ in the momentum equation.

##  2 Variational formulation and main result

###  2.1 Notation

The initial positions of the solid bodies are characterized through sets
$S_{0}^{i}\subset\Omega$, $i=1,...,m$, onto which we impose the conditions

$\displaystyle\emptyset\neq S_{0}^{i}\ \text{is open and
connected,}\quad|\partial
S_{0}^{i}|=0,\quad\overline{S}_{0}^{i}\bigcap\overline{S}_{0}^{j}=\emptyset\quad\quad\forall
i,j=1,...,m,\ i\neq j.$ (18)

Since the motion of the bodies is rigid, we can associate to each body an
isometry $X^{i}(t,\cdot):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$,
$t\in[0,T]$, such that its position $S^{i}(t)$ at an arbitrary time
$t\in[0,T]$ can be expressed through the set-valued function

$\displaystyle S^{i}:\ [0,T]\rightarrow 2^{\mathbb{R}^{3}},\quad
S^{i}(t):=X^{i}\left(t;S_{0}^{i}\right).$

In particular, with the notation

$\displaystyle X(t;\cdot):\
S_{0}:=\bigcup_{i=1}^{m}S_{0}^{i}\rightarrow\mathbb{R}^{3},\quad
X(t;\cdot)|_{S^{i}_{0}}:=X^{i}(t;\cdot)\quad\quad\forall i=1,...,m,\
t\in[0,T],$

the solid region at the time $t$ is given by $S(t):=X(t;S_{0})$. In order to
connect the motion of the bodies to the velocity field $u$ we require $u$ to
be compatible with the system $\\{S_{0}^{i},X^{i}\\}_{i=1}^{m}$, i.e. we
require the existence of rigid velocity fields $u^{s^{i}}(t,\cdot)$,
$i=1,...,m$, such that

$\displaystyle u(t,x)=u^{s^{i}}(t,x)\quad\text{for a.a. }t\in[0,T]\ \text{and
a.a. }x\in S^{i}(t)$ (19)

and $X^{i}$ is the unique Carathéodory solution (c.f. [38, Theorem 1.45]) to
the initial value problem

$\displaystyle\frac{dX^{i}(t;x)}{dt}=u^{s^{i}}\left(t;X^{i}(t;x)\right),\quad
X^{i}(0;x)=x,\quad x\in\mathbb{R}^{3}.$ (20)

Finally we denote by $T(S)$ and $Y(S)$ the test function spaces

$\displaystyle T(S):=$ $\displaystyle\left\\{\phi\in\mathcal{D}(Q):\
\mathbb{D}(\phi)=0\ \text{in a neighbourhood of
}\overline{Q}^{s}(S)\right\\},$ (21) $\displaystyle Y(S):=$
$\displaystyle\left\\{b\in\mathcal{D}(Q):\ \operatorname{curl}b=0\ \text{in a
neighbourhood of }\overline{Q}^{s}(S)\right\\}$ (22)

for our variational formulations of the momentum equation and the induction
equation, respectively, in Definition 2.1 below.

###  2.2 Weak solutions

We are now in the position to present our variational formulation of the
system (2)–(17). With a slight abuse of notation we will write here and in the
following sections $\sigma=\sigma^{f}>0$, since the quantities containing
$\sigma^{s}=0$ are not visible in this weak formulation due to the non-
conductivity of the solids.

###### Definition 2.1

Let $T>0$, let $\Omega\subset\mathbb{R}^{3}$ be a bounded domain and let
$S_{0}=\bigcup_{i=1}^{m}S_{0}^{i}$, where $S_{0}^{i}\subset\Omega$ for
$i=1,...,m\in\mathbb{N}$ satisfy the conditions (18). Assume
$\nu,\lambda,a,\gamma,\sigma,\mu\in\mathbb{R}$ to satisfy

$\displaystyle\nu,a,\sigma,\mu>0,\quad\nu+\lambda\geq
0,\quad\gamma>\frac{3}{2}.$ (23)

Moreover, consider some external forces $g,J\in L^{\infty}((0,T)\times\Omega)$
and some initial data $0\leq\rho_{0}\in L^{\gamma}(\Omega)$, $(\rho u)_{0}\in
L^{1}(\Omega)$, $B_{0}\in L^{2}(\Omega)$ such that

$\displaystyle\frac{\left|(\rho u)_{0}\right|^{2}}{\rho_{0}}\in
L^{1}(\Omega),\quad(\rho u)_{0}=0\ \ \text{a.e. in }\left\\{x\in\Omega:\
\rho_{0}(x)=0\right\\},\quad\operatorname{div}B_{0}=0\ \ \text{in
}\mathcal{D}^{\prime}(\Omega),\quad\left.B_{0}\cdot\text{n}\right|_{\partial\Omega}=0.$
(24)

Then the system (2)–(17) is said to admit a weak solution if there exists a
function

$\displaystyle X:\ [0,T]\times
S_{0}\rightarrow\mathbb{R}^{3},\quad\left.X(t;\cdot)\right|_{S_{0}^{i}}=X^{i}(t;\cdot)\quad\quad\forall
i=1,...,m,\ t\in[0,T],$ (25)

where each $X^{i}(t;\cdot):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ denotes an
isometry, and if there exist functions

$\displaystyle 0\leq\rho$ $\displaystyle\in
L^{\infty}\left(0,T;L^{\gamma}(\Omega;\mathbb{R})\right)\bigcap
C\left([0,T];L^{1}(\Omega;\mathbb{R})\right),$ (26) $\displaystyle u$
$\displaystyle\in\left\\{\phi\in
L^{2}\left(0,T;H_{0}^{1,2}\left(\Omega;\mathbb{R}^{3}\right)\right):\
\mathbb{D}(\phi)=0\ \text{in }Q^{s}(S)\right\\},$ (27) $\displaystyle B$
$\displaystyle\in\left\\{b\in
L^{\infty}\left(0,T;L^{2}\left(\Omega;\mathbb{R}^{3}\right)\right)\bigcap
L^{2}\left(0,T;H_{\operatorname{div}}^{1,2}\left(\Omega;\mathbb{R}^{3}\right)\right):\
\operatorname{curl}b=0\ \text{in }Q^{s}(S),\
b\cdot\operatorname{n}|_{\partial\Omega}=0\right\\},$ (28)

where $S=S(\cdot)=X(\cdot,S_{0})$, such that $\rho$ and $u$, extended by $0$
in $\mathbb{R}^{3}\setminus\Omega$, satisfy the continuity equation and its
renormalized form,

$\displaystyle\partial_{t}\rho+\operatorname{div}\left(\rho u\right)=0\quad$
$\displaystyle\text{in
}\mathcal{D}^{\prime}\left((0,T)\times\mathbb{R}^{3}\right),$ (29)
$\displaystyle\partial_{t}\zeta(\rho)+\operatorname{div}\left(\zeta\left(\rho\right)u\right)+\left[\zeta^{\prime}\left(\rho\right)\rho-\zeta\left(\rho\right)\right]\operatorname{div}u=0\quad$
$\displaystyle\text{in
}\mathcal{D}^{\prime}\left((0,T)\times\mathbb{R}^{3}\right),$ (30)

for any

$\displaystyle\zeta\in
C^{1}\left([0,\infty)\right):\quad\left|\zeta^{\prime}(r)\right|\leq
cr^{\lambda_{1}}\quad\forall r\geq 1,\quad\quad\text{where }c>0,\
\lambda_{1}>-1,$ (31)

such that the momentum equation and the induction equation,

$\displaystyle-\int_{0}^{T}\int_{\Omega}\rho u\cdot\partial_{t}\phi\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left(\rho u\otimes
u\right):\mathbb{D}(\phi)+a\rho^{\gamma}\operatorname{div}\phi-2\nu\mathbb{D}(u):\mathbb{D}(\phi)$
$\displaystyle-\lambda\operatorname{div}u\operatorname{div}\phi+\rho
g\cdot\phi+\frac{1}{\mu}\left(\operatorname{curl}B\times B\right)\cdot\phi\
dxdt,$ (32) $\displaystyle-\int_{0}^{T}\int_{\Omega}B\cdot\partial_{t}b\
dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}B+u\times
B+\frac{1}{\sigma}J\right]\cdot\operatorname{curl}b\ dxdt,$ (33)

are satisfied for any $\phi\in T(S)$ and any $b\in Y(S)$, such that the
initial conditions

$\displaystyle\rho(0)=\rho_{0},\quad\quad\left(\rho u\right)(0)=\left(\rho
u\right)_{0},\quad\quad B(0)=B_{0},$ (34)

hold true, where the latter two equations are to be understood in the sense
that

$\displaystyle\lim_{t\rightarrow 0+}\int_{\Omega}\rho(t,x)u(t,x)\cdot\phi(x)\
dx=\int_{\Omega}(\rho u)_{0}(x)\cdot\phi(x)\ dx,\quad\lim_{t\rightarrow
0+}\int_{\Omega}B(t,x)\cdot b(x)\ dx=\int_{\Omega}B_{0}(x)\cdot b(x)\ dx$ (35)

for all $\phi,b\in\mathcal{D}(\Omega)$ with $\mathbb{D}(\phi)=0$ and
$\operatorname{curl}b=0$ in a neighbourhood of $\overline{S}_{0}$,
respectively, and, finally, such that the system
$\\{S_{0}^{i},X^{i}\\}_{i=1}^{m}$ is compatible with the velocity field $u$.

In this definition the compatibility of the velocity field $u$ and the system
of rigid bodies leads to some vivid consequences for the solids. First of all,
while the bodies are able to touch each other or the domain boundary, the
possibility of interpenetrations are ruled out, c.f. [17, Lemma 3.1, Corollary
3.1]. Moreover, even though the density does not satisfy a transport equation
in the case of a compressible fluid, it still travels along the
characteristics of $u$ in the solid part of the domain, c.f. [17, Lemma 3.2].
For definiteness we present these results in the following lemma.

######  Lemma 2.1 ([17])

Let $\Omega\subset\mathbb{R}^{3}$ and $S_{0}^{i}\subset\mathbb{R}^{3}$,
$i=1,...,m\in\mathbb{N}$ be bounded domains and let further $u\in
L^{2}(0,T;H_{0}^{1,2}(\Omega))$ be extended by $0$ outside of $\Omega$.
Moreover, assume $u$ to be compatible with the system
$\\{S_{0}^{i},X^{i}\\}_{i=1}^{m}$ where each
$X^{i}(t):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, $t\in[0,T]$, $i=1,...,m$,
denotes an isometry. Then it holds:

  * (i)

If, for $i\neq j\in\\{1,...,m\\}$, there exists $\tau\in[0,T]$ such that
$S^{i}(\tau)\bigcap S^{j}(\tau)\neq\emptyset$ then $X^{i}(t)=X^{j}(t)$ for all
$t\in[0,T]$. Further, if there exists $\tau\in[0,T]$ such that
$S^{i}(\tau)\not\subset\Omega$, then $X^{i}(t)=\operatorname{Id}$ for all
$t\in[0,T]$.

  * (ii)

If $\rho\in L^{\infty}(0,T;L^{\gamma}(\Omega))$, $\gamma>1$, extended by $0$
outside of $\Omega$, satisfies

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)=0\quad\quad\text{in
}\mathcal{D}^{\prime}\left((0,T)\times\mathbb{R}^{3}\right)$

then

$\displaystyle\rho\left(t,X^{i}(t;x)\right)=\rho\left(0,x\right)\quad\quad\text{for
all }t\in[0,T],\quad i=1,...,m\quad\text{and a.a. }x\in S_{0}^{i}.$ (36)

A detailed proof of the assertions (i) and (ii) is given in [17, Lemma 3.1,
Corollary 3.1] and [17, Lemma 3.2], respectively. The first part of assertion
(i) can be shown directly from the fact that $\\{S_{0}^{i},X^{i}\\}$ and
$\\{S_{0}^{j},X^{j}\\}$ are compatible with the same velocity field $u$, the
second part then follows by regarding $\mathbb{R}^{3}\setminus\Omega$ as a
rigid body with the associated rigid velocity field $0$. The proof of the
assertion (ii) is achieved via a regularization of $\rho$ with respect to the
spatial variable and a subsequent application of the regularization method by
DiPerna and Lions, c.f. [11], to the continuity equation (36) on compact
subsets of the solid time-space domain.

###  2.3 Main result

Our main result yields existence of the weak solutions as introduced in
Definition 2.1.

###### Theorem 2.1

Let $T>0$, assume $\Omega\subset\mathbb{R}^{3}$ to be a simply connected
domain of class $C^{2}\bigcup C^{0,1}$ and assume $S_{0}^{i}\subset\Omega$,
$i=1,...,m\in\mathbb{N}$ to be domains of class $C^{2}\bigcup C^{0,1}$ which
satisfy the conditions (18). Assume moreover the coefficients
$\sigma,\mu,\nu,\lambda,a,\gamma\in\mathbb{R}$ to satisfy the conditions (23)
and the data $g,J\in L^{\infty}((0,T)\times\Omega)$, $\rho_{0}\in
L^{\gamma}(\Omega)$, $(\rho u)_{0}\in L^{1}(\Omega)$ and $B_{0}\in
L^{2}(\Omega)$ to satisfy the conditions (24). Then the system (2)–(17) admits
a weak solution $\\{X,\rho,u,B\\}$ in the sense of Definition 2.1 which in
addition satisfies the energy inequality

$\displaystyle\int_{\Omega}\frac{1}{2}\rho(t)\left|u(t)\right|^{2}+\frac{a}{\gamma-1}\rho^{\gamma}(t)+\frac{1}{2\mu}\left|B(t)\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}2\nu|\mathbb{D}\left(u\right)|^{2}+\lambda\left|\operatorname{div}u\right|^{2}+\frac{1}{\sigma\mu^{2}}\left|\operatorname{curl}B\right|^{2}\
dxd\tau$ $\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\frac{\left|\left(\rho
u\right)_{0}\right|^{2}}{\rho_{0}}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{1}{2\mu}\left|B_{0}\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}\rho g\cdot
u+\frac{1}{\sigma\mu}J\cdot\operatorname{curl}B\ dxd\tau$ (37)

for almost all $t\in[0,T]$.

The remainder of the article is devoted to the proof of Theorem 2.1. In the
following section we begin with an outline of the main ideas by introducing
the approximation method on which the proof is based.

##  3 Approximate system

The biggest challenge in the extension of the proof in the incompressible case
in [1] to the compressible case lies in the correct construction of the
approximate problem. We fix five parameters $\alpha,\epsilon,\eta>0$,
$n\in\mathbb{N}$ and $\Delta t>0$ and introduce an approximation which
consists of five different approximation levels, each of which corresponds to
one of the parameters. The approximate system is chosen such that it is easy
to solve; a solution to the original system is obtained by passing to the
limit in all of the approximation levels. The first three approximation
levels, associated to $\alpha,\epsilon$ and $\eta$, respectively, correspond
to the approximation used in [17] for the purely mechanical problem: On the
$\alpha$\- and $\epsilon$-levels, the system is regularized through the
addition of an artificial pressure term and multiple further regularization
terms. The $\eta$-level consists of a penalization method which allows us to
pass from a fluid-only system to a system containing both a fluid and rigid
bodies. The fourth level, indexed by $n$, describes a Galerkin method used for
solving the approximate momentum equation. Finally, on the fifth level,
associated to $\Delta t$, the induction equation is discretized with respect
to the time variable while the mechanical part of the problem is split up into
a series of time-dependent problems on the small intervals between the
discrete times. In the following we present the complete approximate system on
the highest approximation level and subsequently give a more explicit
description of each included approximation level and its purpose.

Let $\beta>0$ be sufficiently large such that it satisfies in particular
$\beta>\max\\{4,\gamma\\}$ as in [17, Section 6]. Let $V_{n}$,
$n\in\mathbb{N}$, denote the Galerkin space spanned by the first $n$
eigenfunctions of the Lamé equation in $\Omega$ which constitute an
orthonormal basis of $L^{2}(\Omega)$ and an orthogonal basis of
$H_{0}^{1,2}(\Omega)$, c.f. [37, Lemma 4.33]. Then, provided that the
approximate system has already been solved up to the (discrete) time
$(k-1)\Delta t$ for some $k=1,...,\frac{T}{\Delta t}$, the approximate problem
on the interval $[(k-1)\Delta t,k\Delta t]$ consists of finding a solution

$\displaystyle\rho_{\Delta t,k}$ $\displaystyle\in\left\\{\psi\in
C\left([(k-1)\Delta t,k\Delta
t];C^{2,\frac{1}{2}}\left(\overline{\Omega}\right)\right)\bigcap
C^{1}\left([(k-1)\Delta t,k\Delta
t];C^{0,\frac{1}{2}}\left(\overline{\Omega}\right)\right):\
\left.\nabla\psi\cdot\text{n}\right|_{\partial\Omega}=0\right\\},$ (38)
$\displaystyle u_{\Delta t,k}$ $\displaystyle\in C\left([(k-1)\Delta t,k\Delta
t];V_{n}\right),$ (39) $\displaystyle B_{\Delta t}^{k}$ $\displaystyle\in
Y^{k}\left(S_{\Delta t}\right):=\left\\{b\in H^{2,2}\left(\Omega\right):\
b\cdot\text{n}|_{\partial\Omega}=0,\ \text{div}b=0,\ \text{curl}b=0\ \text{in
}S_{\Delta t}(k\Delta t)\bigcap\Omega\right\\}$ (40)

to the system

$\displaystyle\partial_{t}\rho_{\Delta
t,k}+\operatorname{div}\left(\rho_{\Delta t,k}u_{\Delta t,k}\right)=$
$\displaystyle\epsilon\Delta\rho_{\Delta t,k}\quad\text{in }[(k-1)\Delta
t,k\Delta t]\times\Omega,$ (41) $\displaystyle\int_{(k-1)\Delta t}^{k\Delta
t}\int_{\Omega}\partial_{t}\left(\rho_{\Delta t,k}u_{\Delta
t,k}\right)\cdot\phi\ dxdt=$ $\displaystyle\int_{(k-1)\Delta t}^{k\Delta
t}\int_{\Omega}\left(\rho_{\Delta t,k}u_{\Delta t,k}\otimes u_{\Delta
t,k}\right):\mathbb{D}(\phi)+\left(a\rho_{\Delta
t,k}^{\gamma}+\alpha\rho_{\Delta t,k}^{\beta}\right)\operatorname{div}\phi$
$\displaystyle-2\nu\left(\chi_{\Delta t}^{k-1}\right)\mathbb{D}\left(u_{\Delta
t,k}\right):\mathbb{D}(\phi)-\lambda\left(\chi_{\Delta
t}^{k-1}\right)\operatorname{div}\left(u_{\Delta
t,k}\right)\operatorname{div}\phi$ $\displaystyle+\rho_{\Delta
t,k}g\cdot\phi+\frac{1}{\mu}\left(\text{curl}B_{\Delta t}^{k-1}\times
B_{\Delta t}^{k-1}\right)\cdot\phi-\epsilon\left|u_{\Delta
t,k}\right|^{2}u_{\Delta t,k}\cdot\phi$ $\displaystyle-\epsilon\left(\nabla
u_{\Delta t,k}\nabla\rho_{\Delta t,k}\right)\cdot\phi\ dxdt,$ (42)
$\displaystyle-\int_{\Omega}\frac{B_{\Delta t}^{k}-B_{\Delta t}^{k-1}}{\Delta
t}\cdot b\ dx=$
$\displaystyle\int_{\Omega}\left[\frac{1}{\sigma\mu}\operatorname{curl}B_{\Delta
t}^{k}-\tilde{u}_{\Delta t}^{k-1}\times B_{\Delta t}^{k-1}\right.$
$\displaystyle+\left.\frac{\epsilon}{\mu^{2}}\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{2}\operatorname{curl}B_{\Delta t}^{k}-\frac{1}{\sigma}J_{\Delta
t}^{k}\right]\cdot\operatorname{curl}b$
$\displaystyle+\epsilon\operatorname{curl}\left(\operatorname{curl}B_{\Delta
t}^{k}\right)\cdot\operatorname{curl}\left(\operatorname{curl}b\right)\ dx$
(43)

for all $\phi\in C([(k-1)\Delta t,k\Delta t];V_{n})$ and

$\displaystyle b\in W^{k}\left(S_{\Delta t}\right):=$
$\displaystyle\bigg{\\{}b\in H^{2,2}\left(\Omega\right):\
b\cdot\text{n}|_{\partial\Omega}=0,\ \text{curl}\ b=0\ \text{in }S_{\Delta
t}(k\Delta t)\bigcap\Omega\bigg{\\}},$ (44)

which in addition satisfies the initial conditions

$\displaystyle\rho_{\Delta t,k}((k-1)\Delta t;x)$ $\displaystyle=\rho_{\Delta
t,k-1}((k-1)\Delta t;x),\quad\rho_{\Delta t,1}(0;x)=\rho_{0}(x),\quad
x\in\Omega,$ (45) $\displaystyle u_{\Delta t,k}((k-1)\Delta t;x)$
$\displaystyle=u_{\Delta t,k-1}((k-1)\Delta t;x),\quad u_{\Delta
t,1}(0;x)=u_{0}(x),\quad x\in\Omega,$ (46) $\displaystyle B_{\Delta t}^{0}(x)$
$\displaystyle=B_{0}(x),\quad x\in\Omega.$ (47)

Before we proceed with the explanation of the different approximation levels
in (38)–(47), let us clarify the notation introduced in this system: For the
definition of the set $S_{\Delta t}(k\Delta t)$ in (40) and (44), we first
denote by

$\displaystyle O^{i}:=(S^{i}_{0})_{\delta}:=\left\\{x\in S^{i}_{0}:\
\operatorname{dist}\left(x,\partial S_{0}^{i}\right)>\delta\right\\}$

the $\delta$-kernel of the initial domain $S^{i}_{0}$ of the $i$-th body,
where $\delta>0$ is chosen sufficiently small such that for all $i=1,...,m$
the $\delta$-neighbourhood $(O^{i})^{\delta}=\\{x\in\mathbb{R}^{3}:\
\operatorname{dist}(x,(O^{i})^{\delta})<\delta\\}$ of $O^{i}$ coincides with
$S^{i}_{0}$. Such $\delta>0$ exists due to the $C^{2}$-regularity of
$S^{i}_{0}$, c.f. [40, Proposition 2.1]. Then we set

$\displaystyle O:=\bigcup_{i=1}^{m}O^{i}.$

Moreover, we denote by $X_{\Delta t,k}$ the unique solution to the initial
value problem

$\displaystyle\frac{d}{dt}X_{\Delta t,k}(t;x)=$ $\displaystyle
R_{\delta}\left[u_{\Delta t,k}\right]\left(t,X_{\Delta t,k}(t;x)\right),\quad
t\in[(k-1)\Delta t,k\Delta t],$ (48) $\displaystyle X_{\Delta t,k}((k-1)\Delta
t;x)$ $\displaystyle=X_{\Delta t,k-1}((k-1)\Delta t;x),\quad X_{\Delta
t,1}(0;x)=x,\quad x\in\mathbb{R}^{3},$ (49)

where $R_{\delta}[u_{\Delta t,k}](t,\cdot):=u_{\Delta
t,k}(t,\cdot)*\Theta_{\delta}(\cdot)$ and $\Theta_{\delta}$ denotes a radially
symmetric and non-increasing mollifier with respect to the spatial variable.
With this notation at hand we define the domain $S^{i}_{\Delta t}(t)$ of the
$i$-th approximate solid at an arbitrary time $t\in[(k-1)\Delta t,k\Delta
t]\subset[0,T]$ by

$\displaystyle O_{\Delta t}^{i}(t):=X_{\Delta t,k}\left(t;O^{i}\right),\quad
S^{i}_{\Delta t}(t):=\left(O_{\Delta
t}^{i}(t)\right)^{\delta}=\left\\{x\in\mathbb{R}^{3}:\
\operatorname{dist}\left(x,O^{i}_{\Delta t}(t)\right)<\delta\right\\}.$ (50)

Consequently, the approximate solid region at time $t\in[0,T]$ is given by

$\displaystyle S_{\Delta t}(t):=\bigcup_{i=1}^{m}S^{i}_{\Delta t}(t),$

which in particular defines the set $S_{\Delta t}(k\Delta t)$ in (40) and
(44). We note that, by construction, $S_{\Delta t}(t)$ can be an arbitrary
subset of $\mathbb{R}^{3}$, while the corresponding approximate solid time-
space domain $Q^{s}(S_{\Delta t})$, defined according to (1), is always a
subset of the bounded domain $Q$. For later use we further remark that
$S_{\Delta t}^{i}(t)$, as the $\delta$-neighbourhood of a bounded set,
satisfies the cone condition and thus has the property

$\displaystyle\left|\partial S^{i}_{\Delta t}(t)\right|=0\quad\text{for all
}t\in[0,T],\ i=1,...,m.$ (51)

Next, for the definition of the variable viscosity coefficients
$\nu(\chi_{\Delta t}^{k-1})$ and $\lambda(\chi_{\Delta t}^{k-1})$ in the
momentum equation (42) we denote the signed distance function of arbitrary
sets $U\subset\mathbb{R}^{3}$ by

$\displaystyle\textbf{db}_{U}(x):=\operatorname{dist}\left(x,\overline{\mathbb{R}^{3}\setminus
U}\right)-\operatorname{dist}\left(x,\overline{U}\right).$

Further we introduce the signed distance function of the approximate solid
area,

$\displaystyle\chi_{\Delta t}(t,x):=\textbf{db}_{S_{\Delta
t}(t)}(x),\quad\quad\chi_{\Delta t}^{k-1}(t):=\chi_{\Delta t}((k-1)\Delta
t,\cdot)\quad\text{for }t\in[(k-1)\Delta t,k\Delta t].$

Choosing a convex function $H\in C^{\infty}(\mathbb{R})$ such that

$\displaystyle H(z)=0\quad\text{for }z\in(-\infty,0],\quad\quad
H(z)>0\quad\text{for }z\in(0,+\infty),$ (52)

we then define the variable viscosity coefficients by

$\displaystyle\nu\left(\chi_{\Delta
t}^{k-1}\right):=\nu+\frac{1}{\eta}H\left(\chi_{\Delta
t}^{k-1}\right),\quad\quad\lambda\left(\chi_{\Delta
t}^{k-1}\right):=\lambda+\frac{1}{\eta}H\left(\chi_{\Delta t}^{k-1}\right).$
(53)

Finally, in the induction equation (43) the function $\tilde{u}_{\Delta
t}^{k-1}$ is defined by

$\displaystyle\tilde{u}_{\Delta
t}^{k-1}(x):=\left\\{\begin{array}[]{ll}\frac{1}{\Delta
t}\displaystyle\int_{(k-2)\Delta t}^{(k-1)\Delta t}u_{\Delta
t,k-1}(\tau,x)d\tau&\text{if }k\geq 2,\\\\[14.22636pt] u_{\Delta
t}^{0}(x)&\text{if }k=1.\end{array}\right.$ (56)

while the discretized external force $J_{\Delta t}^{k}$ is defined by

$\displaystyle J_{\Delta t}^{k}:=J_{\omega}(k\Delta t),\quad
J_{\omega}(t):=\int_{0}^{T}\theta_{\omega}\left(t+\omega\frac{T-2t}{T}-s\right)J(s)\
ds,$ (57)

for another mollifier $\theta_{\omega}:\mathbb{R}\to\mathbb{R}$ and a suitable
choice of $\omega=\omega(\Delta t),$ $\omega(\Delta t)\rightarrow 0$ for
$\Delta t\rightarrow 0$. We are now in the position to discuss the several
approximation levels and the reasons why they are required. We start from the
highest level.
The $\Delta t$-level constitutes the level which contains most of the
difficulties. It is here where the main novelties of our proof enter, compared
to the incompressible setting in [1]. The situation presents itself in the
following way: On the one hand, the dependence of the test functions (21),
(22) for the induction equation and the momentum equation on the solution of
the system hinders the effort to solve all of the equations in the system
simultaneously. While we can deal with the test functions in the momentum
equation by means of a penalization method (c.f. the $\eta$-level below), the
same does not work in case of the induction equation. This suggests to
decouple the system by the use of a classical time discretization, e.g. via
the Rothe method ([38, Section 8.2]). In this way, at each fixed discrete time
we can first determine a velocity field and, from this, the position of the
approximate solid. This in turn determines the test functions (44) and solving
the discretized induction equation (43) becomes a routine matter. On the other
hand, however, the various functions evaluated at different discrete times in
a fully discretized system complicate the derivation of a meaningful energy
inequality. The author could not find a way to transfer several of the
techniques known for the continuous compressible Navier-Stokes system (c.f.
[37, Sections 7.6.5, 7.6.6, 7.7.4.2]) - in particular, the proof of the non-
negativity of the density - to the discrete case and it did not seem to be
possible to derive the uniform bounds required for the limit passage with
respect to $\Delta t\rightarrow 0$.
Our solution to this dilemma consists of considering, instead of a strictly
discretized system, a hybrid system in which the induction equation (43) is
indeed discretized by the Rothe method, while the continuity equation and the
momentum equation are solved as continuous equations on the small intervals
between each pair of consecutive discrete times, c.f (41) and (42). Through
this, the solution dependence of the test functions in the induction equation
can be handled as in the fully discrete system, while the mechanical part of
the energy inequality - with the density bounded away from zero - can be
derived as in the strictly continuous case. Moreover, under the consideration
of piecewise linear interpolants of the discrete functions, the discrete
induction equation (43) also leads to a continuous energy estimate, which can
be combined with the mechanical estimate to obtain the full energy inequality,
c.f. Section 5.1. A hybrid approximation scheme was already used in our proof
in the incompressible case [1]. In that case, however, the major part of the
system could be discretized in time while only the transport equation for the
characteristic function of the solid domain had to be treated as a continuous
problem on small time intervals. The idea for the latter procedure, in turn,
stems from [26].
The Galerkin method carried out on the $n$-level is used to solve the
continuous momentum equation (41) on the small time intervals from the $\Delta
t$-level by a standard procedure. The Galerkin-regularity of the velocity
field furthermore helps us during the limit passage with respect to $\Delta
t\rightarrow 0$, c.f. (91) below.
After letting $n$ tend to $\infty$ we find ourselves in the same situation as
in the approximation of the exclusively mechanical system in [17, Section 6].
Indeed, the remaining three approximation levels correspond directly to the
three level approximation scheme used in that article. Hence, for the
mechanical part of our problem we can follow exactly the strategy used
therein. Moreover, the limit passages in the induction equation from here on
do not contain any new difficulties anymore. Consequently, after the limit
passage in $n$ the rest of the proof will become a routine matter.
The penalization method on the $\eta$-level - c.f. Section 7 \- is the same as
the one used for the fluid-rigid bodies system in [17] and was, before that,
also used for example for the corresponding two dimensional problem in [40].
The idea behind it is to approximate the entirety of the fluid and the rigid
bodies by a fluid in the whole domain with viscosity tending to infinity in
the later solid regions. Mathematically this is implemented through the
variable viscosity coefficients (53). Due to the choice of the function $H$ in
(52) these coefficients blow up in the approximate solid region once we let
$\eta$ tend to $0$ and, thanks to the energy inequality, this will cause the
limit velocity field $u$ to coincide with a rigid velocity field in each body.
Moreover, the positions $S^{i}(t)$ of the bodies in the $\eta$-limit are
determined through the flow curves of $R_{\delta}[u]$, c.f. (48) and (50).
This regularized velocity field has the useful property that, for any domain
$U\subset\mathbb{R}^{3}$, it holds

$\displaystyle\mathbb{D}(u(t,\cdot))=0\quad\text{in
}U\quad\quad\Rightarrow\quad\quad
R_{\delta}[u](t,\cdot)=u(t,\cdot)\quad\text{in }U_{\delta}=\left\\{x\in U:\
\operatorname{dist}\left(x,\partial U\right)>\delta\right\\},$ (58)

c.f. [17, Remark 6.1]. Hence $R_{\delta}[u]$ coincides with $u$ itself in the
sets $O^{i}(t)\subset S^{i}(t)$, in which $\mathbb{D}(u)=0$. Consequently, the
rigid velocity fields coinciding with $u$ in $S^{i}(t)$ also coincide, in
$O^{i}(t)$, with the velocity field $R_{\delta}[u]$ which determines the
motion of the bodies. In particular, this shows that the bodies $S^{i}(t)$ are
indeed rigid.
On the $\epsilon$-level, the continuity equation (41) is regularized through
the additional Laplacian $\epsilon\Delta\rho_{\Delta t,k}$. The additional
quantity $\epsilon(\nabla u_{\Delta t,k}\nabla\rho_{\Delta t,k})$ in the
momentum equation (42) ensures that the energy inequality is preserved under
this regularization. This procedure, which is classical in the theory of the
compressible Navier-Stokes equations, is what guarantees us the non-negativity
of the density, c.f. [37, Section 7.3.8]. The other regularization term
$\epsilon|u_{\Delta t,k}|^{2}u_{\Delta t,k}$ in (42) is needed in the time
discrete level where, as opposed to the continuous case, the mixed terms from
the momentum equation and the induction equation do not annihilate each other
in the energy inequality, which prevents a direct application of the Gronwall
lemma. The quantity $\epsilon|u_{\Delta t,k}|^{2}u_{\Delta t,k}$ can be used
to control the velocity part of these mixed terms. The $4$-double-curl
$\epsilon\operatorname{curl}(\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{2}\operatorname{curl}B_{\Delta t}^{k})$ in the induction
equation (43) fulfills, as in the incompressible setting in [1], the same
purpose for the magnetic part of the mixed terms so that we are able to derive
uniform bounds from the energy inequality nevertheless, c.f. Section 5.1. We
remark that this control of the mixed terms was also the motivation for the
definition of the velocity field (56) in the discrete induction equation:
Indeed, defining this quantity as a mean value of the velocity field obtained
from the momentum equation on the intervals $[(k-2)\Delta t,(k-1)\Delta t]$,
we can absorb it into the left-hand side of the energy inequality thanks to
the above-mentioned regularization terms. If instead the term was defined,
more intuitively, as a pointwise evaluation of $u_{\Delta t,k}$, we would not
be able to handle it. The last regularization term in (43), the $4$-th curl of
$B_{\Delta t}^{k}$, enables us to express the induction equation via some
weakly continuous operator on $Y^{k}(S_{\Delta t})$. Seeing that this operator
is moreover coercive, we will then be able to infer the existence of
$B_{\Delta t}^{k}$, c.f. Section 4.2.
Finally, on the $\alpha$-level, the artificial pressure term
$\alpha\rho_{\Delta t,k}^{\beta}$ is added to the momentum equation (42).
Again this method is already well-known from the general existence theory for
the compressible Navier-Stokes system, c.f. [37, Section 7.3.8]. The
artificial pressure gives us an additional amount of integrability of the
density and its gradient, required to pass to the limit in the term
$\epsilon(\nabla u_{\Delta t,k}\nabla\rho_{\Delta t,k})$ from the
$\epsilon$-level, c.f. [37, Section 7.8.2]. It furthermore simplifies the
limit passage with respect to $\epsilon\rightarrow 0$, since the additional
integrability allows for the use of the regularization technique by DiPerna
and Lions, c.f. [37, Lemma 6.8, Lemma 6.9].

##  4 Existence of the approximate solution

We begin the proof of Theorem 2.1 by showing the existence of a solution to
the approximate problem (38)–(47) on the highest approximation level.

###  4.1 Existence of the density and velocity

The existence of the density and the velocity field on the Galerkin level can
be shown by classical methods, c.f. for example [37, Section 7.7]. More
precisely, the continuity equation (41) and the momentum equation (42) can be
solved simultaneously by means of a fixed point argument: For fixed $w\in
C([(k-1)\Delta t,k\Delta t];V_{n})$ we consider the Neumann problem

$\displaystyle\ \quad\quad\quad\partial_{t}\rho+\operatorname{div}\left(\rho
w\right)=\epsilon\Delta\rho\quad\quad\text{in }[(k-1)\Delta t,k\Delta
t]\times\Omega,$ (59)
$\displaystyle\left.\nabla\rho\cdot\text{n}\right|_{\partial\Omega}=0,\quad\rho((k-1)\Delta
t,\cdot)=\rho_{\Delta t,k-1}((k-1)\Delta t;\cdot)\quad\quad\text{in }\Omega,$
(60) $\displaystyle\ \quad\quad\quad 0<\underline{\rho}\leq\rho_{\Delta
t,k-1}((k-1)\Delta t;x)\leq\overline{\rho}<\infty,\quad\quad\text{in }\Omega.$
(61)

It is well known (c.f. Lemma 10.1 in the appendix) that (59)–(61) admits a
unique solution

$\displaystyle\rho=\rho(w)\in C\left([(k-1)\Delta t,k\Delta
t];C^{2,\frac{1}{2}}\left(\overline{\Omega}\right)\right)\bigcap
C^{1}\left([(k-1)\Delta t,k\Delta
t];C^{0,\frac{1}{2}}\left(\overline{\Omega}\right)\right),$

which satisfies the estimate

$\displaystyle 0<\underline{\rho}\exp\left(-\int_{(k-1)\Delta
t}^{t}\left\|w(\tau)\right\|_{V_{n}}\
d\tau\right)\leq\rho(w)(t,x)\leq\overline{\rho}\exp\left(\int_{(k-1)\Delta
t}^{t}\left\|w(\tau)\right\|_{V_{n}}\ d\tau\right)<\infty$ (62)

for all $(t,x)\in[0,T]\times\overline{\Omega}$. Further, we consider a
linearized version of the momentum equation (42). Given $w\in C([(k-1)\Delta
t,k\Delta t];V_{n})$ and the associated solution $\rho(w)$ to the Neumann
problem (59)–(61), we seek $u\in C([(k-1)\Delta t,k\Delta t];V_{n})$ such that

$\displaystyle\int_{\Omega}\partial_{t}\left(\rho(w)u\right)\cdot\phi\ dx=$
$\displaystyle\int_{\Omega}\left(\rho(w)w\otimes
u\right):\mathbb{D}(\phi)+\left(a\rho^{\gamma}(w)+\alpha\rho^{\beta}(w)\right)\operatorname{div}\phi$
$\displaystyle-2\nu\left(\chi_{\Delta
t}^{k-1}\right)\mathbb{D}\left(u\right):\mathbb{D}(\phi)-\lambda\left(\chi_{\Delta
t}^{k-1}\right)\operatorname{div}\left(u\right)\operatorname{div}\phi$
$\displaystyle+\rho(w)g\cdot\phi+\frac{1}{\mu}\left(\text{curl}B_{\Delta
t}^{k-1}\times B_{\Delta t}^{k-1}\right)\cdot\phi$
$\displaystyle-\epsilon\left(\nabla
u\nabla\rho(w)\right)\cdot\phi-\epsilon\left|w\right|^{2}u\cdot\phi\
dx\quad\quad\text{in }[(k-1)\Delta t,k\Delta t],$ $\displaystyle u((k-1)\Delta
t,\cdot)=$ $\displaystyle u_{\Delta t,k-1}((k-1)\Delta
t;\cdot)\quad\quad\quad\quad\quad\quad\quad\quad\ \text{in }\Omega,$

for all $\phi\in C([(k-1)\Delta t,k\Delta t];V_{n})$. Under exploitation of
the fact that, by (62), $\rho(w)$ is bounded away from $0$ and the linearity
of the problem, it follows from classical methods that this problem admits a
unique solution $u=u(w)\in C([(k-1)\Delta t,k\Delta t];V_{n})$. We can thus
define an operator

$\displaystyle\mathbb{T}:C\left([(k-1)\Delta t,k\Delta
t];V_{n}\right)\rightarrow C\left([(k-1)\Delta t,k\Delta
t];V_{n}\right),\quad\mathbb{T}(w):=u.$

It is easy to see that $\mathbb{T}$ is continuous and compact and, by an
energy estimate, fixed points of $s\mathbb{T}$ for $s\in[0,1]$ are bounded in
$C([(k-1)\Delta t,k\Delta t];V_{n})$, uniformly with respect to $s$. Under
these conditions the Schaefer fixed point theorem (see [14, Section 9.2.2,
Theorem 4]) tells us that $\mathbb{T}$ possesses a fixed point $u_{\Delta
t,k}\in C([(k-1)\Delta t,k\Delta t];V_{n})$, which constitutes the desired
solution to the initial value problem (42), (46). Furthermore, by
construction, the associated density $\rho_{\Delta t,k}:=\rho(u_{\Delta t,k})$
is the desired solution to the corresponding initial value problem (41), (45)
for the density.

###  4.2 Existence of the magnetic induction

The existence of the magnetic induction is obtained as in the incompressible
case, c.f. [1, Section 3]. We equip the space $Y^{k}(S_{\Delta t})$ with the
norm $\|\cdot\|_{H^{2,2}(\Omega)}$ and express the identity (43) through the
equation

$\left\langle AB_{\Delta t}^{k},b\right\rangle_{\left(Y^{k}(S_{\Delta
t})\right)^{*}\times Y^{k}(S_{\Delta t})}=\left\langle
f,b\right\rangle_{\left(Y^{k}(S_{\Delta t})\right)^{*}\times Y^{k}(S_{\Delta
t})}\quad\forall b\in Y^{k}\left(S_{\Delta t}\right),$ (63)

where the operator $A:Y^{k}(S_{\Delta t})\rightarrow(Y^{k}(S_{\Delta t}))^{*}$
and the right-hand side $f\in(Y^{k}(S_{\Delta t}))^{*}$ are defined by

$\displaystyle\left\langle A(B),b\right\rangle_{\left(Y^{k}(S_{\Delta
t})\right)^{*}\times Y^{k}(S_{\Delta t})}:=$
$\displaystyle\int_{\Omega}\frac{B}{\Delta t}\cdot
b+\epsilon\text{curl}\left(\text{curl}B\right)\cdot\text{curl}\left(\text{curl}b\right)$
$\displaystyle+\left[\frac{1}{\sigma\mu}\text{curl}B+\frac{\epsilon}{\mu^{2}}\left|\text{curl}B\right|^{2}\text{curl}B\right]\cdot\text{curl}b\
dx,$ $\displaystyle\left\langle f,b\right\rangle_{\left(Y^{k}(S_{\Delta
t})\right)^{*}\times Y^{k}(S_{\Delta t})}:=$
$\displaystyle\int_{\Omega}\frac{B_{\Delta t}^{k-1}}{\Delta t}\cdot
b+\left(\tilde{u}_{\Delta t}^{k-1}\times B_{\Delta
t}^{k-1}\right)\cdot\text{curl}b+\frac{1}{\sigma}J_{\Delta
t}^{k}\cdot\operatorname{curl}b\ dx$

for any $B,b\in Y^{k}(S_{\Delta t})$. The operator $A$ is weakly continuous
and coercive and consequently surjective from $Y^{k}(S_{\Delta t})$ onto
$(Y^{k}(S_{\Delta t}))^{*}$, c.f. [22, Theorem 1.2]. In particular, there
exists a function $B_{\Delta t}^{k}\in Y^{k}(S_{\Delta t})$ which satisfies
(63) and thus the induction equation (43) for all $b\in Y^{k}(S_{\Delta t})$.
Finally, by the Helmholtz decomposition, see [43, Theorem 4.2], we infer that
(43) does not only hold for $b\in Y^{k}(S_{\Delta t})$ but also for the (not
divergence-free) test functions $b\in W^{k}(S_{\Delta t})$. Altogether, we
have shown the following result:

###### Proposition 4.1

Let $n\in\mathbb{N}$, $\Delta t,\eta,\epsilon,\alpha>0$ such that
$\frac{T}{\Delta t}\in\mathbb{N}$ and let $\beta>\max\\{4,\gamma\\}$ be
sufficiently large. Assume the conditions of Theorem 2.1 to be satisfied.
Moreover, assume that

$\displaystyle\rho_{0}\in$ $\displaystyle
C^{2,\frac{1}{2}}\left(\overline{\Omega}\right),\quad\quad(\rho u)_{0}\in
C^{2}\left(\overline{\Omega}\right),\quad\quad u_{0}:=P_{n}\left(\frac{(\rho
u)_{0}}{\rho_{0}}\right)\in V_{n},\quad\quad B_{0}\in H^{2,2}(\Omega),$
$\displaystyle
0<\alpha\leq\rho_{0}\leq\alpha^{-\frac{1}{2\beta}},\quad\quad\left.\nabla\rho_{0}\cdot\operatorname{n}\right|_{\partial\Omega}=0,\quad\quad\operatorname{div}B_{0}=0,\quad\quad\left.B_{0}\cdot\operatorname{n}\right|_{\partial\Omega}=0,$

where $P_{n}$ denotes the orthogonal projection of $L^{2}(\Omega)$ onto
$V_{n}$. Finally, for all $k=1,...,\frac{T}{\Delta t}$, let $J_{\Delta t}^{k}$
be defined by (57). Then, for each $k=1,...,\frac{T}{\Delta t}$, there exist
functions

$\displaystyle 0\leq\rho_{\Delta t,k}$ $\displaystyle\in\bigg{\\{}\psi\in
C\left([(k-1)\Delta t,k\Delta
t];C^{2,\frac{1}{2}}\left(\overline{\Omega}\right)\right)\bigcap
C^{1}\left([(k-1)\Delta t,k\Delta
t];C^{0,\frac{1}{2}}\left(\overline{\Omega}\right)\right):\
\left.\nabla\psi\cdot\operatorname{n}\right|_{\partial\Omega}=0\bigg{\\}},$
$\displaystyle\quad\quad\quad u_{\Delta t,k}\in C\left([(k-1)\Delta t,k\Delta
t];V_{n}\right),\quad\quad B_{\Delta t}^{k}\in Y^{k}\left(S_{\Delta
t}\right),$

which satisfy the continuity equation (41), the momentum equation (42) for all
test functions $\phi\in C\left([(k-1)\Delta t,k\Delta t];V_{n}\right)$ and the
induction equation (43) for all test functions $b\in W^{k}(S_{\Delta t})$ as
well as the initial conditions (45)–(47).

##  5 Limit passage in the time discretization

We continue by passing to the limit with respect to $\Delta t\rightarrow 0$.
As in [1], we first need to assemble the functions constructed in Section 4,
defined up to now only on small time intervals or in discrete time points, to
functions defined on the whole time interval $[0,T]$. More precisely, for
functions $f_{\Delta t,k}$, defined on $[(k-1)\Delta t,k\Delta t]\times\Omega$
for $k=1,...,\frac{T}{\Delta t}$, we introduce the assembled functions

$\displaystyle f_{\Delta t}(t,\cdot):=f_{\Delta t,k}(t,\cdot)\quad\quad\forall
t\in[(k-1)\Delta t,k\Delta t),\quad k=1,...,\frac{T}{\Delta t}$ (64)

while for discrete functions $h_{\Delta t}^{k}$, defined on $\Omega$ for
$k=0,...,\frac{T}{\Delta t}$, we introduce the piecewise affine and piecewise
constant interpolants

$\displaystyle h_{\Delta t}(t)$ $\displaystyle:=\left(\frac{t}{\Delta
t}-(k-1)\right)h_{\Delta t}^{k}+\left(k-\frac{t}{\Delta t}\right)h_{\Delta
t}^{k-1}\ \ \ $ $\displaystyle\forall$ $\displaystyle t\in[(k-1)\Delta
t,k\Delta t),\quad k=1,...,\frac{T}{\Delta t},$ (65)
$\displaystyle\overline{h}_{\Delta t}(t)$ $\displaystyle:=h_{\Delta t}^{k}\ \
\ $ $\displaystyle\forall$ $\displaystyle t\in[(k-1)\Delta t,k\Delta t),\quad
k=0,...,\frac{T}{\Delta t},$ (66) $\displaystyle\overline{h}^{\prime}_{\Delta
t}(t)$ $\displaystyle:=h_{\Delta t}^{k-1}\ \ \ $ $\displaystyle\forall$
$\displaystyle t\in[(k-1)\Delta t,k\Delta t),\quad k=1,...,\frac{T}{\Delta
t}.$ (67)

Moreover, in order to derive a suitable energy inequality in Section 5.1
below, we also introduce a piecewise affine interpolation of the square of the
$L^{2}(\Omega)$-norm,

$\displaystyle h_{\Delta t,\|\cdot\|}(t):=\left(\frac{t}{\Delta
t}-(k-1)\right)\left\|h_{\Delta
t}^{k}\right\|_{L^{2}(\Omega)}^{2}+\left(k-\frac{t}{\Delta
t}\right)\left\|h_{\Delta
t}^{k-1}\right\|_{L^{2}(\Omega)}^{2},\quad\quad\forall t\in[(k-1)\Delta
t,k\Delta t)$

for any $k=1,...,\frac{T}{\Delta t}.$ Since, by Proposition 4.1, the functions
$\rho_{\Delta t,k}$ and $u_{\Delta t,k}$ satisfy the continuity equation (41),
the momentum equation (42), and the initial conditions (45)–(47), it follows
from the definition of $\rho_{\Delta t}$ and $u_{\Delta t}$ in (64) as well as
of $B_{\Delta t}$, $\overline{B}_{\Delta t}$ and
$\overline{B}^{\prime}_{\Delta t}$ in (65)–(67) that these functions solve the
continuity equation

$\displaystyle\partial_{t}\rho_{\Delta t}+\operatorname{div}\left(\rho_{\Delta
t}u_{\Delta t}\right)=\epsilon\rho_{\Delta t}\quad\text{a.e. in
}[0,T]\times\Omega,$ (68)

the momentum equation

$\displaystyle\int_{0}^{T}\int_{\Omega}\partial_{t}\left(\rho_{\Delta
t}u_{\Delta t}\right)\cdot\phi\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left(\rho_{\Delta t}u_{\Delta
t}\otimes u_{\Delta t}\right):\mathbb{D}(\phi)+\left(a\rho_{\Delta
t}^{\gamma}+\alpha\rho_{\Delta t}^{\beta}\right)\operatorname{div}\phi$
$\displaystyle-2\nu\left(\overline{\chi}_{\Delta
t}^{\prime}\right)\mathbb{D}(u_{\Delta
t}):\mathbb{D}(\phi)-\lambda\left(\overline{\chi}_{\Delta
t}^{\prime}\right)\operatorname{div}u_{\Delta
t}\operatorname{div}\phi-\epsilon\left|u_{\Delta t}\right|^{2}u_{\Delta
t}\cdot\phi$ $\displaystyle+\rho_{\Delta
t}g\cdot\phi+\frac{1}{\mu}\left(\operatorname{curl}\overline{B}_{\Delta
t}^{\prime}\times\overline{B}_{\Delta
t}^{\prime}\right)\cdot\phi-\epsilon\left(\nabla u_{\Delta
t}\nabla\rho_{\Delta t}\right)\cdot\phi\ dxdt,$ (69)

for any $\phi\in C([0,T];V_{n})$ and the initial conditions

$\displaystyle\rho_{\Delta t}(0)=\rho_{0},\quad u_{\Delta t}(0)=u_{0},\quad
B_{\Delta t}(0)=B_{0}.$ (70)

Furthermore, from $X_{\Delta t,k}$ being the unique solution to the initial
value problem (48), (49), it follows that $X_{\Delta t}$ is the unique
solution to

$\displaystyle\frac{dX_{\Delta t}(t;x)}{dt}=R_{\delta}\left[u_{\Delta
t}\right]\left(t,X_{\Delta t}(t;x)\right),\quad\text{for }t\in[0,T],\quad\quad
X_{\Delta t}(0;x)=x,\quad\text{for }x\in\mathbb{R}^{3}.$ (71)

Finally, we consider functions

$\displaystyle b\in L^{4}\left(0,T;H_{0}^{2,2}(\Omega)\right)\quad\text{such
that}\quad b(\tau)\in W^{l}\left(S_{\Delta t}\right)\quad\text{for a.a.
}\tau\in[(l-1)\Delta t,l\Delta t],\ l=1,...,\frac{T}{\Delta t}$ (72)

and realize that, after a density argument, the discrete induction equation
(43) at time $k\Delta t$ can be tested by $b(t)$ for almost all
$t\in[(k-1)\Delta t,k\Delta t]$. After integration over $[(k-1)\Delta
t,k\Delta t]$ and summation over $k$ this yields

$\displaystyle\int_{0}^{T}\int_{\Omega}\partial_{t}B_{\Delta t}\cdot b\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}\overline{B}_{\Delta
t}+\overline{\tilde{u}}^{\prime}_{\Delta t}\times\overline{B}_{\Delta
t}^{\prime}+\frac{1}{\sigma}\overline{J}_{\Delta
t}-\frac{\epsilon}{\mu^{2}}\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{2}\operatorname{curl}\overline{B}_{\Delta
t}\right]\cdot\operatorname{curl}b$
$\displaystyle-\epsilon\operatorname{curl}\left(\operatorname{curl}\overline{B}_{\Delta
t}\right)\cdot\operatorname{curl}\left(\operatorname{curl}b\right)\ dxdt.$
(73)

###  5.1 Energy inequality on the $\Delta t$-level

In contrast to the incompressible setting in [1] we have to combine the
discrete induction equation (43) with the continuous Navier-Stokes equations
(41), (42) in a suitable way in order to derive an energy inequality at the
$\Delta t$-level. We pick an arbitrary $t\in(0,T]$ and choose
$k\in\left\\{1,...,\frac{T}{\Delta t}\right\\}$, $\xi\in[0,\Delta t)$ such
that $t=k\Delta t-\xi$. For the magnetic part of the energy inequality we test
the induction equation (43) by $\frac{1}{\mu}B_{\Delta t}^{k}$, which leads to

$\displaystyle\frac{1}{2\mu}\partial_{t}B_{\Delta
t,\|\cdot\|}(t)+\frac{1}{\mu}\int_{\Omega}\frac{1}{\sigma\mu}\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{2}+\frac{\epsilon}{\mu^{2}}\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{4}+\epsilon\left|\Delta B_{\Delta t}^{k}\right|^{2}\ dx$
$\displaystyle=$ $\displaystyle\frac{1}{2\mu\Delta t}\left(\left\|B_{\Delta
t}^{k}\right\|_{L^{2}(\Omega)}^{2}-\left\|B_{\Delta
t}^{k-1}\right\|_{L^{2}(\Omega)}^{2}\right)+\frac{1}{\mu}\int_{\Omega}\frac{1}{\sigma\mu}\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{2}+\frac{\epsilon}{\mu^{2}}\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{4}+\epsilon\left|\Delta B_{\Delta t}^{k}\right|^{2}\ dx$
$\displaystyle\leq$
$\displaystyle\frac{1}{\mu}\int_{\Omega}\left(\tilde{u}_{\Delta t}^{k-1}\times
B_{\Delta t}^{k-1}\right)\cdot\operatorname{curl}B_{\Delta
t}^{k}+\frac{1}{\sigma}J_{\Delta t}^{k}\cdot\operatorname{curl}B_{\Delta
t}^{k}\ dx.$ (74)

Since corresponding estimates hold true also for all time indices
$l=1,...,k-1$ we can integrate (discretely) over the interval $[0,t]$, which
yields

$\displaystyle\frac{1}{2\mu}B_{\Delta t,\|\cdot\|}(t)+\Delta
t\sum_{l=1}^{k-1}\frac{1}{\mu}\int_{\Omega}\frac{1}{\sigma\mu}\left|\operatorname{curl}B_{\Delta
t}^{l}\right|^{2}+\frac{\epsilon}{\mu^{2}}\left|\operatorname{curl}B_{\Delta
t}^{l}\right|^{4}+\epsilon\left|\Delta B_{\Delta t}^{l}\right|^{2}\ dx$
$\displaystyle+\frac{\Delta
t-\xi}{\mu}\int_{\Omega}\frac{1}{\sigma\mu}\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{2}+\frac{\epsilon}{\mu^{2}}\left|\operatorname{curl}B_{\Delta
t}^{k}\right|^{4}+\epsilon\left|\Delta B_{\Delta t}^{k}\right|^{2}\ dx$
$\displaystyle\leq$ $\displaystyle\frac{1}{2\mu}B_{\Delta
t,\|\cdot\|}(0)+\Delta
t\sum_{l=1}^{k-1}\frac{1}{\mu}\int_{\Omega}\left(\tilde{u}_{\Delta
t}^{l-1}\times B_{\Delta t}^{l-1}\right)\cdot\operatorname{curl}B_{\Delta
t}^{l}+\frac{1}{\sigma}J_{\Delta t}^{l}\cdot\operatorname{curl}B_{\Delta
t}^{l}\ dx$ $\displaystyle+\frac{\Delta
t-\xi}{\mu}\int_{\Omega}\left(\tilde{u}_{\Delta t}^{k-1}\times B_{\Delta
t}^{k-1}\right)\cdot\operatorname{curl}B_{\Delta
t}^{k}+\frac{1}{\sigma}J_{\Delta t}^{k}\cdot\operatorname{curl}B_{\Delta
t}^{k}\ dx$ $\displaystyle\leq$
$\displaystyle\frac{1}{2\mu}\int_{\Omega}\left|B_{0}\right|^{2}\ dx+\Delta
t\sum_{l=1}^{k-1}\left[\frac{c}{\epsilon}\left\|B_{\Delta
t}^{l-1}\right\|_{L^{2}(\Omega)}^{2}+\frac{\epsilon}{8}\left\|\tilde{u}_{\Delta
t}^{l-1}\right\|_{L^{4}(\Omega)}^{4}+\frac{\epsilon}{8\mu^{3}}\left\|\operatorname{curl}B_{\Delta
t}^{l}\right\|_{L^{4}(\Omega)}^{4}+\frac{c}{\epsilon^{\frac{1}{3}}}\left\|J_{\Delta
t}^{l}\right\|_{L^{\frac{4}{3}}(\Omega)}^{\frac{4}{3}}\right.$
$\displaystyle\left.+\frac{\epsilon}{8\mu^{3}}\left\|\operatorname{curl}B_{\Delta
t}^{l}\right\|_{L^{4}(\Omega)}^{4}\right]+\left(\Delta
t-\xi\right)\left[\frac{c}{\epsilon}\left\|B_{\Delta
t}^{k-1}\right\|_{L^{2}(\Omega)}^{2}+\frac{\epsilon}{8}\left\|\tilde{u}_{\Delta
t}^{k-1}\right\|_{L^{4}(\Omega)}^{4}+\frac{\epsilon}{8\mu^{3}}\left\|\operatorname{curl}B_{\Delta
t}^{k}\right\|_{L^{4}(\Omega)}^{4}\right.$
$\displaystyle\left.+\frac{c}{\epsilon^{\frac{1}{3}}}\left\|J_{\Delta
t}^{k}\right\|_{L^{\frac{4}{3}}(\Omega)}^{\frac{4}{3}}+\frac{\epsilon}{8\mu^{3}}\left\|\operatorname{curl}B_{\Delta
t}^{k}\right\|_{L^{4}(\Omega)}^{4}\right]$ (75)

under exploitation of Hölder’s and Young’s inequalities in the last estimate.
On the right-hand side of (75) we further estimate, due to the definition of
$\tilde{u}_{\Delta t}^{l-1}$ in (56) and Jensen’s inequality,

$\displaystyle\left\|\tilde{u}_{\Delta
t}^{l-1}\right\|_{L^{4}(\Omega)}^{4}\leq\frac{1}{\Delta
t}\int_{\Omega}\int_{(l-2)\Delta t}^{(l-1)\Delta t}\left|u_{\Delta
t,l-1}(\tau)\right|^{4}\ d\tau dx.$

Moreover, a direct calculation yields

$\displaystyle\Delta t\sum_{l=1}^{k-1}\frac{c}{\epsilon}\left\|B_{\Delta
t}^{l-1}\right\|_{L^{2}(\Omega)}^{2}+\left(\Delta
t-\xi\right)\frac{c}{\epsilon}\left\|B_{\Delta
t}^{k-1}\right\|_{L^{2}(\Omega)}^{2}\leq$ $\displaystyle\frac{c\Delta
t}{\epsilon}\left(\frac{\left\|B_{0}\right\|_{L^{2}(\Omega)}^{2}}{2}+\frac{\left\|B_{\Delta
t}^{k-1}\right\|_{L^{2}(\Omega)}^{2}}{2}+\sum_{l=2}^{k-1}\left\|B_{\Delta
t}^{l-1}\right\|_{L^{2}(\Omega)}^{2}\right)$ $\displaystyle=$
$\displaystyle\frac{c}{\epsilon}\int_{0}^{(k-1)\Delta t}B_{\Delta
t,\|\cdot\|}(\tau)\ d\tau\leq\frac{c}{\epsilon}\int_{0}^{t}B_{\Delta
t,\|\cdot\|}(\tau)\ d\tau.$ (76)

Hence, absorbing the $\|\operatorname{curl}B_{\Delta
t}^{l}\|_{L^{4}(\Omega)}^{4}$-terms in (75) into the left-hand side and
expressing the sums as integrals, we end up with

$\displaystyle\frac{1}{2\mu}B_{\Delta
t,\|\cdot\|}(t)+\frac{1}{\mu}\int_{0}^{t}\int_{\Omega}\frac{1}{\sigma\mu}\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{2}+\frac{6\epsilon}{8\mu^{2}}\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{4}+\epsilon\left|\Delta\overline{B}_{\Delta t}\right|^{2}\ dxd\tau$
$\displaystyle\leq$
$\displaystyle\frac{1}{2\mu}\int_{\Omega}\left|B_{0}\right|^{2}\
dx+\frac{c}{\epsilon^{\frac{1}{3}}}\int_{0}^{T}\int_{\Omega}\left|\overline{J}_{\Delta
t}\right|^{\frac{4}{3}}\ dxd\tau+\frac{\epsilon}{8}\Delta
t\int_{\Omega}\left|u_{0}\right|^{4}\
dx+\frac{\epsilon}{8}\int_{0}^{(k-1)\Delta t}\int_{\Omega}\left|u_{\Delta
t}\right|^{4}\ d\tau dx$
$\displaystyle+\frac{c}{\epsilon}\int_{0}^{t}B_{\Delta t,\|\cdot\|}(\tau)\
d\tau.$ (77)

For the mechanical part of the energy inequality we test the continuity
equation (68) by $\frac{1}{2}|u_{\Delta t}|^{2}$, the momentum equation (69)
by $u_{\Delta t}$, add up the resulting equations and obtain

$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{\Delta t}(t)\left|u_{\Delta
t}(t)\right|^{2}+\frac{a}{\gamma-1}\rho_{\Delta
t}^{\gamma}(t)+\frac{\alpha}{\beta-1}\rho_{\Delta t}^{\beta}(t)\
dx+\int_{0}^{t}\int_{\Omega}2\nu\left(\overline{\chi}_{\Delta
t}^{\prime}\right)|\mathbb{D}\left(u_{\Delta t}\right)|^{2}$
$\displaystyle+\lambda\left(\overline{\chi}_{\Delta
t}^{\prime}\right)\left|\operatorname{div}u_{\Delta
t}\right|^{2}+a\epsilon\gamma\rho_{\Delta
t}^{\gamma-2}\left|\nabla\rho_{\Delta
t}\right|^{2}+\alpha\epsilon\beta\rho_{\Delta
t}^{\beta-2}\left|\nabla\rho_{\Delta t}\right|^{2}+\epsilon\left|u_{\Delta
t}\right|^{4}\ dxd\tau$ $\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}\
dx+\int_{0}^{t}\int_{\Omega}\rho_{\Delta t}g\cdot u_{\Delta
t}+\frac{1}{\mu}\left(\text{curl}\overline{B}_{\Delta
t}^{\prime}\times\overline{B}_{\Delta t}^{\prime}\right)\cdot u_{\Delta t}\
dxd\tau$ $\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}\
dx+\int_{0}^{t}\int_{\Omega}\rho_{\Delta t}g\cdot u_{\Delta t}\ dxd\tau$
$\displaystyle+\frac{\epsilon}{8}\int_{0}^{t}\int_{\Omega}\left|u_{\Delta
t}\right|^{4}\ dxd\tau+\frac{c}{\epsilon}\int_{0}^{t}B_{\Delta
t,\|\cdot\|}(\tau)d\tau+\int_{0}^{t}\int_{\Omega}\frac{\epsilon}{8\mu^{3}}\left|\operatorname{curl}\overline{B}_{\Delta
t}^{\prime}\right|^{4}\ dxd\tau,$ (78)

where the last estimate uses Hölder’s inequality, Young’s inequality and the
same inequality (76) as in the estimate of the corresponding term in the
induction equation. Adding the inequality (78) to the inequality (77) and
absorbing multiple terms from the right-hand side into the left-hand side, we
finally get the energy inequality

$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{\Delta t}(t)\left|u_{\Delta
t}(t)\right|^{2}+a\frac{\rho_{\Delta
t}^{\gamma}(t)}{\gamma-1}+\frac{\alpha\rho_{\Delta t}^{\beta}(t)}{\beta-1}\
dx+\frac{1}{2\mu}B_{\Delta
t,\|\cdot\|}(t)+\int_{0}^{t}\int_{\Omega}2\nu\left(\overline{\chi}_{\Delta
t}^{\prime}\right)|\mathbb{D}\left(u_{\Delta t}\right)|^{2}$
$\displaystyle+\lambda\left(\overline{\chi}_{\Delta
t}^{\prime}\right)\left|\operatorname{div}u_{\Delta
t}\right|^{2}+a\epsilon\gamma\rho_{\Delta
t}^{\gamma-2}\left|\nabla\rho_{\Delta
t}\right|^{2}+\alpha\epsilon\beta\rho_{\Delta
t}^{\beta-2}\left|\nabla\rho_{\Delta
t}\right|^{2}+\frac{3\epsilon}{4}\left|u_{\Delta t}\right|^{4}\ dxdt$
$\displaystyle+\frac{1}{\mu}\int_{0}^{t}\int_{\Omega}\frac{1}{\sigma\mu}\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{2}+\frac{5\epsilon}{8\mu^{2}}\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{4}+\epsilon\left|\Delta\overline{B}_{\Delta t}\right|^{2}\ dxd\tau$
$\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}\
dx+\frac{1}{2\mu}\int_{\Omega}\left|B_{0}\right|^{2}\
dx+\frac{\epsilon}{8}\Delta t\int_{\Omega}\left|u_{0}\right|^{4}\
dx+\frac{c}{\epsilon^{\frac{1}{3}}}\int_{0}^{T}\int_{\Omega}\left|\overline{J}_{\Delta
t}\right|^{\frac{4}{3}}\ dxd\tau$
$\displaystyle+\int_{0}^{t}\int_{\Omega}\rho_{\Delta t}g\cdot u_{\Delta t}\
dxd\tau+\int_{0}^{t}B_{\Delta t,\|\cdot\|}(\tau)d\tau$ $\displaystyle\leq$
$\displaystyle c+\int_{0}^{t}\int_{\Omega}\rho_{\Delta t}g\cdot u_{\Delta t}\
dxd\tau+\int_{0}^{t}B_{\Delta t,\|\cdot\|}(\tau)d\tau\quad\quad\forall
t\in[0,T],$ (79)

where the constant $c>0$ is independent of $\Delta t$ and $t$. In particular,
by use of the Gronwall Lemma and the estimates for the solution to the Neumann
problem for the density, c.f. Lemma 10.1, we find a constant $c>0$,
independent of $\Delta t$, such that the following bounds hold true:

$\displaystyle\left\|u_{\Delta t}\right\|_{C([0,T];V_{n})}+\left\|\rho_{\Delta
t}\right\|_{C([0,T];C^{2}(\overline{\Omega}))}+\left\|\partial_{t}\rho_{\Delta
t}\right\|_{C(\overline{Q})}\leq$ $\displaystyle c,$ (80)
$\displaystyle\left\|B_{\Delta
t}\right\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\left\|B_{\Delta
t}\right\|_{L^{2}(0,T;H^{2,2}(\Omega))}+\left\|\operatorname{curl}B_{\Delta
t}\right\|_{L^{4}(Q)}\leq$ $\displaystyle c,$ (81)
$\displaystyle\left\|\overline{B}_{\Delta
t}\right\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\left\|\overline{B}_{\Delta
t}\right\|_{L^{2}(0,T;H^{2,2}(\Omega))}+\left\|\operatorname{curl}\overline{B}_{\Delta
t}\right\|_{L^{4}(Q)}\leq$ $\displaystyle c,$ (82)
$\displaystyle\left\|\overline{B}_{\Delta
t}^{\prime}\right\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\left\|\overline{B}_{\Delta
t}^{\prime}\right\|_{L^{2}(0,T;H^{2,2}(\Omega))}+\left\|\operatorname{curl}\overline{B}_{\Delta
t}^{\prime}\right\|_{L^{4}(Q)}\leq$ $\displaystyle c.$ (83)

The bounds for the magnetic induction in $L^{\infty}(0,T;L^{2}(\Omega))$ in
(81)–(83) follow from the choice $t=k\Delta t$, $k=1,...,\frac{T}{\Delta t}$
in the energy inequality (79), for which it holds $B_{\Delta
t,\|\cdot\|}(t)=\|B_{\Delta t}^{k}\|_{L^{2}(\Omega)}^{2}$. For a bound of the
time derivative of $u_{\Delta t}$ we introduce the operator

$\displaystyle\mathcal{M}_{\rho_{\Delta t}(t)}:V_{n}\rightarrow
V_{n},\quad\quad\left\langle M_{\rho_{\Delta
t}(t)}v^{1},v^{2}\right\rangle:=\int_{\Omega}\rho_{\Delta t}(t)v^{1}\cdot
v^{2}\ dx\quad\quad\forall v^{1},v^{2}\in V_{n}$

and denote

$\displaystyle\mathcal{N}\left(\rho_{\Delta t},u_{\Delta t}\right):=$
$\displaystyle-\operatorname{div}\left(\rho_{\Delta t}u_{\Delta t}\otimes
u_{\Delta t}\right)-a\nabla\rho_{\Delta t}^{\gamma}-\alpha\nabla\rho_{\Delta
t}^{\beta}+\operatorname{div}\left(2\nu\left(\chi_{\Delta
t}^{k-1}\right)\mathbb{D}\left(u_{\Delta t}\right)\right)$ $\displaystyle+$
$\displaystyle\nabla\left(\lambda\left(\chi_{\Delta
t}^{k-1}\right)\operatorname{div}\left(u_{\Delta t}\right)\right)+\rho_{\Delta
t}g+\frac{1}{\mu}\left(\text{curl}B_{\Delta t}^{k-1}\times B_{\Delta
t}^{k-1}\right)-\epsilon\nabla u_{\Delta t}\nabla\rho_{\Delta
t}-\epsilon\left|u_{\Delta t}\right|^{2}u_{\Delta t}.$

Due to the uniform bound (80) the solution $\rho_{\Delta t}$ to the Neumann
problem (68), (70) is bounded away from $0$ uniformly with respect to $\Delta
t$, c.f. the estimate (163) in Lemma 10.1 below. Consequently,
$\mathcal{M}_{\rho_{\Delta t}}$ is invertible and from the momentum equation
(69) it is possible to derive the representation

$\displaystyle\partial_{t}u_{\Delta t}(t)=$
$\displaystyle\mathcal{M}_{\rho_{\Delta
t}(t)}^{-1}\mathcal{M}_{\partial_{t}\rho_{\Delta
t}(t)}\mathcal{M}_{\rho_{\Delta
t}(t)}^{-1}\left[P_{n}\left(\rho_{0}u_{0}\right)+\int_{0}^{t}P_{n}\left(\mathcal{N}\left(\rho_{\Delta
t},u_{\Delta t}\right)\right)dt\right]$
$\displaystyle+\mathcal{M}_{\rho_{\Delta
t}(t)}^{-1}\left[P_{n}\left(\mathcal{N}\left(\rho_{\Delta t},u_{\Delta
t}\right)\right)(t)\right].$

This together with the uniform bound of $\rho_{\Delta t}$ away from $0$ and
the uniform bounds (80)–(83) leads to the estimate

$\displaystyle\left\|\partial_{t}u_{\Delta
t}\right\|_{L^{2}(0,T;V_{n})}^{2}\leq c.$ (84)

For more details on the derivation of (84) we refer to [37, Section 7.7] and
in particular [37, Section 7.7.4.1]. The bounds (80)–(83) and (84) and the
Aubin-Lions Lemma imply the existence of functions

$\displaystyle 0\leq\rho$ $\displaystyle\in\left\\{\psi\in
C\left([0,T];H^{1,2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{2,2}(\Omega)\right):\ \partial_{t}\psi\in
L^{2}\left(Q\right),\
\left.\nabla\psi\cdot\operatorname{n}\right|_{\partial\Omega}=0\right\\},$
(85) $\displaystyle u$ $\displaystyle\in\left\\{\phi\in
C\left([0,T];V_{n}\right):\ \partial_{t}\phi\in
L^{2}\left(0,T;V_{n}\right)\right\\},$ $\displaystyle B$
$\displaystyle\in\left\\{b\in L^{\infty}\left(0,T;L^{2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{2,2}_{\operatorname{div}}(\Omega)\right):\
b\cdot\text{n}|_{\partial\Omega}=0\right\\}$ (86)

such that, after the extraction of a subsequence, it holds

$\displaystyle\rho_{\Delta t}\rightharpoonup\rho\ \ \quad$
$\displaystyle\text{in }L^{2}\left(0,T;H^{2,2}(\Omega)\right),$
$\displaystyle\rho_{\Delta t}$ $\displaystyle\rightarrow\rho\ \quad\text{in
}C\left([0,T];H^{1,2}(\Omega)\right),$ (87)
$\displaystyle\partial_{t}\rho_{\Delta t}\rightharpoonup\partial_{t}\rho\quad$
$\displaystyle\text{in }L^{2}\left(Q\right),\quad$ $\displaystyle\rho_{\Delta
t}^{\frac{\beta}{2}}$
$\displaystyle\rightharpoonup\rho^{\frac{\beta}{2}}\quad\text{in
}L^{2}\left(0,T;H^{1,2}(\Omega)\right),$ (88) $\displaystyle\quad u_{\Delta
t}\rightarrow u\ \ \quad$ $\displaystyle\text{in
}C\left([0,T];V_{n}\right),\quad$ $\displaystyle\partial_{t}u_{\Delta t}$
$\displaystyle\rightharpoonup\partial_{t}u\quad\text{in
}L^{2}\left(0,T;V_{n}\right),$ (89) $\displaystyle B_{\Delta t},\
\overline{B}_{\Delta t},\ \overline{B}^{\prime}_{\Delta
t}\buildrel\ast\over{\rightharpoonup}B\quad$ $\displaystyle\text{in
}L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right),\quad$ $\displaystyle
B_{\Delta t},\ \overline{B}_{\Delta t},\ \overline{B}^{\prime}_{\Delta t}$
$\displaystyle\rightharpoonup B\quad\ \ \text{in
}L^{2}\left(0,T;H^{2,2}\left(\Omega\right)\right).$

Here, the fact that the weak limits of the different interpolants coincide
follows from Lemma 10.2. Moreover the boundary conditions of the limit
functions in (85) and (86) follow directly from the corresponding boundary
conditions on the $\Delta t$-level, c.f. Proposition 4.1 and the definition of
the space $Y^{k}(S_{\Delta t})$ in (40). Furthermore, the external force
$\overline{J}_{\Delta t}$, discretized via (57), converges to its original
time-dependent counterpart,

$\displaystyle\overline{J}_{\Delta t}\rightarrow J\quad\text{in
}L^{q}(Q)\quad\quad\quad\forall 1\leq q<\infty.$

Finally, $X_{\Delta t}$, as the solution to the initial value problem (71),
satisfies the conditions of Lemma 10.3 which tells us that

$\displaystyle X_{\Delta t}\rightarrow X\quad$ $\displaystyle\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right),$
$\displaystyle\chi_{\Delta
t}\rightarrow\chi:=\textbf{db}_{S(\cdot)}(\cdot)\quad$ $\displaystyle\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right),$ (90)

where $S(t):=(X(t;O))^{\delta}$ and $X$ represents the solution to

$\displaystyle\frac{dX(t;x)}{dt}=R_{\delta}\left[u\right]\left(t,X(t;x)\right),\quad\text{for
}t\in[0,T],\quad\quad X(0;x)=x,\quad\text{for }x\in\mathbb{R}^{3}.$

###  5.2 Continuity equation

Due to the convergences (87)–(89) of the density and the velocity we can pass
to the limit in the continuity equation (68) and infer that the limit
functions $\rho$ and $u$ solve the initial value problem

$\displaystyle\partial_{t}\rho+\nabla\left(\rho u\right)=$
$\displaystyle\epsilon\Delta\rho\quad\text{a.e. in
}Q,\quad\quad\rho(0)=\rho_{0}.$

###  5.3 Induction equation

We first show convergence of the quantity $\overline{\tilde{u}}_{\Delta
t}^{\prime}$ from the mixed term in the discrete induction equation (73). We
fix an arbitrary point $(t,x)\in Q$ and, for each sufficiently small $\Delta
t>0$, we choose $k_{\Delta t}\in\\{2,...,\frac{T}{\Delta t}\\}$ such that
$t\in[(k_{\Delta t}-1)\Delta t,k_{\Delta t}\Delta t)$. It holds

$\displaystyle\left|\overline{\tilde{u}}_{\Delta
t}^{\prime}(t,x)-u(t,x)\right|=$ $\displaystyle\left|\frac{1}{\Delta
t}\int_{(k_{\Delta t}-2)\Delta t}^{(k_{\Delta t}-1)\Delta t}u_{\Delta
t,k_{\Delta t}-1}(\tau,x)d\tau-u(t,x)\right|$ $\displaystyle\leq$
$\displaystyle\sup_{\tau\in[(k_{\Delta t}-2)\Delta t,(k_{\Delta t}-1)\Delta
t]}\left|u_{\Delta t}(\tau-\Delta t,x)-u(t,x)\right|\rightarrow 0$ (91)

due to the uniform convergence (89) of $u_{\Delta t}$. Moreover, the uniform
bound of $u_{\Delta t}$ in (80) shows equiintegrability of
$|\overline{\tilde{u}}_{\Delta t}^{\prime}-u_{\Delta t}|^{q}$ for any $1\leq
q<\infty$. This together with the pointwise convergence (91) gives us the
conditions for the Vitali convergence theorem and we infer that

$\displaystyle\overline{\tilde{u}}^{\prime}_{\Delta t}\rightarrow
u\quad\text{in }L^{q}(Q)\quad\forall 1\leq q<\infty.$ (92)

Further, due to the uniform bounds (82), we find a function $z\in
L^{\frac{4}{3}}(Q)$ such that, possibly after the extraction of a suitable
subsequence, it holds

$\displaystyle\epsilon\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{2}\operatorname{curl}\overline{B}_{\Delta t}\rightharpoonup\epsilon
z\quad\text{in }L^{\frac{4}{3}}(Q).$ (93)

We remark that, since the quantity $\epsilon z$ will vanish from the system as
soon as we let $\epsilon$ tend to zero, there is no need to specify the form
of the limit function $z$ in (93). Now we test (73) by an arbitrary function
$b\in Y(S)$. This is possible, since $b$ is curl-free in an open neighbourhood
of $\overline{Q}^{s}(S)$ and so, by the uniform convergence (90) of the signed
distance function, it also satisfies

$\displaystyle\operatorname{curl}b(t)=0\quad\text{a.e. in }S_{\Delta
t}(k\Delta t)\quad\text{for a.a. }t\in[(k-1)\Delta t,k\Delta t],\
k=1,...,\frac{T}{\Delta t}$ (94)

for all sufficiently small $\Delta t>0$. Hence, for all such $t$, $k$ and
$\Delta t$ it holds $b(t)\in W^{k}(S_{\Delta t})$ and thus $b$ is indeed an
admissible test function in (73), c.f. (72). Using the strong convergence (92)
of $\overline{\tilde{u}}^{\prime}_{\Delta t}$, we can pass to the limit with
respect to $\Delta t\rightarrow 0$ in the resulting equation and obtain the
induction equation in the limit,

$\displaystyle-\int_{0}^{T}\int_{\Omega}B\cdot\partial_{t}b\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}B+u\times
B+\frac{1}{\sigma}J-\frac{\epsilon}{\mu^{2}}z\right]\cdot\operatorname{curl}b$
$\displaystyle-\epsilon\operatorname{curl}\left(\operatorname{curl}B\right)\cdot\operatorname{curl}\left(\operatorname{curl}b\right)\
dxdt\quad\quad\forall b\in Y(S).$

###  5.4 Momentum equation

In order to pass to the limit in the momentum equation, it remains to show
convergence of the piecewise constant Lorentz force. This is achieved by the
same arguments as in the incompressible case, c.f. [1, Section 4.2]: The
uniform bounds (82) and (83) allow us to extract suitable subsequences and
find functions $z_{1},\ z_{2}\in L^{2}(Q)$ with the properties that

$\displaystyle\operatorname{curl}\overline{B}^{\prime}_{\Delta
t}\times\overline{B}^{\prime}_{\Delta t}\rightharpoonup z_{1}\quad$
$\displaystyle\text{in }L^{2}(Q),\quad$
$\displaystyle\operatorname{curl}\overline{B}_{\Delta
t}\times\overline{B}^{\prime}_{\Delta t}$ $\displaystyle\rightharpoonup
z_{2}\quad\text{in }L^{2}(Q).$ (95)

Our goal is to identify the limit functions $z_{1}$ and $z_{2}$ as

$\displaystyle z_{1}=z_{2}=\operatorname{curl}B\times B\quad\text{a.e. in }Q.$
(96)

Since, according to (51), it holds

$\displaystyle\left|Q\right|=\left|Q^{s}(S)\bigcup Q^{f}(S)\right|,$

it is sufficient to do so in $Q^{s}(S)$ and $Q^{f}(S)$. Focussing at first on
the fluid part of the domain, we consider any $c,d\in[0,T]$ and any ball
$U\subset\Omega$ such that $\overline{I\times U}:=\overline{(c,d)\times
U}\subset Q^{f}(S)$. In particular, any function $b\in
L^{4}(I;H_{0}^{2,2}(U))$, extended by $0$ outside of $\overline{I\times U}$,
is curl-free in an open neighbourhood of $\overline{Q}^{s}(S)$. Hence, the
uniform convergence (90) of the signed distance function implies that, for
sufficiently small $\Delta t>0$, any such $b$ also satisfies the curl-free
condition (94) in $S_{\Delta t}(k\Delta t)$, $k=1,...,\frac{T}{\Delta t}$, and
therefore constitutes an admissible test function in the induction equation
(73). Using it as such we deduce the uniform bound

$\displaystyle\left\|\frac{\overline{B}^{\prime}_{\Delta
t}(\cdot)-\overline{B}^{\prime}_{\Delta t}(\cdot-\Delta t)}{\Delta
t}\right\|_{L^{\frac{4}{3}}\left(a+\Delta
t,d;H^{-2,2}(U)\right)}\leq\left\|\frac{\overline{B}_{\Delta
t}(\cdot)-\overline{B}_{\Delta t}(\cdot-\Delta t)}{\Delta
t}\right\|_{L^{\frac{4}{3}}(a,d;H^{-2,2}(U))}\leq c.$

Since, by (83), $\overline{B}^{\prime}_{\Delta t}$ is bounded uniformly in
$L^{2}(Q)$, we can thus apply the discrete Aubin-Lions Lemma [12, Theorem 1]
and conclude that

$\displaystyle\overline{B}^{\prime}_{\Delta t}\rightarrow B\quad\text{in
}L^{2}\left(I;H^{-1,2}(U)\right).$

This strong convergence implies that

$\displaystyle z_{1}=z_{2}=\operatorname{curl}B\times B\quad\text{a.e. in
}I\times U\ \text{and hence in }Q^{f}(S).$ (97)

In order to show the equation (96) also in the solid region, we consider
another arbitrary pair of an interval $I\subset(0,T)$ and a ball
$U\subset\Omega$, this time satisfying $\overline{I\times U}\subset Q^{s}(S)$.
From the fact that $\operatorname{curl}B_{\Delta t}^{k}=0$ in $S_{\Delta
t}^{k}\bigcap\Omega$ for $k=1,...,\frac{T}{\Delta t}$, and the uniform
convergence (90) of the signed distance function it follows that, for any
sufficiently small $\Delta t>0$,

$\displaystyle\operatorname{curl}\overline{B}^{\prime}_{\Delta
t}=\operatorname{curl}\overline{B}_{\Delta t}=0\quad\text{a.e. in }I\times U.$

Thus, letting $\Delta t$ tend to zero, we infer that

$\displaystyle z_{1}=z_{2}=0=\operatorname{curl}B=\operatorname{curl}B\times
B\quad\text{a.e. in }I\times U\ \text{and hence in }Q^{s}(S).$ (98)

In combination with the corresponding equality (97) in $Q^{f}(S)$, we have
therefore shown the desired identification (96) of $z_{1}$ and $z_{2}$. We
remark that (98) moreover shows that the solid region $Q^{s}(S)$ in the limit
is again insulating. Next we exploit the uniform convergence (90) of the
signed distance function together with the definition (53) of the variable
viscosity coefficients as smooth functions of $\overline{\chi}_{\Delta
t}^{\prime}$ to infer that

$\displaystyle\nu\left(\overline{\chi}^{\prime}_{\Delta
t}\right)\rightarrow\nu\left(\chi\right)\quad\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right),\quad\quad\lambda\left(\overline{\chi}^{\prime}_{\Delta
t}\right)\rightarrow\lambda\left(\chi\right)\quad\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right).$ (99)

This, in combination with the convergence of the Lorentz force, c.f. (95),
(96) and the convergences (87)–(89) of $\rho_{\Delta t}$ and $u_{\Delta t}$,
allows us to pass to the limit in the momentum equation (69) and infer the
equation

$\displaystyle\int_{0}^{T}\int_{\Omega}\partial_{t}\left(\rho
u\right)\cdot\phi\ dxdt$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left(\rho u\otimes
u\right):\mathbb{D}(\phi)+\left(a\rho^{\gamma}+\alpha\rho^{\beta}\right)\operatorname{div}\phi-2\nu\left(\chi\right)\mathbb{D}\left(u\right):\mathbb{D}(\phi)-\lambda\left(\chi\right)\operatorname{div}\left(u\right)\operatorname{div}\phi$
$\displaystyle+\rho g\cdot\phi+\frac{1}{\mu}\left(\text{curl}B\times
B\right)\cdot\phi-\epsilon\left|u\right|^{2}u\cdot\phi-\epsilon\left(\nabla
u\nabla\rho\right)\cdot\phi\ dxdt\quad\forall\phi\in
C\left([0,T];V_{n}\right).$

###  5.5 Limit passage in the energy inequality

We choose an arbitrary $t\in(0,T]$ and $k\in\\{1,...,\frac{T}{\Delta t}\\}$,
such that $t=k\Delta t-\xi$ for some $\xi\in[0,\Delta t)$. Subsequently, we
sum the inequality (74) over $l=1,...,k$ and add the first inequality in (78)
to obtain

$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{\Delta t}(t)\left|u_{\Delta
t}(t)\right|^{2}+\frac{a}{\gamma-1}\rho_{\Delta
t}^{\gamma}(t)+\frac{\alpha}{\beta-1}\rho_{\Delta
t}^{\beta}(t)+\frac{1}{2\mu}\left|\overline{B}_{\Delta t}(t)\right|^{2}\ dx$
$\displaystyle+\int_{0}^{t}\int_{\Omega}2\nu\left(\overline{\chi}_{\Delta
t}^{\prime}\right)|\mathbb{D}\left(u_{\Delta
t}\right)|^{2}+\lambda\left(\overline{\chi}_{\Delta
t}^{\prime}\right)\left|\operatorname{div}u_{\Delta
t}\right|^{2}+a\epsilon\gamma\rho_{\Delta
t}^{\gamma-2}\left|\nabla\rho_{\Delta
t}\right|^{2}+\alpha\epsilon\beta\rho_{\Delta
t}^{\beta-2}\left|\nabla\rho_{\Delta t}\right|^{2}+\epsilon\left|u_{\Delta
t}\right|^{4}$
$\displaystyle+\frac{1}{\sigma\mu^{2}}\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{2}+\frac{\epsilon}{\mu^{3}}\left|\operatorname{curl}\overline{B}_{\Delta
t}\right|^{4}+\frac{\epsilon}{\mu}\left|\Delta\overline{B}_{\Delta
t}\right|^{2}\ dxd\tau$ $\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}+\frac{1}{2\mu}\left|B_{0}\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}\rho_{\Delta t}g\cdot u_{\Delta
t}+\frac{1}{\mu}\left(\text{curl}\overline{B}_{\Delta
t}^{\prime}\times\overline{B}_{\Delta t}^{\prime}\right)\cdot u_{\Delta t}$
$\displaystyle+\frac{1}{\mu}\left(\overline{\tilde{u}}_{\Delta
t}^{\prime}\times\overline{B}_{\Delta
t}^{\prime}\right)\cdot\operatorname{curl}\overline{B}_{\Delta
t}+\frac{1}{\sigma\mu}\overline{J}_{\Delta
t}\cdot\operatorname{curl}\overline{B}_{\Delta t}\ dxd\tau+c\left(\Delta
t\right)^{\frac{1}{2}}.$ (100)

Here, on the left-hand side we drop the term $\rho_{\Delta
t}^{\gamma-2}\left|\nabla\rho_{\Delta t}\right|^{2}$ and pass to the limit by
exploiting, in particular, the weak lower semicontinuity of norms and the
strong convergence (99) of the variable viscosity coefficients. Moreover, on
the right-hand side of (100) we can carry out the limit passage by exploiting
the convergence (95), (96) of the Lorentz force. Altogether we obtain

$\displaystyle\int_{\Omega}\frac{1}{2}\rho(t)\left|u(t)\right|^{2}+\frac{a}{\gamma-1}\rho^{\gamma}(t)+\frac{\alpha}{\beta-1}\rho^{\beta}(t)+\frac{1}{2\mu}\left|B(t)\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}2\nu\left(\chi\right)|\mathbb{D}\left(u\right)|^{2}+\lambda\left(\chi\right)\left|\operatorname{div}u\right|^{2}$
$\displaystyle+\alpha\epsilon\beta\rho^{\beta-2}\left|\nabla\rho\right|^{2}+\epsilon\left|u\right|^{4}+\frac{1}{\sigma\mu^{2}}\left|\operatorname{curl}B\right|^{2}+\frac{\epsilon}{\mu}\left|\Delta
B\right|^{2}+\frac{\epsilon}{\mu^{3}}\left|z\right|^{4}\ dxd\tau$
$\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}+\frac{1}{2\mu}\left|B_{0}\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}\rho g\cdot
u+\frac{1}{\sigma\mu}J\cdot\operatorname{curl}B\ dxd\tau\quad\text{for
}t\in[0,T].$

Hence we have proved

###### Proposition 5.1

Let $n\in\mathbb{N}$, $\eta,\epsilon,\alpha>0$, $\beta>\max\\{4,\gamma\\}$
sufficiently large and let the assumptions of Theorem 2.1 be satisfied.
Furthermore, let

$\displaystyle\rho_{0}\in$ $\displaystyle
C^{2,\frac{1}{2}}\left(\overline{\Omega}\right),\quad\quad(\rho u)_{0}\in
C^{2}\left(\overline{\Omega}\right),\quad\quad u_{0}:=P_{n}\left(\frac{(\rho
u)_{0}}{\rho_{0}}\right)\in V_{n},\quad\quad B_{0}\in H^{2,2}(\Omega),$
$\displaystyle
0<\alpha\leq\rho_{0}\leq\alpha^{-\frac{1}{2\beta}},\quad\quad\left.\nabla\rho_{0}\cdot\operatorname{n}\right|_{\partial\Omega}=0,\quad\quad\operatorname{div}B_{0}=0,\quad\quad\left.B_{0}\cdot\operatorname{n}\right|_{\partial\Omega}=0.$

Under these conditions there exist functions
$X_{n}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, $z_{n}\in
L^{\frac{4}{3}}((0,T)\times\Omega)$ and

$\displaystyle 0\leq\rho_{n}$ $\displaystyle\in\left\\{\psi\in
C\left([0,T];H^{1,2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{2,2}(\Omega)\right):\ \partial_{t}\psi\in
L^{2}\left(Q\right),\
\left.\nabla\psi\cdot\operatorname{n}\right|_{\partial\Omega}=0\right\\},$
(101) $\displaystyle u_{n}$ $\displaystyle\in\left\\{\phi\in
C\left([0,T];V_{n}\right):\ \partial_{t}\phi\in
L^{2}\left(0,T;V_{n}\right)\right\\},$ $\displaystyle B_{n}$
$\displaystyle\in\left\\{b\in L^{\infty}\left(0,T;L^{2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{2,2}_{\operatorname{div}}(\Omega)\right):\
\operatorname{curl}b=0\ \text{in }Q^{s}(S_{n}),\
b\cdot\operatorname{n}|_{\partial\Omega}=0\right\\}$ (102)

for $S_{n}=S_{n}(\cdot)=(X_{n}(\cdot;O))^{\delta}$, which satisfy

$\displaystyle\frac{dX_{n}(t;x)}{dt}=R_{\delta}\left[u_{n}\right]\left(t,X_{n}(t;x)\right),\quad$
$\displaystyle\text{for }t\in[0,T],\quad\quad X_{n}(0;x)=x,\quad\text{for
}x\in\mathbb{R}^{3},$
$\displaystyle\partial_{t}\rho_{n}+\operatorname{div}\left(\rho_{n}u_{n}\right)$
$\displaystyle=\epsilon\Delta\rho_{n}\quad\text{a.e. in }Q$ (103)

and

$\displaystyle\int_{0}^{T}\int_{\Omega}\partial_{t}\left(\rho_{n}u_{n}\right)\cdot\phi\
dxdt=$ $\displaystyle\int_{0}^{T}\int_{\Omega}\left(\rho_{n}u_{n}\otimes
u_{n}\right):\mathbb{D}(\phi)+\left(a\rho_{n}^{\gamma}+\alpha\rho_{n}^{\beta}\right)\operatorname{div}\phi-2\nu\left(\chi_{n}\right)\mathbb{D}(u_{n}):\mathbb{D}(\phi)$
$\displaystyle-\lambda\left(\chi_{n}\right)\operatorname{div}u_{n}\operatorname{div}\phi+\rho_{n}g\cdot\phi+\frac{1}{\mu}\left(\operatorname{curl}B_{n}\times
B_{n}\right)\cdot\phi-\epsilon\left|u_{n}\right|^{2}u_{n}\cdot\phi$
$\displaystyle-\epsilon\left(\nabla u_{n}\nabla\rho_{n}\right)\cdot\phi\
dxdt,$ (104) $\displaystyle-\int_{0}^{T}\int_{\Omega}B_{n}\cdot\partial_{t}b\
dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}B_{n}+u_{n}\times
B_{n}+\frac{1}{\sigma}J-\frac{\epsilon}{\mu^{2}}z_{n}\right]\cdot\operatorname{curl}b$
$\displaystyle-\epsilon\operatorname{curl}\left(\operatorname{curl}B_{n}\right):\operatorname{curl}\left(\operatorname{curl}b\right)\
dxdt,$ (105)

where $\chi_{n}(t,x):=\textbf{db}_{S_{n}(t)}(x)$, for any $\phi\in
C([0,T];V_{n})$ and any $b\in Y(S_{n})$. Further, these functions satisfy the
initial conditions

$\displaystyle\rho_{n}(0)=\rho_{0},\quad\quad u_{n}(0)=u_{0},\quad\quad
B_{n}(0)=B_{0},$

of which the latter identity can be understood in the sense of (35), as well
as the energy inequality

$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{n}(t)\left|u_{n}(t)\right|^{2}+\frac{a}{\gamma-1}\rho_{n}^{\gamma}(t)+\frac{\alpha}{\beta-1}\rho_{n}^{\beta}(t)+\frac{1}{2\mu}\left|B_{n}(t)\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}2\nu\left(\chi_{n}\right)|\mathbb{D}\left(u_{n}\right)|^{2}$
$\displaystyle+\lambda\left(\chi_{n}\right)\left|\operatorname{div}u_{n}\right|^{2}+\alpha\epsilon\beta\rho_{n}^{\beta-2}\left|\nabla\rho_{n}\right|^{2}+\epsilon\left|u_{n}\right|^{4}+\frac{1}{\sigma\mu^{2}}\left|\operatorname{curl}B_{n}\right|^{2}+\frac{\epsilon}{\mu}\left|\Delta
B_{n}\right|^{2}+\frac{\epsilon}{\mu^{3}}\left|z_{n}\right|^{\frac{4}{3}}\
dxd\tau$ $\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}+\frac{1}{2\mu}\left|B_{0}\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}\rho_{n}g\cdot
u_{n}+\frac{1}{\sigma\mu}J\cdot\operatorname{curl}B_{n}\ dxd\tau$ (106)

for almost all $t\in[0,T]$.

##  6 Limit passage in the Galerkin method

Next, we pass to the limit with respect to $n\rightarrow\infty$. Using
Lebesgue interpolation, we infer from the energy inequality (106) the
existence of a constant $c(\epsilon,\alpha)>0$, independent of $n$, such that

$\displaystyle\left\|\rho_{n}\right\|_{L^{\frac{5}{3}\beta}(Q)}\leq\left\|\rho_{n}\right\|^{\frac{2}{5}}_{L^{\infty}\left(0,T;L^{\beta}(\Omega)\right)}\left\|\rho_{n}\right\|^{\frac{3}{5}}_{L^{\beta}\left(0,T;L^{3\beta}(\Omega)\right)}\leq
c(\epsilon,\alpha).$ (107)

Moreover, from the classical $L^{p}$-$L^{q}$ regularity results for parabolic
equations, c.f. [37, Lemma 7.37, Lemma 7.38, Section 7.8.2], we infer that
$\rho_{n}$ as the solution to the regularized continuity equation (103)
satisfies the estimates

$\displaystyle\epsilon\left\|\nabla\rho_{n}\right\|_{L^{r}(Q)}+\epsilon\left\|\partial_{t}\rho_{n}\right\|_{L^{\tilde{r}}(Q)}+\epsilon^{2}\left\|\Delta\rho_{n}\right\|_{L^{\tilde{r}}(Q)}\leq
c$ (108)

for

$\displaystyle
r:=\frac{10\beta-6}{3\beta+3}>2,\quad\tilde{r}:=\frac{5\beta-3}{4\beta}>1\quad\quad\forall\beta>6$
(109)

and a constant $c>0$ independent of $n$, $\eta$ and $\epsilon$. The uniform
bounds (107),(108) and the energy inequality (106), together with the Aubin-
Lions Lemma, allow us to extract suitable subsequences and find functions
$z\in L^{\frac{4}{3}}(Q)$ and

$\displaystyle u\in$ $\displaystyle L^{2}(0,T;H_{0}^{1,2}(\Omega)),$ (110)
$\displaystyle 0\leq\rho\in$ $\displaystyle\bigg{\\{}\psi\in
L^{\infty}\left(0,T;L^{\beta}(\Omega)\right)\bigcap
L^{r}\left(0,T;W^{1,r}(\Omega)\right)\bigcap
L^{\tilde{r}}\left(0,T;W^{2,\tilde{r}}(\Omega)\right):$
$\displaystyle\partial_{t}\psi\in
L^{\tilde{r}}\left((0,T)\times\Omega\right),\
\left.\nabla\psi\cdot\operatorname{n}\right|_{\partial\Omega}=0\bigg{\\}}\subset
C\left([0,T];L^{2}(\Omega)\right),$ (111) $\displaystyle B\in$
$\displaystyle\bigg{\\{}b\in L^{\infty}\left(0,T;L^{2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{2,2}_{\operatorname{div}}(\Omega)\right):\
b\cdot\text{n}|_{\partial\Omega}=0\bigg{\\}}$ (112)

with the properties that

$\displaystyle\rho_{n}\rightarrow\rho\quad$ $\displaystyle\text{in
}L^{\beta}\left(Q\right),\quad\quad$ $\displaystyle\rho_{n}$
$\displaystyle\rightarrow\rho\quad\ \ \text{in
}L^{2}\left(0,T;H^{1,2}(\Omega)\right)$ (113)
$\displaystyle\rho_{n}\rightharpoonup\rho\quad$ $\displaystyle\text{in
}L^{\tilde{r}}\left(0,T;W^{2,\tilde{r}}(\Omega)\right),\quad\quad$
$\displaystyle\partial_{t}\rho_{n}$
$\displaystyle\rightharpoonup\partial_{t}\rho\quad\text{in
}L^{\tilde{r}}\left(Q\right),$ (114) $\displaystyle u_{n}\rightharpoonup
u\quad$ $\displaystyle\text{in }L^{2}(0,T;H^{1,2}(\Omega)),\quad\quad$
$\displaystyle B_{n}$ $\displaystyle\buildrel\ast\over{\rightharpoonup}B\quad\
\ \text{in }L^{\infty}(0,T;L^{2}(\Omega)),$ (115) $\displaystyle
B_{n}\rightharpoonup B\quad$ $\displaystyle\text{in
}L^{2}(0,T;H^{2,2}(\Omega)),\quad\quad$ $\displaystyle z_{n}$
$\displaystyle\rightharpoonup z\quad\quad\text{in
}L^{\frac{4}{3}}\left(Q\right).$

The boundary conditions of the limit functions in (110)–(112) follow directly
from the corresponding boundary conditions (101) and (102) of $\nabla\rho_{n}$
and $B_{n}$ and the fact that $u_{n}\in V_{n}$ vanishes on $\partial\Omega$
for all $n\in\mathbb{N}$. Finally, the initial value problem (71), solved by
$X_{n}$, yields that the conditions of Lemma 10.3 are satisfied. Hence

$\displaystyle X_{n}\rightarrow X\quad$ $\displaystyle\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right),$
$\displaystyle\chi_{n}\rightarrow\chi:=\textbf{db}_{S(\cdot)}(\cdot)\quad$
$\displaystyle\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right),$ (116)

where $S(t):=(X(t;O))^{\delta}$ and $X$ denotes the unique solution to the
initial value problem

$\displaystyle\frac{dX(t;x)}{dt}=R_{\delta}\left[u\right]\left(t,X(t;x)\right),\quad\text{for
}t\in[0,T],\quad\quad X(0;x)=x,\quad\text{for }x\in\mathbb{R}^{3}.$

###  6.1 Continuity equation

The convergences (113)–(115) of $\rho_{n}$ and $u_{n}$ allow us to pass to the
limit in the continuity equation (103). Consequently, the limit functions
$\rho$ and $u$ satisfy the continuity equation

$\displaystyle\partial_{t}\rho+\operatorname{div}\left(\rho u\right)=$
$\displaystyle\epsilon\Delta\rho\quad\text{a.e. in }Q,\quad\rho(0)=\rho_{0}.$

###  6.2 Induction equation

At this stage - as well as in the later sections - the limit passage in the
induction equation does not differ from the incompressible case, c.f. [1]. For
the convenience of the reader, we present the arguments here: We begin by
making sure that test functions $b\in Y(S,X)$ for the limit equation are also
admissible on the $n$-level. To this end we fix some arbitrary $\omega>0$.
Then the uniform convergence (116) of the signed distance function implies the
existence of $N=N(\omega)>0$, such that for any function $b\in Y(S)$ with
$\operatorname{curl}b=0$ in an $\omega$-neighbourhood of $\overline{Q}^{s}(S)$
it also holds

$\displaystyle\operatorname{curl}b=0\quad\text{in an
}\frac{\omega}{2}\text{-neighbourhood of
}\overline{Q}^{s}(S_{n})\quad\text{and thus}\quad b\in
Y\left(S_{n}\right)\quad\forall n\geq N.$

In particular, for any interval $I\subset(0,T)$ and any ball $U\subset\Omega$,
such that $\overline{I\times U}\subset Q^{f}(S)$ there exists $N=N(I\times
U)>0$ such that

$\displaystyle b\in Y\left(S_{n}\right)\quad\quad\forall n\geq N,\
b\in\mathcal{D}(I\times U),$ (117)

where $b$ has been extended by $0$ outside of $I\times U$. Next, an
interpolation between $L^{\infty}(0,T;L^{2}(\Omega))$ and
$L^{2}(0,T;L^{6}(\Omega))$ provides a uniform bound of $B_{n}$ in $L^{3}(Q)$.
In combination with the bounds of $u_{n}$ and $\operatorname{curl}B_{n}$ in
$L^{2}(Q)$ this yields the existence of functions $z_{3},z_{4}\in
L^{\frac{6}{5}}(Q)$ such that, for a suitable subsequence, it holds

$\displaystyle\operatorname{curl}B_{n}\times B_{n}\rightharpoonup
z_{3}\quad\text{in }L^{\frac{6}{5}}\left(Q\right),\quad\quad u_{n}\times
B_{n}\rightharpoonup z_{4}\quad\text{in }L^{\frac{6}{5}}\left(Q\right).$ (118)

With the aim of identifying $z_{3}$ and $z_{4}$ we pick an arbitrary interval
$I\subset(0,T)$ and an arbitrary ball $U\subset\Omega$ with the property
$\overline{I\times U}\subset Q^{f}(S)$. From (117) we know that for any
sufficiently large $n\in\mathbb{N}$ the induction equation (105) may be tested
by all functions of the form $\psi b$, where $\psi\in\mathcal{D}(I)$,
$b\in\mathcal{D}(U)$. Doing so, we obtain the dual estimate

$\displaystyle\left\|\partial_{t}\int_{U}B_{n}\cdot
bdx\right\|_{L^{\frac{4}{3}}(I)}\leq c$

with a constant $c>0$ depending on $b$ but not on $n$. This, together with the
Arzelá-Ascoli theorem, yields

$\displaystyle B_{n}\rightarrow B\quad\text{in
}C_{\text{weak}}\left(\overline{I};L^{2}(U)\right)\quad\text{and hence in
}L^{2}\left(I;H^{-1,2}(U)\right).$ (119)

Combining this with the weak convergences of $u_{n}$ and $B_{n}$ in
$L^{2}(0,T;H^{1,2}(\Omega))$ we realize that

$\displaystyle z_{3}=\operatorname{curl}B\times B,\quad z_{4}=u\times
B\quad\quad\text{a.e. in }Q^{f}(S).$

Due to the uniform convergence (116) of the signed distance function and the
fact that $\operatorname{curl}B_{n}=0$ in $Q^{s}(S_{n})$, the former one of
these identities also holds true in $Q^{s}(S)$,

$\displaystyle z_{3}=0=\operatorname{curl}B=\operatorname{curl}B\times
B\quad\text{a.e. in }Q^{s}(S).$

Since moreover test functions $b\in Y(S)$ are curl-free in $Q^{s}(S)$, we end
up with

$\displaystyle z_{3}=\operatorname{curl}B\times B,\quad
z_{4}\cdot\operatorname{curl}b=u\times
B\cdot\operatorname{curl}b\quad\quad\text{a.e. in }Q$ (120)

for any $b\in Y(S)$. The relations (118) and (120) allow us to pass to the
limit in the mixed term of the induction equations and hence, letting $n$ tend
to infinity in (105), we see that

$\displaystyle-\int_{0}^{T}\int_{\Omega}B\cdot\partial_{t}b\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}B+u\times
B+\frac{1}{\sigma}J-\frac{\epsilon}{\mu^{2}}z\right]\cdot\operatorname{curl}b$
$\displaystyle-\epsilon\operatorname{curl}\left(\operatorname{curl}B\right)\cdot\operatorname{curl}\left(\operatorname{curl}b\right)\
dxdt$

for any $b\in Y(S)$. Finally, for any $b\in\mathcal{D}(\Omega)$ with
$\operatorname{curl}b=0$ in a neighbourhood of $\overline{S}_{0}$, we can
argue similarly as in the derivation of the $C_{\text{weak}}$-convergence
(119) of $B_{n}$ to deduce that

$\displaystyle\int_{\Omega}B_{n}(\cdot,x)\cdot b(x)\
dx\rightarrow\int_{\Omega}B(\cdot,x)\cdot b(x)\ dx\quad\text{in
}C\left(\left[0,t_{0}\right]\right),$

for some small $t_{0}=t_{0}(b)>0$. This yields the initial condition
$B(0)=B_{0}$ in the sense of (35).

###  6.3 Momentum equation

By the same methods as for the (purely mechanical) compressible Navier-Stokes
system, c.f. [37, Section 7.8.2], we derive strong convergence of the momentum
function in the Galerkin limit: Recalling that $P_{n}$ denotes the orthogonal
projection of $L^{2}(\Omega)$ onto $V_{n}$, we test the momentum equation
(104) — after a density argument — by $P_{n}(\phi)$ for an arbitrary function
$\phi\in L^{\frac{5\beta-3}{\beta-3}}(0,T;H_{0}^{2,2}(\Omega))$. Since in
particular

$\displaystyle\left|\int_{0}^{T}\int_{\Omega}\epsilon\left|u_{n}\right|^{2}u_{n}\cdot
P_{n}(\phi)\
dxdt\right|\leq\epsilon\left\|\left|u_{n}\right|^{3}\right\|_{L^{\frac{4}{3}}(Q)}\left\|P_{n}(\phi)\right\|_{L^{4}(Q)}\leq\epsilon\left\|u_{n}\right\|^{3}_{L^{4}(Q)}\left\|\phi\right\|_{L^{\frac{5\beta-3}{\beta-3}}(0,T;H^{2,2}(\Omega))},$

this results in the dual estimate

$\displaystyle\left\|\partial_{t}P_{n}\left(\rho_{n}u_{n}\right)\right\|_{L^{\tilde{r}}\left(0,T;H^{-2,2}(\Omega)\right)}=\left\|\partial_{t}P_{n}\left(\rho_{n}u_{n}\right)\right\|_{L^{\frac{5\beta-3}{4\beta}}\left(0,T;H^{-2,2}(\Omega)\right)}\leq
c.$

Consequently, from the Aubin-Lions Lemma, we conclude that

$\displaystyle P_{n}(\rho_{n}u_{n})\rightarrow\rho u\quad\text{in
}L^{2}\left(0,T;H^{-1,2}(\Omega)\right)\quad\text{and
hence}\quad\rho_{n}u_{n}\rightarrow\rho u\quad\text{in
}L^{2}\left(0,T;H^{-1,2}(\Omega)\right).$

Therefore, as $u_{n}$ converges weakly in $L^{2}(0,T;H^{1,2}(\Omega))$, we
infer that, for example,

$\displaystyle\rho_{n}u_{n}\otimes u_{n}\rightharpoonup\rho u\otimes
u\quad\text{in }L^{2}\left(0,T;L^{\frac{6\beta}{4\beta+3}}(\Omega)\right).$

Moreover, the bound of $u_{n}$ in $L^{4}(Q)$ implies the existence of some
$\tilde{z}\in L^{\frac{4}{3}}(Q)$ such that, for a chosen subsequence, it
holds

$\displaystyle\epsilon\left|u_{n}\right|^{2}u_{n}\rightharpoonup\epsilon\tilde{z}\quad\text{in
}L^{\frac{4}{3}}(Q).$

Using further the uniform convergence (116) of the signed distance function
for passing to the limit in the variable viscosity coefficients and the
relations (118), (120) which allow us to pass to the limit in the Lorentz
force, we can now let $n$ tend to infinity in the momentum equation (104) and
infer that

$\displaystyle-\int_{0}^{T}\int_{\Omega}\rho u\cdot\partial_{t}\phi\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left(\rho u\otimes
u\right):\mathbb{D}(\phi)+\left(a\rho^{\gamma}+\alpha\rho^{\beta}\right)\operatorname{div}\phi-2\nu\left(\chi\right)\mathbb{D}(u):\mathbb{D}(\phi)$
$\displaystyle-\lambda\left(\chi\right)\operatorname{div}u\operatorname{div}\phi+\rho
g\cdot\phi+\frac{1}{\mu}\left(\operatorname{curl}B\times
B\right)\cdot\phi-\epsilon\tilde{z}\cdot\phi-\epsilon\left(\nabla
u\nabla\rho\right)\cdot\phi\ dxdt$ (121)

for any $\phi\in C_{0}^{1}([0,T];V_{N})$ with fixed $N\in\mathbb{N}$. Since
$\bigcup_{n=1}^{\infty}V_{n}$ is dense in $H_{0}^{1,2}(\Omega)$, we finally
conclude that (121) also holds true for any $\phi\in\mathcal{D}(Q)$. Using
further the weak lower semicontinuity of norms to pass to the limit in both
the energy inequality (106) and the uniform bounds (108) we have proved the
following proposition:

###### Proposition 6.1

Let $\eta,\epsilon,\alpha>0$, $\beta>\max\\{4,\gamma\\}$ sufficiently large
and let the assumptions of Theorem 2.1 be satisfied. Furthermore, let

$\displaystyle\rho_{0}\in
C^{2,\frac{1}{2}}\left(\overline{\Omega}\right),\quad\quad(\rho u)_{0}\in
C^{2}$ $\displaystyle\left(\overline{\Omega}\right),\quad\quad B_{0}\in
H^{2,2}(\Omega),$ (122) $\displaystyle
0<\alpha\leq\rho_{0}\leq\alpha^{-\frac{1}{2\beta}},\quad\quad\left.\nabla\rho_{0}\cdot\operatorname{n}\right|_{\partial\Omega}=0,$
$\displaystyle\quad\quad\operatorname{div}B_{0}=0,\quad\quad\left.B_{0}\cdot\operatorname{n}\right|_{\partial\Omega}=0.$
(123)

Under these conditions there exist functions
$X_{\eta}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, $u_{\eta}\in
L^{2}(0,T;H_{0}^{1,2}(\Omega))$, $z_{\eta},\tilde{z}_{\eta}\in
L^{\frac{4}{3}}(Q)$ and

$\displaystyle 0\leq\rho_{\eta}\in$ $\displaystyle\bigg{\\{}\psi\in
L^{\infty}\left(0,T;L^{\beta}(\Omega)\right)\bigcap
L^{r}\left(0,T;W^{1,r}(\Omega)\right)\bigcap
L^{\tilde{r}}\left(0,T;W^{2,\tilde{r}}(\Omega)\right):$
$\displaystyle\partial_{t}\psi\in L^{\tilde{r}}\left(Q\right),\
\left.\nabla\psi\cdot\operatorname{n}\right|_{\partial\Omega}=0\bigg{\\}},$
(124) $\displaystyle B_{\eta}\in$ $\displaystyle\left\\{b\in
L^{\infty}\left(0,T;L^{2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{2,2}_{\operatorname{div}}(\Omega)\right):\
\operatorname{curl}b=0\ \text{in }Q^{s}\left(S_{\eta}\right),\
b\cdot\operatorname{n}|_{\partial\Omega}=0\right\\},$

for $r>2$, $\tilde{r}>1$ as in (109) and
$S_{\eta}=S_{\eta}(\cdot)=(X_{\eta}(\cdot;O))^{\delta}$, which satisfy

$\displaystyle\frac{dX_{\eta}(t;x)}{dt}=R_{\delta}\left[u_{\eta}\right]\left(t,X_{\eta}(t;x)\right),\quad$
$\displaystyle\text{for }t\in[0,T],\quad\quad X_{\eta}(0;x)=x,\quad\text{for
}x\in\mathbb{R}^{3},$
$\displaystyle\partial_{t}\rho_{\eta}+\operatorname{div}\left(\rho_{\eta}u_{\eta}\right)$
$\displaystyle=\epsilon\Delta\rho_{\eta}\quad\text{a.e. in }Q$ (125)

and

$\displaystyle-\int_{0}^{T}\int_{\Omega}\rho_{\eta}u_{\eta}\cdot\partial_{t}\phi\
dxdt=$ $\displaystyle\int_{0}^{T}\int_{\Omega}\left(\rho_{\eta}u_{\eta}\otimes
u_{\eta}\right):\mathbb{D}(\phi)+\left(a\rho_{\eta}^{\gamma}+\alpha\rho_{\eta}^{\beta}\right)\operatorname{div}\phi-2\nu\left(\chi_{\eta}\right)\mathbb{D}(u_{\eta}):\mathbb{D}(\phi)$
$\displaystyle-\lambda\left(\chi_{\eta}\right)\operatorname{div}u_{\eta}\operatorname{div}\phi+\rho_{\eta}g\cdot\phi+\frac{1}{\mu}\left(\operatorname{curl}B_{\eta}\times
B_{\eta}\right)\cdot\phi-\epsilon\tilde{z}_{\eta}\cdot\phi$
$\displaystyle-\epsilon\left(\nabla u_{\eta}\nabla\rho_{\eta}\right)\cdot\phi\
dxdt,$ $\displaystyle-\int_{0}^{T}\int_{\Omega}B_{\eta}\cdot\partial_{t}b\
dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}B_{\eta}+u_{\eta}\times
B_{\eta}+\frac{1}{\sigma}J-\frac{\epsilon}{\mu^{2}}z_{\eta}\right]\cdot\operatorname{curl}b$
$\displaystyle-\epsilon\operatorname{curl}\left(\operatorname{curl}B_{\eta}\right):\operatorname{curl}\left(\operatorname{curl}b\right)\
dxdt,$ (126)

where $\chi_{\eta}(t,x):=\textbf{db}_{S_{\eta}(t)}(x)$, for any
$\phi\in\mathcal{D}(Q)$ and any $b\in Y(S_{\eta})$. Further, these functions
satisfy the initial conditions

$\displaystyle\rho_{\eta}(0)=\rho_{0},\quad\quad\left(\rho_{\eta}u_{\eta}\right)(0)=\left(\rho
u\right)_{0},\quad\quad B_{\eta}(0)=B_{0},$

of which the latter identity can be understood in the sense of (35), as well
as the energy inequality

$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{\eta}(t)\left|u_{\eta}(t)\right|^{2}+\frac{a}{\gamma-1}\rho_{\eta}^{\gamma}(t)+\frac{\alpha}{\beta-1}\rho_{\eta}^{\beta}(t)+\frac{1}{2\mu}\left|B_{\eta}(t)\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}2\nu\left(\chi_{\eta}\right)|\mathbb{D}\left(u_{\eta}\right)|^{2}$
$\displaystyle+\lambda\left(\chi_{\eta}\right)\left|\operatorname{div}u_{\eta}\right|^{2}+\alpha\epsilon\beta\rho_{\eta}^{\beta-2}\left|\nabla\rho_{\eta}\right|^{2}+\epsilon\left|\tilde{z}_{\eta}\right|^{\frac{4}{3}}+\frac{1}{\sigma\mu^{2}}\left|\operatorname{curl}B_{\eta}\right|^{2}+\frac{\epsilon}{\mu}\left|\Delta
B_{\eta}\right|^{2}+\frac{\epsilon}{\mu^{3}}\left|z_{\eta}\right|^{\frac{4}{3}}\
dxd\tau$ $\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}+\frac{1}{2\mu}\left|B_{0}\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}\rho_{\eta}g\cdot
u_{\eta}+\frac{1}{\sigma\mu}J\cdot\operatorname{curl}B_{\eta}\ dxd\tau$ (127)

for almost all $t\in[0,T]$ and the estimate

$\displaystyle\epsilon\left\|\partial_{t}\rho_{\eta}\right\|_{L^{\tilde{r}}(Q)}+\epsilon^{2}\left\|\Delta\rho_{\eta}\right\|_{L^{\tilde{r}}(Q)}+\epsilon\left\|\nabla\rho_{\eta}\right\|_{L^{r}(Q)}\leq
c$ (128)

with a constant $c>0$ independent of $\eta$ and $\epsilon$.

From this point on, the remainder of the proof of the main result is straight
forward: In the mechanical part of the problem we can follow precisely the
arguments from [17, Sections 7–9], the additional Lorentz force (c.f. [41])
and regularization term in the momentum equation do not cause any essential
further difficulties. In the induction equation, each limit passage from now
on can be carried out as in the incompressible case in [1] and thus
essentially as in Section 6.2. However, for the convenience of the reader, we
will sketch the main arguments for the remaining three limit passages in the
following sections.

##  7 Limit passage in the penalization method

We continue by passing to the limit with respect to $\eta\rightarrow 0$.
Exactly as in the limit passage with respect to $n\rightarrow\infty$ in
Section 6 we can, due to the energy inequality (127) and the uniform bound
(128), extract suitable subsequences and find functions $z,\tilde{z}\in
L^{\frac{4}{3}}(Q)$ and

$\displaystyle u\in$ $\displaystyle L^{2}(0,T;H_{0}^{1,2}(\Omega)),$ (129)
$\displaystyle 0\leq\rho\in$ $\displaystyle\bigg{\\{}\psi\in
L^{\infty}\left(0,T;L^{\beta}(\Omega)\right)\bigcap
L^{r}\left(0,T;W^{1,r}(\Omega)\right)\bigcap
L^{\tilde{r}}\left(0,T;W^{2,\tilde{r}}(\Omega)\right):$
$\displaystyle\partial_{t}\psi\in L^{\tilde{r}}\left(Q\right),\
\left.\nabla\psi\cdot\operatorname{n}\right|_{\partial\Omega}=0\bigg{\\}}\subset
C\left([0,T];L^{2}(\Omega)\right),$ (130) $\displaystyle B\in$
$\displaystyle\left\\{b\in L^{\infty}\left(0,T;L^{2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{2,2}_{\operatorname{div}}(\Omega)\right):\
\operatorname{curl}b=0\ \text{in }Q^{s}(S),\
b\cdot\text{n}|_{\partial\Omega}=0\right\\},$ (131)

such that

$\displaystyle\rho_{\eta}\rightarrow\rho\ \ \ $ $\displaystyle\text{in
}L^{\beta}\left(Q\right),\quad\quad\quad$ $\displaystyle\rho_{\eta}$
$\displaystyle\rightarrow\rho\quad\ \ \ \text{in
}L^{2}\left(0,T;H^{1,2}(\Omega)\right),$ (132)
$\displaystyle\rho_{\eta}\rightharpoonup\rho\ \ \ $ $\displaystyle\text{in
}L^{\tilde{r}}\left(0,T;W^{2,\tilde{r}}(\Omega)\right),\quad\quad\quad$
$\displaystyle\partial_{t}\rho_{\eta}$
$\displaystyle\rightharpoonup\partial_{t}\rho\ \ \ \ \text{in
}L^{\tilde{r}}\left(Q\right),$ (133) $\displaystyle u_{\eta}\rightharpoonup u\
\ \ $ $\displaystyle\text{in
}L^{2}\left(0,T;H^{1,2}(\Omega)\right),\quad\quad\quad$ $\displaystyle
B_{\eta}$ $\displaystyle\buildrel\ast\over{\rightharpoonup}B\ \ \quad\text{in
}L^{\infty}\left(0,T;L^{2}(\Omega)\right),$ (134) $\displaystyle
B_{\eta}\rightharpoonup B\ \ \ $ $\displaystyle\text{in
}L^{2}\left(0,T;H^{2,2}(\Omega)\right),\quad\quad\quad$ $\displaystyle
z_{\eta}$ $\displaystyle\rightharpoonup z\quad\quad\text{in
}L^{\frac{4}{3}}\left(Q\right),$
$\displaystyle\tilde{z}_{\eta}\rightharpoonup\tilde{z}\ \quad$
$\displaystyle\text{in }L^{\frac{4}{3}}\left(Q\right).$

The boundary conditions of the limit functions in (129)–(131) follow directly
from the corresponding boundary conditions on the $\eta$-level, see
Proposition 6.1. The set-valued function $S$ in (131) is defined by
$S:=S(\cdot):=(X(\cdot;O))^{\delta}$ where $X$, given by

$\displaystyle X_{\eta}\rightarrow X\quad\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right),$

denotes the solution to the initial value problem

$\displaystyle\frac{dX(t;x)}{dt}=R_{\delta}[u]\left(t,X(t;x)\right),\quad
X(0;x)=x\quad\quad\forall x\in\mathbb{R}^{3},$ (135)

c.f. Lemma 10.3. In particular, this lemma also implies that

$\displaystyle\chi_{\eta}\rightarrow\chi:=\textbf{db}_{S(\cdot)}(\cdot)\quad\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right).$

###  7.1 Continuity equation

Making use of the convergences (132)–(134) of $\rho_{\eta}$ and $u_{\eta}$, we
can pass to the limit in the continuity equation (125) and ensure that

$\displaystyle\partial_{t}\rho+\operatorname{div}\left(\rho u\right)=$
$\displaystyle\epsilon\Delta\rho\quad\text{a.e. in }Q,\quad\rho(0)=\rho_{0}.$
(136)

Moreover, this pointwise identity can be renormalized by multiplying it by
$\zeta^{\prime}(\rho_{\epsilon})$ for an arbitrary convex function $\zeta\in
C^{2}([0,+\infty))$. Since $\zeta^{\prime\prime}\geq 0$, this yields

$\displaystyle\partial_{t}\zeta\left(\rho\right)+\operatorname{div}\left(\zeta\left(\rho\right)u\right)+\left[\zeta^{\prime}\left(\rho\right)\rho-\zeta\left(\rho\right)\right]\operatorname{div}u-\epsilon\Delta\zeta\left(\rho\right)=-\epsilon\zeta^{\prime\prime}\left(\rho\right)\left|\nabla\rho\right|^{2}\leq
0$ (137)

almost everywhere in $Q$. This relation will turn out useful in the limit
passage with respect to $\epsilon\rightarrow 0$ in Section 8.

###  7.2 Induction equation

For the limit passage in the induction equation we can argue exactly as in the
limit passage with respect to $n\rightarrow\infty$ in Section 6.2 to show
strong convergence of $B_{\eta}$ in the fluid domain. Hence, we can pass to
the limit in (126) and obtain the identity

$\displaystyle-\int_{0}^{T}\int_{\Omega}B\cdot\partial_{t}b\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}B+u\times
B-\frac{\epsilon}{\mu^{2}}z+\frac{1}{\sigma}J\right]\cdot\operatorname{curl}b$
$\displaystyle-\epsilon\operatorname{curl}\left(\operatorname{curl}B\right)\cdot\operatorname{curl}\left(\operatorname{curl}b\right)\
dxdt\quad\quad\forall b\in Y(S).$ (138)

Moreover, the initial condition $B(0)=B_{0}$ also follows as in Section 6.2.

###  7.3 Momentum equation and compatibility of the velocity field

From the uniform bounds given by the energy inequality (127) we further infer
the existence of $z_{5}\in L^{\frac{6}{5}}(Q)$ such that, for a chosen
subsequence,

$\displaystyle\rho_{\eta}u_{\eta}\otimes\eta\rightharpoonup z_{5}\quad\text{in
}L^{\frac{6}{5}}(Q).$

For the limit passage in the momentum equation we need to identify $z_{5}$:
The choice of test functions $\phi\in T(S)$, which satisfy
$\mathbb{D}(\phi)=0$ in a neighbourhood of $\overline{Q}^{s}(S)$, allows us to
control the variable viscosity coefficients $\nu(\chi_{\eta})$ and
$\lambda(\chi_{\eta})$ in the momentum equation (126), since these remain
bounded in the fluid region according to their definition in (52), (53). This
enables us to deduce strong convergence of the momentum function
$\rho_{\eta}u_{\eta}$ in the fluid domain similarly to the strong convergence
(119) of the magnetic induction in the Galerkin limit. Indeed, we fix an
arbitrary interval $I\subset(0,T)$ and an arbitrary ball $U\subset\Omega$ such
that $\overline{I\times U}\subset Q^{f}(S)$ and deduce from the momentum
equation, for any $\Phi\in\mathcal{D}(U)$, the dual estimate

$\displaystyle\left\|\partial_{t}\int_{U}\rho_{\eta}u_{\eta}\cdot\Phi
dx\right\|_{L^{\min\left(\frac{6}{5},\frac{5\beta-3}{4\beta}\right)}(I)}\leq
c$

for a constant $c>0$ depending on $\Phi$ but not on $\eta$. From the Arzelá-
Ascoli theorem it follows that

$\displaystyle\rho_{\eta}u_{\eta}\rightarrow\rho u\quad\text{in
}C_{\text{weak}}\left(\overline{I};L^{\frac{2\beta}{\beta+1}}(U)\right)\quad\text{and
hence in }L^{2}\left(I;H^{-1,2}(U)\right),$

which implies that

$\displaystyle z_{5}:\mathbb{D}(\phi)=(\rho u\otimes
u):\mathbb{D}(\phi)\quad\text{a.e. in }Q.$

for all test functions $\phi\in T(S)$. Letting $\eta$ tend to $0$ in (126) we
thus obtain

$\displaystyle-\int_{0}^{T}\int_{\Omega}\rho u\cdot\partial_{t}\phi\ dxdt=$
$\displaystyle\int_{0}^{T}\int_{\Omega}\left(\rho u\otimes
u\right):\mathbb{D}(\phi)+\left(a\rho^{\gamma}+\alpha\rho^{\beta}\right)\operatorname{div}\phi-2\nu\mathbb{D}\left(u\right):\mathbb{D}(\phi)-\lambda\operatorname{div}u\operatorname{div}\phi$
$\displaystyle+\rho g\cdot\phi+\frac{1}{\mu}\left(\text{curl}B\times
B\right)\cdot\phi-\epsilon\tilde{z}\cdot\phi-\epsilon\left(\nabla
u\nabla\rho\right)\cdot\phi\ dxdt\quad\quad\forall\phi\in T(S).$ (139)

Moreover, since $\nu(\chi_{\eta})$ and $\lambda(\chi_{\eta})$ blow up in the
solid part of the domain, the energy inequality (127) shows that

$\displaystyle\mathbb{D}(u)=0\quad\text{a.e. in }Q^{s}(S).$

Hence, there are rigid velocity fields $u^{s^{i}}$ which coincide with $u$
almost everywhere in the $\delta$-neighbourhoods
$S^{i}(t):=(X(t;O^{i}))^{\delta}$ of the sets $X(t;O^{i})$. Consequently, due
to the property (58) of the regularized velocity field $R_{\delta}[u]$, we can
replace $R_{\delta}[u]$ in the initial value problem (135) by $u^{s^{i}}$ for
$x\in O^{i}$. The combination of the latter two conditions at first yields
compatibility (c.f. (19), (20)) of $u$ with the system
$\\{O^{i},X^{i}\\}_{i=1}^{m}$, where each $X^{i}(t)$ denotes an isometry which
coincides with $X(t)$ in $O^{i}$. However, the fact that each $X^{i}(t)$ is an
isometry implies that

$\displaystyle
S^{i}(t)=\left(X^{i}\left(t;O^{i}\right)\right)^{\delta}=X^{i}\left(t;S^{i}_{0}\right)\quad\quad\text{and
thus }\quad\quad u(t,\cdot)=u^{s^{i}}(t,\cdot)\quad\text{a.e. in
}X^{i}\left(t;S^{i}\right),$

Consequently, we infer that $u$ is even compatible with
$\\{S_{0}^{i},X^{i}\\}_{i=1}^{m}$. Finally, the initial condition $(\rho
u)(0)=(\rho u)_{0}$, in the sense of (35), follows by similar arguments as the
initial condition for the magnetic induction in Section 6.2.

###  7.4 Energy inequality

We drop, among other non-negative terms, the variable parts of the viscosity
coefficients on the left-hand side of the energy inequality (127) and pass to
the limit to see that

$\displaystyle\int_{\Omega}\frac{1}{2}\rho(t)\left|u(t)\right|^{2}+\frac{a}{\gamma-1}\rho^{\gamma}(t)+\frac{\alpha}{\beta-1}\rho^{\beta}(t)+\frac{1}{2\mu}\left|B(t)\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}2\nu|\mathbb{D}(u)|^{2}+\lambda\left|\operatorname{div}u\right|^{2}+\epsilon\left|\tilde{z}\right|^{\frac{4}{3}}$
$\displaystyle+\frac{1}{\sigma\mu^{2}}\left|\operatorname{curl}B\right|^{2}+\frac{\epsilon}{\mu}\left|\Delta
B\right|^{2}+\frac{\epsilon}{\mu^{3}}\left|z\right|^{\frac{4}{3}}\ dxd\tau$
$\displaystyle\leq$
$\displaystyle\int_{\Omega}\frac{1}{2}\rho_{0}\left|u_{0}\right|^{2}+\frac{a}{\gamma-1}\rho_{0}^{\gamma}+\frac{\alpha}{\beta-1}\rho_{0}^{\beta}+\frac{1}{2\mu}\left|B_{0}\right|^{2}\
dx+\int_{0}^{t}\int_{\Omega}\rho g\cdot
u+\frac{1}{\sigma\mu}J\cdot\operatorname{curl}B\ dxd\tau$ (140)

for almost all $t\in[0,T]$.

##  8 Limit passage in the regularization terms

The next step is the limit passage with respect to $\epsilon\rightarrow 0$.
The energy inequality (140) yields the bounds needed for Corollary 10.1, which
implies the existence of isometries
$X^{i}(t):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ such that

$\displaystyle X^{i}_{\epsilon}\rightarrow X^{i}\quad$ $\displaystyle\text{in
}C\left([0,T];C_{\text{loc}}\left(\mathbb{R}^{3}\right)\right).$

We write $X:[0,T]\times S_{0}\rightarrow\mathbb{R}^{3}$,
$X(t)|_{S_{0}^{i}}:=X^{i}(t)$ and $S=S(\cdot):=X(\cdot;S_{0})$. Further,
testing the continuity equation (136) by $\rho_{\epsilon}$, we see that

$\displaystyle\epsilon^{\frac{1}{2}}\left\|\nabla\rho_{\epsilon}\right\|_{L^{2}(Q)}\leq
c$

for a constant $c>0$ independent of $\epsilon$. This, together with the energy
inequality (140), yields the existence of functions $0\leq\rho\in
L^{\infty}(0,T;L^{\beta}(\Omega))$ and

$\displaystyle u\in$ $\displaystyle\left\\{\phi\in
L^{2}\left(0,T;H_{0}^{1,2}\left(\Omega\right)\right):\ \mathbb{D}(\phi)=0\
\text{in }Q^{s}(S)\right\\},$ (141) $\displaystyle B\in$
$\displaystyle\left\\{b\in L^{\infty}\left(0,T;L^{2}(\Omega)\right)\bigcap
L^{2}\left(0,T;H^{1,2}_{\operatorname{div}}(\Omega)\right):\
\operatorname{curl}b=0\ \text{in }Q^{s}(S),\
b\cdot\text{n}|_{\partial\Omega}=0\right\\}$ (142)

such that for certain extracted subsequences it holds

$\displaystyle\rho_{\epsilon}\buildrel\ast\over{\rightharpoonup}\rho\ \ \ $
$\displaystyle\text{in
}L^{\infty}\left(0,T;L^{\beta}(\Omega)\right),\quad\quad\quad$ $\displaystyle
u_{\epsilon}$ $\displaystyle\rightharpoonup u\ \ \ \ \text{in
}L^{2}\left(0,T;H^{1,2}(\Omega)\right),$ $\displaystyle
B_{\epsilon}\buildrel\ast\over{\rightharpoonup}B\ \ \ $ $\displaystyle\text{in
}L^{\infty}\left(0,T;L^{2}(\Omega)\right),\ \ \ \ \ \ $ $\displaystyle
B_{\epsilon}$ $\displaystyle\rightharpoonup B\ \ \ \text{in
}L^{2}\left(0,T;H^{1,2}(\Omega)\right),$
$\displaystyle\epsilon\nabla\rho_{\epsilon},\ \epsilon\Delta
B_{\epsilon}\rightarrow 0\quad$ $\displaystyle\text{in
}L^{2}\left(Q\right),\quad\quad$ $\displaystyle\epsilon z_{\epsilon},\
\epsilon\tilde{z}_{\epsilon}$ $\displaystyle\rightarrow 0\quad\ \text{in
}L^{\frac{4}{3}}\left(Q\right).$ (143)

The boundary conditions of the limit functions in (141) and (142) follow
directly from the corresponding boundary conditions of the velocity field and
the magnetic induction in (129) and $\eqref{-290}$ on the $\epsilon$-level.
Moreover, the velocity field $u$ is compatible with the system
$\\{S^{i}_{0},X^{i}\\}_{i=1}^{m}$, c.f. Corollary 10.1.

###  8.1 Continuity equation

Similarly to the strong convergence (119) of the magnetic induction, we deduce
from the continuity equation (136) that $\rho_{\epsilon}$ even converges to
$\rho$ in $C_{\text{weak}}([0,T];L^{\beta}(\Omega))$. This, together with the
vanishing artificial viscosity term, c.f. (143), is sufficient to pass to the
limit in (136) and to obtain

$\displaystyle\partial_{t}\rho+\operatorname{div}\left(\rho
u\right)=0\quad\text{in }\mathcal{D}^{\prime}\left(Q\right).$ (144)

In fact, as at this stage of the approximation it holds $\rho\in L^{2}(Q)$,
$u\in L^{2}(0,T;H_{0}^{1,2}(\Omega))$, we can use the regularization procedure
by DiPerna and Lions, c.f. [37, Lemma 6.8, Lemma 6.9], to see that $\rho$ and
$u$, extended by $0$ outside of $\Omega$, even satisfy the renormalized
continuity equation (30), (31). This in turn implies that $\rho\in
C([0,T];L^{1}(\Omega))$, c.f. [37, Lemma 6.15], and $\rho$ satisfies the
initial condition $\rho(0)=\rho_{0}$.

###  8.2 Induction equation

The regularization terms in the induction equation vanish as $\epsilon$ tends
to $0$ according to the convergences (143). Apart from that we can argue
exactly as in the Galerkin limit in Section 6.2 to pass to the limit with
respect to $\epsilon\rightarrow 0$ in (138) and to infer that

$\displaystyle-\int_{0}^{T}\int_{\Omega}B\cdot\partial_{t}b\
dxdt=\int_{0}^{T}\int_{\Omega}\left[-\frac{1}{\sigma\mu}\operatorname{curl}B+u\times
B+\frac{1}{\sigma}J\right]\cdot\operatorname{curl}b\ dxdt\quad\quad\forall
b\in Y(S).$ (145)

###  8.3 Momentum equation

In order to pass to the limit in the pressure terms, we first consider an
arbitrary compact set $K\subset Q^{f}(S)$. Denoting by $\mathcal{B}_{\Omega}$
the Bogovskii operator in $\Omega$ (c.f. [37, Section 3.3.1.2]), we test the
momentum equation (139) by

$\displaystyle\phi_{\epsilon}(t,x):=\Phi(t,x)\mathcal{B}_{\Omega}\left[\rho_{\epsilon}(t,\cdot)-\frac{1}{|\Omega|}\int_{\Omega}\rho_{\epsilon}(t,y)\
dy\right](t,x),$ (146)

where $\Phi\in\mathcal{D}(Q^{f}(S))$ is a cut-off function equal to $1$ in
$K$. This procedure leads to a bound of $\rho_{\epsilon}$ in $L^{\beta+1}(K)$
uniformly in $\epsilon$, c.f. [17, Lemma 8.1] and the references therein.
These bounds in turn allow us to find $z_{6}\in
L^{\frac{\gamma+1}{\gamma}}(K)$, $z_{7}\in L^{\frac{\beta+1}{\beta}}(K)$ such
that

$\displaystyle\rho_{\epsilon}^{\gamma}\rightharpoonup z_{6}\quad\text{in