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

HAL
open science

# A well-posedness result for the compressible two-fluid model with density-dependent viscosity

Sagbo Marcel Zodji

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

## To cite this version:

Sagbo Marcel Zodji. A well-posedness result for the compressible two-fluid model with density-dependent viscosity. 2023. hal-04586675

HAL Id: hal-04586675

https://hal.science/hal-04586675

Preprint submitted on 24 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A well-posedness result for the compressible two-fluid model with density-dependent viscosity

Sagbo Marcel Zodji

► To cite this version:

Sagbo Marcel Zodji. A well-posedness result for the compressible two-fluid model with density-dependent viscosity. Journal of Mathematical Fluid Mechanics, 2023, 27 (3), pp.50. 10.1007/s00021-025-00954-y. hal-04586675v2

HAL Id: hal-04586675

https://hal.science/hal-04586675v2

Submitted on 22 Jul 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A WELL-POSEDNESS RESULT FOR THE COMPRESSIBLE TWO-FLUID MODEL WITH DENSITY-DEPENDENT VISCOSITY

SAGBO MARCEL ZODJI

**Abstract.** In this paper, we study a system of PDE's describing the motion of two compressible viscous fluids occupying the whole space $\mathbb{R}^d$ ($d \in \{2, 3\}$). We prove a local well-posedness result with discontinuous initial data: we require the initial density to be Hölder continuous on both sides of a hyper-surface across which it can experience large jumps. The initial velocity belongs to the Sobolev space $H^1(\mathbb{R}^d)$ and the divergence of the initial stress tensor to belong to $L^2(\mathbb{R}^d)$. The later assumption expresses somehow the continuity of the stress tensor. This result is more general than the one by Tani [26] by allowing low regular initial data, and it can be used to serve as building blocks for global solutions.

## 1. INTRODUCTION AND MAIN RESULT

1.1. **Presentation of the model.** In this paper, we are concerned with the following system of PDE's describing the evolution of two compressible fluids with density-dependent viscosity in the whole space $\mathbb{R}^d$, $d \in \{2, 3\}$,

$$ (1.1) \qquad \left\{ \begin{array}{l} \partial_t \rho_{\pm} + \text{div}(\rho_{\pm} u) = 0, \\ \partial_t (\rho u^j) + \text{div}(\rho u^j u) = \text{div} \, \Pi^j. \end{array} \right. $$

Above, $\rho_{\pm} = \rho_{\pm}(t, x) \ge 0$ is the density of the phase fluid labelled $\pm$, $\rho := \rho_+ + \rho_-$ is the total density of the mixture and $u = u(t, x) \in \mathbb{R}^d$ is the single velocity of the mixture. The stress tensor $\Pi$ is given by

$$ (1.2) \qquad \Pi^{jk} = 2\mu \mathbb{D}^{jk} u + (\lambda \text{div } u - P + \tilde{P}) \delta^{jk} \quad \text{with} \quad 2\mathbb{D}^{jk} u = \partial_j u^k + \partial_k u^j, $$

where $P = P(\rho_+, \rho_-)$, $\mu = \mu(\rho_+, \rho_-)$, $\lambda = \lambda(\rho_+, \rho_-)$ are respectively pressure, dynamic and kinematic viscosity of the mixture and are given $\mathcal{C}^2$ function of $\rho_{\pm}$ and finally $\tilde{P}$ is a positive equilibrium pressure. We refer to [13] for the derivation of the system (1.1) for the more general case of a mixture of two heat conducting fluids. It is assumed that the two pure components of the mixture are immiscible, so that each fluid is implicitly assumed to occupy its own domain and that $\rho_+\rho_- = 0$. In particular they are separated by a sharp interface $C$. The relevant physical situation for multi-phase flows is exactly where the total density $\rho$ experiences discontinuity across the interface $C$. Since the two densities $\rho_{\pm}$ of the mixture are assumed to be different, we see that the class of piecewise regular functions is suitable for the study of this system. This motivates the following notation.

**Notation 1.1.** Let $C = \{\gamma(s) \in \mathbb{R}^d, s \in \mathbb{T}^{d-1}\}$ be a hyper-surface that is the boundary of an open, simply connected set $D$ such that

$$ \mathbb{R}^d = D \cup C \cup (\mathbb{R}^d \setminus \bar{D}). $$

We define the set of piecewise $\alpha$-Hölder functions with respect to $\gamma$ as follows:

$$ \dot{\mathcal{C}}_{pw,\gamma}^{\alpha}(\mathbb{R}^d) := \dot{\mathcal{C}}^{\alpha}(\bar{D}) \cap \dot{\mathcal{C}}^{\alpha}(\mathbb{R}^d \setminus D) $$

and the non-homogeneous Hölder space is:

$$ \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d) := L^{\infty}(\mathbb{R}^d) \cap \dot{\mathcal{C}}_{pw,\gamma}^{\alpha}(\mathbb{R}^d). $$

This space strictly contains the Hölder space $\mathcal{C}^{\alpha}(\mathbb{R}^d)$. In particular, there exists functions in $\mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)$ that have discontinuities across the interface $C$. Moreover, we define the following space: for $I \subset \mathbb{R}$,

$$ L^p(I, \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)) := \left\{ g = g(t, x) : \left\{ \begin{array}{ll} \displaystyle \int_I \|g(t)\|_{\mathcal{C}_{pw,\gamma(t)}^{\alpha}(\mathbb{R}^d)}^p dt < \infty & \text{if } 1 \le p < \infty, \\ \displaystyle \sup_{t \in I} \|g(t)\|_{\mathcal{C}_{pw,\gamma(t)}^{\alpha}(\mathbb{R}^d)} < \infty & \text{if } p = \infty, \end{array} \right. \right\} $$

<footnote>Date: December 12, 2023.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 76N10, 35A02, 76T10, 35Q30.</footnote>

<footnote>Key words and phrases. Two-fluid model, Compressible Navier-Stokes equations, Density-dependent viscosity, Free boundary problem, Lagrangian formulation.</footnote>
# A WELL-POSEDNESS RESULT FOR THE COMPRESSIBLE TWO-FLUID MODEL WITH DENSITY-DEPENDENT VISCOSITY

SAGBO MARCEL ZODJI

**Abstract.** In this paper, we study a system of PDEs describing the motion of two compressible viscous fluids occupying the whole space $\mathbb{R}^d$ ($d \in \{2, 3\}$). The two phases of the mixture are separated by a $\mathcal{C}^{1+\alpha}$-regular sharp interface $\mathcal{C}$ across which the density can experience jumps. We prove the existence of a unique local-in-time solution assuming that the initial density is $\alpha$-Hölder continuous on both sides of $\mathcal{C}$. The initial velocity belongs to the Sobolev space $H^1(\mathbb{R}^d)$, and the divergence of the initial stress tensor belongs to $L^2(\mathbb{R}^d)$. The later assumption expresses somehow the continuity of the stress tensor. This result is more general than the one by Tani [32], as it allows for less regular initial data and furthermore it can serve as a building block for the construction of global-in-time solutions.

## 1. INTRODUCTION AND MAIN RESULT

1.1. **Presentation of the model.** In this paper, we are concerned with the following system of PDE describing the evolution of two compressible fluids with density-dependent viscosity in the whole space $\mathbb{R}^d$ ($d \in \{2, 3\}$),

$$ (1.1) \qquad \begin{cases} \partial_t \rho_{\pm} + \operatorname{div}(\rho_{\pm} u) = 0, \\ \partial_t (\rho u^j) + \operatorname{div}(\rho u^j u) = \operatorname{div} \Pi^j. \end{cases} $$

Above, $\rho_{\pm} = \rho_{\pm}(t, x) \ge 0$ is the density of the phase fluid labelled $\pm$, $\rho := \rho_+ + \rho_-$ is the total density of the mixture and $u = u(t, x) \in \mathbb{R}^d$ is the single velocity of the mixture. The stress tensor $\Pi$ is given by

$$ (1.2) \qquad \Pi^{jk} = 2\mu \mathbb{D}^{jk} u + (\lambda \operatorname{div} u - P + \tilde{P}) \delta^{jk} \quad \text{with} \quad \mathbb{D}^{jk} u = \frac{1}{2} (\partial_j u^k + \partial_k u^j), $$

where $P = P(\rho_+, \rho_-)$, $\mu = \mu(\rho_+, \rho_-)$, $\lambda = \lambda(\rho_+, \rho_-)$ are respectively pressure, dynamic and kinematic viscosity of the mixture and are given $\mathcal{C}^2$-regular functions of $\rho_{\pm}$, and finally $\tilde{P}$ is a positive equilibrium pressure. We refer to [16] for the derivation of (1.1) for the more general case of heat conducting fluids mixture.

It is assumed that the two pure components of the mixture are immiscible, so that each fluid is implicitly assumed to occupy its own domain and that $\rho_+\rho_- = 0$. In particular they are separated by a sharp interface $\mathcal{C}$. The relevant physical situation for multi-phase flows is exactly where the total density $\rho$ experiences discontinuity across the interface $\mathcal{C}$. Since the two densities $\rho_{\pm}$ of the mixture are assumed to be different, we see that the class of piecewise regular functions is suitable for the study of this system. This motivates the following definitions.

### 1.1.1. Piecewise regularity.

(1) Throughout this paper $\alpha \in (0, 1)$, and we assume that $\mathcal{C}$ is a compact $\mathcal{C}^{1+\alpha}$-regular manifold of dimension $d-1$, which is the boundary of an open, bounded and simply connected domain $D \subset \mathbb{R}^d$. It is defined as follows.

There exists a finite number $J \in \mathbb{N}^*$ of open bounded sets $V_j \subset \mathbb{R}^{d-1}$, $j \in [1, J]$ and maps

$$ \gamma_j: V_j \to \mathcal{C} \subset \mathbb{R}^d, \quad \mathcal{C} = \bigcup_{j=1}^{J} \gamma_j(V_j), $$

such that $\gamma_j \in \mathcal{C}^{1+\alpha}(V_j)$ are homeomorphisms on their images. For $(j, k) \in [1, J] \times [1, J]$ such that $\gamma_j(V_j) \cap \gamma_k(V_k) \ne \emptyset$, the map

$$ \gamma_j^{-1} \circ \gamma_k: \gamma_k^{-1}(\gamma_j(V_j) \cap \gamma_k(V_k)) \to \gamma_j^{-1}(\gamma_j(V_j) \cap \gamma_k(V_k)) $$

<footnote>Date: July 22, 2025.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 76N10, 35A02, 76T10, 35Q30.</footnote>

<footnote>Key words and phrases. Two-fluid model, Compressible Navier-Stokes equations, Density-dependent viscosity, Free boundary problem, Lagrangian formulation.</footnote>
for a time dependent hyper-surface $\mathcal{C} = \mathcal{C}(t)$. For $I = (0, T)$, the corresponding space will sometimes be shortly denoted by $L_T^p \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)$ for finite $T$ and $L^p \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)$ for $T = \infty$. Alternatively, the local-in-time space is defined as follows:

$$L_{\text{loc}}^{p}([0, \infty), \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^{d})) := \{g \in L^{p}(I, \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^{d})), \quad \forall I \subset [0, \infty) \text{ bounded}\}.$$

A suitable way to rewrite the system (1.1) is to introduce the so-called fractional density or concentration $c$ defined by

$$c := \frac{\rho_{+}}{\rho}.$$

From the first two equations of (1.1), one checks easily, assuming some regularity, that the concentration solves the transport equation

$$(1.3) \qquad \partial_t c + u \cdot \nabla c = 0.$$

Let us notice that the domain occupied by the phase $\pm$ can be defined by:

$$D_{\pm}(t) := \{x \in \mathbb{R}^d : \rho_{\pm}(t, x) > 0\}.$$

As well, one can write the domain occupied by the one phase fluid with the help of the concentration as follows:

$$(1.4) \qquad D_{+}(t) = \{x \in \mathbb{R}^d : c(t, x) = 1\}, \quad D_{-}(t) = \{x \in \mathbb{R}^d : c(t, x) = 0\}$$

and consequently the interface is:

$$\mathcal{C}(t) := \partial D_{+}(t) = \partial D_{-}(t).$$

Assume that the velocity is regular, for instance Lipschitz, then it admits a continuous flow map $X(t)$ that helps together with (1.3) and (1.4) to obtain that the initial interface is just transported by the fluid

$$\mathcal{C}(t) = X(t)\mathcal{C}(0).$$

On the other hand, by summing the mass equation on the density of the one phase $\rho_{\pm}$ (1.1)$_1$, one obtains that the total density solves a mass equation with velocity $u$. Besides, the density $\rho_{\pm}$ of the one phase fluid can be written in terms of the total density and the concentration as follows:

$$\rho_{+} = c\rho \quad \text{and} \quad \rho_{-} = (1-c)\rho,$$

thus, the pressure can be written as function of the total density and the concentration $P = P(\rho, c)$ such that $P_{+}(\rho) = P(\rho, 1)$ and $P_{-}(\rho) = P(\rho, 0)$ as well as the viscosity. Hence, formally $(c, \rho, u)$ solves:

$$(1.5) \qquad \begin{cases} \partial_t c + u \cdot \nabla c = 0, \\ \partial_t \rho + \text{div}(\rho u) = 0, \\ \partial_t (\rho u) + \text{div}(\rho u \otimes u) + \nabla P(\rho, c) = \text{div}(2\mu(\rho, c)\mathbb{D}u) + \nabla(\lambda(\rho, c) \text{div } u). \end{cases}$$

Conversely, if $(c, \rho, u)$ solves (1.5), one can set

$$\rho_{+} = c\rho \quad \text{and} \quad \rho_{-} = (1-c)\rho$$

and checks easily that $\rho_{\pm}$ solve the equations (1.1)$_1$. Just by replacing

$$\rho = \rho_{+} + \rho_{-} \quad \text{and} \quad c = \frac{\rho_{+}}{\rho_{+} + \rho_{-}}$$

in the pressure and the viscosity, one notices that $(\rho_{+}, \rho_{-}, u)$ solves (1.1). In fact, we obtain that $(\rho_{+}, \rho_{-}, u)$ solves (1.1) if and only if $(c, \rho, u)$ solves (1.5) owing to some change of pressure and viscosity law, at least for classical solutions. The equivalence also holds true for weak solutions à la Lions due to the result in [27, Lemma 2.5], which is a version of DiPerna-Lions theory for renormalized solutions. Thus, we will only work with the system (1.5). The equations (1.5) are supplemented with initial data

$$(1.6) \qquad c|_{t=0} = c_0, \quad \rho|_{t=0} = \rho_0 \quad \text{and} \quad u|_{t=0} = u_0,$$

with the following behavior at infinity

$$(1.7) \qquad \lim_{|x| \to \infty} \rho_0(x) = \tilde{\rho}, \quad \lim_{|x| \to \infty} P(\rho_0(x), c_0(x)) = \tilde{P}, \quad \lim_{|x| \to \infty} \mu(\rho_0(x), c_0(x)) = \tilde{\mu} \quad \text{and} \quad \lim_{|x| \to \infty} \lambda(\rho_0(x), c_0(x)) = \tilde{\lambda},$$

where $\tilde{\rho}$, $\tilde{P}$, $\tilde{\mu}$ and $\tilde{\lambda}$ are positive constants.

The purpose of this paper is to construct a local-in-time solution $(c, \rho, u)$ of the equations (1.5) with discontinuous initial data and with density-dependent viscosity: the initial density $\rho_0$ is assumed to be Hölder continuous
is a $\mathscr{C}^{1+\alpha}$-diffeomorphism. Furthermore, these mappings provide an analytical expression for the outward normal vector field $n_x$ on $\mathcal{C}$. We introduce the following notations:

$$ \|\mathcal{C}\|_{\text{Lip}} := \sup_j \|\nabla \gamma_j\|_{L^\infty(V_j)}, \quad \text{and} \quad \|\mathcal{C}\|_{\mathcal{C}^{1+\alpha}} := \sup_j \|\nabla \gamma_j\|_{\mathcal{C}^{\alpha}(V_j)}. $$

Following [11, page 6], we define the quantity:

$$ (1.3) \qquad |\mathcal{C}|_* := \left( \inf_j \left( \inf_{\substack{s \neq s' \\ s, s' \in V_j}} \frac{|\gamma_j(s) - \gamma_j(s')|}{|s - s'|} \right) \right)^{-1} \in (0, \infty). $$

For $(j, k) \in [1, d] \times [1, d]$, we define the second-order Riesz operator as $\mathcal{R}_{jk} = -\partial_j(-\Delta)^{-1}\partial_k$. The piecewise Hölder estimate in Lemma A.1 for the operator $\mathcal{T} = \mathcal{R}_{jk}$ (respectively, $\mathcal{T} = \mathcal{R}_{jk}\mathcal{R}_{lm}$) involves a polynomial denoted by $\mathfrak{P}^{jk}$ (respectively, $\mathfrak{P}^{jklm}$). Next, we consider a polynomial $\mathfrak{P}$ that satisfies:

$$ \mathfrak{P}^{jk}(a), \mathfrak{P}^{jklm}(a) \le \mathfrak{P}(a), \text{ for all } a > 0. $$

Finally, we introduce the following quantity:

$$ (1.4) \qquad \mathfrak{P}_{\mathcal{C}} := (1 + |\mathcal{C}|)\mathfrak{P}(\|\mathcal{C}\|_{\text{Lip}} + |\mathcal{C}|_*) \|\mathcal{C}\|_{\mathcal{C}^{1+\alpha}}. $$

As noted in [20, Section 3.1] (see also [10, Remark 2, page 35]), there exists a function $\varphi: \mathbb{R}^d \to \mathbb{R} \in \mathscr{C}^{1+\alpha}$ such that

$$ D = \{x \in \mathbb{R}^d : \varphi(x) > 0\}, \quad \text{and} \quad |\nabla \varphi|_{\text{inf}} := \inf_{x \in \partial D} |\nabla \varphi(x)| > 0. $$

We then define:

$$ (1.5) \qquad \ell_{\varphi} = \min \left\{ 1, \left( \frac{|\nabla \varphi|_{\text{inf}}}{\|\nabla \varphi\|_{\mathcal{C}^{\alpha}}} \right)^{1/\alpha} \right\}. $$

(2) Since

$$ \mathbb{R}^d = D \cup \mathcal{C} \cup (\mathbb{R}^d \setminus \bar{D}), $$

we define the space of piecewise $\alpha$-Hölder continuous functions with respect to $\mathcal{C}$ as follows:

$$ \mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d) := \mathcal{C}^{\alpha}(\bar{D}) \cap \mathcal{C}^{\alpha}(\mathbb{R}^d \setminus D), \text{ with } \|g\|_{\mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)} := \|g\|_{\mathcal{C}^{\alpha}(\bar{D})} + \|g\|_{\mathcal{C}^{\alpha}(\mathbb{R}^d \setminus D)} $$

and the non-homogeneous space

$$ \mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d) := L^{\infty}(\mathbb{R}^d) \cap \mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d) \text{ with } \|g\|_{\mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)} := \|g\|_{L^{\infty}(\mathbb{R}^d)} + \|g\|_{\mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)}. $$

This space strictly contains the Hölder space $\mathcal{C}^{\alpha}(\mathbb{R}^d)$.

(3) Given a function $g \in \mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)$, the jump $\llbracket g \rrbracket$, and the average $\langle g \rangle$ of $g$ on $\mathcal{C}$ are defined as follows: for all $\sigma \in \mathcal{C}$,

$$ \llbracket g \rrbracket(\sigma) := \lim_{r \to 0} [g(\sigma + rn_x(\sigma)) - g(\sigma - rn_x(\sigma))] \text{ and } \langle g \rangle(\sigma) := \frac{1}{2} \lim_{r \to 0} [g(\sigma + rn_x(\sigma)) + g(\sigma - rn_x(\sigma))]. $$

Above $n_x$ denotes the outward normal vector of $\mathcal{C}$.

(4) Consider a time-dependent interface $\mathcal{C} = \mathcal{C}(t)$, with local parametrizations $\gamma_j = \gamma_j(t, s)$, $j \in [1, J]$ defined on $I \times V_j \subset \mathbb{R} \times \mathbb{R}^{d-1}$. We assume that $\gamma_j \in \mathcal{C}(I, \mathcal{C}^{1+\alpha}(V_j))$, and for all $t \in I$, $\mathcal{C}(t)$ is a $\mathcal{C}^{1+\alpha}$-regular manifold of dimension $d-1$, that forms the boundary of an open, bounded and simply connected domain $D(t)$. We introduce the following space:

$$ L^p(I, \mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)) := \left\{ g = g(t, x): \begin{cases} \displaystyle \int_I \|g(t)\|_{\mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)}^p dt < \infty & \text{if } 1 \le p < \infty, \\ \displaystyle \sup_{t \in I} \|g(t)\|_{\mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)} < \infty & \text{if } p = \infty \end{cases} \right\}. $$

For $I = (0, T)$, the corresponding space will sometimes be shortly denoted by $L_T^p \mathcal{C}_{\text{pw},\mathcal{C}}^{\alpha}(\mathbb{R}^d)$.
on both sides of a suitable interface $C(0)$ through which it is discontinuous. The initial velocity belongs to the Sobolev space $H^1(\mathbb{R}^d)$ with the following compatibility condition:

$$ (1.8) \qquad \mathrm{div}(\Pi)|_{t=0} \in L^2(\mathbb{R}^d) $$

which expresses somehow the continuity of the stress tensor. This result is more general than the one by Tani [26] and it can be used to serve as building blocks for global solutions.

Unfortunately, we do not have enough information in order to prove existence of global weak solution for the two-fluid model, even in the case of small initial data. However, we are able to prove a global well-posedness result in the spirit of [11] for a simplified situation which is the Navier-Stokes equations with density-dependent viscosity. To keep the paper a reasonable size, we plan to address this issue in the near future.

## 1.2. Review of known results.
The study of the two-phase model (1.1) goes back to Tani's work [26] where the author showed existence of a local-in-time solution for the system (1.1) including temperature. Like the pioneering work of Nash [22], his proof relies on a change in Lagrangian coordinates, followed by the construction of an approximate solution using the kernel associated to the linearized equations. The solution is constructed in a piecewise Hölder spaces: the density is $C^{1+\alpha}$ and the velocity is $C^{2+\alpha}$, on both sides of an interface $C$ which is $C^{2+\alpha}$. To our knowledge, this result is the only one concerning interfaces with Hölder regularity.

In the literature, other results mostly focus on interfaces with Sobolev regularity. As an example, Denisova in [6] has proved the local well-posedness of the system (1.1) including surface tension with different constant viscosity for each component of the mixture and with Sobolev-Slobodetskii interface $C \in W_2^{3/2}$. Her proof also relies on a change in Lagrangian coordinates and the density and the velocity are piecewise $W_2^{1+l}$ (for $l \in (1/2, 1)$) regular on both sides of the interface $C$.

In contrast to the aforementioned results based on a change in Lagrangian coordinates, there are other works [14, 16, 15] by Jang et al. that rely on a change into flattened coordinates. They considered the two-phase problem in a cylinder, taking into account scenarios with and without surface tension. This problem involves two fluids, with one fluid on top of another. In this context, there are two interfaces to consider: the interface between the lower and upper fluids, and the interface between the upper fluid and the atmosphere. Both interfaces are assumed to be represented by graphs. Additionally, the viscosity of the two fluids are assumed to be constant and different from each other. In [14], they successfully constructed local-in-time solution in Sobolev framework, in [15], they studied the Rayleigh-Taylor instability which occurs when the heavier fluid is on top of a lighter fluid. They successfully proved that when the surface tension is below a critical value, then the system is nonlinearly unstable. Finally in [16], they proved that if the heavier fluid is on top of the lighter one and the surface tension is sufficiently large, then, the Rayleigh-Taylor instability is prevented and in the converse configuration, the system is stable around the equilibrium state with non-negative surface tension. Moreover, they studied the zero surface tension limit problem.

On the other hand, there are also results [25, 29] that address the two-phase problem involving a compressible fluid and an incompressible fluid. In [25], Shibata has proved the local well-posedness of this two-phase model, including temperature and surface tension, with an almost flat interface. He worked on an exterior domain and with a viscosity that depends on both density and temperature. The initial density of the compressible fluid has a Sobolev $W^{1,p}$ regularity and the initial velocity and temperature have the Besov regularity $B_{q,p}^{2(1-1/p)}$, on both sides of the interface, for some $1 < p, q < \infty$ that verify:

$$ \frac{2}{p} + \frac{d}{q} < 1. $$

The same Compressible-Incompressible model is considered in [29] by Yan and Zhao with surface tension in horizontal periodic domain. The viscosity of the fluids involved are different from each other and assumed to be constant. They successfully achieved a global well-posedness result by making the assumption that the initial density of the compressible fluid is close to a constant state in $H^4$ with respect to the domain occupied by the compressible fluid, and the initial velocity is small in $H^4$ on both sides of the interface. To our knowledge, this result stands as the only one to address the existence of strong global solution of the two-phase problem.

Finally, we may notice that Kubo et al. in [18] have addressed the local well-posedness of the Compressible-Compressible model with different approach. Their methods rely on proving the existence of a $\mathcal{R}$-bounded operator and maximal regularity of the linearised system. The interface is uniformly $W_r^{2-1/r}$ for $r > d$, the density has a $W^{1,r}$ regularity on both sides of the interface and the initial velocity belongs to the Besov space $B_{r,p}^{2(1-1/p)}$ for some $p > 2$ on both sides of the initial interface. More recently, in [19], Kuo and Shibata have constructed local-in-time solution of the one phase model in half-space with non-slip boundary conditions. As the result above, the solution is also formulated, in some maximal regularity Besov space framework.
1.1.2. Reformulation of the system. A suitable way to rewrite the system (1.1) is to introduce the so-called fractional density or concentration c:

$$c := \frac{\rho_+}{\rho}.$$

From the first two equations of (1.1), one checks easily, assuming some regularity, that the concentration solves the transport equation

$$(1.6) \qquad \partial_t c + u \cdot \nabla c = 0.$$

Let us notice that the domain occupied by the phase $\pm$ can be defined by:

$$D_{\pm}(t) := \{x \in \mathbb{R}^d : \rho_{\pm}(t, x) > 0\}.$$

As well, one can write the domain occupied by the one phase fluid with the help of the concentration as follows:

$$(1.7) \qquad D_{+}(t) = \{x \in \mathbb{R}^d : c(t, x) = 1\}, \quad D_{-}(t) = \{x \in \mathbb{R}^d : c(t, x) = 0\}$$

and consequently the interface is:

$$C(t) := \partial D_{+}(t) = \partial D_{-}(t).$$

Assume that the velocity is regular, for instance Lipschitz, then it admits a continuous flow map $X(t)$ that helps together with (1.6) and (1.7) to obtain that the initial interface is just transported by the fluid

$$C(t) = X(t)C(0).$$

On the other hand, summing up the mass equation on the density of the one phase $\rho_{\pm}$ (1.1)$_1$, one obtains that the total density solves a mass equation with velocity $u$. Besides, the density $\rho_{\pm}$ of the one phase fluid can be written in terms of the total density and the concentration as follows:

$$\rho_{+} = c\rho \quad \text{and} \quad \rho_{-} = (1-c)\rho,$$

thus, the pressure can be written as function of the total density and the concentration, $P = P(\rho, c)$, such that $P_+(\rho) = P(\rho, 1)$ and $P_-(\rho) = P(\rho, 0)$ as well as the viscosity. Hence, formally $(c, \rho, u)$ solves:

$$(1.8) \qquad \begin{cases} \partial_t c + u \cdot \nabla c = 0, \\ \partial_t \rho + \operatorname{div}(\rho u) = 0, \\ \partial_t (\rho u) + \operatorname{div}(\rho u \otimes u) + \nabla P(\rho, c) = \operatorname{div}(2\mu(\rho, c)\mathbb{D}u) + \nabla(\lambda(\rho, c) \operatorname{div} u). \end{cases}$$

Conversely, if $(c, \rho, u)$ solves (1.8), one can set

$$\rho_{+} = c\rho \quad \text{and} \quad \rho_{-} = (1-c)\rho$$

and checks easily that $\rho_{\pm}$ solve the equations (1.1)$_1$. Just by replacing

$$\rho = \rho_{+} + \rho_{-} \quad \text{and} \quad c = \frac{\rho_{+}}{\rho_{+} + \rho_{-}}$$

in the pressure and the viscosity, one infers that $(\rho_+, \rho_-, u)$ solves (1.1). In fact, we obtain that $(\rho_+, \rho_-, u)$ solves (1.1) if and only if $(c, \rho, u)$ solves (1.8) owing to some change of pressure and viscosity law, at least for classical solutions. The equivalence also holds true for weak solutions à la Lions (we mean weak solutions which satisfy the classical energy inequality associated with (1.1), see [34]) due to the result in [34, Lemma 2.5], which is part of the DiPerna-Lions theory. Therefore, we will only work with the system (1.8).

The equations (1.8) are supplemented with initial data

$$(1.9) \qquad c|_{t=0} = c_0, \quad \rho|_{t=0} = \rho_0 \quad \text{and} \quad u|_{t=0} = u_0.$$

Given the equilibrium state of the pressure $\tilde{P}$ introduced in (1.2), the equilibrium state of the density, $\tilde{\rho}$, is defined by the following relation:

$$(1.10) \qquad \tilde{P} = P(\tilde{\rho}, 1).$$

Next, we define the equilibrium states of the viscosities:

$$\tilde{\mu} := \mu(\tilde{\rho}, 1) \quad \text{and} \quad \tilde{\lambda} := \lambda(\tilde{\rho}, 1).$$

We recall that $\mu$, $\lambda$ and $P$ are $\mathcal{C}^2$-regular, with $P$ and $\lambda$ being non-negative, while $\mu(\rho, c) > 0$ for $\rho > 0$.

The purpose of this paper is to construct a unique local-in-time solution $(c, \rho, u)$ of the equations (1.8) with discontinuous initial data and with density-dependent viscosity: the initial density $\rho_0$ is assumed to be Hölder
As explained, for instance in [9], a different approach to tracking the interface is to construct weak solutions for the full model in a class that allows the study of its dynamics. To the best of our knowledge, the only available results for weak solutions pertain to fluid mixture with constant viscosity that are the same for both components. For instance, we refer to the work by Vasseur, Wen, and Yu [27] in a three-dimensional bounded domain where they obtained a weak solution with finite initial energy. The main difficulty in their proof lies in the compactness of an approximate sequence of the pressure that depends on two variables, in contrast to the single-variable case [21, 7]. They successfully overcome the difficulty by considering the following pressure law

