Relaxation of the Boussinesq system and applications to the Rayleigh-Taylor instability

We consider the evolution of two incompressible fluids with homogeneous densities $\rho_-<\rho_+$ subject to gravity described by the inviscid Boussinesq equations and provide the explicit relaxation of the associated differential inclusion. The existence of a subsolution to the relaxation allows one to conclude the existence of turbulently mixing solutions to the original Boussinesq system. As a specific application we investigate subsolutions emanating from the classical Rayleigh-Taylor initial configuration where the two fluids are separated by a horizontal interface with the heavier fluid being on top of the lighter. It turns out that among all self-similar subsolutions the criterion of maximal initial energy dissipation selects a linear density profile and a quadratic growth of the mixing zone. The subsolution selected this way can be extended in an admissible way to exist for all times. We provide two possible extensions with different long-time limits. The first one corresponds to a total mixture of the two fluids, the second corresponds to a full separation with the lighter fluid on top of the heavier. There is no motion in either of the limit states.


Introduction
We investigate two incompressible fluids with homogeneous densities 0 < ρ − < ρ + under the influence of gravity modelled by the Euler equations in Boussinesq approximation div v = 0, ∂ t ρ + div(ρv) = 0. (1.1) The equations are considered on a bounded domain Ω ⊂ R n and a time interval [0, T ), T > 0. The function ρ : Ω × [0, T ) → R is the normalized fluid density, i.e. ρ ∈ {±1} a.e., v : Ω × [0, T ) → R n is the velocity field and p : Ω × [0, T ) → R the pressure of the fluid. Furthermore, e n ∈ R n denotes the nth coordinate vector, g > 0 the gravitational constant and is the Atwood number. The incompressibility condition is complemented by the no-penetration boundary condition where ν denotes the exterior unit normal of the boundary of Ω, which is assumed to be sufficiently smooth. We will mostly consider (1.1), (1.2) with the unstable interface as initial data, i.e., ρ(x, 0) = sign(x n ), v(x, 0) = 0, x ∈ Ω. (1.3) This initial data is a classical instance of the Rayleigh-Taylor instability which occurs whenever a lighter fluid or gas is accelerated into a heavier one -a situation which appears in various research areas and applications, see [1,2,39,40] for an overview. Its linear (in)stability analysis goes back to Rayleigh [32] and Taylor [36].
In the present paper we consider the inviscid and indiffusive Boussinesq system (1.1). More generally one can also add different sorts of diffusion terms in the momentum balance and/or the mass balance. In dimension 2 global well-posedness results of sufficiently regular solutions and for different types of diffusion terms have been established in [7,15,22,24], while finite time singularity formation for equation (1.1), i.e. without any diffusive terms, has recently been shown in [18].
Local well-posedness statements for (1.1) considered in different sufficiently regular classes can be found in [8,14,18]. Note however that the horizontal interface (1.3) does not belong to these classes. On the other hand if one adds the diffusion terms −ν∆v, ν > 0 and −µ∆ρ, µ > 0 to the equations, local wellposedness has also been established for L p initial data [4].
Contrary to local well-posedness the articles [3,10] address the question of nonuniqueness of solutions. More precisely, [3] shows the existence of wild solutions for (1.1) considered with g = 0, i.e. for the homogeneous Euler equations augmented by a transport equation for a passive tracer. In [10] the effect of the Coriolis force and a diffusion term in the continuity equation is added to (1.1) and nonuniqueness of weak solutions to given initial data is proven. Moreover, the nonuniqueness of admissible weak solutions is also shown in [10] by construction of a suitable initial velocity field v 0 to a given initial density ρ 0 ∈ L ∞ (Ω) ∩ C 2 (Ω). The non-uniqueness results [3,10] both rely on the method of convex integration introduced by De Lellis and Székelyhidi to the context of fluid dynamics [16,17].

Heuristic outline of results
In this article we also address the inviscid Boussinesq system by means of convex integration, but our main goal here, as in [21] for the Euler equations (1.4), is to investigate nonlinear instability aspects of the Rayleigh-Taylor configuration (1.3) by providing the existence of solutions to the problem (1.1), (1.2), (1.3) that reflect a turbulent mixing. The oscillatory behaviour of solutions obtained by convex integration has also been utilized as an instance of turbulent mixing in the context of the Kelvin-Helmholtz instability [27,34] and the Muskat problem for the incompressible porous media equation [5,6,11,20,23,26,29,33]. Compared to [3,10] this requires the explicit knowledge of the relaxation associated with (1.1), which will be given here.
Using this relaxation we construct solutions to (1.1), (1.2), (1.3) considered on an n-dimensional quader Ω = (0, 1) n−1 × (−L, L), n ≥ 2, which at any time t > 0 with 1 3 gAt 2 ≤ L are turbulently mixing in the space region U (t) := x ∈ Ω : |x n | < 1 3 gAt 2 . The solutions all have a underlying self-similar subsolution in common, whose ρ-component is linear inside the mixing zone, i.e. ρ sub (x, t) = 3xn gAt 2 for x ∈ U (t). The subsolution and hence the growth rate of the mixing zone 1 3 gAt 2 is selected uniquely and independently of the dimension by asking for maximal initial energy dissipation among all self-similar subsolutions. In particular the induced solutions are admissible with respect to the initial energy. The usage of maximal energy dissipation is motivated by the entropy rate admissibility criterion for hyperbolic conservation laws [13] and has been investigated in [27,34] in the context of Euler subsolutions emanating from vortex-sheet initial data, as well as in [9,19] for compressible Euler systems. As in [27] we focus here on maximal initial dissipation, i.e. the selection applies to small times.
Beyond small times, we show that the subsolutions can be extended in an admissible way past the time where the mixing zone hits the boundary. In fact we provide two possible extensions which in the long-time limit converge to two different states of the fluid: the fully mixed, isotropic state without any turbulent motion, this can be seen as the two different fluids forming now a single homogeneous fluid at rest, and the demixed stationary configuration where the two fluids are also at rest, but completly separated with the heavier fluid below the lighter. In other words this shows the existence of turbulent heteroclinic solutions emanating from the unstable interface configuration.

