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

In this paper we investigate two types of relaxation processes quantitatively in the context of small data global-in-time solutions for compressible one-velocity multi-fluid models. First, we justify the pressure-relaxation limit from a one-velocity Baer-Nunziato system to a Kapila model as the pressure-relaxation parameter tends to zero, in a uniform manner with respect to the time-relaxation parameter associated to the friction forces modeled in the equation of the velocity. This uniformity allows us to further consider the time-relaxation limit for the Kapila model. More precisely, we show that the diffusely time-rescaled solution of the Kapila system converges to the solution of a two-phase porous media type system as the time-relaxation parameter tends to zero. For both relaxation limits, we exhibit explicit convergence rates. Our proof of existence results are based on an elaborate low-frequency and high-frequency analysis via the Littlewood-Paley decomposition and it includes three main ingredients: a refined spectral analysis for the linearized problem to determine the threshold of frequencies 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 the renormalized energy estimates in the high-frequency region to cancel higher-order nonlinear terms. To show the convergence rates, we discover several auxiliary unknowns that reveal better structures. In conclusion, our approach may be applied to a class of non-symmetric partially dissipative hyperbolic system with rough coefficients which do not have any time-integrability, in the context of overdamping phenomenon. It extends the latest results of Danchin and the first author [15, 16].


Models and motivations
Multi-phase flows have been used to simulate a wide range of physical mixing phenomenon, from engineering to biological systems (cf.[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: where the unknowns α ± = α ± (t, x) ∈ [0, 1], ρ ± = ρ ± (t, x) ≥ 0 and u = u(t, x) ∈ R d stand for the volume fractions, the densities and the common velocity of two fluids (denoted by + and −), respectively, which satisfy The two positive constants ε and τ 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 ± (s) = A ± s γ± with constants A ± > 0, 1 ≤ γ − < γ + . (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 onedimensional 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 multifluid 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.
Then, we further investigate the time-relaxation limit of System (K) as τ → 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/τ ) and define the following charge of variables Under the diffusive scaling (1.2), System (K) becomes As τ → 0, one then expects that (β τ ± , ̺ τ ± , v τ ) converges to some limit (β ± , ̺ ± , v) which is the solution of a new two-phase system with ̺ = β + ̺ + + β − ̺ − and Π = P + (̺ + ) = P − (̺ − ).Inserting Darcy's law (1.3) 2 into (1.3) 1 , we derive the following two-phase system in porous media: (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 System (BN) to System (K).In [15], the smallness condition on initial data employed to justify their global well-posedness result depends on min{τ, 1 τ } (due to the overdamping phenomenon that will be explained below) and therefore does not allow to further investigate the limit when τ → 0.
The main results of this article are the quantitative justification of the pressure-relaxation limit from System (BN) to System (K) as ε → 0 uniformly in τ and the time-relaxation limit from System (K τ ) to System (PM) as τ → 0. Consequently, a new two-phase flow system in porous media (PM) is rigorously derived from Systems (K τ ) and (BN τ ), which implies that Kapila and Baer-Nunziato systems considered in our paper can be viewed, for ε and τ small enough, as hyperbolic approximations of (PM).
On the other hand, it is natural to ask what happens for System (BN) as τ tends to 0 first.To investigate this process, we introduce a diffusive scaling similar to (1.2) as follows Under such scaling (1.6), System (BN) becomes − and Π ε,τ = β ε,τ + P + (̺ ε,τ + ) + β ε,τ − P − (̺ ε,τ − ).The crucial observation is that the parameter τ now also appears under the pressure-relaxation term in the equation of the volume fractions.This reveals that as τ → 0, the two pressures in System (BN τ ) should converge to a common pressure, and the solutions of System (BN) should converge to the solutions of 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 situation 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-3.2for 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.The Figure 1 summarizes the limit processes that we tackle in this article.

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 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 System (BN) results of solutions for Systems (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 . We agree that C b ([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 + ; X) and for the associated norm.We will keep the same notations for multi-component functions, namely for f : [0, T ] → X m with m ∈ N. F and F −1 stand for the Fourier transform and its inverse, and define the operator Λ σ := F −1 (|ξ| σ F (•)).Finally, let o(1) denote a generic constant which can be sufficiently small.

Main results
Our first theorem concerns the uniform, in both relaxation parameters ε and τ , global well-posedness of System (BN) in a critical regularity framework.
. There exists a positive constant c 0 independent of τ and ε such that if then the Cauchy problem of System (BN) with the initial data (α ±,0 , ρ ±,0 , u 0 ) has a unique global solution ), ), ). (1.8) Moreover, the following uniform estimate holds: ) where C > 0 is a generic constant.
Remark 1.1.It should be emphasized that the regularity and decay-in-τ properties of the effective flux ρ ε,τ u ε,τ τ + ∇P ε,τ is better than the one verified by the solution (α ε,τ ± , ρ ε,τ ± , u ε,τ ).This is consistent with Darcy's law and plays key role in the justification of the time-relaxation limit.
Next, we present the rigorous justifications of the pressure-relaxation limit for System (BN) to System (K) as ε → 0 uniform with respect to τ , and further the time-relaxation limit for System (K τ ) to System (PM) as τ → 0, with explicit convergence rates.
) be defined by (1.6).Then under the assumptions of Theorem 1.4, there is a generic constant C 3 such that

