\documentclass[JEP,XML,SOM,Unicode]{cedram}
\datereceived{2016-12-01}
\dateaccepted{2017-06-09}
\dateepreuves{2017-06-15}

\multlinegap0pt
\def\mto{\mathchoice{\longmapsto}{\mapsto}{\mapsto}{\mapsto}}
\def\hto{\mathrel{\lhook\joinrel\to}}
\newcommand{\Psfrac}[2]{(\sfrac{#1}{#2})}
\newcommand{\psfrac}[2]{\sfrac{(#1)}{#2}}
\newcommand{\spfrac}[2]{\sfrac{#1}{(#2)}}
\newcommand{\ppfrac}[2]{\sfrac{(#1)}{(#2)}}

\usepackage{mathrsfs}
\usepackage{dsfont}

\numberwithin{equation}{section}
\renewcommand{\labelenumi}{(\roman{enumi})}

\theoremstyle{plain}
\newtheorem{thm}{Theorem}[section]
\newtheorem{prop}[thm]{Proposition}
\newtheorem{lemma}[thm]{Lemma}
\newtheorem{cor}[thm]{Corollary}

\theoremstyle{definition}
\newtheorem*{Def}{Definition}
\newtheorem{rmk}[thm]{Remark}
\newtheorem*{Rmk}{Remark}
\newtheorem{emp}[thm]{Example}

\def\R{\mathbb R}
\def\N{\mathbb N}
\def\Z{\mathbb Z}
\let\d\partial
\let\ep\varepsilon
\let\wt\widetilde
\DeclareMathOperator{\divv}{div}\let\div\divv

\def\cC{\mathcal C}
\def\cE{\mathcal E}
\def\sC{\mathscr C}
\def\sD{\mathscr D}
\def\sE{\mathscr E}
\def\cF{\mathcal F}
\def\cD{\mathcal D}
\def\sH{\mathscr H}
\def\sK{\mathscr K}
\def\cM{\mathcal M}
\def\cP{\mathcal P}
\def\cS{\mathcal S}
\def\cT{\mathcal T}
\def\cX{\mathcal X}

\def\hB{\dot{B}}
\def\hD{\dot{D}}
\def\hE{\dot{E}}
\def\hF{\dot{F}}
\def\hR{\dot{R}}
\def\hS{\dot{S}}
\def\hT{\dot{T}}
\def\hcT{\dot{\mathcal T}}
\def\hsC{\dot{\mathscr C}}
\def\hsE{\dot{\mathscr E}}
\def\tsC{\widetilde{\mathscr C}}
\def\hDelta{\dot{\Delta}}

\def\tf{\widetilde{f}}
\def\tX{\widetilde{X}}

\def\tDelta{\widetilde{\Delta}}
\def\thDelta{\widetilde{\dot{\Delta}}}

\begin{document}
\frontmatter
\title[Density patches in the inhomogeneous Navier-Stokes equations]{On the persistence of\\ Hölder regular patches of density for\\ the inhomogeneous Navier-Stokes equations}

\author[\initial{R.} \lastname{Danchin}]{\firstname{Raphaël} \lastname{Danchin}}
\address{Université Paris-Est, LAMA (UMR 8050), UPEMLV, UPEC, CNRS\\
61 avenue du Général de Gaulle, 94010 Créteil Cedex 10}
\email{danchin@u-pec.fr}
\urladdr{http://perso-math.univ-mlv.fr/users/danchin.raphael/}

\author[\initial{X.} \lastname{Zhang}]{\firstname{Xin} \lastname{Zhang}}
\address{Université Paris-Est, LAMA (UMR 8050), UPEMLV, UPEC, CNRS\\
61 avenue du Général de Gaulle, 94010 Créteil Cedex 10}
\email{xin.zhang@univ-paris-est.fr}

\thanks{R.D.\ is partially supported by ANR-15-CE40-0011. X.Z.\ is granted by the \emph{Réseau de Recherche Doctoral de Mathématiques de l'Île de France} (RDM-IdF)}

\keywords{Inhomogeneous Navier-Stokes equations, $\mathcal{C}^{1, \varepsilon}$ density patch, striated regularity}

\subjclass{35K59, 76D05}

\begin{abstract}
In our recent work dedicated to the Boussinesq equations \cite{DanZhX2017}, we established the persistence of solutions with piecewise constant temperature along interfaces with Hölder regularity. We here address the same question for the inhomogeneous Navier-Stokes equations satisfied by a viscous incompressible and inhomogeneous fluid. We prove that, indeed, in the slightly inhomogeneous case, patches of densities with $\mathcal{C}^{1, \varepsilon}$ regularity propagate for all time.
Our result follows from the conservation of Hölder regularity along vector fields moving with the flow. The proof of that latter result is based on commutator estimates involving para-vector fields, and multiplier spaces. The overall analysis is more complicated than in \cite{DanZhX2017}, since the coupling between the mass and velocity equations in the inhomogeneous Navier-Stokes equations is \emph{quasilinear} while it is linear for the Boussinesq equations.
\end{abstract}

\alttitle{Persistance de la régularité höldérienne des poches de densité pour les équations de Navier-Stokes inhomogène}

\altkeywords{Équations de Navier-Stokes inhomogène, poches de densité, régularité $\mathcal{C}^{1, \varepsilon}$, régularité stratifiée}

\begin{altabstract}
Dans notre travail récent consacré aux équations de Boussinesq \cite{DanZhX2017}, on a établi la persistance de solutions avec température constante par morceaux le long d'interfaces à régularité höldérienne. On aborde ici la même question pour les équations de Navier-Stokes inhomogène satisfaites par un liquide visqueux incompressible à densité variable. On démontre que, dans le cas légèrement non homogène, les poches de densité avec régularité $\mathcal{C}^{1, \varepsilon}$ se propagent pour tout temps.
Notre résultat est conséquence de la conservation de la régularité höldérienne le long des champs de vecteurs transportés par le flot de la solution. La preuve de ce dernier résultat repose sur des estimations de commutateur mettant en jeu des para-champs et des espaces de multiplicateurs. L'analyse est plus compliquée que dans \cite{DanZhX2017}, dans la mesure où le couplage entre les équations de la masse et de la vitesse dans les équations de Navier-Stokes inhomogène est quasilinéaire alors qu'il est linéaire pour les équations de Boussinesq.
\end{altabstract}
\maketitle

\enlargethispage{.5\baselineskip}%
\vspace*{-1.5\baselineskip}%
\tableofcontents
\mainmatter

\section*{Introduction}
We are concerned with the following \emph{inhomogeneous incompressible Navier-Stokes equations} in the whole space $\R^N$ with $N\geq2$:
\begin{equation}\tag{{INS}}\label{eq:INSmu}
\begin{cases}
\d_t\rho+u\cdot\nabla\rho=0,\\
\rho(\d_tu+u\cdot\nabla u)-\mu\Delta u+\nabla P=0,\\
\div u=0,\\
(\rho,u)|_{t=0} = (\rho_0, u_0).
\end{cases}
\end{equation}
Above, the unknowns
$(\rho, u, P) \in \R_+ \times \R^{N} \times \R$ stand for the density, velocity vector field and pressure, respectively, and the so-called viscosity coefficient $\mu$ is a positive constant.

A number of recent works have been dedicated to the mathematical analysis of the system \eqref{eq:INSmu}. In particular, it is well-known that
if $\rho_0$ is positive and bounded, and $\sqrt{\rho_0}\, u_0$ is in $L^2(\R^N)$,
then the system \eqref{eq:INSmu} admits a global weak solution with finite energy (see \cite{AKM1990}
and the references therein). That result has been extended by {J.\,Simon} in \cite{Sim1990} to the case $\rho_0\geq0$, and
by {P.-L.\,Lions} in \cite{LionsP1996} to viscosity coefficients depending on the density.

In \cite{LionsP1996}, {P.-L.\,Lions} raised the so-called \emph{density patch problem}.
It may be stated as follows: assume that $\rho_0 =\mathds{1}_{\mathcal{D}_0}$ for some domain $\mathcal{D}_0$. Can we find conditions on $u_0$ that ensure that
\begin{equation}\label{eq:rhott}
\rho(t)=\mathds{1}_{\mathcal{D}_t}\quad\hbox{for all }\ t\geq0
\end{equation}
for some domain $\mathcal{D}_t$
\emph{with the same regularity as the initial one}~?

Whenever $\sqrt{\rho_0}\, u_0$ is in $L^2(\R^N)$,
the renormalized solutions theory in \cite{DPL1989} by {R.\,DiPerna and P.-L.\,Lions} for transport equations
ensures that the global weak solutions mentioned above have a
volume preserving generalized flow $\psi$ and that we do have \eqref{eq:rhott} with
$\mathcal{D}_t$ being the image of $\cD_0$ by $\psi(t,\cdot)$.
However, without assuming more on $u_0$, it is very unlikely that one can get information on
the persistence of regularity of $\cD_t$ for positive times.

The present paper aims at making one more step toward solving {Lions'} question, by considering the case where
\begin{equation}\label{eq:rho0}
\rho_0 = \eta_1 \mathds{1}_{\mathcal{D}_0}
+ \eta_2 \mathds{1}_{\mathcal{D}_0^c},
\end{equation}
for some simply connected bounded domain $\mathcal{D}_0$ of
class $\cC^{1,\ep}$.

Our goal is to find as general as possible conditions on $u_0$, that guarantee that for all time $t\geq0$, the domain
$\cD_t:= \psi(t,\cD_0)$ remains $\cC^{1,\ep}$, and
the density reads
\begin{equation}\label{eq:rhot}
\rho(t,\cdot) = \eta_1 \mathds{1}_{\mathcal{D}_t} + \eta_2 \mathds{1}_{\mathcal{D}_t^c}.
\end{equation}

Several recent works give a partial answer to that issue if $|\eta_1-\eta_2|$ is small enough.
Indeed, the paper by the first author with {P.\,B.\,Mucha} \cite{DanM2012}
ensures that if $\rho_0$ is given by \eqref{eq:rho0} and $u_0$ belongs to the critical Besov space $\hB^{\Psfrac{N}{p} -1}_{p,1}(\R^N)$
(see the definition below in \eqref{def:besov}), then \eqref{eq:rhot} is fulfilled
for small time (and for all time if $u_0$ is small) and the $\cC^1$ regularity is preserved.
Likewise, according to the work \cite{HPZ2013}
by {J.\,Huang, M.\,Paicu and P.\,Zhang} (see also \cite{DanZhP2014}), one
may solve \eqref{eq:INSmu} if the initial density is close (for the $L^\infty$ norm) to some positive constant and $u_0$ belongs to
\[
\hB^{\Psfrac{N}{p} -1}_{p,r}(\R^N)\cap \hB^{\Psfrac{N}{p} + \delta-1}_{p,r}(\R^N)
\]
for some
$1<r<\infty$.
As the flow of the corresponding solution is $\cC^{1,\delta'}$ for all $\delta'<\delta$, one can deduce that
if $\rho_0$ is given by \eqref{eq:rho0} then the $\cC^{1,\ep}$ regularity of the boundary
is preserved \emph{provided that $\ep<\delta$} (as the flow need not be in
$\cC^{1,\delta}$).

Finally, as noticed in \cite{DanM2013} then improved by {M.\,Paicu, P.\,Zhang and Z.\,Zhang} \cite{PZZ2013},
in the 2D case, if working \emph{within the energy framework}, then one may avoid the smallness condition on the density
and solve \eqref{eq:INSmu} globally if $\rho_0$ and $u_0$ just satisfy
\begin{equation}\label{cdt:u0PZZ}
0< \eta_1 \leq \rho_0 \leq \eta_2, \quad
u_0 \in H^{s}(\mathbb{R}^2)
\quad \hbox{for some}\,\,\, s>0.
\end{equation}
As the constructed velocity field therein admits a $\cC^1$ flow, one can readily deduce that,
if $\rho_0$ is given by \eqref{eq:rho0} with $\mathcal{D}_0 \subset \mathbb{R}^2$ then the $\cC^1$ regularity of the boundary is preserved.

The common point between the above works is that the hypotheses on $u_0$
do not take into account the non-isotropic structure of $\rho_0$. Consequently, the maximal regularity that can be propagated
for the patch is limited by the overall regularity of the initial velocity.
In two recent papers \cite{LiaoZh2016a,LiaoZh2016b} devoted to the 2D case (see also \cite{LiaoL2016} for the 3-D case),
{X.\,Liao and P.\,Zhang} pointed out that \emph{only tangential regularity along the boundary of $\cD_0$} was needed
to propagate high Sobolev regularity of the patch.
They followed {J.-Y.\,Chemin}'s approach in his work \cite{Che1993} dedicated to the vortex patches problem for the 2-D incompressible Euler equations, and characterized the regularity of the boundary of the domain by means of one (or several) tangent vector fields that
evolve according to the flow of the velocity field.

More precisely, assume with no loss of generality that $\partial\cD_0$ coincides with the level set $f_0^{-1}(\{0\})$
of some (at least $\cC^1$) function $f_0:\R^N\to\R$ that does not degenerate in a neighborhood of $\d\cD_0$, namely there exists some open neighborhood
$V_0$ of $\cD_0$ such that
\begin{equation}\label{eq:f0}
\cD_0=f_0^{-1}(\{0\})\cap V_0\quad\hbox{and }\
\nabla f_0\ \hbox{ does not vanish on }\ V_0.
\end{equation}
Then ${\mathcal{D}_t}$ coincides with $f_t^{-1}(\{0\})$, where $f_t \equiv f(t,\cdot):=f_0\circ \psi^{-1}_t$ with $\psi_t:=\psi(t,\cdot)$ and $\psi$ being the solution of the (integrated) ordinary differential equation:
\begin{equation}\label{eq:flow}
\psi(t,x) =x+\int_0^t u\big(\tau,\psi(\tau,x)\big)\,d\tau.
\end{equation}
Now, the tangent vector field $X_t:=\nabla^\perp f_t$
coincides with the evolution of $X_0:=\nabla^\perp f_0$ along the flow of $u$, namely:\footnote{For any vector field $Y = Y^{k}(x) \d_k$ and function $f$ in $\cC^1(\R^N;\R)$,
we denote by $\d_Y f$ the \emph{directional derivative} of $f$ along $Y$, that is,
with the Einstein summation convention, $\d_Y f:= Y^k \d_k f=Y\cdot \nabla f$.}
\begin{equation}\label{eq:Xt}
X(t,\cdot):= (\d_{X_0} \psi_t)\circ \psi^{-1}_t,
\end{equation}
and thus satisfies the transport equation
\begin{equation}\label{eq:X}
\begin{cases}
\d_t X +u\cdot \nabla X = \d_{X} u, \\
X|_{t=0}=X_0.
\end{cases}
\end{equation}
Consequently, the problem of persistence of regularity for the patch reduces to that of
the vector field $X$ solution to \eqref{eq:X}.

In their first paper \cite{LiaoZh2016a}, {X.\,Liao and P.\,Zhang} justified that heuristics
in the case where the jump $|\eta_1-\eta_2|$ is small enough, and $u_0 \in \big(W^{1,p}(\R^2)\big)^2$ for $ 2<p<4$.
Their proof was essentially based on weighted $L_p-L_q$ estimates for the velocity and allowed to propagate Sobolev regularity $W^{k,p}$
of the boundary, with $k$ large enough (in particular the boundary
is at least $\cC^{2,\ep}$ for some $\ep>0$).
In a second paper \cite{LiaoZh2016b}, after revisiting the approach of \cite{PZZ2013} (that is Sobolev spaces $H^s$ with $s>0$ and thus finite energy framework),
{X.\,Liao and P.\,Zhang} succeeded in proving a similar result for \emph{general} positive $\eta_1$ and $\eta_2$
in \eqref{eq:rho0}. The corresponding level set function $f_0$ has to be in $W^{2+k,p}(\R^2)$ for
some integer number $k\geq 1$ and $p\in{}]2,4[$, hence $\cD_0$ is still at least $\cC^{2,\ep}$.
As regards the initial velocity field $u_0$, it has to satisfy
the following \emph{striated regularity} property along the vector field $X_0:=\nabla^\perp f_0$:
$$
\big(\d_{X_0}^\ell u_0 \in B^{s+\ep \psfrac{k-\ell}{k} }_{2,1} (\R^2) \big)^2 \hbox{ for all }\ \ell\in\{0,\dots,k\}\
$$
with $0<s<1-\ep$ and $(s,p)$ in $]0,1[{}\times{}]2,\min\{4,2\slash(1-s)\}[$.

In the present paper, we propose a simpler approach that allows
to propagate just $\cC^{1,\ep}$ Hölder regularity (for all $\ep\in{}]0,1[$), within a \emph{critical} regularity framework.
By~critical, we mean that the solution space that we shall consider has the same scaling invariance by time and space dilations as
\eqref{eq:INSmu} itself, namely:
\begin{equation}\label{eq:critical}
(\rho, u, P)(t,x) \to (\rho, \lambda u, \lambda^2 P)(\lambda^2 t , \lambda x)\quad\hbox{and}\quad
(\rho_0, u_0)(x) \to (\rho_0, \lambda u_0)(\lambda x).
\end{equation}
That framework is by now classical
for the homogeneous Navier-Stokes equations (that is $\rho$ is a positive constant in \eqref{eq:INSmu}) in the whole space $\R^N$ (see \eg \cite{BCD2011,L-R2016} and the references therein).
As observed by the first author in \cite{Dan2003} (see also {H.\,Abidi} in \cite{Ab2007} and {H.\,Abidi and M.\,Paicu} in \cite{AbP2007}),
working in a suitable critical functional framework is still relevant in the inhomogeneous situation.

\subsubsection*{Acknowledgements}
We are grateful to the referees for pointing out the work \cite{GG-J2016} by {F.\,Gancedo and E.\,García-Juárez}
that has been posted on arXiv a few weeks after we submitted our paper.
There, in the 2D case, the $\cC^{1,\varepsilon}$ regularity of the interface is propagated
whenever the density is given by \eqref{eq:rho0} with $\eta_1>0$ and $\eta_2>0$,
and~$u_0$ is in~$H^{\ep'}$ for some $\ep'>\ep$.
That result relies on the nice observation that in that particular case,
the flow of the solution constructed in \cite{PZZ2013} has $\cC^{1,\ep}$ regularity
\emph{in all directions}.
In the present paper, we are able to consider densities which are not piecewise constant
(like in \eqref{eq:rho0} with a smooth variable $\eta_1$ for example),
and velocity fields that have just critical regularity
but we need the density to be close to some positive constant.

Part of this work has been performed during a visit of the second author in Peking University. Zhifei Zhang and Chao Wang are warmly thanked for their welcome and stimulating discussions during this period.

\section{Results}

Before stating our main result, we need to introduce a few notations.
First, we recall the definition of Besov spaces: following \cite[\S 2.2]{BCD2011},
we consider two smooth radial functions $\chi$ and $\varphi$ supported in $\{ \xi \in \R^N:|\xi| \leq {4}/{3} \}$ and $\{ \xi \in \R^N: {3}/{4} \leq |\xi| \leq {8}/{3} \}$, respectively, and satisfying
\begin{align} \label{eq:L-P_decomp}
\sum_{j \in \mathbb{Z} } \varphi(2^{-j}\xi) =1,\quad \forall ~\xi \in \R^N \setminus \{0\},\quad
\chi(\xi) + \sum_{j \geq 0 } \varphi(2^{-j}\xi) = 1,\quad\forall ~\xi \in \R^N .
\end{align}
Then we define Fourier truncation operators as follows:
\begin{equation*}
\hDelta_j:= \varphi(2^{-j}D),\,\, \hS_j:= \chi(2^{-j} D),\,\,\forall j\in \Z; \quad\Delta_j:= \varphi(2^{-j}D), \,\,\forall j \geq 0, \quad\Delta_{-1}: = \chi(D).
\end{equation*}
For all triplet $(s,p,r) \in \R \times [1,\infty]^2$,
the homogeneous Besov space $\hB^s_{p,r}(\R^N)$ (just denoted by $\hB^s_{p,r}$ if the
value of the dimension is clear from the context) is defined~by
\begin{equation}\label{def:besov}
\hB^{s}_{p,r}(\R^N):= \big\{u \in \mathcal{S}'_h(\R^N): \|u\|_{\hB^s_{p,r}}:=\big\|2^{js}\|\hDelta_j u \|_{L^p} \big\|_{\ell^r(\mathbb{Z}) } <\infty \big\},
\end{equation}
where $\mathcal{S}'_h(\R^N)$ is the subspace of tempered distributions $\mathcal{S}'(\R^N)$ defined by
\begin{equation*}
\mathcal{S}'_h(\R^N):=\bigl\{u \in \mathcal{S}'(\R^N): \lim_{j\to -\infty} \hS_j u = 0 \,\bigr\}\cdotp
\end{equation*}
We shall also use sometimes the following inhomogeneous Besov spaces:
\begin{equation}\label{def:inhbesov}
B^{s}_{p,r}(\R^N):= \big\{u \in \mathcal{S}'(\R^N): \|u\|_{B^s_{p,r}}:=\big\|2^{js}\|\Delta_j u \|_{L^p} \big\|_{\ell^r(\N \cup \{-1\}) } <\infty \big\}.
\end{equation}
Throughout, we adopt the common notation $b^s_{p,r}(\R^N)$ to denote
$B^{s}_{p,r}(\R^N)$ or $\hB^{s}_{p,r}(\R^N)$.

It is well-known that Sobolev or Hölder spaces
belong to the Besov spaces hierarchy.
For instance $\dot B^s_{2,2}(\R^N)$ coincides with the homogeneous Sobolev space $\dot H^s(\R^N)$,
and we have \begin{equation}\label{spc:Holder}
B^s_{\infty,\infty}(\R^N) = \cC^{0,s}(\R^N) = L^{\infty}(\R^N) \cap \hB^s_{\infty,\infty}(\R^N)\quad \mbox{if} \quad s\in{}]0,1[.
\end{equation}
To emphasize that connection, we shall often use the notation
$\hsC^{s}:= \hB^{s}_{\infty,\infty}$ (or~$\sC^{s}:=\nobreak B^{s}_{\infty,\infty}$) for \emph{any} $s\in \R$.

When investigating evolutionary equations in critical Besov spaces, it is wise to use
the following \emph{tilde homogeneous Besov spaces} first introduced by {J.-Y.\,Chemin} in~\cite{Che1999}:
for any $T\in{}]0,+\infty]$ and $(s,p,r,\gamma) \in \R \times [1,+\infty]^3$, we set\footnote{For $T\in{}]0,+\infty[$, $p\in[1,+\infty]$ and $E$ a Banach space, the notation $L^p_T(E)$
denotes the space of $L^p$ functions on $]0,T[$ with values in $E$, and $L^p(\R_+; E)$
corresponds to the case $T=+\infty$. We keep the same notation
for vector or matrix-valued functions.}
\begin{equation*}
\widetilde{L}^\gamma_T \big(\hB^{s}_{p,r} \big):= \Bigl\{ u \in \mathcal{S}'(]0,T[{}\times \R^N): \lim_{j \to -\infty} \hS_j u =0\ \mbox{in}\ L^\gamma_T(L^\infty)\text{ and }\|u\|_{\widetilde{L}^\gamma_T (\hB^{s}_{p,r})} <\infty\Bigr\},
\end{equation*}
where
\begin{equation*}
\|u\|_{\widetilde{L}^\gamma_T (\hB^{s}_{p,r})}:=\big\|2^{js}\|\hDelta_j u\|_{L^\gamma_T(L^p)} \big\|_{\ell^r(\Z)} <\infty .
\end{equation*}
The index $T$ will be omitted if equal to $+\infty$, and we shall denote
$$
\wt\cC_b(\R_+;\dot B^s_{p,r}):=\wt L^\infty(\R_+;\dot B^s_{p,r})\cap \cC(\R_+;\dot B^s_{p,r}).
$$
We also need to introduce the following spaces
for $(\sigma, p, T) \in \R \times [1,\infty] \times{}]0,\infty] $:
\begin{equation*}
\hE^\sigma_p (T):=\big\{ (v, \nabla Q): v\in \wt\cC_b \big([0,T[; \hB^{\Psfrac{N}{p}-1+\sigma}_{p,1}),
\: (\d_t v, \nabla^2 v , \nabla Q) \in L^1_T \big(\hB^{\Psfrac{N}{p}-1+\sigma}_{p,1}\big) \big\} ,
\end{equation*}
endowed with the norm
\begin{equation*}
\|(v,\nabla Q)\|_{\hE^\sigma_p (T)}:= \|v\|_{\wt{L}^\infty_T(\hB^{\Psfrac{N}{p}+\sigma-1}_{p,1})} + \|(\d_t v, \nabla^2 v, \nabla Q)\|_{L^1_T (\hB^{\Psfrac{N}{p}+\sigma-1}_{p,1})}.
\end{equation*}
For notational simplicity, we shall omit $\sigma$ or $T$ in the notation $\hE^\sigma_p(T)$ whenever $\sigma$ is zero or $T=\infty$. For instance, $\hE_p:= \hE^0_p (\infty)$.

Finally, we shall make use of \emph{multiplier spaces} associated to pairs $(E, F)$ of Banach spaces
included in the set of tempered distributions. The definition goes as follows:
\begin{Def}
Let $E$ and $F$ be two Banach spaces embedded in $\cS'(\R^N)$.
The \emph{multiplier space} $\cM(E \to F)$ (simply denoted by $\cM(E)$ if $E=F$)
is the set of those functions $\varphi$ satisfying $\varphi u\in F$ for all $u$ in $E$ and,
additionally,
\begin{equation}\label{eq:def_norm_multiplier}
\|\varphi\|_{\cM(E \to F)}:= \sup_{ \substack{u \in E \\ \|u\|_{_E} \leq 1} } \|\varphi u\|_F<\infty.
\end{equation}
\end{Def}
It goes without saying that $\|\cdot\|_{\cM(E\to F)}$ is a norm on $\cM(E\to F)$
and that one may restrict the supremum in \eqref{eq:def_norm_multiplier}
to any \emph{dense} subset of $E$.


The following result that has been proved in \cite{DanM2012}
is the starting point of our analysis:\footnote{As the viscosity coefficient $\mu$ will be fixed once and for all, we shall set it to $1$
for notational simplicity. Likewise, we shall assume the
reference density at infinity to be $1$.}
\begin{thm}\label{thm:wellposed_INS}
Let $p\in [1,2N[$ and $u_0$ be a divergence-free vector field with coefficients in $\hB^{\Psfrac{N}{p}-1}_{p,1}$. Assume that $\rho_0$ belongs to the multiplier space $\cM \big(\hB^{\Psfrac{N}{p}-1}_{p,1}\big)$. There exist two constants $c$ and
$C$ depending only on $p$ and on $N$ such that if
\begin{equation*}
\|\rho_0 -1 \|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1})} + \|u_0\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \leq c
\end{equation*}
then the system \eqref{eq:INSmu} in $\R^N$ with $N\geq2$ has a unique solution $(\rho, u, \nabla P)$ satisfying
$$
\rho \in L^{\infty}\Big(\R_+; \cM \big(\hB^{\Psfrac{N}{p}-1}_{p,1} \big)\Big)\quad\hbox{and}\quad (u,\nabla P) \in \dot E_p.$$
Furthermore, the following inequality is fulfilled:
\begin{equation}\label{est:up}
\|u\|_{\wt L^\infty(\R^+;\dot B^{\Psfrac Np-1}_{p,1})}+\|\d_tu,\nabla^2u,\nabla P\|_{L^1(\R^+;\dot B^{\Psfrac Np-1}_{p,1})} \leq C\|u_0\|_{\dot B^{\Psfrac Np-1}_{p,1}}.
\end{equation}
\end{thm}

By classical embedding, having $\nabla^2u$ in $L^1(\R_+;\dot B^{\Psfrac Np-1}_{p,1})$ implies that $\nabla u$ is
in $L^1(\R_+;\cC_b)$. Therefore the flow $\psi$ of $u$ defined by \eqref{eq:flow} is in $\cC^1$.
Now, it has been observed in \cite{DanM2012} that for any uniformly $\cC^1$ bounded domain $\cD_0$,
the function $\mathds{1}_{\cD_0}$ belongs to
$\cM \big(\hB^{s}_{p,1}\big)$ whenever $-1+\sfrac1p<s<\sfrac1p$.
Therefore, one may deduce from Theorem \ref{thm:wellposed_INS} that if
$\rho_0$ is given by \eqref{eq:rho0}, if $u_0$ is in $\dot B^{\Psfrac Np-1}_{p,1}$ for some $N-1<p<2N$
and if
$$\|u_0\|_{\dot B^{\Psfrac Np-1}_{p,1}}+|\eta_2-\eta_1|\quad\hbox{is small enough}$$
then the system \eqref{eq:INSmu} admits a unique global solution $(\rho,u,\nabla P)$ with
$(u,\nabla P)$ in $\dot E_p$ and $\rho$ given by \eqref{eq:rhot} with $\cD_t$
in $\cC^1$ for all time $t\geq0$.

The (parabolic type) gain of regularity for $u$ pointed out in Theorem \ref{thm:wellposed_INS} is optimal, as well as the embedding of $\dot B^{\sfrac Np}_{p,1}$ in the set of continuous bounded functions.
Therefore, one cannot expect the flow of $u$ given by Theorem \ref{thm:wellposed_INS}
to be in any Hölder space~$\cC^{1,\alpha}$ for some $\alpha>0$,
which prevents our propagating more than $\cC^1$ regularity. Following \cite{LiaoZh2016a,LiaoZh2016b},
it is natural to make an additional tangential regularity assumption on $u_0$.
This motivates the following
general result of persistence of geometric structures for \eqref{eq:INSmu}.
\begin{thm}\label{thm:main_dINS}
Let $\ep$ be in $]0,1[$ and $p$ satisfy
\begin{equation}\label{cdt:p}
\sfrac N2<p < \min \big\{ \spfrac{N}{1-\ep}, 2N\big\}\cdotp
\end{equation}
Let $u_0$ be a divergence-free vector field with coefficients in $\hB^{\Psfrac{N}{p}-1}_{p,1}$. Assume that the initial density $\rho_0$ is bounded and belongs to the multiplier space
\[
\cM \big(\hB^{\Psfrac{N}{p}-1}_{p,1}\big) \cap \cM (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}).
\]
There exists a constant $c$ depending only on $p$ and $N$ such that if
\begin{equation}\label{cdt:u0_small}
\|\rho_0 -1 \|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1}) \cap \cM (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})\cap L^\infty} + \|u_0\|_{{\hB}^{\Psfrac{N}{p}-1}_{p,1}} \leq c,
\end{equation}
then the system \eqref{eq:INSmu} in $\R^N$ has a unique global solution $(\rho, u, \nabla P)$ with
$$\rho\in L^\infty\Big(\R_+ ; L^\infty\cap \cM \big(\hB^{\Psfrac{N}{p}-1}_{p,1}\big) \cap \cM (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}) \Big) \quad\hbox{and}\quad (u,\nabla P)\in \hE_p.$$
Moreover, for any vector field $X_0$ with $\cC^{0,\ep}$ regularity
(assuming in addition that $\ep>2-\sfrac Np$ if $\div X_0\not\equiv0$),
if the following conditions are fulfilled
$$\d_{X_0} \rho_0 \in \cM \big(\hB^{\Psfrac{N}{p}-1}_{p,1}\to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1}\big) \quad\hbox{and}\quad
\d_{X_0} u_0\in\hB^{\Psfrac{N}{p}+\ep-2}_{p,1},$$
then the system \eqref{eq:X} in $\R^N$ has a unique global
solution $X\in C_w(\R_+;\cC^{0,\ep})$, and we have
\begin{equation*}
\d_X \rho \in L^\infty \Big(\R_+ ; \cM \big(\hB^{\Psfrac{N}{p}-1}_{p,1}\to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1}\big) \Big)\quad\hbox{and}\quad
(\d_X u, \d_X\nabla P)\in \hE^{\ep-1}_p.
\end{equation*}
\end{thm}
Some comments are in order:
\begin{itemize}
\item The divergence-free property on $X_0$ is conserved during the evolution because if one takes the divergence of \eqref{eq:X}, then we get, remembering that $\div u =0$,
\begin{equation}\label{eq:divX}
\begin{cases}
\d_t \div X + u\cdot \nabla \div X = 0,\\[1ex]
\div X |_{t=0}=\div X_0.
\end{cases}
\end{equation}
\item In the case $\div X_0 \not= 0$, the additional constraint on $(\ep,p)$
is due to the fact that the product of a general $\cC^{0,\ep}$ function with a $\dot B^{\Psfrac Np-2}_{p,1}$
distribution need not be defined if the sum of regularity coefficients, namely $\ep+\Psfrac Np-2$, is negative.
\item The vector field $X$ given by \eqref{eq:X} has
the Finite Propagation Speed Property.
Indeed, from the definitions of the flow and of the space $\dot E_p$, and from the embedding of
$\dot B^{\sfrac Np}_{p,1}(\R^N)$ in $\cC_b(\R^N)$, we readily get\footnote{All over the paper, we agree that $A \lesssim B$ means $A \leq C B$ for some harmless ``constant'' $C$,
the meaning of which may be guessed from the context.}
for all $t\geq0$ and $x\in\R^N$,
$$
\bigl|\psi(t,x)-x\bigr|\lesssim \sqrt t\,\|u\|_{L^2_t(\dot B^{\sfrac Np}_{p,1})}\leq C\sqrt t\,\|u_0\|_{\dot B^{\Psfrac Np-1}_{p,1}}.$$
Therefore, if the initial vector field $X_0$ is supported in the set $K_0$
then $X(t)$ is supported in some set $K_t$ such that
\begin{equation*}
\mathrm{diam}(K_t) \leq \mathrm{diam}(K_0) + C\sqrt t\,\|u_0\|_{\dot B^{\Psfrac Np-1}_{p,1}}.
\end{equation*}
\item
One can prove a similar result (only local in time) for \emph{large} $u_0$ in $\dot B^{\Psfrac Np-1}_{p,1}$. Moreover, we expect our method to be appropriate for handling Hölder regularity $\cC^{k,\ep}$ if making suitable assumptions on $\d_{X_0}^j\rho_0$ and $\d_{X_0}^j u_0$ for $j=0,\dots,k$. We refrained from writing out here this generalization to keep the presentation as short as possible.
\end{itemize}

In the density patch situation (that is if $\rho_0$ is given by \eqref{eq:rho0})
the condition on $\d_{X_0}\rho_0$ is trivially satisfied as the derivative of the density along
any continuous vector field that is tangent to $\partial\cD_0$, vanishes.
This implies the following statement of propagation
of Hölder regularity of density patches for \eqref{eq:INSmu} in the plane:
\begin{thm}\label{thm:N2}
Let $\cD_0$ be a simply connected bounded domain of $\R^2$ satisfying \eqref{eq:f0} for some
function $f_0\in \cC^{1,\ep}(\R^2;\R)$ with $\ep$ in $]0,1[$. There exists a constant $\eta_0$ depending only on $\cD_0$
and such that if \begin{equation}\label{eq:density_patch}
\rho_0:= (1+\eta) \mathds{1}_{\mathcal{D}_0} + \mathds{1}_{\mathcal{D}^c_0}
\quad\hbox{with}\quad \eta\in{}]-\eta_0,\eta_0[
\end{equation}
and if the divergence free vector-field $u_0\in\cS'_h(\R^2)$
has vorticity $\omega_0:=\d_1u_0^2-\d_2u_0^1$ given by
\begin{equation}\label{eq:vorticity}
\omega_0= \overline{\omega}_0 + \wt\omega_0\, \mathds{1}_{\mathcal{D}_0}\ \hbox{ with }\
\div (\overline{\omega}_0 \nabla^{\perp}f_0) =0
\quad \hbox{and} \quad \int_{\R^2} \omega_0 \,dx=0
\end{equation}
for some small enough compactly supported functions $(\overline{\omega}_0, \wt\omega_0)$ in
$L^p(\R^2) \times \cC^{\ep'}(\R^2)$ with $0<\ep'<\ep$ and $1<p< 2\slash (2-\ep)$,
then the system \eqref{eq:INSmu} has a unique solution $(\rho, u,\nabla P)$ with the properties
listed in Theorem \ref{thm:wellposed_INS}.
Moreover, if we denote by $\psi(t,\cdot)$ the flow of $u$ then for all $t\geq0$, we have
\begin{equation}\label{eq:density_patcht}
\rho(t,\cdot):= (1+\eta) \mathds{1}_{\mathcal{D}_t} + \mathds{1}_{\mathcal{D}^c_t}
\ \hbox{ with }\ \mathcal{D}_t:= \psi (t,\mathcal{D}_0),
\end{equation}
and $\mathcal{D}_t$ remains a simply connected bounded domain of class $\cC^{1,\ep}$.
\end{thm}

\begin{rmk}
Of course, one can take $\overline{\omega}_0 \equiv 0$ or $\wt\omega_0 \equiv 0$ in the above statement.

The zero average condition guarantees that $u_0$ belongs to some homogeneous Besov space
$\dot B^{\Psfrac2p-1}_{p,1}$ (so as to apply Theorem \ref{thm:main_dINS}).
It is not needed if we have constant vortex pattern in $\R^2$ (see Theorem \ref{thm:N22}) or if the dimension $N=3$ (see Theorem \ref{thm:N3}).
\end{rmk}

\begin{rmk}
We imposed the particular structure of the vorticity in Theorem \ref{thm:N2} just to give an
explicit example for which \eqref{eq:rhot} with regularity $\cC^{1,\ep}$ holds true. It goes without saying
that one can consider a much more general class of initial velocities: according to
Theorem \ref{thm:main_dINS}, it suffices that $u_0$ satisfies the smallness condition of Theorem \ref{thm:wellposed_INS}
and that $\div(\nabla^\perp f_0\otimes u_0)$ is
in $\dot B^{\Psfrac2p+\ep-2}_{p,1}$ for some $1<p<\min\{4,2/(1-\ep)\}$.
In other words, we just need ``$u_0$ to have $\ep$ more regularity in
the direction that is tangential to the patch of density.''
This is of course satisfied if $u_0$ vanishes in a neighborhood of~$\cD_0$.
However, one may consider much more singular examples like
the case where $u_0$ is compactly supported and behaves locally near some $x_0 \notin \d \cD_0$,
like the function $|x-x_0|^{-1} (-\log |x-x_0|)^{-(1+\delta)}$ with $\delta>0$.
\end{rmk}
\begin{rmk}
Similar results, only local in time, hold true for large $u_0$ with critical regularity.
\end{rmk}
We end this section with a short presentation of the main ideas of the proof of Theorem \ref{thm:main_dINS}.
{}From Theorem \ref{thm:wellposed_INS}, we have
a global solution $(\rho,u,\nabla P)$ such that $\rho\in L^\infty\bigl(\R_+ ; \cM \big(\hB^{\Psfrac{N}{p}-1}_{p,1}\big)\bigr)$ and $(u,\nabla P)\in\hE_p$.
As already explained, our main task is to prove that $X(t,\cdot)$ remains in $\cC^{0,\ep}$ for all time. Now,
in light of \eqref{eq:X}, we have
\begin{equation*}
X(t,x)= X_0\big(\psi^{-1}_t(x)\big) + \int^t_0 \d_X u \big(t', \psi_{t'} \big(\psi^{-1}_t(x) \big)\big)dt'.
\end{equation*}
Therefore, because $\psi_t$ is a $\cC^1$ diffeomorphism of $\R^N$, it suffices
to show that $\d_Xu$ is in $L^1_\mathrm{loc}(\R_+;\cC^{0,\ep})$. To this end, it is natural to look for a suitable evolution
equation for $\d_Xu$. Since \eqref{eq:X} means that $[D_t, \d_X]=0$, where $D_t: = \d_t + u \cdot \nabla$ stands for the material derivative associated to~$u$,
differentiating the momentum equation of \eqref{eq:INSmu} along $X$ yields
\begin{equation}\label{eq:dXu}
\rho D_t \d_X u + \d_X \rho D_t u- \d_X \Delta u + \d_X \nabla P =0.
\end{equation}
Even though \eqref{eq:dXu} has some similarities with the Stokes system,
it is not clear that it does have the same smoothing properties, as its
coefficients have very low regularity.
One of the difficulties lies in the product of
the discontinuous function $\rho$ with $D_t\d_Xu$, as having
only $\d_X u$ in $\cC^{0,\ep}$ suggests that $D_t\d_Xu$ has \emph{negative} regularity.
At the same time, the term with $\d_X\rho$ is harmless as, owing to $[D_t, \d_X]=0$ and to the mass equation, we have
\begin{equation}\label{eq:dXrho}D_t\d_X\rho=0.\end{equation}
Hence any (reasonable) regularity assumption for $\d_{X_0}\rho_0$ persists through the evolution.

Our strategy is to assume that $\rho$ belongs to the multiplier space corresponding to the space
to which $D_t\d_X u$ is expected to belong. As the flow is $\cC^1$,
propagating multiplier information is straightforward (see Lemma \ref{lemma:multiplier}).
This new viewpoint spares us the tricky energy estimates
and iterated differentiation along vector fields (requiring higher regularity of the patch)
that were the cornerstone of the work by {X.\,Liao and P.\,Zhang}.
In fact, \emph{under the smallness assumption \eqref{cdt:u0_small}}
which, unfortunately, forces the fluid to have small density variations, we succeed
in closing the estimates using \emph{only one} differentiation along $X$.
This makes the proof rather elementary and allows us to propagate low Hölder regularity.

Whether one can differentiate terms like $\Delta u$ or $\nabla P$
along $X$ within our critical regularity framework is not totally clear, though. Therefore,
as in our recent work~\cite{DanZhX2017}
dedicated to the incompressible Boussinesq system, we shall replace differentiation along vector-fields by \emph{para-differentiation}.

Let us briefly recall how it works. Fix some suitably large integer $N_0$ and
introduce the following \emph{paraproduct} and \emph{remainder} operators (after {J.-M. Bony} in \cite{Bo1981}):
$$\hT_u v:= \sum_{j\in \Z} \hS_{j-N_0}u \hDelta_j v\quad\hbox{and}\quad
\hR(u,v) \equiv \sum_{j\in\Z} \hDelta_{j} u \thDelta_{j} v:= \sum_{\substack{j\in\Z\\|j-k|\leq N_0}} \hDelta_{j} u \hDelta_{k} v.
$$
Then any product may be formally decomposed as follows:
\begin{equation}\label{eq:bony}
uv = \hT_u v + \hT_v u + \hR(u,v).
\end{equation}
To overcome the problem with the definition of $\d_X\Delta u$ and $\d_X\nabla P$, we shall
change the vector field $X$ to the
para-vector field operator $\hcT_X \cdot:= \dot{T}_{X^k} \d_k \cdot$ which, in our regularity framework, will turn out to be
\emph{the principal part} of operator $\d_X$. Indeed, for any pair $(X,f)$, the decomposition \eqref{eq:bony} ensures that
\begin{equation*}
(\hcT_X -\d_X) f = \hT_{\d_k f} X^k + \d_k \hR(f, X^k) - \hR(f,\div X).
\end{equation*}
Therefore, taking advantage of classical continuity results for operators $\hT$ and $\hR$ (see~\cite{BCD2011}), we discover that
\begin{equation}\label{es:dXf_TXf}
\| (\hcT_X -\d_X)f \|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|f\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \|X\|_{\hsC^{\ep}}
\end{equation}
whenever $(\ep,p)\in{}]0,1[{}\times[1,+\infty]$ fulfills:
\begin{equation}\label{cdt:p_X}
\sfrac Np \in \,]\,1-\ep,\,2 \,[\ \hbox{ if }\ \div X=0, \quad\hbox{and }\ \sfrac Np \in \,]\,2-\ep,\,2\,[
\ \hbox{ otherwise.}
\end{equation}
In our situation, we will apply \eqref{es:dXf_TXf} with $f=\nabla P$ or $\Delta u$, which are in $L^1 \big(\R_+;\dot B^{\Psfrac Np-1}_{p,1}(\R^N)\big)$.

Now, incising the term $\d_X u$ by the scalpel $\hcT_X $ in \eqref{eq:dXu} and applying $\hcT_X $ to the third equation of \eqref{eq:INSmu} yield
\begin{equation} \label{eq:dXINS}
\begin{cases}
\rho D_t \hcT_X u-\Delta \hcT_X u + \nabla \hcT_X P = g,\\
\div \hcT_X u = \div (\hT_{\d_k X} u^k -\hT_{\div X} u),\\
\hcT_X u|_{t=0} = \hcT_{X_0} u_0
\end{cases}
\end{equation}
with
\begin{multline}\label{eq:g_dINS}
g:= - \rho [\hcT_{X}, D_t]u +[\hcT_X,\Delta] u - [\hcT_X, \nabla] P + (\d_X-\hcT_X) (\Delta u - \nabla P)
\\ -\d_X \rho D_t u + \rho (\hcT_X-\d_X)D_t u.
\end{multline}
This surgery leading to \eqref{eq:dXINS} is effective for three reasons.
First, all the commutator terms in \eqref{eq:g_dINS} are under control (see the Appendix). Second,
$D_t\d_Xu$ and $D_t\hcT_Xu$ are in the same Besov space,
and the multiplier type regularity on the density that was pointed out before is thus appropriate. Last, the condition \eqref{cdt:u0_small} ensures that $(\hcT_X -\d_X)u$ is
a remainder term. Of course, the divergence free condition need not be satisfied by $\hcT_X u$, but one can further modify
\eqref{eq:dXINS} so as to enter in the standard maximal regularity theory.
Then, under the smallness condition \eqref{cdt:u0_small}, one can close the estimates involving $\d_Xu$ or $\d_X\rho$, globally in time.


The rest of the paper unfolds as follows.
In the next section, we show that Theorem~\ref{thm:main_dINS} entails a general (but not so explicit) result of persistence of Hölder regularity for patches of density in any dimension.
We shall then obtain Theorem \ref{thm:N2}, and an analogous result in dimension $N=3$.
Section \ref{sect:3} is devoted to the proof of our general result of all-time persistence of striated regularity (Theorem \ref{thm:main_dINS}).
Some technical results pertaining to commutators and multiplier spaces are postponed in appendix.

\section{The density patch problem}

This section is devoted to the proof of results of persistence of regularity for patches of constant densities,
taking Theorem \ref{thm:main_dINS} for granted.
Throughout this section we shall use repeatedly the fact (proved in \eg \cite[Lem.\,A.7]{DanM2012}) that for any (not necessarily bounded) domain $\cD$ of $\R^N$ with uniform $\cC^{1}$ boundary, we have
\begin{equation*}
\mathds{1}_{\mathcal{D}} \in \cM\big(\hB^{s}_{p,r} (\R^N)\big) \quad \mbox{whenever} \quad (s,p,r) \in{}]\Psfrac{1}{p}-1,\sfrac{1}{p}[{}\times{}]1,\infty[{}\times [1,\infty].
\end{equation*}
{}From that property, we deduce that if $(\ep, p) \in{}]0,1[{}\times{}]N-1, \frac{N-1}{1-\ep}[ , $ then the density~$\rho_0$ given by \eqref{eq:density_patch}
belongs to $\cM\big(\hB^{\Psfrac{N}{p}-1}_{p,r}(\R^N)\big) \cap \cM\big(\hB^{\Psfrac{N}{p}+\ep-2}_{p,r}(\R^N)\big)$.


As a start, let us give a result of persistence of regularity,
under rather general hypotheses.
\begin{prop}\label{prop:NN} Assume that $\rho_0$ is given by \eqref{eq:density_patch}
with small enough $\eta$ and
some~$\cC^{1,\ep}$ domain $\cD_0$ of $\R^N$ satisfying \eqref{eq:f0}.
Let $u_0$ be a small enough divergence free vector field with coefficients in
$\dot B^{\Psfrac Np-1}_{p,1}$ for some
\begin{equation}\label{cdt:u_01}
N-1<p<\min\big\{\ppfrac{N-1}{1-\ep},2N\big\}\cdotp\end{equation}
Consider a family $(X_{\lambda,0})_{\lambda\in\Lambda}$
of $\cC^{0,\ep}$ divergence free vector fields tangent
to $\cD_0$ and such that $\d_{X_{\lambda,0}}u_0 \in \dot B^{\Psfrac Np+\ep-2}_{p,1}$ for
all $\lambda\in\Lambda$.

Then the unique solution $(\rho,u,\nabla P)$ of \eqref{eq:INSmu} given by Theorem \ref{thm:wellposed_INS}
satisfies the following additional properties:
\begin{itemize}
\item $\rho(t,\cdot)$ is given by \eqref{eq:density_patcht},
\item all the time-dependent vector fields $X_{\lambda}$ solutions to \eqref{eq:X} with initial data
$X_{\lambda,0}$ are in $L^\infty_\mathrm{loc}(\R_+;\cC^{0,\ep})$ and remain tangent to the patch for all time.
\end{itemize}
\end{prop}
\begin{proof}
Assumptions \eqref{eq:density_patch} and \eqref{cdt:u_01} guarantee that
$\rho_0$ is in
\[
\cM(\hB^{\Psfrac Np-1}_{p,1}) \cap \cM(\hB^{\Psfrac Np-2+\ep}_{p,1}),
\]
and that
\eqref{cdt:u0_small} is fulfilled if $\eta$ and $u_0$ are small enough.
Of course, $\d_{X_{\lambda, 0}}\rho_0\equiv0$ for all $\lambda\in\Lambda$ because
the vector fields $X_{\lambda,0}$ are tangent to the boundary.
Therefore, applying Theorem \ref{thm:main_dINS} ensures that
all the vector fields $X_{\lambda}$ are in $L^\infty_\mathrm{loc}(\R_+;\cC^{0,\ep})$.
Now, if we consider a level set function $f_0$ in $\cC^{1,\ep}$ associated to $\cD_0$ as in \eqref{eq:f0},
then $f_t:=f_0\circ\psi_t$ is associated to the transported domain
$\cD_t=\psi_t(\cD_0)$, and we have
\begin{equation}\label{eq:Df}
D_t \nabla f=-\nabla u\cdot\nabla f\quad\hbox{with }\ (\nabla u)_{ij}=\d_iu^j.
\end{equation}
Therefore, as $X_\lambda$ satisfies \eqref{eq:X}, we have
$$D_t(X_\lambda\cdot\nabla f)=(D_tX_\lambda)\cdot\nabla f+X_\lambda\cdot(D_t\nabla f)=0,$$
which ensures that $X_\lambda$ remains tangent to the patch for all time.
\end{proof}


\noindent{\it Example.}
As a consequence of Bony decomposition and of $\div X_{\lambda,0}\equiv0$, we have
\begin{equation*}
\d_{X_{\lambda,0}} u_0 = \hcT_{X_{\lambda,0}} u_{0} +\hT_{ \d_k u_{0}}X^k_{\lambda,0} +\div \hR(u_{0},X_{\lambda,0}).
\end{equation*}
Hence, if $u_0 \in \hB^{\Psfrac Np-1}_{p,1} \cap \hB^{\Psfrac Np+\ep-1}_{p,1}$ with $p$ satisfying \eqref{cdt:u_01},
then the conditions of Proposition \ref{prop:NN} are fulfilled.
In fact, the additional regularity $\hB^{\Psfrac Np+\ep-1}_{p,1}$ of $u_0$ implies that the flow $\psi_t$ is in $\cC^{1,\ep}$, because the solution $(u,\nabla P)$ lies in $\hE_q^\ep$ for some $q> (N-1)\slash (1-\ep)$. This is a consequence of the following result that may be obtained along the lines of the proof of
Theorem \ref{thm:main_dINS}.
\begin{prop}\label{prop:NNsmooth}
If the initial data $(\rho_0, u_0)$ are as in Theorem \ref{thm:wellposed_INS}
and if in addition $u_0$ is in $\hB^{\Psfrac{N}{q}+\ep-1}_{q,1}$ and
\begin{equation*}
\|\rho_0-1\|_{\cM\big(\hB^{\Psfrac{N}{q}+\ep-1}_{q,1}\big)} \leq c \quad
\hbox{for small constant}\,\,\, c,\,\,\, 0<\ep<1
\,\,\, \hbox{and}\,\,\, \frac{N-1}{1-\ep}< q \leq \infty,
\end{equation*}
then, beside the properties listed in Theorem \ref{thm:wellposed_INS}, the unique global solution $(\rho, u, \nabla P)$ of the system \eqref{eq:INSmu} satisfies
$$\rho\in L^\infty\big(\R_+ ; \cM \big(\hB^{\Psfrac{N}{q}+\ep-1}_{q,1}\big)\big) \quad\hbox{and}\quad (u,\nabla P)\in \hE_q^{\ep}.$$
\end{prop}

\subsection{The two-dimensional case}

Here we prove Theorem \ref{thm:N2}. As a start, we have to show that if the vorticity $\omega_0$
is given by \eqref{eq:vorticity} then $u_0$ is in $\dot B^{\Psfrac 2p-1}_{p,1}(\R^2)$.
This will be achieved by using the
fact that $u_0$ can be computed from $\omega_0$ by means
of the following \emph{Biot-Savart law}:
\begin{equation}\label{eq:BS}
u_0= K_2\ast \omega_0,
\quad \hbox{with} \quad
K_2(x):= \frac{1}{2\pi |x|^2}
\begin{pmatrix}
-x_2\\ x_1
\end{pmatrix} \cdotp
\end{equation}
Recall that $\omega_0$ is in $L^p(\R^2)$ and is supported in some ball $B(0,R)$.
Now, on the one hand, one may write for all $x$ in $B(0,2R)$,
$$
|u_0(x)|\leq\frac1{2\pi}\int_{B(0,3R)} |\omega_0(x-y)|\,\frac{dy}{|y|},
$$
which, by convolution inequalities and our choice of $p$, implies that $|u_0| \mathds{1}_{B(0,2R)}$ is in $L^r(\R^2)$ for any $r$ satisfying
\begin{equation}\label{cdt:r}
r \in \big[1,\spfrac{2p}{2-p}\big[ \subset \big[1, \spfrac{2}{1-\ep}\big[.
\end{equation}

On the other hand, for $|x|\geq2R$, owing to the zero average condition for $\omega_0$, we have
$$
u_0(x)= \frac1{2\pi} \int_{|y|\leq R} \big(K_2(x-y) - K_2(x)\big) \omega_0(y)\, dy.
$$
Therefore, by computing $K_2(x-y)-K_2(x)$, it is not difficult to see that we have for some constant $C_R$ depending only on $R$,
$$
|u_0(x)|\leq\frac{C_R}{|x|^2}\,\|\omega_0\|_{L^1}\quad\hbox{for all }\ x\ \hbox{ such that }\ |x|\geq2R.
$$
Then putting the two information together, we get $u_0$ in $L^r$ for all $r$ given by \eqref{cdt:r}.

Next, let us write that
$$
u_0=\dot S_0 u_0+(\mathrm{Id}-\dot S_0)u_0.
$$
To handle the first term, we infer from
the embedding of $L^r$ in $\dot B^{\sfrac 2p-\sfrac2r}_{p,\infty}$ for all $1<r<p<2$,
$$
\|\dot S_0 u_0\|_{\dot B^{\Psfrac 2p-1}_{p,1}}
\lesssim \|\dot S_0u_0\|_{\dot B^{\sfrac 2p-\sfrac2r}_{p,\infty}}
\lesssim \|\dot S_0 u_0\|_{L^r}
\lesssim \|u_0\|_{L^r}
\lesssim \|\omega_0\|_{L^p}.
$$
As regards the high frequency part of $u_0$, the Fourier multiplier $(\mathrm{Id}-\dot S_0)\nabla^\bot
(-\Delta)^{-1}$ is homogeneous of degree $-1$ away from a neighborhood of $0$, which yields
\begin{equation}\label{es:highw0}
\begin{split}
\|(\mathrm{Id}-\dot S_0)u_0\|_{\dot B^{\Psfrac 2p-1}_{p,1}}
=\|(\mathrm{Id}-\dot S_0)\nabla^\perp(-\Delta)^{-1}\omega_0\|_{\dot B^{\Psfrac 2p-1}_{p,1}}\\
\lesssim\|(\mathrm{Id}-\dot S_0)\omega_0\|_{\dot B^{\Psfrac2p-2}_{p,1}}
\lesssim \|\omega_0\|_{L^p}.
\end{split}
\end{equation}


Next, consider the divergence free vector field
$X_0=\nabla^\perp f_0$, where $f_0$ is given by~\eqref{eq:f0} and is (with no loss of generality)
compactly supported. If it is true that
\begin{equation}\label{eq:dX0}
\d_{X_0} u_0\in \dot B^{\Psfrac 2p-2+\ep}_{p,1},
\end{equation}
then one can apply Proposition \ref{prop:NN} which ensures that the transported vector field~$X_t$
remains in $\cC^{0,\ep}$ for all $t\geq0$. Now, it is classical that we have $X_t=(\nabla f_t)^\perp$
with $f_t=f_0\circ\psi_t^{-1}$. Hence $\cD_t$ has a $\cC^{1,\ep}$ boundary.


Let us establish \eqref{eq:dX0}.
First note that
\begin{equation}\label{spc:X}
X_0 \in \sC^{\ep}_c \hto b^{\ep}_{p,\infty} \cap b^{\alpha}_{p,r} \quad \hbox{provided }\,\,\,
1\leq r \leq \infty \,\,\hbox{and}\,\,
-\sfrac{2}{p'} <\alpha <\ep
\end{equation}
due to Proposition \ref{prop:besov-holder} and Proposition \ref{prop:equiv_besov}. Now, \eqref{es:dXf_TXf} ensures that
\begin{equation}\label{eq:BS0}
\|\hcT_{X_0}u_0-\d_{X_0}u_0\|_{\dot B^{\Psfrac 2p+\ep-2}_{p,1}}\lesssim \|u_0\|_{\dot B^{\Psfrac 2p-1}_{p,1}}
\|X_0\|_{\hsC^\ep} \quad \hbox{for any}\,\,\, 1<p<\spfrac{2}{1-\ep} \cdotp
\end{equation}
Then thanks to \eqref{eq:BS}, we obtain
$$
\hcT_{X_0}u_0=\hcT_{X_0}(-\Delta)^{-1}\nabla^\perp\omega_0
=(-\Delta)^{-1}\nabla^\perp\hcT_{X_0}\omega_0+[\hcT_{X_0},(-\Delta)^{-1}\nabla^\perp]\omega_0,
$$
whence using Lemma \ref{lemma:ce1} and \eqref{es:highw0},
\begin{equation}\label{eq:BS1}
\| \hcT_{X_0}u_0
-(-\Delta)^{-1}\nabla^\perp\hcT_{X_0}\omega_0\|_{\dot B^{\Psfrac{2}{p} +\ep-2}_{p,1}}
\lesssim \|X_0\|_{\hsC^{\ep}}\|\omega_0\|_{L^p}.
\end{equation}
Next, we notice that
$$\hcT_{X_0}\omega_0-\div(X_0\omega_0)=-\div\bigl(\hT_{\omega_0}X_0+\hR(\omega_0,X_0)\bigr).
$$
Therefore, taking advantage of standard continuity results for $\hT$ and $\hR$,
we have
\begin{equation}\label{eq:BS2}
\|\hcT_{X_0}\omega_0-\div(X_0\omega_0)\|_{\dot B^{\Psfrac{2}{p} +\ep-3}_{p,1}}\lesssim \|\omega_0\|_{L^p}
\|X_0\|_{\hsC^{\ep}},\quad \hbox{ as }\, 1<p<\spfrac{2}{2-\ep} \cdotp
\end{equation}
Since $\div(X_0 \,\overline{\omega}_0) \equiv 0$ by assumption, it is sufficient to study $\div(X_0\,\wt\omega_0\,\mathds{1}_{\mathcal{D}_0})$.
Recall that $\wt\omega_0 \in \sC^{\ep'}_c$ for some $0<\ep'<\ep$,
and that $\div X_0=0$. As we have $\d_{X_0}\mathds{1}_{\mathcal{D}_0}=0$,
Corollary \ref{cor:dXfg_3D} implies that
\begin{equation*}
\div(X_0\,\wt\omega_0\,\mathds{1}_{\mathcal{D}_0})
\in \dot B^{\Psfrac{2}{p}+\ep-3}_{p,1}.
\end{equation*}
Putting \eqref{eq:BS0}, \eqref{eq:BS1} and \eqref{eq:BS2} together, we conclude that \eqref{eq:dX0} is fulfilled, which completes the proof of Theorem \ref{thm:N2}. \qed
\medskip

If dropping off the zero average condition for the function $\omega_0$ in Theorem \ref{thm:N2}, then the corresponding initial velocity field $u_0$
cannot be in $L^r(\R^2)$ for any $r\in{}]1,2]$. Still, one can get a similar statement in the particular case where
$(\overline{\omega}_0, \wt \omega_0) \equiv (0,\eta')$
for some small enough $\eta'$.
Indeed, from \eqref{eq:BS} and Hardy-Littlewood-Sobolev inequality, we deduce that $u_0$ belongs to all spaces $L^r(\R^2)$ with $r\in{}]2,\infty[$.
Repeating the first part of the proof of Theorem \ref{thm:N2} thus yields
$u_0 \in \hB^{\Psfrac{2}{p}-1}_{p,1}(\R^2)$ for any $2<p<\infty$. Now, as $\omega_0$ is bounded and
compactly supported, it is in $\hB^{\Psfrac{2}{q} + \ep-2}_{q,1}(\R^2)$ for any $0<\ep<1$ and $1<q<\infty$,
which implies that $u_0 \in \hB^{\Psfrac{2}{q} + \ep-1}_{q,1}$.
Hence, applying Proposition \ref{prop:NNsmooth}, and using the fact that the flow of the solution
constructed therein is in $\cC^{1,\ep}$, we conclude to the following generalization of of \cite[Rem.\,1.1]{LiaoZh2016a}.
\begin{thm}\label{thm:N22}
Let $\cD_0$ satisfy \eqref{eq:f0} for some $\ep$ in $]0,1[$.
There exists a constant $\eta_0$ depending only on $\cD_0$ so that for all
$\eta,\eta' \in{}]-\eta_0,\eta_0[$
if
\begin{equation*}
\rho_0:= (1+\eta) \mathds{1}_{\mathcal{D}_0} + \mathds{1}_{\mathcal{D}^c_0},
\end{equation*}
and if the divergence free vector-field $u_0$ in $W^{1,p}(\R^2)$ for some $p>2$ is given by
\begin{equation*}
u_0:=(-\Delta)^{-1}\nabla^\perp (\eta' \mathds{1}_{\mathcal{D}_0}),
\end{equation*}
then the system \eqref{eq:INSmu} has a unique solution $(\rho, u,\nabla P)$ with the properties listed in Proposition \ref{prop:NNsmooth} for some suitable
In addition, \eqref{eq:rhot} is fulfilled for all $t\geq0$, and $\mathcal{D}_t$ remains a simply connected bounded domain of class $\cC^{1,\ep}$.
\end{thm}

\subsection{The three-dimensional case}

As another application of Proposition \ref{prop:NN}, one can generalize Theorem \ref{thm:N2}
to the three-dimensional case. Our result reads as follows.
\begin{thm}\label{thm:N3}
Let $\cD_0$ be a $\cC^{1,\ep}$ simply connected bounded domain of $\R^3$ with
$\ep\in{}]0,1[$. Let $\rho_0$ be given by \eqref{eq:density_patch}
with small enough $\eta$. Assume that the initial velocity $u_0$ has
coefficients in $\cS'_h(\R^3)$ and
vorticity\footnote{For any point $Y \in \R^3$, we set $X \wedge Y:= (X^2 Y^3-X^3Y^2,\, X^3 Y^1-X^1Y^3,\, X^1 Y^2-X^2Y^1)$, where~$X$ stands for an element of $\R^3$ or for the $\nabla$ operator.}
\begin{equation*}
\Omega_0:=\nabla\wedge u_0=\wt{\Omega}_0 \mathds{1}_{\mathcal{D}_0},
\end{equation*}
for some small enough $\wt{\Omega}_0$ in $ \cC^{0,\delta}(\R^3;\R^3)$ ($\delta \in{}]0,\ep[$) with $\div \wt{\Omega}_0 =0$ and \hbox{$\wt{\Omega}_0 \cdot \vec{n}_{_{\cD_0}} |_{\d \cD_0} \equiv 0$} (here $\vec{n}_{_{\cD_0}}$ denotes the outwards unit normal of the domain $\cD_0$).

There exists a unique solution $(\rho, u,\nabla P)$ to the system \eqref{eq:INSmu} with the properties
listed in Theorem \ref{thm:wellposed_INS} for some suitable $p$ satisfying
\begin{equation}\label{eq:p3}
2<p<\min\Big\{\spfrac2{1-\ep}, 6\Big\}\cdotp
\end{equation}
Furthermore, for all $t\geq0$, we have \eqref{eq:density_patcht}
and $\mathcal{D}_t$ remains a simply connected bounded domain of class $\cC^{1,\ep}$.
\end{thm}
\begin{proof} With no loss of generality, one may assume that $\wt\Omega_0$ is compactly supported.
Like in the 2D case, we first have to check that $u_0$ fulfills the assumptions of Proposition \ref{prop:NN}.
As it is divergence free and decays at infinity (recall that $u_0\in\cS'_h$), it is
given by the Biot-Savart~law:
\begin{equation}\label{eq:BS3}
u_0=(-\Delta)^{-1}\nabla\wedge\Omega_0, \quad \hbox{with}\,\,\,\Omega_0 = \wt{\Omega}_0 \, \mathds{1}_{\mathcal{D}_0}.
\end{equation}
We claim that $u_0$ belongs to
$\hB^{\Psfrac{3}{p}-1}_{p,1}$ for some $p$ satisfying \eqref{eq:p3}. Indeed,
the characteristic function of any bounded domain with $\cC^1$ regularity
belongs to all Besov spaces $B^{\sfrac1q}_{q,\infty}$ with $1\leq q\leq\infty$ (see \eg \cite{Tri1992}).
Hence combining Proposition \ref{prop:equiv_besov} and the embedding \eqref{es:nonhom_hom} gives
\begin{equation}\label{embd_D0}
\mathds{1}_{\mathcal{D}_0} \in \ \cE' \cap B^{\sfrac{1}{q}}_{q,\infty} \hto \hB^{\Psfrac{3}{q}-2}_{q,1} \quad \hbox{for any} \,\,\, q \in{}]1,\infty[.
\end{equation}
Now, using Bony's decomposition and standard continuity results
for operators $\hR$ and~$\hT$, we discover that
\begin{equation*}
\wt{\Omega}_0 \in \sC^{\delta}_c \hto \cM \big(\hB^{\Psfrac{3}{q}-2}_{q,1}\big) \quad \hbox{for any} \,\,\, q \in\big ]\sfrac{3}{2},\, \spfrac{3}{2-\delta}\big[.
\end{equation*}

Hence the definition of Multiplier space and \eqref{embd_D0} yield
\begin{equation}\label{cdt:q0}
\Omega_0 = \wt{\Omega}_0\, \mathds{1}_{\mathcal{D}_0}\in \hB^{\Psfrac{3}{q}-2}_{q,1} \quad \hbox{for any} \,\,\, q \in \big]\sfrac{3}{2},\, \spfrac{3}{2-\delta}\big[.
\end{equation}
As $u_0$ is in $\cS'_h$ and $(-\Delta^{-1})^{-1}\nabla\wedge$ in \eqref{eq:BS3} is a homogeneous multiplier of degree $-1$, one can conclude that $$u_0 \in \hB^{\Psfrac{3}{q}-1}_{q,1} \hto \hB^{\Psfrac{3}{p}-1}_{p,1}, \quad \mbox{for any $p\geq q$.}$$
Note that for any $\delta$ in $]0,1[$, one can find some $p$ satisfying the above conditions and \eqref{eq:p3} altogether.


Next, consider some (compactly supported) level set function $f_0$ associated to $\d\cD_0$, and the three $\cC^{0,\ep}$ vector-fields $X_{k,0}:=e_k\wedge \nabla f_0$ with $(e_1,e_2,e_3)$ being the canonical basis of $\R^3$. It is clear that those vector-fields are divergence free and tangent to $\d\cD_0$. Let us check that we have $\d_{X_{k,0}}u_0\in\dot B^{\Psfrac 3p-2+\ep}_{p,1}$ for some
$p$ satisfying \eqref{eq:p3}.
As in the two-dimensional case, this will follow from Biot-Savart law and the special structure of $\Omega_0$. Indeed, from \eqref{es:dXf_TXf} and $\div X_{k,0}=0$, we have
$$
\|\hcT_{X_{k,0}}u_0-\d_{X_{k,0}}u_0\|_{\dot B^{\Psfrac 3p+\ep-2}_{p,1}}
\lesssim \|u_0\|_{\dot B^{\Psfrac3p-1}_{p,1}}\|X_0\|_{\hsC^\ep}, \,\,\, \forall \,p \in \big]\sfrac{3}{2}, \spfrac{3}{1-\ep}\big[.
$$
Then \eqref{eq:BS3} yields
$$
\hcT_{X_{k,0}}u_0=\hcT_{X_{k,0}}(-\Delta)^{-1}\nabla\wedge\Omega_0
=(-\Delta)^{-1}\nabla\wedge\hcT_{X_{k,0}}\Omega_0+[\hcT_{X_{k,0}},(-\Delta)^{-1}\nabla\wedge]\,\Omega_0.
$$
Thanks to Lemma \ref{lemma:ce1} and homogeneity of $(-\Delta^{-1})^{-1}\nabla\wedge$, it is thus sufficient to verify
that $\hcT_{X_{k,0}}\Omega_0$ belongs to $\hB^{\Psfrac{3}{p}+\ep-3}_{p,1}$ for some $p$ satisfying \eqref{eq:p3}. In fact, from the decomposition
$$\hcT_{X_{k,0}}\Omega_0-\div(X_{k,0}\Omega_0)
=-\div\bigl(\hT_{\Omega_0}X_{k,0}+\hR(\Omega_0,X_{k,0})\bigr),
$$
and continuity results for $\hR$ and $\hT$, we get
\begin{equation*}
\|\hcT_{X_{k,0}}\Omega_0-\div(X_{k,0}\Omega_0)\|_{\dot B^{\Psfrac{3}{q}+\ep-3}_{q,1}}\lesssim \|\Omega_0\|_{\hB^{\Psfrac{3}{q}-2}_{q,1}} \|X_{k,0}\|_{\hsC^\ep}, \,\,\, \forall\, q \in \big]\sfrac{3}{2}, \spfrac{3}{2-\ep}\big[.
\end{equation*}
Thus, remembering \eqref{cdt:q0} and $0<\delta < \ep$, we have to choose some $p$ satisfying \eqref{eq:p3},
such that the following standard embedding holds
\begin{equation}\label{cdt:q1}
\dot B^{\Psfrac{3}{q}+\ep-3}_{q,1} \hto \dot B^{\Psfrac{3}{p}+\ep-3}_{p,1} \,\,
\hbox{for some}\,\,\, q \in \big]\sfrac{3}{2}, \spfrac{3}{2-\delta}\big[ \,\,\,\hbox{with}\,\,\, q \leq p.
\end{equation}

Now, because $\d_{X_{k,0}}\mathds{1}_{\mathcal{D}_0}\equiv0$ and
$\wt\Omega_0$ is in $\cC^{\delta}$,
Corollary \ref{cor:dXfg_3D} yields for all $0 < \delta_{\star}<\delta$,
$$\d_{X_{k,0}} \Omega_0=\div(X_{k,0}\otimes\Omega_0)=\div(X_{k,0}\otimes
\wt\Omega_0\,\mathds{1}_{\mathcal{D}_0})
\in \dot B^{\delta_{\star}-1}_{q,1}\ \hbox{ for all }\ q\geq1.$$
One can thus conclude that $\d_{X_{k,0}}u_0\in\dot B^{\Psfrac3p-2+\ep}_{p,1}$ for any index $p$ satisfying
$p\geq q$ with $q$ satisfying the condition \eqref{cdt:q1} and $\Psfrac3q+\ep-2=\delta^*\in{}]0,\delta[$.

As one can require in addition $p$ to fulfill \eqref{eq:p3},
Proposition \ref{prop:NN} applies with the family $(X_{k,0})_{1\leq k\leq3}$.
Denoting by $(X_k)_{1\leq k\leq3}$
the corresponding family of divergence free vector fields in $\cC^{0,\ep}$ given by \eqref{eq:X} with initial data $X_{0,k}$,
and introducing $Y_1:=X_3\wedge X_1$, $Y_2:=X_3\wedge X_1$
and $Y_3=X_1\wedge X_2$, we discover that for $\alpha=1,2,3$,
\begin{equation}\label{eq:Y}
\begin{cases}
\d_t Y_\alpha +u\cdot \nabla Y_\alpha = - \nabla u\cdot Y_\alpha, \\[3pt]
(Y_\alpha)|_{t=0}=\d_\alpha f_0\, \nabla f_0.
\end{cases}
\end{equation}
From \eqref{eq:Df}, it is clear that the time-dependent
vector field $\bigl(\d_\alpha f_0(\psi^{-1}_t) \big)\,\nabla f_t$ also satisfies \eqref{eq:Y},
hence we have, by uniqueness, $Y_\alpha (t,\cdot) = \big((\d_\alpha f_0) (\psi^{-1}_t) \big)\nabla f_t$.
So finally,
$$
\big|\nabla f_0\circ \psi_t^{-1}\bigr|^2\,\nabla f_t=\sum_{\alpha=1}^3 Y_\alpha(t,\cdot)\,\d_\alpha f_0\circ\psi_t^{-1}.
$$
As $\psi_t^{-1}$ is $\cC^1$ and as both $Y_\alpha$ and $\nabla f_0$
are in $\cC^{0,\ep}$, one can
conclude that $ \nabla f_t$ is $\cC^{0,\ep}$ in some neighborhood of $\d\cD_0$.
Therefore $\cD_t$ remains of class $\cC^{1,\ep}$ for all time.
\end{proof}

\begin{rmk}
In contrast with the 2D case, one cannot consider constant vortex patterns for
the condition $\wt{\Omega}_0 \cdot \vec{n}_{_{\cD_0}} |_{\d \cD_0} \equiv 0$ is not fulfilled.
One can define directly~$u_0$ through
$u_0=(-\Delta)^{-1}\nabla\wedge e$, where $e$ is a constant vector of $\R^3$
(as we did for the Boussinesq system in \cite{DanZhX2017}), but then, $\nabla\wedge u_0$ does not coincides with $e$.
\end{rmk}

\section{The proof of persistence of striated regularity}\label{sect:3}

That section is devoted to the proof of Theorem \ref{thm:main_dINS}.
The first step is to apply Theorem \ref{thm:wellposed_INS}. From it, we get a unique global solution $(\rho,u,\nabla P)$ with
$\rho\in \cC_b\big(\R_+;\cM(\dot B^{\Psfrac Np-1}_{p,1})\big)$ and $(u,\nabla P)\in \dot E_p$,
satisfying \eqref{est:up}.
Because the product of functions maps $\dot B^{\Psfrac Np-1}_{p,1}\times\dot B^{\sfrac Np}_{p,1}$
to $\dot B^{\Psfrac Np-1}_{p,1}$, we deduce that
$D_tu=\d_tu+u\cdot\nabla u$ is also bounded by the right-hand side of \eqref{est:up}.
So finally,
\begin{equation}\label{es:u_INS}
\|(u, \nabla P)\|_{\hE_p} + \|D_t u\|_{L^1_t(\hB^{\Psfrac{N}{p}-1}_{p,1})}
\leq C\|u_0\|_{\hB^{\Psfrac{N}{p}-1}_{p,1} }.
\end{equation}

In order to complete the proof of the theorem, it
is only a matter of showing that the additional multiplier and striated regularity properties are conserved
for all positive times. We shall mainly concentrate on the proof of a priori estimates
for the corresponding norms, just
explaining at the end how a suitable regularization
process allows to make it rigorous.

\subsection{Bounds involving multiplier norms}

As already pointed out in the introduction, because $\nabla u$ is in $L^1(\R_+;\dot B^{\sfrac Np}_{p,1})$
and $\dot B^{\sfrac Np}_{p,1}$ is embedded in $\cC_b$, the flow $\psi$ of $u$ is~$\cC^1$ and we have
for all $t\geq0$, owing to \eqref{est:up},
\begin{equation}\label{es:dflow}
\|\nabla \psi_t^{\pm1}\|_{L^\infty}\leq \exp\biggl(\int_0^t\|\nabla u\|_{L^\infty}\,d\tau\biggr)
\leq C
\end{equation}
for a suitably large universal constant $C$.

Now, from the mass conservation equation and \eqref{eq:dXrho}, we gather that
$$
\rho(t,\cdot)=\rho_0\circ\psi_t^{-1}\quad\hbox{and}\quad
(\d_X\rho)(t,\cdot)=(\d_{X_0}\rho_0)\circ\psi_t^{-1}.
$$
Hence $\|\rho(t,\cdot)\|_{L^\infty}$ is time independent, and
Lemma \ref{lemma:multiplier} (keeping in mind the condition \eqref{cdt:p}) guarantees
that for all $t\in\R_+$,
\begin{align}\label{es:1-rho_1}
\|\rho(t) -1 \|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1})} &\leq C \|\rho_0 -1 \|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1})},\\\label{es:1-rho_2}
\|\rho(t) -1 \|_{\cM (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} &\leq C \|\rho_0 -1 \|_{\cM (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})},\\\label{es:dXrho}
\quad\|(\d_X \rho)(t)\|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1} \to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} &\leq C \|\d_{X_0} \rho_0\|_{ \cM(\hB^{\Psfrac{N}{p}-1}_{p,1} \to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}.
\end{align}

\subsection{ Estimates for the striated regularity}
Recall that $\hcT_X u$ satisfies the Stokes-like system \eqref{eq:dXINS}.
As $\hcT_X u$ need not be divergence free, to enter into the standard theory, we set
$$ v:= \hcT_X u - w\quad\hbox{with}\quad w:= \hT_{\d_k X} u^k-\hT_{\div X}u.$$
Denoting
$\widetilde{g}:= g- \rho u\cdot \nabla \hcT_X u -(\rho\d_t w - \Delta w)$ with $g$ defined in \eqref{eq:g_dINS}, we see that $v$ satisfies:
\begin{equation} \label{eq:v_dINS} \tag{S}
\begin{cases}
\rho \d_t v - \Delta v + \nabla \hcT_X P = \widetilde{g},\\
\div v = 0,\\
v|_{t=0}=v_0.
\end{cases}
\end{equation}

We shall decompose the proof of a priori estimates for striated regularity into three steps.
The first one is dedicated to bounding $\widetilde{g}$ (which mainly requires the commutator estimates
of the appendix). In the second step, we take advantage of the smoothing effect of the heat
flow so as to estimate $v$.
In the third step, we revert to $\hcT_X u$ and eventually bound $X$.

\subsubsection*{First step: bounds of ${\wt g}$} Recall that $\widetilde{g}:= g- \rho u\cdot \nabla \hcT_X u -(\rho\d_t w - \Delta w)$ with
\begin{multline*}
g= - \rho [\hcT_{X}, D_t]u +[\hcT_X,\Delta] u - [\hcT_X, \nabla] P + (\d_X-\hcT_X) (\Delta u - \nabla P)\\
-\d_X \rho D_t u + \rho (\hcT_X-\d_X)D_t u.
\end{multline*}

The first term of $g$ may be bounded according to Proposition \ref{prop:ce_dINS} and to the definition of multiplier spaces. We get, under assumption \eqref{cdt:p_X},
\begin{multline}\label{es:g1_dINS}
\|\rho [\hcT_X, D_t ]u\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \Big(\|u\|_{\hsC^{-1}}\|\hcT_X u\|_{ \hB^{\Psfrac{N}{p}+\ep}_{p,1}} \\
+ \|u\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}} \|\hcT_{X} u\|_{\hsC^{\ep-2}}
+ \|u\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}} \|u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \|X\|_{\hsC^{\ep}} \Big).
\end{multline}

Next, thanks to the commutator estimates in Lemma \ref{lemma:ce1}, we have
\begin{align}\label{es:g2_dINS}
\|[\hcT_X,\Delta] u\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} &\lesssim \|\nabla X\|_{\hsC^{\ep-1}}\|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}},
\\
\label{es:g3_dINS}
\|[\hcT_X, \nabla ]P\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} &\lesssim \|\nabla X\|_{\hsC^{\ep-1}} \|\nabla P\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}}.
\end{align}
Bounding the fourth term of $g$ stems from \eqref{es:dXf_TXf}: we have
\begin{equation}\label{es:g4_dINS}
\| (\hcT_X -\d_X) (\Delta u -\nabla P) \|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|(\Delta u ,\nabla P) \|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \|X\|_{\hsC^{\ep}}.
\end{equation}
Then the definition of multiplier spaces yields
\begin{equation}\label{es:g5_dINS}
\|\d_X \rho D_t u\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|\d_X \rho\|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1} \to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \|D_t u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}}.
\end{equation}
Finally, using again \eqref{es:dXf_TXf}
and the definition of multiplier spaces, we may write
\begin{equation}\label{es:g6_dINS}
\| \rho (\hcT_X - \d_X) D_t u\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \| \rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \|X\|_{\hsC^\ep} \|D_t u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}}.
\end{equation}
Putting together \eqref{es:g1_dINS} -- \eqref{es:g6_dINS} and integrating with respect to time,
we end up with
\begin{multline}\label{es:g_dINS}
\|g\|_{L^1_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\\[-3pt]
\lesssim \!
\int^t_0\! \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \big(\|u\|_{\hsC^{-1}}\|\hcT_X u\|_{ \hB^{\Psfrac{N}{p}+\ep}_{p,1}} \!+\! \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \|\hcT_{X}u\|_{\hsC^{\ep-2}} \big)\,dt' \\[-3pt]
\shoveleft{+ \! \int_0^t \|X\|_{\hsC^{\ep}} \Big(\! \big(\|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \|u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}}
+ \|D_t u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \big)\|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}}\\[-3pt]
\shoveright{+ \| (\nabla^2 u, \nabla P)\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \Big)dt'} \\[-3pt]
+ \int_0^t\|\d_X \rho\|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1} \to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \| D_t u \|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \,dt'.
\end{multline}

