A New Approach to the Rayleigh–Taylor Instability

In this article we consider the inhomogeneous incompressible Euler equations describing two fluids with different constant densities under the influence of gravity as a differential inclusion. By considering the relaxation of the constitutive laws we formulate a general criterion for the existence of infinitely many weak solutions which reflect the turbulent mixing of the two fluids. Our criterion can be verified in the case that initially the fluids are at rest and separated by a flat interface with the heavier one being above the lighter one—the classical configuration giving rise to the Rayleigh–Taylor instability. We construct specific examples when the Atwood number is in the ultra high range, for which the zone in which the mixing occurs grows quadratically in time.


Introduction
We study the mixing of two different density perfect incompressible fluids subject to gravity, when the heavier fluid is on top. In this setting an instability known as the Rayleigh-Taylor instability forms on the interface between the fluids which eventually evolves into turbulent mixing. For an overview of the investigation of this phenomenon originating in the work of Rayleigh [30] in 1883 we refer to the articles [1,3,4,8,38,39].
The mathematical model (see for example Section 6.4 of [26]) is given by the inhomogeneous incompressible Euler equations ∂ t (ρv) + div (ρv ⊗ v) + ∇ p = −ρge n , div v = 0, ∂ t ρ + div (ρv) = 0, (1.1) which we consider on a bounded domain ⊂ R n , n 2 and a time interval [0, T ). Here ρ : × [0, T ) → R denotes the fluid density, v : × [0, T ) → R n is the velocity field, respectively p : × [0, T ) → R is the pressure, g > 0 is the gravitational constant and e n is the n-th Euclidean coordinate vector. Compared to the homogenous density case, ρ ≡ 1, the solvability of the Cauchy problem of (1.1) for a general non-constant initial density distribution is more delicate even in the planar case; see Section 6.4 of [26]. Results concerning the local wellposedness have only been obtained under sufficiently strong regularity assumptions on the initial density; see [13][14][15]37] and references therein. However, since we are interested in the mixing of two different fluids, our initial data does not fall into the classes considered in [13][14][15]37].
More precisely, we consider (1.1) together with initial data v 0 : → R n , ρ 0 : → R satisfying div v 0 = 0 and ρ 0 ∈ { ρ − , ρ + } almost everywhere (1.2) with two fixed values ρ + > ρ − > 0. In fact our main focus lies on the flat unstable initial configuration v 0 ≡ 0 and ρ 0 (x) = ρ + when x n > 0, ρ − when x n 0, (1.3) giving rise to the Rayleigh-Taylor instability. The linear stability analysis of the flat interface has already been investigated in the article of Rayleigh [30] and for example can also be found in [2]. Regarding the nonlinear analysis, to the best of our knwoledge there has been so far no existence result of mixing solutions for the case of the discontinuous initial data (1.3).
In the spirit of the results by De Lellis and the 3rd author [16,17], for the homogeneous incompressible Euler equations, we develop a convex integration strategy for the inhomogeneous Euler system to prove the existence of weak solutions for the Cauchy problem (1.1), (1.3). Similarly to other unstable interface problems that have recently been attacked by means of convex integration, like the Kelvin-Helmholtz instability in [34] or the Muskat problem for the incompressible porous media equation in [11,33], we can interpret the "wild" behaviour of the weak solutions obtained this way as turbulent mixing. More precisely, we prove the existence of solutions with the following properties: For ρ + > ρ − > 0 define the Atwood number A = ρ + −ρ − ρ + +ρ − and the quadratic functions c ± : R → R, Let T > 0 and = (0, 1) × (−c − (T ), c + (T )) ⊂ R 2 . . The initial value problem (1.1), (1.3) has infinitely many weak admissible solutions (ρ, v) ∈ L ∞ ( × (0, T ); R × R 2 ) with ρ ∈ {ρ − , ρ + } almost everywhere and such that For the precise definition of weak admissible solutions we refer to Definitions 2.1, 2.2.
We would like to point out that the infinitely many weak solutions differ only in their turbulent fine structure, while they all have a continuous coarse grained density profileρ in common. The profileρ(x, t) =ρ(x 2 , t) can be seen as an x 1 -average of the solutions, cf. Remark 2.5 (a) below and [7,Remark 1.2], and is found as the entropy solution of a conservation law ∂ tρ + gt∂ x 2 G(ρ) = 0, which up to the factor t shows similarities to the conservation law appearing in Otto's relaxation for the incompressible porous media equation [29]. Further details and the explicit profileρ can be found in Section 6.
The condition that the density ratio ρ + /ρ − is larger than 4+2 √ 10 3 2 ≈ 11.845, implies that the Atwood number is in the so-called (for example [5]) "ultra high" range (0.845, 1). The main obstruction in establishing a similar result to Theorem 1.1 for a density ratio outside of this range comes from the fact that in our approach the weak admissibility condition on the solutions associated with the profileρ reduces to an algebraic inequality for the density ratio; see Section 6 for further details. The ultra high regime has been of great interest to the physics and numerics communities recently, as it has many applications in fields such as inertial confinement fusion, astrophysics or meteorology (see for example [5,18,25]). For instance the Atwood number for mixing hydrogen and air is 0.85 (see [25]).
A higher Atwood number implies higher turbulence, and compared to the low Atwood regime, one can not use the Boussinesq approximation (see for example [4,8,24]) to accurately model the phenomena. Compared to the homogeneous density case, where the turbulence is only due to mixing in momentum, here it is due to mixing both in momentum and in density, this "double mixing" is reflected also in our relaxation given in Section 2.
We note that up to our knowledge, our result is the first rigorous result leading to existence of weak solutions with quadratic growth in time for the mixing zone. It is also of interest that both numerical simulations and physical experiments predict a growth rate of the mixing zone like αAgt 2 , but there is considerable disagreement about the value of the constant α and its possible dependence on A (see [5,18,25]).
In future work we plan to further study the possibility of constructing solutions with different mixing zone growth rates, to investigate the optimality of the growth rates c ± in Theorem 1.1, and to explore more precisely their relation to the values from experiments and simulations.
Concerning convex integration as a tool in the investigation of unstable interface problems we have already mentioned the papers [11,33,34]. While [11] shows the non-uniqueness of solutions to the incompressible porous media equation, the paper [33] provides the full relaxation of the equation allowing to establish sharp linear bounds for the growth of the mixing zone in the Muskat problem. The knowledge of the relaxation also opened the door to further investigations of the Muskat interface problem, see [6,7,22]. We already mentioned the different relaxation approach for the incompressible porous media equation via gradient flow in [29], the unique solution of this relaxation approach turned out to be recovered as a subsolution in [33].
Another classical instability in fluid dynamics is the Kelvin-Helmholtz instability generated by vortex sheet initial data. Regarding this instability solutions with linearly growing mixing zone have been constructed in [34] based on the computations of the relaxation of the homogeneous Euler equations in [17].
There have also been some recent convex integration results for the compressible Euler [9,20,21] and the inviscid Boussinesq equation [10]. The approach used for the compressible Euler equations ultimately relies on reducing the problem to having a finite partition of incompressible and homogeneous fluids. In [27] the convex hull of the isentropic compressible Euler system has been computed, but so far not used for the construction of weak solutions via convex integration. In the Boussinesq approximation the influence of density variations is neglected in the lefthand side of the momentum equation (1.1). Moreover, the result in [10] addressing the existence of infinitely many weak solutions to a given initial configuration requires the initial density to be of class C 2 and the obtained weak solutions to this prescribed initial data are not admissible in the sense that they violate the energy inequality. We would like to point out that so far there have been no convex integration results relying on the full relaxation of the compressible Euler equations nor the inhomogeneous incompressible Euler equations, the latter will be done in this paper.
The paper is organized as follows: in Section 2 we present our main results, one regarding the convex integration of the inhomogeneous incompressible Euler equations regardless of initial data, and one regarding the existence of appropriate subsolutions in the case of a flat initial interface.
In Section 3 we prove that through an appropriate change of coordinates, which in fact corresponds to the way how actual experiments investigating the Rayleigh-Taylor instability are carried out [18,31,32], problem (1.1) can be recast as a differential inclusion. The differential inclusion fits in a modified version of the Tartar framework of convex integration, adapted from [17,33] to simultaneously handle the absence of the pressure from the set of constraints and the dependence of the set of constraints on (x, t) due to the prescribed energy density function.
In Section 4 we prove the ingredients of the topological framework, most importantly we calculate the -convex hull of the set of constraints, which forms the core of this paper.
In Section 5 we conclude the proof of our main convex integration result, while in Section 6 we construct appropriate subsolutions having the growth rates presented in Theorem 1.1.