Difficulties and strategies
The first difficulty concerning the study of System (BN) are 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 (cf.[44]).
Additionally, 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 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 and the corresponding initial data so that the Cauchy problem of System (BN) subject to the initial data (α ±,0 , ρ ±,0 , u 0 ) is reformulated as where F ε,τ i = F ε,τ i (y, w, r) (i = 0, 1, 2, 3, 4) are the nonlinear terms Fi (i = 0, 1, 2, 3) are the constants and .26)In this formulation, the equation (1.23) 1 is purely transport and the linear part of subsystem (1.23) 2 - (1.23) 4 is partially dissipative and satisfies the "Shizuta-Kawashima" stability condition.Thus, we will estimate the undamped unknown y ε,τ and the dissipative components (w ε,τ , r ε,τ , u ε,τ ) separately.We emphasize here that due to the double parameters ε, τ and the lack of time-integrability of G ε,τ i , the dissipative structures of subsystem (1.23) 2 -(1.23) 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 ε, τ .
• In the high-frequency region |ξ| >> 1 ε , by Taylor's expansion near |εξ| −1 << 1, the real eigenvalue λ 1 and the conjugated complex eigenvalues λ i (i = 2, 3) The above spectral analysis suggests us to separate the whole frequencies into two parts |ξ| 1 τ and |ξ| 1 τ so as to capture the qualitative properties of solutions for System (1.23).Indeed, the time-decay rates (determined by λ 2 ) achieve the fastest rate in the low-frequency region |ξ| 1 τ .Moreover this region recover the whole frequency-space when τ → 0, as expected from the well-known overdamping phenomenon which will be mentioned below.To this end, the threshold J τ between these two regions is used in the definition of the hybrid Besov spaces in next section.
It should be noted that λ 2 and λ 3 exhibit similar behaviors to the eigenvalues of the compressible Euler equations with damping.Indeed, to study System (1.23), one considers the following simplified system of damped Euler type with rough coefficients: (1.27) The well-known spectral analysis for the linear Euler part of System (1.27) implies that the frequency space shall be separated into the low-frequency region |ξ| 1 τ and the high-frequency region |ξ| 1 τ to recover the uniform estimates and optimal regularity of solutions.Formally, this implies that as τ → 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 1 τ gets larger, the decay rates of r τ do not necessarily increase and on the contrary follow min{τ, 1 τ }, 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 . 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 d 2 + 1 is called critical for initial data of general hyperbolic systems since Ḃ 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 < d 2 + 1.
Nevertheless, the methods developed in [18][19][20] are not applicable in the current situation to derive estimates which are uniform with respect to the relaxation parameter τ .This is mainly due to the complex form of the total pressure P ε,τ in the velocity equation (BN) 3 and the fact that one can not expect any time integrability property on R + for the purely transported unknown y ε,τ , which generally leads to a lack of time integrability on R + for (α ε,τ ± − ᾱ± , ρ ε,τ ± − ρ± ) (see Remarks 3.2-3.3),and thus for G τ 3 , G τ 0 in System (1.27).In addition, we can not find a a rescaling to reduce the proof to the case τ = 1 and then recover the corresponding uniform estimates with respect to τ 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 ε, τ 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) 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][19][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 -time dissipation properties in high frequencies (cf.Lemma 3.2).Moreover, we also establish the uniform L 2 -in-time estimates at Ḃ d 2 +1 -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 ε, τ satisfying 0 < ε ≤ τ < 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.

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 χ(ξ) with compact supported in B(0, 4  3 ) and χ(ξ) = 1 in B(0, 3  4 ) such that For any j ∈ Z, the homogeneous dyadic blocks ∆j and the low-frequency cut-off operator Ṡj are defined by ∆j u From now on, we use the shorthand notation With the help of these dyadic blocks, the homogeneous Besov space Ḃs for s ∈ R is defined by We denote the Chemin-Lerner type space L ̺ (0, T ; Ḃs ) for s ∈ R and T > 0: By the Minkowski inequality, it holds that where is the usual Lebesgue-Besov norm.In order to perform our analysis on the low and high frequencies regions, we set the threshold for suitable negative integer k (to be determined).Denote the following notations for p ∈ [1, ∞] and For any u ∈ S ′ h , we also define the low-frequency part u ℓ and the high-frequency part u h by It is easy to check for any 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-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.
Due to the Bernstein inequalities, the Besov spaces have many useful properties: Lemma 2.2.The following properties hold: • For any s ∈ R and q ≥ 2, we have the following continuous embeddings: • Ḃ d 2 is continuously embedded in the set of continuous functions decaying to 0 at infinity.
Then the space Ḃs1 ∩ Ḃs2 is a Banach space and satisfies weak compact and Fatou properties: If u k is a uniformly bounded sequence of Ḃs1 ∩ Ḃs2 , then an element u of Ḃs1 ∩ Ḃs2 and a subsequence u n k exist such that 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: (2.4) The following commutator estimates will used to control some nonlinearities in high frequencies.
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 Ḃ0 .
Moreover, we note that Therefore, applying (2.3), (2.6) and the embedding 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 (cf.[2][Page 157] for example).
We have the regularity estimates of the damped transport equation.Since it can be directly shown by the commutator estimates (2.5) and Grönwall's inequality as in [15,23], we omit the proof for brevity. ) where C > 0 is a constant independent of T and λ * .

Analysis of the linearized system
We now consider the linearized problem associated to (1.23), which reads where h i (i = 1, ..., 5) are given positive constants and 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 ε, τ , which improves the result in [15] without the uniformity with respect to τ .As explained before, the threshold J τ between low and high frequencies given by (2.1) is the key to our analysis.
Proposition 3.1.Let d ≥ 2, 0 < ε ≤ τ < 1, T > 0, and the threshold J τ be given by (2.1).Assume There exists a constant c > 0 independent of T , ε and τ such that if then for t ∈ (0, T ), the solution (y, w, r, u) of the Cauchy problem (3.1) satisfies X (t) := (y, w, r, u) ) where C 0 > 1 is a universal constant, and V(t) is denoted by (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 And direct produce law (2.4) for the equation (3.1) 1 gives that Similarly, we also get from (2.4) and (3.1) 3 that ) u and + u . (3.8) The conclusion of the proof will follow from Lemmas 3.1-3.3given and proven in the next three subsections.
Then making use of the Grönwall inequality and the smallness assumption (3.2) of Z(t), we obtain the uniform a-priori estimate (3.3).

