Nonlinear Rayleigh-Taylor Instability for Nonhomogeneous Incompressible Viscous Magnetohydrodynamic Flows

We investigate the nonlinear instability of a smooth Rayleigh-Taylor steady-state solution (including the case of heavier density with increasing height) to the three-dimensional incompressible nonhomogeneous magnetohydrodynamic (MHD) equations of zero resistivity in the presence of a uniform gravitational field. We first analyze the linearized equations around the steady-state solution. Then we construct solutions of the linearized problem that grow in time in the Sobolev space $H^k$, thus leading to the linear instability. With the help of the constructed unstable solutions of the linearized problem and a local well-posedness result of smooth solutions to the original nonlinear problem, we establish the instability of the density, the horizontal and vertical velocities in the nonlinear problem. Moreover, when the steady magnetic field is vertical and small, we prove the instability of the magnetic field. This verifies the physical phenomenon: instability of the velocity leads to the instability of the magnetic field through the induction equation.


Introduction
This paper is concerned with nonlinear instability of a smooth Rayleigh-Taylor (RT) steadystate solution to the following three-dimensional (3D) nonhomogeneous incompressible magnetohydrodynamic (MHD) equations with zero resistivity (i.e. without magnetic diffusivity) in the presence of a uniform gravitational field (see, for example, [2,25,26,28] on the derivation of the equations): (1.1) Here the unknowns ρ := ρ(t, x), v := v(t, x), M := M(t, x) and p := p(t, x) denote the density, velocity, magnetic field and pressure of the incompressible fluid, respectively; µ > 0 stands for the coefficient of shear viscosity, g > 0 for the gravitational constant, e 3 = (0, 0, 1) for the vertical unit vector, and −ρge 3 for the gravitational force.
In the system (1.1) the equation (1.1) 1 is the continuity equation, while (1.1) 2 describes the balance law of momentum. It is well-known that the electromagnetic field is governed by the Maxwell equations. In MHD, the displacement current can be neglected [25,26]. As a consequence, (1.1) 3 is called the induction equation. As for the constraint div M = 0, it can be seen just as a restriction on the initial value of M since (div M) t = 0. We remark that, the resistivity in (1.1) 3 is zero, which arises in the physics regime with negligible electrical resistance, see [4]. In addition, if M ≡ 0, the system (1.1) reduces to the incompressible Navier-Stokes equations in the presence of a uniform gravitational field.
We point out that by virtue of the condition (1. 3), there is at least a region in which the RT density profile has larger density with increasing x 3 (height), thus leading to the classical RT instability as shown in Theorem 1.1 below. Now, we denote the perturbation to the RT steady state by then, (̺, u, q) satisfies the perturbed equations          ̺ t + u · ∇(̺ +ρ) = 0, (̺ +ρ)u t + (̺ +ρ)u · ∇u + ∇q + g̺e 3 = µ∆u + (∇ × N) × (N +M), N t = ∇ × (u × (N +M)), divu = 0, divN = 0. (1.7) The RT instability is well-known as gravity-driven instability in fluid dynamics when heavy fluid is on top of light one. The linear instability for an incompressible fluid was first introduced by Rayleigh in 1883 [31]. The analogue of the RT instability arises when fluids are electrically conducting and a magnetic field is present, and the growth of the instability will be influenced by the magnetic field due to the generated electromagnetic induction and the Lorentz force. Some authors have extended the partial results concerning the RT instability to the case of MHD fluids by overcoming difficulties induced by presence of the magnetic field. For example, Kruskal and Schwarzchild in 1954 first showed that a horizontal magnetic field has no effect on the development of the linear RT instability [24]. Then the influence of a vertical magnetic field was investigated by Hide in [15] where the effect of finite viscosity and resistivity was included and his analysis was encumbered with many parameters. By a variational approach, Hwang in 2008 studied the nonlinear RT instability of (1.4)-(1.6) for the inviscid case (i.e. µ = 0) in a 2D periodic domain [16].
To our best knowledge, however, it is still open mathematically whether there exists an unstable solution to the nonlinear RT problem (1.4)-(1.6) of 3D viscous MHD fluids. The aim of this article is to rigorously verify the instability for the nonlinear RT problem (1.4)-(1.6). Moreover, the impact of the magnetic filed on the instability will be analyzed, for example, we shall show that if the steady magnetic field is vertical and small, then the magnetic field is unstable, thus verifying the physical phenomenon: instability of the velocity leads to the instability of the magnetic field through the induction equation. The main result of this paper reads as follows. , similarly to that in the derivation of (2.6), we can also deduce a ODE corresponding toM = (M 1 , M 2 , M 3 ). The resulting ODE, however, possesses the terms iψ ′ and iψ ′′′ which destroy the good variational structure, and consequently, one could not directly construct the growing mode solutions to the linearized problem. In addition, Theorem 1.1 also holds when M = N = 0, hence the weak RT instability for incompressible viscous fluids given in [20] can be further shown to be strongly RT unstable as in (1.10) in the case of the horizontally periodic domain. Remark 1.3. In [17] Hwang and Guo proved the nonlinear RT instability for 2D nonhomogeneous incompressible inviscid flows (i.e. µ = 0 andM = N = 0 in (1.4), (1.5)) with boundary condition u·n| ∂Ω ′ = 0 where Ω ′ = {(x 1 , x 2 ) ∈ R 2 | −l < x 2 < m} and n denotes the outer normal vector to ∂Ω ′ . Later, Hwang [16] further investigated the nonhomogeneous incompressible inviscid MHD fluid on a periodic domain, and get the instability of the norm (̺, u, N) L 2 (Ω ′ ) . Our result is more precise than those in [16,17] in the sense that Theorem 1.1 reveals that the vertical velocity induces the instability of the density and horizontal velocity; and moreover, we can show the instability of the magnetic field for the caseM = Me 3 . This mathematically verifies the physical phenomenon: instability of the velocity further leads to the instability of the magnetic field through the induction equation.
) and is finite for Rρ ′ dx < 0 (i.e.,ρ(+∞) <ρ(−∞)), see Proposition 2.1 for the detailed proof. This means that any vertical steady magnetic field cannot restrain growth of the nonlinear RT instability for the case Rρ ′ dx > 0. Now it is still not clear to us that whether a sufficiently large vertical steady magnetic field has impact on growth of the nonlinear RT instability for the case Rρ ′ dx < 0 due to some technical difficulties. However, if we consider Ω = (2πLT) 2 × (−l, m), then, similarly to the derivation of (3.71) in [32], we can show the stability of the velocity for any classical solution of the linearized problem satisfying boundary conditions u| x 3 =−l = u| x 3 =m = 0, provided the vertical steady magnetic field is sufficiently large. This result does not contradict that in [16] where for any vertical steady magnetic field, the nonlinear RT instability for a nonhomogeneous incompressible inviscid MHD flow in a periodic domain (the vertical direction is also periodic) is shown. These mathematical results reveal that the domain and boundary conditions of the velocity do have impact on the instability of MHD flows.
The proof of Theorem 1.1 is divided into four steps given in Sections 2-4: (i) First, we make an ansatz to seek the "normal mode" solutions of the linearized equations (1.7), which grow exponentially in time by the factor e λ(ξ)t with ξ ∈ R 2 being the horizontal spatial frequency and λ(ξ) > 0. This reduces the equations to a system of ODEs defined on R with λ(ξ) > 0 for some ξ (see (2.4), (2.5)). All such points ξ constitute a solvable domain. Because of presence of the magnetic field, the solvable domain is not a ball for the horizontal case as in [12], resulting in some difficulties in constructing a solvable domain. In order to circumvent such difficulties, similarly to [21,32], we introduce the notions of the critical frequency function S(ξ) and critical frequency constant |ξ| M vc to define the solvable domain A g (cf. (2.14), (2.15)). Thus, by careful constructing the solutions and adapting the modified variational method in [12], we can also solve the ODEs for any given ξ ∈ A g and obtain a normal mode with λ(ξ) > 0, thus leading to a mechanism for the global linear RT instability. Using the normal modes, we can further construct real-valued solutions of the linearized problem (1.5)-(1.7) that grow in time, when measured in H k (Ω) for any k ≥ 0. In particular, the density, horizontal velocity and vertical velocity in the solutions of the linearized problem are not zero, this fact will play a key role in the nonlinear instability in (1.10). (ii) In Section 3, we state a local well-posedness result of the perturbed problem (1.4)-(1.6), which will be proved in Section 5. Then we derive the integrand form of Gronwall's inequality of high-order energy estimate E(̺, u, N) for the perturbed problem, and this makes the escape time occurring before break-down of the classical solutions. Since the equilibrium state of the magnetic fieldM is a no-zero vector, we shall introduce a simple technique to deal with the terms includingM in the energy estimates, see Subsection 3.4. (iii) Finally, in Section 4, with the help of the results established in Sections 2-3, we adapt a careful bootstrap argument as in [11] to establish the instability of the nonlinear problem. We mention that although the approach in [11] has been also used to treat the instability of other problems (see [10,11,19,33] for examples), but Duhamel's principle in the standard bootstrap argument can not be directly applied to our problem, since the nonlinear terms in (1.4) do not satisfy the compatibility condition of divergence-free. To circumvent this obstacle, we shall use some specific energy estimates to replace Duhamel's principle. Moreover, we can also find in the proof that Λ is indeed a sharp exponential growth rate for general solutions to the linearized problem (see Remark 4.1).
We end this section by briefly reviewing some of the previous results on the nonlinear RT instability for two layer incompressible fluids separated by a free interface (stratified fluids), where the RT steady-state solution is a denser fluid lying above a lighter one separated by a free interface. When the densities of two layer fluids are two constants, Wang and Tice [33] proved the (local) existence of nonlinear unstable solutions in a horizontally periodic domain T 2 × (−b, 1), where the instability term is described by the sum of L 2 -norm of the velocity and the moving internal interface. Prüss and Simonett used the C 0 -semigroup theory and the Henry instability theorem to show the existence of nonlinear unstable solutions in the Sobolev-Slobodeckii spaces in R 3 [30], where the instability term is described by the sum of (see [30,Theorem 1.2] for details). When densities of two layer fluids are variable, to our best knowledge, the (local) existence of solutions to the nonlinear problem (1.4)-(1.6) still remains open, consequently, the strong nonlinear instability is still open. For compressible fluids, there are very few results on the nonlinear RT instability. Guo and Tice proved the instability of immiscible compressible inviscid fluids in the frame of Lagrangian coordinates under the existence assumption of solutions [13]. This is in some sense a compressible analogue to the local illposedness of the RT problem for incompressible fluids given in [6].
Finally we mention related results on the instability for stratified MHD fluids. Wang [32] introduced the critical number of stratified MHD fluids (denoted by M s c ) to investigate the linear RT instability of stratified MHD fluids in a infinite slab domain R 2 × (−l, l). Later, Jiang et al. [21] further showed the weak nonlinear RT instability of stratified MHD fluids and found that M s c = gl(̺ + − ̺ − )/2, where ̺ + > ̺ − , and ̺ + resp. ̺ − denotes the density of the upper-resp. lower-layer fluid. Obviously, M s c → ∞ as the height 2l → ∞. Hence, the critical number of stratified MHD fluids in R 3 is infinite. This means that the vertical magnetic field may have no impact on growth of the RT instability in R 3 .