Statement of Results
Let ⊂ R n be a bounded domain and T > 0. Our notion of solution to equation (1.1) on × [0, T ) is as follows: Note that the definition of v being weakly divergence-free includes the no-flux boundary condition. Moreover, the last condition automatically holds true when we deal with smooth solutions of (1.1), because then the density is transported along the flow associated with v, but for weaker notions of solutions this property does not necessarily need to be true, see for example [28]. Furthermore, given a weak solution, the (in general distributional) pressure p is determined up to a function depending only on time, as in the case of the homogeneous Euler equations, see [36].
As in the homogeneous case, one can associate with a weak solution (ρ, v) an energy density function E ∈ L 1 ( × (0, T )) given by Furthermore, for smooth solutions of (1.1) one can show that t → E(x, t) dx is constant. For weak solutions this necessarily does not need to be true. As in the case of the homogeneous Euler equations or hyperbolic conservation laws, in order to not investigate physically irrelevant solutions we require our weak solutions to be admissible with respect to the initial energy.
One main contribution of the present article is the relaxation of equation (1.1) viewed as a differential inclusion. For the formulation of the relaxation we need the linear system considered on × (0, T ) and with z = (ρ, v, u, S, P) taking values in the space Here S n×n 0 denotes the space of symmetric n × n matrices with trace 0. We will also write S n×n for the space of symmetric matrices, id ∈ S n×n for the identity and λ max (S), λ min (S) for the maximal, minimal resp., eigenvalue of S ∈ S n×n .
As usual, equations (2.1) will be complemented by a set of pointwise constraints. Let e : × (0, T ) → R + be a given function and define for (x, t) ∈ × (0, T ) the sets as well as the sets U (x,t) by requiring for z ∈ U (x,t) the following four inequalities to hold: where Note that by the definition of K (x,t) in (2.2) and by recalling that S has vanishing trace, every solution of (2.1) taking values in K (x,t) almost everywhere is a solution to the inhomogeneous Euler equations (1.1) with ρ ∈ { ρ − , ρ + } and associated energy which is equivalent to saying that Conversely, if we have a solution (ρ, v, p) of (1.1) with ρ ∈ { ρ − , ρ + } almost everywhere, we can introduce the variables u = ρv, S = ρv ⊗ v − 1 n ρ |v| 2 id, P = p + 1 n ρ |v| 2 to see that z = (ρ, v, u, S, P) will satisfy system (2.1) while pointwise taking values z( Since the pressure P does not play a role in the set of constraints K (x,t) , it is convenient to consider the following projection: for z = (ρ, v, u, S, P) ∈ Z we denote (2.6) Using the linear system (2.1) and the definition of U (x,t) we define relaxed solutions to (1.1) in the following way: Definition 2.3. (Subsolutions) Let e : ×[0, T ) → R + be a bounded function. We say that z = (ρ, v, u, S, P) : × (0, T ) → Z is a subsolution of (1.1) associated with e and the initial data (ρ 0 , v 0 ) ∈ L ∞ ( ) × L 2 ( ; R n ) satisfying (1.2) iff π(z) ∈ L ∞ ( × (0, T ); π(Z )), P is a distribution, z solves (2.1) in the sense that v is weakly divergence-free (including the weak no-flux boundary condition), We call U the mixing zone of z. Moreover, the subsolution is called admissible provided that We now can state the following criterion for the existence of infinitely many weak solutions: ) and satisfies (2.8) with strict inequality for every t ∈ (0, T ], then infinitely many of the induced weak solutions are admissible in the sense of Definition 2.2.

