# Existence of solutions for an interaction problem between a bubble and a compressible viscous fluid

Fabien Lespagnol, Matthieu Hillairet

<figure><img src="image_1.png" /></figure>

## To cite this version:

Fabien Lespagnol, Matthieu Hillairet. Existence of solutions for an interaction problem between a bubble and a compressible viscous fluid. 2026. hal-05514263

HAL Id: hal-05514263

https://hal.science/hal-05514263v1

Preprint submitted on 17 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Existence of solutions for an interaction problem between a bubble and a compressible viscous fluid

FABIEN LESPAGNOL¹ AND MATTHIEU HILLAIRET¹

¹ IMAG, ANGUS, Univ Montpellier, Inria, CNRS, Montpellier, France

## Abstract

In this paper, we study the dynamics of a finite number of spherical bubbles in a compressible fluid within a bounded open domain of $\mathbb{R}^3$. The fluid-bubble interaction is described by a system of nonlinear partial differential equations (PDEs) and ordinary differential equations (ODEs) coupling the fluid's density, velocity and pressure to the bubble's translational, rotational and radial velocities. We prove the existence of weak solutions for this model until the collision or collapse of the bubbles. The formulation of the fluid-bubble system, along with the techniques used for the existence proof, is inspired by penalization methods developed for fluid-solid interaction. The main contribution of this work is the addition of a radial expansion-contraction mode in the bubble motion, which introduces new nonlinear terms in the momentum equations that need to be treated carefully in the compactness arguments.

## 1 Introduction

Let $m$ be the number of bubbles immersed in the fluid flow, for each $i \in \{1, \dots, m\}$ and time $t \in \mathbb{R}_+$, we denote by $B_i(t)$ and $\mathcal{F}(t)$ the time-dependent domains occupied by bubble $i$ and the fluid, respectively. We assume that the bubbles move inside a fixed, bounded, smooth domain $\Omega \subset \mathbb{R}^3$ and

$$ \mathcal{F}(t) = \Omega \setminus \bigcup_{i=1}^{m} \overline{B_i(t)}. $$

The fluid flow is assumed to satisfy the isentropic compressible Navier-Stokes equations. We denote by $\rho_f$, $u_f$ and $p_f$ its density, velocity and pressure, respectively. The system reads

$$ \partial_t \rho_f + \text{div}(\rho_f u_f) = 0, \quad t \in \mathbb{R}_+, \quad x \in \mathcal{F}(t), \qquad (1.1) $$

$$ \partial_t (\rho_f u_f) + \text{div}(\rho_f u_f \otimes u_f) = \text{div}(\mathbb{T}_f(u_f, p_f)) - \rho_f g, \quad t \in \mathbb{R}_+, \quad x \in \mathcal{F}(t), \qquad (1.2) $$

with $-\rho_f g$ the gravitational force. We assume the fluid is Newtonian with stress tensor $\mathbb{T}_f(u_f, p_f)$ given by the constitutive law

$$ \mathbb{T}_f(u_f, p_f) \stackrel{\text{def}}{=} 2\mu_f (\mathbb{D}(u_f) - \frac{1}{3} \text{div}(u_f) \mathbb{I}_3) + (\nu_f \text{div}(u_f) - p_f) \mathbb{I}_3, $$

where

$$ \mathbb{D}(u_f) \stackrel{\text{def}}{=} \frac{1}{2} (\nabla u_f + (\nabla u_f)^\mathrm{T}) $$

denotes the symmetric part of the velocity gradient, $\mathbb{I}_3$ the identity matrix in $\mathbb{R}^{3 \times 3}$, and $\mu_f > 0$, $\nu_f \ge 0$ the shear and bulk viscosity coefficients. We further assume that the flow is in the barotropic regime and the relation between $p_f$ and $\rho_f$ is given by the constitutive law

$$ p_f \stackrel{\text{def}}{=} p_f(\rho_f) = a_f \rho_f^{\gamma_f} \quad \text{with } a_f > 0, \gamma_f > \frac{3}{2}. \qquad (1.3) $$

We now formulate the equations describing the kinematics and dynamics of the bubbles. This model, introduced in [24, 25], is based on the assumptions that the bubbles are spherical and behave as a compressible viscous fluid with infinite shear viscosity. In particular, denoting by $\mathcal{B}(x, R)$ the ball of center $x \in \mathbb{R}^3$ and radius $R \in \mathbb{R}_+$, for each $i \in \{1, \dots, m\}$ and time $t \in \mathbb{R}_+$, there exists $R_i(t) \in \mathbb{R}_+$ and $x_i(t) \in \mathbb{R}^3$ such that

$$ \mathcal{B}_i(t) = \mathcal{B}(x_i(t), R_i(t)). $$

Setting $\mathcal{B}_{i,0} \stackrel{\text{def}}{=} \mathcal{B}_i(0)$, $R_{i,0} \stackrel{\text{def}}{=} R_i(0)$, $x_{i,0} \stackrel{\text{def}}{=} x_i(0)$, the motion of the bubble $\mathcal{B}_i$ is then given by

$$ \mathcal{B}_i(t) = \eta_i[t](\mathcal{B}_{i,0}) \quad \forall t \in \mathbb{R}_+, $$
where

$$ \eta_i[t] : \begin{cases} \mathbb{R}^3 \to \mathbb{R}^3 \\ \hat{x} \mapsto x_i(t) + \frac{R_i(t)}{R_{i,0}} \mathcal{O}_i(t) (\hat{x} - x_{i,0}) \end{cases} \qquad (1.4) $$

and $\mathcal{O}_i : \mathbb{R}_+ \to SO(3)$ is a rotation matrix satisfying $\mathcal{O}_i(0) = \mathbb{I}_3$. The Lagrangian velocity of bubble $i$ is then given for $(t, \hat{x}) \in \mathbb{R}_+ \times \mathcal{B}_{i,0}$ by

$$ \hat{u}_i(t, \hat{x}) \stackrel{\text{def}}{=} \partial_t \eta_i[t](\hat{x}) = \dot{x}_i(t) + \frac{\dot{R}_i(t)}{R_{i,0}} \mathcal{O}_i(t)(\hat{x} - x_{i,0}) + \frac{R_i(t)}{R_{i,0}} \dot{\mathcal{O}}_i(t)(\hat{x} - x_{i,0}). \qquad (1.5) $$

From (1.5), we can deduce the following expression for the Eulerian velocity field $u_i$ of bubble $i$:

$$ u_i(t, x) \stackrel{\text{def}}{=} (\partial_t \eta_i[t]) ((\eta_i[t])^{-1}(x)) = V_i(t) + \omega_i(t) \times (x - x_i(t)) + \frac{1}{3} \Lambda_i(t) (x - x_i(t)), \quad t \in \mathbb{R}_+, \quad x \in \mathcal{B}_i(t), \qquad (1.6) $$

where

$$ V_i \stackrel{\text{def}}{=} \dot{x}_i, \quad \frac{1}{3} \Lambda_i \stackrel{\text{def}}{=} \frac{\dot{R}_i}{R_i} \qquad (1.7) $$

and $\omega_i$ is the unique vector-valued function satisfying

$$ \dot{\mathcal{O}}_i \mathcal{O}_i^T x = \omega_i \times x \quad \forall x \in \mathbb{R}^3. \qquad (1.8) $$

The quantities $V_i$, $\omega_i$ and $\Lambda_i$ then correspond to the translational, rotational and radial components of the velocity of bubble $i$, respectively. Regarding the bubble dynamics, we assume that the density of each bubble $\rho_i$ is uniform in $\mathcal{B}_i$ so that $\rho_i$ is related to the bubble radius $R_i$ and mass $m_i$ by

$$ \rho_i = \left( \frac{R_{i,0}}{R_i} \right)^3 = \frac{3m_i}{4\pi R_i^3}. \qquad (1.9) $$

Moreover, the center of mass of each bubble coincides with its geometric center, that is

$$ \frac{1}{m_i} \int_{\mathcal{B}_i} \rho_i y dy = \frac{3}{4\pi R_i^3} \int_{\mathcal{B}_i} y dy = x_i. \qquad (1.10) $$

We also introduce the inertia tensor of the bubble $i$ defined by

$$ \mathcal{J}_i \stackrel{\text{def}}{=} \int_{\mathcal{B}_i} \rho_i (|y - x_i|^2 \mathbb{I}_3 - (y - x_i) \otimes (y - x_i)) dy = \frac{2}{5} m_i R_i^2 \mathbb{I}_3, \qquad (1.11) $$

and a mechanical parameter $K_i : \mathbb{R}_+ \to \mathbb{R}_+$ analogous to $\mathcal{J}_i$, quantifying the contribution of radius variations to kinetic energy, defined by

$$ K_i \stackrel{\text{def}}{=} \frac{1}{9} \int_{\mathcal{B}_i} \rho_i |y - x_i|^2 dy = \frac{m_i}{15} R_i^2. \qquad (1.12) $$

The dynamics equations of the bubble $i \in \{1, \dots, m\}$ read

$$ \frac{\mathrm{d}}{\mathrm{d}t}(m_i V_i) = - \int_{\partial \mathcal{B}_i} \mathbb{T}_f(u_f, p_f) \mathrm{nds} - \int_{\mathcal{B}_i} \rho_i g dy, \qquad (1.13) $$

$$ \frac{\mathrm{d}}{\mathrm{d}t}(\mathcal{J}_i \omega_i) = - \int_{\partial \mathcal{B}_i} (s - x_i) \times \mathbb{T}_f(u_f, p_f) \mathrm{nds}, \qquad (1.14) $$

$$ \begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t}(K_i \Lambda_i) - \frac{1}{3} ((\mathcal{J}_i \omega_i) \cdot \omega_i + K_i |\Lambda_i|^2) = & \\ - \frac{1}{3} \int_{\partial \mathcal{B}_i} (s - x_i) \cdot \mathbb{T}_f(u_f, p_f) \mathrm{nds} + (\nu_i \Lambda_i - p_i + \frac{\kappa_i}{R_i}) |\mathcal{B}_i|, & \end{aligned} \qquad (1.15) $$

where $\nu_i$ is the bulk viscosity coefficient, $\kappa_i$ a surface-tension coefficient and $p_i$ is the bubble pressure satisfying

$$ p_i \stackrel{\text{def}}{=} p_i(\rho_i) = a_i \rho_i^{\gamma_i} \text{ with } a_i > 0, \gamma_i > \frac{3}{2}. \qquad (1.16) $$

To close the system, we impose the no-slip condition at each fluid-bubble interface and homogeneous Dirichlet boundary conditions on $\partial \Omega$:

$$ u_f - u_i = 0 \text{ on } \partial \mathcal{B}_i, \qquad (1.17) $$

$$ u_f = 0 \text{ on } \partial \Omega, \qquad (1.18) $$
and prescribe the initial data:

$$
\begin{align}
\rho_f(0, x) &= \rho_{f,0}(x), \quad (\rho_f u_f)(0, x) = q_{f,0}(x), \quad x \in \mathcal{F}(0), \nonumber \\
\rho_i(0, x) &= \rho_{i,0}, \quad u_i(0, x) = V_{i,0} + \omega_{i,0} \times (x - x_{i,0}) + \frac{\Lambda_{i,0}}{3}(x - x_{i,0}), \quad x \in \mathcal{B}_i(0). \tag{1.19}
\end{align}
$$

**Remark 1.1.** Note that by substituting equations (1.7)-(1.8) and (1.9)-(1.12) into (1.13)-(1.15), we obtain the dynamical equa-
tions governing the evolution of the bubble $\mathcal{B}_i$ in terms of its center of mass $x_i$, its radius $R_i$ and the rotational velocity $\omega_i$:

$$
\begin{align*}
\frac{d}{dt} (m_i \dot{x}_i) &= - \int_{\partial \mathcal{B}_i} \mathbb{T}_f (u_f, p_f) \mathrm{nds} - m_i g, \\
\frac{d}{dt} \left( \frac{2m_i}{5} R_i^2 \omega_i \right) &= - \int_{\partial \mathcal{B}_i} (s - x_i) \times \mathbb{T}_f (u_f, p_f) \mathrm{nds}, \\
\frac{d}{dt} \left( \frac{m_i}{5} \dot{R}_i R_i \right) &- \frac{m_i R_i^2}{3} \left( \frac{2}{5} |\omega_i|^2 + \frac{3}{5} \left| \frac{\dot{R}_i}{R_i} \right|^2 \right) = \\
&\quad - \frac{1}{3} \int_{\partial \mathcal{B}_i} (s - x_i) \cdot \mathbb{T}_f (u_f, p_f) \mathrm{nds} - \left( 3 \nu_i \frac{\dot{R}_i}{R_i} - p_i \left( \frac{3m_i}{4\pi R_i^3} \right) + \frac{\kappa_i}{R_i} \right) \frac{4\pi R_i^3}{3}.
\end{align*}
$$

**Remark 1.2.** The second term in the left-hand side of (1.15) arises from the time-dependence of $\mathcal{J}_i$ and $K_i$, which itself results from the evolution of $R_i$. These terms ensure that the left-hand side conserves the kinetic energy. In particular, these additional terms allow to recover formally the following energy equation for the bubble $i$:

$$
\frac{d}{dt} \left( \int_{\mathcal{B}_i} \rho_i \frac{|u_i|^2}{2} + a_i \frac{\rho_i^{\gamma_i}}{\gamma_i - 1} dy \right) + \int_{\mathcal{B}_i} \nu_i |\Lambda_i|^2 dy + \frac{2\pi}{3} \kappa_i \frac{d}{dt} |R_i|^2 = - \int_{\partial \mathcal{B}_i} \mathcal{T}_f(u_f, p_f) n \cdot u_i ds.
$$

**1.1 Weak formulation**

In this paragraph, we formally derive a weak formulation for the system (1.1)-(1.19). We set

$$
\rho \stackrel{\text{def}}{=} \rho_f \mathbb{1}_\mathcal{F} + \sum_{i=1}^{m} \rho_i \mathbb{1}_{\mathcal{B}_i}, \quad u \stackrel{\text{def}}{=} u_f \mathbb{1}_\mathcal{F} + \sum_{i=1}^{m} u_i \mathbb{1}_{\mathcal{B}_i}. \qquad (1.20)
$$

Let $T > 0$, from (1.1), (1.9) and (1.20), we obtain for the continuity equation

$$
\int_0^T \int_\Omega \rho \partial_t \varphi + (\rho u) \cdot \nabla \varphi dy d\tau = 0 \quad (1.21)
$$

for all $\varphi \in \mathcal{D}((0,T) \times \Omega)$. To derive a weak formulation of the momentum equation for the fluid and the bubble dynamics, we first define the bubble and fluid regions as follows:

$$
\overline{Q}_i \stackrel{\text{def}}{=} \{(t, x) \mid t \in [0, T], x \in \overline{\mathcal{B}_i(t)}\}, \quad \overline{Q}_b \stackrel{\text{def}}{=} \bigcup_{i=1}^{m} \overline{Q}_i, \quad Q_f \stackrel{\text{def}}{=} ((0, T) \times \Omega) \setminus \overline{Q}_b.
$$

We also introduce the space of motion and velocity fields corresponding to the kinematics of the bubbles:

$$
\mathcal{M} = \{ \phi_b \mid \exists a \in \mathbb{R}^3, R \in \mathbb{R}, \O \in \mathcal{SO}(3), \phi_b(x) = a + R \O x \quad \forall x \in \mathbb{R}^3 \}.
$$

$$
S = \left\{ \varphi_b \mid \exists V, \omega \in \mathbb{R}^3, \Lambda \in \mathbb{R}, \varphi_b(x) = V + \omega \times x + \frac{\Lambda}{3} x \quad \forall x \in \mathbb{R}^3 \right\}.
$$

Note that the spaces $\mathcal{M}$ and $\mathcal{S}$ inherit the canonical topology of finite-dimensional real vector spaces. Finally, we set the following test functions space $\mathcal{T}(Q_b)$ over $(0, T) \times \Omega$:

$$
\mathcal{T}(Q_b) \stackrel{\text{def}}{=} \left\{ \varphi \in \mathcal{D}((0, T) \times \Omega) \text{ such that there exists } \{\varphi_i\}_{i=1}^m \subset \mathcal{D}(0, T; S) \\
\text{ satisfying for } i \in \{1, \dots, m\}, \varphi = \varphi_i \text{ in an open neighborhood of } \overline{Q_i} \right\}.
$$
Observe that the space of test functions $\mathcal{T}(Q_b)$ depends on the solution through the domains $\mathcal{F}(t)$ and $\mathcal{B}_i(t)$. Combining (1.2) with (1.13)-(1.15) and (1.20), we obtain the following weak formulation for the momentum equation:

$$ 
\begin{aligned}
& \int_0^T \int_{\Omega} (\rho u) \cdot \partial_t \varphi + (\rho u \otimes u) : \mathcal{D}(\varphi) + p(\rho, \{1_{\mathcal{B}_i}\}_{i=1}^m) \mathrm{div} \varphi dyd\tau \\
& = \int_0^T \int_{\Omega} 2\mu \left( \mathcal{D}(u) - \frac{1}{3} \mathrm{div}(u) \mathbb{I}_3 \right) : \left( \mathcal{D}(\varphi) - \frac{1}{3} \mathrm{div}(\varphi) \mathbb{I}_3 \right) dyd\tau + \int_0^T \int_{\Omega} \nu \mathrm{div}(u) \mathrm{div}(\varphi) + \tilde{\kappa}_b \mathrm{div} \varphi - \rho g \cdot \varphi dyd\tau \quad (1.22)
\end{aligned}
$$

for all $\varphi \in \mathcal{T}(Q_b)$, where

$$ 
\begin{aligned}
\mu &\stackrel{\mathrm{def}}{=} \mu_f 1_{\mathcal{F}}, \quad \nu \stackrel{\mathrm{def}}{=} \nu_f 1_{\mathcal{F}} + \sum_{i=1}^m \nu_i 1_{\mathcal{B}_i}, \quad \tilde{\kappa}_b \stackrel{\mathrm{def}}{=} \sum_{i=1}^m \frac{\kappa_i}{R_i} 1_{\mathcal{B}_i}, \\
p(\rho, \{1_{\mathcal{B}_i}\}_{i=1}^m) &\stackrel{\mathrm{def}}{=} a_f \rho^{\gamma_f} \left( 1 - \sum_{i=1}^m 1_{\mathcal{B}_i} \right) + \sum_{i=1}^m a_i \rho^{\gamma_i} 1_{\mathcal{B}_i}.
\end{aligned}
$$

**Remark 1.3.** For all $V, \omega, a, b \in \mathbb{R}^3$, $\Lambda \in \mathbb{R}$ and $x \in \mathbb{R}^3$,

$$ V + \omega \times (x-a) + \frac{\Lambda}{3}(x-a) = \tilde{V} + \tilde{\omega} \times (x-b) + \frac{\tilde{\Lambda}}{3}(x-b) $$

where

$$ \tilde{V} \stackrel{\mathrm{def}}{=} V - \omega \times (a-b) - \frac{\Lambda}{3}(a-b), \quad \tilde{\omega} \stackrel{\mathrm{def}}{=} \omega, \quad \tilde{\Lambda} \stackrel{\mathrm{def}}{=} \Lambda. $$

Consequently, for all $a \in \mathbb{R}^3$, the space $S$ can equivalently be defined by

$$ S = \left\{ \varphi_b \mid \exists V, \omega \in \mathbb{R}^3, \Lambda \in \mathbb{R}, \quad \varphi_b(x) = V + \omega \times (x-a) + \frac{\Lambda}{3}(x-a) \quad \forall x \in \mathbb{R}^3 \right\}. $$

In particular, any test function $\varphi \in \mathcal{T}(Q_b)$ satisfies

$$ \varphi(t, x) = \tilde{V}_i(t) + \tilde{\omega}_i(t) \times (x - x_i(t)) + \frac{\tilde{\Lambda}_i(t)}{3}(x - x_i(t)) \quad \text{on a neighbourhood of } \overline{Q}_i, $$

with $\tilde{V}_i, \tilde{\omega}_i, \tilde{\Lambda}_i \in \mathcal{D}(0, T)$.

**Definition 1.1.** Let $\Omega \subset \mathbb{R}^3$ be a bounded smooth domain. For each $i \in \{1, \dots, m\}$, let $\mathcal{B}_{i,0} \subset \Omega$ be a ball of radius $R_{i,0} \in (0, \infty)$ and center $x_{i,0} \in \mathbb{R}^3$. Let $\rho_{i,0} \in (0, \infty)$, $V_{i,0}, \omega_{i,0} \in \mathbb{R}^3$ and $\Lambda_{i,0} \in \mathbb{R}$. We prescribe the following initial conditions:

$$ \rho_0 \in L^{\gamma_f}(\Omega), \quad \rho_0 \ge 0 \text{ a.e in } \Omega, \quad \rho_0 = \rho_{i,0} \text{ a.e. in } \mathcal{B}_{i,0}, \quad (1.23) $$

$$ q_0 \in L^{\frac{2\gamma_f}{\gamma_f+1}}(\Omega), \quad q_0 1_{\rho_0=0} = 0 \text{ a.e. in } \Omega, \quad \frac{|q_0|^2}{\rho_0} 1_{\rho_0>0} \in L^1(\Omega), \quad (1.24) $$

$$ u_0 = V_{i,0} + \omega_{i,0} \times (x - x_{i,0}) + \frac{\Lambda_{i,0}}{3}(x - x_{i,0}) \text{ a.e. in } \mathcal{B}_{i,0}, \quad q_0 = \rho_0 u_0 \text{ a.e. in } \mathcal{B}_{i,0}. \quad (1.25) $$

A triplet $(\{\mathcal{B}_i\}_{i=1}^m, \rho, u)$ is a bounded energy weak solution of the system (1.1)-(1.2), (1.13)-(1.15), (1.17)-(1.19), associated with initial conditions (1.23)-(1.25) if the following holds:

1. For all $i \in \{1, \dots, m\}$, there exists $\eta_i \in H^1([0, T]; \mathcal{M})$ such that for all $t \in [0, T]$, $\mathcal{B}_i(t)$ is a ball of radius $R_i(t) \in (0, \infty)$ and center $x_i(t) \in \mathbb{R}^3$ given by

$$ \mathcal{B}_i(t) = \eta_i[t](\mathcal{B}_{i,0}). $$

2. The velocity $u$ belongs to $L^2(0, T; H^1(\Omega))$ and satisfies

$$ u(t, x) = (\partial_t \eta_i[t]) ((\eta_i[t])^{-1}(x)) \quad \text{for a.e. } (t, x) \in Q_i. \quad (1.26) $$

3. The density $\rho$ belongs to $L^\infty(0, T; L^{\gamma_f}(\Omega))$ and satisfies

$$ \rho(t, x) = (R_{i,0}/R_i(t))^3 \rho_{i,0} \quad \text{for a.e. } (t, x) \in Q_i. \quad (1.27) $$
4. Provided that they are extended by 0 in $\mathbb{R}^3 \setminus \Omega$, $\rho$ and $u$ satisfy the continuity equation (1.21) for all $\varphi \in D((0,T) \times \mathbb{R}^3)$ together with its renormalized form

