# Quantitative derivation of a two-phase porous media system from the one-velocity Baer-Nunziato and Kapila systems

Timothée Crin-Barat, Ling-Yun Shou, Jin Tan

► To cite this version:

Timothée Crin-Barat, Ling-Yun Shou, Jin Tan. Quantitative derivation of a two-phase porous media system from the one-velocity Baer-Nunziato and Kapila systems. Nonlinearity, 2024, 37 (7), pp.075002. 10.1088/1361-6544/ad3f66. hal-04862724

HAL Id: hal-04862724

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

Submitted on 10 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Quantitative derivation of a two-phase porous media system from the one-velocity Baer–Nunziato and Kapila systems

To cite this article: Timothée Crin-Barat et al 2024 *Nonlinearity* **37** 075002

View the [article online](#) for updates and enhancements.

## You may also like

- *Diffraction Measurement for Grating Formed in Optically Poled Polymeric Materials*
Naoto Tsutsumi, Jun Yamamoto and Wataru Sakai

- *Polymer thin-film distributed feedback tunable lasers*
Vincent Dumarcher, Licinio Rocha, Christine Denis et al.

- *One- and two-photon picosecond excitation dynamics of the singlet states of a tetraphenyl-diamine derivative in solution*
Barbara Paci, Chantal Andraud, Remi Anemian et al.
# Quantitative derivation of a two-phase porous media system from the one-velocity Baer–Nunziato and Kapila systems

Timothée Crin-Barat¹*, Ling-Yun Shou² and Jin Tan³

¹ Chair for Dynamics, Control, Machine Learning and Numerics, Alexander Von Humboldt- Professorship, Department of Mathematics, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany

² School of Mathematics and Key Laboratory of Mathematical MIIT, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, People's Republic of China

³ Laboratoire de Mathématiques AGM, UMR CNRS 8088, Cergy Paris Université, 2 Avenue Adolphe Chauvin, 95302 Cergy-Pontoise Cedex, France

E-mail: timotheecrinbarat@gmail.com, shoulingyun11@gmail.com and jin.tan@cyu.fr

Received 3 April 2023; revised 28 January 2024
Accepted for publication 16 April 2024
Published 20 May 2024

Recommended by Dr Yao Yao

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

CrossMark

## Abstract

We derive a novel two-phase flow system in porous media as a relaxation limit of compressible multi-fluid systems. Considering a one-velocity Baer–Nunziato system with friction forces, we first justify its pressure-relaxation limit toward a Kapila model in a uniform manner with respect to the time-relaxation parameter associated with the friction forces. Then, we show that the diffusely rescaled solutions of the damped Kapila system converge to the solutions of the new two-phase porous media system as the time-relaxation parameter tends to zero. In addition, we also prove the convergence of the Baer–Nunziato system to the same two-phase porous media system as both relaxation parameters tend to zero. For each relaxation limit, we exhibit sharp rates of convergence in a critical regularity setting. Our proof is based on an elaborate low-frequency and high-frequency analysis via the Littlewood–Paley decomposition and includes three main ingredients: a refined spectral analysis

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

Original Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote>* Author to whom any correspondence should be addressed.</footnote>
of the linearized problem to determine the frequency threshold explicitly in terms of the time-relaxation parameter, the introduction of an effective flux in the low-frequency region to overcome the loss of parameters due to the *overdamping phenomenon*, and renormalized energy estimates in the high-frequency region to cancel higher-order nonlinear terms. To justify the convergence rates, we discover several *auxiliary unknowns* allowing us to recover crucial $O(\epsilon)$ bounds.

Keywords: multi-fluid system, pressure-relaxation limit,

overdamping phenomenon, critical regularity,

two-phase flow in porous media, Kapila system, Baer–Nunziato system

Mathematics Subject Classification numbers: 35Q35, 35B40, 76N10, 76T17

# 1. Introduction

## 1.1. Models and motivations

Multi-phase flows have been used to simulate a wide range of physical mixing phenomenon, from engineering to biological systems (see [1, 11, 32, 46] and the references therein). In the present paper, we investigate an inviscid compressible one-velocity Baer–Nunziato system with two different pressure laws in presence of drag forces, which was discussed in the recent work [10] of Bresch and Hillairet:

$$ egin{cases} 
\partial_t \alpha_+ + u \cdot \nabla \alpha_+ = \frac{\alpha_+ \alpha_-}{\epsilon} (P_+(\rho_+) - P_-(\rho_-)), \\
\partial_t (\alpha_\pm \rho_\pm) + \operatorname{div} (\alpha_\pm \rho_\pm u) = 0, \\
\partial_t (\rho u) + \operatorname{div} (\rho u \otimes u) + \nabla P + \frac{\rho u}{\tau} = 0, \quad x \in \mathbb{R}^d, \quad t > 0,
\end{cases} 
\quad (\text{BN}) $$

where the unknowns $\alpha_\pm = \alpha_\pm(t, x) \in [0, 1]$, $\rho_\pm = \rho_\pm(t, x) \ge 0$ and $u = u(t, x) \in \mathbb{R}^d$ stand for the volume fractions, the densities and the common velocity of two fluids (denoted by $+$ and $-$), respectively, which satisfy

$$ \alpha_+ + \alpha_- = 1, \quad \rho = \alpha_+ \rho_+ + \alpha_- \rho_-, \quad P = \alpha_+ P_+(\rho_+) + \alpha_- P_-(\rho_-). $$

The two positive constants $\epsilon$ and $\tau$ are (small) relaxation parameters associated to the pressure-relaxation and time-relaxation limits. Finally, the two pressures $P_+$ and $P_-$ take the gamma-law forms

$$ P_{\pm}(s) = A_{\pm}s^{\gamma_{\pm}} \quad \text{with constants} \quad A_{\pm} > 0, \quad 1 \le \gamma_{\pm} < \gamma_{\pm}. \quad (1.1) $$

The Baer–Nunziato terminology refers to the pressure-relaxation mechanism in the equations of volume fractions. Numerically, such relaxation procedure can simplify its resolution as it reduces the number of constraints by introducing new unknowns: two pressures instead of one. The readers can see [10] and references therein for more discussions on this pressure-relaxation process. Very recently, the one-dimensional version of system (BN) was rigorously derived by Bresch, Burtea and Lagoutiére in [6, 7].

There is an extensive literature on the mathematical analysis of multi-fluid systems. For example, in the one-velocity case, the global existence of weak solutions has been studied in [14, 37, 42, 45, 47, 52], and the global well-posedness and optimal time-decay rates of strong
solutions has been established in the framework of Sobolev spaces [30, 51, 53, 54] and critical Besov spaces [15, 31, 36], etc. We also refer to [9, 12, 13, 26, 35] on the study of multi-fluid systems in the two-velocity case. Complete reviews on multi-fluid systems are presented in [8, 48]. Concerning the study of relaxation problems associated to systems of conservation laws, it can be traced back to the work [16] by Chen, Levermore and Liu. Recently, Giovangigli and Yong in [28, 29] studied a relaxation problem arising in the dynamics of perfect gases out of thermodynamic equilibrium.

At the formal level, the solution $(\alpha_{\pm}^{\varepsilon, \tau}, \rho_{\pm}^{\varepsilon, \tau}, u^{\varepsilon, \tau})$ of the system (BN) tends, as $\varepsilon \to 0$, to some limit $(\alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, u^{\tau})$ that satisfies the so-called one-velocity Kapila system (see [34]):

$$ 
\left\{ 
\begin{array}{l}
\partial_t (\alpha_{\pm}^{\tau} \rho_{\pm}^{\tau}) + \operatorname{div} (\alpha_{\pm}^{\tau} \rho_{\pm}^{\tau} u^{\tau}) = 0, \\
\partial_t (\rho^{\tau} u^{\tau}) + \operatorname{div} (\rho^{\tau} u^{\tau} \otimes u^{\tau}) + \nabla P^{\tau} + \frac{\rho^{\tau} u^{\tau}}{\tau} = 0, \\
P^{\tau} = P_{+}^{\tau}(\rho_{+}) = P_{-}^{\tau}(\rho_{-}),
\end{array}
\right. 
\qquad (K) 
$$

with $\alpha_{\pm}^{\tau} + \alpha_{-}^{\tau} = 1$ and $\rho^{\tau} = \alpha_{+}^{\tau} \rho_{+}^{\tau} + \alpha_{-}^{\tau} \rho_{-}^{\tau}$. The system (K) can be rewritten as classical two-phase fluid models of drift-flux type, see [25, 32, 47] and the references therein. For existence of finite energy weak solutions to the system (K) with viscosities, refer to the recent works [14, 37, 42, 47].

Then, we further investigate the time-relaxation limit of the system (K) as $\tau \to 0$. Inspired by the works [17, 33, 50] concerning the relaxation problems for the compressible Euler system with damping, we introduce a large time-scale $O(1/\tau)$ and define the following charge of variables

$$ (\beta_{\pm}^{\tau}, \varrho_{\pm}^{\tau}, v^{\tau})(s,x) := \left( \alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, \frac{u^{\tau}}{\tau} \right) \left( \frac{s}{\tau}, x \right). \qquad (1.2) $$

Under the diffusive scaling (1.2), the system (K) becomes

$$ 
\left\{ 
\begin{array}{l}
\partial_s (\beta_{\pm}^{\tau} \varrho_{\pm}^{\tau}) + \operatorname{div} (\beta_{\pm}^{\tau} \varrho_{\pm}^{\tau} v^{\tau}) = 0, \\
\tau^2 \partial_s (\varrho^{\tau} v^{\tau}) + \tau^2 \operatorname{div} (\varrho^{\tau} v^{\tau} \otimes v^{\tau}) + \nabla \Pi^{\tau} + \varrho^{\tau} v^{\tau} = 0, \\
\beta_{+}^{\tau} + \beta_{-}^{\tau} = 1,
\end{array}
\right. 
\qquad (K_{\tau}) 
$$

with $\varrho^{\tau} = \beta_{+}^{\tau} \rho_{+}^{\tau} + \beta_{-}^{\tau} \varrho_{-}^{\tau}$ and $\Pi^{\tau} = P_{+}(\varrho_{+}^{\tau}) = P_{-}(\varrho_{-}^{\tau})$. As $\tau \to 0$, one then expects that $(\beta_{\pm}^{\tau}, \varrho_{\pm}^{\tau}, v^{\tau})$ converges to some limit $(\beta_{\pm}, \varrho_{\pm}, v)$ which is the solution of a new two-phase system

$$ 
\left\{ 
\begin{array}{l}
\partial_s (\beta_{\pm} \varrho_{\pm}) + \operatorname{div} (\beta_{\pm} \varrho_{\pm} v) = 0, \\
\nabla \Pi + \varrho v = 0, \\
\beta_{+} + \beta_{-} = 1,
\end{array}
\right. 
\qquad (1.3) 
$$

with $\varrho = \beta_{+} \varrho_{+} + \beta_{-} \varrho_{-}$ and $\Pi = P_{+}(\varrho_{+}) = P_{-}(\varrho_{-})$. Inserting Darcy's law (1.3)$_2$ into (1.3)$_1$, we derive the following two-phase system in porous media:

$$ 
\left\{ 
\begin{array}{l}
\partial_s \beta_{+} + v \cdot \nabla \beta_{+} = \frac{(\gamma_{+} - \gamma_{-}) \beta_{+} \beta_{-}}{\gamma_{+} \beta_{-} + \gamma_{-} \beta_{+}} \operatorname{div} \left( \frac{\nabla \Pi}{\beta_{+} \varrho_{+} + \beta_{-} \varrho_{-}} \right), \\
\partial_s \Pi + v \cdot \nabla \Pi = \frac{\gamma_{+} \gamma_{-} \Pi}{\gamma_{+} \beta_{-} + \gamma_{-} \beta_{+}} \operatorname{div} \left( \frac{\nabla \Pi}{\beta_{+} \varrho_{+} + \beta_{-} \varrho_{-}} \right), \\
\beta_{+} + \beta_{-} = 1, \\
\Pi = P_{+}(\varrho_{+}) = P_{-}(\varrho_{-}).
\end{array}
\right. 
\qquad (\mathrm{PM}) 
$$
The present paper is a follow-up to the paper [15] by Burtea, Crin-Barat and Tan where the authors justified the pressure-relaxation limit for the viscous version of the system (BN) to the system (K). In [15], the smallness condition on initial data employed to justify their global well-posedness result depends on $min\{\tau, \frac{1}{\tau}\}$ (due to the overdamping phenomenon that will be explained below) and therefore does not allow to further investigate the limit when $\tau \to 0$.

The main results of this article are the quantitative justification of the pressure-relaxation limit from the system (BN) to the system (K) as $\varepsilon \to 0$ uniformly in $\tau$ and the time-relaxation limit from the system (K$\tau$) to the system (PM) as $\tau \to 0$. Consequently, a new two-phase flow system in porous media (PM) is rigorously derived from (K) and (BN), which implies that Kapila and Baer-Nunziato systems considered in our paper can be viewed, for $\varepsilon$ and $\tau$ small enough, as hyperbolic approximations of (PM).

For both relaxation limits, we will focus on global-in-time strong solutions being small perturbations of constant equilibrium states. In other words, we consider solutions $(\alpha_{\pm}^{\varepsilon, \tau}, \rho_{\pm}^{\varepsilon, \tau}, u^{\varepsilon, \tau})$ to the system (BN) (resp. $(\alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, u^{\tau})$ to the system (K)) with positive densities and volume fractions which, as $|x| \to \infty$, tend to some thermodynamically stable equilibrium state $(\bar{\alpha}_{\pm}, \bar{\rho}_{\pm}, 0)$ fulfilling

$$0 < \bar{\alpha}_{\pm} < 1, \quad \bar{\alpha}_{+} + \bar{\alpha}_{-} = 1, \quad \bar{\rho}_{\pm} > 0, \quad P_{+}(\bar{\rho}_{+}) = P_{-}(\bar{\rho}_{-}). \quad (1.4)$$

For convenience, we also define the corresponding equilibrium state for the total density and the total pressure as

$$\bar{\rho} := \bar{\alpha}_{+} \bar{\rho}_{+} + \bar{\alpha}_{-} \bar{\rho}_{-}, \quad \bar{P} := P_{+}(\bar{\rho}_{+}) = P_{-}(\bar{\rho}_{-}). \quad (1.5)$$

To achieve our goals, we prove uniform in $\varepsilon$ and $\tau$ (such that $\varepsilon \le \tau$) a priori estimate for the system (BN) which improves the analysis performed in [15] that did not provide uniform-in-$\tau$ estimate. Such estimate allows us to justify a global well-posedness for a class of non-symmetric partially dissipative hyperbolic systems with rough coefficients in the context of overdamping phenomenon, which is not covered by the recent lecture of Danchin [24]. Indeed, our proof generalizes the techniques developed in [5, 18, 19, 21] which cannot be directly applied to the system (BN) due to the complex forms of the total pressure and the lack of symmetry.

On the other hand, it is natural to ask what happens for the system (BN) as $\tau$ tends to 0 first. To investigate this process, we introduce a diffusive scaling similar to (1.2) as follows

$$(\beta_{\pm}^{\varepsilon, \tau}, \varrho_{\pm}^{\varepsilon, \tau}, v^{\varepsilon, \tau})(s, x) := \left(\alpha_{\pm}^{\varepsilon, \tau}, \rho_{\pm}^{\varepsilon, \tau}, \frac{u^{\varepsilon, \tau}}{\tau}\right)\left(\frac{s}{\tau}, x\right). \quad (1.6)$$

Under such scaling (1.6), the system (BN) becomes

$$\left\{
\begin{array}{l}
\partial_s \beta_+^{\varepsilon, \tau} + v^{\varepsilon, \tau} \cdot \nabla \beta_+^{\varepsilon, \tau} = - \frac{\beta_+^{\varepsilon, \tau} \beta_-^{\varepsilon, \tau}}{\varepsilon \tau} (P_+ (\varrho_+^{\varepsilon, \tau}) - P_- (\varrho_-^{\varepsilon, \tau})), \\
\partial_s (\beta_+^{\varepsilon, \tau} \varrho_-^{\varepsilon, \tau}) + \text{div} (\beta_+^{\varepsilon, \tau} \varrho_-^{\varepsilon, \tau} v^{\varepsilon, \tau}) = 0, \\
\tau^2 \partial_s (\varrho_-^{\varepsilon, \tau} v^{\varepsilon, \tau}) + \tau^2 \text{div} (\varrho_-^{\varepsilon, \tau} v^{\varepsilon, \tau} \otimes v^{\varepsilon, \tau}) + \nabla \Pi^{\varepsilon, \tau} + \varrho_-^{\varepsilon, \tau} v^{\varepsilon, \tau} = 0,
\end{array}
\right.
\quad (\text{BN}_\tau)$$

with $\beta_+^{\varepsilon, \tau} + \beta_-^{\varepsilon, \tau} = 1$, $\varrho_-^{\varepsilon, \tau} = \beta_+^{\varepsilon, \tau} \varrho_-^{\varepsilon, \tau} + \beta_-^{\varepsilon, \tau} \varrho_-^{\varepsilon, \tau}$ and $\Pi^{\varepsilon, \tau} = \beta_+^{\varepsilon, \tau} P_+(\varrho_+^{\varepsilon, \tau}) + \beta_-^{\varepsilon, \tau} P_-(\varrho_-^{\varepsilon, \tau})$. The crucial observation is that the parameter $\tau$ now also appears under the pressure-relaxation term in the equation of the volume fractions. This reveals that as $\tau \to 0$, the two pressures in
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Relaxation limits diagram.</figcaption></figure>

the system (BNτ) should converge to a common pressure, and the solutions of the system (BN)
should converge to the solutions of the system (K) regardless of ε. Additionally, in the sequel
of the paper we are only able to justify the limit in the case ε ≪ τ which corresponds to the situ-
ation that the time-scale of the pressure-relaxation is small than the time-scale of the diffusive
relaxation. The condition ε ≪ τ appears in the spectral analysis of the system (see section 1.4)
and is essential for us to close the uniform *a priori* estimate in both low and high frequencies
(see sections 3.1 and 3.2 for the details). But in a formal way, the condition ε ≪ τ seems not
necessary in the limit process τ → 0, so the case ε > τ remains an interesting open problem
figure 1 summarizes the limit processes that we tackle in this article.

1.2. *Outline of the paper*

The rest of the paper is organized as follows. Our main results are stated in section 1.3. In
section 1.4, we first recall a reformulation of the system (BN) from [15] and present an explicit
spectral analysis for the associated linear system, then the difficulties and strategies of proof
are discussed. Section 2 is devoted to some notations and properties of Besov spaces and
Littlewood-Paley decomposition, and the regularity estimates for some linear problems are
stated. In section 3, we establish uniform *a priori* estimate for the linearized problem. Next,
in section 4, we prove the global existence and uniqueness results of solutions for (BN), (K)
and (PM), respectively. Section 5 is devoted to the justification of the relaxation limits with
explicit convergence rates.

**Notations.** We end this section by presenting a few notations. As usual, we denote by *C*
(and sometimes with subscripts) harmless positive constants that may change from line to
line, and *A* ≲ *B* (*A* ≳ *B*) means that both *A* ≪ *CB* (*A* ≳ *CB*), while *A* ~ *B* means that *A* ≲ *B* and
*A* ≳ *B*. For *X* a Banach space, *p* ∈ [1, ∞] and *T* > 0, the notation *L*^p(0, *T*; *X*) or *L*^p_T(X) designates
the set of measurable functions *f*: [0, *T*) → *X* with *t* ↦→ ||*f*(*t*)||_X in *L*^p(0, *T*), endowed with the
norm || · ||_T^p(X) := ||| · |||_X ||_L^p(0, *T*). We agree that *C*<sub>*b*</sub>([0, *T*]; *X*) denotes the set of continuous and
bounded (uniformly in *T*) functions from [0, *T*) to *X*. Sometimes, we use the notation *L*^p(X) to
designate the space *L*^p(R<sub>+</sub>; *X*) and || · ||_L^p(X) for the associated norm. We will keep the same
notations for multi-component functions, namely for *f*: [0, *T*) → *X*<sup>*m*</sup> with *m* ∈ N. ℱ and ℱ<sup>-1</sup>
stand for the Fourier transform and its inverse, and define the operator Λ<sup>σ</sup> := ℱ<sup>-1</sup>(|ξ|<sup>σ</sup> ℱ(·)).
Finally, let *o*(1) denote a generic constant which can be sufficiently small.
### 1.3. Main results

Our first theorem concerns the uniform, in both relaxation parameters $\varepsilon$ and $\tau$, global well-posedness of the system (BN) in a critical regularity framework.

**Theorem 1.1.** Let $d \ge 2$ and $0 < \varepsilon \le \tau \le 1$. Given the constants $\bar{\alpha}_{\pm}, \bar{\rho}_{\pm}$ verifying (1.4) and (1.5), assume that the initial data $(\alpha_{\pm,0}, \rho_{\pm,0}, u_0)$ satisfies $(\alpha_{\pm,0} - \bar{\alpha}_{\pm}, \rho_{\pm,0} - \bar{\rho}_{\pm}, u_0) \in \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}$. There exists a positive constant $c_0$ independent of $\tau$ and $\varepsilon$ such that if

$$ \| (\alpha_{\pm,0} - \bar{\alpha}_{\pm}, \rho_{\pm,0} - \bar{\rho}_{\pm}, u_0) \|_{\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}} \le c_0, \quad (1.7) $$

then the Cauchy problem of the system (BN) with the initial data $(\alpha_{\pm,0}, \rho_{\pm,0}, u_0)$ has a unique global solution $(\alpha_{\pm}^{\varepsilon, \tau}, \rho_{\pm}^{\varepsilon, \tau}, u^{\varepsilon, \tau})$ satisfying

$$ \left\{
\begin{array}{l}
(\alpha_{\pm}^{\varepsilon, \tau} - \bar{\alpha}_{\pm}, \rho_{\pm}^{\varepsilon, \tau} - \bar{\rho}_{\pm}, u^{\varepsilon, \tau}) \in C_b(\mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}), \\
P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau}) \in L^1(\mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}), \\
P^{\varepsilon, \tau} - \bar{P} \in L^1(\mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}) \cap L^2(\mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}), \\
u^{\varepsilon, \tau} \in L^1(\mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}) \cap L^2(\mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}).
\end{array}
\right.
\quad (1.8) $$

Moreover, the following uniform estimate holds:

$$ \begin{aligned}
& \| (\alpha_{\pm}^{\varepsilon, \tau} - \bar{\alpha}_{\pm}, \rho_{\pm}^{\varepsilon, \tau} - \bar{\rho}_{\pm}, u^{\varepsilon, \tau}) \|_{L^{\infty}(\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1})} + \| (\partial_t \alpha_{\pm}^{\varepsilon, \tau}, \partial_t \rho_{\pm}^{\varepsilon, \tau}, \partial_t u^{\varepsilon, \tau}) \|_{L^1(\dot{B}_{\frac{d}{2}}^{\frac{d}{2}})} \\
& \quad + \frac{1}{\varepsilon} \| P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau}) \|_{L^1(\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{\frac{d}{2}})} \\
& \quad + \frac{1}{\sqrt{\varepsilon}} \| P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau}) \|_{L^2(\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{\frac{d}{2}+1})} \\
& \quad + \tau \| P^{\varepsilon, \tau} - \bar{P} \|_{L^1(\dot{B}_{\frac{d}{2}+1})} + \sqrt{\tau} \| P^{\varepsilon, \tau} - \bar{P} \|_{L^2(\dot{B}_{\frac{d}{2}}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}+1})} \\
& \quad + \| u^{\varepsilon, \tau} \|_{L^1(\dot{B}_{\frac{d}{2}}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}+1})} + \frac{1}{\sqrt{\tau}} \| u^{\varepsilon, \tau} \|_{L^2(\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}+1})} \\
& \quad + \left\| \frac{\rho^{\varepsilon, \tau} u^{\varepsilon, \tau}}{\tau} + \nabla P^{\varepsilon, \tau} \right\|_{L^1(\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{\frac{d}{2}})} \le C \| (\alpha_{\pm,0} - \bar{\alpha}_{\pm}, \rho_{\pm,0} - \bar{\rho}_{\pm}, u_0) \|_{\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}}
\end{aligned}
\quad (1.9) $$

where $C > 0$ is a generic constant.

**Remark 1.1.** It should be emphasized that the regularity and decay-in-$\tau$ properties of the effective flux $\frac{\rho^{\varepsilon, \tau} u^{\varepsilon, \tau}}{\tau} + \nabla P^{\varepsilon, \tau}$ is better than the one verified by the solution $(\alpha_{\pm}^{\varepsilon, \tau}, \rho_{\pm}^{\varepsilon, \tau}, u^{\varepsilon, \tau})$. This is consistent with Darcy's law and plays key role in the justification of the time-relaxation limit.

By classical compactness arguments and the uniform estimate (1.9), we obtain the following global well-posedness theorems for (K) and (PM) in the critical regularity framework.

**Theorem 1.2.** Let $d \ge 2$ and $0 < \tau \le 1$. Given the constants $\bar{\alpha}_{\pm}, \bar{\rho}_{\pm}$ verifying (1.4) and (1.5), assume that the initial data $(\alpha_{\pm,0}, \rho_{\pm,0}, u_0)$ satisfies $(\alpha_{\pm,0} - \bar{\alpha}_{\pm}, \rho_{\pm,0} - \bar{\rho}_{\pm}, u_0) \in \dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}+1}$. There exists a positive constant $c_1$ independent of $\tau$ such that if
$$ \| (\alpha_{\pm,0} - \bar{\alpha}_{\pm}, \rho_{\pm,0} - \bar{\rho}_{\pm}, u_0) \|_{\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d}} \le c_1, \quad (1.10) $$

then the Cauchy problem of the system (K) with the initial data $(\alpha_{\pm,0}, \rho_{\pm,0}, u_0)$ admits a unique global solution $(\alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, u^{\tau})$ satisfying

$$ \left\{
\begin{array}{l}
(\alpha_{\pm}^{\tau} - \bar{\alpha}_{\pm}, \rho_{\pm}^{\tau} - \bar{\rho}_{\pm}, u^{\tau}) \in C_b \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d} \right), \\
P^{\tau} - \bar{P} \in L^1 \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}+1}^{d} \right) \cap L^2 \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d} \right), \\
u^{\tau} \in L^1 \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d} \right) \cap L^2 \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d} \right).
\end{array}
\right.
\quad (1.11) $$

Moreover, the following uniform estimate holds:

$$ \begin{aligned}
\| (\alpha_{\pm}^{\tau} - \bar{\alpha}_{\pm}, \rho_{\pm}^{\tau} - \bar{\rho}_{\pm}, u^{\tau}) \|_{L^{\infty}(\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d})} + \| (\partial_t \alpha_{\pm}^{\tau}, \partial_t \rho_{\pm}^{\tau}, \partial_t u^{\tau}) \|_{L^1(\dot{B}_{\frac{d}{2}}^{d})} \\
+ \tau \| P^{\tau} - \bar{P} \|_{L^1(\dot{B}_{\frac{d}{2}+1}^{d})} + \sqrt{\tau} \| P^{\tau} - \bar{P} \|_{L^2(\dot{B}_{\frac{d}{2}}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d})} \\
+ \| u^{\tau} \|_{L^1(\dot{B}_{\frac{d}{2}}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d})} + \frac{1}{\sqrt{\tau}} \| u^{\tau} \|_{L^2(\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d})} \\
+ \left\| \frac{\rho^{\varepsilon, \tau} u^{\tau}}{\tau} + \nabla P^{\tau} \right\|_{L^1(\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}}^{d})} \le C \| (\alpha_{\pm,0} - \bar{\alpha}_{\pm}, \rho_{\pm,0} - \bar{\rho}_{\pm}, u_0) \|_{\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d}},
\end{aligned} \quad (1.12) $$

where $C > 0$ is a generic constant.

**Theorem 1.3.** Let $d \ge 2$. Given the constants $\bar{\alpha}_{\pm}, \bar{\rho}_{\pm}$ verifying (1.4) and (1.5), assume that the initial data $(\beta_{\pm,0}, \varrho_{\pm,0})$ satisfies $(\beta_{\pm,0} - \bar{\alpha}_{\pm}, \varrho_{\pm,0} - \bar{\rho}_{\pm}) \in \dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d}$ and

$$ \| (\beta_{\pm,0} - \bar{\alpha}_{\pm}, \varrho_{\pm,0} - \bar{\rho}_{\pm}) \|_{\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d}} \le c_2, \quad (1.13) $$

for a positive constant $c_2$, then the Cauchy problem of the system (PM) with the initial data $(\beta_{\pm,0}, \varrho_{\pm,0})$ admits a unique global solution $(\beta_{\pm}, \varrho_{\pm})$, which satisfies

$$ \left\{
\begin{array}{l}
\beta_{\pm} - \bar{\alpha}_{\pm} \in C_b \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d} \right), \\
\varrho_{\pm} - \bar{\rho}_{\pm} \in C_b \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d} \right) \cap L^1 \left( \mathbb{R}_+; \dot{B}_{\frac{d}{2}}^{d} \right).
\end{array}
\right.
\quad (1.14) $$

Moreover, the following uniform estimate holds:

$$ \begin{aligned}
\| (\beta_{\pm} - \bar{\alpha}_{\pm}, \varrho_{\pm} - \bar{\rho}_{\pm}) \|_{L^{\infty}(\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d})} + \| (\partial_t \beta_{\pm}, \partial_t \varrho_{\pm}) \|_{L^1(\dot{B}_{\frac{d}{2}}^{d})} \\
+ \| \varrho_{\pm} - \bar{\rho}_{\pm} \|_{L^1(\dot{B}_{\frac{d}{2}+1}^{d} \cap \dot{B}_{\frac{d}{2}+3}^{d})} \\
\le C \| (\beta_{\pm,0} - \bar{\alpha}_{\pm}, \varrho_{\pm,0} - \bar{\rho}_{\pm}) \|_{\dot{B}_{\frac{d}{2}-1}^{d} \cap \dot{B}_{\frac{d}{2}+1}^{d}},
\end{aligned} \quad (1.15) $$

where $C > 0$ is a generic constant.

Next, we present the rigorous justifications of the pressure-relaxation limit for the system (BN) to the system (K) as $\varepsilon \to 0$ uniform with respect to $\tau$, and further the time-relaxation limit for the system (Kτ) to the system (PM) as $\tau \to 0$, with explicit convergence rates.

**Theorem 1.4.** Let $d \ge 2$ and $0 < \varepsilon \le \tau \le 1$. Given the constants $\bar{\alpha}_{\pm}, \bar{\rho}_{\pm}$ verifying (1.4) and (1.5), let $(\alpha_{\pm}^{\varepsilon, \tau}, \rho_{\pm}^{\varepsilon, \tau}, u^{\varepsilon, \tau})$, $(\alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, u^{\tau})$ and $(\beta_{\pm}, \varrho_{\pm})$ be the global solutions to the Cauchy
problems of (BN), (K) and (PM) obtained from Theorems 1.1–1.3 associated to their corresponding initial data $(\alpha_{\pm,0}^{\varepsilon,\tau}, \rho_{\pm,0}^{\varepsilon,\tau}, u_0^{\varepsilon,\tau})$, $(\alpha_{\pm,0}^{\tau}, \rho_{\pm,0}^{\tau}, u_0^{\tau})$ and $(\beta_{\pm,0}, \varrho_{\pm,0})$, respectively.