Bounding the second term of $\widetilde{g}$ is obvious: taking advantage of Bony's decomposition
\eqref{eq:bony}
and remembering that $\Psfrac{N}{p}+\ep > 1 $ and that $\div u =0$, we get
\begin{multline}\label{es:tg1_dINS}
\|\rho u\cdot \nabla \hcT_X u\|_{L^1_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\\[-3pt]
\lesssim \int^t_0 \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \big(\|u\|_{\hsC^{-1}}\|\hcT_{X}u\|_{\hB^{\Psfrac{N}{p}+\ep}_{p,1}}
+\|u\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}\|\hcT_{X}u\|_{\hsC^{\ep-2}} \big) \,dt'.
\end{multline}
To bound the last term of $\wt g$, we use the decomposition
$$
\rho \d_t w-\Delta w= \rho (W_1+W_2)+ W_3,\quad\text{with}\quad
\begin{cases}
W_1:= \hT_{\d_k X}\d_tu^k-\hT_{\div X}\d_tu ,\\
W_2:= \hT_{\d_k \d_tX}u^k-\hT_{\div\d_tX}u , \\
W_3:=\Delta\bigl(\hT_{\div X}u-\hT_{\d_k X}u^k\bigr).
\end{cases}
$$
Continuity results for the paraproduct and the definition of
$\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})$ ensure that
\begin{align}\label{eq:W1}
\|\rho W_1 \|_{L^1_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} & \lesssim \int^t_0 \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \|\nabla X \|_{\hsC^{\ep-1}} \|\d_t u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}}\,dt',\\[-3pt]\label{eq:W2}
\| \rho W_2 \|_{L^1_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} &\lesssim \int^t_0 \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \|\d_t X\|_{\hsC^{\ep-2}}\|u\|_{\hB^{{\Psfrac Np}+1}_{p,1}}\,dt',\\[-3pt]\label{eq:W3}
\| W_3\|_{L^1_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} &\lesssim\!\int^t_0 \| \nabla X \|_{\hsC^{\ep-1}}\|u\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}\,dt'.
\end{align}
To estimate $\d_tX$ in \eqref{eq:W2}, we use the fact that
$$
\d_tX=-u\cdot\nabla X+\d_Xu=-\div(u\otimes X)+\d_Xu.
$$
Hence using \eqref{eq:bony}, and
continuity results for the remainder and paraproduct operators, we get under the condition \eqref{cdt:p_X},
$$\|\d_tX \|_{\hsC^{\ep-2}}\lesssim
\|u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \|X\|_{\hsC^{\ep}}+ \|\d_Xu\|_{\hsC^{\ep-2}}.$$
Therefore, taking advantage of \eqref{es:dXf_TXf} yields
\begin{multline}\label{eq:W4}
\| \rho W_2 \|_{L^1_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\\ \lesssim \int^t_0 \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} (\|X\|_{\hsC^{\ep}}\|u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} + \|\hcT_X u\|_{\hsC^{\ep-2}}) \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}}\, dt'.
\end{multline}
Combining \eqref{eq:W1}, \eqref{eq:W2} and \eqref{eq:W4}, we eventually obtain
\begin{multline}\label{es:tg2_dINS}
\|\rho \d_t w - \Delta w \|_{L^1_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\\[-3pt]
\lesssim \int^t_0 \|\hcT_X u\|_{\hsC^{\ep-2}} \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \,dt'\\[-3pt]
\shoveright{+\int^t_0 \|X\|_{\hsC^{\ep}} \Big(\big(\|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \|u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} +1 \big) \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}}}\\[-3pt]
+ \|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \|\d_t u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \Big) \, dt',
\end{multline}
whence, putting together estimate \eqref{es:g_dINS}, \eqref{es:tg1_dINS} and \eqref{es:tg2_dINS},
\begin{multline}\label{es:tg_dINS}
\|\widetilde{g}\|_{L^1_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\\[-3pt]
\lesssim
\!\int^t_0 \!\|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \big(\|u\|_{\hsC^{-1}}\|\hcT_X u\|_{ \hB^{\Psfrac{N}{p}+\ep}_{p,1}}\! +\! \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \|\hcT_{X}u\|_{\hsC^{\ep-2}} \big)\,dt' \\[-3pt]
\shoveright{+ \int_0^t \|X\|_{\hsC^{\ep}} \big(\|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \|u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} + \|(\d_t u,D_t u)\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \big)\|\rho\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \,dt'}
\\[-3pt]
+ \int_0^t \|X\|_{\hsC^{\ep}} \|(\nabla^2 u, \nabla P)\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \,dt'\\[-3pt]
+ \int_0^t\|\d_X \rho\|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1} \to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \| D_t u \|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \,dt'.
\end{multline}