Brief comparison to experiments
There are numerous works addressing the Rayleigh-Taylor instability in different settings by means of experiments, numerical simulations and theoretical investigations for reduced models. For further reading we simply refer to the references given in the reviews [1,2,39,40]. At this point we only like to quickly compare our solutions for (1.1) with the results of the experiments carried out in [31] at Atwood number A ∼ 7.5 · 10 −4 .
First of all both the experiments and our solutions have a growth rate for the mixing zone like αgAt 2 . While the criterion of maximal initial energy dissipation selects α = 1 3 for our solutions, the actual constant observed in [31] is α = 0.07. Concerning self-similarity it is written in [31]: "The saturation of α at late time to a constant value of 0.07 suggests that the flow reaches self-similarity in these experiments." Another quantity we can easily compare is the ratio between dissipated energy and released potential energy. We will see in Section 4 (Remark 4.2) that up to an arbitrary small error our solutions show a ratio of D P rel. = 1 3 , while the ratio measured in [31] in dimension 3 is D P rel. = 0.49. We conclude that the solutions selected by maximal initial energy dissipation stand the comparison to actual experiments on a qualitative level, but it remains the interesting question if the gap between α = 1 3 vs. α = 0.07 and D P rel. = 1 3 vs. D P rel. = 0.49 can be improved in the future. In particular, it would be an interesting open problem to see if the measured values correspond perhaps to the optimization of some other mathematical quantity (other than the initial energy dissipation).

The role of the energy as a prescribed quantity
As in some other previous works of convex integration in fluid mechanics (e.g. [16,17,21]), there is a microscopic quantity which one has to prescribe in a continuous way in order to implement the convex integration. For instance, in the case of the homogeneous density incompressible Euler equations, this quantity was the kinetic energy 1 2 |v| 2 ; respectively in the case of the inhomogeneous incompressible Euler equations in [21] it was the quantity 1 2 ρ|v + gte n | 2 , which corresponded to the kinetic energy of a transformed system, and which can be seen as the kinetic energy of the original system plus a linear function of the momentum and the density. In our case the appropriate prescribed quantity which gives rise to the subsolutions mentioned before is 1 2 |v| 2 + 1 3 ρgAx n , i.e. the kinetic energy plus a fraction of the potential energy of the system.
In fact, both our convex integration strategy and the construction of our subsolutions from Section 4 can be carried out while prescribing the quantity 1 2 |v| 2 + ǫρgAx n , for any ǫ ∈ [0, 1]. However, setting for instance ǫ = 1, i.e. prescribing the total energy of the system, leads to solutions which are not admissible. The value ǫ = 1 3 is obtained by the process of maximizing the initial energy dissipation.
For further details, see the more detailed discussion in Sections 2.1, 2.2 and the construction in Section 4. [21] In [21] the authors together with L. Székelyhidi have addressed the Rayleigh-Taylor instability for the inhomogeneous incompressible Euler equations (1.4). In the aforementioned paper we obtained the existence of admissible turbulently mixing solutions for a sufficiently high density ratio ρ + ρ − , which translates to the Atwood number A being in the so-called "ultra high" range, A ≥ 0.845, i.e. far away from the Boussinesq range.

Comparison to
The proof there also relied on the explicit computation of the relaxation and convex integration within the Tartar framework. The computations for the convex hull in Section 3.3 resemble the computations done in [21]. While in the aforementioned paper a transformation of (1.4) onto an accelerated domain could be used in order to fit the system exactly into the Tartar framework (by which we mean that the gravity term in the momentum equation disappeared), here we can no longer use this transformation due to the Boussinesq approximation, and instead construct localized plane waves for an inhomogeneous linear system, see Section 3.1.
However, the main difference to [21] is the way subsolutions are constructed and selected. In [21] we reduced the relaxed system to a conservation law, which has some similarities to the conservation law appearing in [30] in a different approach to relax the incompressible porous media equation, and picked the unique entropy solution as our subsolution density profile. The subsolution found this way is selfsimilar and admissible for big enough A. Here instead, we consider the whole zoo of self-similar subsolutions and set up a variational problem whose unique minimizer corresponds to the subsolution maximizing the initial energy dissipation.
As illustrated in the previous subsection, contrary to [21] strictly speaking we do not provide one relaxation of the nonlinear system, but several relaxations, which differ in the amount of allowed turbulent behaviour in the local energy density, see Section 2.2.
The discussion of possible long time limits in Section 4.3 for the subsolution is not part of [21].

Outline of the paper
In Section 2 we formulate our results concerning the relaxation of (1.1) and the investigation of subsolutions in a precise way. Section 3 contains the steps needed to carry out convex integration in the Tartar framework and Section 4 contains the construction and selection of self-similar subsolutions.

Statement of results
Let Ω ⊂ R n be a bounded domain and T > 0. Our notion of solution to system is as follows.
Observe that the definition of v being weakly divergence-free includes the noflux boundary condition. Moreover, for a smooth vectorfield v the condition ρ ∈ {±1} automatically holds true, because then the density is transported along the flow associated with v, but for weaker notions of solutions this property in general is lost, see for example [28]. Furthermore, a (in general distributional) pressure p can be recovered from (ρ, v) as in the case of the homogeneous Euler equations, see [37].
The local energy density function E ∈ L 1 (Ω × (0, T )) associate with a weak solution (ρ, v) reads Indeed, testing a sufficiently smooth solution of (1.1) with v one sees that the total energy Ω E(x, t) dx is independent of t. However, this property in general fails to be true for weak solutions of Euler type equations, see [16] for Euler and [10] for the Boussinesq system. In order to rule out unphysical solutions due to an increase in energy and in view of the weak-strong uniqueness principle in various equations in fluid dynamics [38] we require the solutions to satisfy the following admissibility condition. |v 0 (x)| 2 + ρ 0 (x)gAx n dx for a.e. t ∈ (0, T ).