$$ \partial_t b(\rho) + \operatorname{div}(b(\rho)u) + (b'(\rho)\rho - b(\rho))\operatorname{div}(u) = 0 \text{ in } \mathcal{D}'((0,T) \times \mathbb{R}^3), \quad (1.28) $$

for any $b \in C([0, \infty)) \cap C^1((0, \infty))$ satisfying

$$ |b'(z)| \leq cz^{-\lambda_0}, \quad z \in (0,1], \quad |b'(z)| \leq cz^{\lambda_1}, \quad z \geq 1, \quad (1.29) $$

with $\lambda_0 < 1$, $-1 < \lambda_1 \leq \frac{s(\gamma_f)}{2} - 1$, $s(\gamma_f) = \frac{5\gamma_f - 3}{3}$.

5. The weak form of the momentum equation (1.22) holds for all $\varphi \in \mathcal{T}(Q_b)$.

6. The initial conditions are satisfied in the following sense:

$$ \lim_{t \to 0^+} \int_{\Omega} \rho(t)\varphi d\tau = \int_{\Omega} \rho_0\varphi d\tau \quad \forall \varphi \in \mathcal{D}(\Omega), \quad (1.30) $$

$$ \lim_{t \to 0^+} \int_{\Omega} (\rho u)(t) \varphi d\tau = \int_{\Omega} q_0 \varphi d\tau \quad \forall \varphi \in \mathcal{T}_0(Q_b), \quad (1.31) $$

where

$$ \mathcal{T}_0(Q_b) \stackrel{\text{def}}{=} \left\{ \begin{array}{l} \varphi \in \mathcal{D}(\Omega) \text{ such that there exists } \{\varphi_i\}_{i=1}^m \subset \mathcal{S} \\ \text{satisfying for } i \in \{1, \dots, m\}, \quad \varphi = \varphi_i \text{ in an open neighborhood of } \overline{\mathcal{B}_{i,0}} \end{array} \right\}. $$

7. The following energy inequality holds for almost every $t \in (0, T)$:

$$ \int_{\Omega} \left( \frac{1}{2} \rho(t, \cdot) |u(t, \cdot)|^2 + P(\rho(t, \cdot), \{1_{\mathcal{B}_i(t)}\}_{i=1}^m) \right) dy + \int_0^t \int_{\Omega} \left( 2\mu |\mathcal{D}(u) - \frac{1}{3} \operatorname{div}(u) \mathbb{I}_3|^2 + \nu |\operatorname{div}(u)|^2 \right) dy d\tau \\ + \sum_{i=1}^m \frac{2\pi}{3} \kappa_i |R_i(t)|^2 \leq - \int_0^t \int_{\Omega} \rho g \cdot u dy d\tau + \int_{\Omega} \left( \frac{1}{2} \frac{|q_0|^2}{\rho_0} 1_{\rho_0 > 0} + P(\rho_0, \{1_{\mathcal{B}_{i,0}}\}_{i=1}^m) \right) dy + \sum_{i=1}^m \frac{2\pi}{3} \kappa_i |R_{i,0}|^2, \quad (1.32) $$

where $P(\rho, \{1_{\mathcal{B}_i}\}_{i=1}^m)$ is the potential energy given by

$$ P(\rho, \{1_{\mathcal{B}_i}\}_{i=1}^m) \stackrel{\text{def}}{=} \frac{a_f \rho^{\gamma_f}}{\gamma_f - 1} \left( 1 - \sum_{i=1}^m 1_{\mathcal{B}_i} \right) + \sum_{i=1}^m \frac{a_i \rho^{\gamma_i}}{\gamma_i - 1} 1_{\mathcal{B}_i}. $$

We now present the main result of this document.

**Theorem 1.1.** Let $\Omega \subset \mathbb{R}^3$ be a bounded smooth domain and for $i \in \{1, \dots, m\}$, $\mathcal{B}_{i,0} \subset \Omega$ a ball of radius $R_{i,0} \in (0, \infty)$ and center $x_{i,0} \in \mathbb{R}^3$. Let $g \in L^\infty((0, T) \times \Omega)$, assume there exists $\sigma > 0$ such that

$$ \operatorname{dist}(\mathcal{B}_{i,0}, \partial\Omega) > 2\sigma, \quad \operatorname{dist}(\mathcal{B}_{i,0}, \mathcal{B}_{j,0}) > 2\sigma, \quad \forall i, j \in \{1, \dots, m\}, \quad i \neq j. $$

Then there exists $T > 0$ such that a bounded energy weak solution to (1.1)-(1.2), (1.13)-(1.15), (1.17)-(1.19) in the sense of Definition 1.1 exists on $[0, T]$. Moreover, such weak solution satisfies for all $t \in [0, T]$,

$$ R_i(t) > 0, \quad \operatorname{dist}(\mathcal{B}_i(t), \partial\Omega) \geq \sigma, \quad \operatorname{dist}(\mathcal{B}_i(t), \mathcal{B}_j(t)) \geq \sigma, \quad \forall i, j \in \{1, \dots, m\}, \quad i \neq j. \quad (1.33) $$

From now one, for the sake of simplicity, we prove Theorem 1.1 in the case of a single bubble immersed in a linearly viscous compressible fluid ($m = 1$). Due to the property (1.33) satisfied by the solution, the proof of Theorem 1.1 can then be easily extended to the case of several bubbles ($m > 1$).

**Remark 1.4.** Following the terminology of [14], we shall say that the velocity $u$, respectively the density $\rho$, is compatible with the system $\{\overline{\mathcal{B}}_i, \eta_i\}_{i=1}^m$, if it satisfies the condition (1.26), respectively (1.27), of Definition 1.1.
## 1.2 State of the art and discussion

The study of fluid-bubble interaction systems is a classical area of research in fluid mechanics. As discussed in the first part of the paper, the model we study is based on two prior works [24, 25] dealing with the derivation of averaged models for multiphase bubbly flows. More specifically, we aim to establish the existence of weak solutions for the microscopic model derived in [24]. In this model, the bubble dynamics (1.13)-(1.15) are derived considering a two phase compressible Navier-Stokes flow and assuming an infinite shear viscosity in the bubble phase. A common alternative approach to this framework is to couple incompressible Navier-Stokes equations for the fluid and compressible Euler equations for the bubble, imposing continuity of the stress at the interface. In this formulation, the pressure inside the bubble is usually assumed homogeneous and, in contrast with the model considered here, both viscosity and inertia effects of the internal fluid are neglected. The resulting free-surface relations read

$$ \left\{ \begin{array}{ll} p_i(t) = c_i(t)|\mathcal{B}_i(t)|^{-\gamma_b}, & t > 0, \\ \mathcal{T}(u_f, p_f)n = -p_i(t)n, & t > 0, \quad x \in \partial\mathcal{B}_i(t), \\ \partial\mathcal{B}_i(t) = \eta_f[t](\partial\mathcal{B}_{i,0}), & t > 0, \end{array} \right. \qquad (1.34) $$

where $\eta_f$ denotes the fluid flow map. Notably, contrary to our model, these relations allow tangential discontinuity of the velocity at the interface. Under radial symmetry of the initial data, an important subclass of solutions is described by the Rayleigh-Plesset equations [36, 34]. This model has been extensively used and analyzed in bubble dynamics, both theoretically and numerically, see [35] for further details. Variants of this model also exist, incorporating effects such as thermal evolution or modifications of the gas law [1]. The existence of weak solutions for a relaxed version of the interface conditions (1.34) and spherical bubbles, has been proven in [6]. It should be noted that, although the relations (1.34) do not impose any geometric constraint on the bubble shape, most of the available results deal with spherical bubbles.

Because of the no-slip conditions at the fluid-bubble interface and the description of the bubble dynamics by a finite set of ODEs, the fluid-bubble interaction problem we study in this paper shares many features with fluid-solid systems. More specifically, we can consider the bubble as an immersed structure with a velocity constrained to a finite number of modes. However, unlike fluid-solid systems, the bubble kinematics include an additional expansion-compression mode associated with variations in the bubble volume. From this perspective, our system is also related to the study of elastic structures, particularly when deformation is constrained to a finite number of modes.

A vast body of literature addresses the existence of weak and strong solutions for the interaction between rigid bodies and a viscous fluid governed by Navier-Stokes equations, for the most part in the incompressible setting [28, 39, 26, 9, 8, 19, 27, 38, 15, 40, 41, 21, 7] but also in the compressible setting [8, 14, 20, 23, 22, 37, 29, 32].

Regarding weak solutions, for a finite number of rigid bodies in a viscous incompressible fluid, their existence in dimension 2 and 3, until collision, was first established in [8], including an additional assumption of sufficiently small initial data in 3D. These results were then generalized to both incompressible and compressible settings in 2D and 3D, assuming this time only the absence of collision, in [9]. In the incompressible case, the previous results were extended for global in time solutions (including potential collisions) in [38] for 2D and [15] for 3D. In the compressible framework, results allowing collisions in 2D and 3D were given in [14].

Analogous results for the interaction between a viscous fluid and an elastic structure exist both in the incompressible [4] and compressible framework [3]). The case of an elastic structure described by a finite number of modes have also been treated in the incompressible framework, in [10] for the existence of weak solutions and [5] for strong solutions.

Our approach to proving the existence of weak solutions for the fluid-bubble system is mostly inspired by techniques from the theory of rigid bodies [38, 14], where a penalization method is employed by replacing the solid with a fluid of very high viscosity. However, unlike the $H^1$-type penalization used in the aforementioned works, we adopt an $L^2$-type penalization, by penalizing the fluid velocity within the bubble domain. Similar penalization strategies for rigid body motion have been considered as a numerical approach to deal with fluid-structure interaction problems [2] or to take account of slip at the fluid-solid interface [21, 32].

The idea of the proof is to introduce an approximate scheme which combines the theory of compressible fluids introduced by P.L Lions [31] and Feireisl [16] to get the strong convergence of the density (renormalized continuity equations, effective viscous flux, artificial pressure) together with an additional penalization term in the bubble domain to approximatively enforce the constraint on the fluid velocity.

The main novelty of our approach with respect to the rigid-body framework is that the velocity field associated to the motion of the bubbles is not divergence-free. As a result, the volume of the bubble is not constant in time and our analysis must address not only the possibility of collisions, either between bubbles or with the external boundary, but also the potential collapse of the bubbles. In particular, the equations for the kinematics and the dynamics of the bubbles are only defined as long as their radius is strictly positive. Another consequence of the variation of the volume of the bubbles is that contrary to the rigid body framework, the analysis of the non-linear terms, the advection term and the pressure term, needs to be performed in the whole domain and not only in the fluid part.

The rest of the paper is organized as follows. We introduce three levels of approximation schemes in Section 2. In Section 3, we describe some results on the transport equation, which are needed in all the levels of approximation. Section 4.1 and Sec-
tion 4.2 are dedicated to the construction and convergence analysis of a Faedo-Galerkin scheme associated to a finite dimensional approximation level. We discuss the limiting system associated to the high penalization term in Section 4.3 and vanishing viscosity in Section 4.4. Section 5 is devoted to the main part: the existence of weak solutions for the initial fluid-bubble system.

## 2 Approximate solutions

In this section, we state the existence results for the different levels of the approximation scheme. We begin with the $\delta$-level approximation. This step consists of adding an artificial pressure so that

$$p_{\delta} \stackrel{\text{def}}{=} p_{\delta}(\rho_{\delta}, 1_{\mathcal{B}_{\delta}}) = a_{f} \rho_{\delta}^{\gamma_{f}} (1 - 1_{\mathcal{B}_{\delta}}) + a_{b} \rho_{\delta}^{\gamma_{b}} 1_{\mathcal{B}_{\delta}} + \delta \rho_{\delta}^{\beta} \quad (2.1)$$

with

$$\delta > 0, \quad \beta \geq \max\{8, 2\gamma_f, 2\gamma_b\}. \quad (2.2)$$

We also introduce the corresponding potential energy:

$$P_{\delta} \stackrel{\text{def}}{=} P_{\delta}(\rho_{\delta}, 1_{\mathcal{B}_{\delta}}) = \frac{a_{f} \rho_{\delta}^{\gamma_{f}}}{\gamma_{f} - 1} (1 - 1_{\mathcal{B}_{\delta}}) + \frac{a_{b} \rho_{\delta}^{\gamma_{b}}}{\gamma_{b} - 1} 1_{\mathcal{B}_{\delta}} + \frac{\delta \rho_{\delta}^{\beta}}{\beta - 1}.$$

The $\delta$-level approximation problem reads as follows: find a triplet $(\mathcal{B}_{\delta}, \rho_{\delta}, u_{\delta})$ such that

1. There exists $\eta_{\delta} \in H^1([0, T]; \mathcal{M})$ such that, for all $t \in [0, T]$, $\mathcal{B}_{\delta}(t)$ is a ball of radius $R_{\delta}(t) \in (0, \infty)$ and center $x_{\delta}(t) \in \mathbb{R}^3$ given by

$$\mathcal{B}_{\delta}(t) = \eta_{\delta}[t](\mathcal{B}_{0}).$$

2. The velocity $u_{\delta} \in L^2(0, T; H^1(\Omega))$ and the density $\rho_{\delta} \in L^{\infty}(0, T; L^{\gamma_f}(\Omega))$ are compatible with the system $\{\overline{\mathcal{B}}_{\delta}, \eta_{\delta}\}$ in the sense of Remark 1.4.

3. The pair $(\rho_{\delta}, u_{\delta})$ satisfies the continuity equation (1.21) for all $\varphi \in \mathcal{D}((0, T) \times \mathbb{R}^3)$ provided that it is extended by $(0, 0)$ in $(\mathbb{R}^3 \setminus \Omega) \times (\mathbb{R}^3 \setminus \Omega)$.

4. Let

$$\overline{Q}_{\delta} \stackrel{\text{def}}{=} \{(t, x) \mid t \in [0, T], x \in \overline{\mathcal{B}}_{\delta}(t)\}. \quad (2.3)$$

For every test function $\varphi \in \mathcal{T}(Q_{\delta})$, the following weak form of the momentum equation holds:

$$-\int_{0}^{T} \int_{\Omega} \rho_{\delta} (u_{\delta} \cdot \partial_{t}\varphi + u_{\delta} \otimes u_{\delta} : \mathbb{D}(\varphi)) dy d\tau + \int_{0}^{T} \int_{\Omega} 2\mu_{\delta} \left( \mathbb{D}(u_{\delta}) - \frac{1}{3} \mathrm{div}(u_{\delta}) \mathbb{I}_{3} \right) : \left( \mathbb{D}(\varphi) - \frac{1}{3} \mathrm{div}(\varphi) \mathbb{I}_{3} \right) dy d\tau \\ + \int_{0}^{T} \int_{\Omega} (\nu_{\delta} \mathrm{div}(u_{\delta}) - p_{\delta}(\rho_{\delta}, 1_{\mathcal{B}_{\delta}})) \mathrm{div}(\varphi) dy d\tau = - \int_{0}^{T} \int_{\Omega} \rho_{\delta} g \cdot \varphi dy d\tau + \int_{0}^{T} \int_{\Omega} 1_{\mathcal{B}_{\delta}} \frac{\kappa_{b}}{R_{\delta}} \mathrm{div}(\varphi) dy d\tau, \quad (2.4)$$

where the viscosity coefficients are given by

$$\mu_{\delta} = (1 - 1_{\mathcal{B}_{\delta}})\mu_{f}, \quad \nu_{\delta} = (1 - 1_{\mathcal{B}_{\delta}})\nu_{f} + 1_{\mathcal{B}_{\delta}}\nu_{b}. \quad (2.5)$$

5. The initial conditions are satisfied in the sense of (1.30)-(1.31).

A weak solution of problem (1.1)-(1.17), in the sense of Definition 1.1, is obtained as a weak limit of solutions $(\mathcal{B}_{\delta}, \rho_{\delta}, u_{\delta})$ to the $\delta$-level approximation as $\delta \to 0$. The existence result of the approximate system reads:

**Proposition 2.1.** Let $\Omega \subset \mathbb{R}^3$ be a bounded smooth domain and $\sigma, M \in (0, \infty)$. Let $\mathcal{B}_0 \subset \Omega$ be a ball of radius $R_0 \in (0, \infty)$ and center $x_0 \in \mathbb{R}^3$, $\rho_{b,0} \in (0, \infty)$, $V_0, \omega_0 \in \mathbb{R}^3$ and $\Lambda_0 \in \mathbb{R}$. Assume that

$$\begin{align}
\rho_0 \in L^\beta(\Omega), \quad \rho_0 \ge 0 \text{ a.e. in } \Omega, \quad \rho_0 = \rho_{b,0} \text{ a.e. in } \mathcal{B}_0, \nonumber \\
q_0 \in L^{\frac{2\beta}{\beta+1}}(\Omega), \quad q_0 1_{\{\rho_0=0\}} = 0 \text{ a.e. in } \Omega, \quad \frac{|q_0|^2}{\rho_0} 1_{\rho_0>0} \in L^1(\Omega), \tag{2.6} \\
u_0 = V_0 + \omega_0 \times (x-x_0) + \frac{\Lambda_0}{3}(x-x_0) \text{ a.e. in } \mathcal{B}_0. \nonumber
\end{align}$$

Let $g \in L^\infty((0, T) \times \Omega)$ and $\delta \in (0, \infty)$. Assume that

$$\mathrm{dist}(\mathcal{B}_0, \partial\Omega) > 2\sigma,$$
and the initial energy satisfies

$$E_{\delta}(\rho_0, q_0, \mathbb{1}_{\mathcal{B}_0}) \stackrel{\text{def}}{=} \int_{\Omega} \left( \frac{1}{2} \frac{|q_0|^2}{\rho_0} \mathbb{1}_{\rho_0 > 0} + P_{\delta}(\rho_0, \mathbb{1}_{\mathcal{B}_0}) \right) dy + \frac{2}{3} \pi \kappa_b |R_0|^2 \leq M. \quad (2.7)$$

Then there exists $T > 0$, depending only on $(\sigma, M)$, such that the $\delta$-level approximation problem admits a bounded energy weak solution $(\mathcal{B}_\delta, \rho_\delta, u_\delta)$ on $[0, T]$, which satisfies, for almost every $t \in (0, T)$,

$$E_{\delta}(\rho_{\delta}(t), (\rho_{\delta}u_{\delta})(t), \mathbb{1}_{\mathcal{B}_{\delta}}(t)) + \int_{0}^{t} \int_{\Omega} \left( 2\mu_{\delta} \left| \mathbb{D}(u_{\delta}) - \frac{1}{3} \operatorname{div}(u_{\delta}) \mathbb{I}_{3} \right|^{2} + \nu_{\delta} \left| \operatorname{div}(u_{\delta}) \right|^{2} \right) dy d\tau \\ \leq - \int_{0}^{t} \int_{\Omega} \rho_{\delta} g \cdot u_{\delta} dy d\tau + E_{\delta}(\rho_{0}, q_{0}, \mathbb{1}_{\mathcal{B}_{0}}). \quad (2.8)$$

Moreover,

$$R_{\delta}(t) \geq \frac{R_{0}}{2}, \quad \operatorname{dist}(\mathcal{B}_{\delta}(t), \partial\Omega) \geq \frac{3\sigma}{2} \quad \forall t \in [0, T],$$

and provided that $\rho_\delta$ and $u_\delta$ are extended by 0 in $\mathbb{R}^3 \setminus \Omega$, they satisfy the renormalized continuity equation (1.28) for all $b$ satisfying (1.29) with $\lambda_0 < 1$, $-1 < \lambda_1 \le \frac{\beta}{2} - 1$.

To prove Proposition 2.1, we introduce the $\varepsilon$-level approximation problem, in which a dissipative term is added to the continuity equation:

$$\partial_t \rho_{\varepsilon} + \operatorname{div}(\rho_{\varepsilon} u_{\varepsilon}) = \varepsilon \Delta \rho_{\varepsilon} \text{ in } \mathcal{D}'((0,T) \times \Omega), \quad \nabla \rho_{\varepsilon} \cdot n = 0 \text{ on } \partial \Omega, \quad (2.9)$$

with $\varepsilon > 0$. To preserve the energy estimates, adding the dissipative term in (2.9) also requires including

$$\varepsilon \int_{0}^{T} \int_{\Omega} (\nabla u_{\varepsilon} \nabla \rho_{\varepsilon}) \varphi dy d\tau$$

in the momentum equation. The $\varepsilon$-level approximation reads as follows: find a triplet $(\mathcal{B}_\varepsilon, \rho_\varepsilon, u_\varepsilon)$ such that

1. There exists $\eta_\varepsilon \in H^1([0, T]; \mathcal{M})$ such that, for all $t \in [0, T]$, $\mathcal{B}_\varepsilon(t)$ is a ball of radius $R_\varepsilon(t) \in (0, \infty)$ and center $x_\varepsilon(t) \in \mathbb{R}^3$ given by $\mathcal{B}_\varepsilon(t) = \eta_\varepsilon[t](\mathcal{B}_0)$.

2. The velocity field $u_\varepsilon \in L^2(0, T; H^1(\Omega))$ and the density $\rho_\varepsilon \in L^\infty(0, T; L^\beta(\Omega)) \cap L^2(0, T; H^1(\Omega))$ are compatible with the system $\{\overline{\mathcal{B}}_\varepsilon, \eta_\varepsilon\}$ in the sense of Remark 1.4.

3. The pair $(\rho_\varepsilon, u_\varepsilon)$ satisfies the regularized continuity equation (2.9).

4. Let

$$\overline{Q}_{\varepsilon} \stackrel{\text{def}}{=} \{(t, x) \mid t \in [0, T], x \in \overline{\mathcal{B}_{\varepsilon}(t)}\}. \quad (2.10)$$

For every test function $\varphi \in \mathcal{T}(Q_\varepsilon)$, the following weak form of the momentum equation holds:

$$- \int_{0}^{T} \int_{\Omega} \rho_{\varepsilon} (u_{\varepsilon} \cdot \partial_{t} \varphi + u_{\varepsilon} \otimes u_{\varepsilon} : \mathbb{D}(\varphi)) dy d\tau + \int_{0}^{T} \int_{\Omega} 2\mu_{\varepsilon} \left( \mathbb{D}(u_{\varepsilon}) - \frac{1}{3} \operatorname{div}(u_{\varepsilon}) \mathbb{I}_{3} \right) : \left( \mathbb{D}(\varphi) - \frac{1}{3} \operatorname{div}(\varphi) \mathbb{I}_{3} \right) dy d\tau \\ + \int_{0}^{T} \int_{\Omega} (\nu_{\varepsilon} \operatorname{div}(u_{\varepsilon}) - p_{\delta}(\rho_{\varepsilon}, \mathbb{1}_{\mathcal{B}_{\varepsilon}})) \operatorname{div}(\varphi) dy d\tau + \varepsilon \int_{0}^{T} \int_{\Omega} (\nabla u_{\varepsilon} \nabla \rho_{\varepsilon}) \cdot \varphi dy d\tau \\ = - \int_{0}^{T} \int_{\Omega} \rho_{\varepsilon} g \cdot \varphi dy d\tau + \int_{0}^{T} \int_{\Omega} \mathbb{1}_{\mathcal{B}_{\varepsilon}} \frac{\kappa_{b}}{R_{\varepsilon}} \operatorname{div}(\varphi) dy d\tau, \quad (2.11)$$

where the viscosity coefficients are given by

$$\mu_{\varepsilon} = (1 - \mathbb{1}_{\mathcal{B}_{\varepsilon}})\mu_{f}, \quad \nu_{\varepsilon} = (1 - \mathbb{1}_{\mathcal{B}_{\varepsilon}})\nu_{f} + \mathbb{1}_{\mathcal{B}_{\varepsilon}}\nu_{b}. \quad (2.12)$$

5. Initial data are satisfied in the sense of (1.30)-(1.31).

A weak solution of the $\delta$-level approximation, in the sense of Proposition 2.1, is obtained as a weak limit of solutions $(\mathcal{B}_\delta, \rho_\delta, u_\delta)$ to the $\varepsilon$-level approximation as $\varepsilon \to 0$. The existence result of the approximate system reads:
**Proposition 2.2.** Let $\Omega \subset \mathbb{R}^3$ be a bounded smooth domain and $\sigma$, $M$, $\bar{\rho}$, $\underline{\rho} \in (0, \infty)$. Let $\mathcal{B}_0 \subset \Omega$ be a ball of radius $R_0 \in (0, \infty)$ and center $x_0 \in \mathbb{R}^3$, $V_0$, $\omega_0 \in \mathbb{R}^3$ and $\Lambda_0 \in \mathbb{R}$. Assume that

$$ \rho_0 \in W^{1,\infty}(\Omega), \quad 0 < \underline{\rho} \le \rho_0 \le \bar{\rho}, \quad q_0 \in L^2(\Omega), \\ u_0 = V_0 + \omega_0 \times (x-x_0) + \frac{\Lambda_0}{3}(x-x_0) \text{ a.e. in } \mathcal{B}_0, \quad q_0 = \rho_0 u_0 \text{ a.e. in } \mathcal{B}_0. $$

Let $g \in L^\infty((0,T) \times \Omega)$ and $\delta, \varepsilon \in (0, \infty)$. Assume that

$$ \mathrm{dist}(\mathcal{B}_0, \partial\Omega) > 2\sigma, $$

and that the initial energy satisfies

$$ \tilde{E}_{\delta}(\rho_0, q_0) \stackrel{\text{def}}{=} \int_{\Omega} \left( \frac{1}{2} |q_0|^2 + \delta \frac{\rho_0^{\beta}}{\beta - 1} \right) dy + \frac{2}{3} \kappa_b \pi |R_0|^2 \le M. \quad (2.13) $$

Then there exists $T > 0$, depending only on $(\sigma, M, \underline{\rho}, \bar{\rho})$, such that the $\varepsilon$-level approximation problem admits a bounded energy weak solution $(\mathcal{B}_\varepsilon, \rho_\varepsilon, u_\varepsilon)$, which satisfies, for almost every $t \in (0, T)$,

$$ \begin{aligned} \tilde{E}_{\delta}(\rho_{\varepsilon}(t), (\rho_{\varepsilon}u_{\varepsilon})(t)) &+ \int_{0}^{t} \int_{\Omega} \left( 2\mu_{\varepsilon} \left| \mathbb{D}(u_{\varepsilon}) - \frac{1}{3} \mathrm{div}(u_{\varepsilon}) \mathbb{I}_{3} \right|^{2} + \nu_{\varepsilon} \left| \mathrm{div}(u_{\varepsilon}) \right|^{2} \right) \mathrm{d}y\mathrm{d}\tau + \delta\varepsilon\beta \int_{0}^{t} \int_{\Omega} (\rho_{\varepsilon})^{\beta-2} |\nabla\rho_{\varepsilon}|^{2} \mathrm{d}y\mathrm{d}\tau \\ &\le a_{f} \int_{0}^{t} \int_{\Omega} (1 - \mathbb{1}_{\mathcal{B}_{\varepsilon}}) \rho_{\varepsilon}^{\gamma_{f}} \mathrm{div}(u_{\varepsilon}) \mathrm{d}y\mathrm{d}\tau + a_{b} \int_{0}^{t} \int_{\Omega} \mathbb{1}_{\mathcal{B}_{\varepsilon}} \rho_{\varepsilon}^{\gamma_{b}} \mathrm{div}(u_{\varepsilon}) \mathrm{d}y\mathrm{d}\tau - \int_{0}^{t} \int_{\Omega} \rho_{\varepsilon}g \cdot u_{\varepsilon} \mathrm{d}y\mathrm{d}\tau + \tilde{E}_{\delta}(\rho_{0}, q_{0}). \end{aligned} \quad (2.14) $$

Moreover,

$$ R_{\varepsilon}(t) \geq \frac{R_0}{2}, \quad \mathrm{dist}(\mathcal{B}_{\varepsilon}(t), \partial\Omega) \geq \frac{3\sigma}{2} \quad \forall t \in [0, T], $$

and the density $\rho_\varepsilon$ satisfies

$$ \partial_t \rho_\varepsilon, \Delta \rho_\varepsilon \in L^{\frac{5\beta-3}{4\beta}}((0,T) \times \Omega), \quad \sqrt{\varepsilon} \|\nabla \rho_\varepsilon\|_{L^2((0,T) \times \Omega)} + \|\rho_\varepsilon\|_{L^{\beta+1}((0,T) \times \Omega)} \le c, \quad (2.15) $$

where c is a positive constant uniformly independent of $\varepsilon$.

To prove Proposition 2.2, we consider an additional level of approximation. In this n-level scheme, a penalization term is added to the momentum equation, which is designed to approximate the motion of the bubble. We seek a triplet $(\mathcal{B}_n, \rho_n, u_n)$ satisfying

1. For all $t \in [0, T]$, $\mathcal{B}_n(t)$ is a ball of radius $R_n(t) \in (0, \infty)$ and center $x_n(t) \in \mathbb{R}^3$ such that

$$ \chi_n \stackrel{\text{def}}{=} \mathbb{1}_{\mathcal{B}_n} \quad (2.16) $$

belongs to $\mathcal{C}([0, T]; L^p(\mathbb{R}^3))$ for all $1 \le p < \infty$.

2. The velocity field $u_n \in L^2(0, T; H^1(\Omega))$ and the density function $\rho_n \in L^\infty(0, T; L^\beta(\Omega)) \cap L^2(0, T; H^1(\Omega))$ satisfy the regularized continuity equation (2.9).

3. Let for all $t \in [0, T]$, $\Pi_n(t) : L^2(\Omega) \to S$ defined for any $\varphi \in L^2(\Omega)$ by

$$ (\Pi_n(t)\varphi)(x) \stackrel{\text{def}}{=} \frac{3}{4\pi R_n(t)^3} \int_{\Omega} \chi_n \varphi \mathrm{d}y + \frac{15}{8\pi R_n(t)^5} \left( \int_{\Omega} \chi_n(y-x_n(t)) \times \varphi \mathrm{d}y \right) \times (x-x_n(t)), \\ + \frac{5}{4\pi R_n(t)^5} \left( \int_{\Omega} \chi_n(y-x_n(t)) \cdot \varphi \mathrm{d}y \right) (x-x_n(t)) \quad \forall x \in \mathbb{R}^3, \quad (2.17) $$

with

$$ R_n = \left( \frac{3}{4\pi} \int_{\mathbb{R}^3} \chi_n \mathrm{d}y \right)^{\frac{1}{3}}, \quad x_n = \frac{3}{4\pi R_n^3} \left( \int_{\mathbb{R}^3} \chi_n y \mathrm{d}y \right). $$
For every test function $\varphi \in \mathcal{D}((0, T) \times \Omega)$, the following weak form of the momentum equation holds:

$$ -\int_0^T \int_\Omega \rho_n (u_n \cdot \partial_t \varphi + u_n \otimes u_n : \mathcal{D}(\varphi)) \, dyd\tau + \int_0^T \int_\Omega 2\mu_n \left( \mathcal{D}(u_n) - \frac{1}{3} \mathrm{div}(u_n) \mathbb{I}_3 \right) : \left( \mathcal{D}(\varphi) - \frac{1}{3} \mathrm{div}(\varphi) \mathbb{I}_3 \right) \, dyd\tau \\ + \int_0^T \int_\Omega (\nu_n \mathrm{div}(u_n) - p_\delta(\rho_n, \chi_n)) \mathrm{div}(\varphi) \, dyd\tau + \varepsilon \int_0^T \int_\Omega (\nabla u_n \nabla \rho_n) \cdot \varphi \, dyd\tau \\ + n \int_0^T \int_\Omega \chi_n (u_n - \Pi_n u_n)(\varphi - \Pi_n \varphi) \, dyd\tau = -\int_0^T \int_\Omega \rho_n g \cdot \varphi \, dyd\tau + \int_0^T \int_\Omega \frac{\kappa_b}{R_n} \chi_n \mathrm{div}(\varphi) \, dyd\tau, \quad (2.18) $$

where the viscosity coefficients are given by

$$ \mu_n = (1 - \chi_n)\mu_f + n\chi_n, \quad \nu_n = (1 - \chi_n)\nu_f + \chi_n\nu_b. \quad (2.19) $$

4. The color function $\chi_n$ satisfies

$$ \begin{gathered} \partial_t \chi_n + \tilde{\Pi}_n u_n \cdot \nabla \chi_n = 0 \quad \text{in } \mathcal{D}'((0, T) \times \mathbb{R}^3), \\ \chi_n(0, x) = \mathbb{1}_{\mathcal{B}_0}(x) \quad \text{for a.e. } x \in \mathbb{R}^3, \end{gathered} \quad (2.20) $$

where for $t \in [0, T]$, $\tilde{\Pi}_n(t) : L^2(\Omega) \to \mathcal{S}$ is defined for $\varphi \in L^2(\Omega)$ by

$$ (\tilde{\Pi}_n(t)\varphi)(x) \stackrel{\text{def}}{=} \frac{3}{4\pi R_n(t)^3} \int_{\Omega} \chi_n \varphi \, dy + \frac{5}{4\pi R_n(t)^5} \left( \int_{\Omega} \chi_n (y - x_n(t)) \cdot \varphi \, dy \right) (x - x_n(t)) \quad \forall x \in \mathbb{R}^3. \quad (2.21) $$

5. The initial conditions are satisfied in the following sense:

$$ \lim_{t \to 0^+} \int_{\Omega} \rho_n(t) \varphi \, d\tau = \int_{\Omega} \rho_0 \varphi \, d\tau \quad \forall \varphi \in \mathcal{D}(\Omega), \quad (2.22) $$

$$ \lim_{t \to 0^+} \int_{\Omega} (\rho_n u_n)(t) \varphi \, d\tau = \int_{\Omega} q_0 \varphi \, d\tau \quad \forall \varphi \in \mathcal{D}(\Omega). \quad (2.23) $$

A weak solution of the $\varepsilon$-level approximation problem in the sense of Proposition 2.2 will be obtained as a weak limit of the solution $(\mathcal{B}_n, \rho_n, u_n)$ of the $n$-level approximation problem as $n \to \infty$. The existence result of the approximate system reads:

**Proposition 2.3.** Let $\Omega \subset \mathbb{R}^3$ be a bounded smooth domain and $\sigma, M, \bar{\rho}, \underline{\rho} \in (0, \infty)$. Let $\mathcal{B}_0$ be a ball of radius $R_0 \in (0, \infty)$ and center $x_0 \in \mathbb{R}^3$. Assume that

$$ \rho_0 \in W^{1,\infty}(\Omega) \text{ satisfies } 0 < \underline{\rho} \le \rho_0 \le \bar{\rho}, \quad q_0 \in L^2(\Omega). \quad (2.24) $$

Let $g \in L^\infty((0, T) \times \Omega)$, $\delta, \varepsilon \in (0, \infty)$ and $n \ge 0$. Assume that

$$ \mathrm{dist}(\mathcal{B}_0, \partial\Omega) > 2\sigma, $$

and that the initial energy satisfies

$$ \tilde{E}_{\delta}(\rho_0, q_0) = \int_{\Omega} \left( \frac{1}{2} |q_0|^2 + \delta \frac{\rho_0^{\beta}}{\beta - 1} \right) \le M. $$

Then there exists $T > 0$, depending only on $(\sigma, M)$, such that the $n$-level approximation problem admits a bounded energy weak solution $(\mathcal{B}_n, \rho_n, u_n)$, which satisfies, for almost every $t \in (0, T)$,

$$ \begin{aligned} \tilde{E}_{\delta}(\rho_n(t), (\rho_n u_n)(t)) + \int_0^t \int_{\Omega} \left( 2\mu_n \left| \mathcal{D}(u_n) - \frac{1}{3} \mathrm{div}(u_n) \mathbb{I}_3 \right|^2 + \nu_n |\mathrm{div}(u_n)|^2 \right) \, dyd\tau + \delta\varepsilon \int_0^t \int_{\Omega} (\rho_n)^{\beta-2} |\nabla \rho_n|^2 \, dyd\tau \\ + n \int_0^t \int_{\Omega} \chi_n |u_n - \Pi_n u_n|^2 \, dyd\tau \le a_f \int_0^t \int_{\Omega} (1-\chi_n) \rho_n^{\gamma_f} \mathrm{div}(u_n) \, dyd\tau + a_b \int_0^t \int_{\Omega} \chi_n \rho_n^{\gamma_b} \mathrm{div}(u_n) \, dyd\tau \\ - \int_0^t \int_{\Omega} \rho_n g \cdot u_n \, dyd\tau + \int_0^t \int_{\Omega} \chi_n \frac{\kappa_b}{R_n} \mathrm{div}(u_n) \, dyd\tau + \tilde{E}_{\delta}(\rho_0, q_0). \end{aligned} \quad (2.25) $$

Moreover,

$$ R_n(t) \ge \frac{R_0}{2}, \quad \mathrm{dist}(\mathcal{B}_n(t), \partial\Omega) \ge \frac{3\sigma}{2} \quad \forall t \in [0, T]. $$
**Remark 2.1.** Note that, in addition to the penalization term in the momentum equation (2.18), we introduce in (2.19) a shear viscosity inside the bubble domain equal to n. Formally, in the limit $n \to \infty$, this corresponds to imposing

$$D(u_{\varepsilon}) = \frac{1}{3} \mathrm{div}(u_{\varepsilon}) \mathbb{I}_3. \qquad (2.26)$$

Since this constraint is already satisfied by the bubble velocity, it is redundant from a purely theoretical standpoint. However, from a convergence perspective, it may be interesting to examine whether this $H^1$-type penalization improves the convergence of the penalized model toward the exact one. Nevertheless, as discussed in [24], the constraint (2.26) alone is not sufficient to guarantee the preservation of the spherical shape of the bubbles. In particular, to recover the bubble velocity, it would be necessary to additionally penalize a fourth compatible mode by introducing the following term in the momentum equation (2.18):

$$n \int_0^T \int_{\Omega} \Pi_{s,n}(t) u_n \cdot \Pi_{s,n}(t) \varphi dy d\tau,$$

where, for $t \in [0, T]$, the operator $\Pi_{s,n}(t) : L^2(\Omega) \to \mathbb{R}^3$ is defined for $\varphi \in L^2(\Omega)$ by

$$\Pi_{s,n}(t)(\varphi) \stackrel{\mathrm{def}}{=} \int_{\Omega} \left( (y - x_n(t)) \otimes (y - x_n(t)) - \frac{1}{2} |y - x_n(t)|^2 \right) \varphi dy.$$

To prove Proposition 2.3, we need a final level of approximation, the N-level approximation, obtained via a Faedo-Galerkin approximation scheme. Let $\{\psi_i\}_{i \ge 1} \subset \mathcal{D}(\bar{\Omega})$ be an orthonormal basis of $L^2(\Omega)$ dense in $\{v \in W^{1,p}(\Omega) | v = 0 \text{ on } \partial\Omega\}$ for any $1 \le p \le \infty$. We set

$$X_N = \mathrm{span}\{\{\psi_i\}_{1 \le i \le N}\}. \qquad (2.27)$$

Since $X_N$ is a finite-dimensional space, norms on $X_N$ induced by $W^{k,p}$ norms, $k \in \mathbb{N}$, $1 \le p \le \infty$ are equivalent. The task is to find a triplet $(\mathcal{B}_N, \rho_N, u_N)$ satisfying:

1. For all $t \in [0, T]$, $\mathcal{B}_N(t)$ is a ball of radius $R_N(t) \in (0, \infty)$ and center $x_N(t) \in \mathbb{R}^3$ such that

$$\chi_N \stackrel{\mathrm{def}}{=} \mathbb{1}_{\mathcal{B}_N} \qquad (2.28)$$

belongs to $\mathcal{C}([0, T]; L^p(\Omega))$ for all $1 \le p < \infty$.

2. There exists $\{\alpha_i\}_{i=1}^N \subset \mathcal{C}([0, T])$ such that

$$u_N \stackrel{\mathrm{def}}{=} \sum_{i=1}^N \alpha_i \psi_i.$$

3. The velocity field $u_N \in \mathcal{C}([0, T]; X_N)$ and the density function $\rho_N \in L^2(0, T; H^2(\Omega)) \cap H^1(0, T; L^2(\Omega))$ satisfy the regularized continuity equation (2.9).

4. For all $\varphi \in \mathcal{D}(0, T; X_N)$, the following weak formulation of the momentum equation holds:

$$- \int_0^T \int_{\Omega} \rho_N (u_N \cdot \partial_t \varphi + u_N \otimes u_N : \mathcal{D}(\varphi)) dy d\tau + \int_0^T \int_{\Omega} 2\mu_N \left( \mathcal{D}(u_N) - \frac{1}{3} \mathrm{div}(u_N) \mathbb{I}_3 \right) : \left( \mathcal{D}(\varphi) - \frac{1}{3} \mathrm{div}(\varphi) \mathbb{I}_3 \right) dy d\tau \\ + \int_0^T \int_{\Omega} (\nu_N \mathrm{div}(u_N) - p_{\delta}(\rho_N, \chi_N)) \mathrm{div}(\varphi) dy d\tau + \varepsilon \int_0^T \int_{\Omega} (\nabla u_N \nabla \rho_N) \cdot \varphi dy d\tau \\ + n \int_0^T \int_{\Omega} \chi_N (u_N - \Pi_N u_N) \cdot (\varphi - \Pi_N \varphi) dy d\tau = - \int_0^T \int_{\Omega} \rho_N g \cdot \varphi dy d\tau + \int_0^T \int_{\Omega} \chi_N \frac{\kappa_b}{R_N} \mathrm{div}(\varphi) dy d\tau. \quad (2.29)$$

where $\Pi_N$ is defined as in (2.17) with $\chi_n$ replaced by $\chi_N$ and the viscosity coefficients are given by

$$\mu_N \stackrel{\mathrm{def}}{=} (1 - \chi_N)\mu_f + n\chi_N, \quad \nu_N \stackrel{\mathrm{def}}{=} (1 - \chi_N)\nu_f + \chi_N\nu_b. \qquad (2.30)$$

5. The color function $\chi_N$ satisfies

$$\partial_t \chi_N + \tilde{\Pi}_N u_N \cdot \nabla \chi_N = 0 \text{ in } \mathcal{D}'((0, T) \times \mathbb{R}^3), \\ \chi_N(0, x) = \mathbb{1}_{\mathcal{B}_0}(x) \text{ for a.e. } x \in \mathbb{R}^3,$$

where $\tilde{\Pi}_N$ is defined as in (2.21) with $\chi_n$ replaced by $\chi_N$.
6. The initial conditions are satisfied in the sense (2.22) - (2.23).

A weak solution of the n-level approximation problem, in the sense of Proposition 2.3, is obtained as a weak limit of solutions $(\mathcal{B}_N, \rho_N, u_N)$ of the Faedo-Galerkin scheme as $N \to \infty$. The existence result of the approximate system reads:

**Proposition 2.4.** Let $\Omega \subset \mathbb{R}^3$ be a bounded smooth domain and $\sigma$, $M$, $\bar{\rho}$, $\underline{\rho} \in (0, \infty)$. Let $\mathcal{B}_0$ be a ball of radius $R_0 \in (0, \infty)$ and center $x_0 \in \mathbb{R}^3$. Assume that

$$ \rho_0 \in W^{1,\infty}(\Omega), \quad 0 \le \underline{\rho} \le \rho_0 \le \bar{\rho}, \quad u_0 \in X_N, \quad q_0 = \rho_0 u_0 \text{ a.e. in } \Omega. \quad (2.31) $$

Let $g \in L^\infty((0,T) \times \Omega)$, $\delta, \varepsilon \in (0, \infty)$, $n \ge 0$ and $N \ge 1$. Assume that

$$ \mathrm{dist}(\mathcal{B}_0, \partial\Omega) > 2\sigma, $$

and that the initial energy satisfies

$$ \tilde{E}_{\delta}(\rho_0, q_0) \stackrel{\text{def}}{=} \int_{\Omega} \left( \frac{1}{2} \frac{|q_0|^2}{\rho_0} + \delta \frac{\rho_0^{\beta}}{\beta - 1} \right) dy \le M. $$

Then, there exists $T > 0$, depending only on $(\sigma, M, \bar{\rho}, \underline{\rho})$, such that the Faedo-Galerkin scheme admits a bounded energy weak solution $(\mathcal{B}_N, \rho_N, u_N)$, which satisfies, for all $t \in (0, T)$,

$$ \begin{aligned} \tilde{E}_{\delta}(\rho_N(t), (\rho_N u_N)(t)) + \int_0^t \int_{\Omega} \left( 2\mu_N \left| \mathcal{D}(u_N) - \frac{1}{3} \mathrm{div}(u_N) \mathbb{I}_3 \right|^2 + \nu_N |\mathrm{div}(u_N)|^2 \right) dy d\tau + \delta \varepsilon \int_0^t \int_{\Omega} \rho_N^{\beta-2} |\nabla \rho_N|^2 dy d\tau \\ + n \int_0^t \int_{\Omega} \chi_N |u_N - \Pi_N u_N|^2 dy d\tau \le a_f \int_0^t \int_{\Omega} (1 - \chi_N) \rho_N^{\gamma_f} \mathrm{div}(u_N) dy d\tau + a_b \int_0^t \int_{\Omega} \chi_N \rho_N^{\gamma_b} \mathrm{div}(u_N) dy d\tau \\ - \int_0^t \int_{\Omega} \rho_N g \cdot u_N dy d\tau + \int_0^t \int_{\Omega} \chi_N \frac{\kappa_b}{R_N} \mathrm{div}(u_N) dy d\tau + \tilde{E}_{\delta}(\rho_0, q_0) \end{aligned} \quad (2.32) $$

Moreover,

$$ R_N(t) \ge \frac{R_0}{2}, \quad \mathrm{dist}(\mathcal{B}_N(t), \partial\Omega) \ge \frac{3\sigma}{2} \quad \forall t \in [0, T]. $$

## 3 A nonlinear transport equation

This section is devoted to the nonlinear transport equation

$$ \begin{cases} \partial_t \chi + \tilde{\Pi} u \cdot \nabla \chi = 0 \text{ in } \mathcal{D}'((0,T) \times \mathbb{R}^3), \\ \chi|_{t=0} = \mathbb{1}_{\mathcal{B}_0} \text{ in } \mathbb{R}^3, \end{cases} \quad (3.1) $$

where $\mathcal{B}_0$ is a ball of radius $R_0 \in (0, \infty)$ and center $x_0 \in \mathbb{R}^3$, and $\tilde{\Pi}$ is defined analogously to $\tilde{\Pi}_n$ in (2.21), with $\chi_n$ replaced by $\chi$. We first address the regular case.

**Proposition 3.1** (Well-posedness). Let $u \in \mathcal{C}([0, T]; \mathcal{D}(\bar{\Omega}))$. The following statements hold.

* If $T$ satisfies

$$ T \|u\|_{L^\infty(0,T;L^2(\Omega))} \le R_0 \sqrt{\frac{\pi}{5}} R_0, \quad (3.2) $$

then there exists a unique solution $\chi$ of (3.1), which satisfies

$$ \chi \in L^\infty((0,T) \times \mathbb{R}^3) \cap \mathcal{C}([0,T], L^p(\mathbb{R}^3)) $$

for all $1 \le p < \infty$.

* There exist $x_b, R_b \in \mathcal{C}^1([0, T])$ such that for all $t \in [0, T]$,

$$ \chi(t, \cdot) = \mathbb{1}_{\mathcal{B}(t)} \text{ where } \mathcal{B}(t) = \mathcal{B}(x_b(t), R_b(t)) \text{ and } R_b(t) \ge \frac{R_0}{2}. $$
• For all $t \in [0, T]$, $\mathcal{B}(t) = \tilde{\eta}[t](\mathcal{B}_0)$, where $\tilde{\eta} \in \mathcal{C}^1([0, T]; \mathcal{M})$ is the propagator associated with $\tilde{\Pi}u$ given by

$$ \tilde{\eta}[t](x) \stackrel{\text{def}}{=} x_b(t) + \frac{R_b(t)}{R_0}(x - x_0) \quad \forall (t, x) \in [0, T] \times \mathbb{R}^3. \qquad (3.3) $$

Proof of proposition 3.1. We can assume $u \in \mathcal{C}([0, T]; \mathcal{D}(\mathbb{R}^3))$ without loss of generality. Assume now that we have found a solution $\chi$ to (3.1). Then we can see (3.1) as a linear transport equation with prescribed velocity field $\tilde{\Pi}u \in \mathcal{C}([0, T], \mathcal{S})$ expressed by

$$ (\tilde{\Pi}u)(t, x) = V(t) + \frac{\Lambda(t)}{3}(x - x_b(t)) \quad \forall (t, x) \in [0, T] \times \mathbb{R}^3. $$

By the method of characteristics, we obtain

$$ \chi(t, \tilde{\eta}[t](x)) = \mathbb{1}_{\mathcal{B}_0}(x) \quad \forall (t, x) \in [0, T] \times \mathbb{R}^3, \qquad (3.4) $$

where $\tilde{\eta}[t](x)$ is the flow associated to $\tilde{\Pi}u$ given by

$$ \tilde{\eta}[t](x) \stackrel{\text{def}}{=} x_b(t) + \frac{R_b(t)}{R_0}(x - x_0) \quad \forall (t, x) \in [0, T] \times \mathbb{R}^3, $$

and $(x_b, R_b)$ is the unique solution of the ODE

$$ \begin{cases} \dot{x}_b = V, & \dot{R}_b = \frac{\Lambda}{3} R_b, \\ x_b(0) = x_0, & R_b(0) = R_0. \end{cases} $$

In particular, for all $t \in [0, T]$,

$$ \chi(t, \cdot) = \mathbb{1}_{\mathcal{B}(x_b(t), R_b(t))}. $$

We then substitute (3.3) and (3.4) in the expression for $\tilde{\Pi}u$ to obtain for $(t, x) \in [0, T] \times \mathbb{R}^3$,

$$ (\tilde{\Pi}u)(t, x) = \frac{3}{4\pi R_0^3} \int_{\mathcal{B}_0} \mathbb{1}_{\Omega}(\tilde{\eta}[t](y))u(t, \tilde{\eta}[t](y))dy + \frac{5}{4\pi R_0^4 R_b(t)} \left( \int_{\mathcal{B}_0} \mathbb{1}_{\Omega}(\tilde{\eta}[t](y))(y - x_0) \cdot u(t, \tilde{\eta}[t](y))dy \right) (x - x_b(t)). $$

In particular, this implies that $X_b \stackrel{\text{def}}{=} (x_b, R_b)$ must satisfy the nonlinear ODE

$$ \dot{X}_b = f(t, X_b) \text{ with } X_b(0) = (x_0, R_0), \qquad (3.5) $$

where

$$ f: \begin{cases} [0, T] \times (\mathbb{R}^3 \times (0, \infty)) \to \mathbb{R}^3 \times (0, \infty) \\ (t, X) \mapsto (f_x(t, X), f_R(t, X)) \end{cases} $$

is defined, for $X = (x, R) \in \mathbb{R}^3 \times (0, \infty)$, by

$$ f_x(t, X) \stackrel{\text{def}}{=} \frac{3}{4\pi R_0^3} \int_{\mathcal{B}_0} \mathbb{1}_{\Omega}(\tilde{\eta}_X(y))u(t, \tilde{\eta}_X(y))dy, \quad f_R(t, X) \stackrel{\text{def}}{=} \frac{5}{4\pi R_0^4} \int_{\mathcal{B}_0} \mathbb{1}_{\Omega}(\tilde{\eta}_X(y))(y - x_0) \cdot u(t, \tilde{\eta}_X(y))dy \qquad (3.6) $$

and

$$ \tilde{\eta}_X(y) \stackrel{\text{def}}{=} x + \frac{R}{R_0}(y - x_0) \quad \forall y \in \mathbb{R}^3. \qquad (3.7) $$

By the Cauchy-Lipschitz theorem, there exists a unique maximal $\mathcal{C}^1$ solution to equation (3.5), provided that $f$ is continuous in $(t, X)$ and locally Lipschitz in $X$. The continuity is obvious. Let $t \in [0, T]$ and $X_i \stackrel{\text{def}}{=} (x_i, R_i)$, $i \in \{1, 2\}$, we write

$$ f_x(t, X_1) - f_x(t, X_2) = M_1(t) + M_2(t), $$

with

$$ M_1(t) \stackrel{\text{def}}{=} \frac{3}{4\pi R_0^3} \int_{\mathcal{B}_0} \mathbb{1}_{\Omega}(\tilde{\eta}_{X_1}(y))(u(t, \tilde{\eta}_{X_1}(y)) - u(t, \tilde{\eta}_{X_2}(y))) dy, $$

$$ M_2(t) \stackrel{\text{def}}{=} \frac{3}{4\pi R_0^3} \int_{\mathcal{B}_0} u(t, \tilde{\eta}_{X_2}(y))(\mathbb{1}_{\Omega}(\tilde{\eta}_{X_1}(y)) - \mathbb{1}_{\Omega}(\tilde{\eta}_{X_2}(y))) dy. $$
On the one hand,

$$
\begin{align}
|M_1(t)| &\le \frac{3}{4\pi R_0^3} \int_{\mathcal{B}_0} |x_1 - x_2 + \frac{R_1 - R_2}{R_0}(y - x_0)| \int_0^1 \left| \nabla u \left( t, s(x_1 + \frac{R_1}{R_0}(y - x_0)) + (1-s)(x_2 + \frac{R_2}{R_0}(y - x_0)) \right) \right| ds dy \nonumber \\
&\le \|\nabla u\|_{L^\infty((0,T)\times\Omega)} (|x_1 - x_2| + |R_1 - R_2|). \tag{3.8}
\end{align}
$$

On the other hand, we have

$$
|M_2(t)| \le \frac{3}{4\pi R_0^3} \|u\|_{L^\infty((0,T)\times\Omega)} \int_{\mathcal{B}_0} |1_\Omega(\tilde{\eta}_{X_1}(y)) - 1_\Omega(\tilde{\eta}_{X_2}(y))| \, dy
$$

Since for all $y \in \mathcal{B}_0$,

$$
|\tilde{\eta}_{X_1}(y) - \tilde{\eta}_{X_2}(y)| \le |x_1 - x_2| + |R_1 - R_2| \stackrel{\text{def}}{=} d_{1,2},
$$

we deduce that the integrand is nonzero only if $\tilde{\eta}_{X_1}(y)$ and $\tilde{\eta}_{X_2}(y)$ lie in a $d_{1,2}$-neighborhood $\mathcal{U}$ of $\partial\Omega$. Hence,

$$
\begin{align}
|M_2(t)| &\le \frac{3}{4\pi R_0^3} \|u\|_{L^\infty((0,T)\times\Omega)} \left( \int_{\tilde{\eta}_{X_1}^{-1}(\mathcal{U})} dy + \int_{\tilde{\eta}_{X_2}^{-1}(\mathcal{U})} dy \right), \nonumber \\
&\le \frac{3}{4\pi} \|u\|_{L^\infty((0,T)\times\Omega)} \left( \frac{1}{R_1^3} + \frac{1}{R_2^3} \right) |\mathcal{U}|, \nonumber \\
&\lesssim \|u\|_{L^\infty((0,T)\times\Omega)} \left( \frac{1}{R_1^3} + \frac{1}{R_2^3} \right) (|x_1 - x_2| + |R_1 - R_2|), \tag{3.9}
\end{align}
$$

where the hidden constant is independent of $X$. Following similar arguments used for $f_x$, we deduce that

$$
|f_R(t, X_1) - f_R(t, X_2)| \lesssim \|\nabla u\|_{L^\infty((0,T)\times\Omega)} \left(\frac{1}{R_1^3} + \frac{1}{R_2^3}\right) (|x_1 - x_2| + |R_1 - R_2|). \quad (3.10)
$$

We therefore obtain that $f$ is locally Lipschitz in $X$, and conclude that equation (3.5) admits a unique $C^1$ maximal solution $(x_b, R_b)$. Moreover, this solution can be extended to $[0, T]$ provided that

$$
R_b(t) > 0 \quad \forall t \in [0, T].
$$

From (3.5) and (3.6), we obtain the estimate

$$
\| \dot{R}_b \|_{L^1(0,T)} \le T \| \dot{R}_b \|_{L^\infty(0,T)} \le \sqrt{\frac{5}{\pi}} \frac{T}{2R_0 \sqrt{R_0}} \| u \|_{L^\infty(0,T; L^2(\Omega))}.
$$

Since

$$
\begin{align}
R_b(t) &= \int_0^t \dot{R}_b(s) ds + R_0 \quad \forall t \in [0, T], \nonumber \\
T \|u\|_{L^\infty(0,T;L^2(\Omega))} &\le R_0 \sqrt{\frac{\pi}{5}} R_0, \tag{3.11}
\end{align}
$$

it follows that if

then

$$
R_b(t) \ge \frac{R_0}{2} > 0 \quad \forall t \in [0, T].
$$

This completes the proof.

<figure><img src="image_3.png" /></figure>

We now analyse the stability properties of the nonlinear transport equation (3.1).

**Proposition 3.2** (Strong sequential continuity). Let {$u_k, \chi_k$} be a bounded sequence in $C([0, T]; D(\bar{\Omega})) \times L^\infty((0, T) \times \Omega)$ such that

$$
\star u_k \to u \text{ in } C([0, T]; D(\bar{\Omega})), \tag{3.12}
$$

$$
\star \partial_t \chi_k + \tilde{\Pi}_k u_k \cdot \nabla \chi_k = 0 \text{ in } \mathcal{D}((0, T) \times \mathbb{R}^3), \qquad (3.13)
$$

$$
\chi_k|_{t=0} = \mathbb{1}_{\mathcal{B}_0} \quad \text{a.e. in } \mathbb{R}^3, \quad \mathcal{B}_0 = \mathcal{B}(x_0, R_0), \quad x_0 \in \mathbb{R}^3, \quad R_0 \in (0, \infty), \tag{3.14}
$$

$$
\star \chi_k = \mathbb{1}_{\mathcal{B}_k}, \mathcal{B}_k \stackrel{\text{def}}{=} \mathcal{B}(x_k, R_k), R_k \ge \frac{R_0}{2}. \qquad (3.15)
$$
Then, for all $1 \le p < \infty$,

$\chi_k \to \chi$ weakly-* in $L^\infty((0,T) \times \mathbb{R}^3)$, strongly in $C([0,T], L^p(\mathbb{R}^3))$,

and $(u, \chi)$ is solution of

$$ \partial_t\chi + \tilde{\Pi}u \cdot \nabla\chi = 0 \text{ in } (0,T) \times \mathbb{R}^3, \quad \chi|_{t=0} = 1_{B_0} \text{ in } \mathbb{R}^3. $$

Moreover, there exist $x_b, R_b \in C^1([0,T])$ such that

$$ \chi = 1_{B}, \quad \mathcal{B} = \mathcal{B}(x_b, R_b), \quad R_b \ge \frac{R_0}{2}, $$

and the following additional convergences hold:

$$ \begin{aligned} \tilde{\Pi}_k u_k &\to \tilde{\Pi}u \text{ in } C([0, T]; S), \\ \tilde{\eta}_k &\to \tilde{\eta} \text{ in } C^1([0, T]; \mathcal{M}), \end{aligned} $$

where $\tilde{\eta}$ is the propagator associated with $\tilde{\Pi}u$.

**Proof of Proposition 3.2.** This proposition is equivalent to the stability of the family of ODEs

$$ \dot{X}_k = f_k(t, X_k) \stackrel{\text{def}}{=} (f_{x,k}(t, X_k), f_{R,k}(t, X_k)) \text{ with } X_k(0) = (x_0, R_0), $$

where for $t \in [0, T]$ and $X \in \mathbb{R}^3 \times [R_0/2, \infty)$,

$$ \left\{ \begin{array}{ll} f_{x,k}(t, X) & \stackrel{\text{def}}{=} \frac{3}{4\pi R_0^3} \int_{B_0} 1_{\Omega}(\tilde{\eta}_X[t](y))u_k(t, \tilde{\eta}_X[t](y))dy, \\ f_{R,k}(t, X) & \stackrel{\text{def}}{=} \frac{5}{4\pi R_0^4} \int_{B_0} 1_{\Omega}(\tilde{\eta}_X[t](y))(y - x_0) \cdot u_k(t, \tilde{\eta}_X[t](y))dy. \end{array} \right. \qquad (3.16) $$

and $\tilde{\eta}_X$ is defined in (3.7). From (3.8), (3.9) and (3.10) we obtain, for all $k \in \mathbb{N}$, $t \in [0, T]$, $X_1, X_2 \in \mathbb{R}^3 \times [R_0/2, \infty)$,

$$ |f_k(t, X_1) - f_k(t, X_2)| \lesssim \|\nabla u_k\|_{L^\infty((0,T) \times \Omega)} |X_1 - X_2|, $$

where the hidden constant is independent of $k$, $t$ and $X_1, X_2$. By (3.12), the sequence $\{u_k\}$ is uniformly bounded in $W^{1,\infty}((0,T) \times \Omega)$. Consequently, the $f_k$ are uniformly Lipschitz with respect to $X$ in $[R_0/2, \infty)$. Moreover for $t \in [0, T]$ and $X \in \mathbb{R}^3 \times [R_0/2, \infty)$,

$$ |f_k(t, X) - f(t, X)| \lesssim \|u_k(t, \cdot) - u(t, \cdot)\|_{L^\infty(\Omega)} $$

We deduce that $f_k \to f$ uniformly on $[0, T] \times \mathbb{R}^3 \times [R_0/2, \infty)$. The proof then follows by direct application of standard stability results on ODEs. $\square$

**Proposition 3.3** (Weak sequential continuity). Let $\{u_k, \chi_k\} \in L^2(0, T; H^1(\Omega)) \times L^\infty((0, T) \times \mathbb{R}^3)$ such that

$$ \star u_k \to u \text{ weakly in } L^2(0, T; H^1(\Omega)), \qquad (3.17) $$

$$ \star \partial_t \chi_k + \tilde{\Pi}_k u_k \cdot \nabla \chi_k = 0 \text{ in } \mathcal{D}((0, T) \times \mathbb{R}^3), \quad \chi_k|_{t=0} = 1_{B_0} \text{ in } \mathbb{R}^3, \quad B_0 = \mathcal{B}(x_0, R_0), \quad x_0 \in \mathbb{R}^3, \quad R_0 \in (0, \infty), \qquad (3.18) $$

$$ \star \chi_k(t, \cdot) = 1_{B_k(t)} \text{ with } \mathcal{B}_k(t) = \mathcal{B}(x_k(t), R_k(t)) \text{ and } R_k(t) \ge \frac{R_0}{2} \quad \forall t \in [0, T]. \qquad (3.19) $$

Then, up to a subsequence,

$$ \chi_k \to \chi \text{ weakly-} * \text{ in } L^\infty((0, T) \times \mathbb{R}^3), \text{ strongly in } C([0, T], L_{\text{loc}}^p(\mathbb{R}^3)) \quad \forall 1 \le p < \infty, $$

and $(u, \chi)$ is solution of

$$ \partial_t \chi + \tilde{\Pi}u \cdot \nabla \chi = 0 \text{ in } (0, T) \times \mathbb{R}^3, \quad \chi|_{t=0} = 1_{B_0}. $$

Moreover, there exist $x_b, R_b \in H^1([0, T])$ such that

$$ \chi(t, \cdot) = 1_{B(t)}, \quad \mathcal{B}(t) = \mathcal{B}(x_b(t), R_b(t)), \quad R(t) \ge \frac{R_0}{2} \quad \forall t \in [0, T], $$

and the following additional convergences hold:

$$ \begin{aligned} \tilde{\Pi}_k u_k &\to \tilde{\Pi}u \text{ weakly in } L^2(0, T; S), \\ \tilde{\eta}_k &\to \tilde{\eta} \text{ weakly in } H^1([0, T]; \mathcal{M}), \text{ strongly in } C([0, T], \mathcal{M}). \end{aligned} $$

where $\tilde{\eta}$ is the propagator associated with $\tilde{\Pi}u$.
Proof of proposition 3.3. Since the sequence {$u_k$} is bounded in $L^2(0, T; H^1(\Omega))$, the sequence {$\tilde{\Pi}_k u_k$} is bounded in $L^2(0, T; S)$ and up to a subsequence we have

$$ \tilde{\Pi}_k u_k \to \overline{\tilde{\Pi}u} \text{ weakly in } L^2(0, T; S). $$

Since $\chi_k(0)$ converges strongly in $L^1(\mathbb{R}^3) (= \mathbb{1}_{B_0}$ for all $k \in \mathbb{N})$, applying Di Perna-Lions theory, we obtain that $\chi_k$ converges strongly in $C([0, T]; L_{\text{loc}}^p(\mathbb{R}^3))$ for all $1 \le p < \infty$ to $\bar{\chi}$. It follows from (3.17) that, up to a subsequence $\tilde{\Pi}_k u_k \to \overline{\tilde{\Pi}u}$ weakly in $L^2(0, T; S)$ where $\overline{\tilde{\Pi}}$ is defined as in (2.21) replacing $\chi_n$ by $\bar{\chi}$. In particular, $(\bar{\chi}, \overline{\tilde{\Pi}u})$ is a solution to our transport equation. Eventually, the convergence of the propagators $\tilde{\eta}_k$ from the convergence of $\chi_k$ follows easily. $\square$

A straightforward consequence of Proposition 3.3 is the following corollary valid for sequences {$u_k$} satisfying $\tilde{\Pi}_k u_k = u_k$.

**Corollary 3.1.** Let {$u_k, \eta_k$} a bounded sequence of $L^2(0, T; H^1(\Omega)$ and {$B_k$} given for all $t \in [0, T]$ by $B_k(t) = \eta_k[t](B_0)$ where $B_0$ is ball of radius $R_0 \in (0, \infty)$ and center $x_0 \in \mathbb{R}^3$. Assume that for each $k \in \mathbb{N}$, $u_k$ is compatible with the system {$B_k, \eta_k$} and satisfy

$$ u_k \to u \text{ weakly in } L^2(0, T; H^1(\Omega)). $$

Then, passing to a subsequence as the case may be, we have

$$ \eta_k \to \eta \text{ weakly in } H^1([0, T]; \mathcal{M}), \text{ strongly in } C([0, T]; \mathcal{M}). $$

Moreover, defining $B(t)$ for all $t \in [0, T]$ by

$$ B(t) \stackrel{\text{def}}{=} \eta[t](B_0), $$

the velocity field $u$ is compatible with the system {$B, \eta$}.

**Proof.** This proposition is a direct consequence of Proposition 3.3 taking $\chi_k = \mathbb{1}_{B_k}$. $\square$

The last proposition deals with how densities propagate along the characteristics of the bubble motion. Its formulation and its proof is inspired of [14, Lemme 3.2] in which the case of rigid solids is treated.

**Proposition 3.4.** Let $\gamma > 1$, $\rho \in L^\infty(0, T; L^\gamma(\Omega))$ and $u \in L^2(0, T; H^1(\Omega))$ such that $(\rho, u)$ satisfy the continuity equation (1.21) for all $\varphi \in D((0, T) \times \mathbb{R}^3)$. Let us suppose moreover that $u$ is compatible with the system {$B, \eta$} and there exists $\rho_{b,0} \in (0, \infty)$ such that

$$ \rho_0 = \rho_{b,0} \text{ in } B_0. $$

Then

$$ |B(t)|\rho(t, \eta[t](x)) = |B_0|\rho_{b,0} \text{ for all } t \in [0, T] \text{ and a.e. } x \in B_0. $$

**Proof.** Following the same arguments as in [14, Lemme 3.2], we obtain

$$ \int_{\mathcal{U}} \rho(t, \eta[t](y)) \frac{|B(t)|}{|B_0|} dy = \int_{\mathcal{U}} \rho_{b,0} dy $$

for any open ball $\mathcal{U}$. It follows that

$$ |B(t)|\rho(t, \eta[t](x)) = |B_0|\rho_{b,0} \text{ for all } t \in [0, T] \text{ and a.e. } x \in B_0. $$

$\square$

# 4 Existence proofs of the approximate solutions

In this section, we successively prove Proposition 2.4, 2.3, 2.2 and 2.1. The main contribution with respect to fluid-solid interaction is the convergence analysis of the nonlinear terms arising in the bubble domain, specifically the convective term $\rho u \otimes u : \text{div}(\varphi)$ and the pressure term $p(\rho, \mathbb{1}_B) \text{div}(\varphi)$ appearing in the momentum equation (1.22). These terms vanish naturally for rigid-body motion as $\text{div}(\varphi) = 0$ in the solid domain.
## 4.1 Existence proof of the Faedo-Galerkin approximation

In this subsection, we construct a solution $(B_N, \rho_N, u_N)$ to the problem (2.28)-(2.30). The strategy consists in reformulating the Galerkin problem as a fixed point problem. More precisely, we will look for $u_N$ as a fixed point of the operator

$$ \mathcal{F}_N : \begin{cases} B_{Q,T} \mapsto B_{Q,T}, \\ u \mapsto \tilde{u} \end{cases}, $$

where

$$ B_{Q,T} \stackrel{\text{def}}{=} \{u \in \mathcal{C}([0,T]; X_N), \|u\|_{L^\infty(0,T;L^2(\Omega))} \le Q\}, $$

for some constants $Q, T \in (0, \infty)$. This approach first yields a solution on a time interval $[0, T]$, with $T$ possibly depending on $N$. We will then prove that we can extend this solution to an arbitrary time $\tilde{T}$ independent of $N$. To prove the existence of a fixed point for $\mathcal{F}_N$, we will apply Schauder's fixed point theorem (for the proof see e.g. [11]).

**Theorem 4.1** (Schauder's fixed point theorem). Let $X$ be a real Banach space and $K \subset X$ a closed, bounded and convex subset. Assume that

$$ \mathcal{F}: K \to K $$

is continuous and compact. Then $\mathcal{F}$ has a fixed point in $K$.

We now prove that for a certain $(Q, T)$, $\mathcal{F}_N$ satisfies the hypotheses of Theorem 4.1. According to Proposition 3.1, if $Q$ and $T$ satisfy the condition

$$ T \le \sqrt{\frac{\pi R_0^2 \sqrt{R_0}}{5 Q}}, \qquad (4.1) $$

then for all $u \in B_{Q,T}$, there exists a unique solution to

$$ \partial_t \chi + \tilde{\Pi} u \cdot \nabla \chi = 0 \text{ in } \mathcal{D}'((0,T) \times \mathbb{R}^3), \quad \chi|_{t=0} = \mathbb{1}_{B_0} \text{ in } \mathbb{R}^3, \qquad (4.2) $$

which satisfies, for all $1 \le p < \infty$,

$$ \chi \in \mathcal{C}([0, T]; L_{\text{loc}}^p(\mathbb{R}^3)). $$

Moreover, there exist $x_b, R_b \in \mathcal{C}^1([0, T])$ such that, for all $t \in [0, T]$,

$$ \chi(t, \cdot) = \mathbb{1}_{B(x_b(t), R_b(t))} \quad \text{and} \quad R_b(t) \ge R_0/2. \qquad (4.3) $$

Accordingly, we define

$$ \mu \stackrel{\text{def}}{=} (1-\chi)\mu_f + n\chi, \quad \nu \stackrel{\text{def}}{=} (1-\chi)\nu_f + \chi\nu_b, \quad p_\delta(\rho, \chi) \stackrel{\text{def}}{=} (1-\chi)a_f\rho^{\gamma_f} + \chi a_b\rho^{\gamma_b} + \delta\rho^\beta. $$

From now on, we assume that $Q$ and $T$ satisfies (4.1) and we fix $u \in B_{Q,T}$. Let $\rho$ be the solution to

$$ \partial_t \rho + \text{div}(\rho u) = \varepsilon \Delta \rho \text{ in } (0,T) \times \Omega, \quad \nabla \rho \cdot n = 0 \text{ on } \partial \Omega, \quad \rho|_{t=0} = \rho_0 \text{ in } \Omega, \qquad (4.4) $$

with $0 < \underline{\rho} \le \rho_0 \le \bar{\rho}$. We recall the following classical maximal regularity result for the parabolic problem, (see, e.g. [11]):

**Proposition 4.1.** Let $\Omega$ be a bounded smooth domain, $\underline{\rho}, \bar{\rho} \in (0, \infty)$ and $u \in L^\infty(0, T; W^{1,\infty}(\Omega))$. Assume that $\rho_0 \in W^{1,\infty}(\Omega)$ satisfies $\underline{\rho} \le \rho_0 \le \bar{\rho}$. Then the parabolic problem (4.4) admits a unique solution, which satisfies

$$ \rho \in L^2(0, T; H^2(\Omega)) \cap \mathcal{C}([0, T]; H^1(\Omega)) \cap H^1(0, T; L^2(\Omega)) \qquad (4.5) $$

and

$$ \underline{\rho} \exp \left( - \int_0^t \| \text{div}(u(\tau)) \|_{L^\infty(\Omega)} d\tau \right) \le \rho(t, x) \le \bar{\rho} \exp \left( \int_0^t \| \text{div}(u(\tau)) \|_{L^\infty(\Omega)} d\tau \right) \text{ for all } t \in [0, T] \text{ and a.e. } x \in \Omega. \qquad (4.6) $$

Since $\rho_0 \in W^{1,\infty}(\Omega)$ and $u \in B_{Q,T}$ in (4.4), we may apply Proposition 4.1. We conclude that $\rho$ satisfies (4.5) as well as the estimate (4.6). Since all norms are equivalent on the finite-dimensional space $X_N$, there exists a constant $c(N)$, depending only on $N$, such that for all $\varphi \in X_N$,

$$ c(N)^{-1} \|\varphi\|_{W^{1,\infty}(\Omega)} \le \|\varphi\|_{L^2(\Omega)} \le c(N) \|\varphi\|_{W^{1,\infty}(\Omega)}. $$
In particular, we deduce that

$$ \underline{\rho} \exp(-c(N)QT) \le \rho(t, x) \le \bar{\rho} \exp(c(N)QT) \quad \text{for all } t \in [0, T], \text{ and a.e. } x \in \Omega. \quad (4.7) $$

For $\tilde{u} : [0, T] \to X_N$, we may write

$$ \tilde{u}(t) \stackrel{\text{def}}{=} \sum_{i=1}^{N} g_i(t) \psi_i $$

where the functions $\psi_i$ are defined in (2.27). We then consider the following system of ODEs:

$$ A(t) \frac{d}{dt} \tilde{u}(t) + B(t) \tilde{u}(t) = F(t), \quad \tilde{u}(0) = u_0 \stackrel{\text{def}}{=} \sum_{i=1}^{N} \left( \int_{\Omega} u_0 \cdot \psi_i dy \right) \psi_i, \quad (4.8) $$

where the matrices $A(t) \stackrel{\text{def}}{=} (a_{i,j}(t))_{1 \le i,j \le N}$, $B(t) \stackrel{\text{def}}{=} (b_{i,j}(t))_{1 \le i,j \le N}$ and $F(t) \stackrel{\text{def}}{=} (f_i(t))_{1 \le i \le N}$ are defined by

$$ a_{i,j}(t) \stackrel{\text{def}}{=} \int_{\Omega} \rho \psi_i \cdot \psi_j dy, \quad (4.9) $$

$$ b_{i,j}(t) \stackrel{\text{def}}{=} \int_{\Omega} \rho(u \cdot \nabla \psi_j) \cdot \psi_i dy + \int_{\Omega} 2\mu \left( \mathbb{D}(\psi_i) - \frac{1}{3} \operatorname{div}(\psi_i) \right) \cdot \left( \mathbb{D}(\psi_j) - \frac{1}{3} \operatorname{div}(\psi_j) \right) + \nu \operatorname{div}(\psi_i) \operatorname{div}(\psi_j) dy \\ + \int_{\Omega} \varepsilon (\nabla \psi_j \nabla \rho) \cdot \psi_i dy + n \int_{\Omega} \chi(\psi_i - \Pi(t)\psi_i) \cdot (\psi_j - \Pi(t)\psi_j) dy, \quad (4.10) $$

$$ f_j(t) \stackrel{\text{def}}{=} - \int_{\Omega} \rho g \cdot \psi_j dy + \int_{\Omega} \left( \chi \frac{\kappa_b}{R_b} + p_{\delta}(\rho, \chi) \right) \operatorname{div}(\psi_j) dy, \quad (4.11) $$

and the projector $\Pi$ is defined as in (2.17), with $\chi_n$ replaced with $\chi$. The positive lower bound of $\rho$ in (4.7) implies $A(t) \ge \rho \mathbb{I}_N$ in the sense of symmetric positive matrices. In particular, the matrix $A$ is invertible. Using the regularity of $\chi$ (Proposition 3.1) and $\rho$ (Proposition 4.1), we deduce that $A, B$ and $F$ are continuous on $[0, T]$. By the Cauchy-Lipschitz theorem, the system (4.8) admits a unique solution

$$ \tilde{u} \in C^1([0, T]; X_N). $$

We now prove that there exist $T, Q \in (0, \infty)$ such that $\mathcal{F}_N$ defines a continuous and compact mapping from $B_{Q,T}$ into itself. We first show that, for suitable choice of $Q$ and $T$, $B_{Q,T}$ is stable under $\mathcal{F}_N$. Using the continuity equation (4.4), we obtain the following identities by integration by parts:

$$ \int_0^T \int_{\Omega} \rho \frac{d}{dt} \tilde{u} \cdot \tilde{u} dy d\tau = -\frac{1}{2} \int_0^T \int_{\Omega} \partial_t \rho |\tilde{u}|^2 dy d\tau + \frac{1}{2} \int_{\Omega} (\rho |\tilde{u}|^2)(t) dy - \frac{1}{2} \int_{\Omega} \rho_0 |u_0|^2 dy, \quad (4.12) $$

$$ \int_{\Omega} \rho(u \cdot \nabla) \tilde{u} \cdot \tilde{u} dy = -\frac{1}{2} \int_{\Omega} \operatorname{div}(\rho u) |\tilde{u}|^2 dy \quad (4.13) $$

Multiplying (4.8) by $\tilde{u}$, integrating in time, and combining with (4.12)-(4.13), we obtain the following energy identity:

$$ \int_{\Omega} \frac{1}{2} \rho |\tilde{u}|^2 dy + \int_{0}^{T} \int_{\Omega} 2\mu |\mathbb{D}(\tilde{u}) - \frac{1}{3} \operatorname{div}(\tilde{u}) \mathbb{I}_3|^2 + \nu |\operatorname{div}(\tilde{u})|^2 dy d\tau \\ = - \int_{0}^{T} \int_{\Omega} \rho g \cdot \tilde{u} dy d\tau + \int_{0}^{T} \int_{\Omega} \left( \chi \frac{\kappa_b}{R_b} + p_{\delta}(\rho, \chi) \right) \operatorname{div}(\tilde{u}) dy d\tau + \frac{1}{2} \int_{\Omega} \rho_0 |u_0|^2 dy \quad (4.14) $$

Using Hölder's inequality, we bound the right-hand side of (4.14) by

$$ T \|\rho\|_{L^{\infty}((0,T)\times\Omega)} \|g\|_{L^{\infty}(0,T;L^2(\Omega))} \|\tilde{u}\|_{L^{\infty}(0,T;L^2(\Omega))} \\ + \left( \|\chi(\kappa_b/R_b)\|_{L^1((0,T)\times\Omega)} + \|p(\rho, \chi)\|_{L^1((0,T)\times\Omega)} \right) \|\operatorname{div}(\tilde{u})\|_{L^{\infty}((0,T)\times\Omega)} + \frac{1}{2} \|\sqrt{\rho_0} u_0\|_{L^2(\Omega)}^2. $$

Using the lower bound on $R_b$ from (4.3) and the bounds on $\rho$ from (4.7), we get

$$ \|\chi(\kappa_b/R_b)\|_{L^1((0,T)\times\Omega)} \le (2\kappa_b/R_0)|\Omega|T. $$

and

$$ \|p_{\delta}(\rho, \chi)\|_{L^1((0,T)\times\Omega)} \le a_f \|\rho^{\gamma_f}\|_{L^1((0,T)\times\Omega)} + a_b \|\rho^{\gamma_b}\|_{L^1((0,T)\times\Omega)} + \delta \|\rho^{\beta}\|_{L^1((0,T)\times\Omega)}, \\ \le (a_f \bar{\rho}^{\gamma_f} \exp(\gamma_f c(N)QT) + a_b \bar{\rho}^{\gamma_b} \exp(\gamma_b c(N)QT) + \delta \bar{\rho}^{\beta} \exp(\beta c(N)QT)) |\Omega|T. $$
Combining these bounds and using (4.7), we deduce

$$
\begin{aligned}
\|\tilde{u}\|_{L^{\infty}(0,T;L^2(\Omega))}^2 &\le 2 \frac{\bar{\rho}}{\rho} \exp(2c(N)QT) \|g\|_{L^{\infty}(0,T;L^2(\Omega))} QT \\
&\quad + 2c(N) \frac{|\Omega|}{\rho} \exp(c(N)QT) \left( \frac{2\kappa_b}{R_0} + a_f \bar{\rho}^{\gamma_f} \exp(\gamma_f c(N)QT) + a_b \bar{\rho}^{\gamma_b} \exp(\gamma_b c(N)QT) + \delta \bar{\rho}^\beta \exp(\beta c(N)QT) \right) QT \\
&\quad + \frac{1}{\rho} \exp(c(N)QT) \|\sqrt{\rho_0} u_0\|_{L^2(\Omega)}^2.
\end{aligned}
$$

We now choose

$$
Q \geq \max \left\{ 1, \frac{4}{\rho} \| \sqrt{\rho_0} u_0 \|_{L^2(\Omega)}^2 \right\}, \quad (4.15)
$$

and

$$
T \stackrel{\text{def}}{=} T(N) = \min \left\{ \frac{\log(2)}{2Q}, \frac{\underline{\rho}Q}{16\bar{\rho}\|g\|_{L^\infty(0,T;L^2(\Omega))}}, \frac{\log(2)}{\bar{\gamma}c(N)Q}, \frac{\underline{\rho}Q}{32c(N)|\Omega|((\kappa_b/R_0) + a_f\bar{\rho}^{\gamma_f} + a_b\bar{\rho}^{\gamma_b} + \delta\bar{\rho}^\beta)} \right\}, \quad (4.16)
$$

where $\bar{\gamma} = \max\{\gamma_f, \gamma_b, \beta\}$, and therefore

$$
\mathcal{F}_N(B_{Q,T}) \subset B_{Q,T}.
$$

We now want to prove that $\mathcal{F}_N$ is continuous on $B_{Q,T}$. Let $\{u_k\} \subset B_{Q,T}$ such that $u_k \to u$ in $C([0,T]; D(\bar{\Omega}))$ and $\chi_k$ satisfying (4.2). According to Proposition 3.2, we have

$$
\begin{gather*}
\chi_k \to \chi \text{ weakly-}^* \text{ in } L^\infty((0,T) \times \mathbb{R}^3) \text{ and strongly in } C([0,T]; L^p(\mathbb{R}^3)) \quad \forall 1 \le p < \infty, \\
\tilde{\eta}_k \to \tilde{\eta} \text{ strongly in } C^1([0,T]; \mathcal{M}).
\end{gather*}
$$

We deduce that

$$
a_{i,j}^{k} \rightarrow a_{i,j}, \quad b_{i,j}^{k} \rightarrow b_{i,j}, \quad f_{j}^{k} \rightarrow f_{j} \text{ in } C([0,T]),
$$

and so we obtain, from standard stability results on ODE's,

$$
\mathcal{F}_N(u_k) = \tilde{u}_k \rightarrow \tilde{u} = \mathcal{F}_N(u) \text{ in } C([0, T]; X_N).
$$

We now prove that $\mathcal{F}_N : B_{Q,T} \to B_{Q,T}$ is a compact operator. Let $\tilde{u}$ be the solution to (4.8). Then, for all $t \in [0, T]$,

$$
\left| \frac{d}{dt} \tilde{u}(t) \right| \leq \left| \mathbb{A}^{-1}(t) \left| (\mathbb{B}(t)) \|\tilde{u}(t)\| + |F(t)| \right| \right| \leq \left| \mathbb{A}^{-1}(t) \left| (Q|\mathbb{B}(t)|) + |F(t)| \right| \right|,
$$

so that

$$
\sup_{t \in [0, T]} \left\{ |\tilde{u}(t)| + \left| \frac{d}{dt} \tilde{u}(t) \right| \right\} \leq C(N),
$$

for a constant $c(N)$ depending only on $N$. Consequently,

$$
\sup_{u \in B_{Q,T}} \| \mathcal{F}_N(u) \|_{C^1([0,T]; X_N)} \le C(N)
$$

By the Arzelà-Ascoli theorem, it follows that $\mathcal{F}_N : B_{Q,T} \to B_{Q,T}$ is compact. We can then apply Schauder's fixed point theorem 4.1 to conclude the existence of a fixed point $u_N \in B_{Q,T}$. At this stage, we have thus established the existence of a solution $u_N$ of the finite dimensional problem on an interval $[0, T]$ with $T$ possibly depending on $N$. We denote by $\rho_N$ and $\chi_N$ the corresponding solutions of the regularized continuity and transport equations on $[0, T] \times \mathbb{R}^3$. We also define $\tilde{\eta}_N$ the propagator associated with the velocity field $\tilde{\Pi}_N u$ given explicitly by

$$
\tilde{\eta}_N[t](x) \stackrel{\text{def}}{=} x_N(t) + \frac{R_N(t)}{R_0}(x - x_0) \quad \forall (t, x) \in [0, T] \times \mathbb{R}^3. \tag{4.17}
$$

To extend this solution to an arbitrary time $\tilde{T} > 0$ independent of $N$, we use a standard continuation argument based on energy estimates uniform with respect to $T(N)$. To obtain such uniform bounds on $u_N$, we multiply the momentum equation (2.18) by $u_N$ and derive the following energy identity, valid on $[0, T(N)]$:

$$
\begin{align}
& \frac{d}{dt} \tilde{E}_{\delta}(\rho_N, \rho_N u_N) + \int_{\Omega} \left( 2\mu_N |\mathbb{D}(u_N) - \frac{1}{3} \operatorname{div}(u_N)|^2 + \nu_N |\operatorname{div}(u_N)|^2 \right) \mathrm{d}y \mathrm{d}\tau \nonumber \\
& \qquad + \delta \varepsilon \beta \int_{\Omega} (\rho_N)^{\beta-2} |\nabla \rho_N|^2 \mathrm{d}y \mathrm{d}\tau \nonumber \\
&= a_f \int_{\Omega} (1 - \chi_N) \rho_N^{\gamma_f} \operatorname{div}(u_N) \mathrm{d}y + a_b \int_{\Omega} \chi_N \rho_N^{\gamma_b} \operatorname{div}(u_N) \mathrm{d}y - \int_{\Omega} \rho_N g \cdot u_N \mathrm{d}y \mathrm{d}\tau + \int_{\Omega} \chi_N \frac{\kappa_b}{R_N} \operatorname{div}(u_N) \mathrm{d}y \mathrm{d}\tau, \tag{4.18}
\end{align}
$$
where

$$ \tilde{E}_{\delta}(\rho_N, \rho_N u_N) \stackrel{\text{def}}{=} \int_{\Omega} \left( \frac{1}{2} \rho_N |u_N|^2 + \delta \frac{\rho_N^{\beta}}{\beta - 1} \right) \mathrm{d}y. $$

Since $\beta \ge \max\{2\gamma_f, 2\gamma_b\}$, the terms on the right handside of (4.18) can be estimated using Hölder's inequality and the embeddings

$$ L^{\beta}(\Omega) \hookrightarrow L^{2\gamma_f}(\Omega), \quad L^{\beta}(\Omega) \hookrightarrow L^{2\gamma_b}(\Omega). $$

More precisely, they can be bounded by

$$ a_f |\Omega|^{(\beta-2\gamma_f)/2\beta} \|\rho_N\|_{L^\beta(\Omega)}^{\gamma_f} \|\mathrm{div}(u_N)\|_{L^2(\Omega)} + a_b |\Omega|^{(\beta-2\gamma_b)/2\beta} \|\rho_N\|_{L^\beta(\Omega)}^{\gamma_b} \|\mathrm{div}(u_N)\|_{L^2(\Omega)} \\ + \|\rho_N|u_N|^2\|_{L^1(\Omega)}^{\frac{1}{2}} \|\rho_N\|_{L^\beta(\Omega)}^{\frac{1}{2}} \|g\|_{L^{\frac{2\beta}{\beta-1}}(\Omega)}^{\frac{2\beta}{\beta-1}} + \kappa_b \|\mathrm{div}(u_N)\|_{L^2(\Omega)} \|\chi_N/R_N\|_{L^2(\Omega)} $$

Applying Young's inequality repeatedly, it follows that

$$ \frac{a_f^2}{4\theta} |\Omega|^{(\beta-2\gamma_f)/\beta} \|\rho_N\|_{L^\beta(\Omega)}^{2\gamma_f} + \theta \|\mathrm{div}(u_N)\|_{L^2(\Omega)}^2 + \frac{a_b^2}{4\theta} |\Omega|^{(\beta-2\gamma_b)/\beta} \|\rho_N\|_{L^\beta(\Omega)}^{2\gamma_b} + \theta \|\mathrm{div}(u_N)\|_{L^2(\Omega)}^2 \\ \le a_f^{\frac{2\beta}{\beta-2\gamma_f}} \frac{\beta-2\gamma_f}{4\theta\beta} \left( \frac{6\gamma_f(\beta-1)}{4\theta\delta\beta} \right)^{\frac{2\gamma_f}{\beta-2\gamma_f}} |\Omega| + \frac{\delta}{3(\beta-1)} \|\rho_N\|_{L^\beta(\Omega)}^\beta \\ + a_b^{\frac{2\beta}{\beta-2\gamma_b}} \frac{\beta-2\gamma_b}{4\theta\beta} \left( \frac{6\gamma_b(\beta-1)}{4\theta\delta\beta} \right)^{\frac{2\gamma_b}{\beta-2\gamma_b}} |\Omega| + \frac{\delta}{3(\beta-1)} \|\rho_N\|_{L^\beta(\Omega)}^\beta + 2\theta \|\mathrm{div}(u_N)\|_{L^2(\Omega)}^2 $$

and

$$ \frac{1}{2} \|\rho|u_N|^2\|_{L^1(\Omega)} + \frac{1}{2} \|\rho_N\|_{L^\beta(\Omega)} \|g\|_{L^{\frac{2\beta}{\beta-1}}(\Omega)}^2 + \theta \|\mathrm{div}(u_N)\|_{L^2(\Omega)}^2 + \frac{\kappa_b}{4\theta} \|(\chi_N/R_N)\|_{L^2(\Omega)}^2 \\ \le \frac{1}{2} \|\rho|u_N|^2\|_{L^1(\Omega)} + \frac{\delta}{3(\beta-1)} \|\rho_N\|_{L^\beta(\Omega)}^\beta + \frac{\beta-1}{2\beta} \left( \frac{3(\beta-1)}{2\delta\beta} \right)^{\frac{1}{\beta-1}} \|g\|_{L^{\frac{2\beta}{\beta-1}}(\Omega)}^{\frac{2\beta}{\beta-1}} + \frac{\kappa_b^2}{4\theta} \|\chi_N/R_N\|_{L^2(\Omega)}^2 + \theta \|\mathrm{div}(u_N)\|_{L^2(\Omega)}^2, $$

where $\theta > 0$ is a positive constant to be fixed later. Finally, since $R_N \ge R_0/2$ on $[0, T(N)]$, we have the explicit bound

$$ \frac{\kappa_b^2}{4\theta} \|\chi_N/R_{b,N}\|_{L^2(\Omega)}^2 \le \frac{\kappa_b^2 |\Omega|}{\theta R_0^2}. $$

Setting $\mu_m = \min\{n, \mu_f\}$ and $\nu_m = \min\{\nu_f, \nu_b\}$, we have

$$ \int_{\Omega} 2\mu |\mathbb{D}(u_N) - \frac{1}{3} \mathrm{div}(u_N) \mathbb{I}_3|^2 + \nu |\mathrm{div}(u_N)|^2 \mathrm{d}y \ge \int_{\Omega} 2\mu_m |\mathbb{D}(u_N) - \frac{1}{3} \mathrm{div}(u_N) \mathbb{I}_3|^2 + \nu_m |\mathrm{div}(u_N)|^2 \mathrm{d}y = \int_{\Omega} \mu_m |\nabla u_N|^2 + \left(\frac{\mu_m}{3} + \nu_m\right) |\mathrm{div}(u_N)|^2 \mathrm{d}y $$

and choosing suitable $n$ and $\theta$ such that

$$ \mu_m = \mu_f, \quad \theta \le \frac{1}{6} \left( \nu_m + \frac{\mu_f}{3} \right), $$

we obtain

$$ \frac{\mathrm{d}}{\mathrm{d}t} \tilde{E}_{\delta}(\rho_N, \rho_N u_N) + \mu_f \|\nabla u_N\|_{L^2(\Omega)}^2 \le \tilde{E}_{\delta}(\rho_N, \rho_N u_N) + c(\beta, \delta) \|g\|_{L^{\frac{2\beta}{\beta-1}}(\Omega)}^{\frac{2\beta}{\beta-1}} + c(\gamma_f, \gamma_b, \beta, \delta, \mu_f, \nu_m) |\Omega| + c(\kappa_b, \mu_f, \nu_m) \frac{|\Omega|}{R_0^2}. $$

Consequently, by Gronwall's lemma, for $t \in [0, T(N)]$, we have

$$ E_N(\rho_N(t), \rho_N u_N(t)) + \mu_f \int_0^t \|\nabla u_N\|_{L^2(\Omega)}^2 \mathrm{d}\tau \\ \le e^t \tilde{E}_\delta(\rho_0, q_0) + \int_0^t \left( c(\beta, \delta) \|g\|_{L^{\frac{2\beta}{\beta-1}}(\Omega)}^{\frac{2\beta}{\beta-1}} + c(\gamma_f, \gamma_b, \beta, \delta, \mu_f, \nu_m) |\Omega| + c(\kappa_b, \mu_f, \nu_m) \frac{|\Omega|}{R_0^2} \right) e^{t-\tau} \mathrm{d}\tau \quad (4.19) $$
A direct consequence of (4.19) is the uniform bound

$$ \|\nabla u_N\|_{L^2((0,T(N))\times\Omega)}^2 \le e^{\tilde{T}} \left( \tilde{E}_\delta(\rho_0, q_0) + c(\beta, \delta) \|g\|_{L^{\frac{2\beta}{\beta-1}}((0,\tilde{T})\times\Omega)}^{\frac{2\beta}{\beta-1}} + c(\gamma_f, \gamma_b, \beta, \delta, \mu_f, \nu_m) \tilde{T}|\Omega| + c(\kappa_b, \mu_f, \mu_m) \frac{\tilde{T}|\Omega|}{R_0^2} \right) \stackrel{\text{def}}{=} K, \quad (4.20) $$

for any fixed $\tilde{T} > 0$ and $0 < T(N) \le \tilde{T}$. By equivalence of the $L^\infty$ and $L^2$ norms on $X_N$, it follows that

$$ 0 < \rho \exp(-d(N)K) \le \rho(t, x) \le \bar{\rho} \exp(d(N)K) \text{ for all } t \in [0, T(N)] \text{ and a.e. } x \in \Omega. \quad (4.21) $$

Furthermore, (4.19) implies

$$ \sup_{t \in [0, T(N)]} \int_{\Omega} \rho_N(t) |u_N(t)|^2 dy \le K, $$

which, together with (4.21), yields

$$ \|u_N\|_{L^\infty(0,T(N);L^2(\Omega))} \le \bar{\rho} \exp(d(N)K) K. \quad (4.22) $$

Since (4.22) is independent of $T(N)$, we can iterate the fixed-point argument with updated initial data $(B_N(T_N), \rho_N(T_N), u_N(T_N))$, provided that

$$ \inf_{t \in [0, T(N)]} R_N(T(N)) \ge \frac{R_0}{2} \quad \text{and} \quad \inf_{t \in [0, T(N)]} \text{dist}(B_N(t), \partial\Omega) \ge 2\sigma > 0. \quad (4.23) $$

The first condition in (4.23) holds if

$$ \| \dot{R}_N \|_{L^1(0, T(N))} \le \frac{R_0}{2}. $$

Let $\bar{T} > 0$, recalling that, for all $t \in [0, \bar{T}]$,

$$ \dot{R}_N(t) = \frac{5}{4\pi R_0^4} \int_{B_0} \mathbb{1}_{\Omega}(\tilde{\eta}_N[t](y))(y-x_0) \cdot u_N(t, \tilde{\eta}_N[t](y)) dy, \quad \tilde{\eta}_N[t](x) = x_N(t) + \frac{R_N(t)}{R_0}(x-x_0), $$

we deduce

$$ \| \dot{R}_N \|_{L^1(0, \bar{T})} \le \frac{\sqrt{\bar{T}}}{R_0} \sqrt{\frac{5}{4\pi R_0}} \| u_N \|_{L^2((0, \bar{T}) \times \Omega)}. $$

By Poincaré inequality and estimate (4.19),

$$ \|u_N\|_{L^2((0,\bar{T})\times\Omega)}^2 \le c_p^2 K. $$

Consequently, if

$$ T(N) \le T_1 \stackrel{\text{def}}{=} \frac{\pi R_0^5}{5 c_p^2 K}, \quad (4.24) $$

then the first condition in (4.23) is satisfied. The second condition in (4.23) holds if, for almost every $x \in B_0$,

$$ \left\| \frac{d}{dt} \tilde{\eta}_N[\cdot](x) \right\|_{L^1(0, T(N))} < \text{dist}(B_0, \partial\Omega) - 2\sigma. $$

Let $\bar{T} > 0$, for every $t \in [0, \bar{T}]$ and almost every $x \in B_0$,

$$ \left| \frac{d}{dt} \tilde{\eta}_N[t](x) \right| = \left| \dot{x}_N(t) + \frac{\dot{R}_N(t)}{R_0} (x - x_0) \right| \le \left| \dot{x}_N(t) \right| + \left| \frac{\dot{R}_N(t)}{R_0} \right| |x - x_0| \le \left| \dot{x}_N(t) \right| + \left| \dot{R}_N(t) \right|. $$

Recalling that

$$ \dot{x}_N = \frac{3}{4\pi R_0^3} \int_{B_0} \mathbb{1}_{\Omega}(\tilde{\eta}_N[t](y)) u_N(t, \tilde{\eta}_N[t](y)) dy, \quad \dot{R}_N = \frac{5}{4\pi R_0^4} \int_{B_0} \mathbb{1}_{\Omega}(\tilde{\eta}_N[t](y))(y-x_0) \cdot u_N(t, \tilde{\eta}_N[t](y)) dy, $$

we obtain

$$ \|\dot{x}_N\|_{L^2(0, \bar{T})} \le \frac{\sqrt{\bar{T}}}{R_0} \sqrt{\frac{3}{4\pi R_0}} \|u_N\|_{L^2((0, \bar{T}) \times \Omega)}, \quad \|\dot{R}_N\|_{L^2(0, \bar{T})} \le \frac{\sqrt{\bar{T}}}{R_0} \sqrt{\frac{5}{4\pi R_0}} \|u_N\|_{L^2((0, \bar{T}) \times \Omega)}. $$

Thus, for almost every $x \in B_0$,

$$ \left\| \frac{d}{dt} \tilde{\eta}_N[\cdot](x) \right\|_{L^1(0, \bar{T})} \le c_0 \sqrt{\bar{T}} \|u_N\|_{L^2((0, \bar{T}) \times \Omega)}, $$

with $c_0 = \frac{\sqrt{3}+\sqrt{5}}{R_0\sqrt{4\pi R_0}}$. It follows that if

$$ T(N) \le T_2 = \frac{(\text{dist}(B_0, \partial\Omega) - 2\sigma)^2}{c_0^2 c_p^2 K}, \quad (4.25) $$

then the second condition in (4.23) is satisfied. Finally, we choose $\tilde{T}$ in (4.20) such that $\tilde{T} = \min\{T_1, T_2\}$ and conclude the existence of a solution $u_N$ at least until time $\tilde{T}$.
## 4.2 Convergence of the Faedo-Galerkin scheme and the limiting system

Proposition 2.4 establishes the existence of a solution $(\mathcal{B}_N, \rho_N, u_N)$ to (2.28)-(2.30) on a time interval $[0, T]$ with $T$ independent of $N$. In this section, we prove Proposition 2.3 by passing to the limit in (2.28)-(2.30) as $N \to \infty$, thereby recovering a solution to the $n$-level approximation problem (2.16)-(2.19). Given initial data $(\rho_{0,n}, u_{0,n})$ for the $n$-level approximation problem, we construct initial data $(\rho_{0,N}, u_{0,N})$ for the Faedo-Galerkin scheme in such a way

$$ \rho_{0,N} \to \rho_{0,n} \text{ in } W^{1,\infty}(\Omega), \quad \rho_{0,N}u_{0,N} \to q_{0,n} \text{ in } L^2(\Omega), $$

and

$$ \int_{\Omega} \left( \frac{1}{2} \rho_{0,N} |u_{0,N}|^2 + \delta \frac{\rho_N^{\beta}}{\beta-1} dy \right) \to \int_{\Omega} \left( \frac{1}{2} \frac{|q_{0,n}|^2}{\rho_{0,n}} + \delta \frac{\rho_N^{\beta}}{\beta-1} \right) dy. \quad (4.26) $$

The energy estimate (2.32) yields, up to a subsequence,

$$ u_N \to u_n \text{ weakly in } L^2(0, T; H^1(\Omega)), \quad \rho_N \to \rho_n \text{ weakly-}^* \text{ in } L^\infty(0, T; L^\beta(\Omega)), \quad \nabla \rho_N \to \nabla \rho_n \text{ weakly in } L^2((0, T) \times \Omega). \quad (4.27) $$

It follows from Proposition 3.3 that, up to a subsequence and for all $1 \le p < \infty$,

$$ \chi_N \to \chi_n \text{ weakly-}^* \text{ in } L^\infty((0, T) \times \mathbb{R}^3), \text{ strongly in } C([0, T]; L_{\text{loc}}^p(\mathbb{R}^3)), \quad (4.28) $$

where the pair $(\chi_n, u_n)$ is a solution of (2.20). Moreover, there exist functions $x_n, R_n \in H^1(0, T)$ such that, for all $t \in [0, T]$,

$$ \chi_n(t, \cdot) = \mathbb{1}_{\mathcal{B}(x_n(t), R_n(t))} \text{ and } R_n(t) \ge R_0/2. $$

Finally, using the strong convergence of $\chi_N$ given in (4.28) and the weak convergence of $u_N$ given in (4.27), it follows that

$$ \tilde{\Pi}_N u_N \to \tilde{\Pi}_n u_n \text{ weakly in } L^2(0, T; S). \quad (4.29) $$

Following the same arguments as in [17, Section 2.4], [33, Section 7.8.1], we further get, for all $1 \le p < \frac{4}{3}\beta$,

$$ \rho_N \to \rho_n \text{ in } C([0, T]; L_{\text{weak}}^{\beta}(\Omega)), \text{ strongly in } L^p((0, T) \times \Omega), $$

$$ \rho_N u_N \to \rho_n u_n \text{ weakly in } L^2(0, T; L_{\beta+6}^{\frac{6\beta}{\beta+6}}(\Omega)), \text{ weakly-}^* \text{ in } L^\infty(0, T; L_{\beta+1}^{\frac{2\beta}{\beta+1}}(\Omega)). $$

These convergence results allow us to pass to the limit $N \to \infty$ in the regularized continuity equation (2.9). We now consider the limit in the momentum equation (2.18). The delicate terms are, for $1 \le k \le N$,

$$ A_N(t, \psi_k) \stackrel{\text{def}}{=} \int_{\Omega} \rho_N u_N \otimes u_N : \nabla \psi_k dy, \quad B_N(t, \psi_k) \stackrel{\text{def}}{=} \int_{\Omega} \varepsilon(\nabla u_N \nabla \rho_N) \cdot \psi_k dy, \quad C_N(t, \psi_k) = \int_{\Omega} \rho_N^{\beta} \operatorname{div}(\psi_k) dy, $$

where the functions $\psi_k$ are defined in (2.27). Proceeding similarly to [17, Section 2.4], [33, Section 7.8.2], we get, for $1 \le p < \frac{4}{3}\beta$,

$$ \rho_N u_N \otimes u_N \to \rho_n u_n \otimes u_n \text{ weakly in } L^2(0, T; L_{\beta+3}^{\frac{6\beta}{\beta+3}}(\Omega)). \quad (4.30) $$

$$ \varepsilon \nabla u_N \nabla \rho_N \to \varepsilon \nabla u_n \nabla \rho_n \text{ weakly in } L^2(0, T; L_{\beta+3}^{\frac{5\beta-3}{4\beta}}(\Omega)), \quad (4.31) $$

$$ \rho_N \to \rho_n \text{ strongly in } L^p((0, T) \times \Omega). \quad (4.32) $$

It follows from (4.30)-(4.32) the following weak convergences in $L^1(0, T)$:

$$ A_N(\cdot, \psi_k) \to A_n(\cdot, \psi_k) \stackrel{\text{def}}{=} \int_{\Omega} \rho_n u_n \otimes u_n : \nabla \psi_k dy, $$

$$ B_N(\cdot, \psi_k) \to B_n(t, \psi_k) \stackrel{\text{def}}{=} \int_{\Omega} \varepsilon (\nabla u_n \nabla \rho_n) \cdot \psi_k dy, $$

$$ C_N(\cdot, \psi_k) \to C_n(t, \psi_k) \stackrel{\text{def}}{=} \int_{\Omega} \rho_n^{\beta} \operatorname{div}(\psi_k) dy. $$

Thus, we have recovered (2.18) as a limit of equation (2.29) as $N \to \infty$, establishing the existence of a solution $(\mathcal{B}_n, \rho_n, u_n)$ to system (2.16)-(2.19).

It remains to prove energy inequality (2.25). To this purpose, we follow the same lines as in [17, Section 2.4] [33, Section 7.8.3]. The strong convergence of $\chi_N$ given in (4.28) and $\rho_N$ given in (4.32), and the convergence of $\rho_N u_N \otimes u_N$ given in (4.30) ensure that, up to a subsequence,

$$ \int_{\Omega} \left( \frac{1}{2} \rho_N |u_N|^2 + P_{\delta}(\chi_N, \rho_N) \right) dy \to \int_{\Omega} \left( \frac{1}{2} \rho_n |u_n|^2 + P_{\delta}(\chi_n, \rho_n) \right) dy \text{ in } \mathcal{D}'(0, T). $$
By the lower semicontinuity of convex functionals, the weak convergence of $u_N$ given in (4.27), the strong convergence of $\chi_N$ given in (4.28) and the weak convergence of $\Pi_N u_N$ given in (4.29), we deduce

$$ \begin{aligned} & \int_0^T \int_{\Omega} \left( 2\mu_n \left| \mathbb{D}(u_n) - \frac{1}{3} \mathrm{div}(u_n) \right|^2 + \nu_n \left| \mathrm{div}(u_n) \right|^2 \right) \mathrm{d}y \mathrm{d}\tau \\ & \leq \liminf_{N \to \infty} \int_0^T \int_{\Omega} \left( 2\mu_N \left| \mathbb{D}(u_N) - \frac{1}{3} \mathrm{div}(u_N) \right|^2 + \nu_N \left| \mathrm{div} u_N \right|^2 \right) \mathrm{d}y \mathrm{d}\tau, \end{aligned} \quad (4.33) $$
$$ \int_0^T \int_{\Omega} \chi_n |u_n - \Pi_n u_n|^2 \mathrm{d}y \mathrm{d}\tau \leq \liminf_{N \to \infty} \int_0^T \int_{\Omega} \chi_N |u_N - \Pi_N u_N|^2 \mathrm{d}y \mathrm{d}\tau. $$

Using the strong convergence of $\nabla \rho_N \to \nabla \rho_n$ given in (4.27), strong convergence of $\rho_N$ given in (4.32), and Fatou's lemma, we further obtain

$$ \int_0^T \int_{\Omega} \rho_n^{\beta-2} |\nabla \rho_n|^2 \mathrm{d}y \mathrm{d}\tau \leq \liminf_{N \to \infty} \int_0^T \int_{\Omega} \rho_N^{\beta-2} |\nabla \rho_N|^2 \mathrm{d}y \mathrm{d}\tau. \quad (4.34) $$

Collecting these results, we conclude that the energy inequality (2.25) holds.

## 4.3 High penalization limit in the momentum equation

Proposition 2.3 establishes the existence of a weak solution $(\rho_n, u_n, \mathcal{B}_n)$ to (2.16)-(2.19). In this section, we prove Proposition 2.2 by passing to the limit $n \to \infty$, thereby recovering a weak solution to the $\varepsilon$-level approximation problem (2.9)-(2.12). Given initial data $(\rho_{0,\varepsilon}, u_{0,\varepsilon})$ for the $\varepsilon$-level approximation problem, we consider initial data $(\rho_{0,n}, u_{0,n})$ for the $n$-level problem satisfying

$$ \rho_{0,n} = \rho_{0,\varepsilon}, \quad u_{0,n} = u_{0,\varepsilon}. $$

The energy estimate (2.25) yields, up to a subsequence,

$$ u_n \to u_\varepsilon \text{ weakly in } L^2(0,T; H^1(\Omega)), \quad \rho_n \to \rho_\varepsilon \text{ weakly-}^* \text{ in } L^\infty(0,T; L^\beta(\Omega)), \quad \nabla \rho_n \to \nabla \rho_\varepsilon \text{ weakly in } L^2((0,T) \times \Omega). \quad (4.35) $$

It then follows from Proposition 3.3 that, up to a subsequence and for all $1 \le p < \infty$,

$$ \chi_n \to \mathbb{1}_{\mathcal{B}_\varepsilon} \text{ weakly-}^* \text{ in } L^\infty((0,T) \times \mathbb{R}^3), \text{ strongly in } \mathcal{C}([0,T]; L_{\mathrm{loc}}^p(\mathbb{R}^3)), \quad (4.36) $$

where the pair $(\mathbb{1}_{\mathcal{B}_\varepsilon}, u_\varepsilon)$ satisfies

$$ \partial_t \mathbb{1}_{\mathcal{B}_\varepsilon} + \tilde{\Pi}_\varepsilon u_\varepsilon \cdot \nabla \mathbb{1}_{\mathcal{B}_\varepsilon} = 0 \text{ in } \mathcal{D}'((0,T) \times \mathbb{R}^3), \quad \mathcal{B}_\varepsilon(0) = \mathcal{B}_0. $$

Moreover, there exist functions $x_\varepsilon, R_\varepsilon \in H^1([0,T])$ such that

$$ \mathcal{B}_\varepsilon(t) = \mathcal{B}(x_\varepsilon(t), R_\varepsilon(t)), \quad R_\varepsilon(t) \ge R_0/2 \quad \forall t \in [0, T]. $$

Finally, using the strong convergence of $\chi_n$ given in (4.36) and the weak convergence of $u_n$ given in (4.35), it follows that

$$ \tilde{\Pi}_n u_n \to \tilde{\Pi}_\varepsilon u_\varepsilon \text{ weakly in } L^2(0, T; S). \quad (4.37) $$

The penalization term in the energy estimate (2.25) implies

$$ \sqrt{\chi_n}(u_n - \Pi_n u_n) \to 0 \text{ in } L^2((0, T) \times \Omega). \quad (4.38) $$

Using the strong convergence of $\chi_n$ given in (4.36), together with the weak convergences of $u_n$ and $\Pi_n u_n$ given in (4.35) and (4.37), respectively, we obtain in the limit

$$ \mathbb{1}_{\mathcal{B}_\varepsilon}(u_\varepsilon - \Pi_\varepsilon u_\varepsilon) = 0 \text{ in } (0, T) \times \Omega. $$

It follows that

$$ u_\varepsilon = \Pi_\varepsilon u_\varepsilon \text{ in } Q_\varepsilon, $$

where $Q_\varepsilon$ is the space-time set defined in (2.10). By spherical symmetry, the pair $(u_\varepsilon, \mathbb{1}_{\mathcal{B}_\varepsilon})$ satisfies

$$ \partial_t \mathbb{1}_{\mathcal{B}_\varepsilon} + \Pi_\varepsilon u_\varepsilon \cdot \nabla \mathbb{1}_{\mathcal{B}_\varepsilon} = 0 \text{ in } \mathcal{D}'((0,T) \times \mathbb{R}^3), \quad \mathcal{B}_\varepsilon(0) = \mathcal{B}_0. $$

Consequently, the velocity $u_\varepsilon$ is compatible with the system $\{\overline{\mathcal{B}}_\varepsilon, \eta_\varepsilon\}$, where $\eta_\varepsilon$ denotes the propagator associated with $\Pi_\varepsilon u_\varepsilon$.
Following the same arguments as in [17, Section 2.4], [33, Section 7.8.1], we obtain, for all $1 \le p < \frac{4}{3}\beta$,

$$ \rho_n \to \rho_\varepsilon \text{ in } C([0, T]; L_{\text{weak}}^\beta(\Omega)) \text{ and strongly in } L^p((0, T) \times \Omega) \quad \forall 1 \le p < \frac{4}{3}\beta, \qquad (4.39) $$

$$ \rho_n u_n \to \rho_\varepsilon u_\varepsilon \text{ weakly in } L^2(0, T; L_{\beta+6}^{6\beta}(\Omega)) \text{ and weakly-}^* \text{ in } L^\infty(0, T; L_{\beta+1}^{2\beta}(\Omega)), \qquad (4.40) $$

These convergence results allow us to pass to the limit $n \to \infty$ in the regularized continuity equation (2.9).

We now consider the limit of the momentum equation (2.11). Let $\varphi \in T(Q_\varepsilon)$. Since $\overline{Q}_\varepsilon$ is a compact subset of $(0, T) \times \Omega$, there exists $\sigma_0 > 0$ and $V_b, \omega_b, \Lambda_b \in D(0, T)$ such that, for all $t \in [0, T]$,

$$ \varphi(t, \cdot) = \varphi_b(t, \cdot) \text{ in } \mathcal{B}_\varepsilon^{\sigma_0}(t) \stackrel{\text{def}}{=} \mathcal{B}(x_\varepsilon(t), R_\varepsilon(t) + \sigma_0). $$

where

$$ \varphi_b(t, x) \stackrel{\text{def}}{=} V_b(t) + \omega_b(t) \times (x - x_\varepsilon(t)) + \frac{\Lambda_b(t)}{3}(x - x_\varepsilon(t)) \quad \forall (t, x) \in (0, T) \times \mathbb{R}^3. \qquad (4.41) $$

In particular, for $0 < \sigma \le \sigma_0$, the following identity holds:

$$ \varphi = (1 - \mathbb{1}_{\mathcal{B}_\varepsilon^\sigma})\varphi + \mathbb{1}_{\mathcal{B}_\varepsilon^\sigma}\varphi_b. \qquad (4.42) $$

Moreover, using the strong convergence of $\chi_n$ given in (4.36), it follows that there exists $n_0 \in \mathbb{N}$ such that, for all $n \ge n_0$ and $t \in [0, T]$,

$$ \mathcal{B}_n(t) \subset \mathcal{B}_\varepsilon^{\sigma_0}(t), $$

and consequently

$$ \chi_n \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}} = \chi_n, \quad (1 - \chi_n)(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}) = (1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}), \quad \chi_n(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}) = 0. \qquad (4.43) $$