\subsubsection*{Second step: bounds of ${v}$}

We now want to bound $v$ in
\[
\wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}) \cap L^{1}_t (\hB^{\Psfrac{N}{p}+\ep}_{p,1}),
\]
knowing \eqref{es:tg_dINS}.
This will follow from the smoothing properties of the heat flow.
More precisely, introduce the projector $\mathbb P$ over divergence-free vector fields,
and apply $\mathbb{P}\hDelta_j$ (with $j\in\Z$) to the equation \eqref{eq:v_dINS}. We get
\begin{equation*}
\begin{cases}
\d_t \hDelta_jv - \Delta \hDelta_jv = \mathbb{P}\hDelta_j(\widetilde{g} + (1-\rho)\d_t v)\\[1ex]
\hDelta_jv|_{t=0} = \hDelta_jv_{0}.
\end{cases}
\end{equation*}
Lemma 2.1 in \cite{Che1999} implies that if $p\in [1, \infty]$,
\begin{align*}
\|\hDelta_jv(t)\|_{L^{p}} \leq e^{-ct 2^{2j}} \|\hDelta_jv_{0}\|_{L^{p}} + C \int^{t}_{0} e^{-c(t-t') 2^{2j}}
\|\hDelta_j(\widetilde{g} + (1-\rho)\d_t v)(t')\|_{L^{p}}\,dt'.
\end{align*}
Therefore, taking the supremum over $j\in \Z$, using the fact that
\begin{equation*}
\d_t v = \Delta v + \mathbb{P} \big(\widetilde{g} + (1-\rho)\d_t v\big)
\end{equation*}
and that $\mathbb{P}:\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}\to\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}$,
we find that
\begin{multline}\label{es:v1_dINS}
\|v\|_{\wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}
+ \|v\|_{L^{1}_t(\hB^{\Psfrac{N}{p}+\ep}_{p,1})} +
\|\d_tv\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \\
\lesssim \|v_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} } + \|\wt{g}\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}+ \|(1-\rho)\d_t v\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}.
\end{multline}
The smallness condition \eqref{cdt:u0_small} combined with Inequality \eqref{es:1-rho_2} ensure
that the last term of \eqref{es:v1_dINS} may be absorbed by the left-hand side, and we thus
end up with
\begin{multline*}
\|v\|_{\wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}) \cap L^{1}_t (\hB^{\Psfrac{N}{p}+\ep}_{p,1})} + \|\d_t v\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\\
\lesssim \|v_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} } + \|\wt{g}\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}.
\end{multline*}
Next, we use the fact that by definition of $v_0$,
$$\begin{aligned}
v_0 &= \hcT_{X_0} u_0 - \hT_{\d_k X_0} u_0^k + \hT_{\div X_0}u_0 \\
&= \d_{X_0} u_0- \hT_{\d_k u_0}X^k_0 - \d_k \hR(X^k_0, u_0) + \hR(\div X_0, u_0) - \hT_{\d_k X_0} u_0^k + \hT_{\div X_0}u_0.
\end{aligned}$$
Hence continuity results for the paraproduct yield, under the condition \eqref{cdt:p_X},
\begin{equation*}
\|v_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} } \lesssim \|\d_{X_0} u_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} } +\| X_0\|_{\hsC^{\ep}} \|u_0\|_{\hB^{\Psfrac{N}{p} -1}_{p,1}}.
\end{equation*}
Thus
\begin{multline} \label{es:v_dINS}
\|v\|_{\wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}) \cap L^{1}_t (\hB^{\Psfrac{N}{p}+\ep}_{p,1})} + \|\d_t v\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \\ \lesssim \|\d_{X_0} u_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} } +\| X_0\|_{\hsC^{\ep}} \|u_0\|_{\hB^{\Psfrac{N}{p} -1}_{p,1}}
+ \|\wt{g}\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}.
\end{multline}