The relaxation
Next we will reformulate equation (1.1) as a differential inclusion and state its relaxation. Let S n×n be the set of all symmetric n × n matrices, S n×n 0 ⊂ S n×n the subset of matrices with vanishing trace and id ∈ S n×n be the identity matrix. We also write λ max (S), λ min (S) for the maximal, minimal resp., eigenvalue of S ∈ S n×n , and the trace free part of S is denoted by S • := S − 1 n tr(S) id. Consider on Ω × (0, T ) the linear system , which is affine linear in r. A brief discussion on possible choices of e and some general constraints can be found in Section 2.2 below. Now if z : Ω × (0, T ) → Z is a weak solution of (2.3), (2.4) to some initial data (ρ 0 , v 0 ) as in (2.1), see Definition 2.3 below for the precise definition, and if for almost every (x, t) ∈ Ω × (0, T ) there holds z(x, t) ∈ K (x,t) , then (ρ, v) defines a solution to the original equation (1.1) in the sense of Definition 2.1 for the same initial data and with energy density function given by Conversely, if (ρ, v) with associated pressure p is a weak solution in the sense of Definition 2.1, then z = ρ, v, ρv, (v ⊗ v) • , p + 1 n |v| 2 is a weak solution of (2.3), (2.4) and z pointwise a.e. takes values in the set K (x,t) defined with respect to the function e(x, t)[r] = 1 n |v(x, t)| 2 . For the relaxation of (2.3), (2.5) let Z 0 := { z ∈ Z : ρ ∈ (−1, 1) }, as well as In the course of the article we will show that U (x,t) is the interior of the convex hull of K (x,t) . That in particular means that if (ρ k , v k ) k∈N is a sequence of weak solutions with v k ∈ L ∞ (Ω × (0, T ); R n ) and such that the following convergences hold true ) and 1 n |v k | 2 → e in L ∞ (Ω×(0, T )), then there exists a pressure p, such that (ρ, v, m, σ, p) is a weak solution of (2.3), while pointwise a.e. taking values in With the help of the linear system (2.3) and the sets (2.7) we are ready to formulate the notion of subsolutions to (1.1), as well as our general convex integration result. Doing this the following projection turns out to be convenient: for , and if there exists an open set U ⊂ Ω × (0, T ), such that the two restricted maps U ∋ (x, t) → π(z(x, t)) ∈ π(Z) and U × R ∋ (x, t, r) → e(x, t)[r] ∈ R are continuous, and if there holds z(x, t) ∈ U (x,t) for all (x, t) ∈ U , as well as z(x, t) ∈ K (x,t) for a.e. (x, t) ∈ Ω × (0, T ) \ U . The open set U is called the mixing zone of z, and in analogy to solutions we call the subsolution admissible provided Before formulating our convex integration theorem we like to point out the following observation, which follows from Lemma 8 in [17].

Remark 2.4. Without loss of generality the ρ-component of any subsolution or solution is contained in
More precisely, [17,Lemma 8] gives ρ ∈ C 0 ((0, T ); L 2 w (Ω)), but looking into the proof one sees that the functions in [17, equation (90)] can be uniquely extended Observe also that outside the mixing zone U the components (ρ, v) of a subsolution z already solve the Euler-Boussinesq equation (1.1).
Theorem 2.5. Let z = (ρ, v, m, σ, p) be a subsolution associated with e and initial data Moreover, among these solutions one can find a sequence (ρ k , v k ), k ∈ N, such that Remark 2.6 (Admissibility). Observe that by Remark 2.4 and the assumption on e 1 the integral on the left-hand side in Thm. 2.5 c) defines a continuous function on [0, T ]. Moreover, for a.e. t ∈ (0, T ) the energy difference between the subsolution and the solutions is precisely given by this term, i.e.
for a.e. t ∈ (0, T ). In particular, if the subsolution is admissible with strict inequality in (2.10) for a.e. t ∈ (0, T ), then by a suitable choice of error function δ(t) one sees that property c) implies the admissibility of the induced solutions (ρ sol , v sol ) in the sense of Definition 2.2.
In that sense at every t ∈ [0, T ] the subsolution density ρ(·, t) can be seen as a coarse grained or averaged density of the induced solutions ρ sol (·, t), whose turbulent nature is illustrated by means of the mixing at every time slice property d).
The proof of Theorem 2.5 will be carried out in Section 3 and is based on the convex integration methods introduced by De Lellis and Székelyhidi in [16,17] and its refinements in [6,12]. In particular looking at [6] one could in addition also add the "linearly degraded macroscopic behaviour" to the list of properties of the solutions in Theorem 2.5. Moreover, if one is interested in the notion of admissibility at every time, by which we mean that the inequality in Definition 2.2 holds for all t ∈ [0, T ] instead of a.e. t ∈ (0, T ), one can use the convex integration strategy from [6,17] based on a "shifted grid", which is not used here.