Low-frequency analysis
Motivated by Darcy's law (1.3) 2 , we introduce the following effective flux which undergoes a purely damped effect in the low-frequency region |ξ| ≤ 2 k τ 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 where the higher-order linear terms L i (i = 1, 2, 3) are denoted as and the nonlinear terms R i (i = 1, 2, 3) are defined by (3.11) Now, to establish the estimates in low frequencies to the solutions of System (3.1) uniformly with respect to both ε and τ , 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 τ be given by (2.1).Then for t ∈ (0, T ), the solution (w, r, u) where Z(t), X (t), V(t) and z are defined by (3.2), (3.3), (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 τ .
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.

The Ḃ d 2 -estimates
We first perform Ḃ d 2 -estimates in low frequencies for the heat equation (3.10) 2 .It follows from the regularity estimate in Lemma 2.6 that . (3.13) Applying Lemma 2.7 to the damped equation (3.10) 1 , we get where we used inequality (3.13) to control terms involving r in equation (3.10) 1 .
Combining the above three estimates, we are led to

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 -estimates, and furthermore obtain the control of higher-order )-norms.
Lemma 3.2.Let T > 0, and the threshold J τ be given by (2.1).Then for any t ∈ (0, T ), the solution (w 0 , r 0 , u 0 ) h Proof.To prove of Lemma 3.2, we localize in frequencies for the equations (3.1) 2 -(3.1) with the commutator terms Multiplying (3.32) 3 by u j and integrating the resulting equation by parts, we get (3.34) Thence, we multiply (3.32) 1 by h5+H5 h1+H1 w j and integrate the resulting equation by parts to show Similarly, direct computations on (3.32) 2 yield (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 ∇r j that and taking the L 2 -inner product of (3.32) 2 with div u j that (3.38) In the spirit of the work [3] by Beauchard and Zuazua, for a small constant η * > 0 to be determined, we define the following Lyapunov functional with nonlinear weights as and its dissipation rate One derives from assumption (3.2) and the embedding which together with estimates (3.34)-(3.38)and the fact that 2 −j τ ≤ 1 for any j ≥ J τ − 1 yields the following Lyapunov inequality: (3.40) It follows from the smallness condition (3.39), Bernstein's inequality in Lemma 2.1 and the fact 2 −j τ that , and where one has used the condition ε ≤ τ .Thus, we can choose a sufficiently small constant η * > 0 independent of ε and τ so that Dividing the two sides of (3.40) by L j (t) + η for any η > 0, we have From (3.41) and the embedding Ḃ d 2 ֒→ L ∞ we get after integrating the above inequality over [0, t] and taking the limit as η → 0 that ) . (3.42) According to the commutator estimate (2.5), it follows that + τ r )

Z(t)X (t).
This with inequality (3.42) leads to On the other hand, for any η > 0, we deduce from inequality (3.35) that , which together with (3.43) implies Thanks to inequality (2.2), one has ) . (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.

Recovering the Ḃ d 2 +1 -estimates
As explained in Remark 3.3, we need to establish the uniform )-norm estimate of (w, r, u) which in fact leads to the uniform control of ) at both low and high frequencies as a byproduct.Lemma 3.3.Let T > 0, and the threshold J τ be given by (2.1).Then for any t ∈ (0, T ), the solution (w, r, u) to the linear problem (3.1) 2 -(3.1) 4 satisfies (w, r, u) ) where η ∈ (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-τ of u for L 2 -time type norms.In fact, for any j ∈ Z, by combining inequalities (3.34)-(3.35)together, we get (3.46) Furthermore, from (3.46) we have (w, r, u) ) + (S 1 , S 2 , S 3 ) ) (w, r, u) (3.47) The right-hand side of inequality (3.47) can be estimated as follows.By the commutator estimate (2.5), we have Similarly, it holds j∈Z and ) We conclude from the above estimates that (w, r, u) ) ) Applying Hölder's inequality to the above estimate leads to inequality (3.45).
4 Global well-posedness for the nonlinear problems