Using the decomposition (4.42), the identity

$$ \mathbb{D}(\varphi_b) - \frac{1}{3} \mathrm{div}(\varphi_b) = 0, \qquad (4.44) $$

and the definition of $\mu_n \stackrel{\text{def}}{=} \mu_f(1 - \chi_n) + n\chi_n$, we write

$$ \begin{aligned} \mu_n(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)) &= \mu_n(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}})(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)) + \mu_n \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}(\mathbb{D}(\varphi_b) - \frac{1}{3}\mathrm{div}(\varphi_b)), \\ &= (\mu_f(1 - \chi_n)(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}) + n\chi_n(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}))(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)). \end{aligned} $$

Applying (4.43) yields

$$ \mu_n(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)) = \mu_f(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}})(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)). $$

The weak convergence of $u_n$ given in (4.35) and the strong convergence of $\chi_n$ given in (4.36) therefore imply

$$ \int_0^T \int_{\Omega} 2\mu_n \left( \mathbb{D}(u_n) - \frac{1}{3} \mathrm{div}(u_n) \right) : \left( \mathbb{D}(\varphi) - \frac{1}{3} \mathrm{div}(\varphi) \right) + \nu_n \mathrm{div}(u_n) : \mathrm{div}(\varphi) dy d\tau \qquad (4.45) \\ \xrightarrow{n \to \infty} \int_0^T \int_{\Omega} (1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}) 2\mu_f \left( \mathbb{D}(u_\varepsilon) - \frac{1}{3} \mathrm{div}(u_\varepsilon) \right) : \left( \mathbb{D}(\varphi) - \frac{1}{3} \mathrm{div}(\varphi) \right) dy d\tau + \int_0^T \int_{\Omega} \nu_\varepsilon \mathrm{div}(u_\varepsilon) : \mathrm{div}(\varphi) dy d\tau. $$