Choices for e(x, t)[r]
In order to have inside the mixing zone U of a subsolution a non-empty interior of the convex hull U (x,t) we need In general e(x, t)[r] has to be non-negative a.e., because this expression coincides up to a positive factor with the kinetic energy of the solutions.
Besides the above conditions one can a priori use for e(x, t)[r] any function of the type e(x, t) = e 0 (x, t) + e 1 (x, t)r with e 0 , e 1 continuous on U , but in fact we will only consider such e with With this choice the solutions obtained by Theorem 2.5 will have a kinetic energy a.e. given by This means that besides the continuous part n 2 e 0 (x, t), which can be seen as a non turbulent or averaged part, the kinetic energy density of the solutions absorbs a certain fraction, given by n 2 ε ∈ [0, 1], of the turbulent oscillations in the potential energy density gAx n ρ sol (x, t).
A priori also ε can be a function depending on (x, t), but we will mostly stick to constant ε, except for Section 4.3.

Subsolutions
Our second main result addresses the construction and selection of subsolutions associated with the initial data ρ 0 = sgn(x n ), v 0 ≡ 0. We consider the problem on an n-dimensional box Ω = (0, 1) n−1 × (−L, L), L > 0, n ≥ 2 and focus on self-similar subsolutions. For the precise definition let F denote the set of all In Section 4.1 we will prove the following lemma.
We refer to these subsolutions as self-similar subsolutions. Considering solutions with v ≡ 0, m i ≡ 0, 1 ≤ i ≤ n − 1 and independent of x 1 , . . . , x n−1 reflects the interpretation of the subsolution as an (x 1 , . . . , x n−1 )-averaged solution. Moreover, we will see that the symmetry condition on f is needed for the existence of self-similar subsolutions for the Boussinesq system. In contrast the subsolution constructed in [21] for the Euler system without Boussinesq approximation is also self-similar, but the profile f is not symmetric.
Note that the associated mixing zone is given by U a := { (x, t) : |x n | < a(t) }. Note also that at this point the subsolutions are not necessarily admissible.
In order to investigate the admissibility let z = z f,a,ε be a self-similar subsolution and define the functionẽ f,a,ε : Hence by this definition, Theorem 2.5 c) and Remark 2.6 the subsolution z f,a,ε induces mixing solutions whose total energy Ω E sol (x, t) dx for a.e. t ∈ (0, T ) is arbitrarily close to Note that if z f,a,ε is admissible, then E f,a,ε (t) ≤ E(0) for a.e. t ∈ (0, T ), where E(0) = gAL 2 is the initial energy associated with (1.3). In order to evaluate the initial loss of energy define for k = 0, . . . , 4 the functionals whenever the limits exist. We have the following small time selection of a selfsimilar subsolution.
We will see that for f (y) = y, a(t) = 1 3 gAt 2 and ε = 2 3n there holds as long as a(t) ≤ L, i.e. for all t ∈ 0, 3L gA . Next we will formulate the two statements concerning the extension of the subsolution to all times. We like to emphasize that for the extensions we no longer use a selection criterion, instead the constructions contain several choices and for now are only done to illustrate possible options for the long-time behaviour.
Proposition 2.10. The minimizing subsolution from Theorem 2.9 with o(t 2 ) = 0 can be extended in an admissible manner to Ω × (0, +∞) such that it converges to the fully mixed, isotropic state z ≡ 0 as t → +∞ and such that also the associated kinetic energy n 2 e(x, t)[ρ(x, t)] converges to 0. Proposition 2.11. There exists T end ∈ 3L gA , +∞ such that the minimizing subsolution from Theorem 2.9 with o(t 2 ) = 0 can be extended in an admissible manner to Ω × (0, T end ), and at T end it reaches the stable configuration ρ = −ρ 0 , In fact, both subsolutions are not only admissible, but satisfy the strong energy inequality, which means that the total energy Ω E sub (x, t) dx is monotone decreasing w.r.t. time.
Moreover, in the first case the subsolution satisfies z(x, t) ∈ U (x,t) for every x ∈ Ω and t > 3L gA while the closure of the hull U (x,t) collapses as t → +∞ to the set [−1, 1] × {0} × {0} × {0} × R ⊂ Z due to the decay of kinetic energy. Thus technically the mixing zone is unbounded here. In the second case we have that z(x, T end ) actually is a solution, i.e. z(x, T end ) ∈ K (x,T end ) for a.e. x ∈ Ω. Clearly we can extend this subsolution to all times by z(·, t) = z(·, T end ) for all t > T end .

Convex integration via the Tartar framework
To prove our main result, we will use a version of the Tartar framework, originally introduced in the context of compensated compactness [35], for differential inclusions when the set of nonlinear constraints is not constant (c.f. e.g. [6,12,17]).
The general strategy of convex integration in the Tartar framework relies on the idea that if one can find a weak solutionz of (2.3) which instead of taking values in K (x,t) satisfiesz(x, t) ∈ int K co (x,t) , then one may deduce the existence of (infinitely many) solutions z of (2.3), which are nearz in the weak sense while satisfying z(x, t) ∈ K (x,t) a.e., by adding some specially constructed perturbations toz. The perturbations rely on localized plane waves as basic building blocks.