The Cauchy problem of System (BN)
In this section, we prove the uniform in ε and τ global existence and uniqueness of solutions to the Cauchy problem for (BN) subject to the initial data (α ±,0 , ρ ±,0 , u 0 ).i.e.Theorem 1.1.For simplicity, we omit the superscript concerning the parameters ε and τ in this section.

• 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, ε, τ and time.To this end, we first let X 0 ≤ 1.It follows from (4.3) and the composition estimates in Lemma 2.5 that and similarly, 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 Hence, we are able to employ the uniform a-priori estimate established in Proposition 3.1 to obtain ) . (4.6) Applying (4.3), Lemma 2.3 and 2.5 gives directly where C * 4 > 0 is a universal constant.Combining (4.1), (4.3), (4.5), (4.6) and (4.7) together, we have as long as Thus, the uniform estimate (4.3) holds true for any n ≥ 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 → ∞, 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 System (4.2).In what follows, we show that the strong convergence holds in R + × R d for the whole sequence.Lemma 4.1.There exists a small constant c * 3 ∈ (0, min{1, c * 1 , c * 2 }] and a limit (y, w, r, u) such that if In particular, we have Proof.In order to show (4.8), one needs to perform uniform energy estimates on the error unknown Following the framework in Section 3, we aim to estimate the functional .
• Step 4: Global existence By virtue of the strong convergence properties (4.8)-(4.9),one can pass to the limit as n → ∞ 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.23).In addition, taking advantage of Fatou's property in Lemma 2.2, for all t > 0, we have To prove the time continuity property in (1.8)-(1.9),our proof relies on the uniform bound (4.17) and employs a reasoning analogous to that found in [22].Since ), we shall investigate each equations separately.Recall that the solution (y, w, r, u) satisfies As the right-hand side terms of the components y, r and w belong to . We are left with recovering the time continuity of (y, w, r, u) in . First, for a fixed j ∈ 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) ), for any η > 0, there exists an large integer J * such that, for all t > 0 |j|≥J * Thus, for any time t, t ′ ∈ R + , we have ).Finally, applying the inverse function theorem, we can see that once α ± and ρ ± are determined by (1.21), then ) is the unique global strong solution to the original Cauchy problem of System (BN).Using (4.17), product laws and composite estimates, we are able to verify that (α ± , ρ ± , u) satisfies the properties (1.8)-(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.23).Without loss of generality, as the parameters do not affect our argument for proving uniqueness, we set ε = τ = 1.Let (y 1 , w 1 , r 1 , u 1 ) and (y 2 , w 2 , r 2 , u 2 ) be two solutions to (1.23) subject to the same initial data (y 0 , w 0 , r 0 , u 0 ), satisfying (1.8) and Fi + G i (y j , r j , w j ) > 0 for i = 0, 1, 2, 3 and j = 1, 2. The difference w, r, u)(0, x) = (0, 0, 0, 0), (4.18)where we denoted Here G l i := G i (y l , w l , r l ) and F l 4 := 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, Through the application of the weighted Lyapunov functional method, as outlined in (3.32)-(3.36),we where W i = W i (y 1 , w 1 , r 1 ) > 0, i = 1, 2, 3 are some smooth weight functions depending on Fi + G 1 i > 0, and the commutator terms are given by Integrating (4.20) over [0, t], taking the square root of both sides and summing the resulting estimate over j ∈ Z with the weight 2 d 2 j , we get ( w, r, u) + u + w ( w, r, u) , which together with Young's inequality implies Using the classical commutator estimate in Lemma 2.4 implies In addition, according to standard product laws and composite estimates in Lemmas 2.3 and 2.5, the nonlinear terms ( S 1 , S 2 , S 3 , S 4 ) can be handled as Substituting the above two estimates into (4.19) and (4.21) and then employing Grönwall's inequality, we end up with ( y, w, r, u)(t) Ḃ d 2 = 0 for all t > 0, which concludes the proof of Theorem 1.1.