Finally, using again the identity (4.44) together with the definition of $\mu_\varepsilon \stackrel{\text{def}}{=} \mu_f(1 - \mathbb{1}_{\mathcal{B}_\varepsilon})$, we compute

$$ \begin{aligned} \mu_f(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}})(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)) &= \mu_f(1 - \mathbb{1}_{\mathcal{B}_\varepsilon})(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}})(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)), \\ &= \mu_\varepsilon(1 - \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}})(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)) + \mu_\varepsilon \mathbb{1}_{\mathcal{B}_\varepsilon^{\sigma_0}}(\mathbb{D}(\varphi_b) - \frac{1}{3}\mathrm{div}(\varphi_b)), \\ &= \mu_\varepsilon(\mathbb{D}(\varphi) - \frac{1}{3}\mathrm{div}(\varphi)). \end{aligned} $$

Consequently, the limit in (4.45) can be written as

$$ \int_0^T \int_{\Omega} 2\mu_\varepsilon \left( \mathbb{D}(u_\varepsilon) - \frac{1}{3} \mathrm{div}(u_\varepsilon) \right) : \left( \mathbb{D}(\varphi) - \frac{1}{3} \mathrm{div}(\varphi) \right) + \nu_\varepsilon \mathrm{div}(u_\varepsilon) : \mathrm{div}(\varphi) dy d\tau. $$
Similarly, using the identity (4.43), we obtain that the penalization term vanishes for $n \ge n_0$:

$$ 
\begin{aligned}
n \int_0^T \int_{\Omega} \chi_n(u_n - \Pi_n u_n) \cdot (\varphi - \Pi_n \varphi) &= n \int_0^T \int_{\Omega} 1_{B_{\varepsilon}^{\sigma_0}} \chi_n(u_n - \Pi_n u_n) \cdot (\varphi - \Pi_n \varphi), \\
&= n \int_0^T \int_{\Omega} \chi_n(u_n - \Pi_n u_n) \cdot (\varphi_b - \Pi_n \varphi_b) = 0.
\end{aligned}
$$

We now focus on the convergence of the nonlinear term

$$ \int_0^T \int_{\Omega} \rho_n u_n \otimes u_n : \mathcal{D}(\varphi) dy d\tau. \qquad (4.46) $$

For $0 < \sigma \le \sigma_0$, we write

$$ \int_0^T \int_{\Omega} \rho_n u_n \otimes u_n : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_{\Omega} \rho_{\varepsilon} u_{\varepsilon} \otimes u_{\varepsilon} : \mathcal{D}(\varphi) dy d\tau = \sum_{i=1}^5 M_i(\sigma, n), \qquad (4.47) $$

where

$$ M_1(\sigma, n) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} (1 - 1_{B_{\varepsilon}^{\sigma}}) \rho_n u_n \otimes u_n : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_{\Omega} (1 - 1_{B_{\varepsilon}^{\sigma}}) \rho_{\varepsilon} u_{\varepsilon} \otimes u_{\varepsilon} : \mathcal{D}(\varphi) dy d\tau $$