Localized plane waves
such that the wave cone associated with (2.3) can be written as . This allows us to construct solutions which oscillate in the directionz. Note that the condition (ρ,v) = 0 allows us to exclude the degenerate case when ξ = 0, which would correspond to having only oscillations in time.
Let us define a restricted wave cone which also eliminates oscillations only in space, i.e.
In Lemma 3.2 below we construct localized plane wave-like solutions for (2.3) associated withz ∈ Λ ′ . In order to see that it is enough to consider Λ ′ instead of Λ we first show the following density lemma.
We define the following sequence. For N ≥ 1 let Here and in forthcoming formulas the definition ofσ N andp N is understood in the sense that the symmetric matrix on the right hand side is split into its trace free part and its trace. It is easy to check that ξ, − 1 Recall the definition of the projection π : Z → R × R n × R n × S n×n 0 from (2.8). We write d for the euclidian distance function.
Lemma 3.2. There exists C > 0 such that for anyz ∈ Λ ′ , there exists a sequence z N ∈ C ∞ c (B 1 (0); Z), where B 1 (0) ⊂ R n × R, solving the linear system (2.3) and satisfying Proof. We will construct the desired sequence of solutions as a sum of two sequences z N =ẑ N +z N , whereẑ N will be a localized plane wave for the usual Euler equations determining up to a small deviation v N , σ N and p N , whilez N will take care of ρ N and m N .
For higher dimensions, one may proceed analogously, the details are left to the reader. This concludes the first step of our construction.
Step 2. The potential forρ andm. We will show that there exists a constantC > 0 independent ofz, and a sequencez N ⊂ C ∞ c (B 1 (0); Z) of solutions of (2.3), such that It is clear that C := min C ,Ĉ and z N :=ẑ N +z N then satisfies the properties stated in the lemma. For the existence ofz N observe that for any Ψ ∈ C ∞ (R n ×R; S n×n ) the functioñ z defined byρ := ∂ t div div Ψ,ṽ := gA∂ xn div Ψ − gAe n div div Ψ, satisfies equation (2.3).

Perturbing along sufficiently long segments
In this subsection we prove that the wave cone Λ is large with respect to K (x,t) , in the sense that any two points in K (x,t) can be connected with a Λ-segment. Furthermore, this property automatically implies that any point in the interior of the convex hull of K (x,t) can be perturbed along sufficiently long Λ-segments. The set K (x,t) has been defined in (2.5). For simplicity of notation, for the rest of the subsection we will fix a point (x, t) ∈ Ω × [0, T ) and write K instead of K (x,t) .
Without having any further information on the convex hull of K, we can prove the following geometric lemma solely based on this property. Corollary 3.5. For any z ∈ int(K co ) there existsz ∈ Λ such that [z −z, z +z] ⊂ int(K co ) and |π(z)| ≥ 1 2N d(π(z), π(K)), where N = dim(Z) and d is the Euclidean distance on π(Z).
The proof is the same as those of Lemma 6 from [17], respectively Lemma 4.9 from [21], relying on Carathéodory's theorem and Lemma 3.4 above, therefore we omit it.

The convex hull
We now explicitly compute the full Λ-convex hull associated with the differential inclusion (2.3), (2.5), which in our case turns out to coincide with the usual convex hull. The definition of the Λ-convex hull (K ′ ) Λ of K ′ ⊂ Z can be recalled for example from [25].
Let us again fix (x, t) ∈ Ω × [0, T ) and write K instead of K (x,t) , U instead of U (x,t) and e[r] instead of e(x, t)[r] for r ∈ R. Recall the definition of U in (2.7) and of the functions T ± , Q in (2.6). Proposition 3.6. There holds K Λ = K co = U .
In Lemma 3.8 below we will see that the closure of U splits into Moreover, Lemma 3.11 actually shows that K ′ ± is the Λ-convex hull of the sets The proof of Proposition 3.6 is organized as in the corresponding Section of [21] for the inhomogeneous Euler equation and relies on Lemma 3.8 and 3.11.
Lemma 3.7. The function Q is convex.
Proof. Since Q is defined as the maximal eigenvalue of the matrix M(z), there holds Q(z) = sup where for every fixed ξ ∈ S n−1 the function g ξ : { z ∈ Z : ρ ∈ (−1, 1) } → R is given by We will show that every g ξ is convex, such that Q is convex as a supremum of convex functions. In order to prove the convexity of g ξ , ξ ∈ S n−1 fixed, we write Then it is enough to show that the function g : (−1, 1) × R 2 → R, Let us fix (ρ, x) ∈ (−1, 1) × R 2 and observe that A(ρ) is positive definite. Thus the restricted function g(ρ, ·) is strictly convex, or equivalently D 2 g(ρ, x)[0, y] 2 ≥ 0 for all y ∈ R 2 .
It therefore remains to show that D 2 g(ρ, x)[1, y] 2 ≥ 0 for all y ∈ R 2 . By the positive definiteness of A(ρ) we obtain