Remark 2.5. (a)
The second to last two statements justify to call U the mixing zone and to interpret the subsolution density ρ as a kind of coarse-grained or averaged density profile.
(b) The result carries over to the three-or higher-dimensional case by constructing suitable potentials analoguosly to [16], which is not done here, cf. Lemma 4.1.
The other parts of the proof, for example the computation of the -convex hull in Section 4.2, are carried out in arbitrary dimensions. (c) We will see later that the open set U (x,t) is indeed the convex hull of K (x,t) .
In particular we can conclude that weak limits of solutions are subsolutions in the following sense: Let (ρ k , v k ) k∈N be a sequence of essentially bounded weak solutions of (1.1) and define as before ). Assume further that there exists a continuous bounded function e ∈ C 0 ( × (0, T )), such that e k := 1 n ρ k |v k + gte n | 2 → e in L ∞ ( × (0, T )). Then z supplemented by a possibly distributional P is a weak solution of the linear Our second main result addresses the construction of subsolutions associated with the initial data (1.3). Clearly it only makes sense to consider this initial data on domains satisfying ∩ (R n−1 × {0}) = ∅.

Definition 2.6. (Rayleigh-Taylor subsolution)
We call a subsolution z of (1.1) a Rayleigh-Taylor subsolution (short RT-subsolution) provided the initial data is given by (1.3) and the subsolution is admissible with strict inequality in (2.8) for every t ∈ (0, T ).
If ρ + > 4+2 √ 10 3 2 ρ − , then there exists a RT-subsolution z which only depends on An explicit description of the subsolutions and further discussion can be found after the proof of Theorem 2.7 in Section 6. Observe that by combining Theorems 2.4 and 2.7 we arrive at the statement of Theorem 1.1.

Reformulation as a Differential Inclusion
The proof of Theorem 2.4 will rely on a version of the Tartar framework for differential inclusions (cf for example [12,17,35]), where instead of looking for weak solutions of a nonlinear problem, one looks for weak solutions of a first order linear PDE, satisfying a nonlinear algebraic constraint almost everywhere.
In order to reformulate (1.1) into such a framework, we first observe that one can get rid of the gravity in the momentum equation by considering the system in an accelerated domain. As mentioned earlier, this transformation corresponds to actual Rayleigh-Taylor experiments [18,31,32] where the instability is created by considering the stable configuration (light fluid above heavy fluid) and accelerating the surrounding container downwards.
To make this precise, let ⊂ R n be a bounded domain, T > 0 and set such that for t ∈ (0, T ) the slice is given by Let (μ, w, q) be a weak solution of on D for some suitable initial data satisfying (1.2) and with weak boundary condition for y ∈ ∂D(t). More precisely, the notion of weak solution to (3.1), (3.2) is understood as in Definition 2.1, except that now g = 0, in the momentum and continuity equation resp., and the weak formulation of div w = 0 including the weak boundary condition (3.2) becomes Then if we define y := x + 1 2 gt 2 e n and set Observe also that the transformation (3.4) gives a bijective correspondence between solutions of (1.1) and (3.1). Furthermore, the formal energy associated with (3.1) is given by the term for a function e : × [0, T ) → R + . Then the total energy E(x, t) associated with the original system (1.1) is precisely given by (2.5).
We can now reformulate (3.1) as a differential inclusion by considering on D the system where in analogy to the homogeneous Euler equations e : × (0, T ) → R + is given and for the sake of consistency we have denoted μ ± := ρ ± . We will understand weak solutions of (3.5) in the following sense: Definition 3.1. We say that z : D → Z is a weak solution of (3.5) with initial data π(z 0 ) ∈ L 2 ( ; π(Z )) iff π(z) ∈ L 2 (D; π(Z )), q is a distribution, w satisfies (3.3) and one has ). This way we have arrived at a reformulation of equation (1.1) as a differential inclusion. The process is summarized in the following statement.
and the associated energy E by (2.5).