$$ M_2(\sigma, n) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} (1_{B_{\varepsilon}^{\sigma}} - 1_{B_{\varepsilon}}) \rho_n u_n \otimes u_n : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_{\Omega} (1_{B_{\varepsilon}^{\sigma}} - 1_{B_{\varepsilon}}) \rho_{\varepsilon} u_{\varepsilon} \otimes u_{\varepsilon} : \mathcal{D}(\varphi) dy d\tau, $$

$$ M_3(\sigma, n) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} (1_{B_{\varepsilon}} - \chi_n) \rho_n u_n \otimes u_n : \mathcal{D}(\varphi) dy d\tau, $$

$$ M_4(\sigma, n) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} \chi_n \rho_n u_n \otimes (u_n - \Pi_n u_n) : \mathcal{D}(\varphi) dy d\tau, $$

$$ M_5(\sigma, n) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} \chi_n \rho_n u_n \otimes \Pi_n u_n : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_{\Omega} 1_{B_{\varepsilon}} \rho_{\varepsilon} u_{\varepsilon} \otimes \Pi_{\varepsilon} u_{\varepsilon} : \mathcal{D}(\varphi) dy d\tau. $$

This decomposition separates the contributions arising from the fluid region, the fluid-bubble interface, the approximation of the bubble domain, and the penalization errors. Using the convergence of $u_n$ given in (4.35) together with the convergence of $\rho_n u_n$ given in (4.40), we deduce that both $\rho u_n \otimes u_n$ and $\rho u_{\varepsilon} \otimes u_{\varepsilon}$ belong to $L^2(0, T; L^{\frac{6\beta}{\beta+3}}(\Omega))$, and that the sequence $\{\rho_n u_n \otimes u_n\}$ is uniformly bounded in $L^2(0, T; L^{\frac{6\beta}{\beta+3}}(\Omega))$. Moreover, since for all $1 \le p < \infty$,

$$ 1_{B_{\varepsilon}^{\sigma}} \xrightarrow{\sigma \to 0} 1_{B_{\varepsilon}} \text{ in } C([0, T]; L^p(\Omega)), \qquad (4.48) $$

it follows that

$$ M_2(\sigma, n) \xrightarrow{\sigma \to 0} 0 \text{ uniformly in } n. $$

Similarly, by applying (4.36) in place of (4.48), we obtain

$$ M_3(\sigma, n) \xrightarrow{n \to \infty} 0 \text{ uniformly in } \sigma. $$

Finally, the convergence (4.40) implies that $\rho_n u_n$ is uniformly bounded in $L^\infty(0, T; L^{\frac{2\beta}{\beta+1}}(\Omega))$ and combined with (4.38), this yields

$$ M_4(\sigma, n) \xrightarrow{n \to \infty} 0 \text{ uniformly in } \sigma. $$

Let $\theta > 0$, there exists $(\sigma_1, n_1)$ such that for all $(\sigma, n)$ satisfying $0 \le \sigma \le \min\{\sigma_0, \sigma_1\}$ and $n \ge \max\{n_0, n_1\}$,

$$ |M_2(\sigma, n)| \le \theta/5, \quad |M_3(\sigma, n)| \le \theta/5, \quad |M_4(\sigma, n)| \le \theta/5. $$

The delicate terms are $M_1(\sigma, n)$ and $M_5(\sigma, n)$, as so far we have only established weak convergence for the sequences $\{\rho_n u_n\}$, $\{u_n\}$ and $\{\Pi_n u_n\}$. We first consider the convergence of the term $M_5(\sigma, n)$. Recall that, for all $n \in \mathbb{N}$, there exist functions $V_n, \omega_n, \Lambda_n \in L^2(0, T)$ and $V_{\varepsilon}, \omega_{\varepsilon}, \Lambda_{\varepsilon} \in L^2(0, T)$ such that

$$ \Pi_n u_n = V_n + \omega_n \times (x - x_n) + \frac{\Lambda_n}{3}(x - x_n), \quad \Pi_{\varepsilon} u_{\varepsilon} = V_{\varepsilon} + \omega_{\varepsilon} \times (x - x_{\varepsilon}) + \frac{\Lambda_{\varepsilon}}{3}(x - x_{\varepsilon}), \qquad (4.49) $$

and, by (4.37),

$$ V_n \to V_{\varepsilon}, \quad \omega_n \to \omega_{\varepsilon}, \quad \Lambda_n \to \Lambda_{\varepsilon} \text{ weakly in } L^2(0, T), \quad x_n \to x_{\varepsilon} \text{ in } C([0, T]). \qquad (4.50) $$
We introduce the functions $\tilde{V}_n, \tilde{\omega}_n, \tilde{\Lambda}_n$ and $\tilde{V}_\varepsilon, \tilde{\omega}_\varepsilon, \tilde{\Lambda}_\varepsilon$ in $L^2(0, T)$ defined, for $i \in \{1, 2, 3\}$, by

$$ \begin{aligned} \tilde{V}_n \cdot e_i &= (\rho_n u_n, \chi_n e_i)_\Omega, & \tilde{\omega}_n \cdot e_i &= (\rho_n u_n, \chi_n e_i \times (x-x_n))_\Omega, & \tilde{\Lambda}_n(t) &= (\rho_n u_n, \chi_n(x-x_n))_\Omega, \\ \tilde{V}_\varepsilon \cdot e_i &= (\rho_\varepsilon u_\varepsilon, \mathbb{1}_{\mathcal{B}_\varepsilon} e_i)_\Omega, & \tilde{\omega}_\varepsilon \cdot e_i &= (\rho_\varepsilon u_\varepsilon, \mathbb{1}_{\mathcal{B}_\varepsilon} e_i \times (x-x_\varepsilon))_\Omega, & \tilde{\Lambda}_\varepsilon(t) &= (\rho_\varepsilon u_\varepsilon, \mathbb{1}_{\mathcal{B}_\varepsilon}(x-x_\varepsilon))_\Omega. \end{aligned} $$

where $\{e_i\}$ is the canonical orthonormal basis of $\mathbb{R}^3$. Since, for $n \ge n_0$,

$$ \chi_n \operatorname{div}(\varphi) = \Lambda_b \text{ and } \mathbb{1}_{\mathcal{B}_\varepsilon} \operatorname{div}(\varphi) = \Lambda_b, $$

it follows that

$$ M_5(\sigma, n) = \int_0^T \Lambda_b(V_n \cdot \tilde{V}_n + \omega_n \cdot \tilde{\omega}_n + \Lambda_n \tilde{\Lambda}_n - V_\varepsilon \cdot \tilde{V}_\varepsilon + \omega_\varepsilon \cdot \tilde{\omega}_\varepsilon + \Lambda_\varepsilon \tilde{\Lambda}_\varepsilon) d\tau. $$

Consequently, to establish the convergence of the term $M_5(\sigma, n)$, it is sufficient to prove the strong convergence in time of $(\tilde{V}_n, \tilde{\omega}_n, \tilde{\Lambda}_n)$ to $(\tilde{V}_\varepsilon, \tilde{\omega}_\varepsilon, \tilde{\Lambda}_\varepsilon)$. Let $i \in \{1, 2, 3\}$, the idea is to use

$$ \chi_n e_i, \quad \chi_n e_i \times (x-x_n), \quad \chi_n e_i(x-x_n), $$

as test functions in the momentum equation (2.18). However, since these functions are not regular enough to be used directly as test functions, we introduce $\nu$-regularizations defined by

$$ \tilde{V}_{i,n}^{\nu} \stackrel{\text{def}}{=} (\rho_n u_n, \chi_n^{\nu} e_i)_{\Omega}, \quad \tilde{\omega}_n^{\nu} \cdot e_i = (\rho_n u_n, \chi_n^{\nu} e_i \times (x-x_n))_{\Omega}, \quad \tilde{\Lambda}_n(t) = (\rho_n u_n, \chi_n^{\nu}(x-x_n))_{\Omega}, $$

where $\chi_n^{\nu}$ is a smooth approximation of the color function $\chi_n$ defined by

$$ \chi_n^{\nu} \stackrel{\text{def}}{=} \mathcal{K}^{\nu} \star \mathbb{1}_{\mathcal{B}(x_n, (1+\nu)R_n)}, $$

with $\mathcal{K}^{\nu}$ a standard mollifier in $(0, T) \times \Omega$. In particular, the following properties hold:

$$ \chi_n^{\nu} \in \mathcal{D}((0, T) \times \Omega), \quad \chi_n(\Pi_n(\chi_n^{\nu} e_i) - \chi_n^{\nu} e_i) = 0, \quad \chi_n \xrightarrow{\nu \to 0} \chi_n \text{ in } \mathcal{C}([0, T]; L^p(\Omega)) \text{ uniformly in } n. $$

**Lemma 4.1.** For all $\nu > 0$, there exists a non-decreasing map $\phi_\nu$ such that, for all $i \in \{1, 2, 3\}$,

$$ \tilde{V}_{i, \phi_\nu(n)}^{\nu} \to \tilde{V}_i^{\nu}, \quad \tilde{\omega}_{i, \phi_\nu(n)}^{\nu} \to \tilde{\omega}_i^{\nu}, \quad \tilde{\Lambda}_{i, \phi_\nu(n)}^{\nu} \to \tilde{\Lambda}_i^{\nu} \text{ in } \mathcal{C}([0, T]) $$

**Proof of Lemma 4.1.** We only prove the result for $\tilde{V}_n^{\nu}$, the proof for $\tilde{\omega}_n^{\nu}$ and $\tilde{\Lambda}_n^{\nu}$ follows the same arguments. Let $i \in \{1, 2, 3\}$ and $t_1, t_2 \in [0, T]$, setting $\varphi_{i,n}^{\nu} \stackrel{\text{def}}{=} \chi_n^{\nu} e_i$ as a test function in (2.18), we have on one hand,

$$ \begin{aligned} & ((\rho_n u_n)(t_2) - (\rho_n u_n)(t_1), \varphi_{i,n}^{\nu}(t_2))_{\Omega} = \int_{t_1}^{t_2} \int_{\Omega} \rho_n u_n \otimes u_n : \mathcal{D}(\varphi_{i,n}^{\nu}(t_2)) dy d\tau \\ & - \int_{t_1}^{t_2} \int_{\Omega} \mu_n \left( \mathcal{D}(u_n) - \frac{1}{3} \operatorname{div}(u_n) \mathbb{I}_3 \right) : \left( \mathcal{D}(\varphi_{i,n}^{\nu}(t_2)) - \frac{1}{3} \operatorname{div}(\varphi_{i,n}^{\nu}(t_2)) \mathbb{I}_3 \right) dy d\tau + \int_{t_1}^{t_2} \int_{\Omega} \nu_n \operatorname{div}(u_n) \operatorname{div}(\varphi_{i,n}^{\nu}(t_2)) dy d\tau \\ & + \int_{t_1}^{t_2} \int_{\Omega} p(\rho_n) \operatorname{div}(\varphi_{i,n}^{\nu}(t_2)) dy d\tau - \int_{t_1}^{t_2} \int_{\Omega} \varepsilon (\nabla \rho_n \cdot \nabla) u_n \cdot \varphi_{i,n}^{\nu}(t_2) dy d\tau \\ & - \int_{t_1}^{t_2} \int_{\Omega} \rho_n g \cdot \varphi_{i,n}^{\nu}(t_2) dy d\tau + \int_{t_1}^{t_2} \int_{\Omega} \frac{\kappa_n}{R_n} \operatorname{div}(\varphi_{i,n}^{\nu}(t_2)) dy d\tau. \end{aligned} $$

Using Hölder inequalities together with uniform bounds on $\rho_n u_n$, $u_n$, $\rho_n$ and the regularity of $\varphi_{i,n}^{\nu} \stackrel{\text{def}}{=} \chi_n^{\nu} e_i$, we obtain the estimate

$$ \begin{aligned} & \int_{t_1}^{t_2} \|\rho_n u_n \otimes u_n\|_{L^{\frac{6\beta}{4\beta+3}}(\Omega)} \|\nabla \varphi_{i,n}^{\nu}(t_2)\|_{L^{\frac{6\beta}{2\beta-3}}(\Omega)} + (\mu_f + \nu_f) \|\nabla u_n\|_{L^2(\Omega)} \|\nabla \varphi_{i,n}^{\nu}(t_2)\|_{L^2(\Omega)} + \nu_b \|\nabla u_n\|_{L^2(\Omega)} \|\nabla \varphi_{i,n}^{\nu}(t_2)\|_{L^2(\Omega)} \\ & + a_f \|\rho_n\|_{L^{\beta+1}(\Omega)}^{\gamma_f} \|\nabla \varphi_{i,n}^{\nu}(t_2)\|_{L^{\frac{\beta+1}{\beta+1-\gamma_f}}(\Omega)} + a_b \|\rho_n\|_{L^{\beta+1}(\Omega)}^{\gamma_b} \|\nabla \varphi_{i,n}^{\nu}(t_2)\|_{L^{\frac{\beta+1}{\beta+1-\gamma_b}}(\Omega)} + \delta \|\rho_n\|_{L^{\beta+1}(\Omega)}^{\beta} \|\nabla \varphi_{i,n}^{\nu}(t_2)\|_{L^{\beta+1}(\Omega)} \\ & + \varepsilon \|\nabla \rho_n \cdot \nabla u_n\|_{L^{\frac{5\beta-3}{4\beta}}(\Omega)} \|\varphi_{i,n}^{\nu}(t_2)\|_{L^{\frac{5\beta-3}{\beta-3}}(\Omega)} + \|\rho_n g\|_{L^{\beta+1}(\Omega)} \|\varphi_{i,n}^{\nu}(t_2)\|_{L^{\frac{\beta+1}{\beta}}(\Omega)} + \frac{2\kappa_b}{R_0} \|\nabla \varphi_{i,n}^{\nu}(t_2)\|_{L^1(\Omega)} d\tau, \end{aligned} \quad (4.51) $$
where for $1 \le p < \infty$, by Young's inequality,

$$ \|\varphi_{n,i}^{\nu}\|_{W^{1,\infty}(0,T;L^{p}(\Omega))} \le \|1_{\mathcal{B}(x_n,(1+\nu)R_n)}\|_{L^{\infty}((0,T)\times\Omega)} \|\mathcal{K}^{\nu}\|_{W^{1,1}(\Omega)} \le \|\mathcal{K}^{\nu}\|_{W^{1,1}(\Omega)}. $$

On the other hand,

$$ \begin{aligned} (\rho_n u_n(t_2), \varphi_n^{\nu}(t_2) - \varphi_n^{\nu}(t_1))_{\Omega} &= \left( \rho_n u_n(t_2) \cdot e_i, \int_{t_1}^{t_2} \mathcal{K}^{\nu} * \partial_t 1_{\mathcal{B}(x_n,(1+\nu)R_n)} d\tau \right)_{\Omega} \\ &= - \int_{t_1}^{t_2} \left( \rho_n u_n(t_2) \cdot e_i, \mathcal{K}^{\nu} * (\Pi_n u_n \cdot \nabla 1_{\mathcal{B}(x_n,(1+\nu)R_n)}) \right)_{\Omega} d\tau \\ &= - \int_{t_1}^{t_2} \left( \rho_n u_n(t_2) \cdot e_i, \nabla \mathcal{K}^{\nu} * (\Pi_n u_n 1_{\mathcal{B}(x_n,(1+\nu)R_n)}) \right)_{\Omega} d\tau + \int_{t_1}^{t_2} \left( \rho_n u_n(t_2) \cdot e_i, \mathcal{K}^{\nu} * (\mathrm{div}(\Pi_n u_n) 1_{\mathcal{B}(x_n,(1+\nu)R_n)}) \right)_{\Omega} d\tau. \end{aligned} $$

Consequently, using Holder's and Young's inequality, we obtain

$$ |(\rho_n u_n(t_2), \varphi_n^{\nu}(t_2) - \varphi_n^{\nu}(t_1))_{\Omega}| \le \|\rho_n u_n\|_{L^{\infty}(0,T;L^{\frac{2\beta}{\beta+1}}(\Omega))} \|1_{\mathcal{B}(x_n,(1+\nu)R_n)}\|_{L^{\infty}((0,T)\times\Omega)} \|\mathcal{K}^{\nu}\|_{W^{1,1}(\Omega)} \int_{t_1}^{t_2} \|\Pi_n u_n\|_{W^{1,\frac{2\beta}{2\beta-1}}(\Omega)} d\tau. $$

We deduce that $\tilde{V}_{i,\nu}^n$ is uniformly equicontinuous on $(0,T)$ uniformly in $n$ (although not in $\nu$). The Arzela-Ascoli theorem implies that there exists a non-decreasing function $\phi_\nu: \mathbb{N} \to \mathbb{N}$ and $\tilde{V}_i^\nu \in C([0,T])$ such that

$$ \tilde{V}_{i,\phi_\nu(n)}^{\nu} \xrightarrow{n \to \infty} \tilde{V}_i^{\nu} \text{ in } C([0,T]). $$

This concludes the proof of the lemma.

<figure><img src="image_4.png" /></figure>

The next lemma addresses the convergence of $\{\tilde{V}_n\}$, $\{\tilde{\omega}_n\}$ and $\{\tilde{\Lambda}_n\}$.

**Lemma 4.2.** Up to a subsequence, for all $i \in \{1, 2, 3\}$,

$$ \tilde{V}_{i,n} \to \tilde{V}_{i,\varepsilon}, \quad \tilde{\omega}_{i,n} \to \tilde{\omega}_{i,\varepsilon}, \quad \tilde{\Lambda}_n \to \tilde{\Lambda}_\varepsilon \text{ in } C([0,T]). $$

*Proof of Lemma 4.2.* Similarly to Lemma 4.1, we focus on the proof for $\tilde{V}_n$, as the proof is analogous for $\tilde{\omega}_n$ and $\tilde{\Lambda}_n$. Let $i \in \{1, 2, 3\}$ and consider the regularized sequence $\{\tilde{V}_{i,n}^\nu\}$. By construction, we have

$$ \tilde{V}_{i,n}^{\nu} \xrightarrow{\sigma \to 0} \tilde{V}_{i,n} \text{ in } C([0,T]) \text{ uniformly in } n. $$

In particular, for any $\theta > 0$, there exists $\nu > 0$ sufficiently small such that for all $n \in \mathbb{N}$,

$$ \|\tilde{V}_{i,n}^{\nu} - \tilde{V}_{i,n}\|_{L^{\infty}(0,T)} \le \theta/3. $$

By Lemma 4.1, there exists $n_0(\nu) \in \mathbb{N}$ such that for $m, l \ge n_0(\nu)$,

$$ \|V_{i,\phi_\nu(m)}^{\nu} - V_{i,\phi_\nu(l)}^{\nu}\|_{L^{\infty}(0,T)} \le \theta/3. $$

Consequently, for $m, l \ge n_0(\nu)$,

$$ \|\tilde{V}_{i,\phi_\nu(m)} - \tilde{V}_{i,\phi_\nu(l)}\|_{L^{\infty}(0,T)} \le \|\tilde{V}_{i,\phi_\nu(m)} - \tilde{V}_{i,\phi_\nu(l)}^{\nu}\|_{L^{\infty}(0,T)} + \|\tilde{V}_{i,\phi_\nu(m)}^{\nu} - \tilde{V}_{i,\phi_\nu(l)}^{\nu}\|_{L^{\infty}(0,T)} + \|\tilde{V}_{i,\phi_\nu(l)}^{\nu} - \tilde{V}_{i,\phi_\nu(l)}\|_{L^{\infty}(0,T)} \le \theta. \quad (4.52) $$

Applying (4.52) with $\theta = 1, \frac{1}{2}, \frac{1}{4}, \dots$ and using a standard diagonal argument, we can extract a subsequence $\tilde{V}_{i,\psi(n)}$ that is Cauchy in $C([0,T])$ and hence converges in $C([0,T])$. To identify the limit, we consider $\varphi \in D(0,T)$. Using the weak convergence of $\rho_n u_n$ given in (4.40), the strong convergence of $\chi_n$ given in (4.36), we have

$$ \int_0^T \varphi(\rho_n u_n, \chi_n e_i)_\Omega d\tau \to \int_0^T \varphi(\rho_\varepsilon u_\varepsilon, \chi_\varepsilon e_i)_\Omega d\tau. $$

This identifies the limit as $\tilde{V}_{i,\varepsilon}$ and concludes the proof of the lemma.

<figure><img src="image_5.png" /></figure>
From the convergences of the sequences {$V_n$}, {$\omega_n$} and {$\Lambda_n$} given in (4.50) and Lemma 4.2, we deduce

$$V_n \cdot \tilde{V}_n \to V_\epsilon \cdot \tilde{V}_\epsilon, \quad \omega_n \cdot \tilde{\omega}_n \to \omega_\epsilon \cdot \tilde{\omega}_\epsilon, \quad \Lambda_n \tilde{\Lambda}_n \to \Lambda_\epsilon \tilde{\Lambda}_\epsilon \text{ in } L^2(0,T),$$

and there exists $n_2$ such that for all $n$ satisfying $n \ge \{n_0, n_2\}$,

$$|M_5(\sigma, n)| \le \frac{\theta}{5}.$$

For the last term $M_1(\sigma, n)$, we use a result presented in [14, Section 7] in the context of fluid-solid interaction. Let $\sigma_m = \min\{\sigma_0, \sigma_1\}$, we have

$$\rho_n u_n \otimes u_n \to \rho_\epsilon u_\epsilon \otimes u_\epsilon \text{ in } L^{\frac{6\beta}{4\beta+3}}(Q_{f,\epsilon}^{\sigma_m})$$

where

$$Q_{f,\epsilon}^{\sigma_m} \stackrel{\text{def}}{=} (0,T)\backslash \overline{Q_{\epsilon}^{\sigma_m}}, \quad \overline{Q_{\epsilon}^{\sigma_m}} \stackrel{\text{def}}{=} \{(t,x) \in (0,T) \times \Omega \mid x \in (\overline{B_{\epsilon}^{\sigma_m}(t)})\}.$$