Construction of solutions to the linearized problem
We wish to construct a solution to the linearized equations (1.7) that has growing H k -norm for any k. We will construct such solutions via some synthesis as in [12] by first constructing a growing mode for any but fixed spatial frequency. Moreover, we shall introduce the techniques of the critical frequency function and critical frequency constant to carefully analyze the terms involving with the magnetic field.

Linear growing modes
To start with, we make a growing mode ansatz of solutions, i.e., for some λ > 0, Substituting this ansatz into (1.7), and then eliminatingρ andM by using the first and third equations, we arrive at the time-independent system forṽ = (ṽ 1 ,ṽ 2 ,ṽ 3 ) andp: whereM is given by (1.9). After a straightforward calculation, we find that We fix a spatial frequency ξ = (ξ 1 , ξ 2 ) ∈ R 2 , and define the new unknowns Then, in view of (2.1)-(2.3), we see that ϕ, θ, ψ and λ satisfy the following system of ODEs: Eliminating π from the third equation in (2.4), we obtain the following ODE for ψ Next we use the modified variational method to construct a solution of (2.6), (2.7). This idea can be found in Guo and Tice's paper on compressible viscous stratified flows [12], and has been adapted by other authors to study the instability for other fluid models [5,19,22,32]. To this end, we now fix a non-zero vector ξ ∈ R 2 and s > 0. From (2.6) and (2.7) we get a family of the modified problems coupled with the condition (2.7). We define the energy functional of (2.8) by with an associated admissible set where (2.11) Thus we can find a −λ 2 (depending on ξ) by minimizing In order to emphasize the dependence on s ∈ (0, ∞) we will sometimes write E(ψ, s) := E(ψ) and α(s) := inf ψ∈A E(ψ, s) < +∞.
Before constructing the growth solutions, we shall introduce some preliminary results, which will be used in Subsections 2.2-2.4. Let the critical frequency function S(ξ) for the horizontal case and the critical frequency constant |ξ| M vc for the vertical case be given by the following variational forms S(ξ) := sup > 0 (2.14) for 0 < |Mξ 1 |/|ξ| < M c , and respectively, where M c is given by (1.8), and We remark that it depends on the choice ofρ whether M c becomes infinite or finite. More precisely, we have the following conclusions: Proof.
Then ψ n (x) ∈ H 1 (R) and which implies M c = ∞. We mention that, in such a case, we can infer that the critical frequency constant |ξ| M vc defined by (2.15) is equal to zero. (2) We turn to show the second assertion by contradiction. Assume that M c is infinite, i.e., there exists a sequence of functions {ψ n } ∞ n=1 such that ψ n ∈ H 1 (R), and 0 < R gρ ′ |ψ n | 2 dx R |ψ ′ n | 2 dx → ∞ as n → ∞. Denote Now, using Newton-Leibniz's formula, we immediately get that, for any given ε ∈ (0, 1), there exists a N ε > 0 (may depend on ε), such that, for any n > N ε and for any x ∈ (−l, l), Thus, recalling the conditions Rρ ′ dx < 0 andρ ∈ C 0 0 (R), and using the following relation we infer that there exists a constant c, independent of n and ε, such that |ϕ n (0)| ≤ c for any n > N. (2.18) In fact, if this is not true, then there exists a subsequence, denoted by {ϕ nm (0)} ∞ m=1 , satisfying |ϕ nm (0)| ≥ m, which implies that 1 = R(n m ) → ∞. This is a contradiction.
Finally, making use of (2.16)-(2.18) and Rρ ′ dx < 0, we have that for any n > N ε , Consequently, letting ε → 0, we immediately see that R(n) can be sufficiently small for some n, and this contradicts with R(n) ≡ 1. Hence the second assertion holds.
Proposition 2.3. Let A g be defined by (2.19), then (1) the set A g is symmetric on x-axis and y-axis in R 2 , respectively; (2) there exist countably infinite lattice points of (L −1 Z) 2 belongs to A g ; Proof. The first two assertions obviously hold by virtue of the definition of A g . It suffices to show the last assertion. Here we only give the proof of the horizontal case for the reader's convenience.
. For the first case, noting that |ξ 1 M|/|ξ| ∈ (0, M c ) and |ξ| < S(ξ) hold for any |ξ| > 0 with sufficiently small ξ 1 there exists a sufficiently small disk B δ For the latter case, noting that 0 < |ξ 01 M|/|ξ 0 < M c and |ξ 0 | < S(ξ 0 ), we use the continuity of |ξ 1 M|/|ξ|, |ξ| and S(ξ) as ξ → ξ 0 = 0 to deduce that there also exists a sufficiently small disk . Summing up the previous discussions, we immediately conclude that the set A g is a nonempty open set in R 2 by the definition of open sets.
Proof. The above assertions in fact follow from the definitions (1.8), (2.11), (2.14), (2.15) and (2.19), here we only give the proof of the horizontal case for the reader's convenience. Let ξ ∈ A g . If ξ 1 = 0, then obviously, there exists a ψ 0 , such that If ξ satisfies |Mξ 1 |/|ξ| ∈ (0, M c ) and |ξ| < S(ξ), then by virtue of the definition of (2.14), there also exists a ψ 0 , such that Summing up the above discussions, we see that there exists a ψ 0 , such that E 0 (ψ 0 ) < 0 for ξ ∈ A g . Let ξ ∈ A g ∪ {0}, then ξ can be divided by two cases: For the first case, in view of the definition of (2.14), we find that Finally, for the second case, by the definition of (1.8), we have which yields Summarizing the above discussions, we see that This completes the proof for the horizontal case.