* Let the initial quantities $P_0^{\varepsilon,\tau} - P_0^{\tau}$ and $Y_0^{\varepsilon,\tau} - Y_0^{\tau}$ be denoted by (5.1) and (5.3), respectively. If $d \ge 3$ and

$$ \| (P_+(\rho_{+,0}^{\varepsilon,\tau}) - P_-(\rho_{-,0}^{\varepsilon,\tau}), Y_0^{\varepsilon,\tau} - Y_0^{\tau}, P_0^{\varepsilon,\tau} - P_0^{\tau}, u_0^{\varepsilon,\tau} - u_0^{\tau}) \|_{\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1}} \le \sqrt{\varepsilon\tau}, \quad (1.16) $$

then there exists a universal constant $C_1$ such that the following estimate holds:

$$ \begin{aligned} & \| (\alpha_{\pm}^{\varepsilon,\tau} - \alpha_{\pm}^{\tau}, \rho_{\pm}^{\varepsilon,\tau} - \rho_{\pm}^{\tau}, u_{\pm}^{\varepsilon,\tau} - u_{\pm}^{\tau}) \|_{L^{\infty}(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \\ & \quad + \sqrt{\tau} \| \rho_{\pm}^{\varepsilon,\tau} - \rho_{\pm}^{\tau} \|_{L^2(\dot{B}^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\tau}} \| u^{\varepsilon,\tau} - u^{\tau} \|_{L^2(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \\ & \quad + \| u^{\varepsilon,\tau} - u^{\tau} \|_{L^1(\dot{B}^{\frac{d}{2}-1})} \le C_1 \sqrt{\varepsilon\tau}. \end{aligned} \quad (1.17) $$

* Furthermore, define $(\beta_{\pm}^{\tau}, \varrho_{\pm}^{\tau}, v^{\tau})$ by the diffusive scaling (1.2) and $v$ by Darcy's law (1.3)$_2$. Let the initial quantity $Z_0^{\tau} - Z_0$ be denoted by (5.32). If

$$ \| Z_0^{\tau} - Z_0 \|_{\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}}} + \| \varrho_{\pm,0}^{\tau} - \varrho_{\pm,0} \|_{\dot{B}^{\frac{d}{2}-1}} \le \tau, \quad (1.18) $$

then there exists a universal constant $C_2$ such that the following estimate holds:

$$ \begin{aligned} & \| (\beta_{\pm}^{\tau} - \beta_{\pm}, \varrho_{\pm}^{\tau} - \varrho_{\pm}) \|_{L^{\infty}(\dot{B}^{\frac{d}{2}-1})} \\ & \quad + \| \varrho_{\pm}^{\tau} - \varrho_{\pm} \|_{L^1(\dot{B}^{\frac{d}{2}+1})} + \| v^{\tau} - v \|_{L^1(\dot{B}^{\frac{d}{2}})} \le C_2 \tau. \end{aligned} \quad (1.19) $$

Finally, Theorem 1.4 implies the relaxation limit for the system (BN$\tau$) to the system (PM) as both $\varepsilon, \tau \to 0$.

**Corollary 1.1.** Let $d \ge 3$, $0 < \varepsilon \le \tau \le 1$, and $(\beta_{\pm}^{\varepsilon,\tau}, \varrho_{\pm}^{\varepsilon,\tau}, v^{\varepsilon,\tau})$ be defined by (1.6). Then under the assumptions of Theorem 1.4, there is a generic constant $C_3$ such that

$$ \| (\beta_{\pm}^{\varepsilon,\tau} - \beta_{\pm}, \varrho_{\pm}^{\varepsilon,\tau} - \varrho_{\pm}) \|_{L^{\infty}(\dot{B}^{\frac{d}{2}-1})} \le C_3 (\sqrt{\varepsilon\tau} + \tau). $$

## 1.4. Difficulties and strategies

The first difficulty concerning the study of the system (BN) is its lack of dissipativity and symmetrizability. Indeed, the linearization of (BN) admits the eigenvalue 0 and therefore does not satisfy the well-known 'Shizuta-Kawashima' stability condition for partially dissipative hyperbolic systems (see [44]). Additionally, the system (BN) cannot be written in a conservative form and the entropy naturally associated to (BN) is not positive definite, therefore the notion of entropic variables does not make sense in this case. Therefore, the first crucial step in our analysis is to partially symmetrize the system (BN), by hands. We refer to [13, 27] for
the treatment of non-conservative systems in similar contexts. In our setting, as explained in [15], we define the new unknowns

$$
\left\{
\begin{aligned}
y^{\epsilon, \tau} &:= \frac{\alpha_{+}^{\epsilon, \tau} \rho_{+}^{\epsilon, \tau}}{\alpha_{+}^{\epsilon, \tau} \rho_{+}^{\epsilon, \tau} + \alpha_{-}^{\epsilon, \tau} \rho_{-}^{\epsilon, \tau}} - \frac{\bar{\alpha}_{+} \bar{\rho}_{+}}{\bar{\alpha}_{+} \bar{\rho}_{+} + \bar{\alpha}_{-} \bar{\rho}_{-}}, \\
w^{\epsilon, \tau} &:= \frac{\alpha_{+}^{\epsilon, \tau} \alpha_{-}^{\epsilon, \tau}}{\gamma_{+} \alpha_{-}^{\epsilon, \tau} + \gamma_{-} \alpha_{+}^{\epsilon, \tau}} (P_{+}(\rho_{+}^{\epsilon, \tau}) - P_{-}(\rho_{-}^{\epsilon, \tau})), \\
r^{\epsilon, \tau} &:= P^{\epsilon, \tau} - \bar{P} - (\gamma_{+} - \gamma_{-}) w^{\epsilon, \tau},
\end{aligned}
\right.
\qquad (1.20)
$$

and the corresponding initial data

$$
\left\{
\begin{aligned}
y_0 &:= \frac{\alpha_{+,0} \rho_{+,0}}{\alpha_{+,0} \rho_{+,0} + \alpha_{-,0} \rho_{-,0}} - \frac{\bar{\alpha}_{+} \bar{\rho}_{+}}{\bar{\alpha}_{+} \bar{\rho}_{+} + \bar{\alpha}_{-} \bar{\rho}_{-}}, \\
w_0 &:= \frac{\alpha_{+,0} \alpha_{-,0}}{\gamma_{+} \alpha_{-,0} + \gamma_{-} \alpha_{+,0}} (P_{+}(\rho_{+,0}) - P_{-}(\rho_{-,0})), \\
r_0 &:= \alpha_{+,0} P_{+}(\rho_{+,0}) + \alpha_{-,0} P_{-}(\rho_{-,0}) - \bar{P} - (\gamma_{+} - \gamma_{-}) w_0,
\end{aligned}
\right.
\qquad (1.21)
$$

so that the Cauchy problem of the system (BN) subject to the initial data $(\alpha_{\pm,0}, \rho_{\pm,0}, u_0)$ is
reformulated as

$$
\left\{
\begin{aligned}
& \partial_t y^{\epsilon, \tau} + u^{\epsilon, \tau} \cdot \nabla y^{\epsilon, \tau} = 0, \\
& \partial_t w^{\epsilon, \tau} + u^{\epsilon, \tau} \cdot \nabla w^{\epsilon, \tau} + (\bar{F}_1 + G_1^{\epsilon, \tau}) \operatorname{div} u^{\epsilon, \tau} + (\bar{F}_2 + G_2^{\epsilon, \tau}) \frac{w^{\epsilon, \tau}}{\epsilon} = 0, \\
& \partial_t r^{\epsilon, \tau} + u^{\epsilon, \tau} \cdot \nabla r^{\epsilon, \tau} + (\bar{F}_3 + G_3^{\epsilon, \tau}) \operatorname{div} u^{\epsilon, \tau} = F_4^{\epsilon, \tau} \frac{(w^{\epsilon, \tau})^2}{\epsilon}, \\
& \partial_t u^{\epsilon, \tau} + u^{\epsilon, \tau} \cdot \nabla u^{\epsilon, \tau} + \frac{u^{\epsilon, \tau}}{\tau} + (\bar{F}_0 + G_0^{\epsilon, \tau}) \nabla r^{\epsilon, \tau} + (\gamma_+ - \gamma_-) (\bar{F}_0 + G_0^{\epsilon, \tau}) \nabla w^{\epsilon, \tau} = 0, \\
& (y^{\epsilon, \tau}, w^{\epsilon, \tau}, r^{\epsilon, \tau}, u^{\epsilon, \tau}) (0, x) = (y_0, w_0, r_0, u_0) (x),
\end{aligned}
\right.
\tag{1.22}
$$

where $F_i^{\epsilon, \tau} = F_i^{\epsilon, \tau}(y, w, r) (i = 0, 1, 2, 3, 4)$ are the nonlinear terms

$$
\left\{
\begin{aligned}
F_0^{\epsilon, \tau} &:= \frac{1}{\alpha_+^{\epsilon, \tau} \rho_+^{\epsilon, \tau} + \alpha_-^{\epsilon, \tau} \rho_-^{\epsilon, \tau}}, \\
F_1^{\epsilon, \tau} &:= \frac{(\gamma_+ - \gamma_-) \alpha_+^{\epsilon, \tau} \alpha_-^{\epsilon, \tau}}{\gamma_+ \alpha_-^{\epsilon, \tau} + \gamma_- \alpha_+^{\epsilon, \tau}} (\bar{P} + r^{\epsilon, \tau}) + \frac{\gamma_+^2 \alpha_-^{\epsilon, \tau} + \gamma_-^2 \alpha_+^{\epsilon, \tau}}{\gamma_+ \alpha_-^{\epsilon, \tau} + \gamma_- \alpha_+^{\epsilon, \tau}} w^{\epsilon, \tau}, \\
F_2^{\epsilon, \tau} &:= (\gamma_+ \alpha_-^{\epsilon, \tau} + \gamma_- \alpha_+^{\epsilon, \tau}) (\bar{P} + r^{\epsilon, \tau}) - \frac{(\gamma_+ - \gamma_+^2) (\alpha_-^{\epsilon, \tau})^2 - (\gamma_- - \gamma_-^2) (\alpha_+^{\epsilon, \tau})^2}{\alpha_+^{\epsilon, \tau} \alpha_-^{\epsilon, \tau}} w^{\epsilon, \tau}, \\
F_3^{\epsilon, \tau} &:= \frac{\gamma_+ \gamma_-}{\gamma_+ \alpha_-^{\epsilon, \tau} + \gamma_- \alpha_+^{\epsilon, \tau}} P^{\epsilon, \tau}, \\
F_4^{\epsilon, \tau} &:= \frac{\gamma_+ \gamma_-}{\alpha_+^{\epsilon, \tau} \alpha_-^{\epsilon, \tau}} (1 - \gamma_+ \alpha_-^{\epsilon, \tau} - \gamma_- \alpha_+^{\epsilon, \tau}),
\end{aligned}
\right.
\tag{1.23}
$$
$\bar{F}_i$ ($i = 0, 1, 2, 3$) are the constants

$$
\left\{
\begin{aligned}
\bar{F}_0 &:= \frac{1}{\bar{\alpha}_+ \bar{\rho}_+ + \bar{\alpha}_- \bar{\rho}_-} > 0, \\
\bar{F}_1 &:= \frac{(\gamma_+ - \gamma_-) \bar{\alpha}_+ \bar{\alpha}_-}{\gamma_+ \bar{\alpha}_- + \gamma_- \bar{\alpha}_+} \bar{P} > 0, \\
\bar{F}_2 &:= (\gamma_+ \bar{\alpha}_- + \gamma_- \bar{\alpha}_+) \bar{P} > 0, \\
\bar{F}_3 &:= \frac{\gamma_+ \gamma_-}{\gamma_+ \bar{\alpha}_- + \gamma_- \bar{\alpha}_+} \bar{P} > 0,
\end{aligned}
\right.
\qquad (1.24)
$$

and $G_i^{\epsilon, \tau} = G_i^{\epsilon, \tau}(y, w, r)$ ($i = 0, 1, 2, 3$) are the coefficients

$$ G_i^{\epsilon, \tau} := F_i^{\epsilon, \tau} - \bar{F}_i. \qquad (1.25) $$

In this formulation, the equation (1.22)$_1$ is purely transport and the linear part of the subsystem (1.22)$_2$–(1.22)$_4$ is partially dissipative and satisfies the ‘Shizuta–Kawashima’ stability condition. Thus, we will estimate the undamped unknown $y^{\epsilon, \tau}$ and the dissipative components ($w^{\epsilon, \tau}$, $r^{\epsilon, \tau}$, $u^{\epsilon, \tau}$) separately. We emphasize here that due to the double parameters $\epsilon, \tau$ and the lack of time-integrability of $G_i^{\epsilon, \tau}$, the dissipative structures of the subsystem (1.22)$_2$–(1.22)$_4$ does not fit into the general theorems that can be found in [18, 19, 24, 44, 49, 50], and a new analysis is needed to be developed to obtain the uniform estimates with respect to the two relaxation parameters $\epsilon, \tau$.

In order to understand the behaviors of the solution to (1.22) with respect to $\epsilon, \tau$, we perform a spectral analysis of the linear system for (1.22). For simplicity we set $\bar{F}_i = 1$ ($i = 0, 2, 3$) and $\bar{F}_1 = \gamma_+ - \gamma_-$. In terms of Hodge decomposition, we denote the compressible part $m = \Lambda^{-1} \text{div } u$ and the incompressible part $\Omega = \Lambda^{-1} \nabla \times u$ and rewrite the linear system of (1.22) as

$$ \partial_t \begin{pmatrix} w \\ r \\ m \end{pmatrix} = \mathbb{A} \begin{pmatrix} w \\ r \\ m \end{pmatrix}, \quad \mathbb{A} := \begin{pmatrix} -\frac{1}{\epsilon} & 0 & -(\gamma_+ - \gamma_-)\Lambda \\ 0 & 0 & -\Lambda \\ (\gamma_+ - \gamma_-)\Lambda & \Lambda & -\frac{1}{\tau} \end{pmatrix}, \quad \partial_t \Omega + \frac{1}{\tau}\Omega = 0. $$

The eigenvalues of the matrix $\hat{\mathbb{A}}(\xi)$ satisfy

$$ |\hat{\mathbb{A}}(\xi) - \lambda \mathbb{I}_{3 \times 3}| = \lambda^3 + \left(\frac{1}{\tau} + \frac{1}{\epsilon}\right) \lambda^2 + \left[\frac{1}{\epsilon\tau} + (|\gamma_+ - \gamma_-|^2 + 1) |\xi|^2\right] \lambda + \frac{1}{\epsilon}|\xi|^2 = 0. $$

Under the condition $0 < \epsilon \ll \tau$, the behaviors of $\lambda_i$ ($i = 1, 2, 3$) can be analyzed as follows:

* In the low-frequency region $|\xi| \ll \frac{1}{\tau}$, by Taylor's expansion near $|\tau\xi| \ll 1$ as in [41], all the eigenvalues are real, and we have $\lambda_1 = -\frac{1}{\epsilon} + \frac{1}{\tau}O(|\tau\xi|^2)$, $\lambda_2 = -\tau|\xi|^2 + \frac{1}{\tau}O(|\tau\xi|^3)$ and $\lambda_3 = -\frac{1}{\tau} + \frac{1}{\epsilon}O(|\tau\xi|^2)$.

* In the medium-frequency region $\frac{1}{\tau} \ll |\xi| \ll \frac{1}{\epsilon}$, according to Cardano's formula, $\lambda_1$ is real and $\lambda_i$ ($i = 2, 3$) are conjugated complex, and Re $\lambda_i \lesssim -\frac{1}{\tau}$ holds for all $i = 1, 2, 3$.

* In the high-frequency region $|\xi| \gg \frac{1}{\epsilon}$, by Taylor's expansion near $|\epsilon\xi|^{-1} \ll 1$, the real eigenvalue $\lambda_1$ and the conjugated complex eigenvalues $\lambda_i$ ($i = 2, 3$) satisfy $\lambda_1 = -\frac{1}{|\gamma_+ - \gamma_-|^2 + 1} + \frac{1}{\epsilon}O(\frac{1}{|\epsilon\xi|^2})$ and $\lambda_{2,3} = -\frac{1}{2\tau} - \frac{|\gamma_+ - \gamma_-|^2}{|\gamma_+ - \gamma_-|^2 + 1} + \frac{1}{2\epsilon} \pm \sqrt{|\gamma_+ - \gamma_-|^2 + 1}|\xi|i + (\frac{1}{\tau} + \frac{|\gamma_+ - \gamma_-|^2}{\epsilon})O(\frac{1}{|\epsilon\xi|})$.
<figure><img src="image_5.png" /><figcaption>Figure 2. A graph of overdamping phenomenon for the system (1.26).</figcaption></figure>

The above spectral analysis suggests us to separate the whole frequencies into two parts $|\xi| \lesssim \frac{1}{\tau}$ and $|\xi| \gtrsim \frac{1}{\tau}$ so as to capture the qualitative properties of solutions for the system (1.22). Indeed, the time-decay rates (determined by $\lambda_2$) achieve the fastest rate in the low-frequency region $|\xi| \lesssim \frac{1}{\tau}$. Moreover this region recover the whole frequency-space when $\tau \to 0$, as expected from the well-known overdamping phenomenon which will be mentioned below. To this end, the threshold $J_\tau$ between these two regions is used in the definition of the hybrid Besov spaces in next section.

It should be noted that $\lambda_2$ and $\lambda_3$ exhibit similar behaviors to the eigenvalues of the compressible Euler equations with damping. Indeed, to study the system (1.22), one considers the following simplified system of damped Euler type with rough coefficients:

$$ \left\{ \begin{array}{l} \partial_t r^\tau + (1 + G_3^\tau) \operatorname{div} u^\tau = 0, \\ \partial_t u^\tau + (1 + G_0^\tau) \nabla r^\tau + \frac{u^\tau}{\tau} = 0. \end{array} \right. \qquad (1.26) $$

The well-known spectral analysis for the linear Euler part of the system (1.26) implies that the frequency space shall be separated into the low-frequency region $|\xi| \lesssim \frac{1}{\tau}$ and the high-frequency region $|\xi| \gtrsim \frac{1}{\tau}$ to recover the uniform estimates and optimal regularity of solutions. Formally, this implies that as $\tau \to 0$, the low-frequency region covers the whole frequency space and is therefore be dominant at the limit. We observe here the classical overdamping phenomenon: as the friction coefficient $\frac{1}{\tau}$ gets larger, the decay rates of $r^\tau$ do not necessarily increase and on the contrary follow $\min\{\tau, \frac{1}{\tau}\}$, cf figure 2. For more discussion on the overdamping phenomenon, see Zuazua's sildes [55].

Recently, in [18, 19], the issue concerning the relaxation limit from compressible Euler system with damping toward the porous media equation has been rigorously justified in critical space $\dot{B}_{\frac{d}{2}}^d \cap \dot{B}_{\frac{d}{2}+1}^{d+1}$. The readers also can refer to the work [20] about the relaxation limit for a hyperbolic-parabolic chemotaxis system to a parabolic-elliptic Keller-Segel model. The regularity index $\frac{d}{2} + 1$ is called critical for initial data of general hyperbolic systems since $\dot{B}_{\frac{d}{2}+1}$ is embedded in the set of globally Lipschitz functions. Indeed, It has been observed by many authors that controlling the Lipschitz regularity of solutions for general hyperbolic systems can prevent blow-up in finite time, see e.g. [40, 49]. We also refer to [38, 39] about the ill-posedness results for hyperbolic systems in $H^s$ with $s < \frac{d}{2} + 1$.

Nevertheless, the methods developed in [18–20] are not applicable in the current situation to derive estimates which are uniform with respect to the relaxation parameter $\tau$. This is
mainly due to the complex form of the total pressure $P^{\epsilon, \tau}$ in the velocity equation (BN)$_3$ and the fact that one can not expect any time integrability property on $\mathbb{R}_+$ for the purely transported unknown $y^{\epsilon, \tau}$, which generally leads to a lack of time integrability on $\mathbb{R}_+$ for $(\alpha_{\pm}^{\epsilon, \tau} - \bar{\alpha}_{\pm}, \rho_{\pm}^{\epsilon, \tau} - \bar{\rho}_{\pm})$ (see Remarks 3.2 and 3.3), and thus for $G_3^{\tau}, G_0^{\tau}$ in the system (1.26). In addition, we can not find a rescaling to reduce the proof to the case $\tau = 1$ and then recover the corresponding uniform estimates with respect to $\tau$ thanks to the homogeneity of the Besov norms as in [18, 19]. To overcome these new difficulties, we will keep track of the dependence of $\epsilon, \tau$ and perform elaborate energy estimates with mixed $L^1$-time and $L^2$-time type dissipation. More precisely, in the low-frequency region, we introduce a purely damped mode (effective flux)

$$u^{\epsilon, \tau} + \frac{\tau}{\rho^{\epsilon, \tau}} \nabla P^{\epsilon, \tau}$$

corresponding to Darcy's law (1.3)$_2$ in the low-frequency setting to partially diagonalize the system and capture maximal dissipative structures. In addition, we derive some uniform estimates at a lower regularity level compared to [18–20] (see Lemma 3.1). In the high-frequency setting, due to the lack of symmetry, we need to cancel higher-order terms so as not to lose derivatives. For that, the construction of a Lyapunov functional in the spirit of Beauchard and Zuazua as in [3] with additional nonlinear weights allows us to capture the $L^1$-in-time dissipation properties in high frequencies (cf Lemma 3.2). Moreover, we also establish the uniform $L^2$-in-time estimates at $B_2^{\frac{d}{d+1}-\text{regularity}}$ level to recover the necessary bounds of parameters (refer to Lemma 3.3). Applying these ideas, we obtain uniform estimates in terms of the parameters $\epsilon, \tau$ satisfying $0 < \epsilon \le \tau < 1$ for the linearized problem (see proposition 3.1), which is crucial for our later nonlinear analysis.

Let us finally sketch the proof of the justifications of the strong relaxation limits. In fact, to obtain convergence rates, we will not estimate the differences of solutions between systems directly. The reason shares similarities with the proof of global uniform well-posedness. Roughly speaking, since both pressure-relaxation limit and time-relaxation limit are singular limits, there are singular terms which are only uniformly bounded but not necessarily vanishing in the equations satisfied by the difference of solutions. To overcome these difficulties, we discover some auxiliary unknowns associated with the difference systems, which reveal better structures (cancellations), and then perform error estimates on them for each relaxation limit. More details are presented in sections 5.1 and 5.2.

## 2. Functional framework and tools

In this section, we recall the notations of the Littlewood–Paley decomposition and Besov spaces. The reader can refer to [2, chapter 2] for a complete overview. Choose a smooth radial non-increasing function $\chi(\xi)$ with compact supported in $B(0, \frac{4}{3})$ and $\chi(\xi) = 1$ in $B(0, \frac{3}{4})$ such that

$$\varphi(\xi) := \chi(\xi/2) - \chi(\xi), \quad \sum_{j \in \mathbb{Z}} \varphi(2^{-j} \cdot) = 1, \quad \text{Supp } \varphi \subset \left\{ \xi \in \mathbb{R}^d \mid \frac{3}{4} \le |\xi| \le \frac{8}{3} \right\}.$$

For any $j \in \mathbb{Z}$, the homogeneous dyadic blocks $\hat{\Delta}_j$ and the low-frequency cut-off operator $\hat{S}_j$ are defined by

$$\hat{\Delta}_j u := \mathcal{F}^{-1}(\varphi(2^{-j} \cdot) \mathcal{F}u), \quad \hat{S}_j u := \mathcal{F}^{-1}(\chi(2^{-j} \cdot) \mathcal{F}u).$$
From now on, we use the shorthand notation

$$ \dot{\Delta}_j u = u_j. $$

Let $S'_h$ be the set of tempered distributions on $\mathbb{R}^d$ such that every $u \in S'_h$ satisfies $u \in S'$ and $\lim_{j \to -\infty} \|\dot{S}_j u\|_{L^\infty} = 0$. Then it follows that

$$ u = \sum_{j \in \mathbb{Z}} u_j \quad \text{in } S', \qquad \dot{S}_j u = \sum_{j' \le j-1} u_{j'}, \quad \forall u \in S'_h, $$

With the help of these dyadic blocks, the homogeneous Besov space $\dot{B}^s$ for $s \in \mathbb{R}$ is defined by

$$ \dot{B}^s := \left\{ u \in S'_h \mid \|u\|_{\dot{B}^s} := \sum_{j \in \mathbb{Z}} 2^{js} \|u_j\|_{L^2} < \infty \right\}. $$

We denote the Chemin-Lerner type space $\tilde{L}^\rho(0, T; \dot{B}^s)$ for $s \in \mathbb{R}$ and $T > 0$:

$$ \tilde{L}^\rho(0, T; \dot{B}^s) := \left\{ u \in L^\rho(0, T; S'_h) \mid \|u\|_{\tilde{L}^\rho_T(\dot{B}^s)} := \sum_{j \in \mathbb{Z}} 2^{js} \|u_j\|_{L^\rho_T(L^2)} < \infty \right\}. $$

By the Minkowski inequality, it holds that

$$ \|u\|_{L_T^\rho(\dot{B}^s)} \le \|u\|_{\tilde{L}_T^\rho(\dot{B}^s)} \quad \rho > 1, \qquad \|u\|_{L_T^1(\dot{B}^s)} = \|u\|_{\tilde{L}_T^1(\dot{B}^s)}, $$

where $\| \cdot \|_{L_T^\rho(\dot{B}^s)}$ is the usual Lebesgue-Besov norm.

In order to perform our analysis on the low and high frequencies regions, we set the threshold

$$ J_\tau := -[\log_2 \tau] + k, \qquad (2.1) $$

for suitable negative integer $k$ (to be determined in (3.13)-(3.16) in the proof of Lemma 3.1). Denote the following notations for $p \in [1, \infty]$ and $s \in \mathbb{R}$:

$$ \|u\|_{\dot{B}^s}^{\ell} := \sum_{j \le J_{\tau}} 2^{js} \|u_j\|_{L^2}, \qquad \|u\|_{\dot{B}^s}^{h} := \sum_{j \ge J_{\tau}-1} 2^{js} \|u_j\|_{L^\rho_T(L^2)}, $$

$$ \|u\|_{L_T^\rho(\dot{B}^s)}^{\ell} := \sum_{j \le J_\tau} 2^{js} \|u_j\|_{L^2}, \qquad \|u\|_{\tilde{L}_T^\rho(\dot{B}^s)}^{h} := \sum_{j \ge J_\tau-1} 2^{js} \|u_j\|_{\tilde{L}_T^\rho(L^2)}. $$

For any $u \in S'_h$, we also define the low-frequency part $u^\ell$ and the high-frequency part $u^h$ by

$$ u^{\ell} := \sum_{j \le J_{\tau}-1} u_j, \qquad u^{h} := u - u^{\ell} = \sum_{j \ge J_{\tau}} u_j. $$

It is easy to check for any $s' > 0$ that

$$ \left\{ \begin{array}{l} \|u^\ell\|_{\dot{B}^s} \le \|u\|_{\dot{B}^s}^\ell \le 2^{J_\tau s'} \|u\|_{\dot{B}^{s-s'}}^\ell \le 2^{s'} (2^k \tau^{-1})^{s'} \|u\|_{\dot{B}^{s-s'}}^\ell, \\ \|u^h\|_{\dot{B}^s} \le \|u\|_{\dot{B}^s}^h \le 2^{-(J_\tau-1)s'} \|u\|_{\dot{B}^{s-s'}}^h \le 2^{s'} (2^{-k}\tau)^{s'} \|u\|_{\dot{B}^{s-s'}}^h. \end{array} \right. \qquad (2.2) $$

Next, we state some properties of Besov spaces and related estimates which will be repeatedly used in the rest of paper. The reader can refer to [2, chapters 2 and 3] for more
details. Below, all the properties of Besov norms can be easily extended to the Chemin-Lerner norms.

The first lemma pertains to the so-called Bernstein's inequalities.

**Lemma 2.1.** Let $0 < r < R$, $1 \le p \le q \le \infty$ and $k \in \mathbb{N}$. For any function $u \in L^p$ and $\lambda > 0$, it holds that

$$ 
\left\{ 
\begin{array}{l}
\text{Supp } \mathcal{F}(u) \subset \{\xi \in \mathbb{R}^d \mid |\xi| \le \lambda R\} \Rightarrow \|D^k u\|_{L^q} \lesssim \lambda^{k+d(\frac{1}{p}-\frac{1}{q})} \|u\|_{L^p}, \\
\text{Supp } \mathcal{F}(u) \subset \{\xi \in \mathbb{R}^d \mid \lambda r \le |\xi| \le \lambda R\} \Rightarrow \|D^k u\|_{L^p} \sim \lambda^k \|u\|_{L^p}.
\end{array}
\right. 
$$

Due to the Bernstein inequalities, the Besov spaces have many useful properties:

**Lemma 2.2.** The following properties hold:

* For any $s \in \mathbb{R}$ and $q \ge 2$, we have the following continuous embeddings:

$$ \dot{B}^s \hookrightarrow \dot{H}^s, \quad \dot{B}^{\frac{d}{2}-\frac{d}{q}} \hookrightarrow L^q. $$

* $\dot{B}^d$ is continuously embedded in the set of continuous functions decaying to 0 at infinity.

* For any $\sigma \in \mathbb{R}$, the operator $\Lambda^\sigma$ is an isomorphism from $\dot{B}^s$ to $\dot{B}^{s-\sigma}$.

* Let $s_1 \in \mathbb{R}$ and $s_2 \le \frac{d}{2}$. Then the space $\dot{B}^{s_1} \cap \dot{B}^{s_2}$ is a Banach space and satisfies weak compact and Fatou properties: If $u_k$ is a uniformly bounded sequence of $\dot{B}^{s_1} \cap \dot{B}^{s_2}$, then an element $u$ of $\dot{B}^{s_1} \cap \dot{B}^{s_2}$ and a subsequence $u_{n_k}$ exist such that

$$ \lim_{k \to \infty} u_{n_k} = u \quad \text{in } S' \quad \text{and} \quad \|u\|_{\dot{B}^{s_1} \cap \dot{B}^{s_2}} \lesssim \liminf_{n_k \to \infty} \|u_{n_k}\|_{\dot{B}^{s_1} \cap \dot{B}^{s_2}}. $$

The following Morse-type product estimates in Besov spaces play a fundamental role in our analysis of nonlinear terms.

**Lemma 2.3.** The following statements hold:

* Let $s > 0$. Then $\dot{B}^s \cap L^\infty$ is an algebra and

$$ \|uv\|_{\dot{B}^s} \lesssim \|u\|_{L^\infty} \|v\|_{\dot{B}^s} + \|v\|_{L^\infty} \|u\|_{\dot{B}^s}. \quad (2.3) $$

* Let $s_1, s_2$ satisfy $s_1, s_2 \le \frac{d}{2}$ and $s_1 + s_2 > 0$. Then there holds

$$ \|uv\|_{\dot{B}^{s_1+s_2-\frac{d}{2}}} \lesssim \|u\|_{\dot{B}^{s_1}} \|v\|_{\dot{B}^{s_2}}. \quad (2.4) $$

The following commutator estimates will be used to control some nonlinearities in high frequencies.

**Lemma 2.4.** Let $p \in [1, \infty]$ and $-\frac{d}{2} - 1 \le s \le \frac{d}{2} + 1$. Then it holds that

$$ \sum_{j \in \mathbb{Z}} 2^{js} \| [v, \dot{\Delta}_j] \partial_i u \|_{L^2} \lesssim \|\nabla v\|_{\dot{B}^{\frac{d}{2}}} \|u\|_{\dot{B}^s}, \quad i = 1, \dots, d, \quad (2.5) $$

for the commutator $[A, B] := AB - BA$.

We prove the following lemma about the continuity for composition of multi-component functions. It should be noted that (2.7) will be used to deal with the two-dimensional case in $\dot{B}^0$.
**Lemma 2.5.** Let $m \in \mathbb{N}$, $s > 0$, and $G \in C^\infty(\mathbb{R}^m)$ satisfy $G(0, \dots, 0) = 0$. Then for any $f_i \in \dot{B}^s \cap L^\infty$ ($i = 1, \dots, m$), there exists a constant $C_f > 0$ depending on $\sum_{i=1}^m \|f_i\|_{L^\infty}$, $F$, $s$, $m$ and $d$ such that

$$ \|G(f_1, \dots, f_m)\|_{\dot{B}^s} \le C_f \sum_{i=1}^{m} \|f_i\|_{\dot{B}^s}. \qquad (2.6) $$

In the case $s > -\frac{d}{2}$ and $f_i \in \dot{B}^s \cap \dot{B}^{\frac{d}{2}}$, it holds that

$$ \|G(f_1, \dots, f_m)\|_{\dot{B}^s} \le C_f \left( 1 + \sum_{i=1}^{m} \|f_i\|_{\dot{B}^{\frac{d}{2}}} \right) \sum_{i=1}^{m} \|f_i\|_{\dot{B}^s}. \qquad (2.7) $$

Furthermore, for any $f_1^1, f_2^1 \in \dot{B}^s \cap \dot{B}^{\frac{d}{2}}$, we have

$$ \|G(f_1^1, \dots, f_m^1) - G(f_1^1, \dots, f_m^1)\|_{\dot{B}^s} \le C_f^* \left( 1 + \sum_{i=1}^{m} \|f_i\|_{\dot{B}^s \cap \dot{B}^{\frac{d}{2}}} \right) \sum_{i=1}^{m} \|f_i^1 - f_i^2\|_{\dot{B}^s \cap \dot{B}^{\frac{d}{2}}}. \qquad (2.8) $$

Here $C_f^* > 0$ depends on $\sum_{i=1}^m \|(f_i^1, f_i^2)\|_{L^\infty}$, $F$, $s$, $m$ and $d$.

**Proof.** The estimate (2.6) can be found in [43, pp 387–8]. Then for $-\frac{d}{2} < s \le \frac{d}{2}$, Taylor’s formula implies that there exists a sequence $\tilde{H}_i(f_1, \dots, f_m)$ satisfying $\tilde{H}_i(0, \dots, 0) = 0$ and

$$ G(f_1, \dots, f_m) = \sum_{i=1}^{m} \left( \partial_{f_i} G(0, \dots, 0) + \tilde{H}_i(f_1, \dots, f_m) \right) f_i. $$

This together with the product law (2.4) and the estimate (2.6) yields (2.7).

Moreover, we note that

$$ \begin{aligned} G(f_1^1, \dots, f_m^1) - G(f_1^1, \dots, f_m^1) &= \int_0^1 \frac{d}{ds} G(f_1^1 + s(f_1^1 - f_2^1), \dots, f_m^1 + s(f_m^1 - f_m^1)) ds \\ &= \sum_{i=1}^m (f_i^1 - f_i^2) \partial_{f_i} G(0, \dots, 0) + \sum_{i=1}^m (f_i^1 - f_i^2) \int_0^1 (\partial_{f_i} G(f_1^1 + s(f_1^1 - f_2^1), \dots, f_m^1 + s(f_m^1 - f_m^1)) \\ &\qquad - \partial_{f_i} G(0, \dots, 0)) ds. \end{aligned} $$

Therefore, applying (2.3), (2.6) and the embedding $\dot{B}^{\frac{d}{2}} \hookrightarrow L^\infty$, we get (2.8). $\square$

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

Finally, we give optimal regularity estimates of some linear equations. We mention that such estimates on usual Besov norms can be easily extended to the norms restricted in low or high frequencies. We recall the estimates of the heat equation as follows (see [2, p 157] for example).

**Lemma 2.6.** Let $\mu_* > 0$, $s \in \mathbb{R}$ and $1 \le p \le \infty$. For given time $T > 0$, assume $u_0 \in \dot{B}^s$ and $f \in \tilde{L}^p(0, T; \dot{B}^{s-2+\frac{d}{p}})$. If $u$ solves the problem

$$ \left\{ \begin{array}{l} \partial_t u - \mu_* \Delta u = f, \\ u(0, x) = u_0(x), \end{array} \right. $$
then the following estimate is fulfilled:

$$ \|u\|_{\tilde{L}_t^\infty(\dot{B}^s)} + \mu_*^{\frac{1}{p_1}} \|u\|_{\tilde{L}_t^p(\dot{B}^{s+\frac{2}{p}})} \le C \left( \|u_0\|_{\dot{B}^s} + \mu_*^{\frac{1}{p}-1} \|f\|_{\tilde{L}_t^p(\dot{B}^{s-2+\frac{2}{p}})} \right), \quad t \in (0,T), $$

where $C > 0$ is a constant independent of $T$ and $\mu_*$.

We have the regularity estimates of the damped transport equation. Since it can be directly shown by the commutator estimate (2.5) and Grönwall's inequality as in [15, 23], we omit the proof for brevity.

**Lemma 2.7.** Let $\lambda_* \ge 0$, $p=1$ or $\lambda_* > 0$, $1 \le p \le \infty$. For $-\frac{d}{2} < s \le \frac{d}{2} + 1$ and given time $T > 0$, assume that $u_0 \in \dot{B}^s$, $v \in L^1(0,T; \dot{B}^{\frac{d}{2}+1})$ and $f \in \tilde{L}^p(0,T; \dot{B}^s)$. If $u$ solves the problem

$$ \left\{ \begin{array}{l} \partial_t u + v \cdot \nabla u + \lambda_* u = f, \\ u(0,x) = u_0(x), \end{array} \right. $$

then it holds that

$$ \|u\|_{\tilde{L}_t^\infty(\dot{B}^s)} + \lambda_*^{\frac{1}{p}} \|u\|_{\tilde{L}_t^p(\dot{B}^s)} \le C \exp \left( C \|v\|_{\tilde{L}_t^1(\dot{B}^{\frac{d}{2}+1})} \right) \left( \|u_0\|_{\dot{B}^s} + \lambda_*^{\frac{1}{p}-1} \|f\|_{\tilde{L}_t^p(\dot{B}^s)} \right), \quad t \in (0,T), $$

where $C > 0$ is a constant independent of $T$ and $\lambda_*$.

### 3. Analysis of the linearized system

We now consider the linearized problem associated to (1.22), which reads

$$ \left\{ \begin{array}{l} \partial_t y + v \cdot \nabla y = 0, \\ \partial_t w + v \cdot \nabla w + (h_1 + H_1) \operatorname{div} u + (h_2 + H_2) \frac{w}{\varepsilon} = S_1, \\ \partial_t r + v \cdot \nabla r + (h_3 + H_3) \operatorname{div} u = S_2, \\ \partial_t u + v \cdot \nabla u + \frac{u}{\tau} + (h_4 + H_4) \nabla r + (h_5 + H_5) \nabla w = S_3, \\ (y, w, r, u)(0,x) = (y_0, w_0, r_0, u_0)(x), \end{array} \right. \quad (3.1) $$

where $h_i$ ($i = 1, \dots, 5$) are given positive constants and $H_i = H_i(t,x)$ ($i = 1, \dots, 5$), $S_i = S_i(t,x)$ ($i = 1, 2, 3$) are given smooth functions.

We first establish the following a-priori estimate for solutions of the linear problem (3.1) uniformly with respect to the parameters $\varepsilon, \tau$, which improves the result in [15] without the uniformity with respect to $\tau$. As explained before, the threshold $J_\tau$ between low and high frequencies given by (2.1) is the key to our analysis.

**Proposition 3.1.** Let $d \ge 2$, $0 < \varepsilon \le \tau < 1$, $T > 0$, and the threshold $J_\tau$ be given by (2.1). Assume that $(w_0, r_0, u_0) \in \dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1}$, $S_1, S_2, S_3 \in L^1(0, T; \dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1})$, $H_i \in C([0, T]; \dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1})$ and $\partial_t H_i \in L^1(0, T; \dot{B}^{\frac{d}{2}})$ for $i = 1, 2, \dots, 5$. There exists a constant $c > 0$ independent of $T$, $\varepsilon$ and $\tau$ such that if

$$ \mathcal{Z}(t) := \sum_{i=1}^{5} \|H_i\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1})} \le c, \quad t \in (0,T), \quad (3.2) $$
then for $t \in (0, T)$, the solution $(y, w, r, u)$ of the Cauchy problem (3.1) satisfies

$$ 
\begin{aligned}
\mathcal{X}(t) := & \|(y, w, r, u)\|_{\tilde{L}_t^\infty(\dot{B}_2^{d-1} \cap \dot{B}_2^d)} + \|(\partial_t y, \partial_t w, \partial_t r, \partial_t u)\|_{L_t^1(\dot{B}_2^d)} \\
& + \frac{1}{\epsilon}\|w\|_{L_t^1(\dot{B}_2^{d-1} \cap \dot{B}_2^d)} + \frac{1}{\sqrt{\epsilon}}\|w\|_{\tilde{L}_t^2(\dot{B}_2^{d-1} \cap \dot{B}_2^{d+1})} \\
& + \tau\|r\|_{L_t^1(\dot{B}_2^{d+1} \cap \dot{B}_2^d+2)}^\ell + \|r\|_{L_t^1(\dot{B}_2^{d+1})}^h + \tau\|r\|_{L_t^1(\dot{B}_2^{d+1})} + \sqrt{\tau}\|r\|_{\tilde{L}_t^2(\dot{B}_2^d \cap \dot{B}_2^{d+1})} \\
& + \|u\|_{L_t^1(\dot{B}_2^d \cap \dot{B}_2^{d+1})} + \frac{1}{\sqrt{\tau}}\|u\|_{\tilde{L}_t^2(\dot{B}_2^{d-1} \cap \dot{B}_2^{d+1})} + \frac{1}{\tau}\|u + \tau(h_4 + H_4)\nabla r\|_{L_t^1(\dot{B}_2^{d-1} \cap \dot{B}_2^d)} \\
\le & C_0 \exp \left( C_0 \int_0^t \mathcal{V}(s) \, ds \right) \left( \| (y_0, w_0, r_0, u_0) \|_{\dot{B}_2^{d-1} \cap \dot{B}_2^{d+1}} + \| (S_1, S_2, S_3) \|_{L_t^1(\dot{B}_2^{d-1} \cap \dot{B}_2^d+1)} \right),
\end{aligned}
\quad (3.3) $$

where $C_0 > 1$ is a universal constant, and $\mathcal{V}(t)$ is denoted by

$$ \mathcal{V}(t) := \|v(t)\|_{\dot{B}_2^d \cap \dot{B}_2^{d+1}} + \sum_{i=1}^5 \|\partial_t H_i(t)\|_{\dot{B}_2^d}. \quad (3.4) $$

**Proof.** First, we deal with the purely transport unknown y. By the regularity estimate in lemma 2.7 for the transport equation (3.1)$_1$, it follows that

$$ \|y\|_{\tilde{L}_t^\infty(\dot{B}_2^{d-1} \cap \dot{B}_2^d+1)} \lesssim \exp \left( \int_0^t \|v(s)\|_{\dot{B}_2^{d+1} \cap \dot{B}_2^d} \, ds \right) \|y_0\|_{\dot{B}_2^{d-1} \cap \dot{B}_2^d+1}. \quad (3.5) $$

And direct produce law (2.4) for the equation (3.1)$_1$ gives that

$$ \|\partial_t y\|_{L_t^1(\dot{B}_2^d)} \lesssim \int_0^t \|v(s)\|_{\dot{B}_2^d} \|y(s)\|_{\dot{B}_2^{d+1}} \, ds. \quad (3.6) $$

Similarly, we also get from (2.4) and (3.1)$_3$ that

$$ \begin{aligned}
\|\partial_t r\|_{L_t^1(\dot{B}_2^d)} \lesssim & \int_0^t \|v(s)\|_{\dot{B}_2^d} \|r(s)\|_{\dot{B}_2^{d+1}} \, ds \\
& + \left(1 + \|H_3\|_{\tilde{L}_t^\infty(\dot{B}_2^d)}\right) \|u\|_{L_t^1(\dot{B}_2^{d+1})} + \|S_2\|_{L_t^1(\dot{B}_2^d)},
\end{aligned} \quad (3.7) $$

and

$$ 
\begin{aligned}
\| (\partial_t w, \partial_t u) \|_{L_t^1(\dot{B}_2^d)} \lesssim & \int_0^t \|v(s)\|_{\dot{B}_2^d} \|(w, u)(s)\|_{\dot{B}_2^{d+1}} \, ds + \frac{1}{\tau} \|u + \tau(h_5 + H_5)\nabla r\|_{L_t^1(\dot{B}_2^d)} \\
& + \left(1 + \sum_{i=1}^5 \|H_i\|_{\tilde{L}_t^\infty(\dot{B}_2^d)}\right) \left(\frac{1}{\epsilon}\|w\|_{L_t^1(\dot{B}_2^d)} + \|u\|_{L_t^1(\dot{B}_2^{d+1})}\right) \\
& + \|(S_1, S_2)\|_{L_t^1(\dot{B}_2^d)}.
\end{aligned}
\quad (3.8) $$
The conclusion of the proof will follow from Lemmas 3.1–3.3 given and proven in the next three subsections. Indeed, combining (3.5)–(3.8) and the uniform estimates of $(w, r, u)$ from Lemmas 3.1–3.3 together and taking the constant $\eta > 0$ suitable small in Lemma 3.3, we obtain

$$ \begin{aligned} \mathcal{X}(t) \lesssim & \|(S_1, S_2, S_3)\|_{L_1^1(\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d+1})} + (\sqrt{\mathcal{Z}(t)} + \mathcal{Z}(t)) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) ds \\ & + \|(y_0, w_0, r_0, u_0)\|_{\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d+1}}. \end{aligned} $$

Then making use of the Grönwall inequality and the smallness assumption (3.2) of $\mathcal{Z}(t)$, we obtain the uniform $a$-priori estimate (3.3). □

### 3.1. Low-frequency analysis

Motivated by Darcy's law (1.3)$_2$, we introduce the following effective flux

$$ z := u + \tau (h_4 + H_4) \nabla r, \qquad (3.9) $$

which undergoes a purely damped effect in the low-frequency region $|\xi| \le \frac{2^k}{\tau}$ and allows us to diagonalize the subsystem (3.1)$_2$–(3.1)$_4$ up to some higher-order terms that can be absorbed. Indeed, substituting (3.9) into (3.1), we obtain

$$ \left\{ \begin{array}{l} \partial_t w + \frac{h_2}{\varepsilon} w = L_1 + R_1 + S_1, \\ \partial_t r - h_3 h_4 \tau \Delta r = L_2 + R_2 + S_2, \\ \partial_t z + \frac{z}{\tau} = L_3 + R_3 + S_3, \\ (w, r, z)(0,x) = (w_0, r_0, z_0)(x), \end{array} \right. \qquad (3.10) $$

where the higher-order linear terms $L_i$ ($i = 1, 2, 3$) are denoted as

$$ \left\{ \begin{array}{l} L_1 := h_1 (h_4 \tau \Delta r - \operatorname{div} z), \\ L_2 := -h_3 \operatorname{div} z, \\ L_3 := h_3 h_4 \tau \nabla (h_4 \tau \Delta r - \operatorname{div} z) - h_5 \nabla w, \end{array} \right. $$

and the nonlinear terms $R_i$ ($i = 1, 2, 3$) are defined by

$$ \left\{ \begin{array}{l} R_1 := -v \cdot \nabla w - H_1 \operatorname{div} u + h_1 \tau \operatorname{div} (H_4 \nabla r) - \frac{1}{\varepsilon} H_2 w, \\ R_2 := -v \cdot \nabla r - H_3 \operatorname{div} u + h_3 \tau \operatorname{div} (H_4 \nabla r), \\ R_3 := -v \cdot \nabla u - H_5 \nabla w - \tau \partial_t (H_4 \nabla r) + h_4 \tau \nabla R_2. \end{array} \right. \qquad (3.11) $$

Now, to establish the $\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d}$-estimates in low frequencies to the solutions of system (3.1) uniformly with respect to both $\varepsilon$ and $\tau$, we understand the equations in (3.10) are decoupled. More precisely, we will treat the equations of $w$ and $z$ as damped equations and $r$ as a heat equation, respectively. This viewpoint plays a key role in the proof of the following lemma.
**Lemma 3.1.** Let $T > 0$, and the threshold $J_{\tau}$ be given by (2.1). Then for $t \in (0, T)$, the solution $(w, r, u)$ to the linear problem (3.1)$_2$–(3.1)$_4$ satisfies

$$ \begin{aligned} & \|(w, r, u)\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})}^\ell + \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})}^\ell + \frac{1}{\sqrt{\varepsilon}} \|w\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})}^\ell \\ & \qquad + \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1} \cap \dot{B}^{\frac{d}{2}+2})}^\ell + \sqrt{\tau} \|r\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}} \cap \dot{B}^{\frac{d}{2}+1})}^\ell \\ & \qquad + \|u\|_{L_t^1(\dot{B}^{\frac{d}{2}} \cap \dot{B}^{\frac{d}{2}+1})}^\ell + \frac{1}{\sqrt{\tau}} \|u\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})}^\ell + \frac{1}{\tau} \|z\|_{L_t^1(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})}^\ell \quad (3.12) \\ & \lesssim \|(w_0, r_0, u_0)\|_{\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}}} + \|(S_1, S_2, S_3)\|_{L_t^1(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})}^\ell \\ & \qquad + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds, \end{aligned} $$