Consequently, there exists $n_3 \in \mathbb{N}$, possibly depending on $\sigma_m$, such that for all $n \ge n_3$,

$$|M_5(\sigma_m, n)| \le \theta/5.$$

Finally, we conclude that for $n \ge \max\{n_0, n_1, n_2, n_3\}$,

$$\left| \int_0^T \int_{\Omega} \rho_n u_n \otimes u_n : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_{\Omega} \rho_\epsilon u_\epsilon \otimes u_\epsilon : \mathcal{D}(\varphi) dy d\tau \right| \le \theta,$$

which leads to

$$\int_0^T \int_{\Omega} \rho_n u_n \otimes u_n : \mathcal{D}(\varphi) dy d\tau \to \int_0^T \int_{\Omega} \rho_\epsilon u_\epsilon \otimes u_\epsilon : \mathcal{D}(\varphi) dy d\tau.$$

The convergence results for the pressure terms in the momentum equation and the energy estimate are derived using arguments similar to those presented in Section 4.2.

## 4.4 Vanishing dissipation in the continuity equation and the limiting system

Proposition 2.3 establishes the existence of a weak solution $(\rho_\epsilon, u_\epsilon, \mathcal{B}_\epsilon)$ to (2.10)-(2.12). In this section, we prove Proposition 2.1 by passing to the limit $\epsilon \to 0$ in the $\epsilon$-level approximation problem. A key step of the analysis consists in identifying the pressure corresponding to the limiting system. Given initial data $(\rho_{0,\delta}, u_{0,\delta})$ for the $\delta$-level approximation problem, we construct initial data $(\rho_{0,\epsilon}, u_{0,\epsilon})$ for the $\epsilon$-level approximation problem such that

$$\rho_{0,\epsilon} > 0, \quad \rho_{0,\epsilon} \in W^{1,\infty}(\Omega), \quad \rho_{0,\epsilon} \to \rho_{0,\delta} \text{ in } L^\beta(\Omega), \quad q_{0,\epsilon} \to q_{0,\delta} \text{ in } L^{\frac{2\beta}{\beta+1}}(\Omega),$$

and

$$E_\delta(\rho_{0,\epsilon}, q_{0,\epsilon}, \mathbb{1}_{B_0}) \to E_\delta(\rho_{0,\delta}, q_{0,\delta}, \mathbb{1}_{B_0}),$$

where the definition of $E_\delta$ is given in (2.7).

The energy estimate (2.14) yields, up to a subsequence,

$$u_\epsilon \to u_\delta \text{ weakly in } L^2(0, T; H^1(\Omega)). \qquad (4.53)$$

It then follows from Proposition 3.3 that up to a subsequence and for all $1 \le p < \infty$,

$$\mathbb{1}_{B_\epsilon} \to \mathbb{1}_{B_\delta} \text{ weakly-}^* \text{ in } L^\infty((0,T) \times \Omega), \text{ strongly in } \mathcal{C}([0,T]; L_{\text{loc}}^p(\mathbb{R}^3)), \qquad (4.54)$$

where the pair $(\mathbb{1}_{B_\delta}, u_\delta)$ satisfies

$$\partial_t \mathbb{1}_{B_\delta} + u_\delta \cdot \nabla \mathbb{1}_{B_\delta} = 0 \text{ in } \mathcal{D}'((0,T) \times \mathbb{R}^3), \quad \mathcal{B}_\delta(0) = \mathcal{B}_0.$$

Moreover, the sequences {$u_\epsilon$} and {$\eta_\epsilon$} satisfy the assumptions of Corollary 3.1. It then follows that

$$\eta_\epsilon \to \eta_\delta \text{ weakly in } H^1([0, T]; \mathcal{M}), \text{ strongly in } \mathcal{C}([0, T], \mathcal{M}), \qquad (4.55)$$

and, if we define

$$\mathcal{B}_\delta(t) \stackrel{\text{def}}{=} \eta_\delta[t](\mathcal{B}_0) \quad \forall t \in [0, T],$$
then the velocity field $u_\delta$ is compatible with the system $\{\overline{\mathcal{B}}_\delta, \eta_\delta\}$.

Writing $\mathcal{B}_\delta = \mathcal{B}(x_\delta, R_\delta)$ with $x_\delta, R_\delta \in H^1([0, T])$, we have

$$R_\delta(t) \geq R_0/2 \quad \forall t \in [0, T].$$

From the convergence of $\eta_\epsilon$ given in (4.55), it follows that, for $\epsilon > 0$ sufficiently small,

$$\mathcal{B}_\delta(t) \subset \mathcal{B}_\epsilon^\sigma(t) \quad \forall t \in [0, T].$$

In particular, since $\mathcal{B}_\epsilon$ satisfies $\mathrm{dist}(\mathcal{B}_\epsilon, \partial\Omega) \geq 2\sigma$ in $[0, T]$, we deduce that

$$\mathrm{dist}(\mathcal{B}_\delta(t), \partial\Omega) \geq 3\sigma/2 \quad \forall t \in [0, T].$$

Following the same arguments as in [17, Section 3.3], [33, Section 7.9.1], we get

$$\rho_{\epsilon} \rightarrow \rho_{\delta} \text{ in } \mathcal{C}([0, T]; L_{\text{weak}}^{\beta}(\Omega)), \qquad (4.56)$$

$$\rho_{\epsilon} u_{\epsilon} \rightarrow \rho_{\delta} u_{\delta} \text{ weakly-}^* \text{ in } L^{\infty}(0, T; L^{\frac{2\beta}{\beta+1}}(\Omega)), \qquad (4.57)$$

$$\epsilon \nabla \rho_{\delta} \rightarrow 0 \text{ strongly in } L^2((0, T) \times \Omega). \qquad (4.58)$$

These convergences allow us to pass to the limit in the continuity equation (2.9). In addition, the couple $(\rho_\delta, u_\delta)$ satisfies the renormalized continuity equation (1.28) for all $b$ satisfying (1.29). Finally, by Proposition 3.4, the density $\rho_\delta$ is compatible with the system $\{\overline{\mathcal{B}}_\epsilon, \eta_\epsilon\}$.

We now address the convergence of the momentum equation (2.11). Let $\varphi \in \mathcal{T}(Q_\delta)$. Since $\overline{Q}_\delta$ is a compact subset of $(0, T) \times \Omega$, an argument analogous to that of Section 4.3 yields the existence of $\sigma_0 > 0$ and functions $V_b, \omega_b, \Lambda_b \in \mathcal{D}(0, T)$ such that, for all $t \in [0, T]$,

$$\varphi(t, \cdot) = \varphi_b(t, \cdot) \text{ in } \mathcal{B}_\delta^\sigma(t) \stackrel{\text{def}}{=} \mathcal{B}(x_\delta(t), R_\delta(t) + \sigma),$$

with $\varphi_b$ defined by (4.41). In particular, for $0 < \sigma \leq \sigma_0$, the following identity holds:

$$\varphi = (1 - \mathbb{1}_{\mathcal{B}_\delta^\sigma})\varphi + \mathbb{1}_{\mathcal{B}_\delta^\sigma}\varphi_b. \qquad (4.59)$$

Using once more the strong convergence of $\eta_\epsilon$ given in (4.55), we deduce that there exists $\epsilon_0 > 0$ such that for $0 \leq \epsilon \leq \epsilon_0$,

$$\mathcal{B}_\epsilon(t) \subset \mathcal{B}_\delta^{\sigma_0}(t) \quad \forall t \in [0, T],$$

and consequently

$$\mathbb{1}_{\mathcal{B}_\epsilon} \mathbb{1}_{\mathcal{B}_\delta^{\sigma_0}} = \mathbb{1}_{\mathcal{B}_\epsilon}, \quad (1 - \mathbb{1}_{\mathcal{B}_\epsilon})(1 - \mathbb{1}_{\mathcal{B}_\delta^{\sigma_0}}) = (1 - \mathbb{1}_{\mathcal{B}_\delta^{\sigma_0}}), \quad \mathbb{1}_{\mathcal{B}_\epsilon}(1 - \mathbb{1}_{\mathcal{B}_\delta^{\sigma_0}}) = 0.$$

The delicate term is still

$$\int_{0}^{T} \int_{\Omega} \rho_{\epsilon} u_{\epsilon} \otimes u_{\epsilon} : \mathcal{D}(\varphi) dy d\tau. \qquad (4.60)$$

When details are omitted, the arguments for the convergence of (4.60) are analogous to those presented in Section 4.3. For $0 < \sigma \leq \sigma_0$, we write

$$\int_{0}^{T} \int_{\Omega} \rho_{\epsilon} u_{\epsilon} \otimes u_{\epsilon} : \mathcal{D}(\varphi) dy d\tau - \int_{0}^{T} \int_{\Omega} \rho_{\delta} u_{\delta} \otimes u_{\delta} : \mathcal{D}(\varphi) dy d\tau = \sum_{i=1}^{4} M_{i}(\sigma, \epsilon)$$

where

$$M_1(\sigma, \epsilon) \stackrel{\text{def}}{=} \int_0^T \int_\Omega (1 - \mathbb{1}_{\mathcal{B}_\delta^\sigma}) \rho_\epsilon u_\epsilon \otimes u_\epsilon : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_\Omega (1 - \mathbb{1}_{\mathcal{B}_\delta^\sigma}) \rho_\delta u_\delta \otimes u_\delta : \mathcal{D}(\varphi) dy d\tau,$$

$$M_2(\sigma, \epsilon) \stackrel{\text{def}}{=} \int_0^T \int_\Omega (\mathbb{1}_{\mathcal{B}_\delta^\sigma} - \mathbb{1}_{\mathcal{B}_\delta}) \rho_\epsilon u_\epsilon \otimes u_\epsilon : \mathcal{D}(\varphi) dy d\tau + \int_0^T \int_\Omega (\mathbb{1}_{\mathcal{B}_\delta^\sigma} - \mathbb{1}_{\mathcal{B}_\delta}) \rho_\delta u_\delta \otimes u_\delta : \mathcal{D}(\varphi) dy d\tau,$$

$$M_3(\sigma, \epsilon) \stackrel{\text{def}}{=} \int_0^T \int_\Omega (\mathbb{1}_{\mathcal{B}_\delta} - \mathbb{1}_{\mathcal{B}_\epsilon}) \rho_\epsilon u_\epsilon \otimes u_\epsilon : \mathcal{D}(\varphi) dy d\tau,$$

$$M_4(\sigma, \epsilon) \stackrel{\text{def}}{=} \int_0^T \int_\Omega \mathbb{1}_{\mathcal{B}_\epsilon} \rho_\epsilon u_\epsilon \otimes u_\epsilon : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_\Omega \mathbb{1}_{\mathcal{B}_\delta} \rho_\delta u_\delta \otimes u_\delta : \mathcal{D}(\varphi) dy d\tau.$$
Let $\theta > 0$, there exists $\sigma_1, \varepsilon_1 > 0$ such that for $0 < \sigma \le \min\{\sigma_0, \sigma_1\}$ and $0 < \varepsilon \le \min\{\varepsilon_0, \varepsilon_1\}$,

$$|M_2(\sigma, \varepsilon)| \le \theta/4, \quad |M_3(\sigma, \varepsilon)| \le \theta/4.$$

Arguing as in Lemma 4.2, we get for $i \in \{1, 2, 3\}$:

$$(\rho_{\varepsilon} u_{\varepsilon}, \mathbb{1}_{\mathcal{B}_{\varepsilon}} e_i)_{\Omega} \xrightarrow{\varepsilon \to 0} (\rho_{\delta} u_{\delta}, \mathbb{1}_{\mathcal{B}_{\delta}} e_i)_{\Omega} \text{ in } \mathcal{C}([0, T]), \quad (4.61)$$

$$(\rho_{\varepsilon} u_{\varepsilon}, \mathbb{1}_{\mathcal{B}_{\varepsilon}} e_i \times (x - x_{\varepsilon}))_{\Omega} \xrightarrow{\varepsilon \to 0} (\rho_{\delta} u_{\delta}, \mathbb{1}_{\mathcal{B}_{\delta}} e_i \times (x - x_{\delta}))_{\Omega} \text{ in } \mathcal{C}([0, T]), \quad (4.62)$$

$$(\rho_{\varepsilon} u_{\varepsilon}, \mathbb{1}_{\mathcal{B}_{\varepsilon}}(x - x_{\varepsilon}))_{\Omega} \xrightarrow{\varepsilon \to 0} (\rho_{\delta} u_{\delta}, \mathbb{1}_{\mathcal{B}_{\delta}}(x - x_{\delta}))_{\Omega} \text{ in } \mathcal{C}([0, T]). \quad (4.63)$$

Using the convergence (4.61)-(4.63), we obtain

$$\lim_{\varepsilon \to 0} \int_0^T \int_{\Omega} \mathrm{div}(\varphi_b) \mathbb{1}_{\mathcal{B}_{\varepsilon}} \rho_{\varepsilon} |u_{\varepsilon}|^2 \mathrm{d}y \mathrm{d}\tau = \int_0^T \int_{\Omega} \mathrm{div}(\varphi_b) \mathbb{1}_{\mathcal{B}_{\delta}} \rho_{\delta} |u_{\delta}|^2 \mathrm{d}y \mathrm{d}\tau.$$

Consequently, there exists $\varepsilon_2 > 0$ such that for all $0 < \varepsilon \le \varepsilon_2$,

$$|M_4(\sigma, \varepsilon)| \le \theta/4.$$

Finally, setting $\sigma_m = \min\{\sigma_0, \sigma_1\}$ and following similar arguments to the ones presented in [14, Section 8] and already presented in Section 4.3, we can prove that

$$\rho_{\varepsilon} u_{\varepsilon} \otimes u_{\varepsilon} \to \rho_{\delta} u_{\delta} \otimes u_{\delta} \text{ in } L^2(Q_{f,\delta}^{\sigma_m}) \quad (4.64)$$

where

$$Q_{f,\delta}^{\sigma_m} \stackrel{\text{def}}{=} (0,T) \setminus \overline{Q_{\delta}^{\sigma_m}}, \quad \overline{Q_{\delta}^{\sigma_m}} \stackrel{\text{def}}{=} \{(t,x) \in (0,T) \times \Omega \mid x \in (\overline{\mathcal{B}_{\delta}^{\sigma_m}(t)})\}.$$

Consequently, there exists $\varepsilon_3 > 0$, possibly depending on $\sigma_m$, such that for $0 < \varepsilon \le \varepsilon_3$,

$$|M_5(\sigma_m, \varepsilon)| \le \theta/4.$$

Finally, for $0 < \varepsilon \le \min\{\varepsilon_1, \varepsilon_2, \varepsilon_3\}$,

$$\left| \int_0^T \int_{\Omega} \rho_{\varepsilon} u_{\varepsilon} \otimes u_{\varepsilon} : \mathcal{D}(\varphi) \mathrm{d}y \mathrm{d}\tau - \int_0^T \int_{\Omega} \rho_{\delta} u_{\delta} \otimes u_{\delta} : \mathcal{D}(\varphi) \mathrm{d}y \mathrm{d}\tau \right| \le \theta, \quad (4.65)$$

which leads to

$$\int_0^T \int_{\Omega} \rho_{\varepsilon} u_{\varepsilon} \otimes u_{\varepsilon} : \mathcal{D}(\varphi) \mathrm{d}y \mathrm{d}\tau \to \int_0^T \int_{\Omega} \rho_{\delta} u_{\delta} \otimes u_{\delta} : \mathcal{D}(\varphi) \mathrm{d}y \mathrm{d}\tau. \quad (4.66)$$

We now aim to identify the pressure term. First, we first establish that $\{\rho_{\varepsilon}\}$ is locally uniformly bounded in $L^{\beta+1}(Q_{f,\delta})$, where

$$Q_{f,\delta} = (0,T) \times \Omega \setminus \overline{Q_{\delta}}$$

with $\overline{Q}_{\delta}$ defined by (2.3).

**Lemma 4.3.** For any compact $K_f \subset Q_{f,\delta}$, there exists a constant $c$ independent of $\varepsilon$ such that

$$\|\rho_{\varepsilon}\|_{L^{\beta+1}(K_f)} \le c(K_f, E_{\delta}(\rho_{0,\varepsilon}, q_{0,\varepsilon}, \mathbb{1}_{\mathcal{B}_0}), \|g\|_{L^{\infty}((0,T) \times \mathbb{R}^3)}, \kappa_b). \quad (4.67)$$

*Proof.* The proof follows the same strategy as in [14, Lemma 8.1], [17, Lemma 3.1]. It relies on testing the momentum equation (2.11) with functions of the form

$$\varphi = \psi(t) \mathcal{B} \left( \rho_{\varepsilon} - |K_f|^{-1} \int_{K_f} \rho_{\varepsilon} \mathrm{d}y \right),$$

where $\psi \in \mathcal{D}(0, T)$ and $\mathcal{B}$ denotes the Bogovskii operator on $K_f$.

□
Applying Lemma 4.3, we deduce that

$$ \rho_{\varepsilon}^{\beta} \to \overline{\rho_{\delta}^{\beta}} \text{ weakly in } L^{\frac{\beta+1}{\beta}}(K_f), \quad (4.68) $$

$$ \rho_{\varepsilon}^{\gamma_f} \to \overline{\rho_{\delta}^{\gamma_f}} \text{ weakly in } L^{\frac{\gamma_f+1}{\gamma_f}}(K_f), \quad (4.69) $$

for any compact $K_f \subset Q_{f,\delta}$. Together with (4.53), (4.56), and (4.64), the convergences (4.68)-(4.69) allow us to pass to the limit $\varepsilon \to 0$ in (2.11). Consequently, for any test function $\varphi \in D(Q_{f,\delta})$, we obtain

$$ \int_0^T \int_{\mathbb{R}^3} (\rho_{\delta} u_{\delta}) \partial_t \varphi + (\rho_{\delta} u_{\delta} \otimes u_{\delta}) : \mathbb{D}(\varphi) + (a_f \overline{\rho_{\delta}^{\gamma_f}} + \delta \overline{\rho_{\delta}^{\beta}}) \mathrm{div}(\varphi) dy d\tau = \int_0^T \int_{\mathbb{R}^3} \mathbb{T}_f(u_{\delta}) : \mathbb{D}(\varphi) - \rho g \cdot \varphi dy d\tau. $$

The next step consists in establishing the strong convergence of the density. Specifically, we aim to identify the nonlinear pressure term $(a_f \overline{\rho_{\delta}^{\gamma_f}} + \delta \overline{\rho_{\delta}^{\beta}})$ by proving that $\overline{\rho_{\delta}^{\gamma_f}} = \rho_{\delta}^{\gamma_f}$ and $\overline{\rho_{\delta}^{\beta}} = \rho_{\delta}^{\beta}$. To this end, we introduce the effective viscous flux

$$ a_f \rho_{\varepsilon}^{\gamma_f} + \delta \rho_{\varepsilon}^{\beta} - \left( \frac{4}{3} \mu_f + \nu_f \right) \mathrm{div}(u_{\varepsilon}). $$

which exhibits improved compactness properties compared to the pressure or velocity fields separately.

**Lemma 4.4.** For any test function $\varphi \in D(Q_{f,\delta})$,

$$ \lim_{\varepsilon \to 0} \int_0^T \int_{\mathbb{R}^3} \varphi \left( a_f \overline{\rho_{\varepsilon}^{\gamma_f}} + \delta \overline{\rho_{\varepsilon}^{\beta}} - \left( \nu_f + \frac{4}{3} \mu_f \right) \mathrm{div}(u_{\varepsilon}) \right) \rho_{\varepsilon} dy d\tau = \int_0^T \int_{\mathbb{R}^3} \varphi \left( a_f \overline{\rho_{\delta}^{\gamma_f}} + \delta \overline{\rho_{\delta}^{\beta}} - \left( \nu_f + \frac{4}{3} \mu_f \right) \mathrm{div}(u_{\delta}) \right) \rho_{\delta} dy d\tau $$

*Proof.* The proof of this Lemma can be found in [30, Theorem 5.1, Appendix B], [17, Lemma 3.2], [14, Lemma 8.2], [33, Lemma 7.50]. □

Applying Lemma 4.4 and the monotonicity of the mappings $\rho \to \rho^{\gamma_f}$, $\rho \to \rho^{\beta}$, we get

$$ \overline{\rho_{\delta} \mathrm{div}(u_{\delta})} \geq \rho_{\delta} \mathrm{div}(u_{\delta}) \text{ on } Q_{f,\delta}, \quad (4.70) $$

where $\overline{\rho_{\delta} \mathrm{div}(u_{\delta})}$ is defined by

$$ \rho_{\varepsilon} \mathrm{div}(u_{\varepsilon}) \to \overline{\rho_{\delta} \mathrm{div}(u_{\delta})} \text{ weakly in } L^1((0, T) \times \Omega). $$

Moreover, on any compact $K_b \subset Q_{b,\delta}$, there exists $\varepsilon_0$, possibly depending on $K_b$, such that for $0 < \varepsilon \le \varepsilon_0$,

$$ \mathrm{div}(u_{\delta}) = \Lambda_{\delta}, \mathrm{div}(u_{\varepsilon}) = \Lambda_{\varepsilon} \text{ and } \Lambda_{\delta} \to \Lambda_{\varepsilon} \text{ weakly in } L^2(0, T). \quad (4.71) $$

From the convergence of $\rho_{\varepsilon}$ given in (4.56) and of $\mathrm{div}(u_{\varepsilon})$ given in (4.71), we deduce

$$ \overline{\rho_{\delta} \mathrm{div}(u_{\delta})} = \rho_{\delta} \mathrm{div}(u_{\delta}) \text{ in } Q_{\delta}. $$

Finally, we obtain

$$ \overline{\rho_{\delta} \mathrm{div}(u_{\delta})} \geq \rho_{\delta} \mathrm{div}(u_{\delta}) \text{ on } (0, T) \times \mathbb{R}^3. \quad (4.72) $$

The remaining arguments follow standard results from the theory of compressible fluids and can be found in [30], [17][Section 3.5], [33][Lemma 7.51], [14][Section 8]. For completeness, we outline the main steps. Noting that $b(s) = s \log(s)$ satisfies the condition (1.29) and $\rho_{\delta}$ satisfies the renormalized continuity equation (1.28), we deduce, for all $t \in [0, T]$,

$$ \int_0^t \int_{\Omega} \rho_{\delta} \mathrm{div}(u_{\delta}) dy d\tau = \int_{\Omega} \rho_{0,\delta} \log(\rho_{0,\delta}) dy - \int_{\Omega} \rho_{\delta}(t) \log(\rho_{\delta}(t)) dy. \quad (4.73) $$

Since $\rho_{\varepsilon}$ satisfies the regularized continuity equation (2.9) almost everywhere on $(0, T) \times \Omega$, it follows that, for any convex $b$, globally Lipschitz on $[0, \infty)$,