It turns out that in fact
, which shows the convexity of g. Indeed, differentiating Moreover, a straightforward computation shows The following lemma implies the inclusion K Λ ⊂ K co ⊂ U.
Lemma 3.8. The set U is convex and its closure U splits into U = K ′ − ∪ U 0 ∪ K ′ + . In particular K ⊂ U .
Proof. In Lemma 3.7 we have already shown that Q is a convex function. Using the basic triangle inequality one can directly check that T ± (z) < e[±1] define convex sets. Hence U is convex.
For the stated identity concerning U first of all observe that U 0 ⊂ U . Next we will show that K ′ ± ⊂ U . Let z * ∈ K ′ + for instance and take any z ′ ∈ K with ρ ′ = −1, as well as a sequence (ρ j ) j∈N ⊂ (−1, 1) with ρ j → 1. The element clearly converges to z * as j → +∞. Using z * ∈ K ′ + and z ′ ∈ K, ρ ′ = −1 one sees that as well as T − (z j ) = e[−1]. For the matrix M(z j ) we compute This shows that every z j and therefore also the limit z * is contained in U . The case z * ∈ K ′ − works analoguosly. Thus K ′ − ∪ U 0 ∪ K ′ + ⊂ U. Let now (z j ) j∈N ⊂ U be convergent to some z * ∈ Z. If ρ * ∈ (−1, 1), then it is clear that z * ∈ U 0 ⊂ U . Consider the case ρ * = 1. Since on U there holds |m ± v| < ne[±1](1 ± ρ), (3.10) it follows that m * = v * . Recall that e[±1] ≥ 0 from (2.11). Next we rewrite (3.11) and observe that by (3.10). Therefore lim Thus λ max (M(z j )) < e[ρ j ], j ∈ N and the continuity of the maximal eigenvalue function imply z * ∈ K ′ + . The same procedure again works for the other case ρ * = −1, such that the statement of the Lemma follows.
In terms of Proposition 3.6 it now remains to prove the inclusion U ⊂ K Λ . The proof of this inclusion will rely on the Krein-Milman theorem for Λ-convex sets [25,Lemma 4.16]. For this we discuss the following Λ-directions.
Proof. The proof is a straightforward adaption of Lemma 4.6 (ii),(iii) in [21] and therefore only sketched here. As a nontrivial element (ξ, c) ∈ R n × R in the kernel of M Λ (z(z)) one can take any ξ = 0 contained in the orthogonal complement of v(v) and set c = −m(z) · ξ.
The stated invariances can be verified directly. Note that for T ± it helps to rewrite whereas for M • identity (3.11) is useful.
As in [21] we callz(z) the Muskat direction associated with z, since it generalizes the density perturbation of the Muskat problem introduced in [33]. Also as in [21] we have the following lemma concering Euler type directions preserving the density. ,v = 0, there existsp ∈ R, such that for all λ ∈ R the vectorz λ := (0,v, λv,σ,p) belongs to Λ. Moreover, for all t ∈ R there holds Proof. See the proof of Lemma 4.6 (i),(iv) in [21].
We have the following results concerning Λ-extreme points of U . Recall that π : Z → R × R n × R n × S n×n 0 is the projection from (2.8).
Lemma 3.11. The set π(U) is bounded by a constant depending only on e and the dimension n. Moreover, for every z ∈ U \ K there existsz ∈ Λ \ {0}, such that z ±z ∈ U .
Proof. Let z ∈ U. Clearly |ρ| ≤ 1 and the two inequalities (3.10) imply a bound on v and m in terms of e and n. Using (3.11), (3.12), we obtain that M(z) + σ is also bounded by means of e and n. In consequence we obtain |tr M(z)| ≤ c(e, n). Since the trace is bounded and λ max (M(z)) = Q(z) < e[ρ], using that z ∈ U, we get a corresponding bound on the whole spectrum of M(z). Hence, M(z) + σ and M(z) are both uniformly bounded, and therefore |σ| ≤ c(e, n). This proves that π(U) is bounded.
Next we turn to the perturbation property. Let z ∈ U \ K and recall from Lemma 3.
there exists an Euler type directionz +1 , i.e.m =v, from Lemma 3.10, such that z ±z +1 ∈ K ′ + \ K. The proof is the same as in [17] and [21,Lemma 4.8] and therefore omitted. Similar for z ∈ K ′ − \ K. It remains to look at z ∈ U 0 . Let us first check in which cases we can use the associated Muskat directionz =z(z) from Lemma 3.9. By this Lemma the two inequalities T ± (z + tz(z)) ≤ e[±1] remain true for all t ∈ (−1 − ρ, 1 − ρ). Furthermore, a straightforward computation shows that Q(z) − e[ρ] can be rewritten as Using Lemma 3.9 once more we therefore obtain .  Once again proceeding as in [17] or [21,Lemma 4.8], one may easily prove that there exists such an Euler direction which does not effect the maximal eigenvalue of M(z), i.e. such that Q(z + tz) = Q(z) = e[ρ] for small enough |t|. The last condition needed for z + tz +1 ∈ U follows from the continuity of T + , i.e., for all |t| small enough one has T + (z + tz) − e[+1] < T − (z) − e[−1] ≤ 0.
Proof of Proposition 3.6. From Lemma 3.8 one obtains K Λ ⊂ K co ⊂ U, while Lemma 3.11 implies that the Λ-extreme points of the up to the p-component compact set U are contained in K. The inclusion U ⊂ K Λ follows from the Krein-Milman theorem for Λ-convex sets, cf. [25,Lemma 4.16].

Continuity of constraints
We have the following result regarding the continuity of the nonlinear constraints K  Proof. The boundedness of (x,t)∈U π(K (x,t) ) follows from Lemma 3.11 and the boundedness of e.
First of all observe that by the continuity of e there exists δ ∈ (0, ε) such that Furthermore, from (3.13) it follows that This way we have shown that for any y ′ ∈ B δ (y) and any z ∈ K y there exists z ′ ∈ K y ′ ∩ B cε (z) for some c > 0 depending only on the dimension n. Using the symmetry of this construction, one can similarly prove that for any z ′ ∈ K y ′ there exists z ∈ K y such that |z − z ′ | < cε. As illustrated above we then conclude d H (π(K y ), π(K y ′ )) < cε via [12, Lemma 3.1].