The Cauchy problem of Systems (K) and (PM)
We provide a brief explanation of the proof of the global existence and uniqueness for System (PM).
The proof of the result for 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 System (K) allows us to construct solutions for System (PM) by taking the limit as the relaxation parameter τ → 0.
The following lemma states the uniform estimate verified by the solutions of System (K τ ), which are rescaled from estimate (1.12) obtained in Theorem 1.2 for System (K).
This together with the uniform estimate (1.9) leads to (5.4).
We are going to estimate (δQ, δu).By virtue of (BN), (K) and (5.6)-(5.7),(δQ, δu) satisfies the following equations of damped Euler type with rough coefficients: with the nonlinear terms In order to establish the uniform-in-τ convergence estimates, we follow the ideas in Section 3 to overcome the issue caused by the overdamping phenomenon.
For the third difficult term in δF 4 , we apply (2.2), the product law (2.4) for d ≥ 3 and the fact .
For the tricky commutator term R 3,j , we have ) .
It is also one of the reasons why we need to perform the Ḃ 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 L ∞ .19) Remark 3.2.The above estimate (3.19) for H 4 ∇r arising from two pressures implies that one needs uniform Ḃ d 2 −1 -estimates for low frequencies.Indeed, as H 4 does not have the either L 1 -in-time or L 2 -in-