$$P(\rho_+, \rho_-) = A_+ \rho_+^{\gamma_+} + A_- \rho_-^{\gamma_-}$$

with the following closeness constraints on $\gamma_\pm$:

$$\gamma_{\pm} > 9/5 \quad \text{and} \quad \max \left\{ \frac{3\gamma_{+}}{4}; \gamma_{+} - 1; \frac{3(\gamma_{+} + 1)}{5} \right\} < \gamma_{-} < \min \left\{ \frac{4\gamma_{+}}{3}; \gamma_{+} + 1; \frac{5\gamma_{+}}{3} - 1 \right\}.$$

This result has been improved by Novotný [23], Wen [28] by allowing more general pressure law, Novotný and Pokortný [24] and Novotný et al. [17] for multi-components fluid models. It is important to notice that, although the pressure law of the components of the mixture are different, the viscosity of the components are supposed to be the same constant. The solution that they constructed is too weak in order to track down discontinuity in the density. In one dimension, Bresch, Burtea and Lagoutière in [1] successfully solved the global well-posedness for non common constant viscosity, with numerical considerations in [2]. In this paper, we consider density- dependent viscosity for both components of the mixture, and we achieve the construction of unique weak solution in $\mathbb{R}^d$, ($d \in \{2, 3\}$) that allows the study of the dynamics of the interface.

**1.3. Main result.** We consider the Cauchy problem associated with (1.5) with initial data (1.6). We assume that the initial concentration $c_0$ is the characteristic function of $D(0)^c$ where $D(0)$ is an open, bounded, simply connected domain in $\mathbb{R}^d$ such that its boundary

$$\partial D(0) := C(0) = \{\gamma_0(s) \in \mathbb{R}^d, s \in \mathbb{T}^{d-1}\}$$

with $\gamma_0 \in \mathcal{C}^{1+\alpha}$ verifying the non-self intersection condition:

$$(1.9) \qquad \exists c_{\gamma}(0) > 0, \quad \forall (s, s') \in (\mathbb{T}^{d-1})^2, \quad |\gamma_0(s) - \gamma_0(s')| \ge c_{\gamma}(0)^{-1}|s - s'|.$$

Also the initial density $\rho_0$ and velocity $u_0$ have the following regularity:

$$(1.10) \qquad u_0 \in H^1(\mathbb{R}^d), \quad P(\rho_0, c_0) - \tilde{P} \in L^2(\mathbb{R}^d) \cap \mathcal{C}_{\text{pw}, \gamma_0}^{\alpha}(\mathbb{R}^d).$$

Finally, we assume that the initial data is such that initially the material derivative that is

$$\dot{u} := \partial_t u + (u \cdot \nabla)u$$

belongs to $L^2(\mathbb{R}^d)$, in other words,

$$(1.11) \qquad (\rho \dot{u})|_{t=0} = \text{div}(\Pi)|_{t=0} \in L^2(\mathbb{R}^d).$$

A very technical step requires us to obtain estimate for the $L^{8/d}((0, T), L^4(\mathbb{R}^d))$ norm of $\nabla \dot{u}$ and in order to do so, we need to look at the second material derivative of the velocity $\ddot{u}$, which is

$$\ddot{u} := \partial_t \dot{u} + (u \cdot \nabla)\dot{u}.$$

This variable appears while applying the material derivative $\partial_t \cdot + \text{div}(\cdot u)$ to the momentum equation (1.5)$_3$. Our main result reads as follows.

**Theorem 1.2.** Let $(c_0, \rho_0, u_0)$ be initial data associated with the system (1.5) that satisfies the conditions (1.9), (1.10), (1.11) for some $0 < \alpha < 1$ and

$$0 < \rho_* \le \rho_0 \le \rho^* \quad \text{and} \quad 0 < \mu_* \le \mu(\rho_0, c_0), \lambda(\rho_0, c_0) \le \mu^*.$$

There exists $[\mu] > 0$ depending on the dimension $d \in \{2, 3\}$ and the initial interface $C(0)$ such that if

$$(1.12) \qquad \frac{4}{\mu} \|\mu(\rho_0, c_0) - \tilde{\mu}\|_{\mathcal{C}_{\text{pw}, \gamma_0}^{\alpha}(\mathbb{R}^d)} \le [\mu],$$

then, there exists a time $T > 0$ and an unique solution $(c, \rho, u)$ of the system (1.5) verifying the following:

$$(1) \quad \forall 0 \le t \le T, \quad c(t) = 1_{D(t)^c} \quad \text{with} \quad \partial D(t) := C(t) = \{\gamma(t, s) \in \mathbb{R}^d, s \in \mathbb{T}^{d-1}\} \quad \text{and} \quad \gamma(t) \in \mathcal{C}^{1+\alpha};$$

$$(2) \quad P(\rho, c) - \tilde{P} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d) \cap \mathcal{C}_{\text{pw}, \gamma}^{\alpha}(\mathbb{R}^d));$$

$$(3) \quad u \in \mathcal{C}([0, T], H^1(\mathbb{R}^d)) \cap L^\infty((0, T), \dot{W}^{1,6}(\mathbb{R}^d)) \cap L^{16}((0, T), \dot{W}^{1,8}(\mathbb{R}^d)), \sigma^{r/4} \nabla u \in L^4((0, T), \mathcal{C}_{\text{pw}, \gamma}^{\alpha}(\mathbb{R}^d)) \text{ for } 0 < r < 2;$$
continuous on both sides of a suitable interface $C_0$ through which it is discontinuous. The initial velocity belongs to the Sobolev space $H^1(\mathbb{R}^d)$ with the following compatibility condition:

$$ (1.11) \qquad \mathrm{div}(\Pi)|_{t=0} \in L^2(\mathbb{R}^d). $$

This condition expresses the continuity of the normal component of the stress tensor and does not require smoothness of the density. The parabolic effect of the momentum equations ensures that (1.11) holds true at positive times even for less regular initial data (see Hoff [13]). In particular it does not prevent to extend solutions until blow-up.

Our result is more general than the one by Tani [32] (see Section 1.2 below for more details), as it allows for less regular initial data, interface, and it can serve as building blocks for global solutions.

We do not have enough information in order to prove existence of global weak solution for the two-fluid model, even in the case of small initial data. However, we are able to prove a global well-posedness result in the spirit of [14] for a simplified situation which is the Navier-Stokes equations with density-dependent viscosity. To keep the paper a reasonable size, we plan to address this issue in the near future.

## 1.2. Review of known results.
The study of the two-phase model (1.1) goes back to Tani's work [32] where the author showed existence of a local-in-time solution for the heat conductive variant of (1.1). As in the pioneering work of Nash [26], Tani's proof relies on a change of variables into Lagrangian coordinates, followed by the construction of an approximate solution using the kernel associated to the linearized equations. The solution is constructed in piecewise Hölder spaces: the density is $C^{1+\alpha}$ and the velocity is $C^{2+\alpha}$, on both sides of an interface $C$ which is $C^{2+\alpha}$. To the best of our knowledge, Tani's result is the only one addressing the two-fluid model with Hölder regularity for interfaces.

In the literature, other results mostly focus on interfaces with Sobolev regularity. As an example, Denisova in [8] has proved the local well-posedness of the system (1.1) including surface tension with different constant viscosity for each component of the mixture and with Sobolev-Slobodetskii interface $C \in W^{\frac{5}{2}+l,2}$, with $l \in (1/2, 1)$. Her proof also relies on a change in Lagrangian coordinates and the density and the velocity are piecewise $W^{1+l,2}$ regular on both sides of the interface $C$.

In contrast to the aforementioned results based on a change in Lagrangian coordinates, there are other works [17–19] by Jang et al. that rely on a change into flattened coordinates. The authors considered the two-phase problem in a cylinder, taking into account scenarios with and without surface tension. This problem involves two fluids, with one fluid on top of another. In this context, there are two interfaces to consider: the interface between the lower and upper fluids, and the interface between the upper fluid and the atmosphere. Both interfaces are assumed to be represented by graphs. Additionally, the viscosity of the two fluids are assumed to be constant and different from each other. In [17], they successfully constructed local-in-time solution in Sobolev framework, in [18], they studied the Rayleigh-Taylor instability which occurs when the heavier fluid is on top of a lighter fluid. They successfully proved that when the surface tension is below a critical value, then the system is nonlinearly unstable. Finally in [19], they proved that if the heavier fluid is on top of the lighter one and the surface tension is sufficiently large, then, the Rayleigh-Taylor instability is prevented and in the converse configuration, the system is stable around the equilibrium state with non-negative surface tension. Moreover, they studied the zero surface tension limit problem.

On the other hand, there are also results [31, 36] that address the two-phase problem involving a compressible fluid and an incompressible fluid. In [31], Shibata has proved the local well-posedness of this two-phase model, including temperature and surface tension, with an almost flat interface. The author worked on an exterior domain and with a viscosity that depends on both density and temperature. The initial density of the compressible fluid has a Sobolev $W^{1,p}$ regularity and the initial velocity and temperature have the Besov regularity $B_{q,p}^{2(1-1/p)}$, on both sides of the interface, for some $1 < p, q < \infty$ that verify:

$$ \frac{2}{p} + \frac{d}{q} < 1. $$

The same Compressible-Incompressible model is considered in [36] by Yan and Zhao with surface tension in horizontal periodic domain. The viscosity of the fluids involved are different from each other and assumed to be constant. They successfully achieved a global well-posedness result by making the assumption that the initial density of the compressible fluid is close to a constant state in $H^4$ with respect to the domain occupied by the compressible fluid, and the initial velocity is small in $H^4$ on both sides of the interface. To our knowledge, this result stands as the only one to address the existence of strong global solution of the two-phase problem.

Finally, we may notice that Kubo et al. in [22] have addressed the local well-posedness of the Compressible-Compressible model with different approach. Their methods rely on proving the existence of a $\mathcal{R}$-bounded solution
$$ (4) \quad \dot{u} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)), \quad \sqrt{\sigma} \nabla \dot{u} \in L^\infty((0, T), L^2(\mathbb{R}^d)) \cap L^{8/d}((0, T), L^4(\mathbb{R}^d)); $$

$$ (5) \quad \sqrt{\sigma} \dot{u} \in L^2((0, T) \times \mathbb{R}^d), \quad \sigma \ddot{u} \in L^\infty((0, T), L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)), $$

$$ \text{where } \sigma(t) = \min\{1, t\}. $$

**Remark 1.3.** Because of the roughness of the viscosity, we use a perturbation argument to obtain $L^p$ and piecewise Hölder estimates for the velocity gradient. In fact, as we will see in (3.19), we have the following expression for the velocity gradient:

$$ (1.13) \qquad \begin{aligned} \nabla u = \frac{1}{\mu}(-\Delta)^{-1}\nabla [\mathcal{P} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} - \mathcal{P}(\rho_0 \partial_t u) + \mathcal{P} \operatorname{div} F] \\ + \frac{1}{2\tilde{\mu} + \lambda}(-\Delta)^{-1}\nabla [\mathcal{Q} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} - \mathcal{Q}(\rho_0 \partial_t u) + \mathcal{Q}(\operatorname{div} F)] \end{aligned} $$

where $\mathcal{P}$ is the Leray projector onto the space of divergence-free vector fields and where $\mathcal{Q}$ is the projector onto the space of curl-free vector fields. Thus, we need the smallness condition (1.12) on the dynamic viscosity in order to ensure that the first and third terms on the right-hand side of (1.13) are small in comparison to the left-hand side.

**Remark 1.4.** The solution $(c, \rho, u)$ constructed in the above theorem exhibits the following properties.

(1) There exists a constant $R$ depending on the time $T$, on the norms of the initial data $(c_0, \rho_0, u_0)$ and on the initial interface $\mathcal{C}(0)$ such that:

$$ \sup_{[0,T]} \{\|u\|_{H^1(\mathbb{R}^d)} + \|\dot{u}\|_{L^2(\mathbb{R}^d)}\} \le R. $$

(2) From (2) and (4), one has $\rho \dot{u} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d))$ and one deduces from the momentum equation (1.5) that $\operatorname{div}(\Pi) \in \mathcal{C}([0, T], L^2(\mathbb{R}^d))$ so that the compatibility condition (1.11) is propagated in time.

(3) From (3), one deduces by interpolation that the velocity $u$ is continuous in time and space. On the other hand, from (4), its material derivative, $\dot{u}$ is continuous in space for almost every positive time. Adding the fact that its material derivative $\ddot{u}$ verifies (5), then $\ddot{u}$ is at least continuous across the interface $\mathcal{C}(t)$. Roughly, it follows, that for all $s \in \mathbb{T}^{d-1}$, one has:

$$ [[\dot{u}(t, \gamma(t, s))]] = 0 \implies \frac{d}{dt}[[\dot{u}(t, \gamma(t, s))]] = [[\ddot{u}(t, \gamma(t, s))]] = 0. $$

(4) From the last part of (4), we deduce that, $\nabla u \in L^1((0, T), \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d))$ and, of course, this is the key that helps the propagation the regularity of the initial interface $\mathcal{C}(0)$. In particular $\mathcal{C}(t)$ belongs to $\mathcal{C}^{1+\alpha}$ and verifies the non-self intersection condition (1.9) with constant $c_{\gamma}(t)$ verifying:

$$ (1.14) \qquad \begin{cases} c_{\gamma}(t) \le c_{\gamma}(0) \exp \left( \int_{0}^{t} \|\nabla u(\tau)\|_{L^{\infty}(\mathbb{R}^d)} d\tau \right), \\ \|\nabla \gamma(t)\|_{\mathcal{C}^{\alpha}} \le \left( \|\nabla \gamma_0\|_{\mathcal{C}^{\alpha}} + \|\nabla \gamma_0\|_{L^{\infty}}^{1+\alpha} \int_{0}^{t} \|\nabla u(\tau)\|_{\mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)} d\tau \right) \exp \left( (2+\alpha) \int_{0}^{t} \|\nabla u(\tau)\|_{L^{\infty}(\mathbb{R}^2)} d\tau \right). \end{cases} $$

(5) The dynamic viscosity $\mu(\rho, c)$ remains close to $\tilde{\mu}$:

$$ \frac{2}{\tilde{\mu}} \sup_{[0,T]} \|\mu(\rho, c) - \tilde{\mu}\|_{\mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)} \le [\mu]. $$

(6) We notice that:

$$ (1.15) \qquad \rho \ddot{u} = \operatorname{div}(2\mu(\rho, c)\mathbb{D}\dot{u}) + \nabla(\lambda(\rho, c) \operatorname{div} \dot{u}) + \operatorname{div} \mathcal{F}, $$

where $\mathcal{F}$ is a lower order term that belongs, at least, to $L^2((\tau, T), \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d))$ for all $0 < \tau < T$. This equation is obtained by applying the material derivative $\partial_t \cdot + \operatorname{div}(\cdot u)$ to the momentum equation (1.5)3. Using the regularity of $\ddot{u}$ in (5), the smallness of the dynamic viscosity around $\tilde{\mu}$, (1.15), the regularity of $\mathcal{F}$ and the quantitative Hölder estimate [8, 9], we have that:

$$ \nabla \dot{u} \in L^2((\tau, T), \mathcal{C}_{pw,\gamma}^{\beta}(\mathbb{R}^d)) \quad \text{with} \quad \beta = \begin{cases} \alpha, & \text{in } d=2, \\ \min\{\alpha, 1/2\}, & \text{in } d=3. \end{cases} $$

In particular, $\nabla \dot{u}$ is also continuous on both sides of the interface $\gamma$.
operator and maximal regularity of the linearised system. The interface is uniformly $W^{r,2^{-1/r}}$ for $r > d$, the density has a $W^{1,r}$ regularity on both sides of the interface and the initial velocity belongs to the Besov space $B_{r,p}^{2(1-1/p)}$ for some $p > 2$ on both sides of the initial interface. More recently, in [23], Kuo and Shibata have constructed local-in-time solution of the one phase model in half-space with non-slip boundary conditions. As the result above, the solution is also formulated, in some maximal regularity Besov space framework.

As explained, for instance in [12], a different approach to tracking the interface is to construct weak solutions for the full model in a class that allows the study of its dynamic. To the best of our knowledge, the only available results for weak solutions pertain to fluid mixture with constant viscosity that are the same for both components. For instance, we refer to the work by Vasseur, Wen, and Yu [34] in a three-dimensional bounded domain where they obtained a weak solution with finite initial energy. The main difficulty in their proof lies in the compactness of an approximate sequence of the pressure that depends on two variables, in contrast to the single-variable case [9, 24]. They successfully overcome the difficulty by considering the following pressure law

$$P(\rho_+, \rho_-) = A_+ \rho_+^{\gamma_+} + A_- \rho_-^{\gamma_-}$$

with the following closeness constraints on $\gamma_\pm$:

$$\gamma_\pm > 9/5 \quad \text{and} \quad \max \left\{ \frac{3\gamma_+}{4}; \gamma_+ - 1; \frac{3(\gamma_+ + 1)}{5} \right\} < \gamma_- < \min \left\{ \frac{4\gamma_+}{3}; \gamma_+ + 1; \frac{5\gamma_+}{3} - 1 \right\}.$$

This result has been improved by Novotný [27], Wen [35] by allowing more general pressure law, Novotný and Pokortný [28] and Novotný et al. [21] for multi-components fluid models. It is important to notice that, although the pressure law of the components of the mixture are different, the viscosity of the components are supposed to be the same constant. The solutions that the authors constructed are too weak in order to track down discontinuities in the density. In one dimension, Bresch, Burtea and Lagoutière in [2] successfully solved the global well-posedness for non common constant viscosity, with numerical considerations in [3].

A class of weak solutions with more regularity (for the one phase model) has been proposed in the pioneering work of Hoff [13], see also [29, 33] for bulk viscosity coefficient $\lambda$ depending on the density and Perepelitsa [30] for the case of bounded domains. However, these papers do not address the density patch problem. This issue is explored in Hoff [14] and Hoff and Santos [15], which examine the dynamics of sharp interfaces for the single-phase model, namely the Navier-Stokes equations with constant viscosity. Persistence of $C^{1+\alpha}$-regular interface is achieved in [14] for $d=2$, under the condition that the density is piecewise $\alpha$-Hölder continuous. However, the analysis is restricted to a linear pressure law and small viscosity $\lambda$. In [15] the density is only bounded, the gradient of the velocity is only BMO in space and for interfaces that are initially $C^\alpha$ one can only ensure $C^{\alpha(t)}$-regularity for $t > 0$ where $\alpha(t)$ is decaying exponentially to zero. Although not stated, a similar result i.e. propagation of Hölder regularity with exponential loss also holds in the context of [29, 33] in which the viscosity $\lambda$ is allowed to depend on the density, the analysis being very similar to [15]. However, when the viscosity $\mu$ depends on the density, there is no clear notion of effective flux, the analysis complicates and it is not even clear how can one propagate the $L^\infty$-norm of the density, we refer to an interesting work by Bresch and Jabin [4] for a discussion regarding the complications arising in this case.

In this paper, we consider density-dependent viscosity for both components of the mixture, and we achieve the construction of unique weak solution in $\mathbb{R}^d$, ($d \in \{2, 3\}$) for the system (1.8) which allows the study of the dynamics of the interface.

### 1.3. Main result.
We consider the Cauchy problem associated with equations (1.8) and initial data (1.9). We assume the existence of an open, bounded, and simply connected domain $D_0$ such that:

$$(1.12) \qquad c_0 = \mathbb{1}_{D_0}.$$

Also, we assume that the boundary $C_0 := \partial D_0$ fulfills the conditions outlined in Section 1.1.1 item 1. Additionally, the initial density $\rho_0$ and velocity $u_0$ satisfy:

$$(1.13) \qquad u_0 \in H^1(\mathbb{R}^d), \quad P(\rho_0, c_0) - \tilde{P} \in L^2(\mathbb{R}^d) \cap C_{pw,C_0}^{\alpha}(\mathbb{R}^d), \quad \text{and} \quad \text{div}(\Pi)|_{t=0} \in L^2(\mathbb{R}^d).$$

A technical step requires us to obtain estimate for the $L^{8/d}((0, T), L^4(\mathbb{R}^d))$ norm of $\sigma^{\frac{d}{4}}\nabla \dot{u}$, where

$$\sigma(t) = \min\{1, t\} \quad \text{and} \quad \dot{u} := \partial_t u + (u \cdot \nabla)u$$

is the material acceleration of fluid particles. To achieve this, we need to look at the second material derivative of the velocity, which is

$$\ddot{u} := \partial_t \dot{u} + (u \cdot \nabla) \dot{u}.$$
Let us notice that Theorem 1.2 comes with the following blow-up criterion:

**Theorem 1.5** (Blow-up criterion). Let $(c, \rho, u)$ be the solution constructed in Theorem 1.2 defined up to a maximal time $T^*$. Assume that, for all $0 < t < T^*$ we have:

$$ \limsup_{t \to T^*} \left\{ c_{\gamma}(t) + \|\nabla \gamma(t)\|_{\mathscr{C}^{\alpha}} + \left\| \frac{1}{\rho(t)}, \frac{1}{\mu(\rho(t), c(t))}, \frac{1}{\lambda(\rho(t), c(t))} \right\|_{L^{\infty}(\mathbb{R}^d)} \right\} \\ + \limsup_{t \to T^*} \left\{ \|u(t)\|_{H^1(\mathbb{R}^d)} + \|(\rho \dot{u})(t)\|_{L^2(\mathbb{R}^d)} + \|P(\rho(t), c(t)) - \tilde{P}\|_{\mathscr{C}^{\alpha}_{pw, \gamma(t)}(\mathbb{R}^d)} \right\} < \infty $$

and

$$ \frac{1}{\tilde{\mu}} \sup_{[0,t]} \|\mu(\rho, c) - \tilde{\mu}\|_{\mathscr{C}^{\alpha}_{pw, \gamma}(\mathbb{R}^d)} \leq [\mu] $$

then $T^* = \infty$.

**Remark 1.6.** Our main result can be seen as a first step towards the existence of global solutions with small data which we plan to attack in a forthcoming paper. Indeed, as explained in [11], the construction of discontinuous weak-solutions for compressible models cannot be achieved in the classical way by regularizing the initial data and applying some classical results à la Matsmura-Nishida since this procedure destroys the discontinuities in the density. Some difficult constructions are then needed in order to construct discontinuous weak-solutions, see [11, Section 3]. We believe that since we dispose of powerful functional analysis tools, see Appendix A, an approach where local-well posedness is first proved to be quite natural. With Theorem 1.2 in hand, in order to attack the global existence issue, we do not have to worry anymore about constructing appropriate approximate systems that preserve the main feature of the system (1.5).

**Remark 1.7.** Unfortunately, we do not have enough information in order to prove global existence for the two-fluid model, even in the case of small data. The difficulty lies in the fact that we are not able to estimate the $L^1((0, T), \mathscr{C}^{\alpha}_{pw, \gamma}(\mathbb{R}^d))$ norm of $\nabla u$ uniformly with respect to time. In particular, see (1.14), the characteristic $c_{\gamma}(t)$ and the $\mathscr{C}^{1+\alpha}$ norm of the interface $\gamma(t)$ are growing exponentially with respect to time.

In the study of the Navier-Stokes equations in $d=2$, Hoff [11], Hoff and Santos [12] have observed that the exponential growth of these quantities can be compensated by the exponential decay of the pressure jump across the interface. For the two-fluid model, when the pressure laws of the two phases are different, this exponential decay of the pressure jump does not hold, even when both phases have the same constant viscosity. The minimal requirements seem to be that:

$$ 0 < v_* \le \frac{P(\rho, 1) - P_-(\rho', 0)}{f(\rho, 1) - f(\rho', 0)} \le v^* \quad \text{and} \quad \frac{|\mu(\rho, 1) - \mu(\rho', 0)|}{|f(\rho, 1) - f(\rho', 0)|} + \frac{|\lambda(\rho, 1) - \lambda(\rho', 1)|}{|f(\rho, 1) - f(\rho', 0)|} \le v^* $$

at least, for all $\rho$ close to $\tilde{\rho}(1)$ and $\rho'$ close to $\tilde{\rho}(0)$ where $\tilde{P} = P(\tilde{\rho}(1), 1) = P(\tilde{\rho}(0), 0)$ and where $f$ is defined, up to a constant, by:

$$ f(\rho, c) = \int^{\rho} \frac{2\mu(s, c) + \lambda(s, c)}{s} ds. $$

The only case in which this seems to occur is when the pressure and viscosity laws are proportional. However, we are able to prove a global well-posedness result in the spirit of [11] for a simplified situation, forthcoming very soon.

*Outline of the paper.* The rest of this paper is divided into three sections. In the first one, Section 2, we rewrite system (1.5) in Lagrangian coordinates, and in the second section, Section 3 we study the linearized system associated to (1.5). Finally the study of the full nonlinear system is presented in Section 4.

## 2. CHANGE OF VARIABLES IN LAGRANGIAN COORDINATES

The aim of this section is to write the system (1.5) in Lagrangian coordinates. We know that for a Lipschitz vector field $u$, one can define a flow map $\mathcal{X}$ given by the Duhamel formula:

$$ \mathcal{X}(t, y) = y + \int_{0}^{t} u(s, \mathcal{X}(s, y))ds $$

and adopting the following notation: for all function $g = g(t, x)$ we denote by $\bar{g}$ the function:

$$ \bar{g}(t, y) := g(t, \mathcal{X}(t, y)) $$
This variable appears while applying the material derivative $\partial_t \cdot + \text{div}(\cdot u)$ to the momentum equation (1.8)$_3$. Our main result reads as follows:

**Theorem 1.1.** Let $(c_0, \rho_0, u_0)$ be initial data associated with the system (1.8), satisfying the conditions (1.10), (1.12) and (1.13). Additionally, assume:

$$ (1.14) \qquad \rho_{*,0} := \inf_{x \in \mathbb{R}^d} \rho_0(x) > 0 \quad \text{and} \quad \mu_{*,0} := \inf_{x \in \mathbb{R}^d} \mu(\rho_0(x), c_0(x)) > 0. $$

There exists $[\mu] > 0$ depending only on the dimension $d \in \{2, 3\}$, $\alpha$ and $\tilde{\mu}$ such that if

$$ (1.15) \quad \left[ 1 + \|\lambda(\rho_0, c_0)\|_{\mathcal{P}_{pw,c_0}(\mathbb{R}^d)} + (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \|\|\lambda(\rho_0, c_0)\|\|_{L^\infty(C_0)} \|\mu(\rho_0, c_0) - \tilde{\mu}\|_{\mathcal{C}_{pw,c_0}^\alpha(\mathbb{R}^d)} \right. \\ \left. \qquad + (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \left[ \|\|\mu(\rho_0, c_0)\|\|\|_{L^\infty(C_0)} + \|\|\mu(\rho_0, c_0)\|\|\|_{L^\infty(C_0)} \left\| 1 - \frac{\tilde{\mu}}{\langle \mu(\rho_0, c_0) \rangle} \right\|_{L^\infty(C_0)} \right] \le \frac{[\mu]}{4}, $$

then, there exist a time $T > 0$ and a unique solution $(c, \rho, u)$ of the system (1.8) verifying the following:

$$ (1) \quad \forall 0 < t \le T, \quad c(t) = 1_{D(t)} \quad \text{with } \partial D(t) := C(t) \quad \text{and with local parametrization } \gamma(t) \in \mathcal{C}^{1+\alpha}; $$

$$ (2) \quad P(\rho, c) - \tilde{P} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d) \cap \mathcal{C}_{pw,C}^\alpha(\mathbb{R}^d)); $$

$$ (3) \quad u \in \mathcal{C}([0, T], H^1(\mathbb{R}^d)) \cap L^\infty((0, T), \dot{W}^{1,6}(\mathbb{R}^d)) \cap L^{16}((0, T), \dot{W}^{1,8}(\mathbb{R}^d)), \sigma^{r/4} \nabla u \in L^4((0, T), \mathcal{C}_{pw,C}^\alpha(\mathbb{R}^d)) \text{ for} $$

$$ (1.16) \qquad r = \begin{cases} \max(1/3, 2\alpha) & \text{if } d=2, \\ \max\left(\frac{9}{14}, 2\alpha\right) & \text{if } d=3, \alpha \in (0, 1/2), \\ \frac{2}{3}(5\alpha - 1) & \text{if } d=3, \alpha \in [1/2, 1). \end{cases} $$

$$ (4) \quad \dot{u} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)), \quad \sqrt{\sigma} \nabla \dot{u} \in L^\infty((0, T), L^2(\mathbb{R}^d)), \quad \sigma^{d/4} \nabla \dot{u} \in L^{8/d}((0, T), L^4(\mathbb{R}^d)); $$

$$ (5) \quad \sqrt{\sigma} \ddot{u} \in L^2((0, T) \times \mathbb{R}^d), \quad \sigma \ddot{u} \in L^\infty((0, T), L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)). $$

To prove **Theorem 1.1**, we transform the system (1.8) into Lagrangian coordinates, see (2.7) below, to fix the domain. The existence of a unique solution $(\bar{c}, \bar{\rho}, \bar{u})$ of (2.7) is obtained in **Theorem 2.2** below. The flow of $\bar{u}$ is regular enough to ensure that, back in Eulerian coordinates, $(c, \rho, u)$ possesses all the properties stated in **Theorem 1.1**.

**Remark 1.2.**

(1) The smallness assumption (1.15) arises from the study of the linearized equations associated with the equations resulting from the transformation of (1.8) into Lagrangian coordinates:

$$ \rho_0 v_t - \text{div}(\mu \mathbb{D}v + \lambda \text{div} v I_d) = \text{div} F. $$

To obtain $L^p(\mathbb{R}^d)$ and piecewise Hölder estimates for $\nabla u$, we express (see (3.21)):