Conclusion
We have now collected all the ingredients for the proof of Theorem 2.5, which follows by the known convex integration procedures in the Tartar framework [16,17] and its refinements [6,12]. We refrain from formulating another version of the Tartar framework exactly taylored to our needs and instead only point out the small modifications that need to be done in the existing convex integration theorems in order to conclude Theorem 2.5. We begin with the functional setup. Let D := Ω × (0, T ). Fix a function e : D × [−1, 1] → R, a subsolutionẑ = (ρ,v,m,σ,p) with initial data (ρ 0 , v 0 ) and mixing zone U , as well as an error function δ : [0, T ] → R as stated in Theorem 2.5. Define X 0 to be the set of all functions π(z) = (ρ, v, m, σ), such that • z = (ρ, v, m, σ, p) is a subsolution for e, (ρ 0 , v 0 ) and with the same mixing zone U , in the sense of Definition 2.3, • z =ẑ a.e. on D \ U , • there exists C(z) ∈ (0, 1), such that for all t ∈ [0, T ] there holds Recall from Section 2.1 that e(x, t)[r] = e 0 (x, t)+re 1 (x, t) with L ∞ functions e 0 , e 1 , where e 1 is additionally of class C 0 ([0, T ]; L 2 (Ω)).
It therefore remains to improve the convergence of the ρ-component from ρ k ⇀ ρ weakly in L 2 (D) to ρ k → ρ in C 0 ([0, T ]; L 2 w (Ω)) and to show that the functions (π(z k )) k∈N satisfy (3.14) for all k big enough. However, the improved convergence follows from Remark 3.3 by using cylinders instead of balls in the proof of [12,Lemma 2.4]. Finally, the fact that the sequence (π(z k )) k∈N satisfies property (3.14) for k big enough follows as in Step 3 of the proof of [6, Proposition 3.1]. Indeed, since z ∈ X 0 , it is enough to fix C ′ (z) ∈ (0, 1 − C(z)) and to show for all t ∈ [0, T ] and all k sufficiently large. Since by construction ρ k = ρ outside a compact subset of the mixing zone U , hence outside a set contained in [t 0 , t 1 ] × Ω for some 0 < t 0 < t 1 < T , it is enough to show ∀t ∈ [t 0 , t 1 ] : where f (x, t) := n 2 e 1 (x, t) + gAx n and δ 0 := inf { δ(t) > 0 : t ∈ [t 0 , t 1 ] } > 0. But the latter inequality holds true for big enough k due to the uniform continuity of the map [0, T ] ∋ t → f (·, t) ∈ L 2 (Ω), the uniform bound on ρ k (·, t) L 2 (Ω) and the convergence ρ k → ρ in C 0 ([0, T ]; L 2 w (Ω)). Proof of Theorem 2.5. Having Lemma 3.13 at hand we can prove as in [12] or [21] that J −1 (0) is contained in the set of continuity points of I, where I, J were defined in (3.15), (3.16). Since I is Baire-1, this shows that Observe also that if π(z) ∈ J −1 (0), then (ρ, v) is a weak solution of (1.1), (1.2) satisfying properties a) and b) of Theorem 2.5.
Concering property Thm. 2.5 c), approximation by elements from X 0 with respect to d X shows that any element π(z) from X satisfies

Self-similar subsolutions
In this section we prove Lemma 2.8. Recall the definitions of F and A from (2.13), respectively (2.14). For any choice of a profile f ∈ F and a growth rate a ∈ A one can check that x n a(t) e n for x n ∈ (−a(t), a(t)), as well as σ(x, t) = 0 for |x n | ≥ a(t), p(x, t) = −σ nn (x, t) − gA x n ∈ (−a(t), a(t)) } will be the mixing zone. Concerning the pointwise constraints we definẽ for (x, t) ∈ U andẽ ε (x, t) = εgA |x n | for (x, t) ∈ D \ U . Here δ : D → (0, +∞) is a continuous, even, positive and typically small function guaranteeing the inequalities (4.1) to hold in a strict sense. Indeed the first three conditions in (4.1) hold by definition of ρ andẽ ε . For the last inequality we have . This concludes the proof of Lemma 2.8.

Admissibility and maximal initial energy dissipation
Instead of investigating all admissible subsolutions emanating from Section 4.1, we will focus on the one that is selected by asking for maximal initial energy dissipation.
For (f, a, ε) ∈ F × A × 0, 2 n observe that the total energy at time t > 0 of the induced subsolution can be choosen arbitrarily close to E f,a,ε (t) defined in (2.15), which for admissibility has to be less than the initial energy E(0) = Ω gA |x n | dx.
Using this, the definitions of ρ, m from Section 4.1, the transformation x n = a(t)y and the symmetry of f one sees that the difference of the energies can be computed by the following integrals Concerning the well-definedness observe again that for all f ∈ F the quotient F (y) 1−f (y) has a finite limit as y → 1.
For a given profile f ∈ F and a growth rate a ∈ A one can via the above formula simply check by hands the admissibility of the induced self-similar subsolution.
gA , the choices f (y) = y, a(t) = 1 3 gAt 2 and ε = 2 3n give rise to a subsolution on Ω × (0, T ) with In particular this implies that the subsolution is admissible for small δ(x, t).
Remark 4.2. The released potential energy of the subsolution above at time t ∈ [0, T ) is given by Therefore the ratio between dissipated and released energy is 1 3 . Besides the fact of being a simple example, it turns out that these choices for f, a, ε maximize the initial energy dissipation.
Proof of Theorem 2.9. In the formula for the energy difference let us abbreviate the two terms among which the maximum is taken, i.e., set Estimating the maximum from below by the convex combination Observe that (1 − f (y))y dy > 0, (4.6) such that the required admissibility implies and thereforeȧ(0) = 0, J 1 (f, a, ε) = 0. Since now a(t) = 1 2ä (0)t 2 + o(t 2 ) as t → 0 the admissibility also implies J 2 (f, a, ε) = J 3 (f, a, ε) = 0. This proves the first part of the Theorem.
The lowest order for which the initial energy dissipation rate is not necessarily vanishing is 4. There holds In Lemma 4.3 below we will show that the functionalJ : F × [0, +∞) → R has a unique global minimum in f (y) = y andä(0) = 2 3 gA. It follows that for any (f, a, ε) ∈ F × A × 0, 2 n leading to an admissible subsolution there holds Proof. First of all observe that for fixed f ∈ F the functionJ (f, ·) : [0, +∞) → R has a unique minimum in c 0 (f ) = 2 3 gA I 2 (f ) and it remains to show for any f ∈ F \ {id}. Note that for f = id there holds equality, since I 1 (id) = 1 6 , I 2 (id) = 1 6 . Let us rewrite Since Since also I 1 (f ) is positive, we see that (4.8) holds true provided In order to proveĴ (f ) > 0 for f ∈ F \ {id} we write f = id +ϕ with ϕ = 0, such that Thus in terms of f and the L 2 (0, 1) inner product and norm we havê The next (and last for this subsection) lemma implies thatĴ(f ) > 0 for all f ∈ F \ {id}, which allows us to conclude the proof of Lemma 4.3. Proof. We set F 0 := { f ∈ L 2 (0, 1) : |f | ≤ 1 a.e. }, which is the closure of F 0 with respect to · L 2 (0,1) , and observe thatĴ (f ) ≥ − 4 3 for f ∈ F 0 . Now let (f n ) n∈N ⊂ F 0 be a minimzing sequence forĴ. Since F 0 is bounded and convex there exists f * ∈ F 0 with f n ⇀ f * in L 2 (0, 1) along a subsequence. By the weak lower semicontinuity of the norm and since Thus the minimum ofĴ : F 0 → R is achieved at f * . Now there are two cases to consider: f * ∈ F 0 and f * ∈ F 0 \ F 0 . In the first case f * ∈ F 0 one can check that f * is a critical point ofĴ considered as a map from all of L 2 (0, 1) to R.
It is clear thatĴ : L 2 (0, 1) → R is smooth and a quick computation shows that the gradient is given by Thus for a critical point ofĴ there holds S(f ) = 3 2 and Plugging this identity into the definition of S(f ) one obtains that or equivalently S(f ) ∈ { 0, 1 }. ThusĴ : L 2 (0, 1) → R has exactly two critical points in f = id and f = −2 id, and only f = id is contained in F 0 . Consequently if the minimum ofĴ |F 0 is achieved at f * ∈ F 0 , then f * = id andĴ |F 0 ≥Ĵ(id) = 0.
If we assume that id is not minimizingĴ |F 0 , then any minimizer f * lies in F 0 \ F 0 and satisfiesĴ(f * ) < 0. Without loss of generality we can assume f * ≥ 0 and f * to be nondecreasing, otherwise we replace f * by the monotone increasing rearrangement of |f * |, which only decreasesĴ, cf. (4.10). These two properties together with f * ∈ F 0 \ F 0 imply that there exist f 0 ∈ F 0 and a ∈ [0, 1), such that for a.e. y ∈ (0, 1). Here X denotes the indicator function and for a = 0 this expression is understood as f * = X (0,1) . In a straightforward way one sees that such thatĴ (f * ) = a 2Ĵ (f 0 ).
Since by assumptionĴ(f * ) < 0, this equality implies a ∈ (0, 1) andĴ(f 0 ) <Ĵ(f * ), which tells us that f * can not be a minimizer ofĴ |F 0 . Due to this contradiction we conclude that the infimum ofĴ |F 0 is achieved at f * = id. Finally the strict inequalityĴ(f ) > 0 for f ∈ F 0 \ {id} follows from the fact that id is the only critical point ofĴ : L 2 (0, 1) → R lying in F 0 .