where $\mathcal{Z}(t)$, $\mathcal{X}(t)$, $\mathcal{V}(t)$ and $z$ are defined in (3.2)–(3.4) and (3.9), respectively.

**Remark 3.1.** In [15], the authors obtained the low-frequency estimates by constructing a related Lyapunov functional. However, that method does not lead to the desired estimates which uniform with respect to $\tau$. Moreover, it should be noted that the effective unknown $z$ given by (3.9) enables us to capture the heat-like behavior of the unknown $r$ in low frequencies directly, which is consistent with the parabolic nature of the limiting porous media equations.

**3.1.1. The $\dot{B}^{\frac{d}{2}}$-estimates.** We first perform $\dot{B}^{\frac{d}{2}}$-estimates in low frequencies for the heat equation (3.10)$_2$. It follows from the regularity estimate in lemma 2.6 that

$$ \begin{aligned} & \|r\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})}^\ell + \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+2})}^\ell \\ & \lesssim \|r_0\|_{\dot{B}^{\frac{d}{2}}}^\ell + \|L_2\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|R_2\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|S_2\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \quad (3.13) \\ & \lesssim \|r_0\|_{\dot{B}^{\frac{d}{2}}}^\ell + 2^{J_\tau} \|z\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|R_2\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|S_2\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell. \end{aligned} $$

Applying Lemma 2.7 to the damped equation (3.10)$_1$, we get

$$ \begin{aligned} & \|w\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})}^\ell + \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \\ & \lesssim \|w_0\|_{\dot{B}^{\frac{d}{2}}}^\ell + \|L_1\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|R_1\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|S_1\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \\ & \lesssim \|w_0\|_{\dot{B}^{\frac{d}{2}}}^\ell + \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+2})}^\ell + 2^{J_\tau} \|z\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|R_1\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|S_1\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \\ & \lesssim \|(w_0, r_0)\|_{\dot{B}^{\frac{d}{2}}}^\ell + 2^{J_\tau} \|z\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|(R_1, R_2)\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|(S_1, S_2)\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell, \end{aligned} $$

where we used the inequality (3.13) to control terms involving $r$ in the equation (3.10)$_1$.
Similarly, by virtue of the inequality (3.13) and Lemmas 2.6 and 2.7, we have for the equation (3.10)$_3$ that

$$
\begin{align}
\|z\|_{\tilde{L}_t^\infty(\dot{B}^d_2)}^\ell + \frac{1}{\tau} \|z\|_{L_t^1(\dot{B}^d_2)}^\ell &\lesssim \|z_0\|_{\dot{B}^d_2}^\ell + \|L_3\|_{L_t^1(\dot{B}^d_2)}^\ell + \|(R_3, S_3)\|_{L_t^1(\dot{B}^d_2)}^\ell \nonumber \\
&\lesssim \|z_0\|_{\dot{B}^d_2}^\ell + 2^{J_\tau} \|w\|_{L_t^1(\dot{B}^d_2)}^\ell + \tau^2 2^{J_\tau} \|r\|_{L_t^1(\dot{B}^{d+2}_2)}^\ell \nonumber \\
&\quad + \tau 2^{2J_\tau} \|z\|_{L_t^1(\dot{B}^{d+2}_2)}^\ell + \|(R_3, S_3)\|_{L_t^1(\dot{B}^d_2)}^\ell. \tag{3.15}
\end{align}
$$

Since the threshold $J_\tau$ satisfies the condition (2.1), thus $\tau 2^{J_\tau} \sim 2^k \ll 1$ for suitable negative integer $k$. Due to the condition $\varepsilon \le \tau$ so that $2^{J_\tau} \|w\|_{L_t^1(\dot{B}^d_2)}^\ell \le \frac{2^{k+1}}{\varepsilon} \|w\|_{L_t^1(\dot{B}^d_2)}^\ell$, we have by the inequalities (3.13)–(3.15) that

$$
\begin{align}
\| (w, r, z) \|_{\tilde{L}_t^\infty(\dot{B}^d_2)}^\ell + \tau \|r\|_{L_t^1(\dot{B}^{d+2}_2)}^\ell + \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^d_2)}^\ell + \frac{1}{\tau} \|z\|_{L_t^1(\dot{B}^d_2)}^\ell \tag{3.16} \\
\lesssim \| (w_0, r_0, z_0) \|_{\dot{B}^d_2}^\ell + \|(R_1, R_2, R_3)\|_{L_t^1(\dot{B}^d_2)}^\ell + \|(S_1, S_2, S_3)\|_{L_t^1(\dot{B}^d_2)}^\ell. \nonumber
\end{align}
$$

The terms on the right-hand side of (3.16) can be estimated as follows. First, one derives from the inequality (2.2), the product law (2.4) and the composition estimate (2.6) that

$$
\|z_0\|_{\dot{B}^d_2}^\ell \lesssim \|(r_0, u_0)\|_{\dot{B}^d_2}^\ell + \|H_4(0)\|_{\dot{B}^d_2} \|r_0\|_{\dot{B}^d_2} \lesssim \|(r_0, u_0)\|_{\dot{B}^d_2} \quad (3.17)
$$

By the product law (2.4) again, we also get

$$
\left\{
\begin{array}{l}
\displaystyle \left\| v \cdot \nabla w \right\|_{L_t^1(\dot{B}^d_2)}^\ell \lesssim \int_0^t \|v(s)\|_{\dot{B}^d_2} \|w(s)\|_{\dot{B}^{d+1}_2} \, ds, \\[2ex]
\displaystyle \left\| H_1 \operatorname{div} u \right\|_{L_t^1(\dot{B}^d_2)}^\ell \lesssim \left\| H_1 \right\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \left\| u \right\|_{L_t^1(\dot{B}^{d+1}_2)}, \\[2ex]
\displaystyle \frac{1}{\varepsilon} \left\| H_2 w \right\|_{L_t^1(\dot{B}^d_2)}^\ell \lesssim \left\| H_2 \right\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \frac{1}{\varepsilon} \left\| w \right\|_{L_t^1(\dot{B}^d_2)}.
\end{array}
\right.
\qquad (3.18)
$$

According to (2.2) and (2.3)–(2.4), the tricky nonlinear term $H_4 \nabla r$ in (3.11) can be estimated as

$$
\begin{align}
\tau \| \operatorname{div} (H_4 \nabla r) \|_{L_t^1(\dot{B}^d_2)}^\ell &\lesssim \tau \|H_4 \nabla r^\ell\|_{L_t^1(\dot{B}^{d+1}_2)}^\ell + \tau \|H_4 \nabla r^h\|_{L_t^1(\dot{B}^{d+1}_2)}^\ell \nonumber \\
&\lesssim \tau \|H_4 \nabla r^\ell\|_{L_t^1(\dot{B}^{d+1}_2)}^\ell + \|H_4 \nabla r^h\|_{L_t^1(\dot{B}^d_2)}^\ell \tag{3.19} \\
&\lesssim \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \tau \|r\|_{L_t^1(\dot{B}^{d+2}_2)}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^{d+1}_2)} \tau \|r\|_{L_t^1(\dot{B}^{d+1}_2)}^\ell \nonumber \\
&\quad + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \|r\|_{L_t^1(\dot{B}^{d+1}_2)}^h \nonumber
\end{align}
$$

**Remark 3.2.** The above estimate (3.19) for $H_4 \nabla r$ arising from two pressures implies that one needs uniform $\dot{B}^{d-1}_2$-estimates for low frequencies. Indeed, as $H_4$ does not have the either $L^1$-in-time or $L^2$-in-time integrability property, the product law (2.3) in $\dot{B}^{d+1}_2$ indicates us to
discover the control of $\tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^\ell$, which can not be obtained from the $\dot{B}^{\frac{d}{2}}$-estimates in this section.

**Remark 3.3.** It is also one of the reasons why we need to perform the $\dot{B}^{\frac{d}{2}+1}$-estimates in the both low and high frequencies in the later section 3.3. Indeed, in the low-frequency setting, the uniform $\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})$-norm is not enough to produce the uniform $\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}+1})$-estimates required in (3.19) due to the inclusion (2.2).

Now, one derives from inequalities (3.17)–(3.19) that

$$ \begin{aligned} \|R_1\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell &\lesssim \|v \cdot \nabla w\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|H_1 \operatorname{div} u\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \\ &\quad + \tau \|H_4 \nabla r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^\ell + \frac{1}{\varepsilon} \|H_2 w\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \\ &\lesssim \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds. \end{aligned} \quad (3.20) $$

Similarly, we have

$$ \begin{aligned} \|R_2\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell &\lesssim \|v \cdot \nabla r\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|H_3 \operatorname{div} u\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \tau \|H_4 \nabla r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^\ell \\ &\lesssim \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds. \end{aligned} \quad (3.21) $$

To estimate $R_3$, we notice that (2.2) together with (2.4) implies

$$ \begin{aligned} \tau \|\partial_t (H_4 \nabla r)\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell &\lesssim \|\partial_t (H_4 \nabla r)\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell \\ &\lesssim \int_0^t \|\partial_t H_4(s)\|_{\dot{B}^{\frac{d}{2}}} \|r(s)\|_{\dot{B}^{\frac{d}{2}}} \, ds + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \|\partial_t r\|_{L_t^1(\dot{B}^{\frac{d}{2}})} \\ &\lesssim \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds. \end{aligned} \quad (3.22) $$

and

$$ \|H_5 \nabla w\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \lesssim \frac{1}{\tau} \|H_5 \nabla w\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell \lesssim \|H_5\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}})} $$

where we used the assumption $\varepsilon \le \tau$. Thus, it holds that

$$ \begin{aligned} \|R_3\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell &\lesssim \int_0^t \|v(s)\|_{\dot{B}^{\frac{d}{2}}} \|u(s)\|_{\dot{B}^{\frac{d}{2}+1}} \, ds + \tau \|H_4 \nabla r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^\ell + \|H_5 \nabla w\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \\ &\quad + \tau \|\partial_t (H_4 \nabla r)\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell + \|R_2\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell \\ &\lesssim \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds. \end{aligned} \quad (3.23) $$
We substitute inequalities (3.17), (3.20), (3.21) and (3.23) into inequality (3.16) and use standard interpolation to get

$$ \begin{aligned} & \| (w, r, z) \|_{\tilde{L}_t^\infty (\dot{B}^d_2)}^\ell + \tau \|r\|_{L_t^1 (\dot{B}^{d+2}_2)}^\ell + \sqrt{\tau} \|r\|_{\tilde{L}_t^2 (\dot{B}^{d+1}_2)}^\ell \\ & + \frac{1}{\varepsilon} \|w\|_{L_t^1 (\dot{B}^d_2)}^\ell + \frac{1}{\sqrt{\varepsilon}} \|w\|_{\tilde{L}_t^2 (\dot{B}^d_2)}^\ell + \frac{1}{\tau} \|z\|_{L_t^1 (\dot{B}^d_2)}^\ell \\ & \lesssim \| (w_0, r_0, u_0) \|_{\dot{B}^d_2} + \| (S_1, S_2, S_3) \|_{L_t^1 (\dot{B}^d_2)}^\ell + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds. \end{aligned} \quad (3.24) $$

Thence, we rewrite the form (3.9) and use inequalities (2.2) and (3.19) to obtain the $L_t^1(\dot{B}^d_2)$-estimate of $u$ as follows:

$$ \begin{aligned} \|u\|_{L_t^1(\dot{B}^d_2+1)}^\ell & \lesssim \|z\|_{L_t^1(\dot{B}^d_2+1)}^\ell + \tau \|\nabla r\|_{L_t^1(\dot{B}^d_2+1)}^\ell + \tau \|H_4 \nabla r\|_{L_t^1(\dot{B}^d_2+1)}^\ell \\ & \lesssim \frac{1}{\tau} \|z\|_{L_t^1(\dot{B}^d_2)}^\ell + \tau \|r\|_{L_t^1(\dot{B}^d_2+2)}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \tau \|r\|_{L_t^1(\dot{B}^d_2+2)}^\ell \\ & \quad + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^d_2+1)} \tau \|r\|_{L_t^1(\dot{B}^d_2+1)}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \|r\|_{L_t^1(\dot{B}^d_2+1)}^h. \end{aligned} $$

Similarly, we have

$$ \|u\|_{\tilde{L}_t^\infty(\dot{B}^d_2)}^\ell \lesssim \|(z, r)\|_{\tilde{L}_t^\infty(\dot{B}^d_2)}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \|r\|_{\tilde{L}_t^\infty(\dot{B}^d_2)}, $$