The Ingredients of the Tartar Framework
The general strategy of the Tartar framework relies on the following steps: • finding a wave cone ⊂ Z such that for anyz ∈ , one can construct a localized plane wave associated with (3.5) oscillating in the direction ofz; • calculating the -convex hull of K (x,t) (denoted by K (x,t) ) and proving that one can perturb any element in its interior along sufficiently long -segments, provided that one is far enough from K (x,t) ; • deducing an appropriate set of subsolutions using K (x,t) and proving that it is a bounded, nonempty subset of L 2 (D).
In the following subsections we execute each of the above steps in the case of the differential inclusion (3.5), (3.6). Then we can conclude the proof of Theorem 2.4 in Section 5 by using the Baire category method (see [11,16,17,23,35]).

Localized Plane Waves
We begin with the construction of plane wave-like solutions to (3.5) which are localized in space-time. We consider the following wave cone associated with (3.5): It has the property that forz ∈ there exists η ∈ R n+1 \{0} such that every is a solution of (3.5). In Lemma 4.1 below we localize these solutions by constructing suitable potentials. Note that the condition (μ,m) = 0 serves to eliminate the degenerate case when the first n components of η vanish, that is when one is only allowed to oscillate in time.
Recall the projection operator π defined in (2.6).
Lemma 4.1. There exists C > 0 such that for anyz ∈ , there exists a sequence solving (3.5) and satisfying that Proof. We will only present the proof in the two-dimensional case, higher dimensions can be handled analogously to [16].
We start by observing that for any smooth functions ψ : Let S : R → R be a smooth function, N 1 andz ∈ with (μ,m) = 0. It follows from the definition of that there exists We then treat two cases. Case 1: c = 0 Note that in this case we also have ξ = 0, since ξ = 0 would imply (μ,m) = 0. We then set and we claim that Indeed, using (4.1), one has From here on, the localization is done in the standard fashion (for example as in [11,16]). We fix S(·) = − cos(·) and, for ε > 0, consider Case 2: c = 0 In this case we are not allowed to oscillate in time. However, we have ξ = 0, so we may also assume without loss of generality that |ξ | = 1. On the other hand, We set from where with similar calculations as in Case 1, we obtain that To handle the remaining terms (m,σ ,q), we introduce a different type of potential, as done for the homogeneous Euler equations, for instance in [16], Remark 2.
It can be checked through direct calculation that for any smooth function ω : implies thatD(ω) solves (3.5). Now, if we consider ω of the form for some constants a, b ∈ R, with S as before, we obtain that gives us from where, using (4.4), we get The localization is then done as in Case 1, by considering If ξ 1 = 0, then choosing a = k 3 gives us that However, it is easy to see that for any smooth function θ : R 2+1 → R,D(θ ) = (0, 0, ∇ ⊥ θ, 0, 0) also solves (3.5). Therefore, we may consider the potential given by we obtain that and using (4.3), we get that One may then localize this potential by the usual means in order to conclude the proof of the lemma.

The -Convex Hull
We now turn to the set of pointwise constraints K (x,t) , (x, t) ∈ D defined in (3.6). The -convex hull K (x,t) is defined by saying that z ∈ K (x,t) iff for allconvex functions f : Z → R there holds f (z) sup z ∈K (x,t) f (z ), see [23] for more details. In our case it turns out that the -convex hull is nothing else but the usual convex hull, see Proposition 4.2 below.
For the computation of the hull we drop the (x, t) dependence of the sets K (x,t) and consider a general set of pointwise constraints given by where 0 < μ − < μ + , e ∈ R + are given constants.
as well as the open set Proposition 4.2. The -convex hull of K coincides with the convex hull of K and is given by U , that is, K = K co = U .
The proof of Proposition 4.2 relies on Lemmas 4.4 and 4.8.

Lemma 4.3. The function Q is convex.
Proof. We write where for every fixed ξ ∈ S n−1 the function g ξ : Z 0 → R is given by We will show that every g ξ is convex. As a consequence Q is convex as a supremum of convex functions. In order to do this let us complement ξ ∈ S n−1 to a orthonormal basis (ξ, v 2 , . . . , v n ) of R n . Expressing w and m with respect to this basis one sees that it is enough to show that the function g : Thus the restricted function g(μ, ·) is 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

Now we claim that in fact
, which finishes the proof. Indeed, differentiation of the identity where Moreover, in a straightforward way one can check that which implies that