$$ \partial_t b(\rho_{\varepsilon}) + \mathrm{div}(b(\rho_{\varepsilon})u_{\varepsilon}) + (b'(\rho_{\varepsilon})\rho_{\varepsilon} - b(\rho_{\varepsilon})) \mathrm{div}(u_{\varepsilon}) - \varepsilon \Delta b(\rho_{\varepsilon}) \le 0. $$

In particular taking again $b(s) = s \log(s)$ and integrating on $(0, t) \times \Omega$, we obtain

$$ \int_0^t \int_{\Omega} \rho_{\varepsilon} \mathrm{div}(u_{\varepsilon}) dy d\tau \le \int_{\Omega} \rho_{0,\varepsilon} \log(\rho_{0,\varepsilon}) dy - \int_{\Omega} \rho_{\varepsilon}(t) \log(\rho_{\varepsilon}(t)) dy. \quad (4.74) $$
The relations (4.72), (4.73) together with (4.74) lead to

$$ \limsup_{\varepsilon \to 0} \int_{\Omega} \rho_{\delta}(t) \log(\rho_{\delta}(t)) dy \le \int_{\Omega} \rho_{\delta}(t) \log(\rho_{\delta}(t)) dy \quad \forall t \in [0, T] $$
$$ \rho_{\varepsilon} \to \rho_{\delta} \text{ in } L^{p}((0, T) \times \Omega). \qquad (4.75) $$

Consequently, $\overline{\rho_{\delta}^{\gamma f}} + \delta \overline{\rho_{\delta}^{\beta}} = \overline{\rho_{\delta}^{\gamma f}} + \delta \overline{\rho_{\delta}^{\beta}}$ on $(0, T) \times \Omega$.

Finally, using the following convergences

$$ \begin{aligned} \rho_{\varepsilon} u_{\varepsilon} \otimes u_{\varepsilon} &\to \rho_{\delta} u_{\delta} \otimes u_{\delta} \quad \text{weakly in } L^{2}(0, T; L^{\frac{6\beta}{4\beta+3}}(\Omega)), \\ \rho_{\varepsilon} &\to \rho_{\delta} \quad \text{strongly in } L^{p}((0, T) \times \Omega) \quad \forall 1 \le p < \beta + 1, \end{aligned} $$

we have

$$ \begin{aligned} \int_{\Omega} \rho_{\varepsilon} |u_{\varepsilon}|^2 dy &\to \int_{\Omega} \rho_{\delta} |u_{\delta}|^2 dy \quad \text{weakly in } L^2(0, T), \\ \int_{\Omega} P_{\delta}(\rho_{\varepsilon}, \mathbb{1}_{B_{\varepsilon}}) dy &\to \int_{\Omega} P_{\delta}(\rho_{\delta}, \mathbb{1}_{B_{\delta}}) dy \quad \text{weakly in } L^{\frac{\beta+1}{\beta}}(0, T). \end{aligned} $$

and

$$ \int_{\Omega} (\rho_{\varepsilon} |u_{\varepsilon}|^2 + P_{\delta}(\rho_{\varepsilon}, \mathbb{1}_{B_{\varepsilon}})) dy \xrightarrow{\varepsilon \to 0} \int_{\Omega} (\rho_{\delta} |u_{\delta}|^2 + P(\rho_{\delta}, \mathbb{1}_{B_{\delta}})) dy. $$

Due to the weak lower semicontinuity of the $L^2$ norms, the convergence of $u_{\varepsilon}$ given in (4.53), the strong convergence of $\rho_{\varepsilon}$ given in (4.75), the strong convergence of $\chi_{\varepsilon}$ given in (4.54), we can pass to the limit as $\varepsilon \to 0$ in the other terms of the energy inequality (2.14) to establish the estimate (2.8).

# 5 Proof of the main result

Proposition 2.1 establishes the existence of a weak solution $(B_{\delta}, \rho_{\delta}, u_{\delta})$ to system (2.1)-(2.5). In this section, we analyse the convergence of these solutions as $\delta \to 0$ in order to prove Theorem 1.1. We consider initial data $\rho_0, q_0$ satisfying (1.23)-(1.25). Following the same ideas as in [17, Section 4.0], [33, Section 7.10.7], we construct $\{\rho_{0,\delta}\} \in L^{\beta}(\Omega)$ and $\{q_{0,\delta}\} \in L^{\frac{2\beta}{\beta+1}}(\Omega)$ such that

$$ \rho_{0,\delta} = \rho_{b,0} \text{ in } B_0, \quad \rho_{0,\delta} \to \rho_0 \text{ in } L^{\gamma f}(\Omega), \quad q_{0,\delta} \to q_0 \text{ in } L^{\frac{2\gamma f}{\gamma f+1}}(\Omega), $$

and

$$ E_{\delta}(\rho_{0,\delta}, q_{0,\delta}, \mathbb{1}_{B_0}) \to E(\rho_0, q_0, \mathbb{1}_{B_0}). $$

The energy estimate (2.8) yields, up to a subsequence,

$$ u_{\delta} \to u \quad \text{weakly in } L^2(0, T; H^1(\Omega)). \qquad (5.1) $$

By Proposition 3.3, up to a subsequence and for all $1 \le p < \infty$,

$$ \mathbb{1}_{B_{\delta}} \to \mathbb{1}_{B} \quad \text{weakly-}^* \text{ in } L^{\infty}((0, T) \times \Omega), \quad \text{strongly in } C([0, T]; L_{\text{loc}}^{p}(\mathbb{R}^3)), \qquad (5.2) $$

where the pair $(u, \mathbb{1}_B)$ satisfies

$$ \partial_t \mathbb{1}_B + u \cdot \nabla \mathbb{1}_B \text{ in } D'((0, T) \times \mathbb{R}^3), \quad B(0) = B_0. \qquad (5.3) $$

Moreover, the sequences $\{u_{\delta}\}$ and $\{\eta_{\delta}\}$ satisfy the assumptions of Proposition 3.1. It then follows that

$$ \eta_{\delta} \to \eta \quad \text{weakly in } H^1([0, T]; \mathcal{M}), \quad \text{strongly in } C([0, T]; \mathcal{M}), \qquad (5.4) $$

and, if we define $\mathcal{B}$ by

$$ \mathcal{B}(t) \stackrel{\text{def}}{=} \eta[t](\mathcal{B}_0) \quad \forall t \in [0, T], $$

then the velocity field is compatible with the system $\{\mathcal{B}, \eta\}$.

Writing $\mathcal{B} = \mathcal{B}(x_b, R_b)$ with $R_b, x_b \in H^1(0, T)$, we have

$$ R_b(t) \ge R_0/2 \quad \forall t \in [0, T]. $$
Moreover, from the strong convergence of $\eta_\delta$ in $C([0, T]; \mathcal{M})$, we deduce that for $\delta > 0$ sufficiently small,

$$ \mathcal{B}(t) \subset \mathcal{B}_{\delta}^{\frac{\sigma}{2}}(t) \quad \forall t \in [0, T]. $$

Since $\text{dist}(\mathcal{B}_\delta(t), \partial\Omega) \ge 3\sigma/2$ in $[0, T]$, it follows that

$$ \text{dist}(\mathcal{B}(t), \partial\Omega) \ge \sigma \quad \forall t \in [0, T]. $$

The energy estimate (2.14) also yields uniform bounds for $1_{\mathcal{B}_\delta}\rho_\delta$ in $L^\infty(0, T; L^{\gamma_f}(\Omega))$ and $(1 - 1_{\mathcal{B}_\delta})\rho_\delta$ in $L^\infty(0, T; L^{\gamma_b}(\Omega))$. Since $\rho_\delta$ is constant in $\mathcal{B}_\delta$, the uniform bound in $L^\infty(0, T; L^{\gamma_b}(\Omega))$ implies a uniform bound in $L^\infty((0, T) \times \Omega)$ and in particular in $L^\infty(0, T; L^{\gamma_f}(\Omega))$. We deduce that, up to a subsequence,

$$ \rho_{\delta} \to \rho \text{ weakly-}^* \text{ in } L^{\infty}(0, T; L^{\gamma_f}(\Omega)),. \quad (5.5) $$

Following the same arguments as in [17, Section 4.2], [33, Section 7.10.1], we also deduce from (2.8),

$$ \rho_{\delta} \to \rho \text{ in } C([0, T], L_{\text{weak}}^{\gamma_f}(\Omega)), $$

$$ \rho_{\delta} u_{\delta} \to \rho u \text{ weakly-}^* \text{ in } L^{\infty}(0, T; L^{\frac{2\gamma_f}{\gamma_f+1}}(\Omega)), $$

and $\rho, u$ solve the continuity equation in $D'((0, T) \times \mathbb{R}^3)$.

Similarly to the previous sections, we now focus on the convergence of the momentum equation (2.4). Let $\varphi \in \mathcal{T}(Q_b)$. As in Section 4.3 and Section 4.4, since $\bar{Q}_b$ is a compact subset of $(0, T) \times \Omega$, there exists $\sigma_0 > 0$ and $V_b, \omega_b, \Lambda_b \in \mathcal{D}(0, T)$ such that, for all $t \in [0, T]$,

$$ \varphi(t, \cdot) = \varphi_b(t, \cdot) \text{ in } \mathcal{B}^{\sigma_0}(t) \stackrel{\text{def}}{=} \mathcal{B}(x_b(t), R_b(t) + \sigma_0), $$

with $\varphi_b$ defined in (4.41). In particular, for $0 < \sigma \le \sigma_0$, we have the decomposition:

$$ \varphi = (1 - 1_{\mathcal{B}^\sigma})\varphi + 1_{\mathcal{B}^\sigma}\varphi_b. \quad (5.6) $$

Using once more the strong convergence of $\eta_\delta$ given in (5.4), we deduce that there exists $\delta_0 > 0$ such that for all $0 < \delta \le \delta_0$,

$$ \mathcal{B}_{\delta}(t) \subset \mathcal{B}^{\sigma_0}(t) \quad \forall t \in [0, T], $$

and

$$ 1_{\mathcal{B}_\delta} 1_{\mathcal{B}^{\sigma_0}} = 1_{\mathcal{B}_\delta}, \quad (1 - 1_{\mathcal{B}_\delta})(1 - 1_{\mathcal{B}^{\sigma_0}}) = (1 - 1_{\mathcal{B}^{\sigma_0}}), \quad 1_{\mathcal{B}_\delta}(1 - 1_{\mathcal{B}^{\sigma_0}}) = 0. $$

Similarly to the previous sections, the difficult term is the nonlinear convective term

$$ \int_0^T \int_{\Omega} \rho_{\delta} u_{\delta} \otimes u_{\delta} dy d\tau. \quad (5.7) $$

When details are omitted, the arguments for the convergence of (5.7) are analogous to those presented in Section 4.3. For $0 < \sigma \le \sigma_0$, we write

$$ \int_0^T \int_{\Omega} \rho_{\delta} u_{\delta} \otimes u_{\delta} dy d\tau - \int_0^T \int_{\Omega} \rho u \otimes u dy d\tau = \sum_{i=1}^4 M_i(\sigma, \delta) $$

where

$$ M_1(\sigma, \delta) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} (1 - 1_{\mathcal{B}^\sigma}) \rho_\delta u_\delta \otimes u_\delta : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_{\Omega} (1 - 1_{\mathcal{B}^\sigma}) \rho u \otimes u : \mathcal{D}(\varphi) dy d\tau, $$

$$ M_2(\sigma, \delta) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} (1_{\mathcal{B}^\sigma} - 1_{\mathcal{B}}) \rho_\delta u_\delta \otimes u_\delta : \mathcal{D}(\varphi) dy d\tau + \int_0^T \int_{\Omega} (1_{\mathcal{B}^\sigma} - 1_{\mathcal{B}}) \rho u \otimes u : \mathcal{D}(\varphi) dy d\tau, $$

$$ M_3(\sigma, \delta) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} (1_{\mathcal{B}} - 1_{\mathcal{B}_\delta}) \rho_\delta u_\delta \otimes u_\delta : \mathcal{D}(\varphi) dy d\tau, $$

$$ M_4(\sigma, \delta) \stackrel{\text{def}}{=} \int_0^T \int_{\Omega} 1_{\mathcal{B}_\delta} \rho_\delta u_\delta \otimes u_\delta : \mathcal{D}(\varphi) dy d\tau - \int_0^T \int_{\Omega} 1_{\mathcal{B}} \rho u \otimes u : \mathcal{D}(\varphi) dy d\tau. $$
Let $\theta > 0$, there exists $\sigma_1, \delta_1 > 0$ such that for $0 < \sigma \le \min\{\sigma_0, \sigma_1\}$ and $0 < \delta \le \min\{\delta_0, \delta_1\}$,

$$|M_2(\sigma, \delta)| \le \theta/4, \quad |M_3(\sigma, \delta)| \le \theta/4,$$

For $M_4(\sigma, \delta)$, we can prove once again an analogue of Lemma 4.2, that is for $i \in \{1, 2, 3\}$,

$$(\rho_{\delta} u_{\delta}, 1_{\mathcal{B}_{\delta}} e_{i})_{\Omega} \xrightarrow{\delta \to 0} (\rho u, 1_{\mathcal{B}} e_{i})_{\Omega} \text{ in } \mathcal{C}([0, T]), \quad (5.8)$$

$$(\rho_{\delta} u_{\delta}, 1_{\mathcal{B}_{\delta}} e_{i} \times (x - x_{\delta}))_{\Omega} \xrightarrow{\delta \to 0} (\rho u, 1_{\mathcal{B}} e_{i} \times (x - x_{b}))_{\Omega} \text{ in } \mathcal{C}([0, T]), \quad (5.9)$$

$$(\rho_{\delta} u_{\delta}, 1_{\mathcal{B}_{\delta}}(x - x_{\delta}))_{\Omega} \xrightarrow{\delta \to 0} (\rho u, 1_{\mathcal{B}}(x - x_{b}))_{\Omega} \text{ in } \mathcal{C}([0, T]). \quad (5.10)$$

The proof of (5.8)-(5.10) follows the same arguments as the proof of Lemma 4.1 with the estimate (4.51) replaced by

$$\int_{t_1}^{t_2} \|\rho_{\delta} u_{\delta} \otimes u_{\delta}\|_{L^{\frac{6\gamma_f}{4\gamma_f+3}}(\Omega)} \|\nabla \varphi_{\delta,i}^{\nu}\|_{L^{\frac{6\gamma_f}{2\gamma_f-3}}(\Omega)} + (\frac{4}{3}\mu_f + \nu_f)\|\nabla u_{\delta}\|_{L^2(\Omega)}\|\nabla \varphi\|_{L^2(\Omega)} + \nu_b\|\nabla u_{\delta}\|_{L^2(\Omega)}\|\nabla \varphi\|_{L^2(\Omega)} \\ \|\rho_{\delta}\|_{L^{\gamma_f+\theta}(\Omega)}^{\gamma_f} \|\nabla \varphi_{\delta,i}^{\nu}\|_{L^{\frac{\gamma_f+\theta}{\theta}}(\Omega)} + \delta\|\rho_{\delta}\|_{L^{\beta+\theta}(\Omega)}^{\beta} \|\nabla \varphi_{\delta,i}^{\nu}\|_{L^{\frac{\beta+\theta}{\theta}}(\Omega)} + \|\rho_{\delta}\|_{0,\beta+\theta}\|g\|_{L^{\infty}(\Omega)}\|\varphi_{\delta,i}^{\nu}\|_{L^{\frac{\beta+\theta}{\theta}}(\Omega)} + \frac{2\kappa_b}{R_{b,0}}\|\nabla \varphi_{\delta,i}^{\nu}\|_{L^1(\Omega)}$$

From the convergences (5.8)-(5.10), we obtain

$$\lim_{\delta \to 0} \int_0^T \int_{\Omega} \mathrm{div}(\varphi_b) 1_{\mathcal{B}_\delta} \rho_{\delta} |u_{\delta}|^2 \mathrm{d}y \mathrm{d}\tau \to \int_0^T \int_{\Omega} \mathrm{div}(\varphi_b) 1_{\mathcal{B}} \rho |u|^2 \mathrm{d}y \mathrm{d}\tau,$$

and consequently, there exists $\delta_2 > 0$ such that for all $0 < \delta \le \delta_2$,

$$|M_4(\sigma, \delta)| \le \theta/4.$$

Finally, setting $\sigma_m = \{\sigma_0, \sigma_1\}$ and following similar arguments to the one presented in [14, Section 8] and already presented in Section 4.3-Section 4.4, we get

$$\rho_{\delta} u_{\delta} \otimes u_{\delta} \to \rho u \otimes u \text{ in } L^2(Q_f^{\sigma_m})$$

where
$$Q_f^{\sigma_m} \stackrel{\text{def}}{=} (0, T) \times \Omega \setminus \overline{Q_b^{\sigma_m}}, \quad \overline{Q_b^{\sigma_m}} = \{(t, x) \in (0, T) \times \Omega \mid x \in \overline{B^{\sigma_m}(t)}\}.$$

$$\|\rho_{\delta}\|_{L^{\gamma_f+\theta}(\mathcal{K}_f)} + \delta\|\rho_{\delta}\|_{L^{\beta+\theta}(\mathcal{K}_f)}^{\beta} \le c (\mathcal{K}_f, E_{\delta}(\rho_{0,\delta}, q_{0,\delta}, 1_{\mathcal{B}_0}), \|g\|_{L^{\infty}((0,T)\times\Omega)}) \quad (5.11)$$

for any compact $\mathcal{K}_f \subset Q_f$, where $\theta > 0$ is a positive constant independent of $\delta$. The proof of (5.11) follows the same arguments as in [14, Equation (9.5)]. In particular, we can apply [18, Proposition 2.3], [17, Lemma 4.1], [33, Lemma 7.52] to obtain

$$a_f \rho_{\delta}^{\gamma_f} \to a_f \overline{\rho^{\gamma_f}} \text{ weakly in } L^{\frac{\gamma_f+\theta}{\gamma_f}}(\mathcal{K}_f) \text{ for any compact } \mathcal{K}_f \subset Q_f,$$

where the limit functions satisfy the identity

$$\int_0^T \int_{\mathbb{R}^3} (\rho u) \cdot \partial_{\tau} \varphi + (\rho u \otimes u) : \mathbb{D}(\varphi) + a_f \overline{\rho^{\gamma_f}} \mathrm{div}(\varphi) \mathrm{d}y \mathrm{d}\tau = \int_0^T \int_{\mathbb{R}^3} \mathbb{T}_f(u) : \mathbb{D}(\varphi) - \overline{\rho g} \cdot \varphi \mathrm{d}y \mathrm{d}\tau$$

for all test functions $\varphi \in \mathcal{D}(Q_f)$. Similarly to Section 4.4, the proof will be completed by showing the strong convergence of the density. However, unlike Section 4.4, we do not know a priori that the renormalized equation (1.28) holds, as we do not have $\rho \in L^2(Q_T)$. This difficulty dissapears if $\gamma \ge 9/5$ as $L^{\gamma+\theta}(\gamma)(Q_T) \subset L^2(Q_T)$. To prove the strong convergence of the density, we proceed as in [14, Equation (9.9)], employing [13, Proposition 4.3], [17, Lemma 4.2], [33, Lemma 7.55] to derive an analogue of Lemma 4.4, namely,

$$\lim_{\delta \to 0} \int_0^T \int_{\mathbb{R}^3} \varphi \left( a_f \rho_{\delta}^{\gamma_f} - \left(\frac{4}{3}\mu_f + \nu_f\right) \mathrm{div}(u_{\delta}) \right) T_k(\rho_{\delta}) \mathrm{d}y \mathrm{d}\tau = \int_0^T \int_{\mathbb{R}^3} \varphi \left( a_f \overline{\rho^{\gamma_f}} - \left(\frac{4}{3}\mu_f + \nu_f\right) \mathrm{div}(u) \right) \overline{T_k(\rho)} \mathrm{d}y \mathrm{d}\tau \quad (5.12)$$

for any $\varphi \in \mathcal{D}(Q_f)$, where $T_k(\rho) = \min\{\rho, k\}$ are cut-off operators, and $\overline{T_k(\rho)}$ stands for the weak limit of $T_k(\rho_{\delta})$. As in Section 4.4, the relation (5.12) gives rise to the inequality

$$\overline{T_k(\rho)} \mathrm{div}(u) \ge \overline{T_k(\rho)} \mathrm{div}(u) \text{ on } (0, T) \times \mathbb{R}^3, \quad (5.13)$$
and proceeding as in [12, Proposition 6.1], we can prove that for any compact set $\mathcal{K}_f \subset Q_f$ such that $|\mathcal{K}_f| \le M$, then

$$ \sup_{k \ge 1} \left( \lim_{\delta \to 0} \int_{\mathcal{K}_f} |T_k(\rho_\delta) - T_k(\rho)|^{\gamma_f+1} dy d\tau \right) \le c \left( M, \sup_{\delta>0} \|\nabla u_\delta\|_{L^2(\mathcal{K}_f)} \right). $$

Let $\mathcal{K}_b$ a compact subset of $Q_b$, by Lemma 3.4, we get $\rho_\delta \to \rho$ strongly in $L^1(\mathcal{K}_b)$. In particular,

$$ \sup_{k \ge 1} \left( \lim_{\delta \to 0} \int_{\mathcal{K}_b} |T_k(\rho_\delta) - T_k(\rho)|^{\gamma_f+1} dy d\tau \right) = 0. \qquad (5.14) $$

Let $\mathcal{K} \subset (0, T) \times \Omega$ compact such that $|\mathcal{K}| \le M$. Let $\mathcal{U}(k)$ a neighborhood of $\partial Q_b$ such that $|\mathcal{U}(k)| \le k^{-(\gamma_f+1)}$, we have

$$ \int_{\mathcal{K}} |T_k(\rho_\delta) - T_k(\rho)|^{\gamma_f+1} dy d\tau \le \int_{(\mathcal{K} \setminus \mathcal{U}(k)) \cap Q_f} |T_k(\rho_\delta) - T_k(\rho)|^{\gamma_f+1} + \int_{\mathcal{U}(k)} (2k)^{\gamma_f+1} + \int_{(\mathcal{K} \setminus \mathcal{U}(k)) \cap Q_b} |T_k(\rho_\delta) - T_k(\rho)|^{\gamma_f+1}. $$

Consequently

$$ \sup_{k \ge 1} \left( \lim_{\delta \to 0} \int_{\mathcal{K}} |T_k(\rho_\delta) - T_k(\rho)|^{\gamma_f+1} dy d\tau \right) \le c \left( M, \sup_{\delta>0} \|\nabla u_\delta\|_{L^2(\mathcal{K}_f)} \right) + 2^{\gamma_f+1}. $$

We finally obtain

$$ \sup_{k \ge 1} \left( \lim_{\delta \to 0} \int_{(0,T) \times \mathbb{R}^3} |T_k(\rho_\delta) - T_k(\rho)|^{\gamma_f+1} dy d\tau \right) < \infty. \qquad (5.15) $$

By virtue [12, Proposition 7.1], [33, Lemma 7.57], boundedness of the deflection measure (5.15) implies that $\rho, u$ satisfy the renormalized continuity equation (1.28) where $b$ satisfies (1.29) with $-1 < \lambda_1 \le \frac{s(\gamma_f)}{2} - 1$ and $s(\gamma_f) = \frac{5\gamma_f - 3}{3}$. Based on this result, introducing

$$ L_k(s) = \begin{cases} s \log(s), & s \in [0, k), \\ s \log(k) + s - k, & s \in [k, \infty), \end{cases} $$

we can prove that

$$ \partial_t (\overline{L_k(\rho)} - L_k(\rho)) + \operatorname{div} ((\overline{L_k(\rho)} - L_k(\rho))u) + (\overline{T_k(\rho)} \operatorname{div}(u) - T_k(\rho) \operatorname{div}(u)) = 0 $$

which gives together with (5.13)

$$ \sup_{k \ge 0} \left( \lim_{\delta \to 0} \| T_k(\rho_\delta) - T_k(\rho) \|_{L^{\gamma_f+1}((0,T) \times \mathbb{R}^3)} \right) \le 0. $$

This is enough to deduce that $\rho_\delta \to \rho$ in $\mathcal{C}([0, T], L^1(\mathbb{R}^3))$ for instance. This completes the proof of Theorem 1.1, cf [12, Lemma 4.3, 4.4],[18, Section 8], [17, Section 4.5, 4.6], [33, Lemma 7.60].

## Acknowledgement

The authors acknowledge the support of ANR project MSMΦ ANR-23-CE40-0013-01. The authors would like to thank Hélène Mathis and Nicolas Seguin for fruitful discussions when preparing the paper.

## References

[1] Z. Biro and J. J. L. Velazquez. Analysis of a free boundary problem arising in bubble dynamics. *SIAM Journal on Mathematical Analysis*, 32(1):142–171, Jan. 2000.

[2] C. Bost, G.-H. Cottet, and E. Maitre. Convergence analysis of a penalization method for the three-dimensional motion of a rigid body in an incompressible viscous fluid. *SIAM Journal on Numerical Analysis*, 48(4):1313–1337, Jan. 2010.

[3] M. Boulakia. Existence of weak solutions for an interaction problem between an elastic structure and a compressible viscous fluid. *Journal de Mathématiques Pures et Appliquées*, 84(11):1515–1554, Nov. 2005.

[4] M. Boulakia. Existence of weak solutions for the three-dimensional motion of an elastic structure in an incompressible fluid. *Journal of Mathematical Fluid Mechanics*, 9(2):262–294, Sept. 2006.

[5] M. Boulakia, E. L. Schwindt, and T. Takahashi. Existence of strong solutions for the motion of an elastic structure in an incompressible viscous fluid. *Interfaces and Free Boundaries, Mathematical Analysis, Computation and Applications*, 14(3):273–306, Nov. 2012.
[6] C. Burtea and D. Gérard-Varet. Weak solutions of a viscous model for fluid-bubbles interaction, 2025.

[7] N. V. Chemetov, Š. Nečasová, and B. Muha. Weak-strong uniqueness for fluid-rigid body interaction problem with slip boundary condition. *Journal of mathematical physics*, 60(1), 2019.

[8] B. Desjardins and M. Esteban. On weak solutions for fluid-rigid structure interaction: Compressible and incompressible models. *Communications in Partial Differential Equations*, 25(7-8):263-285, Jan. 1999.

[9] B. Desjardins and M. J. Esteban. Existence of weak solutions for the motion of rigid bodies in a viscous fluid. *Archive for Rational Mechanics and Analysis*, 146(1):59-71, Apr. 1999.

[10] B. Desjardins, M. J. Esteban, C. Grandmont, and P. Le Tallec. Weak solutions for a fluid-elastic structure interaction model. *Revista Matemática Complutense*, 14(2):523-538, 2001.

[11] L. C. Evans. *Partial differential equations*, volume 19. American mathematical society, 2022.

[12] E. Feireisl. On compactness of solutions to the compressible isentropic navier-stokes equations when the density is not square integrable. *Commentationes Mathematicae Universitatis Carolinae*, 42(1):83-98, 2001.

[13] E. Feireisl. Viscous and/or heat conducting compressible fluids. *Handbook of mathematical fluid dynamics*, 1:307-371, 2002.

[14] E. Feireisl. On the motion of rigid bodies in a viscous compressible fluid. *Archive for Rational Mechanics and Analysis*, 167(4):281-308, May 2003.

[15] E. Feireisl. *On the motion of rigid bodies in a viscous incompressible fluid*, pages 419-441. Birkhäuser Basel, 2003.

[16] E. Feireisl. *Dynamics of viscous compressible fluids*, volume 26. Oxford University Press, 2004.

[17] E. Feireisl, A. Novotný, and H. Petzeltová. On the existence of globally defined weak solutions to the navier—stokes equations. *Journal of Mathematical Fluid Mechanics*, 3(4):358-392, Nov. 2001.

[18] E. Feireisl, A. Novotný, and H. Petzeltová. On the domain dependence of solutions to the compressible navier-stokes equations of a barotropic fluid. *Mathematical Methods in the Applied Sciences*, 25(12):1045-1073, July 2002.

[19] C. Grandmont and Y. Maday. Existence for an unsteady fluid-structure interaction problem. *ESAIM: Mathematical Modelling and Numerical Analysis*, 34(3):609-636, May 2000.

[20] S. Guerrero and M. Boulakia. A regularity result for a solid-fluid system associated to the compressible navier-stokes equations. *Annales de l'Institut Henri Poincaré C, Analyse non linéaire*, 26(3):777-813, June 2009.

[21] D. Gérard-Varet and M. Hillairet. Existence of weak solutions up to collision for viscous fluid-solid systems with slip. *Communications on Pure and Applied Mathematics*, 67(12):2022-2076, May 2014.

[22] B. H. Haak, D. Maity, T. Takahashi, and M. Tucnak. Mathematical analysis of the motion of a rigid body in a compressible navier-stokes-fourier fluid. *Mathematische Nachrichten*, 292(9):1972-2017, May 2019.

[23] M. Hieber and M. Murata. The $l^p$-approach to the fluid-rigid body interaction problem for compressible fluids. *Evolution Equations and Control Theory*, 4(1):69-87, Feb. 2015.

[24] M. Hillairet, H. Mathis, and N. Seguin. Analysis of compressible bubbly flows. part i: construction of a microscopic model. *ESAIM: Mathematical Modelling and Numerical Analysis*, 57(5):2835-2863, Sept. 2023.

[25] M. Hillairet, H. Mathis, and N. Seguin. Analysis of compressible bubbly flows. part ii: Derivation of a macroscopic model. *ESAIM: Mathematical Modelling and Numerical Analysis*, 57(5):2865-2906, Sept. 2023.

[26] K.-H. Hoffmann and V. N. Starovoitov. *On a motion of a solid body in a viscous fluid: two-dimensional case*. Fakultät für Mathematik, Technische Universität München, 1996.

[27] K.-H. Hoffmann and V. N. Starovoitov. Zur bewegung einer kugel in einer zähen flüssigkeit. *Doc. Math*, 5(15-21):30, 2000.

[28] N. Judakov. The solvability of the problem of the motion of a rigid body in a viscous incompressible fluid. *Dinamika Splošn. Sredy*, 18(249-253):255, 1974.

[29] O. Kreml, v. Nečasová, and T. Piasecki. Weak-strong uniqueness for the compressible fluid-rigid body interaction. *Journal of Differential Equations*, 268(8):4756-4785, Apr. 2020.
[30] P. Lions. Mathematical topics in fluid dynamics, vol. 2, compressible models oxford university press, 1998.

[31] P.-L. Lions. *Mathematical topics in fluid mechanics: volume 2: compressible models*, volume 2. oxford university press, 1996.

[32] v. Nečasová, M. Ramaswamy, A. Roy, and A. Schlömerkemper. Motion of a rigid body in a compressible fluid with navier-slip boundary condition. *Journal of Differential Equations*, 338:256-320, Nov. 2022.

[33] A. Novotný and I. Straškraba. *Introduction to the Mathematical Theory of Compressible Flow*. Oxford University Press Oxford, June 2004.

[34] M. S. Plesset. The dynamics of cavitation bubbles. *Journal of Applied Mechanics*, 16(3):277-282, Sept. 1949.

[35] M. S. Plesset and A. Prosperetti. Bubble dynamics and cavitation. *Annual review of fluid mechanics*, 9:145-185, 1977.

[36] L. Rayleigh. Viii. on the pressure developed in a liquid during the collapse of a spherical cavity. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, 34(200):94-98, Aug. 1917.

[37] A. Roy and T. Takahashi. Stabilization of a rigid body moving in a compressible viscous fluid. *Journal of Evolution Equations*, 21(1):167-200, Apr. 2020.

[38] J. A. San Martín, V. Starovoitov, and M. Tucsnak. Global weak solutions for the two-dimensional motion of several rigid bodies in an incompressible viscous fluid. *Archive for Rational Mechanics and Analysis*, 161(2):113-147, Feb. 2002.

[39] D. Serre. Chute libre d'un solide dans un fluide visqueux incompressible. existence. *Japan Journal of Applied Mathematics*, 4(1):99-110, Feb. 1987.

[40] T. Takahashi. Analysis of strong solutions for the equations modeling the motion of a rigid-fluid system in a bounded domain. *Advances in Differential Equations*, 8(12), Jan. 2003.

[41] T. Takahashi and M. Tucsnak. Global strong solutions for the two-dimensional motion of an infinite cylinder in a viscous fluid. *Journal of Mathematical Fluid Mechanics*, 6(1):53-77, Mar. 2004.