and

$$ \frac{1}{\sqrt{\tau}} \|u\|_{\tilde{L}_t^2(\dot{B}^d_2)}^\ell \lesssim \frac{1}{\sqrt{\tau}} \|z\|_{\tilde{L}_t^2(\dot{B}^d_2)}^\ell + \sqrt{\tau} \|r\|_{\tilde{L}_t^2(\dot{B}^d_2+1)}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \sqrt{\tau} \|r\|_{\tilde{L}_t^2(\dot{B}^d_2+1)}^\ell. $$

We thus obtain from inequality (3.24) that

$$ \begin{aligned} \|u\|_{\tilde{L}_t^\infty(\dot{B}^d_2)}^\ell + \frac{1}{\sqrt{\tau}} \|u\|_{\tilde{L}_t^2(\dot{B}^d_2)}^\ell + \|u\|_{L_t^1(\dot{B}^d_2+1)}^\ell & \quad (3.25) \\ \lesssim \| (w_0, r_0, u_0) \|_{\dot{B}^d_2} + \| (S_1, S_2, S_3) \|_{L_t^1(\dot{B}^d_2)}^\ell + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds. \end{aligned} $$

### 3.1.2. The $\dot{B}^{d-1}_2$-estimates.
We perform the $\dot{B}^{d-1}_2$-estimates so as to control $\tau \|r\|_{L_t^1(\dot{B}^{d+1}_2)}^\ell$, as explained in Remark 3.2. Arguing similarly as for inequalities (3.13)–(3.16), we have

$$ \begin{aligned} & \| (w, r, z) \|_{\tilde{L}_t^\infty (\dot{B}^d_2 - 1)}^\ell + \tau \|r\|_{L_t^1 (\dot{B}^d_2 + 1)}^\ell + \frac{1}{\varepsilon} \|w\|_{L_t^1 (\dot{B}^d_2 - 1)}^\ell + \frac{1}{\tau} \|z\|_{L_t^1 (\dot{B}^d_2 - 1)}^\ell \\ & \lesssim \| (w_0, r_0, u_0) \|_{\dot{B}^d_2 - 1} + \| (R_1, R_2, R_3) \|_{L_t^1 (\dot{B}^d_2 - 1)}^\ell + \| (S_1, S_2, S_3) \|_{L_t^1 (\dot{B}^d_2 - 1)}^\ell. \end{aligned} \quad (3.26) $$
Direct calculations give

$$
\begin{align}
& \|(R_1, R_2)\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell \nonumber \\
& \lesssim \int_0^t \|v(s)\|_{\dot{B}^{\frac{d}{2}}} \|(w, r)(s)\|_{\dot{B}^{\frac{d}{2}}} \mathrm{d}s + \|(H_1, H_3)\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \|u\|_{L_t^1(\dot{B}^{\frac{d}{2}})} \nonumber \\
& \qquad + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} + \|H_2\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \tag{3.27} \\
& \lesssim \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \mathrm{d}s. \nonumber
\end{align}
$$

By inequalities (3.22), (3.27) and product law (2.4) for $d \ge 2$, the term $R_3$ can be bounded by

$$
\begin{align}
\|R_3\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell &\lesssim \int_0^t \|v(s)\|_{\dot{B}^{\frac{d}{2}}} \|u(s)\|_{\dot{B}^{\frac{d}{2}}} \mathrm{d}s + \|H_5\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}})} \nonumber \\
&\quad + \tau \|\partial_t (H_4 \nabla r)\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell + \|R_2\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell \tag{3.28} \\
&\lesssim \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \mathrm{d}s. \nonumber
\end{align}
$$

Inserting (3.27) and (3.28) into (3.26) and taking advantage of interpolation, we obtain

$$
\begin{align}
\|(w, r, z)\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})}^\ell &+ \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^\ell + \sqrt{\tau} \|r\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}})}^\ell \nonumber \\
&\quad + \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell + \frac{1}{\sqrt{\varepsilon}} \|w\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}-1})}^\ell + \frac{1}{\tau} \|z\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell \tag{3.29} \\
&\lesssim \|(w_0, r_0, u_0)\|_{\dot{B}^{\frac{d}{2}-1}} + \|(S_1, S_2, S_3)\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \mathrm{d}s. \nonumber
\end{align}
$$

This together with inequality (2.2) and the fact that $u = z - \tau(h_4 + H_4)\nabla r$ leads to

$$
\|u\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})}^\ell \lesssim \|(z, r)\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \|r\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})}.
$$

Similarly, one gets

$$
\frac{1}{\sqrt{\tau}} \|u\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}-1})}^\ell \lesssim \frac{1}{\sqrt{\tau}} \|z\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}-1})}^\ell + \sqrt{\tau} \|r\|_{\tilde{L}_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \sqrt{\tau} \|r\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}})},
$$

and

$$
\|u\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \lesssim \frac{1}{\tau} \|z\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell + \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^\ell + \|H_4\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^\ell.
$$

Combining the above three estimates, we are led to

$$
\begin{align}
& \|u\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})}^\ell + \frac{1}{\sqrt{\tau}} \|u\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}-1})}^\ell + \|u\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell \nonumber \\
& \lesssim \|(w_0, r_0, u_0)\|_{\dot{B}^{\frac{d}{2}-1}} + \|(S_1, S_2, S_3)\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^\ell + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \mathrm{d}s. \tag{3.30}
\end{align}
$$
Putting the estimates (3.12), (3.24), (3.25) and (3.29), (3.30) together, we complete the proof of Lemma 3.1.

## 3.2. High-frequency analysis

In this section, we establish some uniform high-frequency estimates of solutions to the linear problem (3.1) in terms of the Lyapunov functional. More precisely, we establish the $\tilde{L}_t^\infty(\dot{B}_2^{d-1} \cap \dot{B}_2^d)$-estimates, and furthermore obtain the control of higher-order $L_t^1(\dot{B}_2^d \cap \dot{B}_2^{d+1})$-norms.

**Lemma 3.2.** Let $T > 0$, and the threshold $J_\tau$ be given by (2.1). Then for any $t \in (0, T)$, the solution $(w, r, u)$ to the linear problem (3.1)$_2$–(3.1)$_4$ satisfies

$$ \begin{aligned} & \|(w, r, u)\|_{\tilde{L}_t^\infty(\dot{B}_2^{d-1} \cap \dot{B}_2^d)}^h + \|(w, r, u)\|_{L_t^1(\dot{B}_2^{d+1})}^h \\ & \quad + \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}_2^{d-1} \cap \dot{B}_2^d)}^h + \frac{1}{\sqrt{\varepsilon}} \|w\|_{\tilde{L}_t^2(\dot{B}_2^{d-1} \cap \dot{B}_2^d)}^h \\ & \quad + \sqrt{\tau} \|r\|_{\tilde{L}_t^2(\dot{B}_2^{d+1})}^h + \frac{1}{\tau} \|z\|_{L_t^1(\dot{B}_2^{d-1} \cap \dot{B}_2^d)}^h \\ & \lesssim \|(w_0, r_0, u_0)\|_{\dot{B}_2^{d+1}}^h + \|(S_1, S_2, S_3)\|_{L_t^1(\dot{B}_2^{d+1})}^h \\ & \quad + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds. \end{aligned} \quad (3.31) $$

**Proof.** To prove Lemma 3.2, we localize (3.1)$_2$–(3.1)$_4$ in frequencies as follows

$$ \left\{ \begin{aligned} & \partial_t w_j + v \cdot \nabla w_j + (h_1 + H_1) \operatorname{div} u_j + (h_2 + H_2) \frac{w_j}{\varepsilon} = \dot{\Delta}_j S_1 + T_j^1, \\ & \partial_t r_j + v \cdot \nabla r_j + (h_3 + H_3) \operatorname{div} u_j = \dot{\Delta}_j S_2 + T_j^2, \\ & \partial_t u_j + v \cdot \nabla u_j + \frac{u_j}{\tau} + (h_4 + H_4) \nabla r_j + (h_5 + H_5) \nabla w_j = \dot{\Delta}_j S_3 + T_j^3, \end{aligned} \right. \quad (3.32) $$

with the commutator terms

$$ \left\{ \begin{aligned} T_j^1 &:= [v, \dot{\Delta}_j] \nabla w + [H_1, \dot{\Delta}_j] \operatorname{div} u + \frac{1}{\varepsilon} [H_2, \dot{\Delta}_j] w, \\ T_j^2 &:= [v, \dot{\Delta}_j] \nabla r + [H_3, \dot{\Delta}_j] \operatorname{div} u, \\ T_j^3 &:= [v, \dot{\Delta}_j] \nabla u + [H_4, \dot{\Delta}_j] \nabla r + [H_5, \dot{\Delta}_j] \nabla w. \end{aligned} \right. \quad (3.33) $$

Multiplying (3.32)$_3$ by $u_j$ and integrating the resulting equation by parts, we get

$$ \begin{aligned} & \frac{d}{dt} \int_{\mathbb{R}^d} \frac{1}{2} |u_j|^2 dx + \int_{\mathbb{R}^d} \frac{1}{\tau} |u_j|^2 dx - \int_{\mathbb{R}^d} ((h_4 + H_4) r_j \operatorname{div} u_j + (h_5 + H_5) w_j \operatorname{div} u_j) dx \\ & \lesssim \| \operatorname{div} v \|_{L^\infty} \| u_j \|_{L^2}^2 + \left( \| \dot{\Delta}_j S_3 \|_{L^2}^2 + \| T_j^3 \|_{L^2}^2 \right) \| u_j \|_{L^2}^2 + \| \nabla H_4 \|_{L^\infty} \| (w_j, r_j) \|_{L^2}^2 \| u_j \|_{L^2}^2. \end{aligned} \quad (3.34) $$

Thence, we multiply (3.32)$_1$ by $\frac{h_5+H_5}{h_1+H_1} w_j$ and integrate the resulting equation by parts to show
$$ \begin{aligned} & \frac{d}{dt} \int_{\mathbb{R}^d} \frac{1}{2} \frac{h_5 + H_5}{h_1 + H_1} |w_j|^2 dx + \int_{\mathbb{R}^d} \left( (h_5 + H_5) w_j \operatorname{div} u_j + \frac{(h_2 + H_2)(h_5 + H_5)}{\varepsilon(h_1 + H_1)} |w_j|^2 \right) dx \\ & \lesssim \left( \left\| \partial_t \left( \frac{h_5 + H_5}{h_1 + H_1} \right) \right\|_{L^\infty} + \left\| \frac{h_5 + H_5}{h_1 + H_1} \right\|_{L^\infty} \| \operatorname{div} v \|_{L^\infty} + \left\| \nabla \left( \frac{h_5 + H_5}{h_1 + H_1} \right) \right\|_{L^\infty} \| v \|_{L^\infty} \right) \|w_j\|_{L^2}^2 \\ & \qquad + \left\| \frac{h_5 + H_5}{h_1 + H_1} \right\|_{L^\infty} \left( \| \dot{\Delta}_j S_1 \|_{L^2} + \| T_j^1 \|_{L^2} \right) \|w_j\|_{L^2}. \end{aligned} \quad (3.35) $$

Similarly, direct computations on (3.32)$_2$ yield

$$ \begin{aligned} & \frac{d}{dt} \int_{\mathbb{R}^d} \frac{1}{2} \frac{h_4 + H_4}{h_3 + H_3} |r_j|^2 dx + \int_{\mathbb{R}^d} (h_4 + H_4) r_j \operatorname{div} u_j dx \\ & \lesssim \left( \left\| \partial_t \left( \frac{h_4 + H_4}{h_3 + H_3} \right) \right\|_{L^\infty} + \left\| \frac{h_4 + H_4}{h_3 + H_3} \right\|_{L^\infty} \| \operatorname{div} v \|_{L^\infty} + \left\| \nabla \left( \frac{h_4 + H_4}{h_3 + H_3} \right) \right\|_{L^\infty} \| v \|_{L^\infty} \right) \|r_j\|^2 \\ & \qquad + \left\| \frac{h_4 + H_4}{h_3 + H_3} \right\|_{L^\infty} \left( \| \dot{\Delta}_j S_2 \|_{L^2} + \| T_j^2 \|_{L^2} \right) \|r_j\|_{L^2}. \end{aligned} \quad (3.36) $$

To derive the cross estimate and capture the dissipative property of $r_j$, we gain by taking the $L^2$-inner product of (3.32)$_3$ with $\nabla r_j$ that

$$ \begin{aligned} & \int_{\mathbb{R}^d} \partial_t u_j \cdot \nabla r_j dx + \int_{\mathbb{R}^d} (h_4 + H_4) |\nabla r_j|^2 dx + \int_{\mathbb{R}^d} \left( (h_5 + H_5) \nabla w_j \cdot \nabla r_j + \frac{1}{\tau} u_j \cdot \nabla r_j \right) dx \\ & \lesssim \left( \|v\|_{L^\infty} \|\nabla u_j\|_{L^2} + \|\dot{\Delta}_j S_3\|_{L^2} + \|T_j^3\|_{L^2} \right) \|\nabla r_j\|_{L^2}, \end{aligned} \quad (3.37) $$

and taking the $L^2$-inner product of (3.32)$_2$ with $\operatorname{div} u_j$ that

$$ \begin{aligned} & \int_{\mathbb{R}^d} u_j \cdot \nabla \partial_t r_j dx - \int_{\mathbb{R}^d} (h_3 + H_3) |\operatorname{div} u_j|^2 dx \\ & \lesssim \left( \|v\|_{L^\infty} \|\nabla r_j\|_{L^2} + \|\dot{\Delta}_j S_2\|_{L^2} + \|T_j^2\|_{L^2} \right) \|\operatorname{div} u_j\|_{L^2}. \end{aligned} \quad (3.38) $$

In the spirit of the work [3] by Beauchard and Zuazua, for a small constant $\eta_* > 0$ to be determined, we define the following Lyapunov functional with nonlinear weights as

$$ \mathcal{L}_j(t) := \int_{\mathbb{R}^d} \left( \frac{h_5 + H_5}{h_1 + H_1} |w_j|^2 + \frac{h_4 + H_4}{h_3 + H_3} |r_j|^2 + |u_j|^2 \right) dx + \frac{\eta_*}{\tau} 2^{-2j} \int_{\mathbb{R}^d} u_j \cdot \nabla r_j dx, $$

and its dissipation rate

$$ \mathcal{H}_j(t) := \int_{\mathbb{R}^d} \left( \frac{1}{\tau} |u_j|^2 + \frac{(h_2 + H_2)(h_5 + H_5)}{\varepsilon(h_1 + H_1)} |w_j|^2 \right) dx \\ + \frac{\eta_*}{\tau} 2^{-2j} \int_{\mathbb{R}^d} \left( (h_4 + H_4) |\nabla r_j|^2 + (h_5 + H_5) \nabla w_j \cdot \nabla r_j + \frac{1}{\tau} u_j \cdot \nabla r_j \right) dx. $$

One derives from assumption (3.2) and the embedding $\dot{B}^{\frac{d}{2}} \hookrightarrow L^\infty$ that

$$ \|H_i\|_{L^\infty(L^\infty)} + \|\nabla H_i\|_{L^\infty(L^\infty)} \lesssim \|H_i\|_{L^\infty(\dot{B}^{\frac{d}{2} \cap \dot{B}^{\frac{d}{2}+1})} \lesssim c \ll 1, \quad (3.39) $$
which together with estimates (3.34)–(3.38) and the fact that $2^{-j} \lesssim \tau \le 1$ for any $j \ge J_{\tau} - 1$ yields the following Lyapunov inequality:

$$ \frac{d}{dt} \mathcal{L}_j(t) + \mathcal{H}_j(t) \lesssim \left( \| \operatorname{div} v \|_{L^\infty} + \|v\|_{L^\infty} + \sum_{i=1}^5 \| \partial_i H_i \|_{L^\infty} \right) \|(r_j, w_j, u_j)\|_{L^2}^2 \\ + \left( \sum_{i=1}^5 \| \partial_i H_i \|_{L^\infty} + \| \dot{\Delta}_j (S_1, S_2, S_3) \|_{L^2} + \| (T_j^1, T_j^2, T_j^3) \|_{L^2} \right) \\ \times \|(r_j, w_j, u_j)\|_{L^2}. \quad (3.40) $$

It follows from the smallness condition (3.39), Bernstein's inequality in Lemma 2.1 and the fact $2^{-j} \lesssim \tau$ that

$$ (1 - \eta_*) \|(w_j, r_j, u_j)\|_{L^2}^2 \lesssim \mathcal{L}_j(t) \lesssim (1 + \eta_*) \|(w_j, r_j, u_j)\|_{L^2}^2, $$

and

$$ \mathcal{H}_j(t) \gtrsim \frac{1}{\tau} \|u_j\|_{L^2}^2 + \frac{1}{\varepsilon} \|w_j\|_{L^2}^2 + \frac{\eta_*}{\tau} 2^{-2j} \left( \|\nabla r_j\|_{L^2} - \|\nabla w_j\|_{L^2}^2 - \frac{1}{\tau^2} \|u_j\|_{L^2} \right) \\ \gtrsim \frac{1}{\tau} (1 - \eta_*) \|u_j\|_{L^2}^2 + \frac{1}{\varepsilon} (1 - \eta_*) \|w_j\|_{L^2}^2 + \frac{\eta_*}{\tau} \|r_j\|_{L^2}, $$

where we used the condition $\varepsilon \le \tau$. Thus, we can choose a sufficiently small constant $\eta_*$ > 0 independent of $\varepsilon$ and $\tau$ so that

$$ \mathcal{L}_j(t) \sim \|(w_j, r_j, u_j)\|_{L^2}^2, \quad \mathcal{H}_j(t) \gtrsim \frac{1}{\tau} \|(w_j, r_j, u_j)\|_{L^2}^2 \gtrsim \frac{1}{\tau} \mathcal{L}_j(t). \quad (3.41) $$

Dividing the two sides of (3.40) by $\sqrt{\mathcal{L}_j(t) + \eta}$ for any $\eta > 0$, we have

$$ \frac{d}{dt} \sqrt{\mathcal{L}_j(t) + \eta} + \frac{1}{\tau} \sqrt{\mathcal{L}_j(t) + \eta} - \frac{\eta}{\tau \sqrt{\mathcal{L}_j(t) + \eta}} \\ \lesssim \left( \| \operatorname{div} v \|_{L^\infty} + \|v\|_{L^\infty} + \sum_{i=1}^5 \| \partial_i H_i \|_{L^\infty} \right) \|(r_j, w_j, u_j)\|_{L^2}^2 \\ + \| \dot{\Delta}_j (S_1, S_2, S_3) \|_{L^2} + \| (T_j^1, T_j^2, T_j^3) \|_{L^2}. $$

From (3.41) and the embedding $\dot{B}_{\frac{d}{2}}^d \to L^\infty$ we get after integrating the above inequality over $[0, t]$ and taking the limit as $\eta \to 0$ that

$$ \tau \|(w, r, u)\|_{\tilde{L}_t^\infty(\dot{B}_{\frac{d}{2}+1}^d)}^h + \|(w, r, u)\|_{L_t^1(\dot{B}_{\frac{d}{2}+1}^d)}^h \\ \lesssim \tau \|(w_0, r_0, u_0)\|_{\dot{B}_{\frac{d}{2}+1}^d}^h + \int_0^t \left( \|v(s)\|_{\dot{B}_{\frac{d}{2}}^d \cap \dot{B}_{\frac{d}{2}+1}^d} + \sum_{i=1}^5 \|\partial_i H_i(s)\|_{\dot{B}_{\frac{d}{2}}^d} \right) \tau \|(w, r, u)(s)\|_{\dot{B}_{\frac{d}{2}+1}^d}^h ds \\ + \tau \sum_{j \ge J_\tau - 1} 2^{j(\frac{d}{2}+1)} \|(T_j^1, T_j^2, T_j^3)\|_{L_t^1(L^2)} + \tau \|(S_1, S_2, S_3)\|_{L_t^1(\dot{B}_{\frac{d}{2}+1}^d)}^h. \quad (3.42) $$
According to the commutator estimate (2.5), it follows that

$$
\begin{aligned}
\tau \sum_{j \ge J_{\tau}-1} 2^{j(\frac{d}{2}+1)} \|(T_j^1, T_j^2, T_j^3)\|_{L_t^1(L^2)} & \\
& \lesssim \int_0^t \|v(s)\|_{\dot{B}^{\frac{d}{2}+1}} \|(w, r, u)(s)\|_{\dot{B}^{\frac{d}{2}+1}} \mathrm{d}s + \sum_{i=1}^4 \|H_i\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}+1})} \\
& \times \left( \frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}})} + \tau \|r\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} + \|u\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} \right) \\
& \lesssim \mathcal{Z}(t) \mathcal{X}(t).
\end{aligned}
$$

This with inequality (3.42) leads to

$$
\begin{aligned}
\tau \|(w,r,u)\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}+1})}^h + \|(w,r,u)\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^h & \\
\lesssim \tau \|(w_0,r_0,u_0)\|_{\dot{B}^{\frac{d}{2}+1}}^h + \tau \|(S_1,S_2,S_3)\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^h & + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \mathrm{d}s. \quad (3.43)
\end{aligned}
$$

On the other hand, for any $\eta > 0$, we deduce from inequality (3.35) that

$$
\begin{aligned}
& \frac{d}{dt} \sqrt{\|w_j\|_{L^2}^2 + \eta} + \frac{1}{\varepsilon} \sqrt{\|w_j\|_{L^2}^2 + \eta} \lesssim 2^j \|u_j\|_{L^2} + \| (\partial_t H_1, \partial_t H_5) \|_{L^\infty} \|w_j\|_{L^2} \\
& \qquad + \| \operatorname{div} v \|_{L^\infty} \|w_j\|_{L^2} + \| v \|_{L^\infty} \|w_j\|_{L^2} + \| \dot{\Delta}_j S_1 \|_{L^2} + \| T_j^1 \|_{L^2},
\end{aligned}
$$

which together with (3.43) implies

$$
\begin{aligned}
\frac{1}{\varepsilon} \|w\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^h & \lesssim \|w_0\|_{\dot{B}^{\frac{d}{2}}}^h + \|u\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^h + \|S_1\|_{L_t^1(\dot{B}^{\frac{d}{2}})}^h + \tau \sum_{j \ge J_{\tau}-1} 2^{j(\frac{d}{2}+1)} \|T_j^1\|_{L_t^1(L^2)} \\
& \quad + \int_0^t \left( \| (\partial_t H_1, \partial_t H_5)(s) \|_{\dot{B}^{\frac{d}{2}}} + \|v(s)\|_{\dot{B}^{\frac{d}{2}} \cap \dot{B}^{\frac{d}{2}+1}} \right) \|w(s)\|_{\dot{B}^{\frac{d}{2}}} \mathrm{d}s \\
& \lesssim \tau \|(w_0, r_0, u_0)\|_{\dot{B}^{\frac{d}{2}+1}}^h + \tau \|(S_1, S_2, S_3)\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}^h + \mathcal{Z}(t) \mathcal{X}(t) + \int_0^t \mathcal{V}(s) \mathcal{X}(s) \mathrm{d}s.
\end{aligned}
$$

Thanks to inequality (2.2), one has

$$
\| (w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})}^h \lesssim \tau \| (w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}})}^h \lesssim \tau \| (w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}+1})}^h. \quad (3.44)
$$

Finally, the remain estimates in (3.31) can be achieved similarly to (3.44). We omit the details
here and complete the proof of Lemma 3.2. $\square$

3.3. Recovering the $\dot{B}^{\frac{d}{2}+1}$-estimates

As explained in Remark 3.3, we need to establish the uniform $L_t^\infty(\dot{B}^{\frac{d}{2}+1})$-norm estimate of $(w, r, u)$ which in fact leads to the uniform control of $\tilde{L}_t^2(\dot{B}^{\frac{d}{2}+1})$-norms for $(\frac{1}{\sqrt{\epsilon}}w, \frac{1}{\sqrt{\tau}}u)$ at both low and high frequencies as a byproduct.
**Lemma 3.3.** Let $T > 0$, and the threshold $J_{\tau}$ be given by (2.1). Then for any $t \in (0, T)$, the solution $(w, r, u)$ to the linear problem (3.1)$_2$–(3.1)$_4$ satisfies

$$ \begin{aligned} & \| (w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}_2^{d/2+1})} + \frac{1}{\sqrt{\varepsilon}} \| w \|_{\tilde{L}_t^2 (\dot{B}_2^{d/2+1})} + \frac{1}{\sqrt{\tau}} \| u \|_{\tilde{L}_t^2 (\dot{B}_2^{d/2+1})} \\ & \lesssim \| (w_0, r_0, u_0) \|_{\dot{B}_2^{d/2+1}} + \| (S_1, S_2, S_3) \|_{L_t^1 (\dot{B}_2^{d/2+1})} \\ & \qquad + (\eta + \sqrt{\mathcal{Z}(t)}) \mathcal{X}(t) + \frac{1}{\eta} \int_0^t \mathcal{V}(s) \mathcal{X}(s) \, ds, \end{aligned} \quad (3.45) $$

where $\eta \in (0, 1)$ is a constant to be chosen.

**Proof.** We perform a $L^2$-in-time type estimates and make use of the decay-in-$\tau$ of $u$ for $L^2$-time type norms. In fact, for any $j \in \mathbb{Z}$, by combining inequalities (3.34) and (3.35) together, we get

$$ \begin{aligned} & \frac{d}{dt} \int_{\mathbb{R}^d} \frac{1}{2} \left( \frac{h_5 + H_5}{h_1 + H_1} |w_j|^2 + \frac{h_4 + H_4}{h_3 + H_3} |r_j|^2 + |u_j|^2 \right) dx + \int_{\mathbb{R}^d} \left( \frac{1}{\tau} |u_j|^2 + \frac{(h_2 + H_2)(h_5 + H_5)}{\varepsilon(h_1 + H_1)} |w_j|^2 \right) dx \\ & \lesssim \| \operatorname{div} v \|_{L^\infty} \| (w_j, r_j, u_j) \|_{L^2}^2 + \left( \sum_{i=1}^5 \| \partial_i H_i \|_{L^\infty} + \| v \|_{L^\infty} \right) \| (w_j, r_j) \|_{L^2}^2 \\ & \qquad + \| \nabla H_4 \|_{L^\infty} \| (w_j, r_j) \|_{L^2} \| u_j \|_{L^2} + \| T_j^1 \|_{L^2} \| w_j \|_{L^2} \\ & \qquad + \| T_j^2 \|_{L^2} \| r_j \|_{L^2} + \| T_j^3 \|_{L^2} \| u_j \|_{L^2} + \| \dot{\Delta}_j (S_1, S_2, S_3) \|_{L^2} \| (w_j, r_j, u_j) \|_{L^2}. \end{aligned} \quad (3.46) $$

Furthermore, from (3.46) we have

$$ \begin{aligned} & \| (w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}_2^{d/2+1})} + \frac{1}{\sqrt{\varepsilon}} \| w \|_{\tilde{L}_t^2 (\dot{B}_2^{d/2+1})} + \frac{1}{\sqrt{\tau}} \| u \|_{\tilde{L}_t^2 (\dot{B}_2^{d/2+1})} \\ & \lesssim \| (w_0, r_0, u_0) \|_{\dot{B}_2^{d/2+1}} + \| v \|_{L_t^1 (\dot{B}_2^{d/2+1})} \| (w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}_2^{d/2+1})} \\ & \qquad + \left( \int_0^t \left( \sum_{i=1}^5 \| \partial_i H_i(s) \|_{\dot{B}_2^{d/2}} + \| v(s) \|_{\dot{B}_2^{d/2}} \right) \| (w, r)(s) \|_{\dot{B}_2^{d/2+1}} ds \right)^{\frac{1}{2}} \| (w, r) \|_{\tilde{L}_t^\infty (\dot{B}_2^{d/2+1})} \\ & \qquad + \left( \| H_4 \|_{\tilde{L}_t^\infty (\dot{B}_2^{d/2+1})} \frac{1}{\sqrt{\tau}} \| u \|_{\tilde{L}_t^2 (\dot{B}_2^{d/2+1})} \sqrt{\tau} \| (w, r) \|_{\tilde{L}_t^2 (\dot{B}_2^{d/2+1})} \right)^{\frac{1}{2}} \\ & \qquad + \sum_{j \in \mathbb{Z}} 2^{j(\frac{d}{2}+1)} \left( \int_0^t \left( \| T_j^1 \|_{L^2} \| w_j \|_{L^2} + \| T_j^2 \|_{L^2} \| r_j \|_{L^2} + \| T_j^3 \|_{L^2} \| u_j \|_{L^2} \right) ds \right)^{\frac{1}{2}} \\ & \qquad + \left( \| (S_1, S_2, S_3) \|_{L_t^1 (\dot{B}_2^{d/2+1})} \| (w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}_2^{d/2+1})} \right)^{\frac{1}{2}}. \end{aligned} \quad (3.47) $$
The right-hand side of inequality (3.47) can be estimated as follows. By the commutator estimate (2.5), we have

$$
\begin{aligned}
& \sum_{j \in \mathbb{Z}} 2^{j(\frac{d}{2}+1)} \left( \int_0^t \|T_j^1\|_{L^2} \|w_j\|_{L^2} ds \right)^{\frac{1}{2}} \\
& \lesssim \left( \|w\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \int_0^t \|v(s)\|_{\dot{B}_2^{\frac{d}{2}+1}} \|w(s)\|_{\dot{B}_2^{\frac{d}{2}+1}} ds \right. \\
& \qquad \left. + \|H_1\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \|w\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \|u\|_{L_t^1} (\dot{B}_2^{\frac{d}{2}+1}) + \|H_2\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \frac{1}{\varepsilon} \|w\|_{\tilde{L}_t^2}^2 (\dot{B}_2^{\frac{d}{2}+1}) \right)^{\frac{1}{2}} \\
& \lesssim \sqrt{\mathcal{Z}(t)} \mathcal{X}(t) + \left( \int_0^t \mathcal{V}(s) \mathcal{X}(s) ds \right)^{\frac{1}{2}} \sqrt{\mathcal{X}(t)}.
\end{aligned}
$$

Similarly, it holds

$$
\begin{aligned}
& \sum_{j \in \mathbb{Z}} 2^{j(\frac{d}{2}+1)} \left( \int_0^t \|T_j^2\|_{L^2} \|r_j\|_{L^2} ds \right)^{\frac{1}{2}} \\
& \lesssim \left( \|r\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \int_0^t \|v(s)\|_{\dot{B}_2^{\frac{d}{2}+1}} \|r(s)\|_{\dot{B}_2^{\frac{d}{2}+1}} ds + \|H_3\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \|r\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \|u\|_{L_t^1} (\dot{B}_2^{\frac{d}{2}+1}) \right)^{\frac{1}{2}} \\
& \lesssim \sqrt{\mathcal{Z}(t)} \mathcal{X}(t) + \left( \int_0^t \mathcal{V}(s) \mathcal{X}(s) ds \right)^{\frac{1}{2}} \sqrt{\mathcal{X}(t)},
\end{aligned}
$$