Lemma 4.4. The set U is convex and U
Proof. For μ ∈ (μ − , μ + ) the two conditions T + (z) < e, T − (z) < e can be rewritten as where c ± = ne . Using the basic triangle inequality one can check that the two conditions in (4.10) define a convex set. By Lemma 4.3 we already know that Q is a convex function. Hence we have shown that U is convex. Now we turn to the characterization of U . Clearly U 0 ⊂ U . Let us show that K + ⊂ U . The inclusion K − ⊂ U can be obtained in the same way. Let z * ∈ K + . Take any z ∈ K with μ = μ − and some sequence (μ j ) j∈N ⊂ (μ − , μ + ) converging to μ + . Define Clearly z j → z * as j → ∞. Since z * ∈ K + and z ∈ K − a short calculation shows Similarly we obtain T − (z j ) = e. In a third, slightly longer computation we plug z j into M, sort with respect to the terms w ⊗ w , w ⊗ w * , w * ⊗ w , w * ⊗ w * and find We conclude Q(z j ) = λ max (M(z j )) e. Hence every z j and therefore also the limit z * is contained in U . So far we know For the other inclusion consider (z j ) j∈N ⊂ U , z j → z * . The interesting case of course is μ * / ∈ (μ − , μ + ), say μ * = μ + . By (4.10) we directly see that m * = μ + w * . Moreover, rewriting and looking at (4.10) yields lim j→∞ M(z j ) = μ + w * ⊗ w * − σ * .
Next we introduce the most important -directions.

Definition 4.5.
Let z ∈ Z 0 . We callz(z) ∈ Z defined bỹ the Muskat direction associated with z. Here the definition ofq andσ is understood as decomposition into trace and traceless part. Moreover, any vector of the form z = (0,w, λw,σ ,q), λ ∈ R is called an Euler direction provided it is contained in the wave cone .
Note that the Euler direction comes from the perturbations used in [16] for the homogeneous incompressible Euler equations, while the Muskat direction is a generalization of the perturbations introduced in [33] for the Muskat problem (hence the name), having the property of conserving the quantity m−μw (μ + −μ)(μ−μ − ) , as seen in the proof of the following Lemma: Lemma 4.6. It holds that and the traceless part M(z t ) • are all independent of t. (iv) T + (z + tz) = T + (z) for all t ∈ R and all Euler directionsz withm = μ −w , as well as T − (z + tz) = T − (z) for all t ∈ R and all Euler directions of the form z = (0,w, μ +w ,σ ,q). [17]. We nonetheless present the short proof here as well. Let (w,σ , λ) ∈ R n × S n×n 0 × R,w = 0, λ = 0 and denote by P ⊥ : R n → R n the orthogonal projection ontow ⊥ := {ξ ∈ R n :w · ξ = 0}. Takeq ∈ R, such that −q is an eigenvalue of the linear map P ⊥ •σ :w ⊥ →w ⊥ , and let ξ ∈w ⊥ \{0} denote a corresponding eigenvector. Furthermore, we choose c ∈ R, such that (id −P ⊥ )σ ξ = −cλw. Then one easily checks that ⎛ ⎝σ +q id λw
Next T ± (z t ) = T ± (z) follows immediately after rewriting It remains to check that the traceless part of M(z) is invariant along the line segment in Muskat direction. Plugging into the definition of M(z) leads us to Thus for the traceless part we get (iv) obviously is true, because m+tm−μ ± (w+tw) = m−μ ± w form = μ ±w .
As a corollary, we obtain that any two points in K can be connected with adirection, up to modifying the pressure, which implies that although the wave cone is not the whole space, it is still quite big (with respect to K ).
Since z i ∈ K , we have for i = 1, 2. Therefore, we obtain that Through a simple calculation one can then show that (4.11) holds for γ = −μ − μ + . If μ 1 = μ 2 , recall that in the proof of Lemma 4.6 (i) a suitable pressureq has been chosen to be an eigenvalue of −P ⊥ •σ :w ⊥ →w ⊥ . But z 1 , z 2 ∈ K in fact implies that P ⊥ •σ vanishes on all ofw ⊥ and we can conclude the statement.
Recall the definition of π : Z → R × R n × R n × S n×n 0 in (2.6).