Solutions to the variational problem
In this seusection we show that a minimizer of (2.13) exists for the case of inf ψ∈A E(ψ, s) < 0 which will be shown to be true for sufficiently small s in Proposition 2.6 below, and that the corresponding Euler-Lagrange equations are equivalent to (2.7), (2.8).
Proposition 2.5. For any fixed s > 0 and ξ with |ξ| = 0, we have (3) letψ be a minimizer and −λ 2 := E(ψ), then the pair (ψ, Proof. We only show the proposition for the horizontal case, the vertical case can be dealt with in the same manner. (1) Noticing that for any ψ ∈ A, we see that E is bounded from below on A by virtue of (1.2). This proves (1).
(2) We proceed to show (2). Let ψ n ∈ A be a minimizing sequence, then E(ψ n ) is bounded. This together with (2.9) and (2.24) implies that ψ n is bounded in H 2 (R). So, there exists a ψ 0 ∈ H 2 (R), such that ψ n → ψ 0 weakly in H 2 (R) and strongly in H 1 loc (R). Moreover, by the lower semi-continuity, (1.2) and the assumption that E(ψ) < 0 for someψ ∈ A, we deduce that Suppose by contradiction that J(ψ 0 ) < 1. By the homogeneity of J we may find an α > 1 so that J(αψ 0 ) = 1, i.e., we may scale up ψ 0 so that αψ 0 ∈ A. From this we deduce that which is a contradiction since αψ 0 ∈ A. Hence J(ψ 0 ) = 1 and ψ 0 ∈ A. This shows that E(ψ) achieves its infinimum on A.
Next, we want to show that there is a fixed point such that λ = s. To this end, we first give some properties of α(s) as a function of s > 0.
Proposition 2.6. The function α(s) defined on (0, ∞) enjoys the following properties: (2) for any ξ ∈ A g , there exist constants c 1 , c 2 > 0 depending on g, M,ρ, µ, and ξ, such that (2.28) Proof. We still give the proof for the horizontal case only, and the vertical case can be dealt with in the same way.
(1) Let {ψ n s 2 } ⊂ A be a minimizing sequence of inf ψ∈A E(ψ, s 2 ). From (2.9) and (2.12) it follows that Hence α(s) is nondecreasing on (0, ∞). Next we shall use this fact to show the continuity of α(s).
In view of (2.24) and the monotonicity of α(s), On the other hand, for any s ∈ I, there exists a minimizing sequence {ψ n s } ⊂ A of inf ψ∈A E(ψ, s), such that |α(s) − E(ψ n s , s)| < 1. Making use of (2.9)-(2.12) and (2.29), we infer that For s i ∈ I (i = 1, 2), we further find that Reversing the role of the indices 1 and 2 in the derivation of the inequality (2.30), we obtain the same boundedness with the indices switched. Therefore, we have On the other hand, Putting (2.31) and (2.32) together, we see that there is a subsequenceψ n 0 ∈ H 2 (R), such thatψ n 0 ≡ 0 and Thus, we have for two positive constants c 1 := c 1 (g, M,ρ, |ξ|) and c 2 := c 2 (g, M, µ,ρ, |ξ|). This completes the proof for the vertical case.
Moreover, in view of Propositions 2.7 and 2.8, we conclude the following existence for the problem (2.6)-(2.7).
We end this subsection by giving additional properties of the solutions established in Theorem 2.1 in terms of λ(ξ), which show that λ is a bounded, continuous function of ξ.
(i) We begin with the proof of the following conclusion: By virtue of Proposition 2.7, for any ξ ∈ A g , there exists a function ψ ξ ∈ A, such that which, together with (2.10), yields where c 3 depends on g, M, µ,ρ, a, b and s.

Construction of a solution to the system (2.4), (2.5)
A solution to (2.6), (2.7) gives rise to a solution of the ODEs (2.4), (2.5) for the growing mode velocity u as well. (2.5), and the solution belongs to H ∞ (R). Moreover ψ ∈ A.
Remark 2.2. We mention that the system (2.4), (2.5) for the vertical case enjoys rotational structure. Hence, we can also use the rotation method as in [20] to construct a solution of (2.4), (2.5), which is simper than the above construction process.
Lemma 2.1. Let k ≥ 0, and D ⊂ A g be a nonempty bounded set satisfying the closureD ⊂ A g . Then for any ξ ∈ D there are positive constants A k , B k , C k and D k , which may depend on g, M, µ,ρ and D, such that

49)
where (ϕ, θ, ψ, π)(ξ) and λ(ξ) be constructed as in Theorem 2.2. Moreover, Proof. We still give the proof for the horizonal case only. Throughout this proof, we denote byc a generic positive constant which may vary from line to line, and may depend on g, M, µ, ρ and D.
We begin with the estimate (2.46). We first rewrite (2.6) as which yields (2.53) Noting that since ψ ∈ A, the inequality (2.50) holds, and there is a constantc, such that