and

$$
\begin{aligned}
& \sum_{j \in \mathbb{Z}} 2^{j(\frac{d}{2}+1)} \left( \int_0^t \|T_j^3\|_{L^2} \|u_j\|_{L^2} ds \right)^{\frac{1}{2}} \\
& \lesssim \left( \|u\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \int_0^t \|v(s)\|_{\dot{B}_2^{\frac{d}{2}+1}} \|u(s)\|_{\dot{B}_2^{\frac{d}{2}+1}} ds + \|(H_4, H_5)\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \|r\|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) \|u\|_{L_t^1} (\dot{B}_2^{\frac{d}{2}+1}) \right)^{\frac{1}{2}} \\
& \lesssim \sqrt{\mathcal{Z}(t)} \mathcal{X}(t) + \left( \int_0^t \mathcal{V}(s) \mathcal{X}(s) ds \right)^{\frac{1}{2}} \sqrt{\mathcal{X}(t)}.
\end{aligned}
$$

We conclude from the above estimates that

$$
\begin{aligned}
& \| (w, r, u) \|_{\tilde{L}_t^\infty} (\dot{B}_2^{\frac{d}{2}+1}) + \frac{1}{\sqrt{\varepsilon}} \|w\|_{\tilde{L}_t^2} (\dot{B}_2^{\frac{d}{2}+1}) + \frac{1}{\sqrt{\tau}} \|u\|_{\tilde{L}_t^2} (\dot{B}_2^{\frac{d}{2}+1}) \\
& \lesssim \| (w_0, r_0, u_0) \|_{\dot{B}_2^{\frac{d}{2}+1}} + \| (S_1, S_2, S_3) \|_{L_t^1} (\dot{B}_2^{\frac{d}{2}+1}) \\
& \qquad + \sqrt{\mathcal{Z}(t)} \mathcal{X}(t) + \left( \int_0^t \mathcal{V}(s) \mathcal{X}(s) ds \right)^{\frac{1}{2}} \sqrt{\mathcal{X}(t)}.
\end{aligned}
$$

Applying Hölder's inequality to the above estimate leads to inequality (3.45). □
# 4. Global well-posedness for the nonlinear problems

## 4.1. The Cauchy problem of the system (BN)

In this section, we prove the uniform in $\varepsilon$ and $\tau$ global existence and uniqueness of solutions to the Cauchy problem for (BN) subject to the initial data $(\alpha_{\pm,0}, \rho_{\pm,0}, u_0)$. i.e. Theorem 1.1. For simplicity, we omit the superscript concerning the parameters $\varepsilon$ and $\tau$ in this section.

**Proof of Theorem 1.1.** Let $(\alpha_{\pm,0}, \rho_{\pm,0}, u_0)$ satisfy the smallness condition (1.7) and denote

$$ \mathcal{X}_0 := \| (\alpha_{\pm,0} - \bar{\alpha}_{\pm}, \rho_{\pm,0} - \bar{\rho}_{\pm}, u_0) \|_{\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d+1}} $$

Let $(y_0, w_0, r_0)$ be the perturbation of $(\alpha_{\pm,0}, \rho_{\pm,0})$ given by (1.21).

### • Step 1: Construction of approximation sequence

For any $n \ge 1$, we define the regularized perturbation

$$ (y_0^n, w_0^n, r_0^n, u_0^n)(x) = (\dot{S}_n y_0, \dot{S}_n w_0, \dot{S}_n r_0, \dot{S}_n u_0)(x), $$

where $\dot{S}_n$ is the low-frequency cut-off operator (see section 2). One can verify that $(y_0^n, w_0^n, r_0^n, u_0^n)$ is smooth and converges to $(y_0, w_0, r_0, u_0)$ strongly in $\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d+1}$ as $n \to \infty$. In addition, due to Lemmas 2.3 and 2.5, there exists a constant $C_0^*$ independent of $n$, $\varepsilon$ and $\tau$ such that

$$ \| (y_0^n, w_0^n, r_0^n, u_0^n) \|_{\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d+1}} \le C_0^* \mathcal{X}_0. \quad (4.1) $$

Set $(y^0, w^0, r^0, u^0) := (0, 0, 0, 0)$. For any $n \ge 0$, we consider the approximate scheme for (1.22) as follows:

$$ \left\{
\begin{aligned}
& \partial_t y^{n+1} + u^n \cdot \nabla y^{n+1} = 0, \\
& \partial_t w^{n+1} + u^n \cdot \nabla w^{n+1} + (\bar{F}_1 + G_1^n) \operatorname{div} u^{n+1} + (\bar{F}_2 + G_2^n) \frac{w^{n+1}}{\varepsilon} = 0, \\
& \partial_t r^{n+1} + u^n \cdot \nabla r^{n+1} + (\bar{F}_3 + G_3^n) \operatorname{div} u^{n+1} = F_4^n \frac{(w^n)^2}{\varepsilon}, \\
& \partial_t u^{n+1} + u^n \cdot \nabla u^{n+1} + \frac{u^{n+1}}{\tau} + (\bar{F}_0 + G_0^n) (\nabla r^{n+1} + (\gamma_+ - \gamma_-) \nabla w^{n+1}) = 0, \\
& (y^{n+1}, w^{n+1}, r^{n+1}, u^{n+1})(0,x) = (y_0^n, w_0^n, r_0^n, u_0^n)(x),
\end{aligned}
\right.
\quad (4.2) $$

with $F_i^n = F_i^{\varepsilon, \tau}(y^n, w^n, r^n)$, $\bar{F}_i$ and $G_i^n = G_i^{\varepsilon, \tau}(y^n, w^n, r^n)$ defined in (1.23)–(1.25), respectively. We define the functional space $\mathbb{E}_t$ associated to the following norm:

$$ \begin{aligned}
\| (y, w, r, u) \|_{\mathbb{E}_t} := \| (y, w, r, u) \|_{\tilde{L}_t^\infty (\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d+1})} + \| (\partial_t y, \partial_t w, \partial_t r, \partial_t u) \|_{L_t^1 (\dot{B}_{2/2}^d)} \\
+ \frac{1}{\varepsilon} \|w\|_{L_t^1 (\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^d)} + \frac{1}{\sqrt{\varepsilon}} \|w\|_{\tilde{L}_t^2 (\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^d)} \\
+ \tau \|r\|_{L_t^1 (\dot{B}_{2/2}^{d+1} \cap \dot{B}_{2/2}^d)^2} + \|r\|_{L_t^1 (\dot{B}_{2/2}^{d+1})}^h + \tau \|r\|_{L_t^1 (\dot{B}_{2/2}^{d+1})} + \sqrt{\tau} \|r\|_{\tilde{L}_t^2 (\dot{B}_{2/2}^d \cap \dot{B}_{2/2}^{d+1})} \\
+ \|u\|_{L_t^1 (\dot{B}_{2/2}^d \cap \dot{B}_{2/2}^d)^+} + \frac{1}{\sqrt{\tau}} \|u\|_{\tilde{L}_t^2 (\dot{B}_{2/2}^{d-1} \cap \dot{B}_{2/2}^{d+1})}.
\end{aligned} $$
For any fixed $n \ge 1$, we assume that $(y^n, w^n, r^n, u^n)$ satisfies

$$ \| (y^n, w^n, r^n, u^n) \|_{E_t} + \frac{1}{\tau} \| u^n + \tau (\bar{F}_0 + G_0^{n-1}) \nabla r^n \|_{L_t^1 (\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^d)} \le 2C_0 C_0^* \chi_0, \quad t > 0, \quad (4.3) $$

where the constants $C_0$ and $C_0^*$ are given by (3.3) and (4.1), respectively. Since the initial data is smooth, by virtue of the classical theorems for the transport equation (4.2)$_1$ and the symmetric hyperbolic system (4.2)$_2$–(4.2)$_4$ (see [2, 4]), there exists a unique global solution $(y^{n+1}, w^{n+1}, r^{n+1}, u^{n+1}) \in C(\mathbb{R}_+; H^s)$ with all $s > \frac{d}{2} + 1$.

### Step 2: Uniform estimate

Our goal is to show that $(y^{n+1}, w^{n+1}, r^{n+1}, u^{n+1})$ also satisfies the estimate (4.3) uniformly in $n, \varepsilon, \tau$ and time. To this end, we first let $\chi_0 \le 1$. It follows from (4.3) and the composition estimates in Lemma 2.5 that

$$ \sum_{i=0}^{4} \|G_i^n\|_{L_t^\infty(\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{d+1})} \le C_2^* \|(y^n, r^n, w^n)\|_{L_t^\infty(\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{d+1})}, \quad (4.4) $$

and similarly,

$$ \sum_{i=0}^{4} \| \partial_t G_i^n \|_{L_t^1 (\dot{B}_{\frac{d}{2}}^d)} \le C_3^* \|(\partial_t y^n, \partial_t w^n, \partial_t r^n)\|_{L_t^1 (\dot{B}_{\frac{d}{2}}^d)}, \quad (4.5) $$

for some universal constants $C_2^*$ and $C_3^*$. According to (4.1) and (4.4), we can justify the condition (3.2) provided that

$$ \chi_0 \le c_1^* := \frac{c}{2C_0 C_0^* C_2^*}. $$

Hence, we are able to employ the uniform a-priori estimate established in proposition 3.1 to obtain

$$ \begin{aligned} & \| (y^{n+1}, w^{n+1}, r^{n+1}, u^{n+1}) \|_{E_t} + \frac{1}{\tau} \| u^{n+1} + \tau (\bar{F}_0 + G_0^n) \nabla r^{n+1} \|_{L_t^1 (\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^d)} \\ & \le C_0 \exp \left( C_0 \int_0^t \left( \|u^n(s)\|_{\dot{B}_{\frac{d}{2}}^d \cap \dot{B}_{\frac{d}{2}}^{d+1}} + \sum_{i=0}^4 \| \partial_t G_i^n(s) \|_{\dot{B}_{\frac{d}{2}}^d} \right) ds \right) \quad (4.6) \\ & \times \left( \| (y_0^n, w_0^n, r_0^n, u_0^n) \|_{\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{d+1}} + \left\| F_4^n \frac{(w^n)^2}{\varepsilon} \right\|_{L_t^1 (\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{d+1})} \right). \end{aligned} $$

By (4.3) and lemmas 2.3 and 2.5, it holds that

$$ \left\| F_4^n \frac{(w^n)^2}{\varepsilon} \right\|_{L_t^1 (\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{d+1})} \le \frac{C_4^*}{\varepsilon} \|w^n\|_{L_t^2 (\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{d+1})}^2, \quad (4.7) $$
where $C_4^* > 0$ is a universal constant. Combining (4.1), (4.3), (4.5), (4.6) and (4.7) together, we have

$$ \begin{aligned} & \| (y^{n+1}, w^{n+1}, r^{n+1}, u^{n+1}) \|_{E_t} + \frac{1}{\tau} \|u^{n+1} + \tau (\bar{F}_0 + G_0^n) \nabla r^{n+1} \|_{L_t^1 (\dot{B}_2^{d-1} \cap \dot{B}_2^d)} \\ & \le C_0 e^{2(1+C_3^*)} C_0^2 C_0^* \chi_0 (C_0^* \chi_0 + C_4^* (2C_0 C_0^* \chi_0)^2) \le 2C_0 C_0^* \chi_0, \end{aligned} $$

as long as

$$ \chi_0 \le c_2^* := \min \left\{ \frac{1}{2(1+C_3^*) C_0^2 C_0^* \log^{\frac{3}{2}}}, \frac{2}{9(2C_0 C_0^*)^2 C_4^*} \right\} $$

such that $e^{2(1+C_3^*)} C_0^2 C_0^* \chi_0 \le \frac{3}{2}$ and $C_4^* (2C_0 C_0^* \chi_0)^2 \le \frac{1}{3} \chi_0$. Thus, the uniform estimate (4.3) holds for any $n \ge 0$.

### Step 3: Strong convergence

The uniform estimate (4.3) enables us to obtain the weak compactness of the approximate sequence. In order to pass the limit in every nonlinear term of (4.2) as $n \to \infty$, one needs to have robust strong compactness in a suitable sense. Classical compact embedding theorem merely gives the strong convergence locally in space-time and up to a subsequence, which is not enough for the system (4.2). In what follows, we show that the strong convergence holds in $\mathbb{R}_+ \times \mathbb{R}^d$ for the whole sequence. □

**Lemma 4.1.** There exists a small constant $c_3^* \in (0, \min\{1, c_1^*, c_2^*\}]$ and a limit $(y, w, r, u)$ such that if $\chi_0 \le c_3^*$, then as $n \to \infty$,

$$ (y^n, w^n, r^n, u^n) \to (y, w, r, u) \quad \text{strongly in } L^\infty(\mathbb{R}_+; \dot{B}_2^{d-1} \cap \dot{B}_2^d). \quad (4.8) $$

In particular, we have

$$ (y^n, w^n, r^n, u^n) \to (y, w, r, u) \quad \text{strongly in } L^\infty(\mathbb{R}_+; L^d \cap L^\infty). \quad (4.9) $$

**Proof.** In order to show (4.8), one needs to perform uniform energy estimates on the error unknown

$$ (\tilde{y}^n, \tilde{w}^n, \tilde{r}^n, \tilde{u}^n) := (y^{n+1} - y^n, w^{n+1} - w^n, r^{n+1} - r^n, u^{n+1} - u^n). $$

Following the framework in section 3, we aim to estimate the functional

$$ \begin{aligned} \tilde{\chi}^n(t) = & \| (\tilde{y}^n, \tilde{w}^n, \tilde{r}^n, \tilde{u}^n) \|_{\tilde{L}_t^\infty (\dot{B}_2^{d-1} \cap \dot{B}_2^d)} + \frac{1}{\varepsilon} \|\tilde{w}^n\|_{L_t^1 (\dot{B}_2^{d-1})} + \frac{1}{\sqrt{\varepsilon}} \|\tilde{w}^n\|_{\tilde{L}_t^2 (\dot{B}_2^{d-1} \cap \dot{B}_2^d)} \\ & + \tau \|\tilde{r}^n\|_{L_t^1 (\dot{B}_2^{d+1})} + \|\tilde{r}^n\|_{L_t^1 (\dot{B}_2^d)} + \sqrt{\tau} \|\tilde{r}^n\|_{\tilde{L}_t^2 (\dot{B}_2^d)} \\ & + \|\tilde{u}^n\|_{L_t^1 (\dot{B}_2^d)} + \frac{1}{\sqrt{\tau}} \|\tilde{u}^n\|_{\tilde{L}_t^2 (\dot{B}_2^{d-1} \cap \dot{B}_2^d)} + \frac{1}{\tau} \|\tilde{u}^n + \tau (\bar{F}_0 + G_0^{n-1}) \nabla \tilde{r}^n\|_{L_t^1 (\dot{B}_2^{d-1})}. \end{aligned} $$
To this matter, one can verify that for any $n \ge 1$, $(\tilde{y}^n, \tilde{w}^{n+1}, \tilde{r}^{n+1}, \tilde{u}^{n+1})$ solves

$$ egin{cases} 
\partial_t \tilde{y}^{n+1} + u^n \cdot \nabla \tilde{y}^{n+1} = \tilde{S}_1^n, \\
\partial_t \tilde{w}^{n+1} + u^n \cdot \nabla \tilde{w}^{n+1} + (\bar{F}_1 + G_1^n) \operatorname{div} \tilde{u}^{n+1} + (\bar{F}_2 + G_2^n) \frac{\tilde{w}^{n+1}}{\varepsilon} = \tilde{S}_2^n, \\
\partial_t \tilde{r}^{n+1} + u^n \cdot \nabla \tilde{r}^{n+1} + (\bar{F}_3 + G_3^n) \operatorname{div} \tilde{u}^{n+1} = \tilde{S}_3^n, \\
\partial_t \tilde{u}^{n+1} + u^n \cdot \nabla \tilde{u}^{n+1} + \frac{\tilde{u}^{n+1}}{\tau} + (\bar{F}_0 + G_0^n) (\nabla \tilde{r}^{n+1} + (\gamma_+ - \gamma_-) \nabla \tilde{w}^{n+1}) = \tilde{S}_4^n, \\
(\tilde{y}^{n+1}, \tilde{w}^{n+1}, \tilde{r}^{n+1}, \tilde{u}^{n+1}) (0, x) = \dot{\Delta}_n (y_0, w_0, r_0, u_0) (x),
\end{cases} 
\quad (4.10) $$

with

$$ egin{cases} 
\tilde{S}_1^n := -\tilde{u}^n \cdot \nabla y^n, \\
\tilde{S}_2^n := -\tilde{u}^n \cdot \nabla w^n - (G_1^n - G_1^{n-1}) \operatorname{div} u^n - (G_2^n - G_2^{n-1}) \frac{w^n}{\varepsilon}, \\
\tilde{S}_3^n := -\tilde{u}^n \cdot \nabla r^n - (G_3^n - G_3^{n-1}) \operatorname{div} u^n + (F_4^n - F_4^{n-1}) \frac{(w^n)^2}{\varepsilon} + F_4^{n-1} \frac{(w^n + w^{n-1}) \tilde{w}^n}{\varepsilon}, \\
\tilde{S}_4^n := -\tilde{u}^n \cdot \nabla u^n - (G_0^n - G_0^{n-1}) (\nabla r^n + (\gamma_+ - \gamma_-) \nabla w^n).
\end{cases} $$

First, employing Lemma 2.7 to (4.10)$_1$ yields

$$ \begin{aligned} \|\tilde{y}^{n+1}\|_{\tilde{L}_t^\infty(\tilde{B}_2^{d-1} \cap \tilde{B}_2^d)} & \lesssim \exp \left( \|u^n\|_{L_t^1(\tilde{B}_2^{d+1})} \right) \left( \|\dot{\Delta}_n y_0\|_{\tilde{B}_2^{d-1} \cap \tilde{B}_2^d} + \|\tilde{S}_1^n\|_{L_t^1(\tilde{B}_2^{d-1} \cap \tilde{B}_2^d)} \right). \end{aligned} \quad (4.11) $$

By similar computations on (4.10)$_2$–(4.10)$_4$ as in Lemma 3.1, we have the low-frequency estimate at the $\tilde{B}_2^{d-1}$ regularity level:

$$ egin{aligned} & \| (\tilde{w}^{n+1}, \tilde{r}^{n+1}, \tilde{u}^{n+1}) \|_{\tilde{L}_t^\infty(\tilde{B}_2^{d-1})}^\ell + \tau \|\tilde{r}^{n+1}\|_{L_t^1(\tilde{B}_2^{d+1})}^\ell \\ & \quad + \sqrt{\tau} \|\tilde{r}^{n+1}\|_{\tilde{L}_t^2(\tilde{B}_2^d)}^\ell + \frac{1}{\varepsilon} \|\tilde{w}^{n+1}\|_{L_t^1(\tilde{B}_2^{d-1})}^\ell + \frac{1}{\sqrt{\varepsilon}} \|\tilde{w}^{n+1}\|_{\tilde{L}_t^2(\tilde{B}_2^{d-1})}^\ell \\ & \quad + \frac{1}{\tau} \|\tilde{u}^{n+1} + \tau (\bar{F}_0 + G_0^n) \nabla \tilde{r}^{n+1}\|_{L_t^1(\tilde{B}_2^{d-1})}^\ell \\ & \lesssim \|\dot{\Delta}_n (w_0, r_0, u_0)\|_{\tilde{B}_2^{d-1}} + \| (\tilde{S}_2^n, \tilde{S}_3^n, \tilde{S}_4^n) \|_{L_t^1(\tilde{B}_2^{d-1})}^\ell \\ & \quad + Z^n(t) \tilde{X}^{n+1}(t) + \int_0^t V^n(s) \tilde{X}^{n+1}(s) ds, \end{aligned} \quad (4.12) $$

with

$$ egin{cases} 
Z^n(t) = \sum_{i=0}^3 \|G_i^n\|_{\tilde{B}_2^{d-1} \cap \tilde{B}_2^{d+1}}, \\ 
V^n(t) = \|u^n\|_{\tilde{B}_2^d \cap \tilde{B}_2^{d+1}} + \sum_{i=0}^3 \|\partial_t G_i^n\|_{\tilde{B}_2^{d-1} \cap \tilde{B}_2^{d+1}}. 
\end{cases} $$
Moreover, as in Lemma 3.2, one can obtain the following estimate in the high-frequency region:

$$
\begin{align}
& \| (\tilde{w}^{n+1}, \tilde{r}^{n+1}, \tilde{u}^{n+1}) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1}) \cap L_t^1 (\dot{B}^{\frac{d}{2}})}^h \nonumber \\
& \quad + \frac{1}{\varepsilon} \|\tilde{w}^{n+1}\|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^h + \frac{1}{\sqrt{\varepsilon}} \|\tilde{w}^{n+1}\|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^h + \sqrt{\tau} \|\tilde{r}^{n+1}\|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}})}^h \nonumber \\
& \quad + \frac{1}{\tau} \|\tilde{u}^{n+1} + \tau (\bar{F}_0 + G_0^n) \nabla \tilde{r}^{n+1}\|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^h \tag{4.13} \\
& \lesssim \|\dot{\Delta}_n (w_0, r_0, u_0)\|_{\dot{B}^{\frac{d}{2}}}^h + \| (\tilde{S}_2^n, \tilde{S}_3^n, \tilde{S}_4^n) \|_{L_t^1 (\dot{B}^{\frac{d}{2}})}^h \nonumber \\
& \quad + \mathcal{Z}^n (t) \tilde{\chi}^{n+1} (t) + \int_0^t \mathcal{V}^n (s) \tilde{\chi}^{n+1} (s) \mathrm{d}s. \nonumber
\end{align}
$$

Finally, following the proof of Lemma 3.3, we also have

$$
\begin{align}
& \| (\tilde{w}^{n+1}, \tilde{r}^{n+1}, \tilde{u}^{n+1}) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}})} + \frac{1}{\sqrt{\varepsilon}} \|\tilde{w}^{n+1}\|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}})} + \frac{1}{\sqrt{\tau}} \|\tilde{u}^{n+1}\|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}})} \nonumber \\
& \lesssim \|\dot{\Delta}_n (w_0, r_0, u_0)\|_{\dot{B}^{\frac{d}{2}}} + \| (\tilde{S}_2^n, \tilde{S}_3^n, \tilde{S}_4^n) \|_{L_t^1 (\dot{B}^{\frac{d}{2}})} \tag{4.14} \\
& \quad + (\tilde{\eta} + \sqrt{\mathcal{Z}^n(t)}) \tilde{\chi}^{n+1}(t) + \frac{1}{\tilde{\eta}} \int_0^t \mathcal{V}^n(s) \tilde{\chi}^{n+1}(s) \mathrm{d}s, \nonumber
\end{align}
$$

where $\tilde{\eta} > 0$ is a constant to be chosen. Combining (4.11)–(4.14) together, we arrive at

$$
\begin{align}
\tilde{\chi}^{n+1}(t) &\lesssim \|\dot{\Delta}_n (y_0, w_0, r_0, u_0)\|_{\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}}} + \| (\tilde{S}_2^n, \tilde{S}_3^n, \tilde{S}_4^n) \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})} \tag{4.15} \\
&\quad + (\tilde{\eta} + \sqrt{\mathcal{Z}^n(t)} + \mathcal{Z}^n(t)) \tilde{\chi}^{n+1}(t) + \left(1 + \frac{1}{\tilde{\eta}}\right) \int_0^t \mathcal{V}^n(s) \tilde{\chi}^{n+1}(s) \mathrm{d}s. \nonumber
\end{align}
$$

Now we estimate the right-hand side of (4.15) as follow. First, due to $\Delta_{n'} \dot{\Delta}_n = 0$ with $|n - n'| \ge 2$, one has

$$
\| \dot{\Delta}_n (y_0, w_0, r_0, u_0) \|_{\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}}} \lesssim \sum_{|n'-n| \le 1} \left( 2^{\frac{d}{2}-1} n' + 2^{\frac{d}{2}n'} \right) \| \dot{\Delta}_{n'} (y_0, w_0, r_0, u_0) \|_{L^2}.
$$

Thence, applying uniform estimate (4.3) leads to

$$
\left\{
\begin{array}{l}
\mathcal{Z}^n(t) \lesssim \|(y^n, w^n, r^n)\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1})} \lesssim \chi_0, \\
\\
\int_0^t \mathcal{V}^n(s) \mathrm{d}s \lesssim \|(\partial_t y^n, \partial_t w^n, \partial_t r^n)\|_{L_t^1(\dot{B}^{\frac{d}{2}})} \lesssim \chi_0.
\end{array}
\right.
$$

Regarding the nonlinear terms on the right-hand side of (4.15), one deduces from (2.4), (2.6)–(2.8) and (4.3) that
$$
\begin{aligned}
\| (\tilde{S}_1^n, \tilde{S}_2^n, \tilde{S}_3^n, \tilde{S}_4^n) \|_{L_t^1(\hat{B}_2^{\frac{d}{2}-1} \cap \hat{B}_2^{\frac{d}{2}})} &\lesssim \| (y^n, w^n, r^n) \|_{\tilde{L}_t^\infty(\hat{B}_2^{\frac{d}{2}} \cap \hat{B}_2^{\frac{d}{2}+1})} \| \tilde{u}^n \|_{L_t^1(\hat{B}_2^{\frac{d}{2}})} \\
&\quad + \| (\tilde{y}^n, \tilde{w}^n, \tilde{r}^n, \tilde{u}^n) \|_{\tilde{L}_t^\infty(\hat{B}_2^{\frac{d}{2}})} \| u^n \|_{L_t^1(\hat{B}_2^{\frac{d}{2}} \cap \hat{B}_2^{\frac{d}{2}+1})} \\
&\quad + \frac{1}{\sqrt{\epsilon}} \| (w^n, w^{n-1}) \|_{\tilde{L}_t^2(\hat{B}_2^{\frac{d}{2}})} \frac{1}{\sqrt{\epsilon}} \| \tilde{w}^n \|_{\tilde{L}_t^2(\hat{B}_2^{\frac{d}{2}-1} \cap \hat{B}_2^{\frac{d}{2}})} \\
&\lesssim \chi_0 \tilde{X}^n(t).
\end{aligned}
$$

Gathering the above estimates into (4.15) and (4.11) and letting both $\tilde{\eta}$ and $\chi_0$ be sufficiently small, we obtain