Lemma 4.8. The projection π(U )
is bounded in terms of e, μ ± , n and hence compact. Moreover, for every z ∈ U \K there existsz ∈ \{0}, such that z ±z ∈ U .
Proof. We first prove that π(U ) is bounded in terms of e, μ − , μ + and the dimension n. Let z ∈ U . Obviously μ ∈ (μ − , μ + ) is bounded. The inequalities (4.10) imply that there exists a constant c = c(e, μ − , μ + , n) > 0, such that (4.12) Adapting the constant when necessary we obtain which then also implies |w| c. Next observe that the matrix M(z) can be rewritten to Hence M(z) + σ is uniformly bounded by (4.12). As a consequence we obtain |tr M(z)| c. This bound on the trace together with λ max (M(z)) = Q(z) < e, due to the fact that z ∈ U , gives us a uniform bound on the whole spectrum of M(z). Therefore M(z) + σ and M(z) are both uniformly bounded. Consequently |σ | c, and π(U ) is compact. Next we show that any z ∈ U \K can be perturbed along a -segment without leaving U . Recall that U = U 0 ∪ K + ∪ K − and K ⊂ K + ∪ K − by Lemma 4.4.
The same reasoning applies also to the case z ∈ K − \K . Now let z ∈ U 0 . If Q(z) < e or if T + (z) = T − (z) we take the Muskat direction z =z(z). Because then T ± (z + tz(z)) = T ± (z) e, t ∈ (μ − − μ, μ + − μ) by Lemma 4.6 (iii). Moreover, a straightforward computation shows and thus by Lemma 4.6 (iii) we have From now on we consider the remaining case Q(z) = e and T + (z) = T − (z). Note that then necessarily λ min (M(z)) < e, because otherwise e = λ max (M(z)) = λ min (M(z)) yields M(z) • = 0 and thus Since T + (z) e, T − (z) e this equality can only hold if T + (z) = T − (z) = e, which is excluded in the case we are considering.
Let us assume T − (z) > T + (z), the other case follows similarly. We consider Euler directions of the formz = (0,w, μ +w ,σ ,q), wherew ∈ R n andσ ∈ S n×n 0 will be chosen later andq =q(w,σ ) by Lemma 4.6 (i). These Euler directions allow us to preserve T − due to Lemma 4.6 (iv), that is, Now we need to guarantee that Q(z + tz) = Q(z) = e for small enough |t| and some choice ofw,σ . As in the cases z ∈ K ± \K we can again assume that the matrix M(z) is diagonal with entries e = λ 1 λ 2 · · · λ n and λ n < e. As before we takew = e n ,m = μ + e n and the uniquely determined pair (σ , α) ∈ S n×n 0 × R satisfying For small enough |t| we therefore conclude that this Euler perturbation does not affect the maximal eigenvalue, that is, Q(z + tz) = Q(z) = e for |t| small. Furthermore, the last condition needed for z + tz ∈ U simply follows by the continuity of T + , that is, for all |t| small enough it holds that

Perturbing Along Sufficiently Long Enough Segments
In this subsection we prove that any element from U is contained in a sufficiently long admissible line segment, similarly to Section 4.3 from [17]. We recall the projection operator π defined in (2.6). We have the following result: where N = dim(Z ) and d denotes the Euclidian distance on π(Z ).
Proof. We proceed as in the proof of Lemma 4.7 from [17]. Since z ∈ U = intK co , it follows from Carathéodory's theorem that it lies in the interior of a simplex in Z spanned by K , that is there exist λ i ∈ (0, 1), We may also assume that the coefficients are ordered such that λ 1 = max i λ i , then for any j > 1 we have Indeed, one may rewrite z ± 1 2 λ j (z j − z 1 ) = i κ i z i , where κ 1 = λ 1 ∓ 1 2 λ j , κ j = λ j ± 1 2 λ j and κ i = λ i for i ∈ {1, j}, such that these coefficients are in (0, 1).
Choose j > 1 such that max i=2,...,N +1 λ i |π(z i ) − π(z 1 )| = λ j |π(z j ) − π(z 1 )|, and letz Then [z −z, z +z] ⊂ intK co and To conclude the proof of the lemma, it would suffice to havez ∈ . While this in general may not be true a priori, we know from Corollary 4.7 that it is true up to changing the pressure inz. However, since the constraints in K , and respectively the inequalities in U do not involve the pressure, this can be done such that z ±z ∈ intK co still remains valid. This concludes the proof.

Continuity of Constraints
We now go back to the (x, t) ∈ D dependent sets of constraints K (x,t) defined in (3.6). We have the following result regarding the continuity of the nonlinear constraints in (4.5), given the continuity of the associated energy. This will allow us to have a set of subsolutions which is bounded in L 2 (D). The proof of Lemma 4.10 is based on the following observation, which can be found in [12] as Lemma 3.1: Lemma 4.11. Suppose A, B ⊂ R l for some l ∈ N are compact sets and r > 0 such that • for any z ∈ B there exists z ∈ A ∩ B r (z).
The boundedness of y∈U π(K y ) follows from Lemma 4.8 and the assumption that the function e is bounded.

Proof of Theorem 2.4
In this section we conclude the proof of Theorem 2.4 by using the Baire category method.

The Baire Category Method
We introduce the notion of subsolution associated with (3.5), (3.6). The set of constraints K (x,t) is understood with respect to a from now on fixed bounded function e : × (0, T ) → R + with (x, t) → e x − 1 2 gt 2 e n , t being continuous on an open set U ⊂ D. Furthermore, for simplicity of notation, in this subsection we will, as in the proof of Lemma 4.10, denote y := (x, t). Definition 5.1. We say that z : D → Z is a subsolution of (3.5) associated with the set of constraints K y , iff it is a weak solution of (3.5) in the sense of Definition 3.1 in D, π(z) is continuous in U , z(y) ∈ K y holds for almost every y ∈ D\U and z(y) ∈ U y = intK co y for any y ∈ U .
Furthermore, among these weak solutions there exists a sequence {z k } k 1 such that π(z k ) converges weakly to π(z 0 ) in L 2 (U ; π(Z )).
The proof is similar to those in [12,33], the only main difference being that one has to carefully track the role of the projection π . However, since the existence of the pressure is implicit in Definition 3.1 due to the use of divergence-free test functions, this can be done without any serious difficulty.
The main building block of the proof is the following perturbation lemma.
To prove Lemma 5.3, we will use the following result which can be found together with its proof as Lemma 2.1 in [12].
Lemma 5.4. Let K ⊂ R n be a compact set. Then for any compact set C ⊂ intK co there exists ε > 0 such that for any compact set K ⊂ R n with d H (K , K ) < ε we have C ⊂ int(K ) co .