\subsubsection*{Third step: bounds for striated regularity}
Remembering that
$$\hcT_X u=v+w\quad\hbox{with}\quad w=\hT_{\d_k X}u^k - \hT_{\div X}u,$$
it is now easy to bound the following quantity:
\begin{equation*}
\sH(t):= \|\hcT_X u\|_{\wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}+\|\hcT_X u\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep}_{p,1})}+\|\nabla \hcT_X P\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}.
\end{equation*}
Indeed, we have
\begin{equation}\label{eq:dPi}
\nabla \hcT_X P= (\mathrm{Id}- \mathbb{P}) (\widetilde{g} -\rho \d_t v),
\end{equation}
and thus $\|\nabla \hcT_X P\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}$
may be bounded by the right-hand side of \eqref{es:v_dINS}. Note also that
continuity results for paraproduct operators guarantee that
$$
\begin{aligned}
\|w\|_{\wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} &\lesssim
\|u\|_{ \wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}-1}_{p,1})}\|X\|_{L^{\infty}_t (\hsC^{\ep})},\\
\|w\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep}_{p,1})} &\lesssim
\int^t_0 \|u\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}} \|\nabla X\|_{\hsC^{\ep-1}}\,dt'.
\end{aligned}
$$
Hence we have
\begin{multline}\label{es:w_dINS1}
\sH(t)\lesssim \|\d_{X_0} u_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} }
+\| X_0\|_{\hsC^{\ep}} \|u_0\|_{\hB^{\Psfrac{N}{p} -1}_{p,1}}
+ \|\wt{g}\|_{L^{1}_t (\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\\
+ \|u\|_{ \wt{L}^{\infty}_t (\hB^{\Psfrac{N}{p}-1}_{p,1})\cap L^1_t(\hB^{\Psfrac Np+1}_{p,1})}
\|X\|_{L_t^\infty(\hsC^\ep)}.
\end{multline}
Because $X$ satisfies \eqref{eq:X}, standard Hölder estimates for transport equations imply that
$$
\| X \|_{L^{\infty}_t (\hsC^{\ep})} \leq \|X_0\|_{\hsC^{\ep}} + \int^t_0 \|\nabla u\|_{L^\infty} \| X \|_{\hsC^{\ep}}\, dt' + \int_0^t\|\d_X u\|_{\hsC^{\ep}}\,dt'.
$$
Now, recall that
$$
\d_Xu-\dot\cT_Xu=\dot T_{\d_ku}X^k+\dot R(\d_ku,X^k).
$$
Hence, using standard continuity results for operators $\dot T$ and $\dot R$, and embedding,
\begin{equation}\label{eq:paraX}
\|\dot\cT_Xu-\d_X u\|_{\hsC^{\ep}}\lesssim
\|\dot\cT_Xu-\d_X u\|_{\hB^{\Psfrac{N}{p}+\ep}_{p,1}}\lesssim
\|\nabla u\|_{\dot B^{\sfrac Np}_{p,1}}\|X\|_{\hsC^\ep}.
\end{equation}
Therefore we have
\begin{equation}\label{es:X_dINS}
\| X \|_{L^{\infty}_t (\hsC^{\ep})}
\leq \|X_0\|_{\hsC^{\ep}} + \int^t_0 \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \| X \|_{\hsC^{\ep}}\, dt' + \|\hcT_X u\|_{L^1_t (\hB^{\Psfrac{N}{p}+\ep}_{p,1})}.
\end{equation}
Then, using \eqref{es:u_INS} and plugging the above inequality in \eqref{es:w_dINS1}, we get
$$\displaylines{
\sH(t)\lesssim \|\d_{X_0} u_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} }
+\| X_0\|_{\hsC^{\ep}} \|u_0\|_{\hB^{\Psfrac{N}{p} -1}_{p,1}}
+ \|\wt{g}\|_{L^{1}_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \hfill\cr\hfill
+\|u_0\|_{\dot B^{\Psfrac Np-1}_{p,1}}\biggl(\|\hcT_X u\|_{L^1_t (\hB^{\Psfrac{N}{p}+\ep}_{p,1})}
+ \int^t_0 \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \| X \|_{\hsC^{\ep}}\, dt'\biggr)\cdotp}
$$
Choosing $c$ small enough in \eqref{cdt:u0_small}, we see that the first term of
the second line may be absorbed by the left-hand side.
Therefore, setting $$\sK(t):= \sH(t)+\| X \|_{L^{\infty}_t (\hsC^{\ep})}$$
and using again \eqref{es:X_dINS} and the smallness of $u_0$,
$$
\sK(t)\lesssim \|\d_{X_0} u_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} }
+\| X_0\|_{\hsC^{\ep}} + \|\wt{g}\|_{L^{1}_t(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}
+ \int^t_0 \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \| X \|_{\hsC^{\ep}}\, dt'.
$$
In order to close the estimates, it suffices to bound $\wt g$ by means of \eqref{es:tg_dINS}.
Then the above inequality becomes, after using \eqref{es:1-rho_2} and \eqref{es:dXrho} (and the fact
that \hbox{$\|\rho_0-1\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}$} is small implies that
$\|\rho_0\|_{\cM(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}$ is of order one),
\begin{multline*}
\sK(t)\lesssim \|\d_{X_0} u_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} }
+\| X_0\|_{\hsC^{\ep}}\\[-3pt]
\shoveright{+\int^t_0 \big(\|u\|_{\hsC^{-1}}\|\hcT_X u\|_{ \hB^{\Psfrac{N}{p}+\ep}_{p,1}}\! +\! \|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \|\hcT_{X}u\|_{\hsC^{\ep-2}} \big)\,dt'}\\[-3pt]
\shoveright{+ \int_0^t \|X\|_{\hsC^{\ep}} \big(\|\nabla u\|_{\hB^{\sfrac{N}{p}}_{p,1}} \|u\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} + \|(\d_t u,D_t u,\nabla^2 u, \nabla P)\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \big) \,dt'}\\[-3pt]
+\|\d_{X_0}\rho_0\|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1} \to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1})} \int_0^t \| D_t u \|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \,dt'.
\end{multline*}
The smallness of $u_0$ and \eqref{est:up} imply that all the terms of the right-hand side (except for the ones pertaining to the data), may be
absorbed by the left-hand side. Therefore using the bounds for $D_tu$ in \eqref{es:u_INS}, we eventually get
\begin{multline}\label{eq:sK}
\sK(t)\lesssim \|\d_{X_0} u_0\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1} }
+\| X_0\|_{\hsC^{\ep}}\\
+\|\d_{X_0} \rho_0\|_{\cM(\hB^{\Psfrac{N}{p}-1}_{p,1} \to \hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\|u_0\|_{\dot B^{\Psfrac Np-1}_{p,1}}.
\end{multline}

{}From \eqref{eq:paraX}, we gather that $\d_Xu$ is bounded by the right-hand side of \eqref{eq:sK}.
Next, in order to control the whole nonhomogeneous Hölder norm of $X$, it suffices to remember that
$$
\|X\|_{\cC^{0,\ep}}=\|X\|_{L^\infty}+\|X\|_{\dot\sC^\ep}
$$
and that Relation \eqref{eq:Xt} together with \eqref{es:dflow} directly yield
$$
\|X_t\|_{L^\infty}\leq\|\d_{X_0}\psi_t\|_{L^\infty}\leq C\|X_0\|_{L^\infty}.
$$
Finally, to estimate $\d_X\nabla P$, we use Inequality \eqref{es:dXf_TXf} and get
$$
\|\d_X\nabla P-\nabla\hcT_XP\|_{L_t^1(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}\lesssim \|X\|_{L_t^\infty(\hsC^\ep)}\|\nabla P\|_{L^1_t(\hB^{\Psfrac{N}{p}-1}_{p,1})}.$$
Therefore $\|\d_X\nabla P\|_{L_t^1(\hB^{\Psfrac{N}{p}+\ep-2}_{p,1})}$
may be bounded like $\sK(t)$.

\subsection{The regularization process}

In all the above computations, we implicitly assumed that
$X$ and $\d_Xu$ were in $L^\infty_\mathrm{loc}(\R_+;\cC^{0,\ep})$ and
$L^1_\mathrm{loc}(\R_+;\cC^{0,\ep})$, respectively.
However, Theorem \ref{thm:wellposed_INS} just ensures continuity of those vector-fields, not
Hölder regularity.

To overcome that difficulty, one may smooth out the initial velocity (not the density, not to
destroy the multiplier hypotheses) by setting for example $u_0^n:=\dot S_nu_0$.
Then the condition \eqref{cdt:u0_small} is satisfied by $(\rho_0,u_0^n)$
and, as in addition
$u_0^n$ belongs to all Besov spaces $\dot B^{\Psfrac Np-1}_{\wt p,r}$ with $\wt p\geq p$ and $r\geq1$,
one can
apply\footnote{That paper concerns the half-space; having the same result
in the whole space setting is much easier.}
\cite[Th.\,1.1]{DanZhP2014} for solving (INS) with initial data $(\rho_0,u_0^n)$.
This provides us with a unique
global solution $(\rho^n,u^n,\nabla P^n)$ which, among others, satisfies
$$\nabla u^n\in L^r(\R_+;\dot B^{\sfrac Np}_{\wt p,r})
\quad\hbox{for all }\ r\in{}]1,\infty[\ \hbox{ and }\
\max\Bigl(p,\frac{Nr}{3r-2}\Bigr)\leq \wt p \leq \frac {Nr}{r-1}\cdotp$$
By taking $r$ sufficiently close to $1$ and using embedding, we see that
this implies that $\nabla u^n$ is in $L^1_\mathrm{loc}(\R_+;\hsC^{0,\delta})$
for all $0<\delta<1$ and thus
the corresponding flow $\psi^n$ is (in~particular) in $\cC^{1,\ep}$. This ensures, thanks to
\eqref{eq:Xt}, that $X^n$ is in $L^\infty_\mathrm{loc}(\R_+;\cC^{0,\ep})$ and thus that
$\d_{X^n}u^n$ is in $L^1_\mathrm{loc}(\R_+;\cC^{0,\ep})$.

{}From the previous steps and the fact that the data $(\rho_0,u_0^n)$ satisfy \eqref{cdt:u0_small} uniformly,
we get uniform bounds for $\rho^n$, $u^n$, $\nabla P^n$ and $X^n$, and standard arguments
thus allow to show that $u^n$ tends to $u$ in $L^1_\mathrm{loc}(\R_+;L^\infty)$ and
thus $(\psi^n-\psi)\to0$ in $L^\infty_\mathrm{loc}(\R_+;L^\infty)$.
Interpolating with the uniform bounds and using standard functional analysis
arguments, one can conclude that
$X^n\to X$ in $L^\infty_\mathrm{loc}(\R_+;\cC^{0,\ep'})$ for all $\ep'<\ep$
(and similar results for $(u^n)_{n\in\N}$) and that all the estimates of the previous steps are satisfied.
The details are left to the reader.\qed

\begin{appendix}

\section{Multiplier spaces}

The following relationship between the nonhomogeneous Besov spaces $B^{s}_{p,r}(\R^N)$
and the homogeneous Besov spaces $\hB^{s}_{p,r}(\R^N)$ for
\emph{compactly supported} functions or distributions has been established in \cite[\S2.1]{DanM2015}.
\begin{prop}\label{prop:equiv_besov}
Let $(p,r) \in [1,\infty]^2$ and $s > -\sfrac{N}{p'}:= -N(1-\sfrac{1}{p})$ (or just $s\geq-\sfrac N{p'}$ if $r=\infty$).
For any $u$ in the set $\cE'(\R^N)$ of compactly supported distributions on $\R^N$, we have
\begin{equation*}
u \in B^{s}_{p,r}(\R^N) \Longleftrightarrow u\in \hB^{s}_{p,r}(\R^N).
\end{equation*}
Moreover, there exists a constant $C=C(s,p,r,N, \mathrm{Supp}\, u)$ such that
\begin{equation*}
C^{-1}\|u\|_{\hB^{s}_{p,r}} \leq \|u\|_{B^{s}_{p,r}} \leq C \|u\|_{\hB^{s}_{p,r}}.
\end{equation*}
\end{prop}
A simple consequence of Proposition \ref{prop:equiv_besov} and of standard embeddings for nonhomogeneous Besov spaces is that for any $(s,p,r)$ as above, we have
\begin{equation}\label{es:nonhom_hom}
\cE'(\R^N) \cap \hB^{s+\delta}_{p,r}(\R^N) \hto \cE'(\R^N) \cap \hB^{s}_{p,r}(\R^N) \quad \mbox{for any } ~\delta>0.
\end{equation}
We also used the following statement:
\begin{prop}\label{prop:besov-holder}
Let $(s,p,r)$ be arbitrary in $\R\times[1,\infty]^2$.
Then for all $u\in B^s_{\infty,r}(\R^N)\cap\cE'(\R^N)$,
we have $u\in B^{s}_{p,r}(\R^N)$ and there exists $C=C(s,p, \mathrm{Supp}\,u)$ such that
$$
\|u\|_{B^s_{p,r}}\leq C\|u\|_{B^s_{\infty,r}}.
$$
\end{prop}
\begin{proof}
Let $u$ be in $B^s_{\infty,r}(\R^N)$ with compact support. Fix some
smooth cut-off function~$\phi$ so that $\phi\equiv1$ on $\mathrm{Supp}\, u$.
Being compactly supported and smooth, $\phi$ belongs to any nonhomogeneous Besov space.
Then, using (the nonhomogeneous version of) the decomposition \eqref{eq:bony} and that $u=\phi u$, we get
$$
u=T_\phi u+ T_u\phi+R(u,\phi).
$$
Because $\phi$ is in $L^p$ and $u$ in $B^s_{\infty,r}$, standard continuity results for the paraproduct
ensure that $T_\phi u$ is in $B^s_{p,r}$. For the second term, we just use that
$u$ is in $\sC^{-|s|-1}$ and $\phi$ in $B^{|s|+1+s}_{p,r}$ hence
$T_u\phi$ is in $B^s_{p,r}$. For the remainder term, we use for instance the fact that $\phi$ is in $\sC^{|s|+1}$.
Putting all those information together completes the proof.
\end{proof}

The following result was the key to bounding the density terms in our study of \eqref{eq:INSmu}.
\begin{lemma}\label{lemma:multiplier}
Let $(s,s_k,p,p_k,r, r_k)\!\in{}]-1,1[^2 \times [1,\infty]^4$ with $k\!=\!1,2$, and \hbox{$Z\!:\!\R^N\!\to\!\R^N$} be a $\cC^1$ measure preserving diffeomorphism such that $DZ$ and $DZ^{-1}$ are bounded. When we consider the homogeneous Besov space $\hB^s_{p,r}(\R^N)$ or $\hB^{s_k}_{p_k,r_k}(\R^N)$, we assume in addition that $s\in{}]-\sfrac N{p'},\sfrac Np[$ and $s_k\in{}]-\sfrac N{p'_k},\sfrac N{p_k}[$ for
$k=1,2$. Then we have:
\begin{enumerate}
\item If $b^{s}_{p,r}(\R^N)$ stands for $B^{s}_{p,r}(\R^N)$ or $\hB^{s}_{p,r}(\R^N)$, then the mapping $u \mto u\circ Z $ is continuous on $b^{s}_{p,r}(\R^N)$:
there is a positive constant $C_{Z,s,p,r}$ such that
\begin{equation}\label{es:composition_besov}
\|u\circ Z \|_{b^{s}_{p,r}} \leq C_{Z,s,p,r} \|u \|_{b^{s}_{p,r}}.
\end{equation}

\item If $b^{s_k}_{p_k, r_k}$ with $k=1,2$, denote the same type of Besov spaces, then the mapping $\varphi \mto \varphi \circ Z $ is continuous on $\cM\big(b^{s_1}_{p_1,r_1}(\R^N) \to b^{s_2}_{p_2,r_2}(\R^N) \big)$, that is
\begin{equation*}
\|\varphi \circ Z\|_{\cM(b^{s_1}_{p_1,r_1}\to b^{s_2}_{p_2,r_2})} \leq
C_{Z^{-1},s_1,p_1,r_1} C_{Z,s_2,p_2,r_2} \|\varphi\|_{\cM(b^{s_1}_{p_1,r_1} \to b^{s_2}_{p_2,r_2})}.
\end{equation*}
\item We have the following equivalence for any $\varphi \in \cE'(\R^N)$,
\begin{equation*}
\varphi \in \cM\big(B^{s_1}_{p_1,r_1}(\R^N) \to B^{s_2}_{p_2,r_2}(\R^N) \big)
\Longleftrightarrow
\varphi \in \cM\big(b^{s_1}_{p_1,r_1}(\R^N) \to b^{s_2}_{p_2,r_2}(\R^N) \big) .
\end{equation*}
Here $b^{s_1}_{p_1,r_1}$ and $b^{s_2}_{p_2,r_2}$ can be different type of Besov spaces but obey our convention on the index $s_k$ for homogeneous Besov space.
\end{enumerate}
\end{lemma}

\begin{proof}
Item (i) in the case $b= \hB$ has been proved in \cite[Lem.\,2.1.1]{DanM2015}.
One may easily modify the proof to handle nonhomogeneous Besov spaces:
use the finite difference characterization of \cite[p.\,98]{Tri1992} if $s>0$, argue
by duality if $s<0$ and interpolate for the case $s=0$. We get
$C_{Z,s,p,r}\approx 1+ \|DZ\|^{s+\sfrac{N}{r}}_{L^{\infty}}$ if $s>0$, and $C_{Z,s,p,r}\approx 1+ \|DZ^{-1}\|^{-s+\sfrac{N}{r'}}_{L^{\infty}}$ if $s<0$.

Part (ii) is immediate according to \eqref{eq:def_norm_multiplier} and \eqref{es:composition_besov}. Indeed we may write:
$$
\begin{aligned}
\|\varphi \circ Z\|_{\cM\big(b^{s_1}_{p_1,r_1}\to b^{s_2}_{p_2,r_2} \big)}
&=\sup_{\|u\|_{b^{s_1}_{p_1,r_1}}\leq 1}
\|(\varphi\circ Z) \, u\|_{b^{s_2}_{p_2,r_2}}\\[-3pt]
&=\sup_{\|u\|_{b^{s_1}_{p_1,r_1}}\leq 1}
\|(\varphi\,(u\circ Z^{-1}))\circ Z\|_{b^{s_2}_{p_2,r_2}}\\[-3pt]
&\leq C_{Z,s_2,p_2,r_2} \sup_{\|u\|_{b^{s_1}_{p_1,r_1}}\leq 1}
\|\varphi\,(u\circ Z^{-1})\|_{b^{s_2}_{p_2,r_2}}\\[-3pt]
&\leq C_{Z,s_2,p_2,r_2} \|\varphi\|_{\cM(b^{s_1}_{p_1,r_1}\to b^{s_2}_{p_2,r_2})}
\sup_{\|u\|_{b^{s_1}_{p_1,r_1}}\leq 1}\|u\circ Z^{-1}\|_{b^{s_1}_{p_1,r_2}}\\[-3pt]
&\leq C_{Z^{-1},s_1,p_1,r_1} C_{Z,s_2,p_2,r_2}
\|\varphi\|_{\cM(b^{s_1}_{p_1,r_1}\to b^{s_2}_{p_2,r_2})}.\end{aligned}
$$
To prove the last item, it suffices to check that if $\varphi$ belongs to \hbox{$\cE'\!\cap\!\cM(B^{s_1}_{p_1,r_1}\!\to\!B^{s_2}_{p_2,r_2})$}, then $\varphi$ is also in the multiplier space between the general type Besov spaces. Take \hbox{$u \in b^{s_1}_{p_1,r_1}$} with compact support, and some smooth and compactly supported nonnegative
cut-off function $\psi$ satisfying $\psi \equiv 1$ on $\mathrm{Supp}\,\varphi$. Then from Proposition \ref{prop:equiv_besov} and \eqref{eq:def_norm_multiplier}, we have
\begin{align*}
\|\varphi u\|_{b^{s_2}_{p_2,r_2}} =\|\varphi \psi u\|_{b^{s_2}_{p_2,r_2}} \lesssim \|\varphi \psi u\|_{B^{s_2}_{p_2,r_2}} &\lesssim \|\varphi\|_{\cM(B^{s_1}_{p_1,r_1} \to B^{s_2}_{p_2,r_2} )} \| \psi u\|_{B^{s_1}_{p_1,r_1}}\\
&\lesssim \|\varphi\|_{\cM(B^{s_1}_{p_1,r_1} \to B^{s_2}_{p_2,r_2} )} \|\psi u\|_{b^{s_1}_{p_1,r_1}}\\
&\lesssim \|\varphi\|_{\cM(B^{s_1}_{p_1,r_1} \to B^{s_2}_{p_2,r_2} )} \|\psi\|_{\cM(b^{s_1}_{p_1,r_1} )}\| u\|_{b^{s_1}_{p_1,r_1}}.
\end{align*}
For the last inequality, we used $\cC_c^\infty \hto \cM(b^{s_1}_{p_1,r_1})$ (see \cite[Cor.\,2.1.1]{DanM2015}).
\end{proof}

\section{Commutator Estimates}

We here recall and prove some commutator estimates that were crucial in this paper. All of them strongly rely on continuity results in Besov spaces for the paraproduct and remainder operators, and on the following classical result (see \eg \cite[\S2.10]{BCD2011}).
\begin{lemma}\label{lemma:ce1}
Let $A:\R^N\setminus\{0\}\to\R$ be a smooth function, homogeneous of degree~$m$. Let $(\ep,s,p,r,r_1,r_2,p_1,p_2) \in{}]0,1[{}\times \R \times [1, \infty]^6$ with $\frac{1}{p}=\frac{1}{p_1}+\frac{1}{p_2}$,
$\frac1r=\frac1{r_1}+\frac1{r_2}$ and
$$s-m+\ep <\sfrac{N}{p} \quad \mbox{or} \quad\{s-m+\ep <\sfrac{N}{p} ~~\mbox{and}~~ r=1\}\cdotp$$

There exists a constant $C$ depending only on $s, \ep, N$ and $A$ such that,
$$\|[\hT_g, A(D)]u\|_{\hB^{s-m+\ep}_{p,r}} \leq C \|\nabla g\|_{\hB^{\ep-1}_{p_1,r_1}}\|u\|_{\hB^{s}_{p_2,r_2}}.$$
\end{lemma}

If the integer $N_0$ in the definition of Bony's paraproduct and remainder is large enough (for instance $N_0 =4$ does), then the following fundamental lemma holds.
\begin{lemma}[Chemin-Leibniz Formula]\label{lemma:chemin_leibniz}
Let \hbox{$(\ep,s_k,p, p_k,p_3, r,r_k)\!\in\!{}]0,1[{}\!\times\!\R\!\times\![1, \infty]^5$} for $k=1,2$ satisfying
\begin{equation*}
\frac{1}{p} = \frac{1}{p_1} + \frac{1}{p_2} +\frac{1}{p_3}\quad\text{and}\quad \frac{1}{r} = \frac{1}{r_1} + \frac{1}{r_2}\cdotp
\end{equation*}

\begin{enumerate}
\item If $s_2 < 0$ and $s_1+s_2+\ep-1 < \sfrac{N}{p}$ or $\{s_1+s_2+\ep-1=\sfrac{N}{p} ~~\mbox{and}~~r=1\}$, then we have
\begin{equation*}
\|\hcT_X \hT_g f - \hT_g \hcT_X f - \hT_{ \hcT_X g}f \|_{\hB^{s_1+s_2+\ep-1}_{p,r}} \leq C \|X\|_{\hB^{\ep}_{p_3,\infty}} \|g\|_{\hB^{s_2}_{p_2,\infty}}
\|f\|_{\hB^{s_1}_{p_1,r}}.
\end{equation*}
The above inequality still holds for $s_2=0$, if one replaces $\|g\|_{\hB^{0}_{p_2,\infty}}$ by $\|g\|_{L^{p_2}}$.
\item If $s_1+s_2+\ep-1 \in{}]0, \sfrac{N}{p}[$ or $\{s_1+s_2+\ep-1=\sfrac{N}{p} ~~\mbox{and}~~r=1\}$, then we have
\begin{equation*}
\|\hcT_X \hR(f,g)-\hR(\hcT_X f, g)-\hR(f, \hcT_X g)\|_{\hB^{s_1+s_2+\ep-1}_{p,r}} \leq C \|X\|_{\hB^{\ep}_{p_3,\infty}} \|f\|_{\hB^{s_1}_{p_1,r_1}}\|g\|_{\hB^{s_2}_{p_2,r_2}}.
\end{equation*}
The above inequality still holds for $s_1+s_2+\ep -1 =0$, $r=\infty$ and $\frac{1}{r_1} + \frac{1}{r_2}=1$.
\end{enumerate}
\end{lemma}

\begin{proof} This is a mere adaptation of \cite{DanZhX2017} to the homogeneous framework.
The proof is based on a generalized Leibniz formula for para-vector field operators which was derived by
{J.-Y.\,Chemin} in \cite{Che1988}. More precisely, define the following Fourier multipliers
\begin{equation*}
\hDelta_{k,j}:= \varphi_{k}(2^{-j}D) \quad\ \,\hbox{ with }\, \varphi_{k}(\xi):= i \xi_k \varphi(\xi) \,\,
\hbox{ for }\ k \in\{1,\dots,N\}\ \hbox{ and } \,\, j \in \Z.
\end{equation*}
Then we have
\begin{align*}
\hcT_X \hT_g f & = \sum_{j\in\Z} (\hS_{j-N_0} g \hcT_{X}\hDelta_{j} f + \hDelta_j f \hcT_X \hS_{j-N_0} g) + \sum_{j\in\Z}(\hT_{1,j} + \hT_{2,j})\\
&= \hT_g \hcT_X f +\hT_{ \hcT_X g}f + \sum_{\substack{j\in\Z\\\alpha=1,\ldots,4}} \hT_{\alpha, j},
\end{align*}
where
\begin{align*}
\hT_{1,j} &:= \sum_{\substack{j \leq j' \leq j+1\\j-N_0-1 \leq j''\leq j'-N_0-1}} 2^{j'} \hDelta_{j''}X^k \big(\hDelta_{k,j'} (\hDelta_j f \hS_{j-N_0} g)-\hDelta_{k,j'}\hDelta_j f \hS_{j-N_0}g \big),\\
\hT_{2,j} &:= \sum_{\substack{j'\leq j- 2\\j'-N_0 \leq j''\leq j-N_0 -2}} 2^{j'} \hDelta_{j''}X^k (\hDelta_j f) \hDelta_{k,j'} \hS_{j-N_0} g,\\
\hT_{3,j} &:= \hS_{j-N_0} g [\hT_{X^k}, \hDelta_j] \d_k f,\\
\hT_{4,j} &:= \hDelta_j f [\hT_{X^k}, \hS_{j-N_0}] \d_k g.
\end{align*}
Bounding $\hT_{1,j}$ and $\hT_{2,j}$ stems from the definition
of Besov norms, and Lemmas 2.99, 2.100 of \cite{BCD2011} allow to bound $\hT_{3,j}$ and $\hT_{4,j}$
provided $\ep<1$.

In order to prove the second item, let us set $$A_{j,j'}:= \big\{j-N_0-1,\dots,
j'-N_0-1\big\} \cup \big\{j'-N_0,\dots, j-N_0-2\big\}\cdotp$$
We have
\begin{align*}
\hcT_X \hR(f, g) & = \sum_{j\in\Z} (\thDelta_{j}g \hcT_{X}\hDelta_{j} f + \hDelta_j f \hcT_X \thDelta_{j} g) + \sum_{j\in\Z}(\hR_{1,j} + \hR_{2,j})\\
&= \hR(\hcT_X f, g)+\hR(f, \hcT_X g) + \sum_{\substack{j\in\Z\\\alpha=1,\ldots,4}} \hR_{\alpha, j},
\end{align*}
where, denoting $\thDelta_j:=\hDelta_{j-N_0}+\cdots+\hDelta_{j+N_0}$,
\begin{align*}
\hR_{1,j} &:= \sum_{\substack{|j'-j| \leq N_0 + 1\\j''\in A_{j,j'}}}\hspace*{-5mm}\mathrm{sgn}(j'-j+1) 2^{j'} \hDelta_{j''}X^k \big(\hDelta_{k,j'} (\hDelta_j f \thDelta_{j} g) - \hDelta_j f \hDelta_{k,j'} \thDelta_{j}g \big)\\[-20pt]
& \hspace{5cm}+ \sum_{\substack{j-1 \leq j' \leq j\\j'-N_0 \leq j'' \leq j-N_0}}\hspace*{-5mm} 2^{j'} \hDelta_{j''}X^k(\hDelta_{k,j'} \hDelta_j f) \thDelta_{j} g,\\
\hR_{2,j} &:= \sum_{\substack{j' \leq j-N_0-2\\j'-N_0 \leq j''\leq j-N_0-2}}\hspace*{-5mm} 2^{j'} \hDelta_{j''}X^k \hDelta_{k,j'} (\hDelta_j f \thDelta_{j} g) ,\\
\hR_{3,j} &:= \thDelta_j g [\hT_{X^k}, \hDelta_j] \d_k f,\\
\hR_{4,j} &:= \hDelta_j f [\hT_{X^k}, \thDelta_j] \d_k g.
\end{align*}
Here again, bounding $\hR_{1,j}$ and $\hR_{2,j}$ follows from the definition of
Besov norms, while Lemma 2.100 of \cite{BCD2011} allows to bound $\hR_{3,j}$ and $\hR_{4,j}$.
\end{proof}

\begin{prop}\label{prop:ce_dINS}
Let $(\ep,p)$ be in $]0,1[{}\times [1,\infty]$.
Consider a pair of vector fields $(X,v)$ in
$$\bigl(L^{\infty}_\mathrm{loc}(\R_+; \hsC^{\ep}) \bigr)^N \times \big(L^{\infty}_\mathrm{loc}(\R_+ ; \hB^{\Psfrac{N}{p}-1}_{p,1})\cap L^{1}_\mathrm{loc}(\R_+ ;\hB^{\Psfrac{N}{p}+1}_{p,1})\big)^N,$$
satisfying $\div v=0$ and the transport equation
\begin{equation}\label{eq:prop_ce_dINS}
\begin{cases}
(\d_t+v\cdot \nabla)X = \d_X v,\\[3pt]
X|_{t=0} = X_0.
\end{cases}
\end{equation}
If in addition
\begin{equation}\label{cdt:p_Xbis}
\sfrac Np > 2-\ep, \quad\hbox{or }\ \sfrac Np >1-\ep \: \hbox{ and }\: \div X\equiv0,
\end{equation}
then there exists a constant $C$ such that:
\begin{multline}\label{es:TX}
\|[\hcT_X, \d_t+v\cdot \nabla ]v\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \leq C(\|X\|_{\hsC^\ep} \|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}} \|v\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \\+\|v\|_{\hsC^{-1}}\|\hcT_{X}v\|_{\hB^{\Psfrac{N}{p}+\ep}_{p,1}}+\|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}\|\hcT_{X}v\|_{\hsC^{\ep-2}}).
\end{multline}
\end{prop}

\begin{proof}
This is essentially the proof of \cite[Prop.\,A.5]{DanZhX2017}. For the reader convenience, we here give a sketch of it.
Because $\div v=0$, we may write
\begin{align*}
[\hcT_X, \d_t+v^{\ell} \d_\ell ]v & = -v^\ell\d_\ell \hT_{X^k}\d_kv- \hT_{\d_t X^k}\d_kv
+ \hT_{X^k}\d_k(v^{\ell}\d_\ell v)\\
&=- \hT_{\d_t X^k}\d_kv + \d_\ell\hcT_{X}(v^{\ell} v)-\hcT_{\d_\ell X}(v^{\ell} v) - v^{\ell} \d_{\ell} \hcT_Xv.
\end{align*}
Hence, decomposing $v^{\ell} v$ according to Bony's decomposition, we discover that
$$
[\hcT_X, \d_t+v^{\ell} \d_\ell ]v = \sum_{\alpha=1}^{\alpha=5} \hR_{\alpha}
$$
\[
\tag*{with}
\begin{aligned}
\hR_1 &:=- \hT_{\d_t X^k}\d_kv, &\hR_2 &:= \d_{\ell}(\hcT_X \hT_{v^{\ell}}v +\hcT_X \hT_{v}v^{\ell}),&
\hR_3 &:= \d_{\ell} \hcT_{X} \hR(v^{\ell},v),\\
\hR_4&:=-\hcT_{\d_\ell X}(v^{\ell} v),&
\hR_5 &:=- v^{\ell} \d_{\ell} \hcT_Xv.
\end{aligned}
\]

It suffices to check that all the terms $\hR_\alpha$ may be bounded by the right-hand side of \eqref{es:TX}.

\subsubsection*{Bound of $\hR_1$}

From the equation \eqref{eq:prop_ce_dINS}, we have
\begin{equation*}
\hR_1 = \hT_{v\cdot \nabla X^k} \d_k v - \hT_{\d_X v^k} \d_k v.
\end{equation*}
Hence using standard continuity results for the paraproduct, we deduce that
$$
\|\hR_1\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|\nabla v\|_{\hB^{\sfrac{N}{p}}_{p,1}}\bigl(\|v \cdot \nabla X\|_{\hsC^{\ep-2}} + \|\d_X v\|_{\hsC^{\ep-2}}\bigr).
$$
Keeping in mind \eqref{cdt:p_Xbis}, the last term may be bounded according to \eqref{es:dXf_TXf},
after using the embedding $\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}(\R^N) \hto \hsC^{\ep-2}(\R^N)$.
We get
\begin{equation*}
\|\d_X v- \hcT_X v\|_{\hsC^{\ep-2}} \lesssim \|\nabla v\|_{\hB^{\Psfrac{N}{p}-2}_{p,1}} \|X\|_{\hsC^\ep}.
\end{equation*}

As for the first term, we use the fact $\div v=0$ and the following decomposition
$$
v\cdot\nabla X=\hcT_{v} X + \hT_{\d_{\ell} X } v^\ell +\d_\ell \hR (v^{\ell},X),
$$
which allow to get, as long as \eqref{cdt:p_Xbis} holds
\begin{equation*}
\|\hR_1\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|\nabla v\|_{\hB^{\sfrac{N}{p}}_{p,1}}\bigl(\|v\|_{\hB^{\Psfrac{N}{p}-1}_{p,1}} \|X\|_{\hsC^{\ep}} + \|\hcT_X v\|_{\hsC^{\ep-2}}\bigr).
\end{equation*}
\subsubsection*{Bound of $\hR_2$}

Due to Lemma \ref{lemma:chemin_leibniz} (i) and continuity of paraproduct operator, we have
\begin{multline*}
\|\hR_2\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|X\|_{\hsC^{\ep}} \|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}} \|v\|_{\hsC^{-1}}\\
+\|v\|_{\hsC^{-1}} \|\hcT_{X}v\|_{\hB^{\Psfrac{N}{p}+{\ep}}_{p,1}}+\|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}\|\hcT_{X}v\|_{\hsC^{\ep-2}}.
\end{multline*}