$$ \tilde{\chi}^{n+1}(t) \lesssim \sum_{|n'-n| \le 1} \left(2^{\frac{d}{2}-1} n' + 2^{\frac{d}{2}n'}\right) \|\dot{\Delta}_{n'}(y_0, w_0, r_0, u_0)\|_{L^2} + \chi_0 \tilde{\chi}^n(t). $$

Summing this over $n \ge 1$ leads to

$$ \sum_{n=1}^{\infty} \tilde{\chi}^{n+1}(t) \lesssim \|(y_0, w_0, r_0, u_0)\|_{\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{\frac{d}{2}}} + \chi_0 \sum_{n=1}^{\infty} \tilde{\chi}^n(t). $$

Given $(y^0, w^0, r^0, u^0) = (0, 0, 0, 0)$ and $(y^1, w^1, r^1, u^1) = (y_0, w_0, r_0, u_0)$, we take sufficiently small $\chi_0$ to have

$$ \sum_{n=0}^{\infty} \tilde{\chi}^n(t) \lesssim \|(y_0, w_0, r_0, u_0)\|_{\dot{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \dot{B}_{\frac{d}{2}}^{\frac{d}{2}}} \quad (4.16) $$

Now we define $(y, w, r, u) := \sum_{n'=0}^{\infty} (\tilde{y}^{n'}, \tilde{w}^{n'}, \tilde{r}^{n'}, \tilde{u}^{n'})$. Thanks to (4.16) and

$$ (y^n, w^n, r^n, u^n) = \sum_{n'=0}^{n} (\tilde{y}^{n'}, \tilde{w}^{n'}, \tilde{r}^{n'}, \tilde{u}^{n'}), $$

it follows that

$$ \| (y^n, w^n, r^n, u^n) - (y, w, r, u) \|_{L_t^\infty(\hat{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \hat{B}_{\frac{d}{2}}^{\frac{d}{2}})} \lesssim \sum_{n' \ge n+1}^{\infty} \tilde{\chi}^{n'}(t) \xrightarrow{n \to \infty} 0. $$

Gathering the embeddings in Lemma 2.2, we get (4.9), which finishes the proof of Lemma 4.1. $\square$

### Step 4: Global existence

Let $\chi_0 \le c_3^*$. By virtue of the strong convergence properties (4.8)–(4.9), one can pass to the limit as $n \to \infty$ in (4.2) and justify that the limit $(y, w, r, u)$, obtained in Lemma 4.1, is indeed a global strong solution to the Cauchy problem (1.22). In addition, taking advantage of Fatou’s property in Lemma 2.2, for all $t > 0$, we have

$$ \begin{aligned} \| (y, w, r, u) \|_{\mathbb{E}_t} + \frac{1}{\tau} \| u + \tau (\bar{F}_0 + G_0) \nabla r \|_{L_t^1(\hat{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \hat{B}_{\frac{d}{2}}^{\frac{d}{2}})} & \lesssim \liminf_{n \to \infty} \left( \| (y^n, w^n, r^n, u^n) \|_{\mathbb{E}_t} + \frac{1}{\tau} \| u^n + \tau (\bar{F}_0 + G_0^{n-1}) \nabla r^n \|_{L_t^1(\hat{B}_{\frac{d}{2}-1}^{\frac{d}{2}} \cap \hat{B}_{\frac{d}{2}}^{\frac{d}{2}})} \right) \lesssim \chi_0. \end{aligned} \quad (4.17) $$
To prove the time continuity property in (1.8) and (1.9), our proof relies on the uniform bound (4.17) and employs a reasoning analogous to that found in [22]. Since $\|(\partial_t y, \partial_t w, \partial_t r, \partial_t u)\|$ lies in $L^1(\mathbb{R}_+; \dot{B}^{\frac{d}{2}})$, one has $(y, w, r, u) \in C_b(\mathbb{R}_+; \dot{B}^{\frac{d}{2}})$. To recover $(y, w, r, u) \in C_b(\mathbb{R}_+; \dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1})$, we shall investigate each equations separately. Recall that the solution $(y, w, r, u)$ satisfies

$$ egin{cases} 
\partial_t y = -u \cdot \nabla y, \\
\partial_t w = -u \cdot \nabla w - (\bar{F}_1 + G_1) \operatorname{div} u - (\bar{F}_2 + G_2) \frac{w}{\varepsilon}, \\
\partial_t r = -u \cdot \nabla r - (\bar{F}_3 + G_3) \operatorname{div} u - F_4 \frac{w^2}{\varepsilon}, \\
\partial_t u = -u \cdot \nabla u - \frac{u}{\tau} - (\bar{F}_0 + G_0) \nabla r - (\gamma_+ - \gamma_-) (\bar{F}_0 + G_0) \nabla w.
\end{cases} $$

As the right-hand side terms of the components $y, r$ and $w$ belong to $L^1(\mathbb{R}_+; \dot{B}^{\frac{d}{2}-1})$, we directly get $(y, r, w) \in C_b(\mathbb{R}_+; \dot{B}^{\frac{d}{2}-1})$. Concerning the equation of $u$, its right-hand side lies in $L^2(\mathbb{R}_+; \dot{B}^{\frac{d}{2}-1})$ thus we can also recover that $u$ belongs to $C_b(\mathbb{R}_+; \dot{B}^{\frac{d}{2}-1})$. We are left with recovering the time continuity of $(y, w, r, u)$ in $\dot{B}^{\frac{d}{2}+1}$. First, for a fixed $j \in \mathbb{Z}$, each $(y_j, w_j, r_j, u_j)$ is continuous in time with values in $L^2$ due to Bernstein's inequality. Now, thanks to $(y, w, r, u) \in L^\infty(\mathbb{R}_+; \dot{B}^{\frac{d}{2}+1})$, for any $\eta > 0$, there exists an large integer $J_*$ such that, for all $t > 0$

$$ \sum_{|j| \ge J_*} 2^{j(\frac{d}{2}+1)} \|(y_j, w_j, r_j, u_j)\|_{L^\infty(L^2)} < \frac{\eta}{2}. $$

Thus, for any time $t, t' \in \mathbb{R}_+$, we have

$$ \begin{aligned} \|(y, w, r, u)(t) - (y, w, r, u)(t')\|_{\dot{B}^{\frac{d}{2}+1}} & \le \sum_{|j|<J_*} 2^{j(\frac{d}{2}+1)} \|\dot{\Delta}_j((y, w, r, u)(t) - (y, w, r, u)(t'))\|_{L^2} \\ & \quad + \sum_{|j| \ge J_*} 2^{j(\frac{d}{2}+1)} \|\dot{\Delta}_j((y, w, r, u)(t) - (y, w, r, u)(t'))\|_{L^2} \\ & \le 2^{J_*} \|(y, w, r, u)(t) - (y, w, r, u)(t')\|_{\dot{B}^{\frac{d}{2}} + \eta} \to_{t \to t'} \eta. \end{aligned} $$

Since $\eta$ is an arbitrary constant, we get $(y, w, r, u) \in C_b(\mathbb{R}_+; \dot{B}^{\frac{d}{2}+1})$. Finally, applying the inverse function theorem, we can see that once $\alpha_\pm$ and $\rho_\pm$ are determined by (1.20), then $(\alpha_\pm, \rho_\pm, u) \in C_b(\mathbb{R}_+; \dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1})$ is the unique global strong solution to the original Cauchy problem of the system (BN). Using (4.17), product laws and composite estimates, we are able to verify that $(\alpha_\pm, \rho_\pm, u)$ satisfies the properties (1.8) and (1.9). To complete the proof, we prove the uniqueness in our regularity framework below.

### Step 5: Uniqueness

The final step is to show the uniqueness of solutions to (BN) belonging to the regularity class (1.8). We emphasize that the proof does not require the smallness of regularity for initial data. It suffices to consider the reformulated system (1.22). Without loss of generality, as the parameters do not affect our argument for proving uniqueness, we set $\varepsilon = \tau = 1$. Let $(y_1, w_1, r_1, u_1)$ and $(y_2, w_2, r_2, u_2)$ be two solutions to (1.22) subject to the same initial data $(y_0, w_0, r_0, u_0)$, satisfying (1.8) and $\bar{F}_i + G_i(y_j, r_j, w_j) > 0$ for $i = 0, 1, 2, 3$ and $j = 1, 2$. The difference

$$ (\tilde{y}, \tilde{w}, \tilde{r}, \tilde{u}) := (y_1 - y_2, w_1 - w_2, r_1 - r_2, u_1 - u_2) $$
solves

$$
\left\{
\begin{aligned}
\partial_t \tilde{y} + u_1 \cdot \nabla \tilde{y} &= \tilde{S}_1, \\
\partial_t \tilde{w} + u_1 \cdot \nabla \tilde{w} + (\bar{F}_1 + G_1^1) \operatorname{div} \tilde{u} + \bar{F}_2 \tilde{w} &= \tilde{S}_2, \\
\partial_t \tilde{r} + u_1 \cdot \nabla \tilde{r} + (\bar{F}_3 + G_3^1) \operatorname{div} \tilde{u} &= \tilde{S}_3, \\
\partial_t \tilde{u} + u_1 \cdot \nabla \tilde{u} + \tilde{u} + (\bar{F}_0 + G_0^1) (\nabla \tilde{r} + (\gamma_+ - \gamma_-) \nabla \tilde{w}) &= \tilde{S}_4, \\
(\tilde{y}, \tilde{w}, \tilde{r}, \tilde{u}) (0, x) &= (0, 0, 0, 0),
\end{aligned}
\right.
\qquad (4.18)
$$

where we denoted

$$
\left\{
\begin{aligned}
\tilde{S}_1 &:= -\tilde{u} \cdot \nabla y_2, \\
\tilde{S}_2 &:= -\tilde{u} \cdot \nabla w_2 - (G_1^1 - G_1^2) \operatorname{div} u_2 - (G_2^1 - G_2^2) w_2 - G_2^1 \tilde{w}, \\
\tilde{S}_3 &:= -\tilde{u} \cdot \nabla r_2 - (G_3^1 - G_3^2) \operatorname{div} u_2 + (F_4^1 - F_4^2) (w_1)^2 - F_4^2 (w_1 + w_2) \tilde{w}, \\
\tilde{S}^4 &:= -\tilde{u} \cdot \nabla u_2 - (G_0^1 - G_0^2) (\nabla r_2 + (\gamma_+ - \gamma_-) \nabla w_2).
\end{aligned}
\right.
$$

Here $G_i^l := G_i(y_l, w_l, r_l)$ and $F_4^l := F_4(y_l, w_l, r_l)$ with $i=0,1,2,3$, and $l=1,2$. Applying Lemma 2.7 to (4.18) implies that, for all $t > 0$,

$$
\|\tilde{y}(t)\|_{\dot{B}^{\frac{d}{2}}} \lesssim \exp \left( \|u_1\|_{L^1(\dot{B}^{\frac{d}{2}+1})} \right) \int_0^t \|\tilde{S}_1^n\|_{\dot{B}^{\frac{d}{2}}} \mathrm{d}s. \quad (4.19)
$$

Through the application of the weighted Lyapunov functional method, as outlined in (3.32)–
(3.36), we obtain

$$
\begin{align}
& \frac{d}{dt} \int_{\mathbb{R}^d} \frac{1}{2} (W_1 |\tilde{w}_j|^2 + W_2 |\tilde{r}_j|^2 + |\tilde{u}_j|^2) \, dx + \int_{\mathbb{R}^d} (W_3 |\tilde{w}_j|^2 + |\tilde{u}_j|^2) \, dx \nonumber \\
& \lesssim (\| (\partial_t y_1, \partial_t w_1, \partial_t r_1, \nabla y_1, \nabla w_1, \nabla r_1) \|_{L^\infty} + \| u_1 \|_{W^{1,\infty}}) \| (\tilde{w}_j, \tilde{r}_j, \tilde{u}_j) \|_{L^2}^2 \tag{4.20} \\
& \quad + (\| (\tilde{T}_{1,j}, \tilde{T}_{2,j}, \tilde{T}_{3,j}) \|_{L^2} + \| \dot{\Delta}_j (\tilde{S}_2, \tilde{S}_3, \tilde{S}_4) \|_{L^2}) \| (\tilde{w}_j, \tilde{r}_j, \tilde{u}_j) \|_{L^2}, \nonumber
\end{align}
$$

where $W_i = W_i(y_1, w_1, r_1) > 0$, $i = 1, 2, 3$ are some smooth weight functions depending on
$\bar{F}_i + G_i^1 > 0$, and the commutator terms are given by

$$
\left\{
\begin{aligned}
\tilde{T}_{1,j} &:= [u_1, \dot{\Delta}_j] \nabla \tilde{w} + [G_1^1, \dot{\Delta}_j] \operatorname{div} \tilde{u}, \\
\tilde{T}_{2,j} &:= [u_1, \dot{\Delta}_j] \nabla \tilde{r} + [G_3^1, \dot{\Delta}_j] \operatorname{div} \tilde{u}, \\
\tilde{T}_{3,j} &:= [u_1, \dot{\Delta}_j] \nabla \tilde{u} + [G_0^1, \dot{\Delta}_j] \nabla \tilde{r} + [G_0^1, \dot{\Delta}_j] \nabla \tilde{w}.
\end{aligned}
\right.
$$

Integrating (4.20) over $[0, t]$, taking the square root of both sides and summing the resulting estimate over $j \in \mathbb{Z}$ with the weight $2^{\frac{d}{2}j}$, we get

$$
\begin{align*}
& \| (\tilde{w}, \tilde{r}, \tilde{u}) \|_{\tilde{L}_t^\infty (\tilde{B}^{\frac{d}{2}})} + \| \tilde{u} \|_{\tilde{L}_t^2 (\tilde{B}^{\frac{d}{2}})} + \| \tilde{w} \|_{\tilde{L}_t^2 (\tilde{B}^{\frac{d}{2}})} \\
& \lesssim \left( \int_0^t \left( \| (\partial_t y_1, \partial_t w_1, \partial_t r_1, \nabla y_1, \nabla w_1, \nabla r_1) \|_{L^\infty} + \| u_1 \|_{W^{1,\infty}} \right) \| (\tilde{w}, \tilde{r}, \tilde{u}) \|_{\tilde{B}^{\frac{d}{2}}} \right. \\
& \qquad \left. + \sum_{j \in \mathbb{Z}} 2^{\frac{j d}{2}} \| (\tilde{T}_{1,j}, \tilde{T}_{2,j}, \tilde{T}_{3,j}) \|_{L^2} + \| (\tilde{S}_2, \tilde{S}_3, \tilde{S}_4) \|_{\dot{B}^{\frac{d}{2}}} \right) \frac{1}{ds} \| (\tilde{w}, \tilde{r}, \tilde{u}) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}})}
\end{align*}
$$
which together with Young's inequality implies

$$ \begin{aligned} & \|(\tilde{w}, \tilde{r}, \tilde{u})(t)\|_{\dot{B}^{\frac{d}{2}}} \\ & \lesssim \int_{0}^{t} \left( \| (\partial_t y_1, \partial_t w_1, \partial_t r_1, \nabla y_1, \nabla w_1, \nabla r_1) \|_{L^{\infty}} + \|u_1\|_{W^{1,\infty}} \| (\tilde{w}, \tilde{r}, \tilde{u}) \|_{\dot{B}^{\frac{d}{2}}} \right. \\ & \left. + \sum_{j \in \mathbb{Z}} 2^{j \frac{d}{2}} \| (\tilde{T}_{1,j}, \tilde{T}_{2,j}, \tilde{T}_{3,j}) \|_{L^2} + \| (\tilde{S}_2, \tilde{S}_3, \tilde{S}_4) \|_{\dot{B}^{\frac{d}{2}}} \right) ds, \quad t > 0. \end{aligned} \quad (4.21) $$

Using the classical commutator estimate in Lemma 2.4 implies

$$ \sum_{j \in \mathbb{Z}} 2^{j \frac{d}{2}} \| (\tilde{T}_{1,j}, \tilde{T}_{2,j}, \tilde{T}_{3,j}) \|_{L^2} \lesssim \| (y_1, r_1, w_1, u_1) \|_{\dot{B}^{\frac{d}{2}+1}} \| (\tilde{w}, \tilde{r}, \tilde{u}) \|_{\dot{B}^{\frac{d}{2}}}. $$

In addition, according to standard product laws and composite estimates in Lemmas 2.3 and 2.5, the nonlinear terms $(\tilde{S}_1, \tilde{S}_2, \tilde{S}_3, \tilde{S}_4)$ can be handled as

$$ \| (\tilde{S}_1, \tilde{S}_2, \tilde{S}_3, \tilde{S}_4) \|_{\dot{B}^{\frac{d}{2}}} \lesssim \| (y_2, w_2, r_2, u_2) \|_{\dot{B}^{\frac{d}{2}+1}} + \| (w_1, w_2) \|_{\dot{B}^{\frac{d}{2}}} \| (\tilde{y}, \tilde{w}, \tilde{r}, \tilde{u}) \|_{\dot{B}^{\frac{d}{2}}}. $$

Substituting the above two estimates into (4.19) and (4.21) and then employing Grönwall's inequality, we end up with $\| (\tilde{y}, \tilde{w}, \tilde{r}, \tilde{u})(t) \|_{\dot{B}^{\frac{d}{2}}} = 0$ for all $t > 0$, which concludes the proof of Theorem 1.1.

## 4.2. The Cauchy problem of (K) and (PM)

We provide a brief explanation of the proof of the global existence and uniqueness for the system (PM). The proof of the result for the system (K) (Theorem 1.2) follows a very similar procedure, so we omit the details here for brevity. The uniformity of the estimate (1.9) for the system (K) allows us to construct solutions for the system (PM) by taking the limit as the relaxation parameter $\tau \to 0$.

The following lemma states the uniform estimate verified by the solutions of the system (K$\tau$), which are rescaled from the estimate (1.12) obtained in Theorem 1.2 for the system (K).

**Lemma 4.2.** Let $(\alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, u^{\tau})$ be the global solution to the Cauchy problem of the system (K) subject to the initial data $(\alpha_{\pm,0}^{\tau}, \rho_{\pm,0}^{\tau}, u_0^{\tau})$ given by Theorem 1.2 and $(\beta_{\pm}^{\tau}, \varrho_{\pm}^{\tau}, v^{\tau})$ be defined by the diffusive scaling (1.2), then it holds that

$$ \begin{aligned} & \| (\beta_{\pm}^{\tau} - \bar{\alpha}_{\pm}, \varrho_{\pm}^{\tau} - \bar{\rho}_{\pm}) \|_{\tilde{L}^{\infty}} (\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1}) \\ & + \| (\Pi^{\tau} - \bar{P}, \varrho_{\pm}^{\tau} - \bar{\rho}_{\pm}) \|_{L^1} (\dot{B}^{\frac{d}{2}+1}) + \| (\Pi^{\tau} - \bar{P}, \varrho_{\pm}^{\tau} - \bar{\rho}) \|_{\tilde{L}^2} (\dot{B}^{\frac{d}{2}} \cap \dot{B}^{\frac{d}{2}+1}) \\ & + \| v^{\tau} \|_{L^1} (\dot{B}^{\frac{d}{2}} \cap \dot{B}^{\frac{d}{2}+1}) + \| v^{\tau} \|_{\tilde{L}^2} (\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1}) + \frac{1}{\tau} \| z^{\tau} \|_{L^1} (\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}}) \\ & \le C \| (\alpha_{\pm,0}^{\tau} - \bar{\alpha}_{\pm}, \rho_{\pm,0}^{\tau} - \bar{\rho}_{\pm}, u_0^{\tau}) \|_{\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}+1}}, \end{aligned} \quad (4.22) $$

with $z^{\tau} := v^{\tau} + \frac{1}{\varrho^{\tau}} \nabla \Pi^{\tau}$ and $C > 0$ a universal constant.

**Proof of Theorem 1.3.** Assume that the initial data $(\beta_{\pm,0}, \varrho_{\pm,0})$ satisfies (1.13). For any $\tau \in (0, 1)$, we define the regularized data as

$$ (\alpha_{\pm,0}^{\tau}, \rho_{\pm,0}^{\tau}) (x) := \dot{S}_{[\frac{1}{\tau}]} (\beta_{\pm,0}, \varrho_{\pm,0}) (x) \quad \text{and} \quad u_0^{\tau} (x) := 0. $$
Hence, by employing Theorem 1.2 we can obtain a sequence $(\alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, u^{\tau})$, which is the global solution to the system (K) subject to the initial data $(\alpha_{\pm,0}^{\tau}, \rho_{\pm,0}^{\tau}, u_{0}^{\tau})$. Taking the diffusive scaling (1.2), one has that $(\beta_{\pm}^{\tau}, \varrho_{\pm}^{\tau}, v^{\tau})$ is the global solution to the system (Kτ) subject to the initial data $(\alpha_{\pm,0}^{\varepsilon}, \rho_{\pm,0}^{\varepsilon}, u_{0}^{\varepsilon}/\tau)$. In view of the uniform estimate (4.22) established in Lemma 4.2, the Aubin-Lions lemma and the cantor diagonal process, there exists a limit $(\beta_{\pm}; \varrho_{\pm})$ such that as $\tau \to 0$, up to a subsequence, $\chi(\beta_{\pm}^{\tau}, \varrho_{\pm}^{\tau})$ converges to $\chi(\beta_{\pm}, \varrho_{\pm})$ in $C([0, T]; \dot{B}^s)$ ($s < \frac{d}{2} + 1$) strongly for any given time $T > 0$ and $\chi \in C_c^{\infty}(\mathbb{R}^d \times [0, T])$. Thus, we can check that $(\beta_{\pm}, \varrho_{\pm})$ solves the system (PM) in the sense of distributions. Furthermore, taking advantage of the Fatou property and the optimal regularity estimate in Lemma 2.6 for the equation of Π, we can conclude (1.14)–(1.15). Finally, the uniqueness can be obtained in a simple fashion. The interested reader may also refer to [19, 20] for more details. □

# 5. Relaxation limits with convergence rates

## 5.1. Pressure-relaxation limit: (BN) to (K)

In this section, we prove Theorem 1.4 related to the convergence rate of the relaxation process between the system (BN) and the system (K). Let $(\alpha_{\pm}^{\varepsilon, \tau}, \rho_{\pm}^{\varepsilon, \tau}, u^{\varepsilon, \tau})$ and $(\alpha_{\pm}^{\tau}, \rho_{\pm}^{\tau}, u^{\tau})$ be the global solutions to the system (BN) with the initial data $(\alpha_{\pm,0}^{\varepsilon, \tau}, \rho_{\pm,0}^{\varepsilon, \tau}, u_{0}^{\varepsilon, \tau})$ and the system (K) with the initial data $(\alpha_{\pm,0}^{\varepsilon, \tau}, \rho_{\pm,0}^{\varepsilon, \tau}, u_{0}^{\varepsilon, \tau})$ given by Theorems 1.1 and 1.2, respectively. Denote the error variables

$$ (\delta\alpha_{\pm}, \delta\rho_{\pm}, \delta u) := (\alpha_{\pm}^{\varepsilon, \tau} - \alpha_{\pm}^{\tau}, u^{\varepsilon, \tau} - u^{\tau}), \\ (\delta\rho, \delta P_{\pm}, \delta P) := (\rho_{\pm}^{\varepsilon, \tau} - \rho_{\pm}^{\tau}, \rho^{\varepsilon, \tau} - \rho^{\tau}, P_{\pm}(\rho_{\pm}^{\varepsilon, \tau}) - P_{\pm}(\rho_{\pm}^{\tau}), P^{\varepsilon, \tau} - P^{\tau}), $$

and the initial data of $\delta P$

$$ \delta P|_{t=0} = P_{0}^{\varepsilon, \tau} - P_{0}^{\tau}, \quad P_{0}^{\varepsilon} := \alpha_{+,0}^{\varepsilon, \tau} P_{+}(\rho_{+,0}^{\varepsilon, \tau}) + \alpha_{-,0}^{\varepsilon, \tau} P_{-}(\rho_{-,0}^{\varepsilon, \tau}), \quad P_{0}^{\tau} := P_{+}(\rho_{+,0}^{\tau}). \quad (5.1) $$

First, to avoid dealing with difficult nonlinearities in the equation of $\delta\alpha_{\pm}$, we work with the following purely transported variable instead of $\delta\alpha_{\pm}$:

$$ \delta Y := \frac{\alpha_{+}^{\varepsilon, \tau} \rho_{+}^{\varepsilon, \tau}}{\rho^{\varepsilon, \tau}} - \frac{\alpha_{+}^{\tau} \rho_{+}^{\tau}}{\rho^{\tau}}. \quad (5.2) $$

with the initial data

$$ \delta Y|_{t=0} = Y_{0}^{\varepsilon, \tau} - Y_{0}^{\tau}, \quad Y_{0}^{\varepsilon, \tau} := \frac{\alpha_{+,0}^{\varepsilon, \tau} \rho_{+,0}^{\varepsilon, \tau}}{\alpha_{+,0}^{\varepsilon, \tau} \rho_{+,0}^{\varepsilon, \tau} + \alpha_{-,0}^{\varepsilon, \tau} \rho_{-,0}^{\varepsilon, \tau}}, \quad Y_{0}^{\tau} := \frac{\alpha_{+,0}^{\tau} \rho_{+,0}^{\tau}}{\alpha_{+,0}^{\tau} \rho_{+,0}^{\tau} + \alpha_{-,0}^{\tau} \rho_{-,0}^{\tau}}. \quad (5.3) $$

**Lemma 5.1.** For $d \ge 3$, under the assumption (1.16), $\delta Y$ satisfies the following estimate:

$$ \|\delta Y\|_{L_t^\infty(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \lesssim \sqrt{\varepsilon\tau} + o(1) \|\delta u\|_{L^1(\dot{B}^{\frac{d}{2}-1})}. \quad (5.4) $$

**Proof.** Since the equation of $\delta Y$ reads

$$ \partial_t \delta Y + u^{\varepsilon, \tau} \cdot \nabla \delta Y = -\delta u \cdot \nabla \frac{\alpha_{+}^{\tau} \rho_{+}^{\tau}}{\rho^{\tau}}, $$
Lemma 2.7 and the product law (2.4) for $d \ge 3$ give

$$ \|\delta Y\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \lesssim \exp \left( \|u^{\varepsilon, \tau}\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} \right) \\ \times \left( \sqrt{\varepsilon\tau} + \|\delta u\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \|\nabla \frac{\alpha_+^\tau \rho_+^\tau}{\rho^\tau}\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})} \right). $$

This together with the uniform estimate (1.9) leads to (5.4). □

We are now ready to estimate $(\delta\alpha_\pm, \delta\rho_\pm, \delta P_\pm, \delta P)$. It is easy to verify that $P^{\varepsilon, \tau}$ satisfies

$$ \begin{aligned} \partial_t P^{\varepsilon, \tau} + u^{\varepsilon, \tau} \cdot \nabla P^{\varepsilon, \tau} = & - (\gamma_+ \alpha_+^{\varepsilon, \tau} P_+ (\rho_+^{\varepsilon, \tau}) + \gamma_- \alpha_-^{\varepsilon, \tau} P_- (\rho_-^{\varepsilon, \tau})) \operatorname{div} u^{\varepsilon, \tau} \\ & - \alpha_+^{\varepsilon, \tau} \alpha_-^{\varepsilon, \tau} ((\gamma_+ - 1) P_+ (\rho_+^{\varepsilon, \tau}) \\ & \qquad \qquad \qquad - (\gamma_- - 1) P_- (\rho_-^{\varepsilon, \tau})) \frac{P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau})}{\varepsilon}. \end{aligned} \quad (5.5) $$

And the equation of $P^\tau$ reads

$$ \partial_t P^\tau + u^\tau \cdot \nabla P^\tau + \frac{\gamma_+ \gamma_- P^\tau}{\gamma_+ \alpha_-^\tau + \gamma_- \alpha_+^\tau} \operatorname{div} u^\tau = 0. \quad (5.6) $$

However it is not suitable to estimate $\delta P$ directly from (5.5) and (5.6) as the decay rate of $P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau})$ can not be faster than $\varepsilon$ in view of (1.9). To overcome this difficulty, we introduce an auxiliary unknown $Q^{\varepsilon, \tau} := P^{\varepsilon, \tau} - \Gamma_1^{\varepsilon, \tau} (P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau}))$ which verifies

$$ \begin{aligned} \partial_t Q^{\varepsilon, \tau} + u^{\varepsilon, \tau} \cdot \nabla Q^{\varepsilon, \tau} + \Gamma_2^{\varepsilon, \tau} \operatorname{div} u^{\varepsilon, \tau} \\ = - \Gamma_3^{\varepsilon, \tau} (P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau})) \operatorname{div} u^{\varepsilon, \tau} \\ + (\partial_t \Gamma_1^{\varepsilon, \tau} + u^{\varepsilon, \tau} \cdot \nabla \Gamma_1^{\varepsilon, \tau}) (P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau})), \end{aligned} \quad (5.7) $$

with

$$ \left\{ \begin{aligned} \Gamma_1^{\varepsilon, \tau} &:= \frac{\alpha_+^{\varepsilon, \tau} \alpha_-^{\varepsilon, \tau} ((\gamma_+ - 1) P_+ (\rho_+^{\varepsilon, \tau}) - (\gamma_- - 1) P_- (\rho_-^{\varepsilon, \tau}))}{\gamma_+ \alpha_-^{\varepsilon, \tau} P_+ (\rho_+^{\varepsilon, \tau}) + \gamma_- \alpha_+^{\varepsilon, \tau} P_- (\rho_-^{\varepsilon, \tau})}, \\ \Gamma_2^{\varepsilon, \tau} &:= \frac{\gamma_+ \gamma_- P_+ (\rho_+^{\varepsilon, \tau}) P_- (\rho_-^{\varepsilon, \tau})}{\gamma_+ \alpha_- P_+ (\rho_+^{\varepsilon, \tau}) + \gamma_- \alpha_+ P_- (\rho_-^{\varepsilon, \tau})}, \\ \Gamma_3^{\varepsilon, \tau} &:= \frac{\alpha_+^{\varepsilon, \tau} \alpha_-^{\varepsilon, \tau} (\gamma_+ P_+ (\rho_+^{\varepsilon, \tau}) - \gamma_- P_- (\rho_-^{\varepsilon, \tau}))}{\gamma_+ \alpha_- P_+ (\rho_+^{\varepsilon, \tau}) + \gamma_- \alpha_+ P_- (\rho_-^{\varepsilon, \tau})}. \end{aligned} \right. $$

With this formulation, it will be possible to derive the $O(\varepsilon)$ bounds for the last term on the right-hand side of (5.7). Define

$$ \delta Q := P^{\varepsilon, \tau} - P^{\tau} - \Gamma_1^{\varepsilon, \tau} (P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau})). \quad (5.8) $$

The next lemma implies that to estimate $(\delta\alpha_\pm, \delta\rho_\pm, \delta P_\pm, \delta P)$, it is sufficient to control $(\delta Y, \delta Q, P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau}))$.
**Lemma 5.2.** For $d \ge 3$, under the assumption (1.16), the following estimates follow:

$$ \left\{ \begin{array}{l} \| (\delta \alpha_{\pm}, \delta \rho_{\pm}, \delta \rho) \|_{\tilde{L}_i^{\infty} (\hat{B}_2^{\frac{d}{2}-1} \cap \hat{B}_2^{\frac{d}{2}})} \\ \lesssim \| (\delta Y, \delta Q, P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau})) \|_{\tilde{L}_i^{\infty} (\hat{B}_2^{\frac{d}{2}-2} \cap \hat{B}_2^{\frac{d}{2}-1})}, \quad (5.9) \\ \| \delta \rho_{\pm} \|_{\tilde{L}_i^2 (\hat{B}_2^{\frac{d}{2}-1})} \lesssim \| (\delta Q, P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau})) \|_{\tilde{L}_i^2 (\hat{B}_2^{\frac{d}{2}-1})}. \end{array} \right. $$

**Proof.** Due to (5.2) and

$$ \delta \rho = (\rho_+^{\varepsilon, \tau} - \rho_-^{\varepsilon, \tau}) \delta \alpha_+ + \alpha_+^{\tau} \delta \rho_+ + \alpha_-^{\tau} \delta \rho_-, \quad (5.10) $$

it holds that

$$ \begin{aligned} \delta Y &= \frac{1}{\rho^{\varepsilon, \tau} \rho^{\tau}} (\rho_+^{\varepsilon, \tau} \rho^{\tau} \delta \alpha_+^{\varepsilon, \tau} + \alpha_+^{\tau} \rho^{\tau} \delta \rho_+ - \alpha_+^{\tau} \rho_+^{\tau} \delta \rho) \\ &= \frac{1}{\rho^{\varepsilon, \tau} \rho^{\tau}} ((\alpha_-^{\tau} \rho_+^{\varepsilon, \tau} \rho_-^{\tau} + \alpha_+^{\tau} \rho_+^{\tau} \rho_-^{\varepsilon, \tau}) \delta \alpha_+ + \alpha_+^{\tau} \alpha_-^{\tau} \rho_-^{\tau} \delta \rho_+ - \alpha_+^{\tau} \alpha_-^{\tau} \rho_+^{\tau} \delta \rho_-). \end{aligned} $$

This implies

$$ \delta \alpha_+ = \frac{1}{\alpha_-^{\tau} \rho_+^{\varepsilon, \tau} \rho_-^{\tau} + \alpha_+^{\tau} \rho_+^{\varepsilon, \tau} \rho_-^{\tau}} (\rho^{\varepsilon, \tau} \rho^{\tau} \delta Y - \alpha_+^{\tau} \alpha_-^{\tau} \rho_-^{\tau} \delta \rho_+ + \alpha_+^{\tau} \alpha_-^{\tau} \rho_+^{\tau} \delta \rho_-). \quad (5.11) $$

Inserting (5.11) into (5.10), we have

$$ \begin{aligned} \delta \rho ={}& \frac{\rho_+^{\varepsilon, \tau} - \rho_-^{\varepsilon, \tau}}{\alpha_-^{\tau} \rho_+^{\varepsilon, \tau} \rho_-^{\tau} + \alpha_+^{\tau} \rho_+^{\varepsilon, \tau} \rho_-^{\tau}} (\rho^{\varepsilon, \tau} \rho^{\tau} \delta Y - \alpha_+^{\tau} \alpha_-^{\tau} \rho_-^{\tau} \delta \rho_+ + \alpha_+^{\tau} \alpha_-^{\tau} \rho_+^{\tau} \delta \rho_-) \\ &+ \alpha_+^{\tau} \delta \rho_+ + \alpha_-^{\tau} \delta \rho_- \end{aligned} \quad (5.12) $$

Moreover, we have

$$ \delta P_{\pm} = \delta \rho_{\pm} \int_{0}^{1} P'_{\pm} (\theta \rho_{\pm}^{\varepsilon, \tau} + (1-\theta) \rho_{\pm}^{\tau}) d\theta, \quad \delta P = \alpha_{+}^{\varepsilon, \tau} (P_{+}^{\varepsilon, \tau} - P_{-}^{\varepsilon, \tau}) + \delta P_{-}. \quad (5.13) $$

Using the previous uniform estimates (1.9) and (4.22), the product laws (2.3)–(2.4) and the composition estimates (2.6)–(2.7), for some constant states $\bar{\Gamma}_i > 0$ ($i = 1, 2, 3$), we have

$$ \sum_{i=1}^{3} \left( \| \Gamma_i - \bar{\Gamma}_i \|_{\tilde{L}_i^{\infty} (\hat{B}_2^{\frac{d}{2}-1} \cap \hat{B}_2^{\frac{d}{2}+1})} + \| \partial_i \Gamma_i \|_{\tilde{L}_i^2 (\hat{B}_2^{\frac{d}{2}})} \right) = o(1). \quad (5.14) $$