$$ (1.17) \quad \begin{aligned} \nabla v = & \frac{1}{\tilde{\mu}}(-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div} F - \nabla^2(-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda}(-\Delta)^{-1} \operatorname{div} \operatorname{div} F \right] \\ & - \frac{1}{\tilde{\mu}}(-\Delta)^{-1} \nabla \mathcal{P} (\rho_0 v_t) + \nabla^2(-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda}(-\Delta)^{-1} \operatorname{div}(\rho_0 v_t) \right] \\ & + \frac{1}{\tilde{\mu}}(-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}v\} - \nabla^2(-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda}(-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}v\} \right]. \end{aligned} $$

where $\mathcal{P}$ is the Leray projector onto the space of divergence-free vector fields. Under the smallness condition (1.15), the last two terms of (1.17) are small in comparison to the left-hand side.

(2) One of the main difficulties is to establish an estimate for the $L^\infty(\mathbb{R}^d)$-norm for $\nabla v$, (see page 17). Except for the third term in the expression (1.17), which is regular in space, all other terms are second-order or fourth-order Riesz transforms of discontinuous functions (see the regularity of $\mu$, $\lambda$ and $F$ in (3.2)-(3.3)). It is well known that Riesz operators are not continuous on $L^\infty(\mathbb{R}^d)$, which makes the estimation of the $L^\infty(\mathbb{R}^d)$-norm of $\nabla u$ more challenging. We employ a quantitative piecewise Hölder estimate for even-order Riesz operators obtained by Gancedo and García-Juárez (see Lemma A.1 or [11, 12]) to derive a piecewise Hölder estimate for $\nabla v$. Roughly speaking, the authors proved that these operators are continuous on $L^q(\mathbb{R}^d) \cap \mathcal{C}_{pw,c_0}^\alpha(\mathbb{R}^d)$, $q \in [1, \infty)$, with a norm that may depend on the regularity of $C_0$ (see Section 1.1.1-item 1).

<footnote n="1">¹We refer to (1.4)-(1.5) for the definition of $\ell_{\varphi_0}$ and $\mathfrak{P}_{c_0}$.</footnote>
then, the flow of the vector field $u$ can be written in term of $\bar{u}$ by:

$$ \mathcal{X}(t, y) = y + \int_{0}^{t} \bar{u}(s, y) ds. $$

And for any $g = g(t, x)$, one has that

$$ (2.1) \qquad \partial_t \bar{g}(t, y) = (\partial_t g + u \cdot \nabla g)(t, \mathcal{X}(t, y)). $$

Let us notice that by chain rule, for any vector field $G$, one has

$$ D\bar{G} = DG \circ \mathcal{X} \cdot D\mathcal{X}, $$

but, since

$$ (2.2) \qquad (D\mathcal{X}(t, y))_{ij} = \delta^{ij} + \int_{0}^{t} \partial_j \bar{u}^i(s, y) ds, $$

it holds that, at least, for small time, the Jacobian matrix $D\mathcal{X}$ of the change of variable is inversible and then

$$ \overline{DG}(t, y) = D\overline{G}(t, y) \cdot (D\mathcal{X}(t, y))^{-1} = \frac{1}{\det D\mathcal{X}(t, y)} D\overline{G}(t, y) \cdot \text{Adj } D\mathcal{X}(t, y). $$

In particular, one has that:

$$ (2.3) \qquad \overline{\partial_j G^i} = \frac{1}{\det D\mathcal{X}} \partial_k \overline{G}^i (\text{Adj } D\mathcal{X})_{kj}. $$

On the other hand, from the Piola formula: for any matrix $A \in \mathcal{M}_{n \times n}(\mathbb{R})$,

$$ \partial_k (\text{Adj}(A))_{kj} = 0 $$

one deduces that:

$$ (2.4) \qquad \overline{\partial_j G^i} = \frac{1}{\det D\mathcal{X}} \partial_k \left[ \overline{G}^i (\text{Adj } D\mathcal{X})_{kj} \right]. $$

The first form (2.3) will be referred as the non conservative form and the second (2.4) is the conservative one and both will be used. From the non conservative form, it turns out that

$$ \overline{\text{div } G} = \frac{1}{\det D\mathcal{X}} \text{Tr}(D\overline{G} \cdot \text{Adj } D\mathcal{X}), $$

and by setting

$$ A(t, y) = (D\mathcal{X}(t, y))^{-1} $$

one has that

$$ (2.5) \qquad \text{div}_A G := \overline{\text{div } G} = D\overline{G} : A = A^T : \nabla\overline{G}. $$

On the other hand, for all matrix $B \in \mathcal{C}^1([0, T], \mathcal{M}_{n \times n}(\mathbb{R}))$, one has:

$$ \begin{align*} \frac{d \det B}{ds}(s) &= \frac{\partial \det B}{\partial b_{ij}} \frac{\partial b_{ij}}{\partial s}(s) \\ &= (-1)^{i+j} M_{ij}(B(s)) \frac{\partial b_{ij}}{\partial s}(s) \\ &= \text{Adj}(B(s))_{ji} \left( \frac{\partial B}{\partial s}(s) \right)_{ij} \\ &= \frac{\partial B(s)}{\partial s} : \text{Adj}(B(s)). \end{align*} $$

In particular, for $B = D\mathcal{X}$, and remembering from (2.2) that

$$ \frac{\partial D\mathcal{X}}{\partial t} = D\bar{u} $$

one deduces that

$$ \frac{\partial \det D\mathcal{X}}{\partial t}(t, y) = D\bar{u}(t, y) : \text{Adj } D\mathcal{X}(t, y). $$

Thus, from the identity (2.5), one has that

$$ (2.6) \qquad \frac{\partial \det D\mathcal{X}}{\partial t}(t, y) = \det D\mathcal{X}(t, y) \text{div}_A u(t, y). $$
**Remark 1.3.** Compared to Tani's work [32], we only require the interface $\mathcal{C}$ to be $\mathcal{C}^{1+\alpha}$-regular and the density to be piecewise $\mathcal{C}^{\alpha}$-Hölder continuous. Furthermore, initially non-Lipschitz velocities are allowed. Indeed, for constant viscosity coefficients, the compatibility condition (1.11) reads:

$$ \tilde{\mu}\Delta u_0 + (\tilde{\mu} + \tilde{\lambda})\nabla \operatorname{div} u_0 = \nabla P(\rho_0) + \operatorname{div}(\Pi)|_{t=0}, $$

which leads to (similar to the derivation of (1.17)):

$$ \begin{aligned} \nabla u_0 = & -\frac{1}{\tilde{\mu}}\nabla(-\Delta)^{-1}\mathcal{P}\operatorname{div}(\Pi)|_{t=0} - \frac{1}{2\tilde{\mu} + \tilde{\lambda}}\nabla(-\Delta)^{-1}(\mathcal{I} - \mathcal{P})\operatorname{div}(\Pi)|_{t=0} \\ & - \frac{1}{2\tilde{\mu} + \tilde{\lambda}}\nabla(-\Delta)^{-1}\nabla(P(\rho_0, c_0) - \tilde{P}). \end{aligned} $$

We observe that, in general, the first two terms are not bounded for $\operatorname{div}(\Pi)|_{t=0} \in L^2(\mathbb{R}^d)$.

**Remark 1.4.** The solution $(c, \rho, u)$ constructed in the above theorem exhibits the following properties.

(1) There exists a constant $R$ depending on the lifespan $T$, on the norms of the initial data and on the initial interface $C_0$ such that:

$$ \sup_{[0,T]} \{\|u\|_{H^1(\mathbb{R}^d)} + \|\dot{u}\|_{L^2(\mathbb{R}^d)}\} \le R. $$

(2) From item 4, we have $\rho\dot{u} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d))$ and from the momentum equation (1.8) we infer that $\operatorname{div}(\Pi) \in \mathcal{C}([0, T], L^2(\mathbb{R}^d))$: the compatibility condition (1.13) is propagated in time.

(3) From item 3 and interpolation we infer that the velocity $u$ is continuous in time and space. Additionally, from item 4, its material derivative, $\dot{u}$ is continuous in space. Furthermore, from item 5 we deduce that $\ddot{u}$ is at least continuous across the interface $\mathcal{C}(t)$. Roughly, for all $j \in [1, J]$, $s \in V_j$, we have:

$$ \|\dot{u}(t, \gamma_j(t, s))\| = 0 \implies \frac{d}{dt} \|\dot{u}(t, \gamma_j(t, s))\| = \|\ddot{u}(t, \gamma_j(t, s))\| = 0. $$

(4) From item 3, we infer that $\nabla u \in L^1((0, T), \mathcal{C}_{pw,C}^{\alpha}(\mathbb{R}^d))$, and hence we have the persistence of the initial interface's regularity. Indeed $\mathcal{C}(t)$ is $\mathcal{C}^{1+\alpha}$-regular manifold of dimension $d-1$, with the following estimates:

$$ \begin{aligned} & |\mathcal{C}(t)|_* \le |\mathcal{C}_0|_* \exp \left( \int_0^t \|\nabla u(\tau)\|_{L^\infty(\mathbb{R}^d)} d\tau \right), \\ & \|\nabla \gamma_j(t)\|_{\mathcal{C}^\alpha} \le \left( \|\nabla \gamma_{0,j}\|_{\mathcal{C}^\alpha} + \|\nabla \gamma_{0,j}\|_{L^\infty}^{1+\alpha} \int_0^t \|\nabla u(\tau)\|_{\dot{\mathcal{C}}_{pw,C(\tau)}^{\alpha}(\mathbb{R}^d)} d\tau \right) \\ & \times \exp \left( (2+\alpha) \int_0^t \|\nabla u(\tau)\|_{L^\infty(\mathbb{R}^d)} d\tau \right). \end{aligned} \quad (1.18) $$

(5) Applying $\partial_t \cdot + \operatorname{div}(\cdot, u)$ to the momentum equations (1.8)$_3$, we find:

$$ (1.19) \qquad \rho \ddot{u} = \operatorname{div}(2\mu(\rho, c)\mathbb{D}\dot{u}) + \nabla(\lambda(\rho, c)\operatorname{div}\dot{u}) + \operatorname{div}\mathcal{F}, $$

where $\mathcal{F}$ is a lower order term that belongs, at least, to $L^2((\tau, T), \mathcal{C}_{pw,C}^{\alpha}(\mathbb{R}^d))$, for all $\tau \in (0, T)$. With the help of item 5, small fluctuation assumption on $\mu(\rho, c)$, quantitative Hölder estimate Lemma A.1, we infer from (1.19) that:

$$ \nabla \dot{u} \in L^2((\tau, T), \mathcal{C}_{pw,C}^{\beta}(\mathbb{R}^d)) \quad \text{with} \quad \beta = \begin{cases} \alpha, & \text{in } d=2, \\ \min\{\alpha, 1/2\}, & \text{in } d=3, \end{cases} $$

for all $\tau \in (0, T)$. In particular, $\nabla \dot{u}$ is also continuous on both sides of $\mathcal{C}$.

**Theorem 1.1** immediately implies:

**Corollary 1.5** (Blow-up criterion). Let $(c, \rho, u)$ be the solution constructed in **Theorem 1.1** defined up to a maximal time $T^*$. If $T^* < \infty$, then we have either:

$$ \limsup_{t \to T^*} \left\{ |\mathcal{C}(t)|_* + \|\mathcal{C}(t)\|_{\operatorname{Lip}} + \|\mathcal{C}(t)\|_{\dot{\mathcal{C}}^{1+\alpha}} + \left\| \frac{1}{\rho(t)}, \frac{1}{\mu(\rho, c)(t)} \right\|_{L^\infty(\mathbb{R}^d)} \right\} \\ + \limsup_{t \to T^*} \left\{ \|u(t)\|_{H^1(\mathbb{R}^d)} + \|(\rho\dot{u})(t)\|_{L^2(\mathbb{R}^d)} + \|P(\rho, c)(t) - \tilde{P}\|_{\mathcal{C}_{pw,C(t)}^{\alpha}(\mathbb{R}^d)} \right\} = \infty, $$
Finally, if we combine (2.1) and the previous identity, one notices that the continuity equation on the density reads:

$$ \frac{\partial(\bar{p} \det D\mathcal{X})}{\partial t}(t, y) = 0 $$

and consequently

$$ \bar{p}(t, y) = \frac{\rho_0(y)}{\det D\mathcal{X}(t, y)}. $$

Let us turn now to the rewriting of the momentum equation in Lagrangian coordinates. First, from what precedes, it is clear that

$$ [\partial_t(\rho u^j) + \operatorname{div}(\rho u^j u)](t, \mathcal{X}(t, y)) = \frac{1}{\det D\mathcal{X}(t, y)} \rho_0(y) \partial_t \bar{u}(t, y) $$

so it only remains to compute the term of the right hand side of the momentum equation. From the conservative form (2.4), one has that:

$$ \begin{aligned} \overline{\det \Pi^j} &= \frac{1}{\det D\mathcal{X}(t, y)} \partial_k \{ (\operatorname{Adj} D\mathcal{X})_{kl} \overline{\Pi_l^j} \} \\ &= \frac{1}{\det D\mathcal{X}(t, y)} \partial_k \{ (\operatorname{Adj} D\mathcal{X})_{kl} \{ 2\mu(\bar{p}, c_0) \mathbb{D}_A^{jl} \bar{u} + (\lambda(\bar{p}, c_0) \operatorname{div}_A \bar{u} - P(\bar{p}, c_0) + \tilde{P}) \delta_{jl} \} \} \end{aligned} $$

where, by the help of the non conservative form (2.3), one has

$$ 2\mathbb{D}_A v = Dv \cdot A + A^T \cdot \nabla v. $$

This ends the transformation of the momentum equation in Lagrangian coordinates. Finally, if $(c, \rho, u)$ solves the system (1.5), then $(\bar{c}, \bar{\rho}, \bar{u})$ solves:

$$ (2.7) \qquad \left\{ \begin{aligned} \partial_t \bar{c} &= 0, \\ \partial_t \{ \bar{p} J \} &= 0, \\ \rho_0 \partial_t \bar{u} &= \operatorname{div} \{ \operatorname{Adj}(D\mathcal{X}) \{ 2\mu(\bar{p}, c_0) \mathbb{D}_A \bar{u} + (\lambda(\bar{p}, c_0) \operatorname{div}_A \bar{u} - P(\bar{p}, c_0) + \tilde{P}) Id \} \} \end{aligned} \right. $$

where $J := |\det D\mathcal{X}|$.

### 3. STUDY OF THE LINEARISED SYSTEM

#### 3.1. Statement of the main result.
The aim of this section is to prove the well-posedness of the linear non-homogeneous parabolic equation

$$ (3.1) \qquad \rho_0 u_t - \operatorname{div}(2\mu \mathbb{D}u + \lambda \operatorname{div} u I_d) = \operatorname{div} F $$

associated to (2.7). The initial density $\rho_0 = \rho_0(x)$ and the viscosity $\mu = \mu(t, x)$, $\lambda = \lambda(t, x)$, are supposed to verify:

$$ (3.2) \qquad \left\{ \begin{aligned} & 0 < \rho_* \le \rho_0 \le \rho^*, \quad 0 < \mu_* \le \mu, \quad \lambda \le \mu^*, \\ & \mu, \lambda \in L^\infty((0,T), \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)), \\ & \mu_t, \lambda_t \in L^\infty((0,T), L^4(\mathbb{R}^d)) \cap L^1((0,T), L^\infty(\mathbb{R}^d)) \cap L^{16}((0,T), L^8(\mathbb{R}^d)), \\ & \mu_{tt}, \lambda_{tt} \in L^1((0,T), L^2(\mathbb{R}^d)), \quad \sqrt{\sigma}(\mu_{tt}, \lambda_{tt}) \in L^{8/d}((0,T), L^4(\mathbb{R}^d)), \end{aligned} \right. $$

for positive constants $\rho_*, \rho^*, \mu_*, \mu^*$ and time $0 < T < \infty$. Above, $\gamma$ is a parametrization of an hyper-surface $C \in \mathcal{C}^{1+\alpha}$ in $\mathbb{R}^d$, $d \in \{2, 3\}$:

$$ C := \{ \gamma(s) \in \mathbb{R}^d : s \in \mathbb{T}^{d-1} \} $$

that verifies the following non self-intersection condition:

$$ \exists c_\gamma > 0, \quad \forall (s, s') \in (\mathbb{T}^{d-1})^2, \quad |\gamma(s) - \gamma(s')| \ge c_\gamma^{-1} |s - s'|. $$

We refer to **Notation 1.1** for the definition of the space of piecewise Hölder continuous functions. The source term $F$ in (3.1) and the initial datum $u_0$ are such that $(F, u_0) \in Y_T(\mu_0, \lambda_0)$ where:

$$ \begin{aligned} Y_T(\mu_0, \lambda_0) := & \{(F, u_0) \in \mathcal{C}([0, T], L^2(\mathbb{R}^d)) \times H^1(\mathbb{R}^d) : F \in L^\infty((0, T), L^6(\mathbb{R}^d)), \\ & F \in L^{16}((0, T), L^8(\mathbb{R}^d)), \sigma^{r/4} F \in L^4((0, T), \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)), \\ & F_t \in L^2((0, T) \times \mathbb{R}^d), \sqrt{\sigma} F_t \in L^{8/d}((0, T), L^4(\mathbb{R}^d)), \\ & \sigma F_{tt} \in L^2((0, T) \times \mathbb{R}^d), \operatorname{div}(2\mu_0 \mathbb{D}u_0 + \lambda_0 \operatorname{div} u_0 I_d - F|_{t=0}) \in L^2(\mathbb{R}^d) \}. \end{aligned} $$
or

$$ (1.20) \quad \limsup_{t \to T^*} \left[ 1 + \|\lambda(\rho, c)(t)\|_{\mathcal{C}_{pw,C(t)}^{\alpha}(\mathbb{R}^d)} + (\mathfrak{P}_{C(t)} + \ell_{\varphi(t)}^{-\alpha}) \|\!\|\!\|\lambda(\rho, c)(t)\!\|\!\|\!\|_{L^{\infty}(C(t))} \right] \|\mu(\rho, c)(t) - \tilde{\mu}\|_{\mathcal{C}_{pw,C(t)}^{\alpha}(\mathbb{R}^d)} \\ + \limsup_{t \to T^*} \left( \mathfrak{P}_{C(t)} + \ell_{\varphi(t)}^{-\alpha} \right) \left[ \|\!\|\!\|\mu(\rho, c)(t)\!\|\!\|\!\|_{L^{\infty}(C(t))} + \|\!\|\!\|\mu(\rho, c)(t)\!\|\!\|\!\|_{L^{\infty}(C(t))} \left\| 1 - \frac{\tilde{\mu}}{\langle\mu(\rho, c)(t)\rangle} \right\|_{L^{\infty}(C(t))} \right] \geq [\mu]. $$

**Remark 1.6.** Our result can be seen as a first step towards the existence of global solutions with small data which we plan to attack in a forthcoming paper. The construction of discontinuous weak-solutions for compressible models with variable viscosity cannot be achieved in the classical way by regularizing the initial data and applying some classical results à la Matsmura-Nishida since this procedure destroys the discontinuities in the density. Some difficult constructions are then needed in order to obtain discontinuous weak-solutions, see [14, Section 3]. We believe that since we dispose of powerful functional analysis tools, see Appendix A, an approach where local-well posedness is first proved to be quite natural. With Theorem 1.1 in hand, in order to attack the global existence issue, we do not have to worry anymore about constructing appropriate approximate systems that preserve the main feature of the system (1.8).

**Remark 1.7.** The difficulty in proving global existence for the two-fluid model lies in the fact that we are not able to estimate the $L^1((0,T), \mathcal{C}_{pw,C}^{\alpha}(\mathbb{R}^d))$ norm of $\nabla u$ uniformly with respect to time. In particular, see (1.18), the characteristic $|C(t)|_*$ and the $\mathcal{C}^{1+\alpha}$-norm of the interface $\gamma_j(t)$ are expected to grow exponentially with respect to time which prevents the validity of the condition (1.20) globally in time.

In the study of the Navier-Stokes equations with constant viscosity for $d=2$, Hoff [14] has observed that the exponential-in-time growth of these quantities can be compensated by the exponential decay of the pressure jump across the interface (in our case we would be needing the exponential-in-time decay of the viscosity jump, see the explosion criterion (1.20)). In order to obtain exponential decay, we use the continuity of the normal component of the stress tensor, we follow the arguments starting from (3.38) below to establish:

$$ \|\!\|\!\| (2\mu(\rho, c) + \lambda(\rho, c)) \operatorname{div} u - P(\rho, c) \|\!\|\!\| = 2 \|\!\|\!\| \mu(\rho, c) \|\!\|\!\| (\langle \operatorname{div} u \rangle - \langle \mathbb{D}^{jk} u \rangle n_x^j n_x^k). $$

Next, we observe that $f = f(\rho, c)$ defined up to a constant as

$$ f(\rho, c) = \int^{\rho} \frac{2\mu(s, c) + \lambda(s, c)}{s} ds $$

satisfies:

$$ \partial_t f(\rho, c) + u \cdot \nabla f(\rho, c) + P(\rho, c) - \tilde{P} = -[(2\mu(\rho, c) + \lambda(\rho, c)) \operatorname{div} u - P(\rho, c)]. $$

We write the above equation along the flow $X$ of the velocity and then take jumps to obtain:

$$ (1.21) \quad \frac{d}{dt} \|f(\rho, c)\|(t, X(t, s)) + \|P(\rho, c)\|(t, X(t, s)) = 2 \|\!\|\!\| \mu(\rho, c) \|\!\|\!\| (\langle \operatorname{div} u \rangle - \langle \mathbb{D}^{jk} u \rangle n_x^j n_x^k) (t, X(t, s)), \quad s \in \mathcal{C}_0. $$

Observe that in the constant viscosity case, the RHS of the above equation vanishes and using the fact that the pressure is strictly increasing one infers in a straightforward manner the exponential-in-time decay of the jump of $f$. For the two-fluid model, when the pressure and viscosity laws of the two phases are different, it is less obvious how to obtain exponential decay of the jump or even if this property holds in general. First of all the RHS of the above identity is no longer 0 and second of all in order to obtain exponential-in-time decay, the jump of $f$, $P$ and $\mu$ need to be "comparable", (this is the case for constant $\mu$ and increasing $P = P(\rho)$). For multi-phase flows the minimal requirements seem to be that:

$$ (1.22) \quad 0 < v_* \le \frac{P(\rho, 1) - P_-(\rho', 0)}{f(\rho, 1) - f(\rho', 0)} \le v^* \quad \text{and} \quad \frac{|\mu(\rho, 1) - \mu(\rho', 0)|}{|f(\rho, 1) - f(\rho', 0)|} + \frac{|\lambda(\rho, 1) - \lambda(\rho', 0)|}{|f(\rho, 1) - f(\rho', 0)|} \le v^*, $$

at least, for all $\rho$ close to $\tilde{\rho}(1)$ and $\rho'$ close to $\tilde{\rho}(0)$ where $\tilde{P} = P(\tilde{\rho}(1), 1) = P(\tilde{\rho}(0), 0)$. This requirement does not seem to occur in general. For instance, it does not hold in the case of two different constant viscosity coefficients say $\mu(\rho, 1) = \tilde{\mu}$, $\mu(\rho, 0) = \mu_-$, where $\mu_- \ne \tilde{\mu}$.

**Outline of the paper.** The remainder of this paper is divided into three sections. In the first one, Section 2, we rewrite system (1.8) in Lagrangian coordinates, see (2.7) below, and in the second one, Section 3, we study the linearized system associated to (2.7). Finally, the study of the full nonlinear system in Lagrangian coordinates (2.7) is presented in Section 4 from which Theorem 1.1 follows.
Above, $\mu_0 := \mu_{|t=0}$, $\lambda_0 := \lambda_{|t=0}$ and $0 < r < 2$ depends on $0 < \alpha < 1$ and the dimension $d \in \{2, 3\}$. This dependence is explicitly given in (3.32). The existence and uniqueness of solution to (3.1) is established in the following space:

$$X_T := \{u \in \mathcal{C}([0, T], H^1(\mathbb{R}^d)) : \begin{array}{l} \nabla u \in L^\infty((0, T), L^6(\mathbb{R}^d)) \cap L^{16}((0, T), L^8(\mathbb{R}^d)), \sigma^{r/4} \nabla u \in L^4((0, T), \mathcal{C}_{pw,\gamma}^\alpha(\mathbb{R}^d)), \\ u_t \in \mathcal{C}([0, T], L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)), \sqrt{\sigma} \nabla u_t \in L^{8/d}((0, T), L^4(\mathbb{R}^d)) \cap L^\infty((0, T), L^2(\mathbb{R}^d)), \\ \sqrt{\sigma} u_{tt} \in L^2((0, T) \times \mathbb{R}^d), \sigma u_{tt} \in L^\infty((0, T), L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)) \end{array}\}.$$

The time weights are inspired by the pioneer work of Hoff [10]. One of the main difficulties is to obtain estimates for the gradient of $u$, of the equations (3.1), allowing to propagate the regularity of the interface $\gamma$. The difficulty stems from the fact that the viscosity do not have any Sobolev regularity.

The main result of this section is stated in the following theorem:

**Theorem 3.1.** Let $0 < T < \infty$, the viscosity $\mu$, $\lambda$ and the initial density $\rho_0$ verify (3.2), $0 < \alpha < 1$. There exists a constant $[\mu]$ depending on the dimension $d$, on $\alpha$ and on the interface $\gamma$, such that if:

$$(3.4) \qquad \frac{1}{\mu} \|\mu - \tilde{\mu}\|_{L^{\infty}((0,T), \mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d))} \le [\mu],$$

then, for all $(F, u_0) \in Y_T(\mu_0, \lambda_0)$, there exists an unique solution of the Cauchy problem:

$$(3.5) \qquad \begin{cases} \rho_0 u_t - \operatorname{div}(\mu \mathbb{D} u) - \nabla(\lambda \operatorname{div} u) = \operatorname{div} F, \\ u|_{t=0} = u_0, \end{cases}$$

in $X_T$. Moreover, the following estimates hold true for some constant $M$ that depends on the upper and lower bounds of the density and viscosity (see first line of (3.2)):

(1) Estimates for $u$ and $u_t$:

$$\begin{align} \sup_{[0,T]} \|u, u_t, \nabla u\|_{L^2(\mathbb{R}^d)}^2 &+ \|\nabla u\|_{L^\infty((0,T), L^6(\mathbb{R}^d))}^2 + \int_0^T \|\nabla u, \nabla u_t, u_t\|_{L^2(\mathbb{R}^d)}^2 \nonumber \\ &\le M \left( \| (F, u_0)\|_{Y_T(\mu_0, \lambda_0)}^2 + \int_0^T \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F\|_{L^2 \cap L^6(\mathbb{R}^d)}^2 \right) \nonumber \\ &\quad \times \exp \left( M \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \right] \right). \tag{3.6} \end{align}$$

(2) Estimates for $u_{tt}$:

$$\begin{align} \int_0^T \|\sqrt{\sigma} u_{tt}, \sigma \nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \|\sqrt{\sigma} \nabla u_t, \sigma u_{tt}\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^2 &\le M \left[ \int_0^{\sigma(T)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|F_t, \sigma F_{tt}\|_{L^2(\mathbb{R}^d)}^2 \right. \nonumber \\ &\qquad \left. + \int_0^T \|\mu_t, \lambda_t, \sigma(\mu_{tt}, \lambda_{tt})\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \right] \nonumber \\ \tag{3.7} &\quad \times \exp \left( M \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \|\mu_{tt}, \lambda_{tt}\|_{L^2(\mathbb{R}^d)} \right] + M \int_0^T \sigma \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)}^2 \right). \end{align}$$

The norms of $u$ involved in this estimate are controlled thanks to (3.6).

**Remark 3.2.** The first estimate (3.6) follows from the classical energy estimate and the so-called Hoff estimates introduced in [10, 11]. These estimates provide $H^1$ regularity for the time derivative of the velocity $u$. The other estimate (3.7) can be seen as a higher Hoff estimate and has been used in [5] in the study of the incompressible Navier-Stokes model.

**Corollary 3.3.** The solution in the above theorem comes with subsequent estimates:

(1) $L^{16}((0,T), L^8(\mathbb{R}^d))$ estimate for the $\nabla u$:

$$(3.8) \qquad \|\nabla u\|_{L^{16}((0,T), L^8(\mathbb{R}^d))} \lesssim \|F\|_{L^{16}((0,T), L^8(\mathbb{R}^d))} + \|F, u_t, \nabla u\|_{L^2 \cap L^\infty((0,T), L^2(\mathbb{R}^d))} + \|\nabla u_t\|_{L^2((0,T) \times \mathbb{R}^d)}.$$
## 2. CHANGE OF VARIABLES IN LAGRANGIAN COORDINATES

The aim of this section is to write the system (1.8) in Lagrangian coordinates. We know that for a Lipschitz vector field $u$, one can define a flow map $\mathcal{X}$ given by the Duhamel formula:

$$ \mathcal{X}(t, y) = y + \int_{0}^{t} u(s, \mathcal{X}(s, y))ds $$

and adopting the following notation: for all function $g = g(t, x)$ we denote by $\bar{g}$ the function:

$$ \bar{g}(t, y) := g(t, \mathcal{X}(t, y)) $$

then, the flow of the vector field $u$ can be written in term of $\bar{u}$ by:

$$ \mathcal{X}(t, y) := \mathcal{X}_{\bar{u}}(t, y) = y + \int_{0}^{t} \bar{u}(s, y)ds. $$

And for any $g = g(t, x)$, one has that

$$ (2.1) \qquad \partial_t \bar{g}(t, y) = (\partial_t g + u \cdot \nabla g)(t, \mathcal{X}_{\bar{u}}(t, y)). $$

We introduce some notations.

**Notation 2.1.** Let us fix a square matrix $B = (b_{ij})_{ij} \in \mathcal{M}_n(\mathbb{R})$, $n \in \mathbb{N}^*$.

* The trace of $B$ is the sum of the diagonal elements of $B$:

$$ \mathrm{Tr}(B) = b_{ii}. $$

The transposed matrix of $B$, denoted $B^T$, is:

$$ B^T = (b_{ji})_{ij} $$

and for $B' = (b'_{ij})_{ij} \in \mathcal{M}_n(\mathbb{R})$, we have:

$$ \mathrm{Tr}(B \cdot B') = b_{ij} b'_{ji} = B : B' = (B')^T : B^T. $$

* For $i, j \in [1, n]$, $M_{ij}(B)$ denotes the determinant of the $(n-1) \times (n-1)$ matrix obtained by removing the $i$-th row and $j$-th column of $B$.

* The adjugate matrix of $B$, denoted as $\mathrm{Adj}(B)$, is given by:

$$ (\mathrm{Adj}(B))_{ij} = (-1)^{i+j} M_{ji}(B). $$

* If $\det(B) \neq 0$, then $B$ is invertible and we have:

$$ B^{-1} = \frac{1}{\det(B)} \mathrm{Adj}(B). $$

If $DG$ denotes the Jacobian matrix of $G$, then chain rule yields:

$$ D\overline{G} = DG \circ \mathcal{X}_{\bar{u}} \cdot D\mathcal{X}_{\bar{u}}, $$

but, since

$$ (2.2) \qquad (D\mathcal{X}_{\bar{u}}(t, y))_{ij} = \delta^{ij} + \int_{0}^{t} \partial_{j}\bar{u}^{i}(s, y)ds, $$

it holds, at least for small time, that the Jacobian matrix of the change of variable is invertible and hence:

$$ \overline{DG}(t, y) = D\overline{G}(t, y) \cdot (D\mathcal{X}_{\bar{u}}(t, y))^{-1} = \frac{1}{\det D\mathcal{X}_{\bar{u}}(t, y)} D\overline{G}(t, y) \cdot \mathrm{Adj} D\mathcal{X}_{\bar{u}}(t, y). $$

In particular, one has that:

$$ (2.3) \qquad \overline{\partial_j G^i} = \frac{1}{\det D\mathcal{X}_{\bar{u}}} \partial_k \overline{G}^i (\mathrm{Adj} D\mathcal{X}_{\bar{u}})_{kj}. $$

On the other hand, from Piola's formula: for any matrix $A \in \mathcal{M}_{n \times n}(\mathbb{R})$,

$$ \partial_k (\mathrm{Adj}(A))_{kj} = 0 $$

we infer that:

$$ (2.4) \qquad \overline{\partial_j G^i} = \frac{1}{\det D\mathcal{X}_{\bar{u}}} \partial_k \overline{G}^i (\mathrm{Adj} D\mathcal{X}_{\bar{u}})_{kj}. $$
$$ (2) \quad L^{8/d}((0,T), L^4(\mathbb{R}^d)) \text{ norm estimate for } \sqrt{\sigma\nabla}u_t: $$

$$ (3.9) \quad \int_0^T \sigma^{4/d} \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \lesssim \int_0^T \sigma^{4/d} \|F_t\|_{L^4(\mathbb{R}^d)}^{8/d} + \int_0^T \sigma^{4/d} \|\nabla u\|_{L^8(\mathbb{R}^d)}^{8/d} \|\mu_t, \lambda_t\|_{L^8(\mathbb{R}^d)}^{8/d} + \|\sigma u_{tt}\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^{8/d} \\ \qquad + \left[ \int_0^T \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathbb{D}u\|_{L^2(\mathbb{R}^d)}^{2} \right]^{4/d} . $$

$$ (3) \text{ Hölder estimate for } \nabla u: $$

$$ (3.10) \quad \int_0^T \sigma^r \|\nabla u\|_{\mathcal{E}_{\mu,\nu}^{\alpha}(\mathbb{R}^d)}^4 \lesssim (1+\mathcal{P}(C)) \int_0^T \sigma^r \|F\|_{\mathcal{E}_{\mu,\nu}^{\alpha}(\mathbb{R}^d)} + \|F, \nabla u, u_t, \sqrt{\sigma}\nabla u_t\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^4 \\ \qquad + \left[ \int_0^T \|F, \nabla u, u_t, \nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \right]^2 + \left[ \int_0^T \sigma^{4/d} \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \right]^{d/2} . $$

**3.2. Sketch of the Proof.** The proof of the above theorem is done in Section 3.3.3 and is based on a homotopy argument. For this purpose, we define for $\theta \in [0, 1]$

$$ \mu^{\theta} = (1 - \theta)\tilde{\mu} + \theta\mu \quad \text{and} \quad \lambda^{\theta} = (1 - \theta)\tilde{\lambda} + \theta\lambda $$

and consider the equation (3.5) where the viscosity $\mu$, $\lambda$ are replaced by $\mu^\theta$ and $\lambda^\theta$ respectively:

$$ (3.11) \quad \begin{cases} \rho_0 u_t - \operatorname{div}(2\mu^\theta \mathbb{D}u + \lambda^\theta \operatorname{div} u I_d) = \operatorname{div} F, \\ u|_{t=0} = u_0. \end{cases} $$

An important remark is that the last condition in the definition of the set $Y_T(\mu_0, \lambda_0)$ complicates the proof. In fact, it is not obvious that if $(F, u_0)$ belongs to $Y_T(\mu_0^\theta, \lambda_0^\theta)$ then it also belongs to $Y_T(\mu_0^{\theta'}, \lambda_0^{\theta'})$. In general, this is not the case unless $F|_{t=0}$, $\nabla u_0$, $\lambda|_{t=0}$ and $\mu|_{t=0}$ have a Sobolev regularity. Of course, this condition does not hold true in our context because $F$ and the viscosity are assumed to be discontinuous across a hyper-surface so that their weak gradient are just a measure. In order to ease the reading, we set $Y_T^\theta := Y_T(\mu_0^\theta, \lambda_0^\theta)$. The principle of homotopy argument is to look for $\theta \in [0, 1]$ such that (3.11) is uniquely solvable in $X_T$. So we define the following set:

$$ (3.12) \quad \mathcal{E}_T := \{\theta \in [0,1]: \forall (F, u_0) \in Y_T^\theta, \exists! u \in X_T \text{ solution of (3.11)}\}. $$

It is obvious that **Theorem 3.1** follows as soon as $1 \in \mathcal{E}$. This will be obtained by proving that $\mathcal{E}$ is at the same time a non-empty, open and closed set of $[0, 1]$. The proof is structured into three steps: In the first one, **Section 3.3.1**, we derive some a priori estimates on a solution $u \in X_T$ of the equation (3.5). In the second step, we prove that $0 \in \mathcal{E}$ and finally, we prove in **Section 3.3.3** that $\mathcal{E}$ is an open and a closed set of $[0, 1]$ from which, **Theorem 3.1** follows. In the following, we will discuss a little bit the main difficulties in proving that $0 \in \mathcal{E}$.

When $\theta = 0$, the equation (3.11) is reduced to the following:

$$ (3.13) \quad \begin{cases} \rho_0 u_t - \tilde{\mu} \Delta u - (\tilde{\mu} + \tilde{\lambda}) \nabla \operatorname{div} u = \operatorname{div} F, \\ u|_{t=0} = u_0 \end{cases} $$

and the aim is to prove that for all $(F, u_0) \in Y_T^0$ there exists an unique solution of (3.13) in $X_T$. This is not obvious since, (3.13) is not really a heat equation and also the initial density is rough. Despite, one can add $(\tilde{\rho} - \rho_0)u_t$ in both sides of (3.13), in order to obtain a heat equation with a right hand side that belongs to $L^2((0, T), L^2(\mathbb{R}^d)) + L^2((0, T), \dot{H}^{-1}(\mathbb{R}^d))$. As in [4], we may require a smallness of $\tilde{\rho} - \rho_0$ in order to obtain a stability for approximate solution. In order to get around this difficulty, the strategy adopted is that, we consider a mollifying sequence of the initial density $\rho_0$, velocity $u_0$ and $F$ respectively $\rho_{0,\delta}$ $u_{0,\delta}$, and $F_\delta$ and consider the equation:

$$ (3.14) \quad \begin{cases} \rho_{0,\delta} u_t - \tilde{\mu} \Delta u - (\tilde{\mu} + \tilde{\lambda}) \nabla \operatorname{div} u = \operatorname{div} F_\delta, \\ u|_{t=0} = u_{0,\delta}. \end{cases} $$

The existence and uniqueness of global solution $u^\delta$ of the above system is based on, a change of unknowns, a study of the heat kernel and a homotopy argument. Notice that the homotopy argument applied in this context seems simpler than the one we are applying in order to deduce the proof of **Theorem 3.1**. Taking the advantage that
The first form (2.3) will be referred as the non-conservative form and the second (2.4) is the conservative one and both will be used. From the non-conservative form, we find

$$ \overline{\det G} = \frac{1}{\det D\mathcal{X}_{\bar{u}}} \mathrm{Tr}(D\overline{G} \cdot \mathrm{Adj} D\mathcal{X}_{\bar{u}}), $$

and by setting

$$ A_{\bar{u}}(t, y) = (D\mathcal{X}_{\bar{u}}(t, y))^{-1} $$

we have

$$ (2.5) \qquad \mathrm{div}_{A_{\bar{u}}} G = \overline{\mathrm{div} G} = D\overline{G} : A_{\bar{u}} = A_{\bar{u}}^{T} : \nabla \overline{G}. $$

On the other hand, for all matrix $B \in \mathscr{C}^1([0, T], \mathcal{M}_n(\mathbb{R}))$, one has:

$$ \begin{aligned} \frac{d \det B}{ds}(s) &= \frac{\partial \det B}{\partial b_{ij}} \frac{\partial b_{ij}}{\partial s}(s) \\ &= (-1)^{i+j} M_{ij}(B(s)) \frac{\partial b_{ij}}{\partial s}(s) \\ &= \mathrm{Adj}(B(s))_{ji} \left( \frac{\partial B}{\partial s}(s) \right)_{ij} \\ &= \frac{\partial B(s)}{\partial s} : \mathrm{Adj}(B(s)). \end{aligned} $$

In particular, for $B = D\mathcal{X}_{\bar{u}}$, and remembering from (2.2) that

$$ \frac{\partial D\mathcal{X}_{\bar{u}}}{\partial t} = D\bar{u} $$

we find

$$ \frac{\partial \det D\mathcal{X}_{\bar{u}}}{\partial t}(t, y) = D\bar{u}(t, y) : \mathrm{Adj} D\mathcal{X}_{\bar{u}}(t, y). $$

Hence, the identity (2.5) yields:

$$ (2.6) \qquad \frac{\partial \det D\mathcal{X}_{\bar{u}}}{\partial t}(t, y) = \det D\mathcal{X}_{\bar{u}}(t, y) \mathrm{div}_{A_{\bar{u}}} u(t, y). $$

Finally, combining (2.1) and the previous identity, we notice that the continuity equation on the density reads:

$$ \frac{\partial(\bar{\rho} \det D\mathcal{X}_{\bar{u}})}{\partial t}(t, y) = 0 $$

and consequently

$$ \bar{\rho}(t, y) = \frac{\rho_0(y)}{\det D\mathcal{X}_{\bar{u}}(t, y)}. $$

Now, let us proceed to rewrite the momentum equations in Lagrangian coordinates. First, from what precedes, it is clear that

$$ [\partial_t(\rho u^j) + \mathrm{div}(\rho u^j u)](t, \mathcal{X}_{\bar{u}}^j(t, y)) = \frac{1}{\det D\mathcal{X}_{\bar{u}}(t, y)} \rho_0(y) \partial_t \bar{u}(t, y) $$

so it only remains to compute the term involving the stress tensor. From the conservative form (2.4), we have:

$$ \begin{aligned} \overline{\mathrm{div}} \Pi^j &= \frac{1}{\det D\mathcal{X}_{\bar{u}}(t, y)} \partial_k \{ (\mathrm{Adj} D\mathcal{X}_{\bar{u}})_{kl} \Pi_l^j \} \\ &= \frac{1}{\det D\mathcal{X}_{\bar{u}}(t, y)} \partial_k \{ (\mathrm{Adj} D\mathcal{X}_{\bar{u}})_{kl} \{ 2\mu(\bar{\rho}, c_0) \mathbb{D}_{A_{\bar{u}}}^{jl} \bar{u} + (\lambda(\bar{\rho}, c_0) \mathrm{div}_{A_{\bar{u}}} \bar{u} - P(\bar{\rho}, c_0) + \tilde{P}) \delta_{jl} \} \} \end{aligned} $$

where, we use the following notation:

$$ 2\mathbb{D}_B v = Dv \cdot B + B^T \cdot \nabla v. $$

This ends the transformation of the momentum equation in Lagrangian coordinates. Summing up, we find that $(\bar{c}, \bar{\rho}, \bar{u})$ solves:

$$ (2.7) \qquad \left\{ \begin{array}{l} \partial_t \bar{c} = 0, \\ \partial_t \{ \bar{\rho} J_{\bar{u}} \} = 0, \\ \rho_0 \partial_t \bar{u} = \mathrm{div} \left\{ \mathrm{Adj}(D\mathcal{X}_{\bar{u}}) \left\{ 2\mu(\bar{\rho}, c_0) \mathbb{D}_{A_{\bar{u}}} \bar{u} + (\lambda(\bar{\rho}, c_0) \mathrm{div}_{A_{\bar{u}}} \bar{u} - P(\bar{\rho}, c_0) + \tilde{P}) I_d \right\} \right\}, \end{array} \right. $$

where $J_{\bar{u}} := |\det D\mathcal{X}_{\bar{u}}|$. We establish the well-posedness of (2.7) in the following theorem.
for any $\delta > 0$, $\rho_{0,\delta}$ is smooth enough and bounded away from vacuum, one has that $u$ solves (3.14) if and only if $v := \rho_{0,\delta}u$ solves the following Cauchy problem:

$$ (3.15) \qquad \begin{cases} v_t^j - \operatorname{div} \left( \frac{\tilde{\mu}}{\rho_{0,\delta}} \mathbb{D}^j v \right) - \partial_j \left( \frac{\tilde{\lambda}}{\rho_{0,\delta}} \operatorname{div} v \right) = \operatorname{div} F_\delta^j - \tilde{\mu} \partial_k \left( \rho_{0,\delta}^{-2} (v^k \partial_j \rho_{0,\delta} + v^j \partial_k \rho_{0,\delta}) \right) - \tilde{\lambda} \partial_j \left( \rho_{0,\delta}^{-2} v \cdot \nabla \rho_{0,\delta} \right), \\ v|_{t=0} = \rho_{0,\delta} u_{0,\delta}. \end{cases} $$

In order to ease computations, we consider the following problem for $\mu'_0, \lambda'_0, b_0 = \mu'_0(x), \lambda'_0(x), b_0(x)$

$$ (3.16) \qquad \begin{cases} v_t - \operatorname{div}(2\mu'_0 \mathbb{D}v + \lambda'_0 \operatorname{div} v I_d) = f + \operatorname{div}(b_0 : v), \\ v|_{t=0} = v_0 \end{cases} $$

and we have the following:

**Theorem 3.4.** Let $0 < T < \infty$ and $f \in L^2((0, T) \times \mathbb{R}^d)$, $0 < \mu_* \le \mu'_0$, $\lambda'_0 \in W^{1,\infty}(\mathbb{R}^d)$, $b_0 \in W^{1,\infty}(\mathbb{R}^d)$, and finally $v_0 \in H^1(\mathbb{R}^d)$. Then there exists an unique solution $v \in \mathcal{C}([0, T], H^1(\mathbb{R}^d))$ with $v_t, \nabla^2 v \in L^2((0, T), L^2(\mathbb{R}^d))$ for the system (3.16). Moreover, there exist constants $\Lambda, K > 0$ that depend on the norms of $\mu'_0, \lambda'_0, b_0$ and the lower and upper bounds $\mu_*, \mu^*$ such that

$$ \sup_{[0,T]} e^{-2\Lambda t} \left\{ \|v(t)\|_{H^1(\mathbb{R}^d)}^2 + \int_0^t \|v_t, \nabla v, \nabla^2 v\|_{L^2(\mathbb{R}^d)}^2 \right\} \le K \left( \|v_0\|_{H^1(\mathbb{R}^d)}^2 + \int_0^T \|f(s)\|_{L^2(\mathbb{R}^d)}^2 ds \right). $$

The proof of Theorem 3.4 is given in Section 3.3.2.

**Remark 3.5.** One deduces from Theorem 3.4, (3.15) and then (3.14) that the Cauchy problem (3.14) admits a unique solution $u^\delta \in \mathcal{C}([0, T], H^1(\mathbb{R}^d))$ with $u_t^\delta, \nabla^2 u^\delta \in L^2((0, T) \times \mathbb{R}^d)$. Moreover, by deriving the equation (3.14) in space, and noticing that $\rho_0^\delta$ and $u_0^\delta$ are smooth, one obtains that for all $s \ge 2$, $u^\delta \in \mathcal{C}([0, T], H^s(\mathbb{R}^d))$, and $u_t, \nabla^2 u^\delta \in L^2((0, T), H^{s-1}(\mathbb{R}^d))$. In the same manner, one can derive the equation (3.14) with respect to time, in order to obtain the equation satisfied by $u_t^\delta$. Of course, $F_t^\delta$ is regular enough and the regularity $u_t^\delta \in L^2((0, T), H^{s-1}(\mathbb{R}^d))$ is sufficient in order to derive the equation in a classical sense. Thus, applying again Theorem 3.4 one deduces from the previous discussion that $u_t^\delta \in \mathcal{C}([0, T], H^s(\mathbb{R}^d))$ and $u_{tt}^\delta, \nabla u_t^\delta \in L^2((0, T), H^{s-1}(\mathbb{R}^d))$. This closes the existence of global smooth solution of the Cauchy problem (3.14).

We continue the sketch of the proof of Theorem 3.1 in the following.

The sequence of solution $(u^\delta)_\delta$ of the system (3.14) verifies the estimates for $u, u_t$, (in (3.6)) $u_{tt}$ (in (3.7)) and the $L^p$ estimates in Corollary 3.3 for $u$ and $u_t$, uniformly with respect to the regularisation parameter $\delta > 0$ with the exception of the piecewise Hölder estimate. Hence, one can let $\delta$ go to zero in order to obtain an solution of the Cauchy problem (3.13). Since $(u^\delta)_\delta$ is bounded in $L^\infty((0, T), H^1(\mathbb{R}^d))$ and its time derivative $(u_t^\delta)_\delta$ is bounded in $L^2((0, T), H^1(\mathbb{R}^d))$, it then turns out that $u \in \mathcal{C}([0, T], H^1(\mathbb{R}^d))$. On the other hand, one recalls that from the second Hoff estimate $(\rho_0 u_t^\delta)_\delta$ is bounded in $L^\infty((0, T), L^2(\mathbb{R}^d))$ and the time derivative of the equation (3.14) implies that $((\rho_0 u_t^\delta)_t)_\delta$ is bounded in $L^2((0, T), H^{-1}(\mathbb{R}^d))$. Then, from [20, Appendix C], $\rho_0 u_t$ belongs to $\mathcal{C}([0, T], L^2(\mathbb{R}^d))$ and since $\rho_0$ is independent from time, one deduces that $u_t$ is continuous in time with values in the weak topology of $L^2(\mathbb{R}^d)$. We will obtain the continuity of $u_t$ in time with values in the strong topology in two steps. In the first step, we notice that $(u_t^\delta)_\delta$ is bounded in $L^\infty((0, T), L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d))$ and $(u_{tt}^\delta)_\delta$ is bounded in $L^2((\tau, T) \times \mathbb{R}^d)$ for all $0 < \tau < T$, and since $H^1(\mathbb{R}^d)$ embeds compactly in $L^2(\mathbb{R}^d)$, one deduces from the Aubin-Lions theorem that $u_t \in \mathcal{C}((0, T], L^2(\mathbb{R}^d))$. It only remains the continuity of $u_t$ at the initial time. At this level, we know that $u$ already solves (3.13), the regularity of $u_t$ and $u_{tt}$ ensures that one can derive (3.13) in time and use $u_t$ as a test function in the weak formulation of the obtained equation in order to obtain:

$$ \frac{1}{2} \int_{\mathbb{R}^d} \rho_0 |u_t(\tau)|^2 + \int_0^\tau 2\tilde{\mu} \int_{\mathbb{R}^d} |\mathbb{D}u_t|_{L^2(\mathbb{R}^d)}^2 + \tilde{\lambda} \int_0^\tau \int_{\mathbb{R}^d} (\operatorname{div} u_t)^2 = \frac{1}{2} \int_{\mathbb{R}^d} \rho_0 |u_t|_{t=0}^2 - \int_0^\tau \int_{\mathbb{R}^d} \partial_t F^{jk} \partial_k u_t^j. $$

From which we deduce that

$$ \limsup_{\tau \to 0} \| \sqrt{\rho_0} u_t(\tau) \|_{L^2(\mathbb{R}^d)} \le \| \sqrt{\rho_0} u_{t|t=0} \|_{L^2(\mathbb{R}^d)}. $$

On the other hand, we have

$$ \| \sqrt{\rho_0} u_{t|t=0} \|_{L^2(\mathbb{R}^d)} \le \liminf_{\tau \to 0} \| \sqrt{\rho_0} u_t(\tau) \|_{L^2(\mathbb{R}^d)}, $$
**Theorem 2.2.** Assume the hypotheses in **Theorem 1.1** hold true for $(c_0, \rho_0, u_0)$, the pressure law, and the viscosity laws $\mu$, $\lambda$. Additionally, suppose that the small fluctuation assumption (1.15) also holds true, and consider $r$ as given in (1.16).

Then, there exist $T > 0$ and a unique solution $\bar{u}$ of the Cauchy problem

$$ (2.8) \qquad \left\{ \begin{array}{l} \rho_0 \partial_t \bar{u} = \text{div} \left\{ \text{Adj}(D\mathcal{X}_{\bar{u}}) \left\{ 2\mu(\rho_0 J_{\bar{u}}^{-1}, c_0) \mathbb{D}_{A_{\bar{u}}} \bar{u} + (\lambda(\rho_0 J_{\bar{u}}^{-1}, c_0) \text{div}_{A_{\bar{u}}} \bar{u} - P(\rho_0 J_{\bar{u}}^{-1}, c_0) + \tilde{P} \right\} I_d \right\}, \\ \bar{u}|_{t=0} = u_0, \end{array} \right. $$

verifying the following:

(1) $\bar{u} \in \mathcal{C}([0, T], H^1(\mathbb{R}^d)) \cap L^\infty((0, T), \dot{W}^{1,6}(\mathbb{R}^d)) \cap L^{16}((0, T), \dot{W}^{1,8}(\mathbb{R}^d))$, $\sigma^{r/4} \nabla \bar{u} \in L^4((0, T), \mathcal{C}^{\alpha}_{pw,C_0}(\mathbb{R}^d))$;

(2) $\partial_t \bar{u} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d))$, $\sqrt{\sigma} \nabla \partial_t \bar{u} \in L^\infty((0, T), L^2(\mathbb{R}^d))$, $\sigma^{d/4} \nabla \partial_t \bar{u} \in L^{8/d}((0, T), L^4(\mathbb{R}^d))$;

(3) $\sqrt{\sigma} \partial_{tt} \bar{u} \in L^2((0, T) \times \mathbb{R}^d)$, $\sigma \partial_{tt} \bar{u} \in L^\infty((0, T), L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d))$.

In Section 3 below, we study the linearized equations associated with (2.8). The proof of the existence of $\bar{u}$ is the purpose of Section 4, and it relies on a fixed point argument. The proof of uniqueness follows a similar argument to that in [6] and is therefore omitted here. It only requires viscosity far from vacuum and the velocity gradient to meet the following conditions:

$$ \nabla \bar{u} \in L^1((0, T), L^\infty(\mathbb{R}^d)), \quad \text{and} \quad \sqrt{\sigma} \nabla \bar{u} \in L^2((0, T), L^\infty(\mathbb{R}^d)). $$

In particular, no small fluctuation assumption on the viscosity $\mu$ is required. Once we have constructed $\bar{u}$, we express $(\bar{c}, \bar{\rho})$ from (2.7)$_{1,2}$, and we finally go back into Eulerian coordinates to obtain that $(c, \rho, u)$ solves (1.8) with initial data $(c_0, \rho_0, u_0)$. We now proceed to the study of the linearized equations associated with (2.8).

### 3. STUDY OF THE LINEARISED SYSTEM

#### 3.1. Statement of the main result.
The aim of this section is to establish the well-posedness of the linear non-homogeneous parabolic equation

$$ (3.1) \qquad \rho_0 u_t - \text{div}(2\mu \mathbb{D}u + \lambda \text{div} u I_d) = \text{div} F $$

associated to (2.7). The initial density $\rho_0 = \rho_0(x)$ and the viscosity $\mu = \mu(t, x)$, $\lambda = \lambda(t, x)$, are supposed to verify:

$$ (3.2) \qquad \left\{ \begin{array}{l} 0 < \rho_* \le \rho_0 \le \rho^*, \quad 0 < \mu_* \le \mu \le \mu^*, \quad 0 \le \lambda \le \mu^*; \\ \mu, \lambda \in L^\infty((0,T), \mathcal{C}^\alpha_{pw,C_0}(\mathbb{R}^d)); \\ \mu_t, \lambda_t \in L^\infty((0,T), L^4(\mathbb{R}^d)) \cap L^1((0,T), L^\infty(\mathbb{R}^d)) \cap L^{16}((0,T), L^8(\mathbb{R}^d)); \\ \sqrt{\sigma}(\mu_t, \lambda_t) \in L^2((0,T), L^\infty(\mathbb{R}^d)), \sigma(\mu_{tt}, \lambda_{tt}) \in L^2((0,T), L^4(\mathbb{R}^d)), \end{array} \right. $$

for positive constants $\rho_*, \rho^*, \mu_*, \mu^*$ and time $0 < T < \infty$. The source term $F$ in (3.1) and the initial datum $u_0$ are such that $(F, u_0) \in Y_T(\mu_0, \lambda_0)$ where:

$$ (3.3) \qquad Y_T(\mu_0, \lambda_0) := \left\{ \begin{aligned} &(F, u_0) \in \mathcal{C}([0, T], L^2(\mathbb{R}^d)) \times H^1(\mathbb{R}^d) : F \in L^\infty((0, T), L^6(\mathbb{R}^d)), \\ &F \in L^{16}((0, T), L^8(\mathbb{R}^d)), \sigma^{r/4} F \in L^4((0, T), \mathcal{C}^\alpha_{pw,C_0}(\mathbb{R}^d)), \\ &F_t \in L^2((0, T) \times \mathbb{R}^d), \sigma^{d/4} F_t \in L^{8/d}((0, T), L^4(\mathbb{R}^d)), \\ &\sigma F_{tt} \in L^2((0, T) \times \mathbb{R}^d), \text{div}(2\mu_0 \mathbb{D}u_0 + \lambda_0 \text{div} u_0 I_d + F|_{t=0}) \in L^2(\mathbb{R}^d) \end{aligned} \right\}. $$

Above, $\mu_0 := \mu|_{t=0}$, $\lambda_0 := \lambda|_{t=0}$ and $0 < r < 8/3$ depends on $0 < \alpha < 1$ and the dimension, see (3.46) below. The existence and uniqueness of solution to (3.1) is established in the following space:

$$ (3.4) \qquad X_T := \left\{ \begin{array}{l} u \in \mathcal{C}([0, T], H^1(\mathbb{R}^d)) : \nabla u \in L^\infty((0, T), L^6(\mathbb{R}^d)) \cap L^{16}((0, T), L^8(\mathbb{R}^d)), \sigma^{r/4} \nabla u \in L^4((0, T), \mathcal{C}^\alpha_{pw,C_0}(\mathbb{R}^d)), \\ u_t \in \mathcal{C}([0, T], L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)), \sqrt{\sigma} \nabla u_t \in L^\infty((0, T), L^2(\mathbb{R}^d)), \sigma^{d/4} \nabla u_t \in L^{8/d}((0, T), L^4(\mathbb{R}^d)), \\ \sqrt{\sigma} u_{tt} \in L^2((0, T) \times \mathbb{R}^d), \sigma u_{tt} \in L^\infty((0, T), L^2(\mathbb{R}^d)) \cap L^2((0, T), \dot{H}^1(\mathbb{R}^d)) \end{array} \right\}. $$

The time weights are inspired by the pioneering work of Hoff [13]. One of the main difficulties is to obtain estimates for $\nabla u$, the solution of (3.1), that allow us to propagate the regularity of the interface $C_0$. The difficulty arises from the fact that the density and viscosity do not have any (global-in-space) Sobolev regularity, and therefore we are unable to differentiate (3.1) with respect to space.
and thus, we deduce that $\sqrt{\rho_0}u_t(\tau)$ convergences to $\sqrt{\rho_0}u_{t|t=0}$ strongly in $L^2(\mathbb{R}^d)$. From the fact that $\rho_0$ is bounded away from vacuum, we finally deduce the continuity of $u_t$ at the initial time. Moreover, this solution is unique, since, the difference $\delta u$ of two solutions $u_1$ and $u_2$ of (3.13) verifies:

$$ \begin{cases} \rho_0 \delta u_t - \tilde{\mu} \Delta \delta u - (\tilde{\mu} + \tilde{\lambda}) \nabla \delta u = 0, \\ \delta u|_{t=0} = 0 \end{cases} $$

and one has, from the classical energy computation, that $\delta u = 0$. Furthermore the gradient of $u$ is given by:

$$ \nabla u = -\frac{1}{\tilde{\mu}}(-\Delta)^{-1}\nabla [\mathcal{P}(\rho_0 u_t) - \mathcal{P}\operatorname{div} F] - \frac{1}{2\tilde{\mu} + \tilde{\lambda}}(-\Delta)^{-1}\nabla [\mathcal{Q}(\rho_0 u_t) - \mathcal{Q}\operatorname{div} F]. $$

where $\mathcal{P}$ is the Leray projector onto the space of divergence-free vector fields and $\mathcal{Q}$ is the projector onto the space of curl-free vector fields. Finally, using the regularity of $u_t$ and the Hölder regularity of $F$, one has the last estimate of $\nabla u$ that is (3.10). All these discussions establish that $\mathcal{E}$ is a non-empty set. The only thing left to prove is that $\mathcal{E}$ is both an open and a closed subset of $[0, 1]$. This objective is the purpose of the final section Section 3.3.3.

### 3.3. Proofs.

#### 3.3.1. A priori estimates.