\subsubsection*{Bound of $\hR_3$}
Applying Lemma \ref{lemma:chemin_leibniz} (ii) and continuity of remainder operator under the condition
$\Psfrac{N}{p} +\ep-1>0$ yields
\begin{equation*}
\|\hR_3\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|X\|_{\hsC^{\ep}} \|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}\|v\|_{\hsC^{-1}} + \|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}} \|\hcT_X v\|_{\hsC^{\ep-2}}.
\end{equation*}

\subsubsection*{Bound of $\hR_4$}
{}From Bony decomposition \eqref{eq:bony}, it is easy to get
\begin{equation*}
\|v^l v\|_{\hB^{\sfrac{N}{p}}_{p,1}} \lesssim \|v\|_{\hsC^{-1}}\|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}.
\end{equation*}
Hence
\begin{equation*}
\|\hR_4\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|\nabla X\|_{\hsC^{\ep-1}}\|v\|_{\hsC^{-1}}\|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}.
\end{equation*}

\subsubsection*{Bound of $\hR_5$}
Applying Bony decomposition and using that $\div v\!=\!0$ and~\hbox{$\frac{N}{p}\!+\!\ep\!>\!1$} give
\begin{equation*}
\|\hR_5\|_{\hB^{\Psfrac{N}{p}+\ep-2}_{p,1}} \lesssim \|v\|_{\hsC^{-1}}\|\hcT_{X}v\|_{\hB^{\Psfrac{N}{p}+\ep}_{p,1}}+\|v\|_{\hB^{\Psfrac{N}{p}+1}_{p,1}}\|\hcT_{X}v\|_{\hsC^{\ep-2}}.
\end{equation*}
Combining the above estimates for all $\hR_{\alpha}$, with $\alpha=1,\dots,5$ yields \eqref{es:TX}.
\end{proof}