Therefore, (5.9) follows from (5.8), (5.11)–(5.14), the product laws (2.3)–(2.4) and the fact $\delta \alpha_+ = -\delta \alpha_-$.

The next lemma pertains to $O(\sqrt{\varepsilon \tau})$ bounds for $P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau})$, which leads to the convergence rate $\sqrt{\varepsilon \tau}$.

**Lemma 5.3.** For $d \ge 3$, under the assumption (1.16), the following estimate is valid:

$$ \begin{aligned} & \| P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau}) \|_{\tilde{L}_i^{\infty} (\hat{B}_2^{\frac{d}{2}-2} \cap \hat{B}_2^{\frac{d}{2}-1})} \\ & \qquad + \frac{1}{\sqrt{\varepsilon}} \| P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau}) \|_{\tilde{L}_i^2 (\hat{B}_2^{\frac{d}{2}-1} \cap \hat{B}_2^{\frac{d}{2}-1})} \lesssim \sqrt{\varepsilon \tau}. \end{aligned} \quad (5.15) $$
**Proof.** It is easy to verify from (BN) that $P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})$ satisfies the damped equation

$$ 
\begin{aligned}
& \partial_t (P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})) + u^{\varepsilon,\tau} \cdot \nabla (P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})) \\
& \qquad + \frac{c_*}{\varepsilon} (P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})) \\
& = ((\gamma_+ \alpha_-^{\varepsilon,\tau} P_+(\rho_+^{\varepsilon,\tau}) + \gamma_- \alpha_+^{\varepsilon,\tau} P_-(\rho_-^{\varepsilon,\tau})) - c_*) \frac{1}{\varepsilon} (P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})) \\
& \qquad - (\gamma_+ P_+(\rho_+^{\varepsilon,\tau}) - \gamma_- P_-(\rho_-^{\varepsilon,\tau})) \operatorname{div} u^{\varepsilon,\tau} := W_1 + W_2,
\end{aligned} 
\quad (5.16) $$

with $c_* := (\gamma_+ \bar{\alpha}_- + \gamma_- \bar{\alpha}_+)^{\bar{P}}$. Then the $L^2$-in-time type estimate in Lemma 2.7 for the damped transport equation (5.16) leads to

$$ 
\begin{aligned}
& \| P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau}) \|_{\bar{L}_t^\infty(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\varepsilon}} \| P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau}) \|_{\bar{L}_t^2(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \\
& \lesssim \exp \left( \| u^{\varepsilon,\tau} \|_{\bar{L}_t^1(\dot{B}^{\frac{d}{2}+1})} \right) \left( \sqrt{\varepsilon \tau} + \sqrt{\varepsilon} \|(W_1, W_2)\|_{\bar{L}_t^2(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \right).
\end{aligned} 
$$

By (1.9) and (2.4), there holds

$$ 
\begin{aligned}
\sqrt{\varepsilon} \|W_1\|_{\bar{L}_t^2(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} & \lesssim \| (\alpha_{\pm}^{\varepsilon,\tau} - \bar{\alpha}_{\pm}, \rho_{\pm}^{\varepsilon,\tau} - \bar{\rho}_{\pm}) \|_{\bar{L}_t^{\infty}(\dot{B}^{\frac{d}{2}})} \frac{1}{\sqrt{\varepsilon}} \|P_+(\rho_+^{\varepsilon,\tau}) \\
& \qquad - P_-(\rho_-^{\varepsilon,\tau})\|_{\bar{L}_t^2(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \\
& \lesssim o(1) \frac{1}{\sqrt{\varepsilon}} \|P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})\|_{\bar{L}_t^2(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})},
\end{aligned} 
$$

and

$$ \sqrt{\varepsilon} \|W_2\|_{\bar{L}_t^2(\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} \lesssim \sqrt{\varepsilon} \|u^{\varepsilon,\tau}\|_{\bar{L}_t^2(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})} \lesssim \sqrt{\varepsilon \tau}. $$

Therefore, we gain (5.15). □

We are going to estimate $(\delta Q, \delta u)$. By virtue of (BN), (K) and (5.6)–(5.7), $(\delta Q, \delta u)$ satisfies the following equations of damped Euler type with rough coefficients:

$$ 
\left\{ 
\begin{aligned}
\partial_t \delta Q + u^{\varepsilon,\tau} \cdot \nabla \delta Q + \Gamma_2^{\varepsilon,\tau} \operatorname{div} \delta u &= \delta F_1, \\
\partial_t \delta u + u^{\varepsilon,\tau} \cdot \nabla \delta u + \frac{1}{\bar{\rho}} \nabla \delta Q + \left( \frac{1}{\rho^{\varepsilon,\tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} + \frac{\delta u}{\tau} &= \delta F_2,
\end{aligned} 
\right. 
\quad (5.17) $$

with the nonlinear terms

$$ 
\left\{ 
\begin{aligned}
\delta F_1 &= -\delta u \cdot \nabla P^{\tau} - \left( \Gamma_2^{\varepsilon,\tau} - \frac{\gamma_+ \gamma_- P^{\tau}}{\gamma_+ \alpha_-^{\varepsilon,\tau} + \gamma_- \alpha_+^{\varepsilon,\tau}} \right) \operatorname{div} u^{\varepsilon,\tau} \\
& \qquad - \Gamma_3^{\varepsilon,\tau} (P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})) \operatorname{div} u^{\varepsilon,\tau} + (\partial_t \Gamma_1^{\varepsilon,\tau} + u^{\varepsilon,\tau} \cdot \nabla \Gamma_1^{\varepsilon,\tau}) (P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau})), \\
\delta F_2 &:= -\delta u \cdot \nabla u^{\tau} - \frac{1}{\rho^{\varepsilon,\tau}} \nabla (\Gamma_1^{\varepsilon,\tau} (P_+(\rho_+^{\varepsilon,\tau}) - P_-(\rho_-^{\varepsilon,\tau}))).
\end{aligned} 
\right. 
$$
In order to establish the uniform-in-$\tau$ convergence estimates, we follow the ideas in section 3 to overcome the issue caused by the overdamping phenomenon.

**Lemma 5.4.** Let $d \ge 3$, $0 < \varepsilon \le \tau \le 1$, and the threshold $J_{\tau}$ be given by (2.1). Then under the assumption (1.16), there holds

$$ \begin{aligned} & \| (\delta Q, \delta u) \|_{\tilde{L}_i^\infty (\dot{B}_2^{d/2-2} \cap \dot{B}_2^{d/2-1})} + \tau \| \delta Q \|_{L_i^1 (\dot{B}_2^d)}^\ell + \| \delta Q \|_{L_i^1 (\dot{B}_2^{d/2-1})}^h \\ & \qquad + \sqrt{\tau} \| \delta Q \|_{\tilde{L}_i^2 (\dot{B}_2^{d/2-1})} + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_i^2 (\dot{B}_2^{d/2-2} \cap \dot{B}_2^{d/2-1})} + \| \delta u \|_{L_i^1 (\dot{B}_2^{d/2-1})} \\ & \lesssim \sqrt{\varepsilon \tau} + o(1) \| \delta Y \|_{\tilde{L}_i^\infty (\dot{B}_2^{d/2-2} \cap \dot{B}_2^{d/2-1})}. \end{aligned} \quad (5.18) $$

**Proof.** As in section 3, we split the proof into three parts:

• Step 1: $\dot{B}_2^{d/2-2}$-estimates in low frequencies

We introduce the new damped mode (effective flux)

$$ \delta z := \delta u + \frac{\tau}{\rho^{\varepsilon, \tau}} \nabla \delta Q + \tau \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau}, $$

so that (5.17) is rewritten as

$$ \left\{ \begin{aligned} & \partial_t \delta Q - \frac{\bar{\Gamma}_2 \tau}{\bar{\rho}} \Delta \delta Q = -\bar{\Gamma}_2 \operatorname{div} z + \delta F_3, \\ & \partial_t \delta z + \frac{\delta z}{\tau} = \frac{\tau}{\bar{\rho}} \nabla \left( \frac{\bar{\Gamma}_2 \tau}{\bar{\rho}} \Delta \delta Q - \bar{\Gamma}_2 \operatorname{div} z \right) + \delta F_4, \end{aligned} \right. \quad (5.19) $$

where $\bar{\Gamma}_2 > 0$ is the constant state of $\Gamma_2^{\varepsilon, \tau}$, and $\delta F_i$ ($i = 3, 4$) is defined by

$$ \left\{ \begin{aligned} \delta F_3 := & -u^{\varepsilon, \tau} \cdot \nabla \delta Q - (\Gamma_2^{\varepsilon, \tau} - \bar{\Gamma}_2) \operatorname{div} \delta u \\ & + \bar{\Gamma}_2 \tau \operatorname{div} \left( \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla \delta Q + \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right) + \delta F_1, \\ \delta F_4 := & -u^{\varepsilon, \tau} \cdot \nabla \delta u + \frac{\tau}{\bar{\rho}} \nabla \delta F_3 \\ & + \tau \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla \partial_t \delta Q - \tau \partial_t \left( \frac{1}{\rho^{\varepsilon, \tau}} \right) \nabla \delta Q + \tau \partial_t \left( \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right) + \delta F_2. \end{aligned} \right. $$

Then by similar arguments used to get (3.13)–(3.16), we deduce from (5.19) and the choice (2.1) of the threshold $J_{\tau}$ that

$$ \begin{aligned} & \| (\delta Q, \delta z) \|_{\tilde{L}_i^\infty (\dot{B}_2^{d/2-2})}^\ell + \tau \| \delta Q \|_{L_i^1 (\dot{B}_2^d)}^\ell + \tau \| \partial_t \delta Q \|_{L_i^1 (\dot{B}_2^{d/2-2})}^\ell \\ & \qquad + \sqrt{\tau} \| \delta Q \|_{\tilde{L}_i^2 (\dot{B}_2^{d/2-1})}^\ell + \frac{1}{\tau} \| \delta z \|_{L_i^1 (\dot{B}_2^{d/2-2})}^\ell \\ & \lesssim \sqrt{\varepsilon \tau} + \| (\delta F_3, \delta F_4) \|_{L_i^1 (\dot{B}_2^{d/2-2})}^\ell. \end{aligned} \quad (5.20) $$
We first estimate $\delta F_3$. From (1.9), (5.14) and the product map $\dot{B}_{\frac{d}{2}-2} \times \dot{B}_{\frac{d}{2}} \to \dot{B}_{\frac{d}{2}-2}$ for $d \ge 3$, one obtains

$$ \begin{aligned} & \|u^{\epsilon, \tau} \cdot \nabla \delta Q + (\Gamma_2^{\epsilon, \tau} - \bar{\Gamma}_2) \operatorname{div} \delta u\|_{L_t^1(\dot{B}_{\frac{d}{2}-2}^d)} \\ & \lesssim \frac{1}{\sqrt{\tau}} \|u^{\epsilon, \tau}\|_{\bar{L}_t^2(\dot{B}_{\frac{d}{2}}^d)} \sqrt{\tau} \|\delta Q\|_{\bar{L}_t^2(\dot{B}_{\frac{d}{2}-1}^d)} + \|\Gamma_2^{\epsilon, \tau} - \bar{\Gamma}_2\|_{\bar{L}_t^\infty(\dot{B}_{\frac{d}{2}}^d)} \|\delta u\|_{L_t^1(\dot{B}_{\frac{d}{2}-1}^d)} \quad (5.21) \\ & \lesssim o(1) \left( \sqrt{\tau} \|\delta Q\|_{\bar{L}_t^2(\dot{B}_{\frac{d}{2}-1}^d)} + \|\delta u\|_{L_t^1(\dot{B}_{\frac{d}{2}-1}^d)} \right). \end{aligned} $$

By virtue of (1.9), (2.2), (2.4), (5.9) and (5.15), we also have

$$ \begin{aligned} & \tau \| \operatorname{div} \left( \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right) \|_{L_t^1(\dot{B}_{\frac{d}{2}-2}^d)}^{\ell} \\ & \lesssim \|\delta \rho\|_{\bar{L}_t^{\infty}(\dot{B}_{\frac{d}{2}-1}^d)} \tau \|P^{\tau} - \bar{P}\|_{L_t^1(\dot{B}_{\frac{d}{2}+1}^d)} \\ & \lesssim o(1) \| (\delta Y, \delta Q) \|_{\bar{L}_t^{\infty}(\dot{B}_{\frac{d}{2}-1}^d)} + \sqrt{\epsilon \tau}. \end{aligned} $$

As in the previous analysis (3.19), the tricky nonlinear term can be estimated as

$$ \begin{aligned} & \tau \| \operatorname{div} \left( \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla \delta Q \right) \|_{L_t^1(\dot{B}_{\frac{d}{2}-2}^d)}^{\ell} \\ & \lesssim \tau \left\| \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla \delta Q^{\ell} \right\|_{L_t^1(\dot{B}_{\frac{d}{2}-1}^d)}^{\ell} + \left\| \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla \delta Q^h \right\|_{L_t^1(\dot{B}_{\frac{d}{2}-2}^d)}^{\ell} \\ & \lesssim o(1) \left( \tau \|\delta Q\|_{L_t^1(\dot{B}_{\frac{d}{2}}^d)}^{\ell} + \|\delta Q\|_{L_t^1(\dot{B}_{\frac{d}{2}-1}^d)}^{h} \right). \end{aligned} $$

Similarly, one can show

$$ \begin{aligned} \|\delta F_1\|_{L_t^1(\dot{B}_{\frac{d}{2}-2}^d)} & \lesssim \frac{1}{\sqrt{\tau}} \|\delta u\|_{\bar{L}_t^2(\dot{B}_{\frac{d}{2}-2}^d)} \sqrt{\tau} \|Q^\tau\|_{\bar{L}_t^2(\dot{B}_{\frac{d}{2}}^d)} \\ & + \left\| \left( \Gamma_2^{\epsilon, \tau} - \frac{\gamma_+ \gamma_- P^\tau}{\gamma_+ \alpha_-^\tau + \gamma_- \alpha_+^\tau} \right) \right\|_{\bar{L}_t^\infty(\dot{B}_{\frac{d}{2}-2}^d)} \|u^\tau\|_{L_t^1(\dot{B}_{\frac{d}{2}+1}^d)} \\ & + \|P_+(\rho_+^{\epsilon, \tau}) - P_-(\rho_-^{\epsilon, \tau})\|_{\bar{L}_t^\infty(\dot{B}_{\frac{d}{2}-2}^d)} \|u^{\epsilon, \tau}\|_{L_t^1(\dot{B}_{\frac{d}{2}+1}^d)} \\ & + \left( \|\partial_t \Gamma_1^{\epsilon, \tau}\|_{L_t^1(\dot{B}_{\frac{d}{2}}^d)} + \|u^{\epsilon, \tau}\|_{L_t^1(\dot{B}_{\frac{d}{2}}^d)} \|\nabla \Gamma_1^{\epsilon, \tau}\|_{\bar{L}_t^\infty(\dot{B}_{\frac{d}{2}}^d)} \right) \\ & \times \|P_+(\rho_+^{\epsilon, \tau}) - P_-(\rho_-^{\epsilon, \tau})\|_{\bar{L}_t^\infty(\dot{B}_{\frac{d}{2}-2}^d)}, \end{aligned} $$

which implies

$$ \|\delta F_1\|_{L_t^1(\dot{B}_{\frac{d}{2}-2}^d)} \lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\bar{L}_t^\infty(\dot{B}_{\frac{d}{2}-2}^d)} + \frac{1}{\sqrt{\tau}} \|\delta u\|_{\bar{L}_t^2(\dot{B}_{\frac{d}{2}-2}^d)} \right) + \sqrt{\epsilon \tau}. \quad (5.22) $$
Therefore, we have

$$ \|\delta F_3\|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}^{\ell} \lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2}-2})} + \|\delta Q\|_{L_t^1(\dot{B}^{\frac{d}{2}})} + \sqrt{\tau} \|\delta Q\|_{\tilde{L}_t^2(\dot{B}^{\frac{d}{2}-1})} \right. \\ \left. + \frac{1}{\sqrt{\tau}} \|\delta u\|_{\tilde{L}_t^1(\dot{B}^{\frac{d}{2}-2})} + \|\delta u\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \right) + \sqrt{\varepsilon\tau}. \quad (5.23) $$

We turn to the estimate of $\delta F_4$. Similar calculations give

$$ \| -u^{\varepsilon, \tau} \cdot \nabla \delta u + \frac{\tau}{\rho^{\varepsilon, \tau}} \nabla \delta F_3 \|_{L_t^1(\dot{B}^{\frac{d}{2}-2})} \lesssim o(1) \|\delta u\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} + \|\delta F_3\|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}, $$

and

$$ \| \tau \partial_t \left( \frac{1}{\rho^{\varepsilon, \tau}} \right) \nabla \delta Q \|_{L_t^1(\dot{B}^{\frac{d}{2}-2})} \lesssim \| \partial_t \rho^{\varepsilon, \tau} \|_{L_t^1(\dot{B}^{\frac{d}{2}})} \| \delta Q \|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2}-1})} \lesssim o(1) \| \delta Q \|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2}-1})}. $$

For the third difficult term in $\delta F_4$, we apply (2.2), the product law (2.4) for $d \ge 3$ and the fact

$$ \tau \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla \partial_t \delta Q = \tau \nabla \left( \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \partial_t \delta Q \right) - \tau \nabla \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \partial_t \delta Q $$

to have

$$ \begin{aligned} \| \tau \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla \partial_t \delta Q \|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}^{\ell} & \lesssim \tau \left\| \nabla \left( \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \partial_t \delta Q \right) \right\|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}^{\ell} \\ & + \tau \left\| \nabla \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \partial_t \delta Q \right\|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}^{\ell} \\ & \lesssim \left\| \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \partial_t \delta Q \right\|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}^{\ell} \\ & + \left\| \nabla \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \partial_t \delta Q \right\|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}^{\ell} \\ & \lesssim \| \rho^{\varepsilon, \tau} - \bar{\rho} \|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2} \cap \dot{B}^{\frac{d}{2}+1})} \| \partial_t \delta Q \|_{L_t^1(\dot{B}^{\frac{d}{2}-2})} \\ & \lesssim o(1) \| \partial_t \delta Q \|_{L_t^1(\dot{B}^{\frac{d}{2}-2})}. \end{aligned} $$

Similarly, the term $\delta F_2$ can be easily estimated as follows:

$$ \begin{aligned} \|\delta F_2\|_{L_t^1(\dot{B}^{\frac{d}{2}-2})} & \lesssim \|\delta u\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \|u^\tau\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} + \|P_+(\rho_+^{\varepsilon, \tau}) - P_-(\rho_-^{\varepsilon, \tau})\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \\ & \lesssim o(1) \|\delta u\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} + \sqrt{\varepsilon\tau} + \varepsilon \\ & \lesssim o(1) \|\delta u\|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} + \sqrt{\varepsilon\tau}. \end{aligned} $$

To bound the term $\tau \partial_t \left( \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\bar{\rho}} \right) \nabla P^\tau \right)$, noticing that

$$ \partial_t \delta \rho = - \operatorname{div} (\delta \rho u^{\varepsilon, \tau} + \rho^\tau \delta u), $$
we use (1.9), (5.9), (5.15) and (2.6)–(2.7) that

$$
\begin{aligned}
\| \tau \partial_t \left( \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\rho^\tau} \right) \nabla P^\tau \right) \|_{L_t^1 (\dot{B}^{\frac{d}{2}-2})} \\
&\lesssim \| \partial_t \delta \rho \|_{L_t^1 (\dot{B}^{\frac{d}{2}-2})} \| \nabla P^\tau \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}})} + \| \delta \rho \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-2})} \tau \| \nabla P^\tau \|_{L_t^1 (\dot{B}^{\frac{d}{2}})} \\
&\lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} + \| \delta u \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})} \right) + \sqrt{\varepsilon \tau}.
\end{aligned}
$$

We thence get

$$
\begin{aligned}
\| \delta F_3 \|_{L_t^1 (\dot{B}^{\frac{d}{2}-2})}^\ell &\lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-2})} + \| \delta Q \|_{L_t^1 (\dot{B}^{\frac{d}{2}})} + \sqrt{\tau} \| \delta Q \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})} \right. \\
&\qquad \left. + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-2})} + \| \delta u \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})} \right) + \sqrt{\varepsilon \tau}.
\end{aligned}
\quad (5.24)
$$

Substituting the above estimates (5.23)–(5.24) into (5.20) and taking advantage of $\delta u = \delta z - \frac{\tau}{\rho^{\varepsilon, \tau}} \nabla \delta Q - \tau(\frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\rho^\tau}) \nabla P^\tau$, we obtain

$$
\begin{aligned}
& \| (\delta Q, \delta z) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-2})}^\ell + \tau \| \delta Q \|_{L_t^1 (\dot{B}^{\frac{d}{2}})}^\ell + \sqrt{\tau} \| \delta Q \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^\ell + \frac{1}{\tau} \| \delta z \|_{L_t^1 (\dot{B}^{\frac{d}{2}-2})}^\ell \\
& + \| \delta u \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-2})}^\ell + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-2})}^\ell + \| \delta u \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^\ell \\
& \lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-2} \cap \dot{B}^{\frac{d}{2}-1})} + \tau \| \delta Q \|_{L_t^1 (\dot{B}^{\frac{d}{2}})}^\ell \right. \\
& \qquad \left. + \| \delta Q \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^h + \sqrt{\tau} \| \delta Q \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^\ell \right. \\
& \qquad \left. + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-2})} + \| \delta u \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})} + \| \partial_t \delta Q \|_{L_t^1 (\dot{B}^{\frac{d}{2}-2})} \right) + \sqrt{\varepsilon \tau}.
\end{aligned}
\quad (5.25)
$$

• Step 2: $\dot{B}^{\frac{d}{2}-2}$-estimates of $(\delta Q, \delta u)$ in high frequencies

Applying $\Delta_j$ to (5.17), one gets

$$
\left\{
\begin{aligned}
& \partial_t \dot{\Delta}_j \delta Q + u^{\varepsilon, \tau} \cdot \nabla \dot{\Delta}_j \delta Q + \Gamma_2^{\varepsilon, \tau} \operatorname{div} \dot{\Delta}_j \delta u = \dot{\Delta}_j \delta F_1 + \delta R_{1,j}, \\
& \partial_t \dot{\Delta}_j \delta u + u^{\varepsilon, \tau} \cdot \nabla \dot{\Delta}_j \delta u + \frac{1}{\rho^{\varepsilon, \tau}} \nabla \dot{\Delta}_j \delta Q + \frac{\dot{\Delta}_j \delta u}{\tau} \\
& \qquad = - \dot{\Delta}_j \left( \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\rho^\tau} \right) \nabla P^\tau \right) + \dot{\Delta}_j \delta F_2 + \delta R_{2,j} + \delta R_{3,j},
\end{aligned}
\right.
$$
with

$$
\left\{
\begin{aligned}
\delta R_{1,j} &:= [u^{\epsilon, \tau}, \dot{\Delta}_j] \nabla \delta Q + [\Gamma_2^{\epsilon, \tau}, \dot{\Delta}_j] \operatorname{div} \dot{\Delta}_j \delta u, \\
\delta R_{2,j} &:= [u^{\epsilon, \tau}, \dot{\Delta}_j] \nabla \delta u, \\
\delta R_{3,j} &:= \left[ \frac{1}{\rho^{\epsilon, \tau}}, \dot{\Delta}_j \right] \nabla \delta Q.
\end{aligned}
\right.
$$

Similarly to the high-frequency analysis in section 3.2, one gains

$$
\begin{align}
& \frac{d}{dt} \int_{\mathbb{R}^d} \left( \frac{1}{\rho^{\epsilon, \tau}} |\dot{\Delta}_j \delta Q|^2 + \Gamma_2^{\epsilon, \tau} |\dot{\Delta}_j \delta u|^2 \right) dx + \frac{1}{\tau} \|\dot{\Delta}_j \delta u\|_{L^2}^2 \nonumber \\
& \lesssim \left\| \left( \operatorname{div} u^{\epsilon, \tau}, \nabla \Gamma_2^{\epsilon, \tau}, \nabla \frac{1}{\rho^{\epsilon, \tau}}, \partial_t \frac{1}{\rho^{\epsilon, \tau}}, \partial_t \Gamma_2^{\epsilon, \tau} \right) \right\|_{\tilde{L}^{\infty}} \|\dot{\Delta}_j \delta Q\|_{L^2} \|\dot{\Delta}_j \delta u\|_{L^2} \tag{5.26} \\
& \qquad + \left\| \dot{\Delta}_j \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right\|_{L^2} \|\dot{\Delta}_j \delta u\|_{L^2} + \| \dot{\Delta}_j (\delta F_1, \delta F_2) \|_{L^2} \|\dot{\Delta}_j (\delta Q, \delta u)\|_{L^2} \nonumber \\
& \qquad + \| \delta R_{1,j} \|_{L^2} \|\dot{\Delta}_j \delta Q\|_{L^2} + \| (\delta R_{2,j}, \delta R_{3,j}) \|_{L^2} \|\dot{\Delta}_j \delta u\|_{L^2}, \nonumber
\end{align}
$$

and the cross term

$$
\begin{align}
& \frac{d}{dt} \int_{\mathbb{R}^d} \dot{\Delta}_j \delta u \cdot \nabla \dot{\Delta}_j \nabla \delta P dx + \int_{\mathbb{R}^d} \left( \frac{1}{\rho^{\epsilon, \tau}} |\nabla \dot{\Delta}_j \delta Q|^2 - \Gamma_2^{\epsilon, \tau} |\operatorname{div} \dot{\Delta}_j \delta u_j|^2 + \frac{1}{\tau} \dot{\Delta}_j \delta u \cdot \nabla \dot{\Delta}_j \nabla \delta P \right) dx \nonumber \\
& \lesssim \left( \|u^{\epsilon, \tau}\|_{\tilde{L}^{\infty}} \|\nabla \dot{\Delta}_j \delta u\|_{L^2}^2 + \left\| \dot{\Delta}_j \left( \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right) \right\|_{L^2} \right. \tag{5.27} \\
& \qquad \left. + \left\| (\dot{\Delta}_j \delta F_2, \delta R_{2,j}, \delta R_{3,j}) \right\|_{L^2} \|\nabla \dot{\Delta}_j \nabla \delta P\|_{L^2} \right. \nonumber \\
& \qquad \left. + \left( \|u^{\epsilon, \tau}\|_{\tilde{L}^{\infty}} \|\nabla \dot{\Delta}_j \delta Q\|_{L^2}^2 + \|\dot{\Delta}_j \delta F_1\|_{L^2}^2 + \|\delta R_{1,j}\|_{L^2} \right) \|\delta u\|_{L^2} \right). \nonumber
\end{align}
$$

For all $j \ge J_\tau$, multiplying (5.27) by a suitable small constant and adding the resulting inequal-
ity and (5.26) together, we can derive the Lyapunov inequality similar to (3.34)–(3.41) and then
show the following $L^1$-in-time type estimate:

$$
\begin{align*}
& \tau \| (\delta Q, \delta u) \|_{\tilde{L}^\infty}^h (\dot{B}^{\frac{d}{2}-1}) + \| (\delta Q, \delta u) \|_{L^1}^h (\dot{B}^{\frac{d}{2}-1}) \\
& \lesssim \sqrt{\tau \epsilon} + \left( \| u^{\epsilon, \tau} \|_{L^1} (\dot{B}^{\frac{d}{2} \cap \dot{B}^{\frac{d}{2}+1}}) + \| (\partial_t \rho^{\epsilon, \tau}, \partial_t \Gamma_2^{\epsilon, \tau}) \|_{L^1} (\dot{B}^{\frac{d}{2}}) \right) \tau \| (\delta Q, \delta u) \|_{\tilde{L}^\infty}^h (\dot{B}^{\frac{d}{2}-1}) \\
& \qquad + \| (\rho^{\epsilon, \tau} - \bar{\rho}, \Gamma_2^{\epsilon, \tau} - \bar{\Gamma}_2) \|_{\tilde{L}^\infty} (\dot{B}^{\frac{d}{2}+1}) \| (\delta Q, \delta u) \|_{L^1}^h (\dot{B}^{\frac{d}{2}-1}) + \tau \| (\delta F_1, \delta F_2) \|_{L^1}^h (\dot{B}^{\frac{d}{2}-1}) \\
& \qquad + \tau \left\| \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right\|_{L^1}^h (\dot{B}^{\frac{d}{2}-1}) \\
& \qquad + \tau \sum_{j \ge J_\tau - 1} 2^{\frac{d}{2}-1} j \| (\delta R_{1,j}, \delta R_{2,j}, \delta R_{3,j}) \|_{L^2}.
\end{align*}
$$
By (1.9), (2.4), (2.5), (5.9) and (5.15), it is easy to show

$$
\tau \left\| \left( \frac{1}{\rho^{\varepsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right\|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^h \\
\lesssim \| \delta \rho \|_{\tilde{L}_t^{\infty} (\dot{B}^{\frac{d}{2}-1})} \tau \| P^{\tau} - \bar{P} \|_{L_t^1 (\dot{B}^{\frac{d}{2}+1})} \lesssim o(1) \| (\delta Y, \delta Q) \|_{\tilde{L}_t^{\infty} (\dot{B}^{\frac{d}{2}-1})} + \sqrt{\varepsilon \tau},
$$

and

$$
\tau \sum_{j \ge J_{\tau}-1} 2^{\frac{d}{2}-1} j \| (\delta R_{1,j}, \delta R_{2,j}) \|_{L^2} \\
\lesssim \| \nabla u^{\varepsilon, \tau} \|_{L_t^1 (\dot{B}^{\frac{d}{2}})} \| (\delta Q, \delta u) \|_{\tilde{L}_t^{\infty} (\dot{B}^{\frac{d}{2}-1})} \lesssim o(1) \| (\delta Q, \delta u) \|_{\tilde{L}_t^{\infty} (\dot{B}^{\frac{d}{2}-1})}.
$$

For the tricky commutator term $R_{3,j}$, we have

$$
\tau \sum_{j \ge J_{\tau}-1} 2^{\frac{d}{2}-1} j \| \mathbf{R}_{3,j} \|_{L_t^1(L^2)} \\
\lesssim \tau \sum_{j \ge J_{\tau}-1} 2^{\frac{d}{2}j} \left\| \left[ \frac{1}{\rho^{\varepsilon, \tau}}, \dot{\Delta}_j \right] \nabla \delta Q^\ell \right\|_{L_t^1(L^2)} + \sum_{j \ge J_{\tau}-1} 2^{\frac{d}{2}-1} j \left\| \left[ \frac{1}{\rho^{\varepsilon, \tau}}, \dot{\Delta}_j \right] \nabla \delta Q^h \right\|_{L_t^1(L^2)} \\
\lesssim \| \nabla \rho^{\varepsilon, \tau} \|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}})} \left( \tau \| \delta Q^\ell \|_{L_t^1(\dot{B}^{\frac{d}{2}})} + \| \delta Q^h \|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \right) \\
\lesssim o(1) \left( \tau \| \delta Q \|_{L_t^1(\dot{B}^{\frac{d}{2}})}^\ell + \| \delta Q \|_{L_t^1(\dot{B}^{\frac{d}{2}-1})}^h \right).
$$