Proof of Theorem 5.2. Let
X 0 = z ∈ L 2 (D; π(Z )) such that z = π(z) for some subsolution z in the sense of Definition 5.1 satisfying z = z 0 on D\U } , and X denote the closure of X 0 with respect to the weak L 2 topology. From Lemma 4.10 and the boundedness of the function e it follows that X 0 is bounded, therefore X is metrizable, denote its metric by d X (·, ·). Also since the existence of the pressure is implicit in Definition 3.1 due to the use of divergence-free test functions, it follows that for any z ∈ X there exists a possibly distributional pressure q such that (z , q ) is indeed a weak solution of (3.5).
We observe that the functional I (z ) = U |z | 2 dy is a Baire-1 function on X . Indeed, setting where χ j ∈ C ∞ c (B 1/j (0)) is the standard mollifying sequence, one obtains that I j is continuous on X and that I j (z ) → I (z ) as j → +∞.
It follows from the Baire category theorem that the set Y = {z ∈ X : I is continuous at z } is residual in X . We claim that for any z ∈ Y it follows that Suppose the contrary. Then J (z ) = ε > 0 for some z ∈ Y , and let z j ∈ X 0 be a sequence which converges to z with respect to d X . Since I is continuous at z , it follows that z j → z strongly in L 2 (U ; π(Z )). Note that J is continuous with respect to the strong L 2 -topology. Therefore we may assume that J (z j ) > ε/2 for all z j .
Since z j ∈ X 0 , there exists some z j : D → Z which is a subsolution in the sense of Definition 5.1 and such that z j = π(z j ). We may then apply Lemma 5.3 to deduce that there exists δ = δ(ε) > 0 and a subsolutionz j such that U |π(z j (y) −z j (y))| 2 dy δ and π(z j −z j ) 0 weakly in L 2 . Since z j = π(z j ) → z and z ∈ Y , we conclude as before π(z j ) → z strongly in L 2 contradicting the fact that π(z j ) and z j are uniformly bounded away from each other. We thus have showed that the set of solutions J −1 (0) is residual in X .
The proof of the mixing property (5.2) follows by another application of the Baire category theorem and is exactly the same as in [6]. For convenience we briefly present it here as well. Let B be an open ball contained in U . The set for any countable union of balls B i ⊂ U . By taking all balls (B i ) i∈N ⊂ U with rational centers and radii we can conclude the statement.

Conclusion
In order to prove our convex integration result for (1.1) we apply a transformation similar to (3.4) to the differential inclusion (3.5), (3.6) and in particular also its relaxation. Recall from Section 3 that for a bounded domain ⊂ R n and T > 0 we defined D = (y, t) ∈ R n × (0, T ) : y − 1 2 gt 2 e n ∈ . Now let z = (μ, w, m, σ, q) be a weak solution of (3.5) with some suitable initial data. Defining again y := x + 1 2 gt 2 e n , as well as S(x, t) = σ (y, t) − gt (m (y, t) ⊗ e n + e n ⊗ m (y, t)) one obtains through lenghty but straightforward calculations that (ρ, v, u, S, P) is a weak solution of (2.1) with the same initial data. Also here the transformation can be inverted in an obvious way, mapping a solution of (2.1) to a solution of (3.5). Furthermore, for a given function e : × [0, T ) → R + the condition z(y, t) ∈ K (y,t) for y = x + 1 2 gt 2 , (x, t) ∈ × (0, T ) and with K (y,t) defined in (3.6) translates to (ρ, v, u, S, P)(x, t) ∈ K (x,t) with K (x,t) defined in (2.2). Similarly, if we define U (y,t) to be the interior of the convex hull of K (y,t) then by Proposition 4.2 the condition z(y, t) ∈ U (y,t) translates to (ρ, v, u, S, P)(x, t) ∈ U (x,t) with U (x,t) defined in (2.3), (2.4). Since the transformation is an affine bijection, we also see that U (x,t) is the interior of the convex hull of K (x,t) .
We have now all pieces together to prove our main result.
In other words z is a subsolution of the differential inclusion (3.5), (3.6) in the sense of Definition 5.1 (with mixing zone U ). Theorem 5.2 therefore provides us with infinitely many solutions of our differential inclusion (3.5), (3.6) which outside of U agree with z and inside U satisfy the mixing property (5.2), as well as with a sequence of solutions such that π(z k ) converges L 2 -weakly to π(z ).
One may then transfer these conclusions to the setting of Theorem 2.4 via Lemma 3.2.
Let us now briefly explain how to establish the admissibility of the obtained solutions, provided that π(z) is in addition of class C 0 ([0, T ]; L 2 ( ; π(Z ))). As before let z be the corresponding transformed subsolution defined on D. Due to an improvement of the Tartar framework as in [7,17] one can show that the induced sequence {π(z k )} k∈N not only converges weakly in L 2 (D) to π(z ), but weakly on every time-slice D(t) uniformly in t ∈ [0, T ]. It is in fact straightforward but quite lengthy to adapt the proof from [17] to our situation, therefore we omit the details, cf. also [7] and in particular Remark 2.3 therein. Transforming z k again to z k we conclude that the associated energies uniformly in t ∈ [0, T ] as k → ∞, recall the definitions (2.5), (2.7). However this does not yet allow us to conclude the admissibility of the induced solutions, since the difference goes to 0 as t 0. Nonetheless, similarly to [7, Definition 2.4] (but a lot less technical for our purposes) we can extend the definiton of the space X 0 , such that the sequence (or any solution obtained by the convex integration scheme) satisfies for all t ∈ [0, T ], k 0. The statement follows.