Exponential growth rate
In this subsection we will construct a linear real-valued solution to the linearized problem (1.7) along with (1.6) which grows in-time in the Sobolev space of order k. (2) For every t > 0, (̺, u, N, q)(t) ∈ H k , and

Nonlinear energy estimates for the perturbed problem
In this section, we derive some nonlinear energy estimates for the perturbed problem (1.4)-(1.6), and the integrand form of Gronwall's inequality of the high-order energy estimates.
First we shall give a local well-posedness result of the perturbed problem. We mention that the local existence of classical solutions and global existence of classical small solutions to the non-resistive MHD equations have been established by many authors, see [7,18,23,27,28] for example. To our best knowledge, there is no existence result for the MHD equations (1.1) in the horizontally periodic domain Ω. However, with the help of the usual approaches in the proof of local existence for fluid dynamical equations and some mathematical techniques to deal with the horizontally periodic domain, we can establish the following local well-posedness result, the proof of which will be offered in Section 5 for the completeness. With Proposition 3.1 in hand, we further derive the integrand form of Gronwall's inequality of the high-order energy estimates. It should be pointed out that the solution (̺, u, N, q) constructed in Proposition 3.1 possesses more regularity, see the proof in Section 5. This additional regularity makes it sense to derive some identities and high-order energy estimates later. In particular, we omit the standard regularization argument in the derivation of some identities, (3.16) for example, we refer to (5.23) or (5.7) in Section 5 for the proof.