*Proof.* This section consists in performing five kind of energy estimates on the solution $u \in X_T$ of the problem (3.5). We notice that the regularity of $u$ is enough in order to use $u$ or $u_t$ as test function in subsequent calculations. Let us begin by the first energy estimate:

**First step:** Estimates for $u$. This energy consists in taking the scalar product of the equation (3.5) with the velocity $u$, in order to obtain the following:

$$ \frac{1}{2} \frac{d}{dt} \|\sqrt{\rho_0} u\|_{L^2(\mathbb{R}^d)}^2 + \int_{\mathbb{R}^d} 2\mu |\mathbb{D}u|^2 + \int_{\mathbb{R}^d} \lambda |\operatorname{div} u|^2 = - \int_{\mathbb{R}^d} F^{jk} \partial_k u^j. $$

We then integrate the above in time in order to obtain the following:

$$ (3.17) \qquad \frac{1}{2} \|\sqrt{\rho_0} u(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \int_{\mathbb{R}^d} 2\mu |\mathbb{D}u|^2 + \int_0^t \int_{\mathbb{R}^d} \lambda |\operatorname{div} u|^2 = \frac{1}{2} \|\sqrt{\rho_0} u_0\|_{L^2(\mathbb{R}^d)}^2 - \int_0^t \int_{\mathbb{R}^d} F^{jk} \partial_k u^j $$

and finally Hölder and Young inequalities yield:

$$ \sup_{[0,T]} \|\sqrt{\rho_0} u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 \lesssim \|\sqrt{\rho_0} u_0\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|F\|_{L^2(\mathbb{R}^d)}^2. $$

**Second Step:** Estimates for $u_t$. We use $u_t$ as a test function in the weak formulation of (3.1) in order to obtain:

$$ \int_0^t \|\sqrt{\rho_0} u_t\|_{L^2(\mathbb{R}^d)}^2 + \int_{\mathbb{R}^d} \mu(t) |\mathbb{D}u(t)|^2 + \frac{1}{2} \int_{\mathbb{R}^d} \lambda(t) (\operatorname{div} u(t))^2 = \int_{\mathbb{R}^d} \mu_0 |\mathbb{D}u_0|^2 + \frac{1}{2} \int_{\mathbb{R}^d} \lambda_0 (\operatorname{div} u_0)^2 \\ - \int_{\mathbb{R}^d} F^{jk}(s) \partial_k u^j(s) \Big|_{s=0}^{s=t} + \int_0^t \int_{\mathbb{R}^d} F_t^{jk} \partial_k u^j + \int_0^t \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u|^2 + \frac{1}{2} \int_0^t \int_{\mathbb{R}^d} \lambda_t (\operatorname{div} u)^2 $$

In the following lines, we will estimate the right hand side in order to use Grönwall's lemma. The third and the fourth terms are:

$$ \left| \int_{\mathbb{R}^d} F^{jk}(s) \partial_k u^j(s) \right|_{s=0}^{s=t} \le \eta \|\nabla u(t)\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \|F(t)\|_{L^2(\mathbb{R}^d)}^2 + \|\nabla u(0)\|_{L^2(\mathbb{R}^d)}^2 + \|F(0)\|_{L^2(\mathbb{R}^d)}^2 $$

for some positive $\eta > 0$ small in order to absorb the norm of the gradient in the left hand side of the previous above equality. Next, the fifth term can be bounded as follows:

$$ \left| \int_0^t \int_{\mathbb{R}^d} F_t^{jk} \partial_k u^j \right| \le \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2. $$

Finally, the last two terms are controlled by:

$$ \int_0^t \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u|^2 + \frac{1}{2} \int_0^t \int_{\mathbb{R}^d} \lambda_t (\operatorname{div} u)^2 \lesssim \int_0^t (\|\mu_t\|_{L^\infty(\mathbb{R}^d)} + \|\lambda_t\|_{L^\infty(\mathbb{R}^d)}) \|\nabla u\|_{L^2(\mathbb{R}^d)}^2. $$

Gathering all of these estimates, and choosing $\eta$ small, one has:

$$ \int_0^t \|\sqrt{\rho_0} u_t\|_{L^2(\mathbb{R}^d)}^2 + \sup_{[0,t]} \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 \lesssim \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2 + \sup_{[0,t]} \|F\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} \|\nabla u\|_{L^2(\mathbb{R}^d)}^2. $$
The main result of this section is stated in the following theorem:

**Theorem 3.1.** Let $T \in (0, \infty)$, the viscosity $\mu$, $\lambda$ and the initial density $\rho_0$ verify (3.2), $0 < \alpha < 1$. There exists a constant $[\mu] > 0$ depending only on the dimension $d$, $\alpha$ and $\tilde{\mu}$ such that if:

$$ (3.5) \qquad \sup_{[0,T]} \left[ 1 + \|\lambda\|_{\dot{\mathcal{C}}_{\text{pw},C_0}^{\alpha}(\mathbb{R}^d)} + (\mathfrak{P}_{C_0} + \ell_{\varphi_0}^{-\alpha}) \|\llbracket \lambda \rrbracket\|_{L^{\infty}(C_0)} \right] \sup_{[0,T]} \|\mu - \tilde{\mu}\|_{\mathcal{C}_{\text{pw},C_0}^{\alpha}(\mathbb{R}^d)} \\ \qquad\qquad\qquad + (\mathfrak{P}_{C_0} + \ell_{\varphi_0}^{-\alpha}) \sup_{[0,T]} \left[ \|\llbracket \mu \rrbracket\|_{L^{\infty}(C_0)} + \|\llbracket \mu \rrbracket\|_{L^{\infty}(C_0)} \left\| \frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle} \right\|_{L^{\infty}(C_0)} \right] \le [\mu], $$

then, for all $(F, u_0) \in Y_T(\mu_0, \lambda_0)$, there exists a unique solution of the Cauchy problem:

$$ (3.6) \qquad \begin{cases} \rho_0 u_t - \operatorname{div}(\mu \mathbb{D} u) - \nabla(\lambda \operatorname{div} u) = \operatorname{div} F, \\ u|_{t=0} = u_0, \end{cases} $$

in $X_T$. Moreover, the following estimates hold true with a constant $C_*$ that depends polynomially on the upper and lower bounds of the density and viscosity (see first line of (3.2)):

* **Estimates for $u$ and $u_t$:**

$$ (3.7) \qquad \sup_{[0,T]} \|u, u_t, \nabla u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla u, \nabla u_t, u_t\|_{L^2(\mathbb{R}^d)}^2 \le C_* \left( \|(F, u_0)\|_{Y_T(\mu_0, \lambda_0)}^2 + \int_0^T \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F\|_{L^2 \cap L^6(\mathbb{R}^d)}^2 \right) \\ \qquad\qquad\qquad \times \exp \left( C_* \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \right] \right). $$

* **Higher integrability estimates for $\nabla u$:**

$$ (3.8) \qquad \begin{cases} \|\nabla u\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \le C_* \|F\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} + C_* \sup_{[0,T]} \|F, u_t, \nabla u\|_{L^2(\mathbb{R}^d)}, \\ \|\nabla u\|_{L^{16}((0,T),L^8(\mathbb{R}^d))} \le C_* (\|F\|_{L^{16}((0,T),L^8(\mathbb{R}^d))} + \|F, u_t, \nabla u\|_{L^2 \cap L^\infty((0,T),L^2(\mathbb{R}^d))} + \|\nabla u_t\|_{L^2((0,T) \times \mathbb{R}^d)}). \end{cases} $$

* **Estimates for $u_{tt}$:**

$$ (3.9) \qquad \int_0^T \|\sqrt{\sigma} u_{tt}, \sigma \nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \|\sqrt{\sigma} \nabla u_t, \sigma u_{tt}\|_{L^\infty((0,T),L^2(\mathbb{R}^d))}^2 \le C_* \left[ \int_0^{\sigma(T)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|F_t, \sigma F_{tt}\|_{L^2(\mathbb{R}^d)}^2 \right. \\ \qquad\qquad\qquad \left. + \int_0^T \|\mu_t, \lambda_t, \sigma(\mu_{tt}, \lambda_{tt})\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \right] \\ \qquad\qquad\qquad \times \exp \left( C_* \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \sigma \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)}^2 \right] \right). $$

The norms of $u$ involved in this estimate are controlled thanks to (3.7)-(3.8).

* $L^{8/d}((0,T), L^4(\mathbb{R}^d))$ norm estimate for $\sigma^{\frac{d}{4}}\nabla u_t$:

$$ (3.10) \qquad \int_0^T \sigma^2 \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \le C_* \int_0^T \sigma^2 \|F_t\|_{L^4(\mathbb{R}^d)}^{8/d} + C_* \int_0^T \sigma^2 \|\nabla u\|_{L^8(\mathbb{R}^d)}^{8/d} \|\mu_t, \lambda_t\|_{L^8(\mathbb{R}^d)}^{8/d} \\ \qquad\qquad\qquad + C_* \|\sigma u_{tt}\|_{L^\infty((0,T),L^2(\mathbb{R}^d))}^2 \int_0^T \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathbb{D} u\|_{L^2(\mathbb{R}^d)}^{\frac{8}{d}-2} $$

* Hölder estimate for $\nabla u$:

$$ (3.11) \qquad \int_0^T \sigma^r \|\nabla u\|_{\mathcal{C}_{\text{pw},C_0}^{\alpha}(\mathbb{R}^d)}^4 \le C_{**} \int_0^T \sigma^r \|F\|_{\mathcal{C}_{\text{pw},C_0}^{\alpha}(\mathbb{R}^d)}^4 + C_{**} \|F, \nabla u, u_t, \sqrt{\sigma} \nabla u_t\|_{L^{\infty}((0,T),L^2(\mathbb{R}^d))}^4 \\ \qquad\qquad\qquad + C_{**} \left[ \int_0^T \|F, \nabla u, u_t, \nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \right]^2 + C_{**} \left[ \int_0^T \sigma^2 \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \right]^{d/2}, $$

with a constant $C_{**}$ that depends polynomially on the same parameters as $C_*$, and on $\alpha$, $\ell_{\varphi_0}^{-\frac{\alpha}{2}}$, $\mathfrak{P}_{C_0}$ and $\|\mu, \lambda\|_{\mathcal{C}_{\text{pw},C_0}^{\alpha}(\mathbb{R}^d)}$.
We turn to the computations of the second energy that consists in deriving the equation (3.5) in time and use $u_t$ as test function in the weak formulation of the obtained equation. By doing so, we obtain:

$$ \begin{aligned} & \frac{1}{2} \int_{\mathbb{R}^d} \rho_0 |u_t(t)|^2 + \int_0^t 2 \int_{\mathbb{R}^d} \mu |\mathbb{D}u_t|^2_{L^2(\mathbb{R}^d)} + \int_0^t \int_{\mathbb{R}^d} \lambda(\mathrm{div} \, u_t)^2 = \frac{1}{2} \int_{\mathbb{R}^d} \rho_0 |u_t|_{t=0}|^2 \\ & \qquad - \int_0^t \int_{\mathbb{R}^d} 2\mu_t \mathbb{D}^j u \partial_k u_t^j - \int_0^t \int_{\mathbb{R}^d} \lambda_t \mathrm{div} \, u \mathrm{div} \, u_t + \int_0^t \int_{\mathbb{R}^d} F_t^{jk} \partial_k u_t^j \\ & \le \int_{\mathbb{R}^d} \rho_0 |u_t|_{t=0}|^2 + \eta \int_0^t \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 \\ & \qquad + C_\eta \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \end{aligned} \quad (3.18) $$

with some positive $\eta$ small related to the lower and upper bounds of the viscosity. We now turn to the estimate of the $L^4(\mathbb{R}^d)$ norm of $\nabla u$ that appears in the above inequality. We begin by observing that from (3.5), one has:

$$ -\tilde{\mu}\Delta u - \nabla\{(\tilde{\mu} + \lambda)\mathrm{div}\,u\} = \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} - \rho_0 u_t + \mathrm{div}\,F $$

and so, we can express:

$$ \begin{aligned} & \nabla u = \nabla P u + \nabla (I - P) u \\ & = \frac{1}{\tilde{\mu}}(-\Delta)^{-1} \nabla [P \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} - P(\rho_0 \partial_t u) + P \mathrm{div}\,F] \\ & + \frac{1}{2\tilde{\mu} + \lambda}(-\Delta)^{-1} \nabla [Q \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} - Q(\rho_0 \partial_t u) + Q(\mathrm{div}\,F)]. \end{aligned} \quad (3.19) $$

Above, $P$ is the Leray projector onto the space of divergence-free vector fields and where $Q$ is the projector onto the space of curl-free vector fields. Since $F \in L^4(\mathbb{R}^d)$ a.e in time, singular integral theory helps us to obtain:

$$ \begin{aligned} \|\nabla u\|_{L^4(\mathbb{R}^d)} \lesssim_d \frac{1}{\tilde{\mu}} \|\mu - \tilde{\mu}\|_{L^\infty(\mathbb{R}^d)} \|\nabla u\|_{L^4(\mathbb{R}^d)} + \|F\|_{L^4(\mathbb{R}^d)} \\ + \frac{1}{\tilde{\mu}} \|(-\Delta)^{-1} \nabla P(\rho_0 u_t), (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{L^4(\mathbb{R}^d)}. \end{aligned} \quad (3.20) $$

Therefore, for small perturbation of the viscosity $\mu$ around the constant state $\tilde{\mu}$, one can absorb the first term of the right hand side in the left one. By Gagliardo-Nirenberg inequality, the last terms above can be bounded by:

$$ \|(-\Delta)^{-1} \nabla P(\rho_0 u_t), (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{L^4(\mathbb{R}^d)} \lesssim \|\rho_0 u_t\|_{L^2(\mathbb{R}^d)}^{d/4} \|\nabla u, F\|_{L^2(\mathbb{R}^d)}^{1-d/4}. $$

In consequence, the last term of (3.18) can be estimated as:

$$ (3.21) \quad \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \lesssim \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F\|_{L^4(\mathbb{R}^d)}^2 + \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F, \nabla u, u_t\|_{L^2(\mathbb{R}^d)}^2. $$

We now add (3.21) to (3.18) and choose $\eta$ related to the lower bound and or upper bounds of the viscosity and density in order to absorb the second term of the right hand side of (3.18) in the left hand side and then obtains, after Gronwall's Lemma:

$$ (3.22) \quad \sup_{[0,T]} \|u_t, \nabla u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla u_t, u_t\|_{L^2(\mathbb{R}^d)}^2 \le M \left( \| (F, u_0) \|_{Y_T(\mu_0, \lambda_0)}^2 + \int_0^T \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F\|_{L^2 \cap L^6(\mathbb{R}^d)}^2 \right) \\ \times \exp \left( M \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \right] \right). $$

We turn to the estimate for the second time derivative of $u$.

Third step: Estimates for $u_{tt}$. They consist in two kinds of estimates and they are obtained from the equation satisfied by $u_t$. In fact, when we derive equations (3.5) in time, we have:

$$ (3.23) \quad \rho_0 u_{tt} - \mathrm{div}(2\mu \mathbb{D}u_t + \lambda \mathrm{div} \, u_t I_d) = \mathrm{div} \, F_t + \mathrm{div}(2\mu_t \mathbb{D}u + \lambda_t \mathrm{div} \, u_t I_d). $$
**Remark 3.2.** The first estimate (3.7) follows from the classical energy estimate and the so-called Hoff estimates introduced in [13, 14]. These estimates provide $H^1$ regularity for the time derivative of the velocity $u$. The other estimate (3.9) can be seen as a higher-order Hoff estimate and has be used in [7] in the study of the incompressible Navier-Stokes model.

The proof of **Theorem 3.1** is the subject of **Section 3.2** below and it is based on the classical method of a priori estimates combined with a homotopy argument.

**3.2. Proof of Theorem 3.1.** The proof of **Theorem 3.1** relies on a priori estimates and a homotopy argument and is structured into three sections: in **Section 3.2.1**, we derive a priori estimates (3.7), (3.8), (3.9), (3.10) and (3.11), whereas **Section 3.2.2-Section 3.2.3** focus on the homotopy argument.

**3.2.1. A priori estimates.** The aim of this section is to derive the a priori estimates (3.7), (3.8), (3.9), (3.10) and (3.11). We observe that the regularity of $u$ is sufficient to use $u$ or $u_t$ as test functions in subsequent calculations. We use $C$ to denote any constant that depends only on the dimension, $\tilde{\mu}$, and $\alpha$. Conversely, $C_*$ denotes a constant that may additionally depend polynomially on the lower and upper bounds of $\rho_0$, $\mu$ and $\lambda$.

**Proof of (3.7).** Here we derive energy estimates for $u$ and $u_t$. We take the scalar product of the equation (3.6) with the velocity $u$. By doing so, we obtain:

$$ \frac{1}{2} \frac{d}{dt} \|\sqrt{\rho_0} u\|_{L^2(\mathbb{R}^d)}^2 + \int_{\mathbb{R}^d} 2\mu |\mathbb{D}u|^2 + \int_{\mathbb{R}^d} \lambda |\mathrm{div} \, u|^2 = - \int_{\mathbb{R}^d} F^{jk} \partial_k u^j. $$

We then integrate the above in time and we have:

$$ (3.12) \quad \frac{1}{2} \|\sqrt{\rho_0} u(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \int_{\mathbb{R}^d} 2\mu |\mathbb{D}u|^2 + \int_0^t \int_{\mathbb{R}^d} \lambda |\mathrm{div} \, u|^2 = \frac{1}{2} \|\sqrt{\rho_0} u_0\|_{L^2(\mathbb{R}^d)}^2 - \int_0^t \int_{\mathbb{R}^d} F^{jk} \partial_k u^j. $$

Using the identity:

$$ (3.13) \quad \|\mathbb{D}w\|_{L^2(\mathbb{R}^d)}^2 = \frac{1}{2} \|\nabla w\|_{L^2(\mathbb{R}^d)}^2 + \frac{1}{2} \|\mathrm{div} \, w\|_{L^2(\mathbb{R}^d)}^2 $$

we deduce that

$$ (3.14) \quad \mu_* \left( \|\nabla w\|_{L^2(\mathbb{R}^d)}^2 + \|\mathrm{div} \, w\|_{L^2(\mathbb{R}^d)}^2 \right) \le \int_{\mathbb{R}^d} 2\mu(x) |\mathbb{D}w(x)|^2 dx. $$

This inequality will be freely used in the forthcoming computations to translate the $L^2(\mathbb{R}^d)$-norm estimate of $\sqrt{\mu}\mathbb{D}w$ into that of $\nabla w$ with $w \in \{u, u_t, u_{tt}\}$. For instance, by applying (3.14) with $w = u$ and using Hölder's and Young's inequalities, we obtain from (3.12):

$$ (3.15) \quad \sup_{[0,T]} \|\sqrt{\rho_0} u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 \le C_* \|\sqrt{\rho_0} u_0\|_{L^2(\mathbb{R}^d)}^2 + C_* \int_0^T \|F\|_{L^2(\mathbb{R}^d)}^2. $$

Next, we use $u_t$ as a test function in the weak formulation of (3.1) in order to obtain:

$$ \int_0^t \|\sqrt{\rho_0} u_t\|_{L^2(\mathbb{R}^d)}^2 + \int_{\mathbb{R}^d} \mu(t) |\mathbb{D}u(t)|^2 + \frac{1}{2} \int_{\mathbb{R}^d} \lambda(t) (\mathrm{div} \, u(t))^2 = \int_{\mathbb{R}^d} \mu_0 |\mathbb{D}u_0|^2 + \frac{1}{2} \int_{\mathbb{R}^d} \lambda_0 (\mathrm{div} \, u_0)^2 \\ - \int_{\mathbb{R}^d} F^{jk}(s) \partial_k u^j(s) \Big|_{s=0}^{s=t} + \int_0^t \int_{\mathbb{R}^d} F_t^{jk} \partial_k u^j + \int_0^t \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u|^2 + \frac{1}{2} \int_0^t \int_{\mathbb{R}^d} \lambda_t (\mathrm{div} \, u)^2. $$

In the following lines, we will estimate the right hand side in order to use Grönwall's lemma. The third and the fourth terms are:

$$ \left| \int_{\mathbb{R}^d} F^{jk}(s) \partial_k u^j(s) \right|_{s=0}^{s=t} \le \eta \|\nabla u(t)\|_{L^2(\mathbb{R}^d)}^2 + \frac{1}{4\eta} \|F(t)\|_{L^2(\mathbb{R}^d)}^2 + \|\nabla u(0)\|_{L^2(\mathbb{R}^d)}^2 + \|F(0)\|_{L^2(\mathbb{R}^d)}^2, $$

for some positive $\eta > 0$ small in order to absorb the norm of the gradient in the left hand side of the above equality. Next, the fifth term can be bounded as follows:

$$ \left| \int_0^t \int_{\mathbb{R}^d} F_t^{jk} \partial_k u^j \right| \le \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2. $$

Finally, the last two terms are controlled by:

$$ \int_0^t \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u|^2 + \frac{1}{2} \int_0^t \int_{\mathbb{R}^d} \lambda_t (\mathrm{div} \, u)^2 \le C_* \int_0^t (\|\mu_t\|_{L^\infty(\mathbb{R}^d)} + \|\lambda_t\|_{L^\infty(\mathbb{R}^d)}) \|\nabla u\|_{L^2(\mathbb{R}^d)}^2. $$
The first energy consists in multiplying the above equation by $\sigma u_{tt}$ whereas, the second one consists in deriving the above equation with respect to time before testing with $\sigma^2 u_{tt}$. Testing the equation (3.23) with $\sigma u_{tt}$ leads to:

$$ \frac{1}{2} \int_0^t \sigma \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \sigma(t) \int_{\mathbb{R}^d} \mu |\mathbb{D}u_t|^2 + \frac{\sigma(t)}{2} \int_{\mathbb{R}^d} \lambda |\text{div } u_t|^2 = \frac{1}{2} \int_0^{\sigma(t)} \int_{\mathbb{R}^d} \mu |\mathbb{D}u_t|^2 + \frac{1}{2} \int_0^{\sigma(t)} \int_{\mathbb{R}^d} \lambda |\text{div } u_t|^2 \\ + \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u_t|^2 + \frac{1}{2} \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t |\text{div } u_t|^2 - 2 \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t \mathbb{D}^{jk} u \partial_k u_{tt}^j \\ - \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t \text{div } u \text{div } u_{tt} - \int_0^t \sigma \int_{\mathbb{R}^d} F_t^{jk} \partial_k u_{tt}^j. \tag{3.24} $$

By Hölder's inequality, the last term can be estimated as:

$$ \left| \int_0^t \sigma \int_{\mathbb{R}^d} F_t^{jk} \partial_k u_{tt}^j \right| \lesssim \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 $$

and on the other hand, the remainder terms are bounded by:

$$ \left| \frac{1}{2} \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t |\text{div } u_t|^2 - \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t \text{div } u \text{div } u_{tt} + \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u_t|^2 - 2 \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t \mathbb{D}^{jk} u \partial_k u_{tt}^j \right| \\ \lesssim \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \\ + \int_0^t \sigma \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} $$

Gathering all of these estimates, we have:

$$ \int_0^t \sigma \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \sigma(t) \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \lesssim \int_0^{\sigma(t)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 \\ + C_\eta \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 + \int_0^t \sigma \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2. \tag{3.25} $$

The second term on the right hand side above will be reduced with the terms coming from the second energy which we turn to, in the following lines.

We first derive (3.23) with respect to time before multiplying the obtained equation by $\sigma^2 u_{tt}$. By doing so, one obtains:

$$ \frac{\sigma^2}{2} \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + 2 \int_0^t \sigma^2 \|\sqrt{\mu} \mathbb{D} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma^2 \|\sqrt{\lambda} \text{div } u_{tt}\|_{L^2(\mathbb{R}^d)}^2 = \int_0^{\sigma(t)} \sigma \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 \\ - 4 \int_0^t \sigma^2 \mu_t \mathbb{D}^{jk} u_t \partial_k u_{tt}^j - 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_t \text{div } u_t \text{div } u_{tt} - \int_0^t \sigma^2 \int_{\mathbb{R}^d} F_{tt}^{jk} \partial_k u_{tt}^j \\ - 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \mu_{tt} \mathbb{D}^{jk} u \partial_k u_{tt}^j - \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_{tt} \text{div } u \text{div } u_{tt}. $$

We estimate the second and the third terms of the right hand side using Young's inequality, as follows:

$$ \left| 4 \int_0^t \sigma^2 \mu_t \mathbb{D}^{jk} u_t \partial_k u_{tt}^j + 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_t \text{div } u_t \text{div } u_{tt} \right| \lesssim \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t \sigma^2 \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2, $$
whereas the last two terms are bounded by

$$ \left| 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \mu_{tt} \mathbb{D}^{jk} u \partial_k u_{tt}^j + \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_{tt} \text{div } u \text{div } u_{tt} \right| \lesssim \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t \sigma^2 \|\mu_{tt}, \lambda_{tt}\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2. $$

Finally, the fourth term of the right hand side is:

$$ \left| \int_0^t \sigma^2 \int_{\mathbb{R}^d} F_{tt}^{jk} \partial_k u_{tt}^j \right| \lesssim \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t \sigma^2 \|F_{tt}\|_{L^2(\mathbb{R}^d)}^2. $$
Gathering all of these estimates, and choosing $\eta$ small, one has:

$$ (3.16) \qquad \int_0^t \|\sqrt{\rho_0} u_t\|_{L^2(\mathbb{R}^d)}^2 + \sup_{[0,t]} \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 \le C_* \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2 + C_* \sup_{[0,t]} \|F\|_{L^2(\mathbb{R}^d)}^2 + C_* \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 \\ \qquad \qquad + C_* \int_0^t \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} \|\nabla u\|_{L^2(\mathbb{R}^d)}^2. $$

We differentiate (3.6) in time and we use $u_t$ as a test function in the weak formulation of the resulting equation. By doing so, we obtain:

$$ (3.17) \qquad \frac{1}{2} \int_{\mathbb{R}^d} \rho_0 |u_t(t)|^2 + \int_0^t \int_{\mathbb{R}^d} 2\mu |\mathbb{D}u_t|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \int_{\mathbb{R}^d} \lambda(\mathrm{div} \, u_t)^2 = \frac{1}{2} \int_{\mathbb{R}^d} \rho_0 |u_t|_{t=0}^2 \\ \qquad - \int_0^t \int_{\mathbb{R}^d} 2\mu_t \mathbb{D}^j k u \partial_k u_t^j - \int_0^t \int_{\mathbb{R}^d} \lambda_t \mathrm{div} \, u \mathrm{div} \, u_t + \int_0^t \int_{\mathbb{R}^d} F_t^{jk} \partial_k u_t^j \\ \qquad \le \frac{1}{2} \int_{\mathbb{R}^d} \rho_0 |u_t|_{t=0}^2 + \eta \int_0^t \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + \frac{C_*}{\eta} \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 \\ \qquad + \frac{C_*}{\eta} \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2, $$

with some positive $\eta$ small related to the lower bound of the viscosity $\mu$. We now turn to the estimate of the $L^4(\mathbb{R}^d)$ norm of $\nabla u$ that appears in the above inequality. We rewrite (3.6) as:

$$ (3.18) \qquad -\tilde{\mu}\Delta u - \nabla\{(\tilde{\mu} + \lambda)\mathrm{div} \, u\} = \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} - \rho_0 u_t + \mathrm{div} \, F. $$

We denote by $\mathcal{P} = I - \mathcal{Q}$ and $\mathcal{Q} = -\nabla(-\Delta)^{-1} \mathrm{div}$ the usual projectors onto the spaces of divergence-free and curl-free vector fields, respectively. First, we apply $\mathcal{P}$ to (3.18) and we have:

$$ (3.19) \qquad \tilde{\mu}\mathcal{P}u = (-\Delta)^{-1}\mathcal{P}\mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u + F\} - (-\Delta)^{-1}\mathcal{P}(\rho_0 u_t). $$

Second, we apply the divergence operator to obtain:

$$ (3.20) \qquad (2\tilde{\mu} + \lambda) \mathrm{div} \, u = (-\Delta)^{-1} \mathrm{div} \, \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u + F\} - (-\Delta)^{-1} \mathrm{div}(\rho_0 u_t) $$

and given that $\mathcal{Q}w = -\nabla(-\Delta)^{-1} \mathrm{div} \, w$, we find:

$$ \mathcal{Q}u = -\nabla(-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \mathrm{div} \, \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u + F\} \right] + \nabla(-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \mathrm{div}(\rho_0 u_t) \right] $$

Summing up, we express:

$$ \begin{aligned} \nabla u &= \nabla \mathcal{P}u + \nabla \mathcal{Q}u \\ &= \frac{1}{\tilde{\mu}} (-\Delta)^{-1} \nabla [\mathcal{P} \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u + F\} - \mathcal{P}(\rho_0 u_t)] + \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \mathrm{div}(\rho_0 u_t) \right] \\ (3.21) \qquad & \quad - \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \mathrm{div} \, \mathrm{div}\{2(\mu - \tilde{\mu})\mathbb{D}u + F\} \right]. \end{aligned} $$

Furthermore, from (3.18), we have:

$$ (-\Delta)^{-1} \nabla(\rho_0 u_t) = \nabla(-\Delta)^{-1} \mathrm{div}(2\mu\mathbb{D}u) + \nabla(-\Delta)^{-1} \nabla(\lambda \mathrm{div} \, u) + \nabla(-\Delta)^{-1} \mathrm{div} \, F $$

and owing to the fact that $\nabla(-\Delta)^{-1} \nabla$ is continuous over $L^p(\mathbb{R}^d)$, $p \in (1, \infty)$, we have:

$$ (3.22) \qquad \|(-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{L^2(\mathbb{R}^d)} \le C_* \|\nabla u, F\|_{L^2(\mathbb{R}^d)}. $$

Additionally, singular integral theory helps us to obtain from (3.21):

$$ (3.23) \qquad \|\nabla u\|_{L^4(\mathbb{R}^d)} \le \frac{C}{\tilde{\mu}} \|\mu - \tilde{\mu}\|_{L^\infty(\mathbb{R}^d)} \|\nabla u\|_{L^4(\mathbb{R}^d)} + C \|F\|_{L^4(\mathbb{R}^d)} + C \|(-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{L^4(\mathbb{R}^d)}. $$

Therefore, under the hypothesis (3.5), we absorb the first term of the right hand side in the left one. With the help of Gagliardo-Nirenberg's inequality and (3.22), the last terms above can be bounded as (recall $d \in \{2, 3\}$):

$$ \|(-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{L^4(\mathbb{R}^d)} \le C_* \|\rho_0 u_t\|_{L^2(\mathbb{R}^d)}^{d/4} \|\nabla u, F\|_{L^2(\mathbb{R}^d)}^{1-d/4}. $$
Gathering all of these estimates, and adding (3.25) before choosing $\eta$ small related to the lower bounds of the viscosity, one finally applies Grönwall's lemma in order to obtain the following:

$$ (3.26) \qquad \int_0^T \|\sigma u_{tt}, \sigma^2 \nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \|\sigma \nabla u_t, \sigma^2 u_{tt}\|_{L^\infty((0,T),L^2(\mathbb{R}^d))}^2 \le M \left[ \int_0^{\sigma(T)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|F_t, \sigma^2 F_{tt}\|_{L^2(\mathbb{R}^d)}^2 \right. \\ \qquad \qquad \qquad \qquad \left. + \int_0^T \|\mu_t, \lambda_t, \sigma(\mu_{tt}, \lambda_{tt})\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \right] \\ \qquad \qquad \qquad \qquad \times \exp \left( M \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \|\mu_{tt}, \lambda_{tt}\|_{L^2(\mathbb{R}^d)} \right] + M \int_0^T \sigma \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)}^2 \right). $$

This ends the computations of the last energy.

Fourth Step: $L^p((0,T), L^q(\mathbb{R}^d))$ norm estimate for $\nabla u$. Starting from (3.19), one has:

$$ \|\nabla u\|_{L^q(\mathbb{R}^d)} \lesssim_{d,q} \frac{1}{\mu} \|\mu - \tilde{\mu}\|_{L^\infty(\mathbb{R}^d)} \|\nabla u\|_{L^q(\mathbb{R}^d)} + \|F\|_{L^q(\mathbb{R}^d)} + \frac{1}{\mu} \|(-\Delta)^{-1} \nabla P(\rho_0 u_t), (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{L^q(\mathbb{R}^d)}. $$

Therefore, for small perturbation of the viscosity $\mu$ around the constant state $\tilde{\mu}$, one can absorb the first term of the right hand side in the left one. We only have to estimate the last terms, by Gagliardo-Nirenberg inequality:

$$ \|(-\Delta)^{-1} \nabla P(\rho_0 u_t), (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{L^6(\mathbb{R}^d)} \lesssim \|\rho_0 u_t\|_{L^2(\mathbb{R}^d)}^{\beta} \|\nabla u, F\|_{L^2(\mathbb{R}^d)}^{1-\beta} $$

with $\beta = d/2 - d/6 \in \{2/3, 1\}$. In consequence, we have:

$$ (3.27) \qquad \|\nabla u\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \lesssim \|F\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} + \sup_{[0,T]} \|F, u_t, \nabla u\|_{L^2(\mathbb{R}^d)}. $$

On the other hand, for $q = 8$, in two space dimension one has:

$$ \|(-\Delta)^{-1} \nabla P(\rho_0 u_t), (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{L^8(\mathbb{R}^2)} \lesssim \|\rho_0 u_t\|_{L^2(\mathbb{R}^2)}^{3/4} \|\nabla u, F\|_{L^2(\mathbb{R}^2)}^{1/4} $$

whereas in three space dimension, Sobolev and interpolation inequalities yield:

$$ \|(-\Delta)^{-1} \nabla P(\rho_0 u_t), (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{L^8(\mathbb{R}^3)} \lesssim \|\rho_0 u_t\|_{L^{24/11}(\mathbb{R}^3)} \lesssim \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{1/8} \|u_t\|_{L^2(\mathbb{R}^3)}^{7/8}, $$

such that,

$$ (3.28) \qquad \|\nabla u\|_{L^{16}((0,T),L^8(\mathbb{R}^d))} \lesssim \|F\|_{L^{16}((0,T),L^8(\mathbb{R}^d))} + \|F, u_t, \nabla u\|_{L^2 \cap L^\infty((0,T),L^2(\mathbb{R}^d))} + \|\nabla u_t\|_{L^2((0,T) \times \mathbb{R}^d)}. $$

Fifth Step: $L^4((0,T) \times \mathbb{R}^d)$ norm estimate for $\nabla u_t$. As above, we begin by recalling that, from (3.23), one has:

$$ \begin{aligned} \nabla u_t &= \nabla P u_t + \nabla Q u_t \\ &= \frac{1}{\mu} (-\Delta)^{-1} \nabla [P \operatorname{div}\{2(\mu - \tilde{\mu})Du_t\} - P(\rho_0 u_{tt}) + P \operatorname{div} F_t + P \operatorname{div}(2\mu_t Du + \lambda_t \operatorname{div} uI_d)] \\ &\quad + \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \nabla [Q \operatorname{div}\{2(\mu - \tilde{\mu})Du_t\} - Q(\rho_0 u_{tt}) + Q \operatorname{div} F_t + Q \operatorname{div}(2\mu_t Du + \lambda_t \operatorname{div} uI_d)] \end{aligned} $$

so it holds that, under the smallness of the viscosity, one has:

$$ \|\nabla u_t\|_{L^4(\mathbb{R}^d)} \lesssim \|F_t\|_{L^4(\mathbb{R}^d)} + \|\nabla u\|_{L^8(\mathbb{R}^d)} \|\mu_t, \lambda_t\|_{L^8(\mathbb{R}^d)} + \|(-\Delta)^{-1} \nabla P(\rho_0 u_{tt}), (-\Delta)^{-1} \nabla Q(\rho_0 u_{tt})\|_{L^4(\mathbb{R}^d)} $$

and we only need to compute

$$ \int_0^T \sigma^{4/d} \|(-\Delta)^{-1} \nabla P(\rho_0 u_{tt}), (-\Delta)^{-1} \nabla Q(\rho_0 u_{tt})\|_{L^4(\mathbb{R}^d)}^{8/d}. $$

By Gagliardo-Nirenberg inequality, one has:

$$ \int_0^T \sigma^{4/d} \|(-\Delta)^{-1} \nabla P(\rho_0 u_{tt}), (-\Delta)^{-1} \nabla Q(\rho_0 u_{tt})\|_{L^4(\mathbb{R}^d)}^{8/d} \lesssim \int_0^T \sigma^{\frac{4}{d}} \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t Du\|_{L^2(\mathbb{R}^d)}^{\frac{8}{d}-2} \|u_{tt}\|_{L^2(\mathbb{R}^d)}^2 \|u_t\|_{L^2(\mathbb{R}^d)}^2. $$

When $d = 2$, one has:

$$ \int_0^T \sigma^2 \|(-\Delta)^{-1} \nabla P(\rho_0 u_{tt}), (-\Delta)^{-1} \nabla Q(\rho_0 u_{tt})\|_{L^4(\mathbb{R}^2)}^4 \lesssim \|\sigma u_{tt}\|_{L^\infty((0,T),L^2(\mathbb{R}^2))}^2 \int_0^T \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t Du\|_{L^2(\mathbb{R}^2)}^2 \|u_t\|_{L^2(\mathbb{R}^2)}^2. $$
In consequence, the last term of (3.17) can be estimated as:

$$ (3.24) \qquad \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \le C_* \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F\|_{L^4(\mathbb{R}^d)}^2 \\ \qquad \qquad + C_* \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F, \nabla u, u_t\|_{L^2(\mathbb{R}^d)}^2. $$

We now use (3.24) in (3.17), choose $\eta$ related to the lower bound of the viscosity and density in order to absorb the second term of the right hand side of (3.17) in the left hand side. Applying Gronwall's Lemma and adding (3.16) we find:

$$ (3.25) \quad \sup_{[0,T]} \|u_t, \nabla u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla u_t, u_t\|_{L^2(\mathbb{R}^d)}^2 \le C_* \left( \| (F, u_0)\|_{Y_T(\mu_0, \lambda_0)}^2 + \int_0^T \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|F\|_{L^2 \cap L^6(\mathbb{R}^d)}^2 \right) \\ \qquad \times \exp \left( C_* \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \right] \right). $$

Finally, (3.7) follows from combining (3.25) and (3.15). We now turn to the proof of (3.8). □

**Proof of (3.8).** Here we derive $L^q(\mathbb{R}^d)$-norm estimate for the velocity gradient. According to (3.21), for $q \in \{6, 8\}$ we have:

$$ \|\nabla u\|_{L^q(\mathbb{R}^d)} \le \frac{C}{\mu} \|\mu - \tilde{\mu}\|_{L^\infty(\mathbb{R}^d)} \|\nabla u\|_{L^q(\mathbb{R}^d)} + C \|F\|_{L^q(\mathbb{R}^d)} + C \|(-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{L^q(\mathbb{R}^d)}. $$

Therefore, under the hypothesis (3.5), we absorb the first term of the right hand side in the left one. We only have to estimate the last term. Using Gagliardo-Nirenberg's inequality and (3.22) we obtain:

$$ \|(-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{L^6(\mathbb{R}^d)} \le C_* \|\rho_0 u_t\|_{L^2(\mathbb{R}^d)}^{\beta} \|\nabla u, F\|_{L^2(\mathbb{R}^d)}^{1-\beta} $$

with $\beta = d/2 - d/6 \in \{2/3, 1\}$. In consequence, we have:

$$ (3.26) \qquad \|\nabla u\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \le C_* \|F\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} + C_* \sup_{[0,T]} \|F, u_t, \nabla u\|_{L^2(\mathbb{R}^d)}. $$

For $q = 8$ and $d = 2$, we have:

$$ \|(-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{L^s(\mathbb{R}^2)} \le C_* \|\rho_0 u_t\|_{L^2(\mathbb{R}^2)}^{3/4} \|\nabla u, F\|_{L^2(\mathbb{R}^2)}^{1/4} $$

whereas for $d = 3$, Sobolev's and interpolation inequalities yield:

$$ \|(-\Delta)^{-1} \nabla P(\rho_0 u_t)\|_{L^s(\mathbb{R}^3)} \le C_* \|\rho_0 u_t\|_{L^{24/11}(\mathbb{R}^3)} \le C_* \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{1/8} \|u_t\|_{L^2(\mathbb{R}^3)}^{7/8}, $$

such that,

$$ (3.27) \qquad \|\nabla u\|_{L^{16}((0,T),L^s(\mathbb{R}^d))} \le C_* \|F\|_{L^{16}((0,T),L^s(\mathbb{R}^d))} + C_* \|F, u_t, \nabla u\|_{L^2 \cap L^\infty((0,T),L^2(\mathbb{R}^d))} + C_* \|\nabla u_t\|_{L^2((0,T) \times \mathbb{R}^d)}. $$

This concludes (3.8) and we turn to the proof of (3.9). □

**Proof of (3.9).** Here we establish estimates for $u_{tt}$. We first differentiate equation (3.6) in time:

$$ (3.28) \qquad \rho_0 u_{tt} - \operatorname{div}(2\mu \mathbb{D}u_t + \lambda \operatorname{div} u_t I_d) = \operatorname{div} F_t + \operatorname{div}(2\mu_t \mathbb{D}u + \lambda_t \operatorname{div} u I_d). $$

Testing (3.28) with $\sigma u_{tt}$ (recall $\sigma(t) = \min\{1, t\}$ and hence $\sigma(0) = 0$) leads to:

$$ \begin{aligned} & \frac{1}{2} \int_0^t \sigma \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \sigma(t) \int_{\mathbb{R}^d} \mu |\mathbb{D}u_t|^2 + \frac{\sigma(t)}{2} \int_{\mathbb{R}^d} \lambda |\operatorname{div} u_t|^2 = \int_0^t \sigma'(t) \int_{\mathbb{R}^d} \mu |\mathbb{D}u_t|^2 + \frac{1}{2} \int_0^t \sigma'(t) \int_{\mathbb{R}^d} \lambda |\operatorname{div} u_t|^2 \\ & \qquad + \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u_t|^2 + \frac{1}{2} \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t |\operatorname{div} u_t|^2 - 2 \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t \mathbb{D}^k u \partial_k u_{tt}^j \\ & \qquad - \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t \operatorname{div} u \operatorname{div} u_{tt} - \int_0^t \sigma \int_{\mathbb{R}^d} F_t^{jk} \partial_k u_{tt}^j. \end{aligned} \quad (3.29) $$

Using Hölder's inequality, the last term can be estimated as:

$$ \left| \int_0^t \sigma \int_{\mathbb{R}^d} F_t^{jk} \partial_k u_{tt}^j \right| \le \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \frac{1}{4\eta} \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 $$
whereas, in $d = 3$, one has:

$$ \begin{aligned} \int_0^T \sigma^{4/3} \|(-\Delta)^{-1} \nabla P(\rho_0 u_{tt}), (-\Delta)^{-1} \nabla Q(\rho_0 u_{tt})\|_{L^4(\mathbb{R}^3)}^{8/3} &\lesssim \int_0^T \sigma^{4/3} \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathbb{D}u\|_{L^2(\mathbb{R}^3)}^{\frac{2}{3}} \|u_{tt}\|_{L^2(\mathbb{R}^3)}^2 \\ &\lesssim \left[ \int_0^T \sigma^2 \|u_{tt}\|_{L^2(\mathbb{R}^3)}^3 \right]^{2/3} \left[ \int_0^T \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathbb{D}u\|_{L^2(\mathbb{R}^3)}^2 \right]^{1/3} . \end{aligned} $$

Finally, one has:

$$ (3.29) \qquad \int_0^T \sigma^{4/d} \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \lesssim \int_0^T \sigma^{4/d} \|F_t\|_{L^4(\mathbb{R}^d)}^{8/d} + \int_0^T \sigma^{4/d} \|\nabla u\|_{L^8(\mathbb{R}^d)}^{8/d} \|\mu_t, \lambda_t\|_{L^8(\mathbb{R}^d)}^{8/d} + \sigma^{8/d} \|\sigma u_{tt}\|_{L^\infty((0,T),L^2(\mathbb{R}^d))}^{8/d} \\ \qquad \qquad \qquad \qquad \qquad + \left[ \int_0^T \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathbb{D}u\|_{L^2(\mathbb{R}^d)}^2 \right]^{4/d} . $$

Sixth step: Piecewise Hölder norm estimate for $\nabla u$. In this step, we apply the quantitative Hölder estimate Proposition A.2 to (3.19), in order to obtain:

$$ \begin{aligned} \|\nabla u\|_{\mathscr{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)} &\lesssim \frac{1}{\mu} \|(\mu - \tilde{\mu})\mathbb{D}u\|_{\mathscr{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)} + \|F\|_{\mathscr{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)} + \eta^{-\frac{d}{2}} (\|\mathbb{D}u\|_{L^2(\mathbb{R}^d)} + \|F\|_{L^2(\mathbb{R}^d)}) \\ &\quad + \|(-\Delta)^{-1} \nabla P(\rho_0 u_t)\|_{\mathscr{C}^{\alpha}(\mathbb{R}^d)} + \|(-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{\mathscr{C}^{\alpha}(\mathbb{R}^d)} + \left\| \frac{1}{\mu} (\mu - \tilde{\mu}) \mathbb{D}u, F \right\|_{L^{\infty}(\mathbb{R}^d)} \mathcal{P}(C). \end{aligned} $$

Thus, there exists a constant $[\mu] > 0$ that depends on $d$ and the regularity of the interface such that if

$$ \frac{1}{\mu} \|\mu - \tilde{\mu}\|_{L^{\infty}((0,T), \mathscr{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d))} \le [\mu] \approx \frac{1}{1 + \mathcal{P}(C)} $$

then the first and last terms of the right hand side above can be absorbed in the left hand side and this will lead to:

$$ \begin{aligned} \|\nabla u\|_{\mathscr{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)} &\lesssim \|\mathbb{D}u\|_{L^2(\mathbb{R}^d)} + \|F\|_{L^2(\mathbb{R}^d)} + (1 + \mathcal{P}(C)) \|F\|_{\mathscr{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)} \\ &\quad + \|(-\Delta)^{-1} \nabla P(\rho_0 u_t)\|_{\mathscr{C}^{\alpha}(\mathbb{R}^d)} + \|(-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{\mathscr{C}^{\alpha}(\mathbb{R}^d)}. \end{aligned} $$

We only have to compute the time norm of the two last term. By interpolation, one has:

$$ (3.30) \qquad \begin{aligned} \int_0^T \sigma^{\frac{d\beta}{4}} \|(-\Delta)^{-1} \nabla P(\rho_0 u_t); (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{L^4(\mathbb{R}^d)}^{4} &\lesssim \int_0^T \sigma^{\frac{d\beta}{4}} \|\rho_0 u_t\|_{L^4(\mathbb{R}^d)}^{4\beta} \|\nabla(-\Delta)^{-1}(\rho_0 u_t)\|_{L^2(\mathbb{R}^d)}^{4(1-\beta)} \\ &\lesssim (\rho^*)^{4\beta} \int_0^T \sigma^{\frac{d\beta}{4}} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^{d\beta} \|u_t\|_{L^2(\mathbb{R}^d)}^{(4-d)\beta} \|F, \nabla u\|_{L^2(\mathbb{R}^d)}^{4(1-\beta)} \\ &\lesssim \left[ \int_0^T \sigma \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^4 \right]^{\frac{d\beta}{4}} \left[ \int_0^T \|u_t\|_{L^2(\mathbb{R}^d)}^4 \right]^{(1-\frac{d}{4})\beta} \left[ \int_0^T \|\nabla(-\Delta)^{-1}(\rho_0 u_t)\|_{L^2(\mathbb{R}^d)}^4 \right]^{1-\beta} \end{aligned} $$

with

$$ \beta = \frac{2d}{4+d}. $$

On the other hand, one has:

$$ (3.31) \qquad \int_0^T \sigma^r \|(-\Delta)^{-1} \nabla P(\rho_0 u_t); (-\Delta)^{-1} \nabla Q(\rho_0 u_t)\|_{\mathscr{C}^{\alpha}(\mathbb{R}^d)}^4 \lesssim \int_0^T \sigma^r \|u_t\|_{L^{d/(1-\alpha)}(\mathbb{R}^d)}^4. $$

In two space dimension, one has:

$$ \begin{aligned} \int_0^T \sigma^\alpha \|u_t\|_{L^{2/(1-\alpha)}(\mathbb{R}^2)}^4 &\lesssim \int_0^T \sigma^\alpha \|\nabla u_t\|_{L^2(\mathbb{R}^2)}^{4\alpha} \|u_t\|_{L^2(\mathbb{R}^2)}^{4(1-\alpha)} \\ &\lesssim \left[ \int_0^T \sigma \|\nabla u_t\|_{L^2(\mathbb{R}^2)}^4 \right]^\alpha \left[ \int_0^T \|u_t\|_{L^2(\mathbb{R}^2)}^4 \right]^{1-\alpha} \end{aligned} $$
and the remaining terms are bounded by:

$$
\begin{aligned}
& \left| \frac{1}{2} \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t |\operatorname{div} u_t|^2 - \int_0^t \sigma \int_{\mathbb{R}^d} \lambda_t \operatorname{div} u \operatorname{div} u_{tt} + \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t |\mathbb{D}u_t|^2 - 2 \int_0^t \sigma \int_{\mathbb{R}^d} \mu_t \mathbb{D}^{jk} u \partial_k u_{tt}^j \right| \\
& \le \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \frac{C_*}{\eta} \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \\
& \qquad + C_* \int_0^t \sigma \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)}.
\end{aligned}
$$

Gathering all of these estimates, we have:

$$
\begin{aligned}
& \int_0^t \sigma \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \sigma(t) \int_{\mathbb{R}^d} \mu |\mathbb{D}u_t|^2 + \frac{\sigma(t)}{2} \int_{\mathbb{R}^d} \lambda |\operatorname{div} u_t|^2 \\
& \le C_* \int_0^t \sigma' \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + 2\eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \frac{C_*}{\eta} \int_0^t \|F_t\|_{L^2(\mathbb{R}^d)}^2 \\
& \qquad + \frac{C_*}{\eta} \int_0^t \|\mu_t, \lambda_t\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 + C_* \int_0^t \sigma \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2.
\end{aligned}
\tag{3.30}
$$

The second term on the right hand side above will be reduced with the terms coming from the second energy which we turn to. We differentiate (3.28) with respect to time before multiplying the resulting equation with $\sigma^2 u_{tt}$. By doing so, we obtain:

$$
\begin{aligned}
& \frac{\sigma^2}{2} \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + 2 \int_0^t \sigma^2 \|\sqrt{\mu} \mathbb{D} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma^2 \|\sqrt{\lambda} \operatorname{div} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 \\
& = \int_0^t \sigma \sigma' \|\sqrt{\rho_0} u_{tt}\|_{L^2(\mathbb{R}^d)}^2 \\
& \quad - 4 \int_0^t \sigma^2 \mu_t \mathbb{D}^{jk} u_t \partial_k u_{tt}^j - 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_t \operatorname{div} u_t \operatorname{div} u_{tt} - \int_0^t \sigma^2 \int_{\mathbb{R}^d} F_{tt}^{jk} \partial_k u_{tt}^j \\
& \quad - 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \mu_{tt} \mathbb{D}^{jk} u \partial_k u_{tt}^j - \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_{tt} \operatorname{div} u \operatorname{div} u_{tt}.
\end{aligned}
$$

We estimate the second and third terms of the right hand side as follows:

$$
\left| 4 \int_0^t \sigma^2 \mu_t \mathbb{D}^{jk} u_t \partial_k u_{tt}^j + 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_t \operatorname{div} u_t \operatorname{div} u_{tt} \right| \le \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \frac{C_*}{\eta} \int_0^t \sigma^2 \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)}^2 \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2,
$$

whereas the last two terms are bounded by

$$
\begin{aligned}
& \left| 2 \int_0^t \sigma^2 \int_{\mathbb{R}^d} \mu_{tt} \mathbb{D}^{jk} u \partial_k u_{tt}^j + \int_0^t \sigma^2 \int_{\mathbb{R}^d} \lambda_{tt} \operatorname{div} u \operatorname{div} u_{tt} \right| \\
& \le \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \frac{C_*}{\eta} \int_0^t \sigma^2 \|\mu_{tt}, \lambda_{tt}\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2.
\end{aligned}
$$

Finally, the fourth term of the right hand side is:

$$
\left| \int_0^t \sigma^2 \int_{\mathbb{R}^d} F_{tt}^{jk} \partial_k u_{tt}^j \right| \le \eta \int_0^t \sigma^2 \|\nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \frac{1}{4\eta} \int_0^t \sigma^2 \|F_{tt}\|_{L^2(\mathbb{R}^d)}^2.
$$

Gathering all of these estimates, adding (3.30) and choosing $\eta$ small related to the lower bounds of the viscosity, we finally apply Grönwall's Lemma, leading to (noticing $\sigma'(t) = 0$ for $t > 1$):

$$
\int_0^T \|\sigma u_{tt}, \sigma^2 \nabla u_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \|\sigma \nabla u_t, \sigma^2 u_{tt}\|_{L^\infty((0,T),L^2(\mathbb{R}^d))}^2 \le C_* \left[ \int_0^{\sigma(T)} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|F_t, \sigma^2 F_{tt}\|_{L^2(\mathbb{R}^d)}^2 \right. \\
\left. \qquad + \int_0^T \|\mu_t, \lambda_t, \sigma \mu_{tt}, \sigma \lambda_{tt}\|_{L^4(\mathbb{R}^d)}^2 \|\nabla u\|_{L^4(\mathbb{R}^d)}^2 \right] \exp \left( C_* \int_0^T \left[ \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)} + \sigma^2 \|\mu_t, \lambda_t\|_{L^\infty(\mathbb{R}^d)}^2 \right] \right).
$$

This concludes (3.9) and we turn to (3.10). □

<figure><img src="image_2.png" /></figure>
whereas in three space dimension, we only have for all $0 < \alpha < 1/2$

$$ \begin{aligned} \int_0^T \sigma^{2\alpha} \|u_t\|_{L^{3/(1-\alpha)}(\mathbb{R}^3)}^4 &\le \int_0^T \sigma^{2\alpha} \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{2(1+2\alpha)} \|u_t\|_{L^2(\mathbb{R}^3)}^{2(1-2\alpha)} \\ &\le \sup_{[0,T]} \|u_t\|_{L^2(\mathbb{R}^3)}^{2(1-2\alpha)} \int_0^T \sigma^{2\alpha} \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{2(1+2\alpha)}. \end{aligned} $$

For $1/2 \le \alpha < 1$, one has:

$$ \begin{aligned} \int_0^T \sigma^{2\alpha} \|u_t\|_{L^{3/(1-\alpha)}(\mathbb{R}^3)}^4 &\lesssim \int_0^T \sigma^{2\alpha} \|\nabla u_t\|_{L^4(\mathbb{R}^3)}^{\frac{16}{3}(\alpha-\frac{1}{2})} \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{\frac{4}{3}(5-4\alpha)} \\ &\lesssim \left[ \|\sqrt{\sigma}\nabla u_t\|_{L^\infty((0,T),L^2(\mathbb{R}^3))}^2 \right]^{\frac{2}{3}(2-\alpha)} \left[ \int_0^T \sigma^{4/3} \|\nabla u_t\|_{L^4(\mathbb{R}^3)}^{8/3} \right]^{2\alpha-1} \left[ \int_0^T \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^2 \right]^{2(1-\alpha)}. \end{aligned} $$

In conclusion, we choose $r$ to be

$$ (3.32) \qquad r = \max \left\{ \frac{d^2}{2(4+d)}, 2\alpha \right\} < 2 $$

and we have:

$$ (3.33) \qquad \begin{aligned} \int_0^T \sigma^r \|\nabla u\|_{\mathcal{E}_{p,w,\gamma}^{\alpha}(\mathbb{R}^d)}^4 &\lesssim (1 + \mathcal{P}(C)) \int_0^T \sigma^r \|F\|_{\mathcal{E}_{p,w,\gamma}^{\alpha}(\mathbb{R}^d)}^4 + \|F, \nabla u, u_t, \sqrt{\sigma}\nabla u_t\|_{L^\infty((0,T),L^2(\mathbb{R}^d))}^4 \\ &\quad + \left[ \int_0^T \|F, \nabla u, u_t, \nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \right]^2 + \left[ \int_0^T \sigma^{4/d} \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \right]^{d/2}. \end{aligned} $$

This completes the first part of the proof of Theorem 3.1. $\square$

We turn to the second step in the proof of Theorem 3.1 that consists in proving that the set $\mathcal{E}$, defined in (3.12), is not empty, that is $0 \in \mathcal{E}$. We have already discussed this fact linearly on page 10¹, and now it only remains to prove Theorem 3.4, which is the purpose of the following section.

### 3.3.2. Proof of Theorem 3.4.
The proof of Theorem 3.4 is based on a homotopy argument and for this purpose, we define

$$ \mu_0^{\theta} = (1 - \theta)\tilde{\mu}' + \theta\mu_0' \quad \text{and} \quad \lambda_0^{\theta} = (1 - \theta)\tilde{\lambda}' + \theta\lambda_0', $$

and we consider the equation (3.16) where the viscosity $\mu_0'$, $\lambda_0'$ are replaced by $\mu_0^{\theta}$, $\lambda_0^{\theta}$ respectively. Also, we define

$$ X_T' := \{v \in \mathcal{C}([0,T], H^1(\mathbb{R}^d)): v_t, \nabla^2 v \in L^2((0,T) \times \mathbb{R}^d)\} $$

and the map:

$$ \Pi_{\theta}: \begin{aligned} X_T' &\to L^2((0,T) \times \mathbb{R}^d) \times H^1(\mathbb{R}^d) \\ v &\mapsto (v_t - \operatorname{div}(2\mu_0^{\theta} \mathbb{D}v + \lambda_0^{\theta} \operatorname{div} vI_d) - \operatorname{div}(b_0 : v), v|_{t=0}) \end{aligned} $$

and finally we consider the set:

$$ \mathcal{E}'_T := \{\theta \in [0,1]: \Pi_{\theta} \text{ is one to one}\}. $$

Obviously, the map $\Pi_{\theta}$ is well-defined and continuous, we will only have to show that $1 \in \mathcal{E}'_T$ to conclude the proof of Theorem 3.4. To obtain $0 \in \mathcal{E}$, one can observe that for $\theta = 0$, the equation (3.16)$_\theta$ is reduced to the Lamé system:

$$ (3.34) \qquad v_t - \tilde{\mu}' \Delta u - (\tilde{\mu}' + \tilde{\lambda}') \nabla \operatorname{div} v = f + \operatorname{div}(b_0 : v) $$

in such a way that $v = \mathcal{P}v + \mathcal{Q}v$ with

$$ \left\{ \begin{array}{l} \partial_t \mathcal{A}v - k \mathcal{A}v = \mathcal{A}f + \mathcal{A} \operatorname{div}(b_0 : v), \\ \mathcal{A}v|_{t=0} = \mathcal{A}v_0 \end{array} \right. $$

where $k = \tilde{\mu}'$ if $\mathcal{A} = \mathcal{P}$ and $k = 2\tilde{\mu}' + \tilde{\lambda}'$ if $\mathcal{A} = \mathcal{Q}$. In the following theorem we address the global well-posedness of (3.34).

<footnote n="1">¹from (3.13)</footnote>
**Proof of (3.10).** Here we derive $L^{8/d}((0, T), L^4(\mathbb{R}^d))$-norm estimate for $\sigma^{\frac{d}{4}}\nabla u_t$. Differentiating (3.19)-(3.20) with respect to time yields:

$$ 
\begin{cases} 
\tilde{\mu} \mathcal{P} u_t = (-\Delta)^{-1} \mathcal{P} \operatorname{div}\{2(\mu - \tilde{\mu})\mathcal{D}u_t + 2\mu_t \mathcal{D}u + F_t\} - (-\Delta)^{-1} \mathcal{P}(\rho_0 u_{tt}), \\ 
(2\tilde{\mu} + \lambda) \operatorname{div} u_t = -\lambda_t \operatorname{div} u + (-\Delta)^{-1} \operatorname{div}\{2(\mu - \tilde{\mu})\mathcal{D}u_t + 2\mu_t \mathcal{D}u + F_t\} - (-\Delta)^{-1} \operatorname{div}(\rho_0 u_{tt}) 
\end{cases} 
$$

and therefore (recall $Qw = -\nabla(-\Delta)^{-1} \operatorname{div} w$):

$$ 
\begin{aligned} 
\nabla u_t &= \nabla \mathcal{P} u_t + \nabla \mathcal{Q} u_t \\ 
&= \frac{1}{\tilde{\mu}} (-\Delta)^{-1} \nabla [\mathcal{P} \operatorname{div}\{2(\mu - \tilde{\mu})\mathcal{D}u_t + 2\mu_t \mathcal{D}u + F_t\} - \mathcal{P}(\rho_0 u_{tt})] + \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div}(\rho_0 u_{tt}) \right] \\ 
&\quad + \nabla^2 (-\Delta)^{-1} \left[ \frac{\lambda_t}{2\tilde{\mu} + \lambda} \operatorname{div} u \right] - \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathcal{D}u_t + 2\mu_t \mathcal{D}u + F_t\} \right]. 
\end{aligned} 
$$