Subsolutions
We now turn to the construction of Rayleigh-Taylor subsolutions. We start by observing that the relaxation inside the mixing zone U ⊂ × (0, T ) given in Definition 2.3 can be equivalently rewritten (in the spirit of [6]) as the system in U . The condition on λ max (A(z)) from (2.4) with u replaced by ρv + f is kept in accordance with Definition 2.3. Indeed, in order to see this, given a subsolution z = (ρ, v, u, S, P) it suffices to set Conversely, given f , it suffices to set u := ρv + f to obtain a subsolution in the sense of Definition 2.3.
Proof of Theorem 2.7. Now let n = 2, T > 0 and ⊂ R 2 be the rectangle stated in the Theorem. In view of the equivalent reformulation above our goal is to find a suitable combination of functions ξ, η and e, such that (6.1) has a solution satisfying the energy inequality (2.8) in a strict sense.
In fact we will look for one-dimensional solutions of (6.1), that is a subsolution z in the sense of Definition 2.3, which is independent of x 1 and satisfies u = u 2 e 2 , ξ = ξ 2 e 2 , η = η 2 e 2 respectively. We further assume v ≡ 0.
If we have chosen ξ , η, then condition (6.2) implies that e in the mixing zone is determined by Note also that under condition (6.2) the denominator will always be positive for t > 0. Outside the mixing zone we will have e = 1 2 ρg 2 t 2 in accordance with (2.5). The last equation in (6.1) then becomes We recall (1.3) and hence that ρ 0 only depends on the sign of x 2 . Using the change of coordinates ρ(x, t) = y(x, gt 2 /2) and interpreting the ξ 2 , η 2 as functions of ρ only, one obtains equivalently Now if G : [ρ − , ρ + ] → R is uniformly strictly convex, then one may consider the unique entropy solution (cf. Section 3.4.4 in [19]) of (6.6) with Rayleigh-Taylor initial data ρ 0 to obtain that  Observe that this already implies that the height of the mixing zone grows (up to a constant) like 1 2 gt 2 , more precisely we will have U = (x, t) ∈ × (0, T ) : It is easy to check that if one is able to choose ξ 2 , η 2 ∈ (−1, 1) such that G is indeed uniformly strictly convex and the above entropy solution exists, then defining with u 2 and S extended by 0 outside U , one truly obtains a subsolution in the sense of Definition 2.3. Indeed the relaxed momentum equation holds by definition of P, and S is chosen in a way, such that the trace free part of A(z) vanishes. In consequence inequality (2.4) reduces to e > (ρ + + ρ − − ρ)u 2 2 2(ρ + − ρ)(ρ − ρ − ) + gtu 2 + 1 2 which holds, since by our reformulation inequalities (2.3) are automatically satisfied for ξ 2 , η 2 ∈ (−1, 1) and e defined in (6.4). Therefore, all that remains to do in order to finish the construction of RTsubsolutions is to find ξ 2 , η 2 : (ρ − , ρ + ) → (−1, 1) such that G is uniformly strictly convex and to assure the admissibility (2.8) (in a strict sense for t > 0) of the associated total energy (2.7). Denoting one has e(x 2 , t) = g 2 t 2ẽ (ρ(x 2 , t)) with e(ρ) := 1 2 By the transformation x 2 = 1 2 gt 2 G (ρ) the desired admissibility (2.8) in the strict sense is then equivalent to We further make the ansatzẽ(ρ ± ) = 1 2 ρ ± , in other words that e is continuious across ∂U . Then partial integration shows that (6.9) is equivalent to Observe that the conditionẽ(ρ ± ) = 1 2 ρ ± requires ξ 2 (ρ + ) = 1, η 2 (ρ − ) = −1. Inspired by the known families of subsolutions for the Muskat problem [33] or the Kelvin-Helmholtz instability [34], it is of interest to investigate the limit case when one is in the boundary of the convex hull, instead of its interior, as this corresponds to the limiting mixing zone growth rates of these families. In our case this means to choose |ξ | = |η| = 1 throughout all of [ρ − , ρ + ], that is ξ 2 ≡ −η 2 ≡ 1. Of course this will not lead to a strict subsolution inside the mixing zone, so we will later consider a slight perturbation in order to be into the interior of the convex hull.
Funding Open Access funding enabled and organized by Projekt DEAL. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.