Estimate of the perturbation density
We first note that by the classical Sobolev embedding results (see [1,Chapter 5]), we have With the help of the above estimates, we can bound the perturbation density ̺. In fact, applying ∂ α to (1.4) 1 , multiplying the resulting identity by ∂ α ̺ in L 2 (Ω), we get where we have defined ∂ β u · ∂ γ ∇̺ for |α| = 3, so that (3.4) makes sense (cf. Lemma 5.1). Using (3.3), Hölder's and Young's inequalities, I 1 can be bounded as follows: On the other hand, we use (3.1), (3.2) and Hölder's inequality to control I 2 as follows.

Estimate of the perturbation magnetic field
We continue to bound the perturbation magnetic field N. Applying ∂ α to (1.4) 3 , multiplying the resulting identity by ∂ α N in L 2 , we obtain where we have used the facts ∇ × (u ×M) =M · ∇u and ∇ × (u × N) = N · ∇u − u · ∇N, and defined Similarly to the derivation of (3.5) and (3.6), we can control J 1 and J 2 as follows.

Estimate of the perturbation velocity
First we restrict the density ρ = ̺ +ρ so that it has a positive lower bound. By virtue of (3.2), there is a constant δ ′ 0 ∈ (0, 1), such that Consequently, in view of (1.4) 1 and divu = 0, we find that for any (t, x) ∈ [0, T ) × Ω by the method of characteristics, where ρ 0 = ̺ 0 +ρ. From now on, we always assume that E ≤ δ 0 ≤ δ ′ 0 < 1. We remark that the upper and lower boundedness of ρ will be repeatedly used below.
We now estimate the time-derivative of the perturbation velocity. Multiplying (1.4) 2 by u in L 2 , utilizing (1.4) 1 and integrating by parts, we obtain where Y 0 := (∇ × N) × (N +M) − g̺e 3 . Applying ∂ t to (1.4) 2 , multiplying the resulting identity by u t in L 2 , we make use of (1.4) 1 , and integrate by parts to deduce Applying ∂ 2 t to (1.4) 2 again, multiplying the resulting identity by u t , employing (1.4) 1 and Adding the above three equalities, we get (3.17) Obviously, Recalling that E ≤ δ < 1, we use (3.1), (3.2), (3.8), (3.9), (3.12) and (3.13) to arrive at To bound Y 2 , we argue in a way similar to that used for (3.19), with additional estimates (3.10) and (3.3), to infer that (3.20) Finally, using (3.12)-(3.14), we obtain (3.21) Putting (3.17)-(3.21) together and using Cauchy-Schwarz's inequality, we conclude and multiplying (3.23) by u t in L 2 , we have Differentiating (3.23) with respect to t and multiplying the resulting equations by u tt in L 2 , we obtain µ 2 Adding the above two equalities, we get Consequently, in view of (3.19), (3.20) and (3.25), using Cauchy-Schwarz's inequality, we get from (3.24) that d dt (3.26) Next, we continue to derive more estimates of higher-order derivatives of the perturbation velocity. Multiplying (1.4) 2 by u t , integrating and using (3.3), one gets while applying ∂ i to (1.4) 2 , multiplying the resulting equations by ∂ i u t in L 2 , and using (3.3), we have (3.28) If we take ∂ i ∂ j to (1.4) 2 and multiply the resulting equations with ∂ i u t in L 2 , we see that In particular, summing up the above three estimates, we conclude provided δ 0 is sufficiently small.
Finally, summing up the above estimates and the existence statement, we arrive at the following conclusion: be constructed in Proposition 3.1, and δ ′ 0 be chose as in (3.15).