Whence, under the hypothesis (3.5) we obtain:

$$ \|\nabla u_t\|_{L^4(\mathbb{R}^d)} \le C_* \|F_t\|_{L^4(\mathbb{R}^d)} + C_* \|\nabla u\|_{L^8(\mathbb{R}^d)} \|\mu_t, \lambda_t\|_{L^8(\mathbb{R}^d)} + C_* \|(-\Delta)^{-1}\nabla(\rho_0 u_{tt})\|_{L^4(\mathbb{R}^d)} $$

and we only need to estimate the last term. Following the computations leading to (3.22), we have from (3.28):

$$ (3.31) \qquad \|(-\Delta)^{-1}\nabla(\rho_0 u_{tt})\|_{L^2(\mathbb{R}^d)} \le C_* \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathcal{D}u\|_{L^2(\mathbb{R}^d)}. $$

and Gagliardo-Nirenberg's inequality implies:

$$ 
\begin{aligned} 
\int_0^T \sigma^2 \|(-\Delta)^{-1} \nabla(\rho_0 u_{tt})\|_{L^4(\mathbb{R}^d)}^{8/d} \le C_* \int_0^T \sigma^2 \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathcal{D}u\|_{L^2(\mathbb{R}^d)}^{\frac{8}{d}-2} \|u_{tt}\|_{L^2(\mathbb{R}^d)}^2 \\ 
\le C_* \|\sigma u_{tt}\|_{L^\infty((0,T), L^2(\mathbb{R}^2))}^2 \int_0^T \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathcal{D}u\|_{L^2(\mathbb{R}^d)}^{\frac{8}{d}-2}. 
\end{aligned} 
$$

In summary, we obtain:

$$ 
\begin{aligned} 
\int_0^T \sigma^2 \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \le C_* \int_0^T \sigma^2 \|F_t\|_{L^4(\mathbb{R}^d)}^{8/d} + C_* \int_0^T \sigma^2 \|\nabla u\|_{L^8(\mathbb{R}^d)}^{8/d} \|\mu_t, \lambda_t\|_{L^8(\mathbb{R}^d)}^{8/d} \\ 
+ C_* \|\sigma u_{tt}\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^2 \int_0^T \|\nabla u_t, F_t, \lambda_t \operatorname{div} u, \mu_t \mathcal{D}u\|_{L^2(\mathbb{R}^d)}^{\frac{8}{d}-2}. 
\end{aligned} 
$$

This ends the proof of (3.10) and we turn to the last a priori estimate. $\square$

**Proof of (3.11).** Here, we translate the bounds of $u_t$ into a piecewise Hölder regularity for the velocity gradient. To achieve this, we will use the results from Appendix A to derive a piecewise Hölder estimate for the discontinuous terms in the expression of the velocity gradient (3.21). We fix $q = 2d/\alpha$.

**Step 1: Terms involving F.** We estimate the terms

$$ \frac{1}{\tilde{\mu}} (-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div} F \quad \text{and} \quad \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div} F \right]. $$

With the help of (A.2) and (A.3), we have:

$$ (3.32) \qquad \begin{cases} \|(-\Delta)^{-1} \nabla^2 F, (-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div} F\|_{L^\infty(\mathbb{R}^d)} \le C \left( \ell_{\varphi_0}^{-\frac{\alpha}{2}} \|F\|_{L^q(\mathbb{R}^d)} + \|F\|_{\mathcal{C}^{\alpha}_{pw,c_0}(\mathbb{R}^d)} \right), \\ \|(-\Delta)^{-1} \nabla^2 F, (-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div} F\|_{\dot{\mathcal{C}}^{\alpha}_{pw,c_0}(\mathbb{R}^d)} \le C(1+\mathfrak{P}_{C_0}) \|F\|_{\mathcal{C}^{\alpha}_{pw,c_0}(\mathbb{R}^d)}, \end{cases} $$

where $\mathfrak{P}_{C_0}$ is defined in (1.4). Similarly, we have:

$$ 
\begin{aligned} 
& \left\| \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div} F \right] \right\|_{L^\infty(\mathbb{R}^d)} \\ 
& \le C \left[ \left( \|\lambda\|_{\dot{\mathcal{C}}^{\alpha}_{pw,c_0}(\mathbb{R}^d)} + \ell_{\varphi_0}^{-\frac{\alpha}{2}} \right) \|F\|_{L^q(\mathbb{R}^d)} + \|(-\Delta)^{-1} \operatorname{div} \operatorname{div} F\|_{\mathcal{C}^{\alpha}_{pw,c_0}(\mathbb{R}^d)} \right] \\ 
& \le C \left[ \left( \|\lambda\|_{\dot{\mathcal{C}}^{\alpha}_{pw,c_0}(\mathbb{R}^d)} + \ell_{\varphi_0}^{-\frac{\alpha}{2}} \right) \|F\|_{L^q(\mathbb{R}^d)} + (1 + \mathfrak{P}_{C_0}) \|F\|_{\mathcal{C}^{\alpha}_{pw,c_0}(\mathbb{R}^d)} \right], 
\end{aligned} 
$$
**Theorem 3.6.** Let $0 < T < \infty$, $v_0 \in H^1(\mathbb{R}^d)$, $b_0 \in W^{1,\infty}(\mathbb{R}^d)$ and $f \in L^2((0,T) \times \mathbb{R}^d)$. Then, the Cauchy problem associated to the equation (3.34) with initial datum $v_0$ admits an unique solution $v \in \mathcal{C}([0,T], H^1(\mathbb{R}^d))$. The solution $v$ verifies the following integral formulation

$$ \mathcal{A}v(t) = \mathcal{T}(kt)\mathcal{A}v_0 + \int_0^t \mathcal{T}(k(t-s))\mathcal{A}\{f + \operatorname{div}(b_0v)\}(s)ds $$

where $(\mathcal{T}(t))_t$ is the heat semi-group. Moreover, there are constants $K, \Lambda$ depending on $b_0$ and $\tilde{\mu}, \tilde{\lambda}$ such that:

$$ \sup_{t \in [0,T]} e^{-\Lambda t} \|v(t)\|_{H^1(\mathbb{R}^d)} \leq K R_T $$

where

$$ R_T := \sup_{[0,T]} \left\| \mathcal{T}(kt)\mathcal{A}v_0 + \int_0^t \mathcal{T}(k(t-s))\mathcal{A}f(s)ds \right\|_{H^1(\mathbb{R}^d)} . $$

As well, $v_t, \nabla^2 v \in L^2((0,T) \times \mathbb{R}^d)$ with the following estimate:

$$ \sup_{[0,T]} e^{-2\Lambda t} \int_0^t \|v_t(s), \nabla^2 v(s)\|_{L^2(\mathbb{R}^d)}^2 ds \leq K^2 R_T^2. $$

The proof of this theorem is done in Appendix B and is based on some appropriate estimate on the heat kernel. From this discussion, it turns out that $\mathcal{E}'_T$ is a non-empty set. The remainder of Theorem 3.4 consists in first deriving some a priori estimates for the solution $v \in X'_T$ of the equation (3.16) in Section 3.3.2. Then, in Section 3.3.2, we prove that $\mathcal{E}'_T$ is at the same time open and closed in much easier way.

A priori estimates. We establish three types estimates for the solution $v \in X'_T$ of the equation (3.16). The first one is the classical energy that consists in multiplying the equation (3.16) by $v$ and integrating in time and in space in order to obtain:

$$ \begin{aligned} \frac{1}{2} \|v(t)\|_{L^2(\mathbb{R}^d)}^2 &+ \int_0^t \int_{\mathbb{R}^d} 2\mu_0' |\mathbb{D}v|^2 + \int_0^t \int_{\mathbb{R}^d} \lambda_0' |\operatorname{div} v|^2 = \frac{1}{2} \|v_0\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \int_{\mathbb{R}^d} fv - \int_0^t \int_{\mathbb{R}^d} b_0 : v \nabla v \\ &\leq \eta \int_0^t \|\nabla v\|_{L^2(\mathbb{R}^d)}^2 + \frac{1}{2} \|v_0\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|f\|_{L^2(\mathbb{R}^d)}^2 + C_\eta \int_0^t (1 + \|b_0\|_{L^\infty(\mathbb{R}^d)}^2) \|v\|_{L^2(\mathbb{R}^d)}^2 \end{aligned} $$

for some positive constant $\eta$ small related to the lower bound of the viscosity, in such way that the first term of the right hand side can be absorbed in the left hand side in order to obtain:

$$ (3.35) \qquad \|v(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|\nabla v\|_{L^2(\mathbb{R}^d)}^2 \lesssim \|v_0\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|f\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t (1 + \|b_0\|_{L^\infty(\mathbb{R}^d)}^2) \|v\|_{L^2(\mathbb{R}^d)}^2. $$

The second estimate is obtained by testing the equation (3.16) with $\Delta v$. The first term is:

$$ \int_{\mathbb{R}^d} v_t \Delta v = - \int_{\mathbb{R}^d} \nabla v_t \cdot \nabla v = - \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} |\nabla v|^2 $$

and the second one is:

$$ \begin{aligned} -2 \int_{\mathbb{R}^d} \partial_k (\mu_0' \mathbb{D}^{jk} v) \Delta v^j & - \int_{\mathbb{R}^d} \partial_j (\lambda_0' \operatorname{div} v) \Delta v^j = - \int_{\mathbb{R}^d} \mu_0' |\Delta v|^2 - \int_{\mathbb{R}^d} \mu_0' \partial_j \operatorname{div} v \Delta v^j - \int_{\mathbb{R}^d} \partial_k \mu_0' \mathbb{D}^{jk} v \Delta v^j \\ & - \int_{\mathbb{R}^d} \lambda_0' |\nabla \operatorname{div} v|^2 - \int_{\mathbb{R}^d} \operatorname{div} v \nabla \lambda_0' \cdot \nabla \operatorname{div} v \\ & = - \int_{\mathbb{R}^d} \mu_0' |\Delta v|^2 - \int_{\mathbb{R}^d} (\mu_0' + \lambda_0') |\nabla \operatorname{div} v|^2 + \int_{\mathbb{R}^d} \operatorname{div} v \partial_j \mu_0' \Delta v^j \\ & - \int_{\mathbb{R}^d} \operatorname{div} v \nabla \mu_0' \cdot \nabla \operatorname{div} v - \int_{\mathbb{R}^d} \partial_k \mu_0' \mathbb{D}^{jk} v \Delta v^j - \int_{\mathbb{R}^d} \operatorname{div} v \nabla \lambda_0' \cdot \nabla \operatorname{div} v. \end{aligned} $$

Gathering all of these computations, one has:

$$ \begin{aligned} \frac{1}{2} \frac{d}{dt} \|\nabla v\|_{L^2(\mathbb{R}^d)}^2 &+ \int_{\mathbb{R}^d} \mu_0' |\Delta v|^2 + \int_{\mathbb{R}^d} (\mu_0' + \lambda_0') |\nabla \operatorname{div} v|^2 = \int_{\mathbb{R}^d} f^j \Delta v^j + \int_{\mathbb{R}^d} \operatorname{div}(b_0 : v) \Delta v^j \\ &+ \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} \operatorname{div} v \partial_j \mu_0' \Delta v^j - \int_{\mathbb{R}^d} \operatorname{div} v \nabla \mu_0' \cdot \nabla \operatorname{div} v - \int_{\mathbb{R}^d} \partial_k \mu_0' \mathbb{D}^{jk} v \Delta v^j - \int_{\mathbb{R}^d} \operatorname{div} v \nabla \lambda_0' \cdot \nabla \operatorname{div} v. \end{aligned} $$
and

$$
\begin{aligned}
\left\| \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div} F \right] \right\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)}
&\le C \|(-\Delta)^{-1} \operatorname{div} \operatorname{div} F \|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
&\quad + C \left( \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \mathfrak{P}_{c_0} \right) \|(-\Delta)^{-1} \operatorname{div} \operatorname{div} F \|_{L^{\infty}(\mathbb{R}^d)} \\
&\le C \left( 1 + \mathfrak{P}_{c_0} + \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \right) \|F\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
&\quad + C \left( \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \mathfrak{P}_{c_0} \right) \ell_{\varphi_0}^{-\frac{\alpha}{2}} \|F\|_{L^q(\mathbb{R}^d)}.
\end{aligned}
$$

In sum, we have:

$$
\begin{equation} \tag{3.33}
\begin{aligned}
\| (-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div} F, \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div} F \right] \|_{\mathcal{E}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)}
&\le C \left( 1 + \mathfrak{P}_{c_0} + \|\lambda\|_{\mathcal{E}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \right) \left( \ell_{\varphi_0}^{-\frac{\alpha}{2}} \|F\|_{L^q(\mathbb{R}^d)} + \|F\|_{\mathcal{E}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \right).
\end{aligned}
\end{equation}
$$

**Step 2: Small terms.** We estimate the terms

$$
\frac{1}{\tilde{\mu}} (-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div}\{2(\mu - \tilde{\mu}) \mathbb{D}u\} \quad \text{and} \quad \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu}) \mathbb{D}u\} \right].
$$

Again, with the help of (A.2), we have:

$$
\begin{equation} \tag{3.34}
\begin{aligned}
\| (-\Delta)^{-1} \nabla^2 \{2(\mu - \tilde{\mu}) \mathbb{D}u \} \|_{L^{\infty}(\mathbb{R}^d)}
&\le C \left( \|\mu\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \ell_{\varphi_0}^{-\frac{\alpha}{2}} \|\mu - \tilde{\mu}\|_{L^{\infty}(\mathbb{R}^d)} \right) \|\nabla u\|_{L^q(\mathbb{R}^d)} \\
&\quad + C \|\mu - \tilde{\mu}\|_{L^{\infty}(\mathbb{R}^d)} \|\nabla u\|_{\mathcal{E}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)}.
\end{aligned}
\end{equation}
$$

Next, we gather Lemma A.2 and (A.3) to obtain:

$$
\begin{aligned}
\|(-\Delta)^{-1} \nabla^2(w)\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} &\le \|(-\Delta)^{-1} \nabla^2(w^e)\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \|(-\Delta)^{-1} \nabla^2(w - w^e)\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
&\le C \left[ \|w^e\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \|w - w^e\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \mathfrak{P}_{c_0} \|w - w^e\|_{L^{\infty}(\mathbb{R}^d)} \right] \\
&\le C \left[ \|w\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + (\ell_{\varphi_0}^{-\alpha} + \mathfrak{P}_{c_0}) \|\!\|w\!\|\!\|_{L^{\infty}(c_0)} \right].
\end{aligned}
$$

Therefore:

$$
(3.35) \quad \| (-\Delta)^{-1} \nabla^2 \{2(\mu - \tilde{\mu}) \mathbb{D}u \} \|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \le C \left[ \| (\mu - \tilde{\mu}) \mathbb{D}u \|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + (\ell_{\varphi_0}^{-\alpha} + \mathfrak{P}_{c_0}) \|\!\| (\mu - \tilde{\mu}) \mathbb{D}u \!\|\!\|_{L^{\infty}(c_0)} \right].
$$

Similarly, we have:

$$
\begin{equation} \tag{3.36}
\begin{aligned}
\left\| \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div} \{2(\mu - \tilde{\mu}) \mathbb{D}u \} \right] \right\|_{L^{\infty}(\mathbb{R}^d)}
&\le C \left( \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \ell_{\varphi_0}^{-\frac{\alpha}{2}} \right) \|\mu - \tilde{\mu}\|_{L^{\infty}(\mathbb{R}^d)} \|\nabla u\|_{L^q(\mathbb{R}^d)} \\
&\quad + C \|(-\Delta)^{-1} \operatorname{div} \operatorname{div} \{2(\mu - \tilde{\mu}) \mathbb{D}u \} \|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
&\le C \left[ \|\mu\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \left( \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \ell_{\varphi_0}^{-\frac{\alpha}{2}} \right) \|\mu - \tilde{\mu}\|_{L^{\infty}(\mathbb{R}^d)} \right] \|\nabla u\|_{L^q(\mathbb{R}^d)} \\
&\quad + C (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \|\!\| (\mu - \tilde{\mu}) \mathbb{D}u \!\|\!\|_{L^{\infty}(c_0)} + C \|\mu - \tilde{\mu}\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \|\nabla u\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)};
\end{aligned}
\end{equation}
$$
Thus, Hölder and Young inequalities yield:

$$ (3.36) $$

$$ \|\nabla v(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|\Delta v\|_{L^2(\mathbb{R}^d)}^2 \lesssim \|\nabla v_0\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|f\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|\nabla v\|_{L^2(\mathbb{R}^d)}^2 \|\nabla \mu_0, \nabla \lambda_0, b_0, \nabla b_0\|_{L^\infty(\mathbb{R}^d)}^2. $$

Finally, we gather (3.35) and (3.36), and we make use of Grönwall's lemma and equation (3.16), in order to obtain:

$$ \sup_{[0,T]} e^{-\Lambda t} \left\{ \|v(t)\|_{H^1(\mathbb{R}^d)}^2 + \int_0^t \|v_t, \nabla v, \nabla^2 v\|_{L^2(\mathbb{R}^d)}^2 \right\} \lesssim \|\nabla v_0\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|f\|_{L^2(\mathbb{R}^d)}^2 $$

$$ \Lambda = C(1 + \|\nabla \mu'_0, \nabla \lambda'_0, b_0, \nabla b_0\|_{L^\infty(\mathbb{R}^d)}^2) $$

and constant $C$ that depends on the lower and upper bounds of the viscosity $\mu'_0, \lambda'_0$. The $L^2((0, T) \times \mathbb{R}^d)$ norm estimate for $v_t$ can be obtained differently by testing the equation (3.16) with $v_t$. This ends the computations of the a priori estimates and we turn to the last step.

Closing argument. In this section, we will show that $\mathcal{E}'_T = [0, 1]$. We have already showed that $\mathcal{E}'_T$ is not empty. Also, we may notice that if $\Pi_\theta$ is invertible then, from the previous subsection, $\Pi_\theta^{-1}$ is linear continuous with a norm that depends on $\Lambda$ and time $T$. In particular, one has:

$$ K_T := \sup_{\theta \in \mathcal{E}'_T} \|\Pi_\theta^{-1}\|_{\mathcal{L}} < \infty. $$

Let us prove that $\mathcal{E}'_T$ is an open set. For this purpose, let $\theta_0 \in \mathcal{E}'_T$, then, $\Pi_{\theta_0}$ is a one to one map and one can write, for all $\theta \in [0, 1]$,

$$ \Pi_{\theta} = \Pi_{\theta_0} + (\theta - \theta_0)(\Pi_1 - \Pi_0) = \Pi_{\theta_0} (I + (\theta - \theta_0)\Pi_{\theta_0}^{-1}(\Pi_1 - \Pi_0)). $$

Thus, all $\theta$ such that $2|\theta - \theta_0|\|\Pi_{\theta_0}^{-1}(\Pi_1 - \Pi_0)\|_{\mathcal{L}(X'_T)} < 1$ belongs to $\mathcal{E}'_T$, therefore $\mathcal{E}'_T$ is an open set of $[0, 1]$. We turn to the last part. Let $(\theta_n)_n$ be a sequence of element of $\mathcal{E}'_T$ that converges to some $\theta \in [0, 1]$. As above, one can write, for some $n_0$

$$ \Pi_{\theta} = \Pi_{\theta_{n_0}} \left( I + (\theta - \theta_{n_0}) \Pi_{\theta_{n_0}}^{-1} (\Pi_1 - \Pi_0) \right). $$

Since $K_T < \infty$, one can choose large $n_0$ such that

$$ 2|\theta - \theta_{n_0}| K_T \|\Pi_1 - \Pi_0\|_{\mathcal{L}} < 1 $$

thus $\Pi_\theta$ is one to one as composition of two one to one operators. Therefore $\theta \in \mathcal{E}'_T$ and consequently $\mathcal{E}'_T$ is a closed set of $[0, 1]$. One then concludes that $1 \in \mathcal{E}'_T$ and **Theorem 3.4** then follows.

### 3.3.3. Remainder of the proof of **Theorem 3.1**.

*Proof.* We recall that **Theorem 3.1** follows as soon as we will prove that the set $\mathcal{E} = [0, 1]$. What remains to be proved is that $\mathcal{H} \mathcal{E}$ is at the same time open and closed. We begin by proving that it is an open subset of $[0, 1]$. So let $\theta_0 \in \mathcal{E}$, the goal is to find some $\eta > 0$ such that

$$ (3.37) \qquad B(\theta_0, \eta) \cap [0, 1] \subset \mathcal{E}. $$

For this purpose, let $\theta \in [0, 1]$, $(F, u_0) \in Y_T^\theta$ and consider the Cauchy problem (3.11) that can be written as:

$$ (3.38) \qquad \begin{cases} \rho_0 u_t - \operatorname{div}(2\mu^{\theta_0} \mathbb{D}u + \lambda^{\theta_0} \operatorname{div} u I_d) = \operatorname{div} F + (\theta - \theta_0) \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u + (\lambda - \tilde{\lambda}) \operatorname{div} u I_d\}, \\ u|_{t=0} = u_0. \end{cases} $$

We will apply a fixed point theorem in order to solve the above Cauchy problem with a small $\eta$ that does not depend on $F$ nor $u_0$. In fact, the proof shows that the constant $\eta$ only depends on the norm of the viscosity and initial density as required in (3.2), as well on $\mathcal{P}(C)^2$. The fixed point theorem will be applied in $X_T(u_0)$ defined as:

$$ X_T(u_0) := \{v \in X_T : v|_{t=0} = u_0\} $$

a complete subset of the Banach space $X_T$. So for $v \in X_T(u_0)$, we consider the following equation:

$$ (3.39) \qquad \begin{cases} \rho_0 u_t - \operatorname{div}(2\mu^{\theta_0} \mathbb{D}u + \lambda^{\theta_0} \operatorname{div} u I_d) = \operatorname{div} F + (\theta - \theta_0) \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}v + (\lambda - \tilde{\lambda}) \operatorname{div} v I_d\}, \\ u|_{t=0} = u_0. \end{cases} $$

<footnote n="2">²see (3.10)</footnote>
and

$$
\begin{align*}
& \left\| \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} \right] \right\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
& \le C \left\| \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} \right\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
& \quad + C (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \left\| \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} \right] \right\|_{L^{\infty}(C_0)} \\
& \le C \left( \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \left\| \left[ \lambda \right] \right\|_{L^{\infty}(C_0)} \right) \|(-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} \|_{L^{\infty}(\mathbb{R}^d)} \\
& \quad + C (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \left\| \left[ (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} \right] \right\|_{L^{\infty}(C_0)} \\
& \quad + C \left\| (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} \right\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
& \le C \left( \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \left\| \left[ \lambda \right] \right\|_{L^{\infty}(C_0)} \right) \left( \|\mu\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \ell_{\varphi_0}^{-\frac{\alpha}{2}} \|\mu - \tilde{\mu}\|_{L^{\infty}(\mathbb{R}^d)} \right) \|\nabla u\|_{L^q(\mathbb{R}^d)} \\
& \quad + C \left[ 1 + \|\lambda\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \left\| \left[ \lambda \right] \right\|_{L^{\infty}(C_0)} \right] \|\mu - \tilde{\mu}\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \|\nabla u\|_{\dot{\mathcal{E}}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \\
& \quad + C (\mathfrak{P}_{c_0} + \ell_{\varphi_0}^{-\alpha}) \left\| \left[ (-\Delta)^{-1} \operatorname{div} \operatorname{div}\{2(\mu - \tilde{\mu})\mathbb{D}u\} \right] \right\|_{L^{\infty}(C_0)} \left\| (\mu - \tilde{\mu})\mathbb{D}u \right\|_{L^{\infty}(C_0)}.
\end{align*}
$$

We now turn to estimating the norms of the jumps involved in the above inequalities. We first notice that the
solution $u$ of (3.1) satisfies:

$$
(3.38) \qquad \left\| 2\mu \mathbb{D}^{jk} u + \lambda \operatorname{div} u \delta^{jk} + F^{jk} \right\| n_x^k = 0,
$$

where $n_x$ is the normal vector field of $C_0$ and $\delta^{jk}$ denotes the Kronecker symbol. We consider an unit vector field $\tau_{x,l}$, $l \in [1, d-1]$ such that $(\tau_{x,1}, \dots, \tau_{x,d-1}, n_x)$ forms a basis of $\mathbb{R}^d$. Given that $u$ is continuous in the whole space and $\nabla u$ is continuous on both sides of $C_0$, then discontinuity in $\nabla u$ may appear only in the normal direction of $C_0$. Namely, we have:

$$
\forall l \in [1, d-1], \quad \llbracket \nabla u \rrbracket \cdot \tau_{x,l} = 0
$$

and there exists a vector field $\vec{a} \in \mathbb{R}^d$ such that:

$$
(3.39) \qquad \llbracket \nabla u \rrbracket = \vec{a} \otimes n_x.
$$

If $\langle g \rangle$ denotes the average of $g$ at the interface, then we rewrite (3.38) as

$$
2\langle\mu\rangle\llbracket\mathbb{D}^{jk}u\rrbracket n_x^k + \langle\lambda\rangle\llbracket\operatorname{div}u\rrbracket n_x^j = -2\llbracket\mu\rangle\langle\mathbb{D}^{jk}u\rangle n_x^k - \llbracket\lambda\rangle\langle\operatorname{div}u\rangle n_x^j - \llbracket F^{jk}\rrbracket n_x^k,
$$

and scalar products with $n_x$ and $\tau_{x,l}$ together with (3.39) lead to:

$$
(3.40) \quad \begin{cases} \langle 2\mu + \lambda \rangle \vec{a} \cdot n_x = -2 \llbracket \mu \rrbracket \langle \mathbb{D}^{jk} u \rangle n_x^k n_x^j - \llbracket \lambda \rrbracket \langle \operatorname{div} u \rangle - \llbracket F^{jk} \rrbracket n_x^k n_x^j, \\ \langle \mu \rangle \vec{a} \cdot \tau_{x,l} = -2 \llbracket \mu \rrbracket \langle \mathbb{D}^{jk} u \rangle n_x^k \tau_{x,l}^j - \llbracket F^{jk} \rrbracket n_x^k \tau_{x,l}^j. \end{cases}
$$

Therefore:

$$
\begin{align*}
\llbracket (\mu - \tilde{\mu}) \mathbb{D} u \rrbracket &= \llbracket \mu \rrbracket \langle \mathbb{D} u \rangle + \langle \mu - \tilde{\mu} \rangle \llbracket \mathbb{D} u \rrbracket \\
&= \llbracket \mu \rrbracket \langle \mathbb{D} u \rangle + \frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle} \langle \mu \rangle \llbracket \mathbb{D} u \rrbracket
\end{align*}
$$

and recalling (3.39), (3.40) with

$$
\vec{a} = \vec{a} \cdot n_x n_x + \sum_l \vec{a} \cdot \tau_{x,l} \tau_{x,l},
$$

we infer:

$$
(3.41) \qquad \begin{aligned}
\| \llbracket (\mu - \tilde{\mu}) \mathbb{D} u \rrbracket \|_{L^{\infty}(C_0)} &\le C_* \|F\|_{L^{\infty}(\mathbb{R}^d)} + C \| \llbracket \lambda \rrbracket \|_{L^{\infty}(C_0)} \left\| \frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle} \right\|_{L^{\infty}(C_0)} \| \nabla u \|_{L^{\infty}(\mathbb{R}^d)} \\
&\quad + C \| \llbracket \mu \rrbracket \|_{L^{\infty}(C_0)} \left( 1 + \left\| \frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle} \right\|_{L^{\infty}(C_0)} \right) \| \nabla u \|_{L^{\infty}(\mathbb{R}^d)}.
\end{aligned}
$$
The regularity of $v$ and $\mu, \lambda$ ensures that the right hand side can be written as $\operatorname{div} \tilde{F}(v)$ with

$$ (3.40) \qquad \tilde{F}(v) := F + (\theta - \theta_0)\{2(\mu - \tilde{\mu})\mathbb{D}v + (\lambda - \tilde{\lambda})\operatorname{div} vI_d\} $$

that verifies $(\tilde{F}(v), u_0) \in Y_T^{\theta_0}$. The last condition in the definition of $Y_T^{\theta_0}$ is valid since $v|_{t=0} = u_0$. Now, owing to the fact that $\theta_0 \in \mathcal{E}$, the underlying equation admits an unique solution $\Phi(v) \in X_T(u_0)$ that verifies all of the estimates in Theorem 3.1 where $F$ is replaced by $\tilde{F}(v)$. In the following we are looking for some small $\theta - \theta_0$ and some $R_T$ that depends on the initial data such that the condition $\|v\|_{X_T} \le R_T$ implies the same for $\Phi(v)$. In order to facilitate the computations and the reading, the estimates will be modulo a constant that depends on, not only of the lower and upper bound of $\mu, \lambda, \rho_0$, but also on all of the norms of $\mu, \lambda$ as required in (3.2).