Beyond small-time behaviour
While the subsolution constructed in the previous subsection focused on minimizing the initial energy dissipation, one could also be interested in the long-time behaviour of such subsolutions. In particular, how can the subsolution be continued after a reaches L, i.e. the mixing zone touches the upper boundary. There are two long-time states which are of interest, namely the one where both the density and the momentum are vanishing everywhere (hence there are no longer two different density fluids, but only one completely mixed fluid), and the configuration −ρ 0 , where the higher density fluid occupies the lower half of the domain, respectively the lower density fluid occupies the upper half (i.e. gravity demixes the two fluids in the long term). We will show that both of these configurations can be achieved.
To conclude the proof of Proposition 2.10, observe that the limit of the subsolution as t → +∞ is identically zero, and δ can be chosen such that the energy of the system also decays to zero in the limit at +∞. Remark 4.5. Since the kinetic energy of the solutions associated with the constructed subsolution goes to 0 as t → +∞, any turbulent motion, in fact any motion, will vanish as t → +∞. Note that one could have made the same construction while keeping ε = 2 3n constant and still have obtained an admissible subsolution. However, the associated energy as t → +∞ would not vanish, which would imply that there is still some turbulence at infinite time.

Demixing in finite time
Let us now construct an example of a different admissible continuation past the time when the mixing zone touches the upper boundary, one where first the density profile is rotated by 180 degrees, and then the mixing zone shrinks until the stable configuration −ρ 0 is reached.
Proof of Proposition 2.11. We will do this in two steps.
We claim that T max < +∞ and r as a function extends continuously to [T 0 , T max ] with r(T max ) = − 1 L . Assume to the contrary that T max = +∞, then r(t) > − 1 L for all t ≥ T 0 . But now integrating (4.12) for t ∈ (T 0 , T max ) and using that I is decreasing, one has the contradiction as t → +∞. Hence T max < +∞ and then necessarily lim t→Tmax r(t) = − 1 L , because the orbit r([T 0 , T max )) is bounded from above due to the monotonicity of r. We therefore setT := T max .
Next due to I 1 L = 1 3 L 5 and (4.12) it is easy to see thatṙ (T 0 ) = −2 gA 3L 3 . Finally, let us show that the associated corrected total energy function E r is decreasing, to conclude the admissibility on [T 0 ,T ] of our subsolution. For this, we first show that one has r(t)ṙ(t) 2 ≤ 4gA 3L 4 , so that inẽ ε indeed the first term under the maximum is selected for x n ≥ 0 and thus (4.11) holds. Once again, this follows from (4.12) by using the monotonicity of I and r: Since the corrected total energy function E r is then given by formula (4.11), we may plug (4.12) into (4.11) to further obtain that E r (t) = 1 9 gAL 2 + 4 9 gAL 3 r(t) + 1 3 gAL 2 = 4 9 gAL 2 (1 + Lr(t)) , which is clearly decreasing since r is decreasing. This concludes the construction for the rotation of the profile.
Step 2: Shrinking of the mixing zone.