Another consequence of Lemma \ref{lemma:chemin_leibniz} is the following estimate of $\div (X f g)$:
\begin{prop} \label{prop:dXfg_ND}
Let $(s,p,r)$ be in $]0,1[{}\times [1,\infty]^2$ and $\eta$ in $]0,1-s[$.
Consider a bounded vector field $X$ and two bounded functions $f,g$ satisfying
\begin{equation*}
X \in \big(\hB^s_{p,r} (\R^N) \cap \sC^{s+\eta} \big)^{N},\,\,\, (f, g)\in \hB^s_{p,r} (\R^N) \times \hB^{-\eta}_{p, r} (\R^N) \,\,\,\hbox{and}\,\,\, \d_X g \in \hB^{s-1}_{p,r} (\R^N).
\end{equation*}
If in addition $\div X $ belongs to $\cM\big(\hB^{s}_{p,r}(\R^N) \to \hB^{s-1}_{p,r}(\R^N)\big)$,
and there exists some $q\in[1,p[$
such that
\begin{equation}\label{cdt:s_pq}
\div X \in \hB^{s_{p,q}}_{q,r}(\R^N)\ \hbox{ with }\ s_{p,q}:= s-1 + N(\sfrac{1}{q} - \sfrac{1}{p})>0,
\end{equation}
then we have $\div (X f g) \in \hB^{s-1}_{p,r}(\R^N)$, and the following estimate holds true:
\begin{multline*}
\|\div (X f g)\|_{\hB^{s-1}_{p,r}}\lesssim \|X\|_{\hB^s_{p,r} \cap \sC^{s+\eta}} \|f\|_{ L^\infty \cap \hB^s_{p,r}}\|g\|_{L^\infty \cap \hB^{-\eta}_{p,r}}
+ \|f\|_{L^\infty} \|\d_X g\|_{\hB^{s-1}_{p,r}} \\
+ \|\div X\|_{\hB^{s_{p,q}}_{q,r} \cap \cM(\hB^{s}_{p,r}\to \hB^{s-1}_{p,r})} \|g\|_{L^\infty} \|f\|_{\hB^{s}_{p,r} \cap L^\infty}.
\end{multline*}
\end{prop}
\begin{proof}
In light of Bony's decomposition \eqref{eq:bony}, and denoting
$\hT'_g f:= \hT_g f + \hR(f,g)$, we can decompose $\div (X fg)$ into
\begin{equation*}
\div (Xfg) = \div \big(\hT'_{fg} X + \hT_X (fg) \big) = \sum_{\alpha=1}^{4} \hF_{\alpha},\text{ where }
\begin{cases}
\hF_1\!:=\!\div (\hT'_{fg} X),&\hspace*{-5pt}\hF_3\!:=\!\hcT_X \hT'_g f,\\
\hF_2\!:=\!\hT_{\div X} (fg),&\hspace*{-5pt}\hF_4\!:=\!\hcT_X \hT_f g.
\end{cases}
\end{equation*}

\subsubsection*{Bound of $\hF_1$} As $s>0$, standard
continuity results for $\hT$ and $\hR$ yield
\begin{equation*}
\|\hF_1\|_{\hB^{s-1}_{p,r}} \lesssim \|\hT'_{fg} X\|_{\hB^{s}_{p,r}} \lesssim \|f\|_{L^\infty}\|g\|_{L^\infty} \|X\|_{\hB^{s}_{p,r}}.
\end{equation*}

\subsubsection*{Bound of $\hF_2$}
Thanks to continuity results for $\hT$, we have for $s <1$,
\begin{equation*}
\|\hF_2\|_{\hB^{s-1}_{p,r}} \lesssim \| \div X\|_{\hB^{s-1}_{p,r}} \|f\|_{L^\infty}\|g\|_{L^\infty}.
\end{equation*}
\subsubsection*{Bound of $\hF_3$}
Because $X$ and $g$ are bounded and $s>0$,
we readily have
\begin{equation*}
\|\hF_3\|_{\hB^{s-1}_{p,r}} \lesssim \|X\|_{L^\infty} \|\hT'_g f\|_{\hB^{s}_{p,r}} \lesssim \|X\|_{L^\infty}\|g\|_{L^\infty} \|f\|_{\hB^{s}_{p,r}}.
\end{equation*}
\subsubsection*{Bound of $\hF_4$}
Because $0<s<s+\eta<1$,
Lemma \ref{lemma:chemin_leibniz} and continuity results for the paraproduct $\hT$ imply that
$$\begin{aligned}
\|\hcT_X \hT_f g\|_{\hB^{s-1}_{p,r}} &\lesssim
\|X\|_{\hsC^{s+\eta}} \|f\|_{L^\infty}\|g\|_{\dot B^{-\eta}_{p,r}}
+ \| \hT_f \hcT_X g \|_{\hB^{s-1}_{p,r}}
+\| \hT_{ \hcT_X f}g \|_{\hB^{s-1}_{p,r}}\\
&\lesssim \|X\|_{\hsC^{s+\eta}} \|f\|_{L^\infty}\|g\|_{\dot B^{-\eta}_{p,r}}
+ \|f\|_{L^\infty} \|\hcT_X g\|_{\hB^{s-1}_{p,r}}
+ \|g\|_{L^\infty} \|\hcT_X f\|_{\hB^{s-1}_{p,r}}.
\end{aligned}$$
To bound the last term, one may use the decomposition
$$\hcT_X f=\div(\hT_Xf)-f\div X+\hT_f\,\div X+\hR(f,\div X).$$
Hence using continuity results for $\hR$ and $\hT$ and the fact that
$(s_{p,q},q)$ satisfies \eqref{cdt:s_pq},
$$
\|\hcT_X f\|_{\dot B^{s-1}_{p,r}}\lesssim
\|f\|_{\dot B^s_{p,r}}\bigl(\|X\|_{L^\infty}+ \|\div X\|_{\cM(\hB^{s}_{p,r}\to \hB^{s-1}_{p,r})}\bigr)
+\|f\|_{L^\infty}\|\div X\|_{ \hB^{s_{p,q}}_{q,r}}.
$$
Finally, to bound the term with $\hcT_Xg$, we use the fact that
$$\d_X g - \hcT_X g = \hT_{\nabla g}\cdot X + \div \hR(X, g) - \hR(\div X, g),$$
whence
\begin{equation}\label{es:dXh}
\|\d_X g - \hcT_X g\|_{\hB^{s-1}_{p,r}} \lesssim \|g\|_{L^\infty} \big(\|X\|_{\hB^s_{p,r}} + \|\div X\|_{\hB^{s_{p,q}}_{q,r}} \big).
\end{equation}
This completes the proof of the proposition.
\end{proof}
Proposition \ref{prop:dXfg_ND} above reveals that the bounded function $g$ may behave like some element in $\cM(\hB^{s-1}_{p,\infty})$ under a suitable additional structure assumption.
If in addition~$g$ has compact support, then one can relax a bit the regularity of $X$ and $f$ to study $\d_X (fg)$,
and get the following generalization of \cite[Lem.\,A.6]{Dan1999}.
\begin{cor} \label{cor:dXfg_3D}
Consider a divergence-free vector field $X$ with coefficients in $\sC^{\ep}$,
and some function $f$ in $\sC^{\ep'}$ with $0<\ep,\ep'<1$.
Let $g \in L^\infty$ be compactly supported and satisfy $\d_X g \in \dot B^{\alpha-1}_{p,r}$ for some $(p,r) \in [1,\infty]^2$
and $\alpha \in{}]0,\min\{\ep,\ep'\}[$.
Then $\div (Xfg) = \d_X (fg) \in \dot B^{\alpha-1}_{p,r}$.
\end{cor}
\begin{proof}
Let $\psi \in \cC^{\infty}_c$ be a cut-off function such that
$\psi\equiv 1$ near $\mathrm{Supp}\, g$. Denote $(\tX, \tf):= (\psi X, \psi f)$.
From Proposition \ref{prop:equiv_besov} and the proof of Proposition \ref{prop:besov-holder}, we know that
\begin{equation*}
(\tX,\tf,g) \in (B^{\ep}_{q,\infty})^N \times B^{\ep'}_{q,\infty} \times B^{-\eta}_{q,r}
\hto (\hB^{\alpha}_{q,1}\cap L^\infty)^{N+1} \times \hB^{-\eta}_{q,r},
\end{equation*}
for any $q \in [1,\infty]$ and some $\eta \in{}]0, \min\{N\slash q' ,\ep-\alpha\}[$.
It is also clear that $\d_X (fg) = \div (\tX \tf g)$ and $\d_{\tX} g =\d_X g$.
Hence applying Proposition \ref{prop:dXfg_ND}
gives the result.
\end{proof}
\end{appendix}

\backmatter
\bibliographystyle{jepplain}
\bibliography{danchin-zhang}
\end{document}