39)
and where Λ * is provided by Theorem 2.3.
We then define where T max denotes the maximal time of existence. Obviously, T * T * * > 0, and furthermore, Then for all t ≤ min{T δ , T * , T * * }, we deduce from the estimate (3.39), and the definitions of T * and T * * that for some constant C 2 := C 2 (δ 0 ) > 0 independent of δ. Then applying Gronwall's inequality to the above estimate, one obtains for any t ≤ min{T δ , T * , T * * }.
is also a linear solution to (1.7) with the initial data (̺ δ 0 , u δ 0 , N δ 0 ) ∈ H ∞ . Moreover, (̺ d , u d , N d ) satisfies the following non-homogenous equations: From the estimates (3.40) and (4.6), it follows that Moreover, we have the following error estimate for ̺ d , u d , N d : There exists a constant C 3 such that Proof. The proof is divided into two steps, in which we omit the superscript "d" in ̺ d , u d , N d for simplicity.
Step 2: Proof of the error estimate (4.9) In what follows, we still denote by C a generic positive constant which may depend on physical parameters and the known functions̺ 0 ,ū 0 , andN 0 .
We differentiate (4.7) 2 withM = Me i (i = 1, 3) in time, multiply the resulting equation by u t and integrate by part over Ω. Then, using the first and third equations in (4.7), we obtain (4.20) Here, utilizing (4.1), (4.6) and (4.8), we see that Noting that u t (0) ≤ Cδ 3 , we integrate (4.20) from 0 to t and use (4.10) to infer that The following inequality follows easily from integration in time and Cauchy-Schwarz's inequality.  On the other hand, Hence, putting (4.21)-(4.23) together, we obtain the differential inequality: An application of Gronwall's inequality then implies that for any t ≥ 0, where we have used the fact that 3Λ * − 2Λ > 0. Thus, making use of (4.21), (4.22) and (4.24), we deduce that Now we use (4.24), (4.25) and (4.7) 1 to arrive at where X = L 2 or H 1 . Similarly to the derivation of (4.26), one obtains Thus, putting the above four estimates together, we immediately get (4.9). This completes the proof of the lemma.