Step 1: Estimates for $\Phi(v)$ and $\partial_t\Phi(v)$. Since $\Phi(v)$ solves the Cauchy problem (3.39), $(\tilde{F}(v), u_0) \in Y_T^{\theta_0}$, and $\theta_0 \in \mathcal{E}$, one has from the a priori estimate the following estimate for $\Phi(v)$:

$$ (3.41) \qquad \sup_{[0,T]} \|\Phi(v), \partial_t\Phi(v), \nabla\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla\Phi(v), \nabla\partial_t\Phi(v), \partial_t\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 \lesssim \|u_0\|_{H^1(\mathbb{R}^d)}^2 + \|\partial_t\Phi(v)|_{t=0}\|_{L^2(\mathbb{R}^d)}^2 \\ + \sup_{[0,T]} \|\tilde{F}(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\partial_t\tilde{F}(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\tilde{F}(v)\|_{L^2 \cap L^6(\mathbb{R}^d)}^2. $$

On the first hand, from (3.40), one deduces that:

$$ (3.42) \qquad \sup_{[0,T]} \|\tilde{F}(v) - F\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\tilde{F}(v) - F\|_{L^2(\mathbb{R}^d)}^2 \lesssim |\theta - \theta_0|^2 \left\{ \sup_{[0,T]} \|\nabla v\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla v\|_{L^2(\mathbb{R}^d)}^2 \right\} \\ \lesssim |\theta - \theta_0|^2 \|v\|_{X_T}^2 $$

and

$$ \int_0^T \|\tilde{F}(v) - F\|_{L^6(\mathbb{R}^d)}^2 \lesssim |\theta - \theta_0|^2 \|v\|_{X_T}^2. $$

On the other hand, from (3.40) one can express the time derivative of $\tilde{F}(v)$ as:

$$ (3.43) \qquad \partial_t \tilde{F}(v) = F_t + (\theta - \theta_0)\{2(\mu - \tilde{\mu})\mathbb{D}v_t + (\lambda - \tilde{\lambda})\operatorname{div} v_t I_d + 2\mu_t \mathbb{D}v + \lambda_t \operatorname{div} v I_d\} $$

and consequently, it holds:

$$ (3.44) \qquad \int_0^T \|\partial_t \tilde{F}(v) - F_t\|_{L^2(\mathbb{R}^d)}^2 \lesssim |\theta - \theta_0|^2 \left\{ \|\nabla v\|_{L^4((0,T) \times \mathbb{R}^d)}^2 + \int_0^T \|\nabla v_t\|_{L^2(\mathbb{R}^d)}^2 \right\} \lesssim |\theta - \theta_0|^2 \|v\|_{X_T}^2. $$

Gathering all of these estimates, one has the following:

$$ (3.45) \qquad \sup_{[0,T]} \|\Phi(v), \partial_t\Phi(v), \nabla\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla\Phi(v), \nabla\partial_t\Phi(v), \partial_t\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 \lesssim \|(F, u_0)\|_{Y_T^\theta}^2 + |\theta - \theta_0|^2 \|v\|_{X_T}^2. $$

Also, one has the following estimate for $\nabla\Phi(v)$:

$$ (3.46) \qquad \|\nabla\Phi(v)\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \lesssim \|\tilde{F}(v)\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} + \sup_{[0,T]} \|\tilde{F}(v), \partial_t\Phi(v), \nabla\Phi(v)\|_{L^2(\mathbb{R}^d)}^2. $$

Owing to the fact that the estimate (3.45) provides a bound for $\partial_t\Phi(v)$ in $L^2((0,T) \times \mathbb{R}^d)$ and for $\Phi(v)$ in $L^\infty((0,T), H^1(\mathbb{R}^d))$ and owing to the estimate (3.42), we only need to estimate the $L^\infty((0,T), L^6(\mathbb{R}^d))$ norm of $\tilde{F}(v)$. In fact, this can be easily obtained from (3.40) as:

$$ \|\tilde{F}(v) - F\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \lesssim |\theta - \theta_0| \|\nabla v\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \lesssim |\theta - \theta_0| \|v\|_{X_T}. $$

Gathering all of these estimates, one has:

$$ (3.47) \qquad \|\nabla\Phi(v)\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \lesssim \|(F, u_0)\|_{Y_T^\theta} + |\theta - \theta_0| \|v\|_{X_T}. $$

We turn to the $L^{16}((0,T), L^8(\mathbb{R}^d))$ estimate for $\nabla\Phi(v)$. For this purpose, we recall from (3.8) that:

$$ \int_0^T \|\nabla\Phi(v)\|_{L^8(\mathbb{R}^d)}^{16} \lesssim \int_0^T \|\tilde{F}(v)\|_{L^8(\mathbb{R}^d)}^{16} + \|\tilde{F}(v), \partial_t\Phi(v), \nabla\Phi(v)\|_{L^2 \cap L^\infty((0,T), L^2(\mathbb{R}^d))}^{16} + \|\nabla\partial_t\Phi(v)\|_{L^2((0,T) \times \mathbb{R}^d)}^{16}. $$

We have already obtained bound for the last two terms of the right hand side above in (3.45) and (3.42) and we have only to estimate the first one. So from (3.40), one has:

$$ \int_0^T \|\tilde{F}(v) - F\|_{L^8(\mathbb{R}^d)}^{16} \lesssim |\theta - \theta_0|^{16} \int_0^T \|\nabla v\|_{L^8(\mathbb{R}^d)}^{16} \lesssim |\theta - \theta_0|^p \|v\|_{X_T}^{16}. $$
Moreover, from (3.20) and (3.40)$_1$, we find (recall $\llbracket \mathrm{div} \, u \rrbracket = a \cdot n_x$):

$$
\begin{align*}
\|\llbracket (-\Delta)^{-1} \mathrm{div} \, \mathrm{div}\{(\mu - \tilde{\mu})\mathbb{D}u\}\rrbracket\|_{L^{\infty}(C_0)} &\le C \|\llbracket \mu \rrbracket\|_{L^{\infty}(C_0)} \left(1 + \left\|\frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle}\right\|_{L^{\infty}(C_0)}\right) \|\nabla u\|_{L^{\infty}(\mathbb{R}^d)} \\
&= \frac{\langle 2\tilde{\mu} + \lambda \rangle}{\langle 2\mu + \lambda \rangle} \langle 2\mu + \lambda \rangle \|\mathrm{div} \, u\| + \|\lambda\| \langle \mathrm{div} \, u \rangle - \|\llbracket (-\Delta)^{-1} \mathrm{div} \, \mathrm{div} \, F\rrbracket\| \\
&= \frac{\langle 2\tilde{\mu} + \lambda \rangle}{\langle 2\mu + \lambda \rangle} (-2\|\mu\| \langle \mathbb{D}^j u \rangle n_x^k n_x^j - \|F^{jk}\| n_x^k n_x^j) \\
&\quad + 2 \frac{\langle \mu - \tilde{\mu} \rangle}{\langle 2\mu + \lambda \rangle} \|\lambda\| \langle \mathrm{div} \, u \rangle - \|\llbracket (-\Delta)^{-1} \mathrm{div} \, \mathrm{div} \, F\rrbracket\|
\end{align*}
$$

and hence (recall (3.32)):

$$
(3.42) \qquad \begin{aligned}
\|\llbracket (-\Delta)^{-1} \mathrm{div} \, \mathrm{div}\{(\mu - \tilde{\mu})\mathbb{D}u\}\rrbracket\|_{L^{\infty}(C_0)} &\le C \|\llbracket \mu \rrbracket\|_{L^{\infty}(C_0)} \left(1 + \left\|\frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle}\right\|_{L^{\infty}(C_0)}\right) \|\nabla u\|_{L^{\infty}(\mathbb{R}^d)} \\
&\quad + C \|\llbracket \lambda \rrbracket\|_{L^{\infty}(C_0)} \left\|\frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle}\right\|_{L^{\infty}(C_0)} \|\nabla u\|_{L^{\infty}(\mathbb{R}^d)} + C \left(\ell_{\varphi_0}^{-\frac{\alpha}{2}} \|F\|_{L^q(\mathbb{R}^d)} + \|F\|_{\mathcal{C}_{\mathrm{pw},C_0}^{\alpha}(\mathbb{R}^d)}\right).
\end{aligned}
$$

Summing up (3.34), (3.35), (3.36), (3.37), (3.41) and (3.42) we find:

$$
\begin{align}
& \|(-\Delta)^{-1} \nabla \mathcal{P} \operatorname{div}\{(\mu - \tilde{\mu})\mathbb{D}u\}, \nabla^2(-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div}\{(\mu - \tilde{\mu})\mathbb{D}u\} \right] \|_{\mathcal{C}_{\mathrm{pw},C_0}^{\alpha}(\mathbb{R}^d)} \nonumber \\
& \le C_{**} \left( \|\nabla u\|_{L^q(\mathbb{R}^d)} + \|F\|_{L^q(\mathbb{R}^d)} + \|F\|_{\mathcal{C}_{\mathrm{pw},C_0}^{\alpha}(\mathbb{R}^d)} \right) + C \left( \mathfrak{P}_{C_0} + \ell_{\varphi_0}^{-\alpha} \right) \|\mu\|_{L^{\infty}(C_0)} \|\nabla u\|_{L^{\infty}(\mathbb{R}^d)} \nonumber \\
& \quad + C \left( 1 + \|\lambda\|_{\mathcal{C}_{\mathrm{pw},C_0}^{\alpha}(\mathbb{R}^d)} + \left( \mathfrak{P}_{C_0} + \ell_{\varphi_0}^{-\alpha} \right) \|\lambda\|_{L^{\infty}(C_0)} \right) \|\mu - \tilde{\mu}\|_{\mathcal{C}_{\mathrm{pw},C_0}^{\alpha}(\mathbb{R}^d)} \|\nabla u\|_{\mathcal{C}_{\mathrm{pw},C_0}^{\alpha}(\mathbb{R}^d)} \nonumber \\
& \quad + C \left( \mathfrak{P}_{C_0} + \ell_{\varphi_0}^{-\alpha} \right) \|\lambda, \mu\|_{L^{\infty}(C_0)} \left\| \frac{\langle \mu - \tilde{\mu} \rangle}{\langle \mu \rangle} \right\|_{L^{\infty}(C_0)} \|\nabla u\|_{L^{\infty}(\mathbb{R}^d)} \tag{3.43}
\end{align}
$$

Above, the constant $C_{**}$ depends polynomially on the lower and upper bounds of the viscosity and density, on $\sup_{[0,T]} \|\mu, \lambda\|_{\mathcal{C}_{\mathrm{pw},C_0}^{\alpha}(\mathbb{R}^d)}$, and on $\ell_{\varphi_0}^{-\frac{\alpha}{2}}$ and $\mathfrak{P}_{C_0}$.

**Step 3: Remaining terms.** We estimate the terms

$$
\frac{1}{\tilde{\mu}} (-\Delta)^{-1} \nabla \mathcal{P}(\rho_0 u_t) \quad \text{and} \quad \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \mathrm{div}(\rho_0 u_t) \right].
$$

The first term above is Hölder continuous in the whole space, given the regularity previously obtained for $u_t$.
Indeed, by the interpolation inequality and (3.22):

$$
\begin{align}
& \int_0^T \sigma^{\frac{d\beta}{4}} \|(-\Delta)^{-1} \nabla \mathcal{P}(\rho_0 u_t); (-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{L^{\infty}(\mathbb{R}^d)}^4 \le C \int_0^T \sigma^{\frac{d\beta}{4}} \|\rho_0 u_t\|_{L^4(\mathbb{R}^d)}^{4\beta} \|\nabla(-\Delta)^{-1}(\rho_0 u_t)\|_{L^2(\mathbb{R}^d)}^{4(1-\beta)} \nonumber \\
& \le C_* \int_0^T \sigma^{\frac{d\beta}{4}} \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^{d\beta} \|u_t\|_{L^2(\mathbb{R}^d)}^{(4-d)\beta} \|F, \nabla u\|_{L^2(\mathbb{R}^d)}^{4(1-\beta)} \nonumber \\
& \le C_* \left[ \int_0^T \sigma \|\nabla u_t\|_{L^2(\mathbb{R}^d)}^4 \right]^{\frac{d\beta}{4}} \left[ \int_0^T \|u_t\|_{L^2(\mathbb{R}^d)}^4 \right]^{(1-\frac{d}{4})\beta} \left[ \int_0^T \|F, \nabla u\|_{L^2(\mathbb{R}^d)}^4 \right]^{1-\beta} \tag{3.44}
\end{align}
$$

with

$$
\beta = \frac{2d}{4+d}.
$$

Given that $\dot{W}^{1, \frac{d}{1-\alpha}}(\mathbb{R}^d) \hookrightarrow \dot{\mathcal{C}}^{\alpha}(\mathbb{R}^d)$, we have:

$$
(3.45) \qquad \int_0^T \sigma^r \|(-\Delta)^{-1} \nabla \mathcal{P}(\rho_0 u_t); (-\Delta)^{-1} \nabla(\rho_0 u_t)\|_{\dot{\mathcal{C}}^{\alpha}(\mathbb{R}^d)}^4 \le C_* \int_0^T \sigma^r \|u_t\|_{L^{d/(1-\alpha)}(\mathbb{R}^d)}^4.
$$
In consequence, one has:

$$ (3.48) \qquad \int_0^T \|\nabla\Phi(v)\|_{L^8(\mathbb{R}^d)}^{16} \lesssim \|(F, u_0)\|_{Y_T^\theta}^p + |\theta - \theta_0|^{16} \|v\|_{X_T}^{16}. $$

Step 2: Estimates for $\partial_{tt}\Phi(v)$ and $L^4$ norm estimate for $\nabla\partial_t\Phi(v)$. As in the previous step, one first recalls that $\Phi(v)$ verifies:

$$ \int_0^T \|\sqrt{\sigma}\partial_{tt}\Phi(v), \sigma\nabla\partial_{tt}\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 + \|\sqrt{\sigma}\partial_t\nabla\Phi(v), \sigma\partial_{tt}\Phi(v)\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^2 \\ \lesssim \int_0^{\sigma(T)} \|\partial_t\nabla\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\partial_t\tilde{F}(v), \sigma\partial_{tt}\tilde{F}(v)\|_{L^2(\mathbb{R}^d)}^2 + \|\nabla\Phi(v)\|_{L^\infty((0,T), L^4(\mathbb{R}^d))}^2. $$

All terms that appear in the right hand side above are controlled in (3.45), (3.47), (3.42) except the $L^2((0,T) \times \mathbb{R}^d)$ norm estimate for $\sigma\partial_{tt}\tilde{F}(v)$ that is:

$$ \int_0^T \sigma^2 \|\partial_{tt}\tilde{F}(v)\|_{L^2(\mathbb{R}^d)}^2 \le \int_0^T \sigma^2 \|F_{tt}\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \sigma^2 \|\partial_{tt}\tilde{F}(v) - F_{tt}\|_{L^2(\mathbb{R}^d)}^2. $$

In order to link the last term in terms of the norm of $v$, we first express $\partial_{tt}\tilde{F}(v) - F_{tt}$, from (3.43), as

$$ \partial_{tt}\tilde{F}(v) - F_{tt} = (\theta - \theta_0)\{2(\mu - \tilde{\mu})\mathbb{D}v_{tt} + (\lambda - \tilde{\lambda})\operatorname{div}v_{tt}I_d + 4\mu_t\mathbb{D}v_t + 2\lambda_t\operatorname{div}v_tI_d + 2\mu_{tt}\mathbb{D}v + \lambda_{tt}\operatorname{div}vI_d\} $$

and then

$$ \int_0^T \sigma^2 \|\partial_{tt}\tilde{F}(v) - F_{tt}\|_{L^2(\mathbb{R}^d)}^2 \lesssim |\theta - \theta_0|^2 \|v\|_{X_T}^2 $$

from which, we deduce that

$$ (3.49) \qquad \int_0^T \|\sqrt{\sigma}\partial_{tt}\Phi(v), \sigma\nabla\partial_{tt}\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 + \|\sqrt{\sigma}\partial_t\nabla\Phi(v), \sigma\partial_{tt}\Phi(v)\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^2 \lesssim \|(F, u_0)\|_{Y_T^\theta}^2 + |\theta - \theta_0|^2 \|v\|_{X_T}^2. $$

This ends the computations of the higher Hoff estimates. Now, we may turn to the $L^{8/d}((0,T), L^4(\mathbb{R}^d))$ norm estimate for $\sqrt{\sigma}\nabla\partial_t\Phi(v)$. Thanks to (3.9), one has the following:

$$ \int_0^T \sigma^{4/d} \|\nabla \partial_t \Phi(v)\|_{L^4(\mathbb{R}^d)}^{8/d} \lesssim \int_0^T \sigma^{4/d} \|\partial_t \tilde{F}(v)\|_{L^4(\mathbb{R}^d)}^{8/d} + \left[ \int_0^T \|\nabla \Phi(v)\|_{L^8(\mathbb{R}^d)}^{16} \right]^{1/(2d)} + \|\sigma \partial_{tt} \Phi(v)\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^{8/d} \\ + \|\nabla \Phi(v)\|_{L^4((0,T) \times \mathbb{R}^d)}^{8/d} + \left[ \int_0^T \|\nabla \partial_t \Phi(v), \partial_t \tilde{F}(v)\|_{L^2(\mathbb{R}^d)}^2 \right]^{4/d}. $$

The higher estimate (3.49) and (3.48), (3.47) and (3.44) control all of the terms that appear in the left hand side above, except the norm of $\sqrt{\sigma}\partial_t\tilde{F}(v)$ in $L^{8/d}((0,T), L^4(\mathbb{R}^d))$. Starting from (3.43), one has:

$$ \int_0^T \sigma^{4/d} \|\partial_t \tilde{F}(v) - F_t\|_{L^4(\mathbb{R}^d)}^{8/d} \lesssim |\theta - \theta_0|^{8/d} \left\{ \int_0^T \sigma^{4/d} \|\nabla v_t\|_{L^4(\mathbb{R}^d)}^{8/d} + \left[ \int_0^T \|\nabla v\|_{L^8(\mathbb{R}^d)}^{16} \right]^{1/(2d)} \right\} \\ \lesssim |\theta - \theta_0|^{8/d} \|v\|_{X_T}^{8/d}. $$

Gathering all of the estimates one obtains the following:

$$ \int_0^T \sigma^{4/d} \|\partial_t \nabla \Phi(v)\|_{L^4(\mathbb{R}^d)}^{8/d} \lesssim \|(F, u_0)\|_{Y_T^\theta}^{8/d} + |\theta - \theta_0|^{8/d} \|v\|_{X_T}^{8/d}. $$

We now turn to the last step:
For $d=2$, we have:

$$ 
\begin{aligned}
\int_0^T \sigma^\alpha \|u_t\|_{L^{2/(1-\alpha)}(\mathbb{R}^2)}^4 &\le C \int_0^T \sigma^\alpha \|\nabla u_t\|_{L^2(\mathbb{R}^2)}^{4\alpha} \|u_t\|_{L^2(\mathbb{R}^2)}^{4(1-\alpha)} \\
&\le C \left[ \int_0^T \sigma \|\nabla u_t\|_{L^2(\mathbb{R}^2)}^4 \right]^\alpha \left[ \int_0^T \|u_t\|_{L^2(\mathbb{R}^2)}^4 \right]^{1-\alpha}
\end{aligned}
$$

whereas for $d=3$, we only have for all $0 < \alpha < 1/2$

$$ 
\begin{aligned}
\int_0^T \sigma^{2\alpha} \|u_t\|_{L^{3/(1-\alpha)}(\mathbb{R}^3)}^4 &\le C \int_0^T \sigma^{2\alpha} \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{2(1+2\alpha)} \|u_t\|_{L^2(\mathbb{R}^3)}^{2(1-2\alpha)} \\
&\le C \sup_{[0,T]} \|u_t\|_{L^2(\mathbb{R}^3)}^{2(1-2\alpha)} \int_0^T \sigma^{2\alpha} \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{2(1+2\alpha)}.
\end{aligned}
$$

For higher regularity, $1/2 \le \alpha < 1$, we have:

$$ 
\begin{aligned}
\int_0^T \sigma^{\frac{2}{3}(5\alpha-1)} \|u_t\|_{L^{3/(1-\alpha)}(\mathbb{R}^3)}^4 &\le C \int_0^T \sigma^{\frac{2}{3}(5\alpha-1)} \|\nabla u_t\|_{L^4(\mathbb{R}^3)}^{\frac{16}{3}(\alpha-\frac{1}{2})} \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^{\frac{4}{3}(5-4\alpha)} \\
&\le C \left[ \|\sqrt{\sigma}\nabla u_t\|_{L^\infty((0,T),L^2(\mathbb{R}^3))}^2 \right]^{\frac{2}{3}(2-\alpha)} \\
&\quad \times \left[ \int_0^T \sigma^2 \|\nabla u_t\|_{L^4(\mathbb{R}^3)}^{8/3} \right]^{2\alpha-1} \left[ \int_0^T \|\nabla u_t\|_{L^2(\mathbb{R}^3)}^2 \right]^{2(1-\alpha)}.
\end{aligned}
$$

In conclusion, choosing $r$ as:

$$ (3.46) \qquad r = \begin{cases} \max(1/3, 2\alpha) & \text{if } d=2, \\ \max\left(\frac{9}{14}, 2\alpha\right) & \text{if } d=3, \alpha \in (0, 1/2) \\ \frac{2}{3}(5\alpha-1) & \text{if } d=3, \alpha \in [1/2, 1). \end{cases} $$

we have, owing to Young's inequality:

$$ (3.47) \qquad \int_0^T \sigma^r \|(-\Delta)^{-1} \nabla \mathcal{P}(\rho_0 u_t), (-\Delta)^{-1} \nabla (\rho_0 u_t)\|_{\mathscr{C}^\alpha(\mathbb{R}^d)}^4 \le C_* \|F, \nabla u, u_t, \sqrt{\sigma} \nabla u_t\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^4 \\ + C_* \left[ \int_0^T \|F, \nabla u, u_t, \nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \right]^2 + C_* \left[ \int_0^T \sigma^2 \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \right]^{d/2}. $$

The second term is discontinuous at the interface, and we have:

$$ \left\| \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div}(\rho_0 u_t) \right] \right\|_{\mathscr{C}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} \le C \left( \|\lambda\|_{\mathscr{C}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \ell_{\varphi_0}^{-\frac{\alpha}{2}} \right) \|(-\Delta)^{-1} \nabla (\rho_0 u_t)\|_{L^q(\mathbb{R}^d)} \\ + (1 + \|\lambda\|_{\mathscr{C}_{\text{pw},c_0}^{\alpha}(\mathbb{R}^d)} + \mathfrak{P}c_0) \|(-\Delta)^{-1} \nabla (\rho_0 u_t)\|_{\mathscr{C}^{\alpha}(\mathbb{R}^d)}. $$

Finally, we interpolate the $L^q(\mathbb{R}^d)$ norm above between $L^2(\mathbb{R}^d)$ and $L^\infty(\mathbb{R}^d)$, and we recall (3.47) to obtain:

$$ (3.48) \qquad \int_0^T \sigma^r \left\| \frac{1}{\tilde{\mu}} (-\Delta)^{-1} \nabla \mathcal{P}(\rho_0 u_t), \nabla^2 (-\Delta)^{-1} \left[ \frac{1}{2\tilde{\mu} + \lambda} (-\Delta)^{-1} \operatorname{div}(\rho_0 u_t) \right] \right\|^4 \\
\le C_{**} \|F, \nabla u, u_t, \sqrt{\sigma} \nabla u_t\|_{L^\infty((0,T), L^2(\mathbb{R}^d))}^4 \\
+ C_{**} \left[ \int_0^T \|F, \nabla u, u_t, \nabla u_t\|_{L^2(\mathbb{R}^d)}^2 \right]^2 + C_{**} \left[ \int_0^T \sigma^2 \|\nabla u_t\|_{L^4(\mathbb{R}^d)}^{8/d} \right]^{d/2}. $$
Step 3: Piecewise Hölder norm estimate for $\nabla\Phi(v)$. As above, we begin by recalling that, thanks to (3.10), the following Hölder estimate holds for $\nabla\Phi(v)$:

$$ \int_0^T \sigma^r \|\nabla\Phi(v)\|_{\mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)}^4 \lesssim (1 + \mathcal{P}(C)) \int_0^T \sigma^r \|\tilde{F}(v)\|_{\mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)}^4 + \left[ \int_0^T \sigma^{4/d} \|\partial_t \nabla\Phi(v)\|_{L^4(\mathbb{R}^d)}^{8/d} \right]^{d/2} \\ + \|\tilde{F}(v), \nabla\Phi(v), \partial_t\Phi(v), \sqrt{\sigma}\nabla\partial_t\Phi(v)\|_{L^\infty((0,T),L^2(\mathbb{R}^d))}^4 \\ + \left[ \int_0^T \|\tilde{F}(v), \nabla\Phi(v), \partial_t\Phi(v), \partial_t\nabla\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 \right]^2. $$

According to the above computations, one only needs to estimate the first term of the right hand side above and from (3.40) one has:

$$ \int_0^T \sigma^r \|\tilde{F}(v) - F\|_{\mathcal{C}_{pw,\gamma}^{\alpha}(\mathbb{R}^d)}^4 \lesssim |\theta - \theta_0|^4 \|v\|_{X_T}^4. $$

Gathering all of these estimates, one has the following:

$$ \|\Phi(v)\|_{X_T} \le C(R_0 + |\theta - \theta_0|\|v\|_{X_T}) $$

where $R_0$ depends on the norm of $(F, u_0) \in Y_T^\theta$ and the constant $C$ that depends not only on the lower upper bounds of the initial density and viscosity but also all of the norm of the viscosity $\mu$, $\lambda$ and initial density $\rho_0$ as required in (3.2), as well on $\mathcal{P}(C)$. So, choosing large $R_T$, for instance $R_T = 2CR_0$ and $|\theta - \theta_0| \le \min\{1, 1/(2C)\}$, one obtains $\|\Phi(v)\|_{X_T} \le R_T$ as soon as $\|v\|_{X_T} \le R_T$. We thus define an application $\Phi$ from the close ball $B_{X_T(u_0)}(0, R_T)$ of $X_T(u_0)$ to itself. It now remains to prove that $\Phi$ is a contraction for some small $|\theta - \theta_0|$. For this purpose, we consider $v, w \in B_{X_T(u_0)}(0, R_T)$, $\delta v = u - w$ and $\delta\Phi(v) = \Phi(v) - \Phi(w)$. It is obvious that $\delta\Phi(v)$ solves the following Cauchy problem:

$$ \left\{ \begin{array}{l} \rho_0 \partial_t \delta \Phi(v) - \operatorname{div}\{2\mu^{\theta_0} \mathbb{D} \delta \Phi(v) + \lambda^{\theta_0} \operatorname{div} \delta \Phi(v)\} = \operatorname{div}\{\tilde{F}(v) - \tilde{F}(w)\}, \\ \delta \Phi(v)|_{t=0} = 0 \end{array} \right. $$

and moreover, one has, $\partial_t\Phi(v)|_{t=0} = 0$. We also recall that

$$ (3.50) \qquad \delta \tilde{F}(v) := \tilde{F}(v) - \tilde{F}(w) = (\theta - \theta_0)\{2(\mu - \tilde{\mu})\mathbb{D}\delta v + (\lambda - \tilde{\lambda})\operatorname{div}\delta v I_d\}. $$

In the following lines, we will repeat the above steps in order to express the norm of $\delta\Phi(v)$ in $X_T$.

Step 4: Estimates for $\delta\Phi(v)$ and $\partial_t\delta\Phi(v)$. In a similar way as in Step 1, one has:

$$ \begin{align} & \sup_{[0,T]} \|\delta\Phi(v), \partial_t\delta\Phi(v), \nabla\delta\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\nabla\delta\Phi(v), \nabla\partial_t\delta\Phi(v), \partial_t\delta\Phi(v)\|_{L^2(\mathbb{R}^d)}^2 \nonumber \\ & \lesssim \sup_{[0,T]} \|\delta\tilde{F}(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\partial_t\delta\tilde{F}(v)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^T \|\delta\tilde{F}(v)\|_{L^2 \cap L^6(\mathbb{R}^d)}^2 \nonumber \\ & \lesssim |\theta - \theta_0|^2 \|\delta v\|_{X_T}^2 \tag{3.51} \end{align} $$

where we have used the expression (3.50) in order to link the norm of $\delta\tilde{F}(v)$ with the norm of $\delta v$. In the similar way as in step 1, a priori estimates give us:

$$ \|\nabla\delta\Phi(v)\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \lesssim \|\delta\tilde{F}(v)\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} + \sup_{[0,T]} \|\delta\tilde{F}(v), \partial_t\delta\Phi(v), \nabla\delta\Phi(v)\|_{L^2(\mathbb{R}^d)}. $$

Next, we use (3.50) in order to estimate the norm of $\delta\tilde{F}(v)$ and from (3.51) we obtain:

$$ (3.52) \qquad \|\nabla\delta\Phi(v)\|_{L^\infty((0,T),L^6(\mathbb{R}^d))} \lesssim |\theta - \theta_0| \|\delta v\|_{X_T}. $$

On the other hand, for the $L^{16}((0,T), L^8(\mathbb{R}^d))$ norm estimate of $\nabla\delta\Phi(v)$, one has:

$$ \int_0^T \|\nabla\delta\Phi(v)\|_{L^8(\mathbb{R}^d)}^{16} \lesssim \int_0^T \|\delta\tilde{F}(v)\|_{L^8(\mathbb{R}^d)}^{16} + \|\delta\tilde{F}(v), \partial_t\delta\Phi(v), \nabla\delta\Phi(v)\|_{L^2 \cap L^\infty((0,T), L^2(\mathbb{R}^d))}^{16} + \|\nabla\partial_t\delta\Phi(v)\|_{L^2((0,T) \times \mathbb{R}^d)}^{16}. $$

We now gather (3.51) and compute the norm of $\delta\tilde{F}(v)$ from (3.50) just as in the first step in order to obtain:

$$ (3.53) \qquad \int_0^T \|\nabla\delta\Phi(v)\|_{L^8(\mathbb{R}^d)}^{16} \lesssim |\theta - \theta_0|^p \|\delta v\|_{X_T}^{16}. $$