For $\delta F_1$ and $\delta F_2$, similar computations give rise to

$$
\begin{align}
\| (\delta F_1, \delta F_2) \|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} &\lesssim \| \delta u \|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \| (P^\tau - \bar{P}, u^\tau) \|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}+1})} \nonumber \\
&\quad + \| \Gamma_2^{\varepsilon, \tau} - \frac{\gamma+\gamma-P^\tau}{\gamma_+\alpha_-^\tau + \gamma_-\alpha_+^\tau} \|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})} \| u^\tau \|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} \nonumber \\
&\quad + \left( \| u^{\varepsilon, \tau} \|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} + \| \partial_t \Gamma_1^{\varepsilon, \tau} + u^{\varepsilon, \tau} \nabla \Gamma_1^{\varepsilon, \tau} \|_{L_t^1(\dot{B}^{\frac{d}{2}})} \right) \nonumber \\
&\quad \times \| P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau}) \|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})} \nonumber \\
&\quad + \| P_+ (\rho_+^{\varepsilon, \tau}) - P_- (\rho_-^{\varepsilon, \tau}) \|_{L_t^1(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})} \nonumber \\
&\lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1})} + \| \delta u \|_{L_t^1(\dot{B}^{\frac{d}{2}-1})} \right) + \sqrt{\varepsilon \tau} + \varepsilon. \tag{5.28}
\end{align}
$$
We thus get

$$ \begin{aligned} & \| (\delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})}^h + \tau \| (\delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})}^h \\ & + \| (\delta Q, \delta u) \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^h + \sqrt{\tau} \| \delta Q \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^h + \frac{1}{\sqrt{\tau}} \| u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-2})}^h \\ & \lesssim o(1) \left( \| (\delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})} + \tau \| \delta Q \|_{L_t^1 (\dot{B}^{\frac{d}{2}})}^\ell \right. \\ & \qquad \left. + \| \delta Q \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^h + \tau \| \delta u \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})} \right) + \sqrt{\epsilon \tau}. \end{aligned} \quad (5.29) $$

### Step 3: $\dot{B}^{\frac{d}{2}-1}$-estimates of $(\delta Q, \delta u)$ in all frequencies

We need to further establish the uniform $\dot{B}^{\frac{d}{2}-1}$-bounds. To this end, owing to (5.26), we obtain the $L^2$-in-time estimate

$$ \begin{aligned} & \| (\delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})} \\ & \lesssim \sqrt{\epsilon \tau} + \left( \| u^{\epsilon, \tau} \|_{L_t^1 (\dot{B}^{\frac{d}{2}})} + \| (\partial_t \rho^{\epsilon, \tau}, \partial_t \Gamma_2^{\epsilon, \tau}) \|_{L_t^1 (\dot{B}^{\frac{d}{2}})} \right)^{\frac{1}{2}} \| (\delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})}^h \\ & + \| (\rho^{\epsilon, \tau} - \bar{\rho}, \Gamma_2^{\epsilon, \tau} - \bar{\Gamma}_2) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}+1})}^{\frac{1}{2}} \left( \sqrt{\tau} \| \delta Q \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^h \right)^{\frac{1}{2}} \left( \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^h \right)^{\frac{1}{2}} \\ & + \| (\delta F_1, \delta F_2) \|_{L_t^1 (\dot{B}^{\frac{d}{2}-1})}^{\frac{1}{2}} \| (\delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})}^{\frac{1}{2}} \\ & + \left\| \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right\|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^{\frac{1}{2}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^{\frac{1}{2}} \\ & + \sum_{j \in \mathbb{Z}} 2^{\frac{d}{2}-1} j \left( \int_0^t \left( \| \delta R_{1,j} \|_{L^2} \| \dot{\Delta}_j \delta Q \|_{L^2} + \| (\delta R_{2,j}, \delta R_{3,j}) \|_{L^2} \| \dot{\Delta}_j u \|_{L^2} \right) ds \right)^{\frac{1}{2}}. \end{aligned} $$

One has

$$ \begin{aligned} & \left\| \left( \frac{1}{\rho^{\epsilon, \tau}} - \frac{1}{\rho^{\tau}} \right) \nabla P^{\tau} \right\|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^{\frac{1}{2}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})}^{\frac{1}{2}} \\ & \lesssim \| \delta \rho \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})} \left( \sqrt{\tau} \| P^{\tau} - \bar{P} \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}+1})} \right)^{\frac{1}{2}} \left( \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})} \right)^{\frac{1}{2}} \\ & \lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\tilde{L}_t^\infty (\dot{B}^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}^{\frac{d}{2}-1})} \right) + \sqrt{\epsilon \tau}. \end{aligned} $$
Concerning the commutator terms, we have

$$
\begin{aligned}
& \sum_{j \in \mathbb{Z}} 2^{\left(\frac{d}{2}-1\right)j} \left( \int_0^t \left( \|\delta R_{1,j}\|_{L^2} \|\dot{\Delta}_j \delta Q\|_{L^2} + \|\delta R_{2,j}\|_{L^2} \|\dot{\Delta}_j u\|_{L^2} \right) \mathrm{d}s \right)^{\frac{1}{2}} \\
& \lesssim \left( \frac{1}{\sqrt{\tau}} \|\nabla u^{\varepsilon, \tau}\|_{\tilde{L}_t^2(\dot{B}_2^{\frac{d}{2}-1})} \|\delta Q\|_{\tilde{L}_t^{\infty}(\dot{B}_2^{\frac{d}{2}-1})} \sqrt{\tau} \|\delta Q\|_{\tilde{L}_t^2(\dot{B}_2^{\frac{d}{2}-1})} \right)^{\frac{1}{2}} \\
& \quad + \left( \frac{1}{\sqrt{\tau}} \|\delta u\|_{\tilde{L}_t^2(\dot{B}_2^{\frac{d}{2}-1})} \|\nabla \Gamma_2^{\varepsilon, \tau}\|_{\tilde{L}_t^{\infty}(\dot{B}_2^{\frac{d}{2}+1})} \sqrt{\tau} \|\delta Q\|_{\tilde{L}_t^2(\dot{B}_2^{\frac{d}{2}-1})} \right)^{\frac{1}{2}} \\
& \quad + \left( \|\nabla u^{\varepsilon, \tau}\|_{L_t^1(\dot{B}_2^{\frac{d}{2}-1})} \|\delta u\|_{\tilde{L}_t^{\infty}(\dot{B}_2^{\frac{d}{2}-1})}^2 \right)^{\frac{1}{2}} \\
& \quad + \left( \|\nabla \rho^{\varepsilon, \tau}\|_{\tilde{L}_t^{\infty}(\dot{B}_2^{\frac{d}{2}-1})} \sqrt{\tau} \|\delta Q\|_{\tilde{L}_t^{\infty}(\dot{B}_2^{\frac{d}{2}-1})} \frac{1}{\sqrt{\tau}} \|\delta u\|_{\tilde{L}_t^2(\dot{B}_2^{\frac{d}{2}-1})} \right)^{\frac{1}{2}} \\
& \lesssim o(1) \left( \| (\delta Q, \delta u) \|_{\tilde{L}_t^{\infty}(\dot{B}_2^{\frac{d}{2}-1})} + \sqrt{\tau} \|\delta Q\|_{\tilde{L}_t^2(\dot{B}_2^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\tau}} \|\delta u\|_{\tilde{L}_t^2(\dot{B}_2^{\frac{d}{2}-1})} \right).
\end{aligned}
$$

Gathering (5.28) and the above three estimates, we have

$$
\begin{aligned}
& \| (\delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}_2^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}_2^{\frac{d}{2}-1})} \\
& \lesssim o(1) \left( \| (\delta Y, \delta Q, \delta u) \|_{\tilde{L}_t^\infty (\dot{B}_2^{\frac{d}{2}-1})} + \| \delta Q \|_{\tilde{L}_t^2 (\dot{B}_2^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}_2^{\frac{d}{2}-1})} \right) + \sqrt{\varepsilon \tau}.
\end{aligned}
\tag{5.30}
$$

• Step 4: Proof of convergence rate

Finally, as required in (5.25), one needs to estimate $\partial_t \delta Q$. We make use of the equation (5.17)$_1$, (5.21) and (5.22) to get

$$
\begin{aligned}
\| \partial_t \delta Q \|_{L_t^1 (\dot{B}_2^{\frac{d}{2}-2})} & \lesssim \| u^{\varepsilon, \tau} \cdot \nabla \delta Q + (\Gamma_2^{\varepsilon, \tau} - \bar{\Gamma}_2) \operatorname{div} \delta u \|_{L_t^1 (\dot{B}_2^{\frac{d}{2}-2})} \\
& \quad + \| \operatorname{div} \delta u \|_{L_t^1 (\dot{B}_2^{\frac{d}{2}-2})} + \| \delta F_1 \|_{L_t^1 (\dot{B}_2^{\frac{d}{2}-2})} \\
& \lesssim o(1) \left( \| (\delta Y, \delta Q) \|_{\tilde{L}_t^\infty (\dot{B}_2^{\frac{d}{2}-1})} + \frac{1}{\sqrt{\tau}} \| \delta u \|_{\tilde{L}_t^2 (\dot{B}_2^{\frac{d}{2}-1})} \right) \\
& \quad + \| \delta u \|_{L_t^1 (\dot{B}_2^{\frac{d}{2}-1})} + \sqrt{\varepsilon \tau}.
\end{aligned}
\tag{5.31}
$$

Combining (5.4), (5.25), (5.29), (5.30) and (5.31) together, we end up with (5.18) which completes the proof of Lemma 5.4. □
## 5.2. Time-relaxation limit: (Kτ) to (PM)

This section is devoted to the proof of (1.19) in Theorem 1.4. Define the error variables

$$ (\delta\beta_{\pm}, \delta\varrho_{\pm}, \delta\varrho, \delta\Pi, \delta v) := (\beta_{\pm}^{\tau} - \beta_{\pm}, \varrho_{\pm}^{\tau} - \varrho_{\pm}, \varrho^{\tau} - \varrho, \Pi^{\tau} - \Pi, v^{\tau} - v). $$

First, similarly to (5.1), instead of $\delta\beta$, we need to estimate the variable

$$ \delta Z := \frac{\beta_{+}^{\tau} \varrho_{+}^{\tau}}{\varrho^{\tau}} - \frac{\beta_{+} \varrho_{+}}{\varrho}, $$

where the initial data of $\delta z$ is

$$ \delta Z|_{t=0} = Z_0 - Z_0, \quad Z_0 := \frac{\alpha_{+,0}^{\tau} \rho_{+,0}^{\tau}}{\alpha_{+,0}^{\tau} \rho_{+,0}^{\tau} + \alpha_{-,0}^{\tau} \rho_{-,0}^{\tau}}, \quad Z_0 := \frac{\beta_{+,0} \varrho_{+,0}}{\beta_{+,0} \varrho_{+,0} + \beta_{-,0} \varrho_{-,0}}. \quad (5.32) $$

Indeed, arguing similarly as in Lemma 5.2, we obtain from (Kτ) and (PM) that

$$ \left\{ \begin{aligned} \delta\beta_+ &= \frac{1}{\beta_-^\tau \varrho_+^\tau \varrho_- + \beta_+ \varrho_+ \varrho_-^\tau} (\varrho^\tau \varrho \delta Z - \beta_+ \beta_- \varrho_- \delta \varrho_+ + \beta_+^\tau \beta_- \varrho_+ \delta \varrho_-), \\ \delta\varrho^\tau &= \frac{\varrho_+^\tau - \varrho_-^\tau}{\beta_- \varrho_+^\tau \varrho_- + \beta_+ \varrho_+ \varrho_-^\tau} (\varrho^\tau \varrho \delta Z - \beta_+ \beta_- \varrho_- \delta \varrho_+ + \beta_+ \beta_- \varrho_+ \delta \varrho_-) \\ &\quad + \beta_+ \delta \varrho_+ + \beta_- \delta \varrho_-^\tau, \\ \delta\Pi &= \delta \varrho_+ \int_0^1 P'_+ (\theta \varrho_+^\tau + (1-\theta) \varrho_+) d\theta = \delta \varrho_- \int_0^1 P'_- (\theta \varrho_-^\tau + (1-\theta) \varrho_-) d\theta, \end{aligned} \right. \quad (5.33) $$

which leads to

$$ \left\{ \begin{aligned} \|\delta\varrho_{\pm}\|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2}-1})} + \|\delta\varrho_{\pm}\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} &\sim \|\delta\Pi\|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2}-1})} + \|\delta\Pi\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})}, \\ \|\delta\beta_{\pm}\|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2}-1})} &\lesssim \|(\delta Z, \delta \Pi)\|_{\tilde{L}_t^{\infty}(\dot{B}^{\frac{d}{2}-1})}. \end{aligned} \right. \quad (5.34) $$

It is therefore sufficient to estimate $(\delta\Pi, \delta v, \delta Z)$ to recover the information on all the error unknowns.

Next, note that $\delta Z$ satisfies the transport equation

$$ \partial_t \delta Z + v^\tau \cdot \nabla \delta Z = -\delta v \cdot \nabla \frac{\beta_+ \varrho_+}{\varrho}. \quad (5.35) $$

Using (2.4), (4.22) and Lemma 2.7, we get

$$ \begin{aligned} \|\delta Z\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2}-1} \cap \dot{B}^{\frac{d}{2}})} &\lesssim \exp \left( \|v^\tau\|_{L_t^1(\dot{B}^{\frac{d}{2}+1})} \right) \|\delta v\|_{L_t^1(\dot{B}^{\frac{d}{2}})} \|\nabla \frac{\beta_+ \varrho_+}{\varrho}\|_{\tilde{L}_t^\infty(\dot{B}^{\frac{d}{2} \cap \dot{B}^{\frac{d}{2}+1}})} \\ &\lesssim o(1) \|\delta v\|_{L_t^1(\dot{B}^{\frac{d}{2}})}. \end{aligned} \quad (5.36) $$

Then, we perform the key estimates of $\delta\Pi$. From (Kτ), it is easy to see

$$ \begin{aligned} \partial_t \Pi^\tau + v^\tau \cdot \nabla \Pi^\tau &= \frac{\gamma_+ \gamma_- \Pi^\tau}{\gamma_+ \beta_-^\tau + \gamma_- \beta_+^\tau} \operatorname{div} \left( \frac{\nabla \Pi^\tau}{\varrho^\tau} \right) - \frac{\gamma_+ \gamma_- \Pi^\tau}{\gamma_+ \beta_-^\tau + \gamma_- \beta_+^\tau} \operatorname{div} z^\tau, \\ z^\tau &:= v^\tau + \frac{\nabla \Pi^\tau}{\varrho^\tau}. \end{aligned} $$
Thence by the above equation and (PM), $\delta\Pi$ satisfies

$$ 
\begin{aligned}
\partial_t \delta\Pi - \bar{c} \Delta \delta\Pi = & -v^\tau \cdot \nabla \delta\Pi - \delta v \cdot \nabla \Pi \\
& + \left( \frac{\gamma_+ \gamma_- \Pi^\tau}{\gamma_+ \beta_+^\tau + \gamma_- \beta_-^\tau} - \frac{\gamma_+ \gamma_- \Pi}{\gamma_+ \beta_+ + \gamma_- \beta_-} \right) \operatorname{div} \left( \frac{\nabla \Pi^\tau}{\rho^\tau} \right) \\
& + \frac{\gamma_+ \gamma_- \Pi}{\gamma_+ \beta_+ + \gamma_- \beta_-} \operatorname{div} \left( \left( \frac{1}{\rho^\tau} - \frac{1}{\rho} \right) \nabla \Pi^\tau \right) - \frac{\gamma_+ \gamma_- \Pi^\tau}{\gamma_+ \beta_-^\tau + \gamma_- \beta_+^\tau} \operatorname{div} z^\tau,
\end{aligned} 
\quad (5.37) $$

with the constant $\bar{c} := \frac{\gamma_+ \gamma_- P_+(\bar{\rho}_+)}{(\gamma_+ \bar{\alpha}_- + \gamma_- \bar{\alpha}_+) \bar{\rho}} > 0$. We mention that the convergence rate $\tau$ is $O(\tau)$ bound comes from the uniform estimate (4.22) of the effective unknown $z^\tau$. Indeed, by Lemma 2.6, the uniform estimate (4.22), the smallness of the initial data (1.7), the product laws (2.3) and the composition estimates (2.6)–(2.7), one obtains

$$ 
\begin{aligned}
\|\delta\Pi\|_{\tilde{L}_t^\infty(\dot{B}^{d-1}_2)} + \|\delta\Pi\|_{L_t^1(\dot{B}^{d+1}_2)} &\lesssim \tau + \|v^\tau\|_{\tilde{L}_t^2(\dot{B}^d_2)} \|\delta\Pi\|_{\tilde{L}_t^2(\dot{B}^{d-1}_2)} \\
&+ \|\delta v\|_{L_t^1(\dot{B}^d_2)} \|\Pi - \bar{P}\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} + \|\delta\Pi\|_{\tilde{L}_t^2(\dot{B}^d_2)} \|\Pi^\tau - \bar{P}\|_{\tilde{L}_t^2(\dot{B}^{d+1}_2)} \\
&+ \left\| \left(\frac{1}{\rho^\tau} - \frac{1}{\rho}\right) \nabla\Pi^\tau \right\|_{L_t^1(\dot{B}^d_2)} + \|z^\tau\|_{L_t^1(\dot{B}^d_2)} \\
&\lesssim o(1) \left( \|\delta Z\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} + \| (\delta\varrho_\pm, \delta\Pi) \|_{\tilde{L}_t^\infty(\dot{B}^{d-1}_2)} + \| (\delta\varrho_\pm, \delta\Pi) \|_{L_t^1(\dot{B}^{d+1}_2)} \right) + \tau,
\end{aligned} 
\quad (5.38) $$

where we have used the key fact

$$ 
\begin{aligned}
\left\| \left( \frac{1}{\rho^\tau} - \frac{1}{\rho} \right) \nabla \Pi^\tau \right\|_{L_t^1(\dot{B}^d_2)} &\lesssim \|\delta Z\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} \|\Pi^\tau - \bar{P}\|_{L_t^1(\dot{B}^{d+1}_2)} \\
&+ \|\delta \varrho_\pm\|_{\tilde{L}_t^2(\dot{B}^d_2)} \|\Pi^\tau - \bar{P}\|_{\tilde{L}_t^2(\dot{B}^{d+1}_2)} \\
&\lesssim o(1) \left( \|\delta Z\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} + \|\delta \varrho_\pm\|_{\tilde{L}_t^2(\dot{B}^d_2)} \right),
\end{aligned} 
\quad (5.39) $$

derived from (2.4), (4.22) and (5.33). Gathering (5.34) and (5.38) together, we get

$$ \|\delta\Pi\|_{\tilde{L}_t^\infty(\dot{B}^{d-1}_2)} + \|\delta\Pi\|_{L_t^1(\dot{B}^{d+1}_2)} \lesssim o(1) \left( \|\delta Z\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} + \|\delta\Pi\|_{L_t^1(\dot{B}^{d+1}_2)} \right) + \tau. \quad (5.40) $$

For the error unknown $\delta v$, in view of (4.22), (5.39) and $\delta v = (\frac{1}{\rho^\tau} - \frac{1}{\rho})\nabla\Pi^\tau - \frac{1}{\rho}\nabla\delta\Pi + z^\tau$, it can be bounded by

$$ 
\begin{aligned}
\|\delta v\|_{L_t^1(\dot{B}^d_2)} &\lesssim \|\delta\Pi\|_{L_t^1(\dot{B}^{d+1}_2)} + \left\| \left(\frac{1}{\rho^\tau} - \frac{1}{\rho}\right) \nabla\Pi^\tau \right\|_{L_t^1(\dot{B}^d_2)} + \|z^\tau\|_{L_t^1(\dot{B}^d_2)} \\
&\lesssim \|\delta Z\|_{\tilde{L}_t^\infty(\dot{B}^d_2)} + \|\delta\Pi\|_{\tilde{L}_t^\infty(\dot{B}^{d-1}_2)} + \|\delta\Pi\|_{L_t^1(\dot{B}^{d+1}_2)} + \tau.
\end{aligned} 
\quad (5.41) $$

The combination of the estimate (5.34) and the inequalities (5.36)–(5.40) gives rise to (1.19), which completes the proof of Theorem 1.4.
**Data availability statement**

No new data were created or analysed in this study.

**Acknowledgments**

The authors are indebted to the anonymous referees for their valuable suggestions and com-
ments on the manuscript. T Crin-Barat has been funded by the Alexander von Humboldt
Professorship program and the Transregio 154 Project Mathematical Modelling, Simulation
and Optimization Using the Example of Gas Networks of the DFG. L-Y Shou is supported
by the National Natural Science Foundation of China (12301275) and the China Postdoctoral
Science Foundation (2023M741694). J Tan is partially supported by the ANR Project BORDS
(ANR-16-CE40-0027-01), by the Labex MMEDII and the CY Initiative of Excellence, Project
CYNA (CY Nonlinear Analysis).

**References**

[1] Baer M R and Nunziato J W 1986 A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials *Int. J. Multiph. Flow* **12** 861–89

[2] Bahouri H, Chemin J-Y and Danchin R 2011 *Fourier Analysis and Nonlinear Partial Differential Equations* (Grundlehren der Mathematischen Wissenschaften vol 343) (Springer)

[3] Beauchard K and Zuazua E 2011 Large time asymptotics for partially dissipative hyperbolic systems *Arch. Ration. Mech. Anal.* **199** 177–227

[4] Benzoni-Gavage S and Serre D 2007 *Multidimensional Hyperbolic Partial Differential Equations: First-Order Systems and Applications (Oxford Mathematical Monographs)* (The Clarendon, Oxford University Press)

[5] Bianchini R, Crin-Barat T and Paicu M 2024 Relaxation approximation and asymptotic stability of stratified solutions to the IPM equation *Arch. Ration. Mech. Anal.* **248** 2

[6] Bresch D, Burtea C and Lagoutiére F 2020 Mathematical justification of a compressible bifluid system with different pressure laws: a continuous approach *Appl. Anal.* **101** 4235–66

[7] Bresch D, Burtea C and Lagoutiére F 2023 Mathematical justification of a compressible bi-fluid system with different pressures laws: a semi-discrete approach and numerical illustrations *J. Comput. Phys.* **490** 112259

[8] Bresch D, Desjardins B, Ghidaglia J M, Grenier E and Hilliairet M 2018 Multifluid models including compressible fluids *Handbook of Mathematical Analysis in Mechanics of Viscous Fluids* ed G Y Novotný (Springer)

[9] Bresch D, Desjardins B, Ghidaglia J-M and Grenier E 2010 Global weak solutions to a generic two-fluid model *Arch. Ration. Mech. Anal.* **196** 599–629

[10] Bresch D and Hillairet M 2015 Note on the derivation of multicomponent flow systems *Proc. Am. Math. Soc.* **143** 3429–43

[11] Bresch D and Hillairet M 2019 A compressible multifluid system with new physical relaxation terms *Ann. ENS* **52** 255–95

[12] Bresch D and Huang X 2011 A multi-fluid compressible system as the limit of weak solutions of the isentropic compressible Navier–Stokes equations *Arch. Ration. Mech. Anal.* **201** 647–80

[13] Bresch D, Huang X and Li J 2012 Global weak solutions to one-dimensional non-conservative viscous compressible two-phase system *Commun. Math. Phys.* **309** 737–55

[14] Bresch D, Mucha P-B and Zatorska E 2019 Finite-energy solutions for compressible two-fluid Stokes system *Arch. Ration. Mech. Anal.* **232** 987–1029

[15] Burtea C, Crin-Barat T and Tan J 2023 Pressure-relaxation limit for a damped one-velocity Baer–Nunziato model to a Kappila model *Math. Models Methods Appl. Sci.* **33** 687–753

[16] Chen G-Q, Levermore C D and Liu T-P 1994 Hyperbolic conservation laws with stiff relaxation terms and entropy *Commun. Pure Appl. Math.* **47** 787–830
[17] Coulombel J-F and Goudon T 2007 The strong relaxation limit of the multidimensional isothermal Euler equations *Trans. Am. Math. Soc.* **359** 637–48

[18] Crin-Barat T and Danchin R 2022 Partially dissipative hyperbolic systems in the critical regularity setting: the multi-dimensional case *J. Math. Pures Appl.* **165** 1–41

[19] Crin-Barat T and Danchin R 2023 Global existence for partially dissipative hyperbolic systems in the $L^p$ framework and relaxation limit *Math. Ann.* **386** 2159–206

[20] Crin-Barat T, He Q and Shou L-Y 2023 The hyperbolic-parabolic chemotaxis system for vasculogenesis: global dynamics and relaxation limit toward a Keller–Segel model *SIAM J. Math. Anal.* **55** 4445–92

[21] Crin-Barat T and Shou L-Y 2023 Diffusive relaxation limit of the multi-dimensional Jin-Xin system *J. Differ. Equ.* **357** 302–31

[22] Danchin R 2000 Global existence in critical spaces for compressible Navier–Stokes equations *Invent. Math.* **141** 579–614

[23] Danchin R 2018 Fourier analysis methods for the compressible Navier–Stokes equations *Handbook of Mathematical Analysis in Mechanics of Viscous Fluids* ed G Y Novotný (Springer)

[24] Danchin R 2023 Partially dissipative systems in the critical regularity setting and strong relaxation limit *EMS Surv. Math. Sci.* **9** 135–92

[25] Evje S and Karlsen K H 2008 Global existence of weak solutions for a viscous two-phase model *J. Differ. Equ.* **245** 2660–703

[26] Evje S, Wang W and Wen H 2016 Global well-posedness and decay rates of strong solutions to a non-conservative compressible two-fluid model *Arch. Ration. Mech. Anal.* **221** 1285–316

[27] Forestier A and Gavrilyuk S 2011 Criterion of hyperbolicity for non-conservative quasilinear systems admitting a partially convex conservation law *Math. Methods Appl. Sci.* **34** 2148–58

[28] Giovangigli V and Yong W-A 2014 Volume viscosity and internal energy relaxation: symmetrization and Chapman–Enskog expansion *Kinet. Relat. Models* **8** 79–116

[29] Giovangigli V and Yong W-A 2018 Volume viscosity and internal energy relaxation: error estimates *Nonlinear Anal. Real World Appl.* **43** 213–44

[30] Guo Z, Yang J and Yao L 2011 Global strong solution for a three-dimensional viscous liquid-gas two-phase flow model with vacuum *J. Math. Phys.* **52** 093102

[31] Hao C and Li H-L 2012 Well-posedness for a multidimensional viscous liquid-gas two-phase flow model *SIAM J. Math. Anal.* **44** 1304–32

[32] Ishii M 2006 *Thermo-Fluid Dynamics of Two-Phase Flow* (Springer)

[33] Junca S and Rascle M 2002 Strong relaxation of the isothermal Euler system to the heat equation *Z. Angew. Math. Phys.* **53** 239–64

[34] Kapila A, Menikoff R, Bdzil J, Son S and Stewart D 2001 Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations *Phys. Fluids* **13** 3002–24

[35] Kraccčmar S, Kwon Y-S, Neccčasová Š and Novotný A 2022 Weak solutions for a bifluid model for a mixture of two compressible noninteracting fluids with general boundary data *SIAM J. Math. Anal.* **54** 818–71

[36] Li H-L and Shou L-Y 2023 Global existence and optimal time-decay rates of the compressible Navier–Stokes–Euler system *SIAM J. Math. Anal.* **55** 1810–46

[37] Li H-L and Shou L-Y 2023 Global existence of weak solutions to the drift-flux system for general pressure laws *Sci. China Math.* **66** 251–84

[38] Li J, Yu Y and Zhu C 2024 Ill-posedness for the Burgers equation in Sobolev spaces *Indian J. Pure Appl. Math.* **55** 189–97

[39] Linares F, Pilod D and Saut J-C 2014 Dispersive perturbations of Burgers and hyperbolic equations: local theory *SIAM J. Math. Anal.* **46** 1505–37

[40] Majda A 1984 *Compressible Fluid Flow and Systems of Conservation Laws in Several Space Variable* (Springer)

[41] Matsumura A and Nishida T 1979 The Cauchy problem for the equations of motion of compressible viscous and heat-conductive fluids *Proc. Jpn. Acad.* **A** **55** 337–42

[42] Novotný A and Pokorný M 2020 Weak solutions for some compressible multicomponent fluid models *Arch. Ration. Mech. Anal.* **235** 355–403

[43] Runst T and Sickel W 1996 *Sobolev Spaces of Fractional Order, Nemytskij Operators and Nonlinear Partial Differential Equations (Nonlinear Analysis and Applications)* (Walter de Gruyter & Co.)

[44] Shizuta S and Kawashima S 1985 Systems of equations of hyperbolic-parabolic type with applications to the discrete Boltzmann equation *Hokkaido Math. J.* **14** 249–75
[45] Vasseur A, Wen H and Yu C 2019 Global weak solution to the viscous two-fluid model with finite energy J. Math. Pures Appl. 125 247-82

[46] Wallis G 1979 *One-Dimensional Two-Fluid Flow* (McGraw-Hill)

[47] Wen H 2021 On global solutions to a viscous compressible two-fluid model with unconstrained transition to single-phase flow in three dimensions Calc. Var. PDE 60 158

[48] Wen H, Yao L and Zhu C 2018 Review on mathematical analysis of some two-phase flow models Acta Math. Sci. 38 1617-36

[49] Xu J and Kawashima S 2014 Global classical solutions for partially dissipative hyperbolic system of balance laws Arch. Ration. Mech. Anal. 211 513-53

[50] Xu J and Wang Z 2013 Relaxation limit in Besov spaces for compressible Euler equations J. Math. Pures Appl. 99 43-61

[51] Yao L, Zhang T and Zhu C 2010 Existence of asymptotic behavior of global weak solutions to a 2D viscous liquid-gas two-phase flow model SIAM J. Math. Anal. 42 1874-97

[52] Yao L and Zhu C 2011 Existence and uniqueness of global weak solution to a two-phase flow model with vacuum Math. Ann. 349 903-28

[53] Zhang Y 2016 Decay of the 3D inviscid liquid-gas two-phase flow model Z. Angew. Math. Phys. 67 54

[54] Zhang Y and Zhu C 2015 Global existence and optimal convergence rates for the strong solutions in $H^2$ to the 3d viscous liquid-gas two-phase flow model J. Differ. Equ. 258 2315-38

[55] Zuazua E 2022 Decay of partially dissipative hyperbolic systems (available at: www.dm.uniba.it/it/ricerca/convegni/2022/acipdif22/decaypartiallydiss-1dfocus.pdf/view)