Proof of the local well-posedness
In this section, we adapt the standard method to briefly show the local well-posedness for the problem (1.4)-(1.6) by three steps, in which we shall exploit some mathematical techniques to deal with the horizontally periodic domain. Firstly we solve the following linearized problem for given v: with initial and boundary conditions: where (u 0 , N 0 ) satisfies the compatibility conditions divu 0 = 0 and divN 0 = 0. Secondly, we use the technique of iteration to construe a sequence of solutions {(̺ k , u k , N k )} ∞ k=1 , in which (̺ k , u k , N k ) solves the above linearized problem with (̺ k , u k , N k ) in place of (̺, u, N) and u k−1 in place of v, and show the uniform boundedness of the solution sequence in some function space. Finally, we further prove that the sequence of solutions is a Cauchy sequence, and thus the limit function is the unique solution of the original problem.
Throughout the rest of this article we shall repeatedly use the abbreviations:

Unique solvability of linearized problems
To show the solvability of the linearized problem (5.1)-(5.2), it suffices to solve the following two problems   with initial and boundary conditions where we have defined that v · ∂ α ∇̺ = 0 and v · ∂ α ∇N = 0 for any |α| = 3.
Proof. Following the proof of [ The strong continuity (5.7) and (5.8) can be shown by adapting the idea in the proof of [29,Lemm 6.9]. Here we give the proof (5.7) for the reader's convenience. Without a generality, we assume that α = (3, 0, 0). Similarly in [29,Lemm 6.9], we can deduce from (5.3) that a.e. in Q T , where S ε is a standard mollifier, and r ε = div(S ε (∂ 3 where the constant c is independent of n. Multiplying (5.3) 1 by φ n ψ ∈ C ∞ 0 (I T ), and integrating the resulting equality, we get Thus, taking ε → 0, and then n → +∞, we immediately get (5.7).
We turn to the existence proof of a unique solution to the second problem (5.5)-(5.6) by employing the Galerkin method, the domain expansion technique and a technique of improving regularity. Lemma 5.1, then there exists a unique classical solution u ∈ V T to the problem (5.5)-(5.6) with an associated pressureq satisfying q ∈ L ∞ (I T , L 6 ), ∇q ∈ C 0 (Ī T , H 2 ), (∇q) t ∈ C 0 (Ī T , L 2 ). (5.11) Proof. We divide the proof of Lemma 5.2 into three steps.
(1) First we solve (5.5)-(5.6) without pressure by the Galerkin method, and the existence of solutions to the Galerkin approximate problem is established by adapting the basic idea from the proof of [14,Theorem 4.3]. Since By the theory of ODEs, we can find the coefficients a m j ∈ C 2 (Ī T ), such that with initial data u m (0) = P m u 0 ∈ U m for each w ∈ U m . Next, we derive uniform estimates for u m . In what follows, we denote by C(· · · ) a generic positive constant depending only on its variables. The notation a b means that a ≤Cb for a universal constantC > 0, which may depend on T ,M, µ, g, ρ,ρ,ρ, and the norms ( Moreover,C is increasing with T , and the norms ̺ t (0) L ∞ , f(0) H 2 and f t (0) L 2 , ifC depends on them.
Taking w = u m in (5.12), we see that Applying Gronwall's inequality, we get which, together with (5.13), implies that (5.14) We differentiate (5.12) in time and take w = u m t to deduce Finally, taking w = u m t in (5.12), we find that Summing up the previous estimates and making use of the estimate we arrive at (5.16) (2) Now we can obtain a strong solution from the above uniform estimates for the approximate solutions. In view of (5.16), up to the extraction of a subsequence, one has u m → u weakly-* in L ∞ (I T , H 1 σ ), u m t → u t weakly-* in L ∞ (I T , L 2 ), ∇u m → ∇u weakly in L 2 (I T , L 2 ), u m t → u t weakly in L 2 (I T , H 1 ), u m → u strongly in C 0 (Ī T , L 2 (K)) for any bounded domain K ⊂ Ω, u(0) = u 0 , and divu t = 0. Thus, if we take limit in (5.12), we get which can be regarded as a weak solution of the following Stokes equations Next we proceed to derive more estimates for u in the Stokes equations by the domain expansion technique (i.e., the classical regularity theory on the Stokes equations). Let L n = (−n, n) and Ω n = (2πLT) 2 × L n . Similarly to [14, Propositions 2.9 and 3.7], we can show that the Stokes problem    −µ∆u n + ∇q n = F in Ω n , divu n = 0, u| x 3 =−2πnL = u| x 3 =2πnL = 0 (5.18) admits a unique strong solution u n with a unique associated pressureq n , such that ∇ k u n 2 H 2 (Ωn) + ∇ k+1q n 2 L 2 (Ωn) ≤ C(Ω n , µ) ∇ k F 2 L 2 (Ωn) for k = 0, 1, q n 2 L 2 (Ωn) ≤ C(Ω n ) F 2 L 2 (Ωn) and Ωnq n dx = 0.
Noting that u n L 6 (Ωn) ≤ C(Ω n ) ∇u n L 2 (Ωn) , we may scale the spatial variables on a cuboid domain (0, 2πnL) 2 ×L n and utilize the horizontally periodic property to get u n L 6 (Ωn) ≤ c ∇u L 2 for some constant c.
This shows that C ∞ 0,σ is dense in H 1 σ .

Uniform estimates
Let (̺ 0 , u 0 , N 0 ) satisfy the assumption in Proposition 3.1, then, in view of Lemmas 5.1 and 5.2, for any given T > 0 and any given v ∈ V T satisfying v(0) = u 0 , the linearized problem (5.1)-(5.2) possesses a unique classical solution ((̺, N), u) ∈ H T × V T . To emphasize such relation between u and v, we can define u := S(v). In addition, we can deduce from (5.1) 2 that where the positive constant C 0 depends on the initial data (̺ 0 , u 0 , N 0 ) and other physical parameters in the perturbed equations, and is independent of v. Now, let us denote We show next that there exists two positive constants T ∈ (0, 1) and κ ≥ 1 depending on (̺ 0 , u 0 , N 0 ) and other physical parameters in the perturbed equations, such that κ 8 T = 1 and S(v) V T ≤ κ for any v ∈ D T κ . (5.33) Throughout this and next subsections, the notation a b means that a ≤ Cb for some constant C > 0 which may depend on (̺ 0 , u 0 , N 0 ) and other physical parameters in the perturbed equations.