1 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

$$\begin{aligned} \begin{aligned} \partial _t(\rho v)+\text {div }(\rho v\otimes v) + \nabla p&= - \rho g e_n,\\ \text {div } v&=0,\\ \partial _t \rho + \text {div }(\rho v)&=0, \end{aligned} \end{aligned}$$
(1.1)

which we consider on a bounded domain \(\Omega \subset {\mathbb {R}}^n\), \(n\geqq 2\) and a time interval [0, T). Here \(\rho :\Omega \times [0,T)\rightarrow {\mathbb {R}}\) denotes the fluid density, \(v:\Omega \times [0,T)\rightarrow {\mathbb {R}}^n\) is the velocity field, respectively \(p:\Omega \times [0,T)\rightarrow {\mathbb {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, \(\rho \equiv 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 well-posedness 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:\Omega \rightarrow {\mathbb {R}}^n\), \(\rho _0:\Omega \rightarrow {\mathbb {R}}\) satisfying

$$\begin{aligned} {{\,\mathrm{div}\,}}v_0=0\text { and } \rho _0\in \left\{ \,\rho _-,\rho _+\,\right\} \text { almost everywhere } \end{aligned}$$
(1.2)

with two fixed values \(\rho _+>\rho _->0\). In fact our main focus lies on the flat unstable initial configuration

$$\begin{aligned} v_0\equiv 0\text { and }\rho _0(x)= \left\{ \begin{array}{ll} \rho _+\text { when }x_n>0,\\ \rho _-\text { when }x_n\leqq 0, \end{array} \right. \end{aligned}$$
(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 \(\rho _+>\rho _->0\) define the Atwood number \( {{\mathcal {A}}}=\frac{\rho _+-\rho _-}{\rho _++\rho _-} \) and the quadratic functions \(c_\pm :{\mathbb {R}}\rightarrow {\mathbb {R}}\),

$$\begin{aligned} c_+(t)=\frac{\rho _++\rho _-}{2\sqrt{\rho _-}(\sqrt{\rho _+}+\sqrt{\rho _-})}{{\mathcal {A}}}gt^2,\quad c_-(t)=\frac{\rho _++\rho _-}{2\sqrt{\rho _+}(\sqrt{\rho _+}+\sqrt{\rho _-})}{{\mathcal {A}}}gt^2. \end{aligned}$$

Let \(T>0\) and \(\Omega =(0,1)\times (-c_-(T),c_+(T))\subset {\mathbb {R}}^2\).

Theorem 1.1

Let \(\frac{\rho _+}{\rho _-}\geqq \left( \frac{4+2\sqrt{10}}{3}\right) ^2\). The initial value problem (1.1), (1.3) has infinitely many weak admissible solutions \((\rho ,v)\in L^\infty (\Omega \times (0,T);{\mathbb {R}}\times {\mathbb {R}}^2)\) with \(\rho \in \{\rho _-,\rho _+\}\) almost everywhere and such that

  1. (i)

    \(\rho (x,t)=\rho _+\), \(v=0\) for \(x_2\geqq c_+(t)\),

  2. (ii)

    \(\rho (x,t)=\rho _-\), \(v=0\) for \(x_2\leqq -c_-(t)\),

  3. (iii)

    for any open ball B contained in \(\left\{ \,(x,t)\in \Omega \times (0,T):x_2\in (-c_-(t),c_+(t))\,\right\} \) it holds that

    $$\begin{aligned} \int _B\rho _+-\rho (x,t)\,\mathrm{d}x\,\mathrm{d}t\cdot \int _B\rho (x,t)-\rho _-\,\mathrm{d}x\,\mathrm{d}t>0. \end{aligned}$$

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 \({\bar{\rho }}\) in common. The profile \({\bar{\rho }}(x,t)={\bar{\rho }}(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

$$\begin{aligned} \partial _t{\bar{\rho }}+gt\partial _{x_2}G({\bar{\rho }})=0, \end{aligned}$$

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 \({\bar{\rho }}\) can be found in Section 6.

The condition that the density ratio \(\rho _+/\rho _-\) is larger than \(\left( \frac{4+2\sqrt{10}}{3}\right) ^2 \approx 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 \({{\bar{\rho }}}\) 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 \(\alpha {\mathcal {A}}g t^2\), but there is considerable disagreement about the value of the constant \(\alpha \) and its possible dependence on \({\mathcal {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_\pm \) 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 left-hand 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 \({{\mathcal {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 (xt) due to the prescribed energy density function.

In Section 4 we prove the ingredients of the topological framework, most importantly we calculate the \(\Lambda \)-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.

2 Statement of Results

Let \(\Omega \subset {\mathbb {R}}^n\) be a bounded domain and \(T>0\). Our notion of solution to equation (1.1) on \(\Omega \times [0,T)\) is as follows:

Definition 2.1

(Weak solutions) Let \((\rho _0,v_0)\in L^\infty (\Omega )\times L^2(\Omega ;{\mathbb {R}}^n)\) such that (1.2) holds almost everywhere in \(\Omega \). We say that \((\rho ,v)\in L^\infty (\Omega \times (0,T))\times L^2(\Omega \times (0,T);{\mathbb {R}}^n)\) is a weak solution to equation (1.1) with initial data \((\rho _0,v_0)\) if for any test functions \(\Phi \in C^\infty _c(\Omega \times [0,T);{\mathbb {R}}^n) \), \(\Psi _1\in C^\infty _c({\overline{\Omega }}\times [0,T)) \), \(\Psi _2\in C^\infty _c(\Omega \times [0,T)) \) such that \(\Phi \) is divergence-free, we have

$$\begin{aligned}&\int _0^T\int _{\Omega } \left[ \rho v \cdot \partial _t \Phi + \langle \rho v\otimes v ,\nabla \Phi \rangle - g \rho \Phi _n \right] \ \mathrm{d}x \ \mathrm{d}t \\&\quad +\int _{\Omega }\rho _0(x)v_0 (x)\cdot \Phi (x,0)\ \mathrm{d}x=0,\\&\int _0^T\int _\Omega v\cdot \nabla \Psi _1\,\mathrm{d}x\,\mathrm{d}t=0,\\&\int _0^T\int _{\Omega } \left[ \rho \partial _t \Psi _2 + \rho v\cdot \nabla \Psi _2 \right] \ \mathrm{d}x \ \mathrm{d}t +\int _{\Omega } \rho _0 (x) \Psi _2(x,0)\ \mathrm{d}x=0, \end{aligned}$$

and if \(\rho (x,t)\in \left\{ \,\rho _-,\rho _+\,\right\} \) for almost every \((x,t)\in \Omega \times (0,T)\).

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 \((\rho ,v)\) an energy density function \(E\in L^1(\Omega \times (0,T))\) given by

$$\begin{aligned} E(x,t):=\frac{1}{2}\rho (x,t)\left| v(x,t)\right| ^2+\rho (x,t)gx_n. \end{aligned}$$

Furthermore, for smooth solutions of (1.1) one can show that \(t\mapsto \int _{\Omega }E(x,t)\,\mathrm{d}x\) 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.

Definition 2.2

(Admissible weak solutions) A weak solution \((\rho ,v)\) in the sense of Definition 2.1 is called admissible provided it satisfies the weak energy inequality

$$\begin{aligned} \int _\Omega E(x,t)\,\mathrm{d}x\leqq \int _\Omega E(x,0)\,\mathrm{d}x\text { for almost every }t\in (0,T). \end{aligned}$$

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

$$\begin{aligned} \begin{aligned} \partial _t u+{{\,\mathrm{div}\,}}S +\nabla P&=-\rho g e_n,\\ \partial _t\rho +{{\,\mathrm{div}\,}}u&=0,\\ {{\,\mathrm{div}\,}}v&=0, \end{aligned} \end{aligned}$$
(2.1)

considered on \(\Omega \times (0,T)\) and with \(z=(\rho ,v,u,S,P)\) taking values in the space \(Z={\mathbb {R}}\times {\mathbb {R}}^n\times {\mathbb {R}}^n\times {{\mathcal {S}}}_0^{n\times n}\times {\mathbb {R}}\). Here \({{\mathcal {S}}}^{n\times n}_0\) denotes the space of symmetric \(n\times n\) matrices with trace 0. We will also write \({{\mathcal {S}}}^{n\times n}\) for the space of symmetric matrices, \({{\,\mathrm{id}\,}}\in {{\mathcal {S}}}^{n\times n}\) for the identity and \(\lambda _{\text {max}}(S),\lambda _{\text {min}}(S)\) for the maximal, minimal resp., eigenvalue of \(S\in {{\mathcal {S}}}^{n\times n}\).

As usual, equations (2.1) will be complemented by a set of pointwise constraints. Let \(e:\Omega \times (0,T)\rightarrow {\mathbb {R}}_+\) be a given function and define for \((x,t)\in \Omega \times (0,T)\) the sets

$$\begin{aligned} {{\mathcal {K}}}_{(x,t)}:= & {} \left\{ z\in Z: \rho \in \{\rho _-,\rho _+\},~u=\rho v,\right. \nonumber \\&\left. ~\rho v\otimes v-S=\left( e(x,t)-\frac{2}{n}\rho gt e_n\cdot v-\frac{1}{n}\rho g^2t^2\right) {{\,\mathrm{id}\,}}\right\} , \end{aligned}$$
(2.2)

as well as the sets \({{\mathcal {U}}}_{(x,t)}\) by requiring for \(z\in {{\mathcal {U}}}_{(x,t)}\) the following four inequalities to hold:

$$\begin{aligned}&\rho _-<\rho<\rho _+,\nonumber \\&\quad \frac{\rho _+}{n}\frac{\left| u-\rho _-v+(\rho -\rho _-)gte_n\right| ^2}{(\rho -\rho _-)^2} <e(x,t),\nonumber \\ \end{aligned}$$
(2.3)
$$\begin{aligned}&\frac{\rho _-}{n}\frac{\left| u-\rho _+v+(\rho -\rho _+) gte_n\right| ^2}{(\rho -\rho _+)^2}<e(x,t),\nonumber \\&\quad \lambda _{\text {max}}\left( A(z)\right) <e(x,t)-\frac{2}{n}gt e_n\cdot u-\frac{1}{n}\rho g^2t^2, \end{aligned}$$
(2.4)

where

$$\begin{aligned} A(z)=\frac{\rho \rho _-\rho _+v\otimes v-\rho _-\rho _+(u\otimes v+v\otimes u)+(\rho _++\rho _--\rho )u\otimes u}{(\rho _+-\rho )(\rho -\rho _-)}-S. \end{aligned}$$

Note that by the definition of \({{\mathcal {K}}}_{(x,t)}\) in (2.2) and by recalling that S has vanishing trace, every solution of (2.1) taking values in \({{\mathcal {K}}}_{(x,t)}\) almost everywhere is a solution to the inhomogeneous Euler equations (1.1) with \(\rho \in \left\{ \,\rho _-,\rho _+\,\right\} \) and associated energy

$$\begin{aligned} E=\frac{1}{2}\rho \left| v\right| ^2+\rho gx_n=\frac{n}{2}e(x,t)-\rho gt e_n\cdot v-\frac{1}{2}\rho g^2t^2+\rho gx_n, \end{aligned}$$
(2.5)

which is equivalent to saying that

$$\begin{aligned} \frac{1}{2}\rho \left| v+gte_n\right| ^2=\frac{n}{2}e(x,t). \end{aligned}$$

Conversely, if we have a solution \((\rho ,v,p)\) of (1.1) with \(\rho \in \left\{ \,\rho _-,\rho _+\,\right\} \) almost everywhere, we can introduce the variables \(u=\rho v\), \(S=\rho v\otimes v-\frac{1}{n}\rho \left| v\right| ^2{{\,\mathrm{id}\,}}\), \(P=p+\frac{1}{n}\rho \left| v\right| ^2\) to see that \(z=(\rho ,v,u,S,P)\) will satisfy system (2.1) while pointwise taking values \(z(x,t)\in {{\mathcal {K}}}_{(x,t)}\), where \({{\mathcal {K}}}_{(x,t)}\) is defined with respect to the function

$$\begin{aligned} e(x,t):=\frac{1}{n}\rho (x,t)\left| v(x,t)+gte_n\right| ^2. \end{aligned}$$

Since the pressure P does not play a role in the set of constraints \({{\mathcal {K}}}_{(x,t)}\), it is convenient to consider the following projection: for \(z=(\rho ,v,u,S,P)\in Z\) we denote

$$\begin{aligned} \pi (z)=(\rho ,v,u,S)\in {\mathbb {R}}\times {\mathbb {R}}^n\times {\mathbb {R}}^n\times \mathcal S_0^{n\times n}. \end{aligned}$$
(2.6)

Using the linear system (2.1) and the definition of \({{\mathcal {U}}}_{(x,t)}\) we define relaxed solutions to (1.1) in the following way:

Definition 2.3

(Subsolutions) Let \(e:\Omega \times [0,T)\rightarrow {\mathbb {R}}_+\) be a bounded function. We say that \(z=(\rho ,v,u,S,P):\Omega \times (0,T)\rightarrow Z\) is a subsolution of (1.1) associated with e and the initial data \((\rho _0,v_0)\in L^\infty (\Omega )\times L^2(\Omega ;{\mathbb {R}}^n)\) satisfying (1.2) iff \(\pi (z)\in L^\infty (\Omega \times (0,T);\pi (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),

$$\begin{aligned}&\int _0^T\int _{\Omega } \left[ u \cdot \partial _t \Phi + \langle S,\nabla \Phi \rangle - g \rho \Phi _n \right] \ \mathrm{d}x \ \mathrm{d}t +\int _{\Omega }\rho _0(x)v_0 (x)\cdot \Phi (x,0)\ \mathrm{d}x=0,\\&\int _0^T\int _{\Omega } \left[ \rho \partial _t \Psi + u\cdot \nabla \Psi \right] \ \mathrm{d}x \ \mathrm{d}t +\int _{\Omega } \rho _0 (x) \Psi (x,0)\ \mathrm{d}x=0, \end{aligned}$$

for any test functions \(\Phi \in C^\infty _c(\Omega \times [0,T);{\mathbb {R}}^n)\), \({{\,\mathrm{div}\,}}\Phi =0\), \(\Psi \in C^\infty _c(\Omega \times [0,T))\), and if there exists an open set \({\mathscr {U}}\subset \Omega \times (0,T)\), such that the maps \((x,t)\mapsto \pi (z(x,t))\) and \((x,t)\mapsto e(x,t)\) are continuous on \({\mathscr {U}}\) with

$$\begin{aligned}&z(x,t)\in {{\mathcal {U}}}_{(x,t)},\text { for all }(x,t)\in {\mathscr {U}},\\&z(x,t)\in {{\mathcal {K}}}_{(x,t)},\text { for almost every }(x,t)\in \Omega \times (0,T){\setminus }{\mathscr {U}}. \end{aligned}$$

We call \({\mathscr {U}}\) the mixing zone of z. Moreover, the subsolution is called admissible provided that

$$\begin{aligned} E_{sub}(x,t):=\frac{n}{2} e(x,t)- gt e_n\cdot u(x,t)-\frac{1}{2}\rho (x,t)g^2t^2+\rho (x,t) gx_n \end{aligned}$$
(2.7)

satisfies

$$\begin{aligned} \int _\Omega E_{sub}(x,t)\, \mathrm{d}x\leqq \int _\Omega E_{sub}(x,0)\,\mathrm{d}x \text { for almost every }t\in (0,T). \end{aligned}$$
(2.8)

We now can state the following criterion for the existence of infinitely many weak solutions:

Theorem 2.4

Let \(n=2\) and \(e:\Omega \times [0,T)\rightarrow {\mathbb {R}}_+\) be bounded. If there exists a subsolution z associated with e in the sense of Definition 2.3, then for the same initial data of the subsolution there exist infinitely many weak solutions in the sense of Definition 2.1, which coincide almost everywhere on \(\Omega \times (0,T){\setminus }{\mathscr {U}}\) with z and whose total energy is given by E defined in (2.5). The solutions are turbulently mixing on \({\mathscr {U}}\) in the sense that for any open ball \(B\subset {\mathscr {U}}\) it holds that

$$\begin{aligned} \int _B\rho _+-\rho (x,t)\,d(x,t)\cdot \int _B\rho (x,t)-\rho _-\,d(x,t)>0. \end{aligned}$$
(2.9)

Among these weak solutions there exists a sequence \(\{z_k\}_{k\geqq 0}\) such that \(\rho _k\rightharpoonup \rho \) in \(L^2({\mathscr {U}})\). If in addition \(\pi (z)\) is in \({{\mathcal {C}}}^0([0,T];L^2(\Omega ;\pi (Z)))\) and satisfies (2.8) with strict inequality for every \(t\in (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 \({\mathscr {U}}\) the mixing zone and to interpret the subsolution density \(\rho \) 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 \(\Lambda \)-convex hull in Section 4.2, are carried out in arbitrary dimensions.

(c) We will see later that the open set \({{\mathcal {U}}}_{(x,t)}\) is indeed the convex hull of \({{\mathcal {K}}}_{(x,t)}\). In particular we can conclude that weak limits of solutions are subsolutions in the following sense: Let \((\rho _k,v_k)_{k\in {\mathbb {N}}}\) be a sequence of essentially bounded weak solutions of (1.1) and define as before \(u_k:=\rho _kv_k\), \(S_k:=\rho _kv_k\otimes v_k-\frac{1}{n}\rho _k\left| v_k\right| ^2{{\,\mathrm{id}\,}}\). Assume that \(z_k':=(\rho _k,v_k,u_k,S_k)\overset{*}{\rightharpoonup } (\rho ,v,u,S)=:z'\) in \(L^\infty (\Omega \times (0,T);{\mathbb {R}}\times {\mathbb {R}}^n\times {\mathbb {R}}^n\times {{\mathcal {S}}}_0^{n\times n})\). Assume further that there exists a continuous bounded function \(e\in {{\mathcal {C}}}^0(\Omega \times (0,T))\), such that \(e_k:=\frac{1}{n}\rho _k\left| v_k+gte_n\right| ^2\rightarrow e\) in \(L^\infty (\Omega \times (0,T))\). Then \(z'\) supplemented by a possibly distributional P is a weak solution of the linear system 2.1 with \((z'(x,t),P(x,t))\in {\overline{{{\mathcal {U}}}}}_{(x,t)}\) for almost every \((x,t)\in \Omega \times (0,T)\), where \({{\mathcal {U}}}_{(x,t)}\) is defined with respect to the function e.

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 \(\Omega \cap ({\mathbb {R}}^{n-1}\times \{0\})\ne \emptyset \).

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\in (0,T)\).

Theorem 2.7

Let \(n=2\), \(\Omega =(0,1)\times (-c_-(T),c_+(T))\), where

$$\begin{aligned} c_-(t)=\frac{1}{2}\left( 1-\sqrt{\frac{\rho _-}{\rho _+}}\right) gt^2,\quad c_+(t)=\frac{1}{2}\left( \sqrt{\frac{\rho _+}{\rho _-}}-1\right) gt^2. \end{aligned}$$

If \(\rho _+>\left( \frac{4+2\sqrt{10}}{3}\right) ^2\rho _-\), then there exists a RT-subsolution z which only depends on \(\frac{x_2}{t^2}\), and at time \(t>0\) the mixing zone \({\mathscr {U}}(t):=\left\{ \,x\in \Omega :(x,t)\in {\mathscr {U}}\,\right\} \) associated with z is \((0,1)\times (-c_-(t),c_+(t))\).

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.

3 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 \(\Omega \subset {\mathbb {R}}^n\) be a bounded domain, \(T>0\) and set

$$\begin{aligned} {\mathscr {D}}=\left\{ \,(y,t)\in {\mathbb {R}}^n\times (0,T):y-\frac{1}{2}gt^2e_n\in \Omega \,\right\} , \end{aligned}$$

such that for \(t\in (0,T)\) the slice is given by

$$\begin{aligned} {\mathscr {D}}(t):=\left\{ \,y\in {\mathbb {R}}^n:(y,t)\in {\mathscr {D}}\,\right\} =\Omega +\frac{1}{2}gt^2e_n. \end{aligned}$$

Let \((\mu ,w,q)\) be a weak solution of

$$\begin{aligned} \begin{aligned} \partial _t(\mu w)+\text {div }(\mu w\otimes w) + \nabla q&= 0,\\ \text {div } w&=0,\\ \partial _t \mu + \text {div }(\mu w)&=0 \end{aligned} \end{aligned}$$
(3.1)

on \({\mathscr {D}}\) for some suitable initial data satisfying (1.2) and with weak boundary condition

$$\begin{aligned} (w(y,t)-gte_n)\cdot \nu _{{\mathscr {D}}(t)}(y)=0 \end{aligned}$$
(3.2)

for \(y\in \partial {\mathscr {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 \(\Omega \times (0,T)\), \(\Omega \times [0,T)\) is replaced by \({\mathscr {D}}\), \({\mathscr {D}}\cup (\Omega \times \{0\})\) resp., and the weak formulation of \({{\,\mathrm{div}\,}}w=0\) including the weak boundary condition (3.2) becomes

$$\begin{aligned}&\int _{{\mathscr {D}}}w\cdot \nabla \Psi \,d(y,t)-\int _0^T \int _{\partial {\mathscr {D}}(t)}\Psi (y,t)gte_n\cdot \nu _{{\mathscr {D}} (t)}(y)\,\mathrm{d}S(y)\,\mathrm{d}t \nonumber \\&\quad =0\text { for all }\Psi \in {{\mathcal {C}}}^\infty (\overline{{\mathscr {D}}}). \end{aligned}$$
(3.3)

Then if we define \(y:=x+\frac{1}{2}gt^2e_n\) and set

$$\begin{aligned} \begin{aligned} \rho (x,t)&=\mu \left( y,t\right) ,\\ v(x,t)&=w\left( y,t\right) -g t e_n,\\ p(x,t)&=q\left( y,t\right) , \end{aligned} \end{aligned}$$
(3.4)

it is straightforward to check that \((\rho ,v)\) is a weak solution of (1.1) on \(\Omega \times (0,T)\) with the same initial data \((\rho _0,v_0)=(\mu _0,w_0)\). 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 \(\frac{1}{2}\mu (y,t)\left| w(y,t)\right| ^2\). Let us write

$$\begin{aligned} \frac{1}{2}\mu (y,t)\left| w(y,t)\right| ^2=\frac{n}{2}e(y-1/2gt^2e_n,t) \end{aligned}$$

for a function \(e:\Omega \times [0,T)\rightarrow {\mathbb {R}}_+\). Then the total energy E(xt) 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 \({\mathscr {D}}\) the system

$$\begin{aligned} \begin{aligned} \partial _t m+{{\,\mathrm{div}\,}}\sigma +\nabla q=0,\\ {{\,\mathrm{div}\,}}w=0,\\ \partial _t\mu +{{\,\mathrm{div}\,}}m=0, \end{aligned} \end{aligned}$$
(3.5)

where \(z:=(\mu ,w,m,\sigma ,q)\) takes values in \(Z={\mathbb {R}}\times {\mathbb {R}}^n\times {\mathbb {R}}^n\times {\mathcal {S}}_0^{n\times n}\times {\mathbb {R}}\), together with the set of pointwise constraints

$$\begin{aligned} K_{(y,t)}=\left\{ z\in Z: \mu \in \{\mu _-,\mu _+\},~m=\mu w,\ \mu w\otimes w-\sigma =e\left( y-\frac{1}{2}gt^2e_n,t\right) {{\,\mathrm{id}\,}}\right\} , \end{aligned}$$
(3.6)

where in analogy to the homogeneous Euler equations \(e:\Omega \times (0,T)\rightarrow {\mathbb {R}}_+\) is given and for the sake of consistency we have denoted \(\mu _\pm :=\rho _\pm \). We will understand weak solutions of (3.5) in the following sense:

Definition 3.1

We say that \(z:{\mathscr {D}}\rightarrow Z\) is a weak solution of (3.5) with initial data \(\pi (z_0)\in L^2(\Omega ;\pi (Z))\) iff \(\pi (z)\in L^2({\mathscr {D}};\pi (Z))\), q is a distribution, w satisfies (3.3) and one has

$$\begin{aligned}&\int _{{\mathscr {D}}} \left[ m \cdot \partial _t \Phi + \langle \sigma ,\nabla \Phi \rangle \right] \ \mathrm{d}x \ \mathrm{d}t +\int _{\Omega }\mu _0(x)w_0 (x)\cdot \Phi (x,0)\ \mathrm{d}x=0,\\&\int _{{\mathscr {D}}} \left[ \mu \partial _t \Psi + m\cdot \nabla \Psi \right] \ \mathrm{d}x \ \mathrm{d}t +\int _{\Omega } \mu _0 (x) \Psi (x,0)\ \mathrm{d}x=0, \end{aligned}$$

for any \(\Phi \in C^\infty _c({\mathscr {D}}\cup (\Omega \times \{0\});{\mathbb {R}}^n)\), \({{\,\mathrm{div}\,}}\Phi =0\), \(\Psi \in C^\infty _c({\mathscr {D}}\cup (\Omega \times \{0\}))\).

This way we have arrived at a reformulation of equation (1.1) as a differential inclusion. The process is summarized in the following statement.

Lemma 3.2

Let \((\rho _0,v_0)\in L^\infty (\Omega )\times L^2(\Omega ;{\mathbb {R}}^n)\) be initial data satisfying (1.2), \(e\in L^1(\Omega \times (0,T);{\mathbb {R}}_+)\) be a prescribed function. If \(z=(\mu ,w,m,\sigma ,q)\) is a weak solution of (3.5) in the sense of Definition 3.1 with initial data \(\mu (\cdot ,0)=\rho _0\), \(w(\cdot ,0)=v_0\) and if \(z(y,t)\in K_{(y,t)}\) for almost every \((y,t)\in {\mathscr {D}}\), then the pair \((\rho ,v)\) defined by (3.4) is a weak solution of (1.1) on \(\Omega \times (0,T)\) with initial data \((\rho _0,v_0)\). Moreover, the (possibly distributional) pressure is given by

$$\begin{aligned} p(x,t):=q(y,t)-\frac{1}{n}\mu (y,t)\left| w(y,t)\right| ^2,\quad y=x+\frac{1}{2}gt^2e_n, \end{aligned}$$

and the associated energy E by (2.5).

4 The Ingredients of the Tartar Framework

The general strategy of the Tartar framework relies on the following steps:

  • finding a wave cone \(\Lambda \subset Z\) such that for any \({\bar{z}}\in \Lambda \), one can construct a localized plane wave associated with (3.5) oscillating in the direction of \({\bar{z}}\);

  • calculating the \(\Lambda \)-convex hull of \(K_{(x,t)}\) (denoted by \(K_{(x,t)}^\Lambda \)) and proving that one can perturb any element in its interior along sufficiently long \(\Lambda \)-segments, provided that one is far enough from \(K_{(x,t)}\);

  • deducing an appropriate set of subsolutions using \(K_{(x,t)}^\Lambda \) and proving that it is a bounded, nonempty subset of \(L^2({\mathscr {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]).

4.1 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):

$$\begin{aligned} \Lambda =\left\{ \,{\bar{z}}\in Z:\ker \begin{pmatrix} {\bar{\sigma }} +{\bar{q}}{{\,\mathrm{id}\,}}&{}\quad {\bar{m}}\\ {\bar{m}}^T &{}\quad {\bar{\mu }}\\ {\bar{w}}^T &{}\quad 0 \end{pmatrix} \ne \{0\},\quad ({\bar{\mu }},{\bar{m}})\ne 0\,\right\} . \end{aligned}$$

It has the property that for \({\bar{z}}\in \Lambda \) there exists \(\eta \in {\mathbb {R}}^{n+1}{\setminus }\{0\}\) such that every \(z(x,t)={\bar{z}}h((x,t)\cdot \eta )\), \(h\in {{\mathcal {C}}}^1({\mathbb {R}})\) is a solution of (3.5). In Lemma 4.1 below we localize these solutions by constructing suitable potentials. Note that the condition \(({\bar{\mu }},{\bar{m}})\ne 0\) serves to eliminate the degenerate case when the first n components of \(\eta \) vanish, that is when one is only allowed to oscillate in time.

Recall the projection operator \(\pi \) defined in (2.6).

Lemma 4.1

There exists \(C>0\) such that for any \({\bar{z}}\in \Lambda \), there exists a sequence

$$\begin{aligned}z_N\in C_c^\infty (B_1(0);Z)\end{aligned}$$

solving (3.5) and satisfying that

  1. (i)

    \(d(z_N,[-{\bar{z}},{\bar{z}}])\rightarrow 0\) uniformly,

  2. (ii)

    \(z_N\rightharpoonup 0\) in \(L^2\),

  3. (iii)

    \(\int \int |\pi (z_N)|^2\, \mathrm{d}x \, \mathrm{d}t\geqq C|\pi ({\bar{z}})|^2.\)

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 \(\psi :{\mathbb {R}}^{2+1}\rightarrow {\mathbb {R}}\), \(\phi :{\mathbb {R}}^{2+1}\rightarrow {\mathcal {S}}^{2\times 2}\), setting \(D(\phi ,\psi )=(\mu ,w,m,\sigma ,q)\) with

$$\begin{aligned} \mu ={{\,\mathrm{div}\,}}{{\,\mathrm{div}\,}}\phi ,\quad w=\nabla ^{\perp } \psi ,\quad m=-\partial _t {{\,\mathrm{div}\,}}\phi ,\quad q=\frac{1}{2}{{\,\mathrm{tr}\,}}\partial _{tt} \phi ,\quad \sigma =\partial _{tt} \phi -q{{\,\mathrm{id}\,}}, \end{aligned}$$

implies that \(D(\phi ,\psi )\) solves (3.5).

Let \(S:{\mathbb {R}}\rightarrow {\mathbb {R}}\) be a smooth function, \(N\geqq 1\) and \({\bar{z}}\in \Lambda \) with \(({\bar{\mu }},{\bar{m}})\ne 0\). It follows from the definition of \(\Lambda \) that there exists

$$\begin{aligned} 0\ne (\xi ,c)\in \ker \begin{pmatrix} {\bar{\sigma }} +{\bar{q}}{{\,\mathrm{id}\,}}&{}\quad {\bar{m}}\\ {\bar{m}}^T &{} \quad {\bar{\mu }}\\ {\bar{w}}^T &{}\quad 0 \end{pmatrix}. \end{aligned}$$
(4.1)

We then treat two cases.

Case 1: \(c\ne 0\)

Note that in this case we also have \(\xi \ne 0\), since \(\xi =0\) would imply \(({\bar{\mu }},{\bar{m}})=0\).

We then set

$$\begin{aligned} \phi _N(x,t)&=\frac{1}{c^2}({\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}}) \frac{1}{N^2}S(N(\xi ,c)\cdot (x,t)),\\ \psi _N(x,t)&=|{\bar{w}}|\frac{\text {sgn}(\xi ^\perp \cdot {\bar{w}})}{|\xi |}\frac{1}{N}S'(N(\xi ,c)\cdot (x,t)), \end{aligned}$$

and we claim that

$$\begin{aligned} D(\phi _N,\psi _N)={\bar{z}}S''(N(\xi ,c)\cdot (x,t)). \end{aligned}$$
(4.2)

Indeed, using (4.1), one has

$$\begin{aligned} {{\,\mathrm{div}\,}}{{\,\mathrm{div}\,}}\phi _N&=\frac{1}{c^2}\xi ^T({\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}}) \xi S''(N(\xi ,c)\cdot (x,t))\\ {}&=\frac{1}{c^2}\xi ^T(-c{\bar{m}}) S''(N(\xi ,c)\cdot (x,t))={\bar{\mu }}S''(N(\xi ,c)\cdot (x,t)),\\ \partial _t{{\,\mathrm{div}\,}}\phi _N&=\frac{1}{c}({\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}}) \xi S''(N(\xi ,c)\cdot (x,t))=-{\bar{m}} S''(N(\xi ,c)\cdot (x,t)),\\ \partial _{tt}\phi _N&=c^2\frac{1}{c^2}({\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}}) S''(N(\xi ,c)\cdot (x,t))=({\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}}) S''(N(\xi ,c)\cdot (x,t)),\\ \nabla ^{\perp }\psi _N&=\xi ^{\perp }|{\bar{w}}|\frac{\text {sgn} (\xi ^\perp \cdot {\bar{w}})}{|\xi |}S''(N(\xi ,c)\cdot (x,t)) ={\bar{w}} S''(N(\xi ,c)\cdot (x,t)). \end{aligned}$$

From here on, the localization is done in the standard fashion (for example as in [11, 16]). We fix \(S(\cdot )=-\cos (\cdot )\) and, for \(\varepsilon >0\), consider \(\chi _\varepsilon \in C_c^\infty (B_1(0))\) satisfying \(|\chi _\varepsilon |\leqq 1\) on \(B_{1}(0)\), \(\chi _\varepsilon =1\) on \(B_{1-\varepsilon }(0)\). It is then straightforward to check that \(z_N=D(\chi _\varepsilon (\phi _N,\psi _N))\) satisfies the conclusions of the lemma.

Case 2: \(c=0\)

In this case we are not allowed to oscillate in time. However, we have \(\xi \ne 0\), so we may also assume without loss of generality that \(|\xi |=1\). On the other hand, (4.1) implies that there exist constants \(k_1,k_2,k_3\in {\mathbb {R}}\) such that

$$\begin{aligned} {\bar{w}}=k_1\xi ^\perp ,\quad {\bar{m}}=k_2\xi ^\perp ,\quad {\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}}=k_3\xi ^\perp \otimes \xi ^\perp . \end{aligned}$$
(4.3)

We set

$$\begin{aligned} \phi _{N}(x,t)={\bar{\mu }}{{\,\mathrm{id}\,}}\frac{1}{N^2}S(N\xi \cdot x),\quad \psi _{N}(x,t)=|{\bar{w}}|\frac{\text {sgn}(\xi ^\perp \cdot {\bar{w}})}{|\xi |}\frac{1}{N}S'(N\xi \cdot x), \end{aligned}$$

from where with similar calculations as in Case 1, we obtain that

$$\begin{aligned} D(\phi _{N},\psi _{N})=({\bar{\mu }},{\bar{w}},0,0,0)S''(N\xi \cdot x). \end{aligned}$$
(4.4)

To handle the remaining terms \(({\bar{m}},{\bar{\sigma }},{\bar{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 \(\omega :{\mathbb {R}}^{2+1}\rightarrow {\mathbb {R}}^{2+1}\), defining \(W={{\,\mathrm{curl}\,}}_{(x,t)}\omega \) and \({\tilde{D}}(\omega )=(0,0,m,\sigma ,q)\) by

$$\begin{aligned} m=-\frac{1}{2}\nabla ^\perp W_3,\quad \sigma +q{{\,\mathrm{id}\,}}=\begin{pmatrix} \partial _2 W_1 &{}\quad \frac{1}{2}(\partial _2W_2-\partial _1W_1)\\ \frac{1}{2}(\partial _2W_2-\partial _1W_1) &{}\quad -\partial _1W_2 \end{pmatrix} \end{aligned}$$

implies that \({\tilde{D}}(\omega )\) solves (3.5).

Now, if we consider \(\omega \) of the form

$$\begin{aligned} \omega _N(x)=(a,b,a)\frac{1}{N^2} S(N \xi \cdot x), \end{aligned}$$

for some constants \(a,b\in {\mathbb {R}}\), with S as before, we obtain that

$$\begin{aligned}&\begin{pmatrix} \partial _2 W_1 &{}\quad \frac{1}{2}(\partial _2W_2-\partial _1W_1)\\ \frac{1}{2}(\partial _2W_2-\partial _1W_1) &{}\quad -\partial _1W_2 \end{pmatrix}=a\xi ^\perp \otimes \xi ^\perp S''(N \xi \cdot x),\\&\quad \nabla ^\perp W_3=(\xi _1b-\xi _2a)\xi ^\perp S''(N \xi \cdot x). \end{aligned}$$

If \(\xi _1\ne 0\), it follows from (4.3) that setting \(a=k_3\), \(b=\frac{-2k_2+k_3\xi _2}{\xi _1}\) gives us

$$\begin{aligned} {\tilde{D}}(\omega _{N})=(0,0,{\bar{m}},{\bar{\sigma }},{\bar{q}})S''(N\xi \cdot x). \end{aligned}$$

from where, using (4.4), we get

$$\begin{aligned} D(\phi _{N},\psi _{N})+{\tilde{D}}(\omega _{N})={\bar{z}} S''(N\xi \cdot x). \end{aligned}$$

The localization is then done as in Case 1, by considering \(z_N=D(\chi _\varepsilon (\phi _{N},\psi _{N}))+{\tilde{D}}(\chi _\varepsilon \omega _{N}).\)

If \(\xi _1=0\), then choosing \(a=k_3\) gives us that

$$\begin{aligned} {\tilde{D}}(\omega _{N})=\left( 0,0,\frac{k_3}{2}\xi _2\xi ^\perp , {\bar{\sigma }},{\bar{q}}\right) S''(N\xi \cdot x). \end{aligned}$$

However, it is easy to see that for any smooth function \(\theta :{\mathbb {R}}^{2+1}\rightarrow {\mathbb {R}}\), \({\hat{D}}(\theta )=(0,0,\nabla ^\perp \theta ,0,0)\) also solves (3.5). Therefore, we may consider the potential given by

$$\begin{aligned} \theta _N(x)=\left( k_2-\xi _2\frac{k_3}{2}\right) \frac{1}{N}S'(N\xi \cdot x), \end{aligned}$$

we obtain that

$$\begin{aligned} \nabla ^\perp \theta _N(x)=\left( k_2-\xi _2\frac{k_3}{2}\right) \xi ^\perp S''(N\xi \cdot x), \end{aligned}$$

and using (4.3), we get that

$$\begin{aligned} D(\phi _{N},\psi _{N})+{\tilde{D}}(\omega _{N})+{\hat{D}}(\theta _N)={\bar{z}} S''(N\xi \cdot x). \end{aligned}$$

One may then localize this potential by the usual means in order to conclude the proof of the lemma. \(\square \)

4.2 The \(\Lambda \)-Convex Hull

We now turn to the set of pointwise constraints \(K_{(x,t)}\), \((x,t)\in {\mathscr {D}}\) defined in (3.6). The \(\Lambda \)-convex hull \(K_{(x,t)}^\Lambda \) is defined by saying that \(z\in K_{(x,t)}^\Lambda \) iff for all \(\Lambda \)-convex functions \(f:Z\rightarrow {\mathbb {R}}\) there holds \(f(z)\leqq \sup _{z'\in K_{(x,t)}}f(z')\), see [23] for more details. In our case it turns out that the \(\Lambda \)-convex hull is nothing else but the usual convex hull, see Proposition 4.2 below.

For the computation of the hull we drop the (xt) dependence of the sets \(K_{(x,t)}\) and consider a general set of pointwise constraints given by

$$\begin{aligned} K=\left\{ \,z=(\mu ,w,m,\sigma ,q)\in Z: \mu \in \{\mu _-,\mu _+\},~m=\mu w,~\mu w\otimes w-\sigma =e{{\,\mathrm{id}\,}}\,\right\} , \end{aligned}$$
(4.5)

where \(0<\mu _-<\mu _+\), \(e\in {\mathbb {R}}_+\) are given constants.

Define \(Z_0:=\left\{ \,z\in Z:\mu \in (\mu _-,\mu _+)\,\right\} \) and \(T_+,T_-,Q:Z_0\rightarrow {\mathbb {R}}\), \(M:Z_0\rightarrow {{\mathcal {S}}}^{n\times n}\),

$$\begin{aligned} M(z)&=\frac{\mu \mu _-\mu _+w\otimes w-\mu _-\mu _+(m\otimes w+w\otimes m)+(\mu _++\mu _--\mu )m\otimes m}{(\mu _+-\mu )(\mu -\mu _-)}-\sigma ,\\ Q(z)&=\lambda _{\text {max}}(M(z)),\quad T_{\pm }(z)=\frac{\mu _\pm }{n}\frac{\left| m-\mu _\mp w\right| ^2}{(\mu -\mu _\mp )^2}, \end{aligned}$$

as well as the open set

$$\begin{aligned} U=\left\{ \,z\in Z:\mu \in (\mu _-,\mu _+),~T_+(z)<e,~T_-(z)<e,~Q(z)<e\,\right\} . \end{aligned}$$
(4.6)

Proposition 4.2

The \(\Lambda \)-convex hull of K coincides with the convex hull of K and is given by \({\overline{U}}\), that is, \(K^\Lambda =K^{co}={\overline{U}}\).

Lemma 4.4 below shows that the closure of U can be written as

$$\begin{aligned} {\overline{U}}=K_-'\cup {\overline{U}}_0\cup K_+', \end{aligned}$$

where

$$\begin{aligned} {\overline{U}}_0&=\left\{ \,z\in Z:\mu \in (\mu _-,\mu _+),~T_+(z)\leqq e,~T_-(z)\leqq e,~Q(z)\leqq e\,\right\} ,\\ K_\pm '&=\{z\in Z:\mu =\mu _\pm ,~m=\mu _\pm w,~\lambda _{\text {max}}(\mu _\pm w\otimes w-\sigma )\leqq e\}. \end{aligned}$$

Moreover, Lemma 4.8 actually shows that \(K_+'\), \(K_-'\) resp., is nothing but the \(\Lambda \)-convex hull of \(K_+:=K\cap \{\mu =\mu _+\}\), \(K_-:=K\cap \{\mu =\mu _-\}\) resp..

Furthermore, notice that if one lets \(\mu _+-\mu _-\rightarrow 0\), one recovers from \({\overline{U}}\) exactly the convex hull of the constraints for the homogeneous Euler equations, cf. [17].

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

$$\begin{aligned} Q(z)=\sup _{\xi \in S^{n-1}}\xi ^TM(z)\xi =\sup _{\xi \in S^{n-1}}\left( g_\xi (z)-\xi ^T\sigma \xi \right) , \end{aligned}$$

where for every fixed \(\xi \in S^{n-1}\) the function \(g_\xi :Z_0\rightarrow {\mathbb {R}}\) is given by

$$\begin{aligned} g_\xi (z)&=\xi ^TM(z)\xi +\xi ^T\sigma \xi \\&=\frac{\mu \mu _-\mu _+(w\cdot \xi )^2-2\mu _-\mu _+(m\cdot \xi ) (w\cdot \xi )+(\mu _++\mu _--\mu )(m\cdot \xi )^2}{(\mu _+-\mu )(\mu -\mu _-)}. \end{aligned}$$

We will show that every \(g_\xi \) is convex. As a consequence Q is convex as a supremum of convex functions. In order to do this let us complement \(\xi \in S^{n-1}\) to a orthonormal basis \((\xi ,v_2,\ldots ,v_n)\) of \({\mathbb {R}}^n\). Expressing w and m with respect to this basis one sees that it is enough to show that the function \(g:(\mu _-,\mu _+)\times {\mathbb {R}}^2\rightarrow {\mathbb {R}}\),

$$\begin{aligned} g(\mu ,x)=\frac{\mu \mu _-\mu _+x_1^2-2\mu _-\mu _+x_1x_2+(\mu _++\mu _--\mu ) x_2^2}{(\mu _+-\mu )(\mu -\mu _-)} \end{aligned}$$

is convex. We write \(g(\mu ,x)=x^TA(\mu )x\) with

$$\begin{aligned} A(\mu ):=\frac{1}{(\mu _+-\mu )(\mu -\mu _-)}\begin{pmatrix} \mu \mu _-\mu _+ &{}\quad -\mu _-\mu _+ \\ -\mu _-\mu _+ &{}\quad \mu _++\mu _--\mu \end{pmatrix}. \end{aligned}$$

Let us fix \((\mu ,x)\in (\mu _-,\mu _+)\times {\mathbb {R}}^2\) and observe that \(A(\mu )\) is positive definite because \(\mu \mu _-\mu _+>0\) and

$$\begin{aligned} \det [(\mu _+-\mu )(\mu -\mu _-)A(\mu )]=\mu _-\mu _+(\mu _+-\mu )(\mu -\mu _-)>0. \end{aligned}$$

Thus the restricted function \(g(\mu ,\cdot )\) is convex, or equivalently \(D^2g(\mu ,x)[0,y]^2\geqq 0\) for all \(y\in {\mathbb {R}}^2\). It therefore remains to show that \(D^2g(\mu ,x)[1,y]^2\geqq 0\) for all \(y\in {\mathbb {R}}^2\). By the positive definiteness of \(A(\mu )\) we obtain

$$\begin{aligned} D^2g(\mu ,x)[1,y]^2&=x^TA''(\mu )x+4y^TA'(\mu )x +2y^TA(\mu )y\\&=2\left( y+A(\mu )^{-1}A'(\mu )x\right) ^TA(\mu )\left( y+A(\mu )^{-1}A'(\mu )x\right) \\&\quad +x^TA''(\mu )x -2x^TA'(\mu )A(\mu )^{-1}A'(\mu )x\\&\geqq x^T\left( A''(\mu ) -2A'(\mu )A(\mu )^{-1}A'(\mu )\right) x. \end{aligned}$$

Now we claim that in fact \(A''(\mu )=2A'(\mu )A(\mu )^{-1}A'(\mu )\), which finishes the proof. Indeed, differentiation of the identity

$$\begin{aligned}&(\mu _+-\mu )(\mu -\mu _-)A(\mu )=\begin{pmatrix} \mu \mu _-\mu _+&{} -\mu _-\mu _+ \\ -\mu _-\mu _+ &{} \mu _++\mu _--\mu \end{pmatrix} \end{aligned}$$

shows that

$$\begin{aligned} (\mu _+-\mu )(\mu -\mu _-)A'(\mu )&=(2\mu -\mu _--\mu _+)A(\mu )+C, \end{aligned}$$
(4.7)
$$\begin{aligned} (\mu _+-\mu )^2(\mu -\mu _-)^2A''(\mu )&=2((\mu _+-\mu )(\mu -\mu _-) +(2\mu -\mu _--\mu _+)^2)A(\mu ) \nonumber \\&\quad +2(2\mu -\mu _--\mu _+)C, \end{aligned}$$
(4.8)

where

$$\begin{aligned} C:=\begin{pmatrix} \mu _-\mu _+ &{}\quad 0\\ 0 &{}\quad -1 \end{pmatrix}. \end{aligned}$$

Moreover, in a straightforward way one can check that

$$\begin{aligned} (\mu _+-\mu )(\mu -\mu _-)\left( A(\mu )C^{-1}\right) ^2+(\mu _-+\mu _+-2\mu ) A(\mu )C^{-1}={{\,\mathrm{id}\,}}_{{\mathbb {R}}^2}, \end{aligned}$$

which implies that

$$\begin{aligned} CA(\mu )^{-1}C=(\mu _+-\mu )(\mu -\mu _-)A(\mu )+(\mu _-+\mu _+-2\mu )C. \end{aligned}$$
(4.9)

Now (4.7)–(4.9) imply the identity \(A''(\mu )=2A'(\mu )A(\mu )^{-1}A'(\mu )\). \(\square \)

Lemma 4.4

The set U is convex and \({\overline{U}}=K_-'\cup {\overline{U}}_0\cup K_+'\). In particular \(K\subset {\overline{U}}\).

Proof

For \(\mu \in (\mu _-,\mu _+)\) the two conditions \(T_+(z)<e\), \(T_-(z)<e\) can be rewritten as

$$\begin{aligned} \begin{aligned} \left| m-\mu _- w\right|&<c_+ (\mu -\mu _-),\\ \left| m-\mu _+ w\right|&<c_- (\mu _+-\mu ), \end{aligned} \end{aligned}$$
(4.10)

where \( c_\pm =\left( \frac{n e}{\mu _\pm }\right) ^{1/2}. \) 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 \({\overline{U}}\). Clearly \({\overline{U}}_0\subset {\overline{U}}\). Let us show that \(K_+'\subset {\overline{U}}\). The inclusion \(K_-'\subset {\overline{U}}\) can be obtained in the same way. Let \(z_*\in K_+'\). Take any \(z'\in K\) with \(\mu '=\mu _-\) and some sequence \((\mu _j)_{j\in {\mathbb {N}}}\subset (\mu _-,\mu _+)\) converging to \(\mu _+\). Define

$$\begin{aligned} z_j=\frac{\mu _+-\mu _j}{\mu _+-\mu _-}z'+\frac{\mu _j-\mu _-}{\mu _+-\mu _-}z_*. \end{aligned}$$

Clearly \(z_j\rightarrow z_*\) as \(j\rightarrow \infty \). Since \(z_*\in K_+'\) and \(z'\in K_-\) a short calculation shows

$$\begin{aligned} T_+(z_j)=\frac{\mu _+}{n}\left| w_*\right| ^2 =\frac{1}{n}{{\,\mathrm{tr}\,}}(\mu _+w_*\otimes w_*-\sigma _*) \leqq \lambda _{\text {max}}(\mu _+w_*\otimes w_*-\sigma _*) \leqq e. \end{aligned}$$

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'\otimes w'\), \(w'\otimes w_*\), \(w_*\otimes w'\), \(w_*\otimes w_*\) and find

$$\begin{aligned} M(z_j)&=\frac{\mu _+-\mu _j}{\mu _+-\mu _-}\big (\mu _-w'\otimes w'-\sigma '\big ) +\frac{\mu _j-\mu _-}{\mu _+-\mu _-}\big (\mu _+w_*\otimes w_*-\sigma _*\big )\\&=\frac{\mu _+-\mu _j}{\mu _+-\mu _-}e{{\,\mathrm{id}\,}}+\frac{\mu _j-\mu _-}{\mu _+-\mu _-} \big (\mu _+w_*\otimes w_*-\sigma _*\big ). \end{aligned}$$

We conclude \(Q(z_j)=\lambda _{\text {max}}(M(z_j))\leqq e\). Hence every \(z_j\) and therefore also the limit \(z_*\) is contained in \({\overline{U}}\). So far we know \(K_-'\cup {\overline{U}}_0\cup K_+'\subset {\overline{U}}\).

For the other inclusion consider \((z_j)_{j\in {\mathbb {N}}}\subset U\), \(z_j\rightarrow z_*\). The interesting case of course is \(\mu _*\notin (\mu _-,\mu _+)\), say \(\mu _*=\mu _+\). By (4.10) we directly see that \(m_*=\mu _+w_*\). Moreover, rewriting

$$\begin{aligned} M(z)=\mu _-\frac{m-\mu w}{\mu -\mu _-}\otimes \frac{m-\mu _+w}{\mu _+-\mu }+\frac{m-\mu _-w}{\mu -\mu _-}\otimes m-\sigma \end{aligned}$$

and looking at (4.10) yields

$$\begin{aligned} \lim _{j\rightarrow \infty } M(z_j)= \mu _+w_*\otimes w_*-\sigma _*. \end{aligned}$$

Thus \(\lambda _{\text {max}}(M(z_j))<e\), \(j\in {\mathbb {N}}\) implies \(z_*\in K_+'\). The case \(\mu _*=\mu _-\) can again be treated by obvious adaptations. Consequently \({\overline{U}}= K_-'\cup {\overline{U}}_0\cup K_+'\). \(\square \)

Next we introduce the most important \(\Lambda \)-directions.

Definition 4.5

Let \(z\in Z_0\). We call \({\tilde{z}}(z)\in Z\) defined by

$$\begin{aligned}&{\tilde{\mu }}=1,\quad {\tilde{w}}(z)=\frac{m-\mu w}{(\mu _+-\mu )(\mu -\mu _-)},\quad {\tilde{m}}(z)=w+(\mu _++\mu _--\mu ){\tilde{w}}(z),\\&{\tilde{\sigma }}(z)+{\tilde{q}}(z){{\,\mathrm{id}\,}}={\tilde{m}}(z)\otimes {\tilde{m}}(z)-\mu _+\mu _-{\tilde{w}}(z)\otimes {\tilde{w}}(z) \end{aligned}$$

the Muskat direction associated with z. Here the definition of \({\tilde{q}}\) and \({\tilde{\sigma }}\) is understood as decomposition into trace and traceless part. Moreover, any vector of the form \({\bar{z}}=(0,{\bar{w}},\lambda {\bar{w}},{\bar{\sigma }},{\bar{q}})\), \(\lambda \in {\mathbb {R}}\) is called an Euler direction provided it is contained in the wave cone \(\Lambda \).

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 \(\frac{m-\mu w}{(\mu _+-\mu )(\mu -\mu _-)}\), as seen in the proof of the following Lemma:

Lemma 4.6

It holds that

  1. (i)

    For any pair \(({\bar{w}},{\bar{\sigma }})\in {\mathbb {R}}^n\times \mathcal S_0^{n\times n}\), \({\bar{w}}\ne 0\), there exists \({\bar{q}}\in {\mathbb {R}}\), such that for \(\lambda \in {\mathbb {R}}{\setminus }\{0\}\) the vector \({\bar{z}}=(0,{\bar{w}},\lambda {\bar{w}},{\bar{\sigma }},{\bar{q}})\) is an Euler direction.

  2. (ii)

    The Muskat directions \({\tilde{z}}(z)\), \(z\in Z_0\) are contained in \(\Lambda \).

  3. (iii)

    For \(z\in Z_0\) define \(z_t:=z+t{\tilde{z}}(z)\), \(t\in (\mu _--\mu ,\mu _+-\mu )\). Then \({\tilde{z}}(z_t)\), \(T_\pm (z_t)\) and the traceless part \(M(z_t)^\circ \) are all independent of t.

  4. (iv)

    \(T_+(z+t{\bar{z}})=T_+(z)\) for all \(t\in {\mathbb {R}}\) and all Euler directions \({\bar{z}}\) with \({\bar{m}}=\mu _-{\bar{w}}\), as well as \(T_-(z+t{\bar{z}})=T_-(z)\) for all \(t\in {\mathbb {R}}\) and all Euler directions of the form \({\bar{z}}=(0,{\bar{w}},\mu _+{\bar{w}},{\bar{\sigma }},{\bar{q}})\).

Proof

(i) This basically has been shown in [17]. We nonetheless present the short proof here as well. Let \(({\bar{w}},{\bar{\sigma }},\lambda )\in {\mathbb {R}}^n\times {\mathcal {S}}_0^{n\times n}\times {\mathbb {R}}\), \({\bar{w}}\ne 0\), \(\lambda \ne 0\) and denote by \(P_\perp :{\mathbb {R}}^n\rightarrow {\mathbb {R}}^n\) the orthogonal projection onto \({\bar{w}}^\perp :=\{\xi \in {\mathbb {R}}^n:\ {{\bar{w}}}\cdot \xi =0\}\). Take \({\bar{q}}\in {\mathbb {R}}\), such that \(-{\bar{q}}\) is an eigenvalue of the linear map \(P_\perp \circ {\bar{\sigma }}:{\bar{w}}^\perp \rightarrow {\bar{w}}^\perp \), and let \(\xi \in {\bar{w}}^\perp {\setminus }\{0\}\) denote a corresponding eigenvector. Furthermore, we choose \(c\in {\mathbb {R}}\), such that \(({{\,\mathrm{id}\,}}-P_\perp ){\bar{\sigma }}\xi =-c\lambda {\bar{w}}\). Then one easily checks that

$$\begin{aligned} \begin{pmatrix} {\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}}&{}\quad \lambda {\bar{w}}\\ \lambda {\bar{w}}^T &{}\quad 0\\ {\bar{w}}^T &{}\quad 0 \end{pmatrix}\begin{pmatrix} \xi \\ c \end{pmatrix}=0. \end{aligned}$$

(ii) Let \(z\in Z_0\), take any element \(\xi \in {\mathbb {R}}^n{\setminus }\{0\}\) with \({\tilde{w}}(z)\cdot \xi =0\) and define \(c:=-{\tilde{m}}(z)\cdot \xi \). Then

$$\begin{aligned} \begin{pmatrix} {\tilde{\sigma }}(z)+{\tilde{q}}(z){{\,\mathrm{id}\,}}&{}\quad {\tilde{m}}(z)\\ {\tilde{m}}(z)^T &{}\quad 1\\ {\tilde{w}}(z)^T &{}\quad 0 \end{pmatrix}\begin{pmatrix} \xi \\ c \end{pmatrix}=({\tilde{m}}(z)\cdot \xi +c)\begin{pmatrix} {\tilde{m}}(z)\\ 1\\ 0 \end{pmatrix}=0. \end{aligned}$$

(iii) Let \(z\in Z_0\), \(t\in (\mu _--\mu ,\mu _+-\mu )\), \(z_t=z+t{\tilde{z}}(z)\). First of all observe that

$$\begin{aligned} (\mu _+-\mu -t)(\mu +t-\mu _-){\tilde{w}}(z_t)&=m+t{\tilde{m}}(z)-(\mu +t)(w+t{\tilde{w}}(z))\\&=m-\mu w +t(\mu _++\mu _--2\mu ){\tilde{w}}(z)-t^2{\tilde{w}}(z)\\&=(\mu _+-\mu -t)(\mu +t-\mu _-){\tilde{w}}(z). \end{aligned}$$

Hence \({\tilde{w}}(z_t)={\tilde{w}}(z)\) and

$$\begin{aligned} {\tilde{m}}(z_t)&=w+t{\tilde{w}}(z)+(\mu _++\mu _--\mu -t){\tilde{w}}(z_t)\\&=w+(\mu _++\mu _--\mu ){\tilde{w}}(z)={\tilde{m}}(z). \end{aligned}$$

The invariances \({\tilde{\sigma }}(z_t)={\tilde{\sigma }}(z)\) and \({\tilde{q}}(z_t)={\tilde{q}}(z)\) then follow by the definition of \({\tilde{\sigma }}\), \({\tilde{q}}\). Thus \({\tilde{z}}(z_t)={\tilde{z}}(z)\).

Next \(T_\pm (z_t)=T_\pm (z)\) follows immediately after rewriting

$$\begin{aligned} T_+(z)=\frac{\mu _+}{n}\left| w+(\mu _+-\mu ){\tilde{w}}(z)\right| ^2,\quad T_-(z)=\frac{\mu _-}{n}\left| w+(\mu _--\mu ){\tilde{w}}(z)\right| ^2. \end{aligned}$$

It remains to check that the traceless part of M(z) is invariant along the line segment in Muskat direction. Plugging

$$\begin{aligned} w&={\tilde{m}}(z)-(\mu _++\mu _--\mu ){\tilde{w}}(z),\\ m&=\mu w +(\mu _+-\mu )(\mu -\mu _-){\tilde{w}}(z)=\mu {\tilde{m}}(z)-\mu _-\mu _+{\tilde{w}}(z) \end{aligned}$$

into the definition of M(z) leads us to

$$\begin{aligned} M(z)&=\mu {\tilde{m}}(z)\otimes {\tilde{m}}(z)-\mu _-\mu _+({\tilde{m}}(z)\otimes {\tilde{w}}(z)+{\tilde{w}}(z)\otimes {\tilde{m}}(z))\\&\quad +\mu _-\mu _+(\mu _++\mu _--\mu ){\tilde{w}}(z)\otimes {\tilde{w}}(z)-\sigma . \end{aligned}$$

Thus for the traceless part we get

$$\begin{aligned} M(z_t)^\circ&=M(z)^\circ +t\big ({\tilde{m}}(z)\otimes {\tilde{m}}(z) -\mu _-\mu _+{\tilde{w}}(z)\otimes {\tilde{w}}(z)\big )^\circ -t{\tilde{\sigma }}(z)=M(z)^\circ . \end{aligned}$$

(iv) obviously is true, because \(m+t{\bar{m}}-\mu _{\pm }(w+t{\bar{w}})=m-\mu _\pm w\) for \({\bar{m}}=\mu _\pm {\bar{w}}\). \(\square \)

As a corollary, we obtain that any two points in K can be connected with a \(\Lambda \)-direction, up to modifying the pressure, which implies that although the wave cone \(\Lambda \) is not the whole space, it is still quite big (with respect to K).

Corollary 4.7

For any \(z_1,z_2\in K\), \(z_1\ne z_2\), one has \(z_2-z_1+(0,0,0,0,q_1-q_2)\in \Lambda .\)

Proof

In the case \(\mu _1\ne \mu _2\) we assume without loss of generality that \(\mu _1=\mu _{-}\) and \(\mu _2=\mu _{+}\). Set \({\bar{z}}=z_2-z_1+(0,0,0,0,q_1-q_2)\), such that \({\bar{q}}=0\). Similarly to (ii) from Lemma 4.6 one can prove that \({\bar{z}}\in \Lambda \) if

$$\begin{aligned} {\bar{\mu }}({\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}})={\bar{m}}\otimes {\bar{m}}+\gamma {\bar{w}}\otimes {\bar{w}}, \end{aligned}$$
(4.11)

for some \(\gamma \in {\mathbb {R}}.\)

Since \(z_i\in K\), we have

$$\begin{aligned} \sigma _i=\mu _i w_i\otimes w_i-e{{\,\mathrm{id}\,}}, \end{aligned}$$

for \(i=1,2.\) Therefore, we obtain that

$$\begin{aligned} {\bar{\sigma }}=({\bar{\sigma }}+{\bar{q}}{{\,\mathrm{id}\,}})=\mu _2 w_2\otimes w_2-\mu _1 w_1\otimes w_1. \end{aligned}$$

Through a simple calculation one can then show that (4.11) holds for \(\gamma =-\mu _{-}\mu _{+}\).

If \(\mu _1=\mu _2\), recall that in the proof of Lemma 4.6 (i) a suitable pressure \({\bar{q}}\) has been chosen to be an eigenvalue of \(-P_\perp \circ {\bar{\sigma }}:{\bar{w}}^\perp \rightarrow {\bar{w}}^\perp \). But \(z_1,z_2\in K\) in fact implies that \(P_\perp \circ {\bar{\sigma }}\) vanishes on all of \({\bar{w}}^\perp \) and we can conclude the statement. \(\square \)

Recall the definition of \(\pi :Z\rightarrow {\mathbb {R}}\times {\mathbb {R}}^n\times {\mathbb {R}}^n\times {{\mathcal {S}}}_0^{n\times n}\) in (2.6).

Lemma 4.8

The projection \(\overline{\pi (U)}\) is bounded in terms of e, \(\mu _\pm \), n and hence compact. Moreover, for every \(z\in {\overline{U}}{\setminus }K\) there exists \({\bar{z}}\in \Lambda {\setminus }\{0\}\), such that \(z\pm {\bar{z}}\in {\overline{U}}\).

Proof

We first prove that \(\pi (U)\) is bounded in terms of \(e,\mu _-,\mu _+\) and the dimension n. Let \(z\in U\). Obviously \(\mu \in (\mu _-,\mu _+)\) is bounded. The inequalities (4.10) imply that there exists a constant \(c=c(e,\mu _-,\mu _+,n)>0\), such that

$$\begin{aligned} \left| m-\mu _-w\right| \leqq c(\mu -\mu _-),\quad \left| m-\mu _+w\right| \leqq c(\mu _+-\mu ). \end{aligned}$$
(4.12)

Adapting the constant when necessary we obtain

$$\begin{aligned} \left| m\right| =\left| \frac{\mu _+}{\mu _+-\mu _-}(m-\mu _-w)-\frac{\mu _-}{\mu _+-\mu _-} (m-\mu _+w)\right| \leqq c, \end{aligned}$$

which then also implies \(\left| w\right| \leqq c\). Next observe that the matrix M(z) can be rewritten to

$$\begin{aligned} M(z)&=-\mu \frac{m-\mu _-w}{\mu -\mu _-}\otimes \frac{m-\mu _+w}{\mu -\mu _+} +\frac{m-\mu _-w}{\mu -\mu _-}\otimes m +m\otimes \frac{m-\mu _+w}{\mu -\mu _+}-\sigma . \end{aligned}$$

Hence \(M(z)+\sigma \) is uniformly bounded by (4.12). As a consequence we obtain \(\left| {{\,\mathrm{tr}\,}}M(z)\right| \leqq c\). This bound on the trace together with \(\lambda _{\text {max}}(M(z))=Q(z)<e\), due to the fact that \(z\in U\), gives us a uniform bound on the whole spectrum of M(z). Therefore \(M(z)+\sigma \) and M(z) are both uniformly bounded. Consequently \(\left| \sigma \right| \leqq c\), and \(\overline{\pi (U)}\) is compact.

Next we show that any \(z\in {\overline{U}}{\setminus }K\) can be perturbed along a \(\Lambda \)-segment without leaving \({\overline{U}}\). Recall that \({\overline{U}}={\overline{U}}_0\cup K_+'\cup K_-'\) and \(K\subset K_+'\cup K_-'\) by Lemma 4.4.

If \(z\in K_+'{\setminus }K\), we can find similarly as in [17] a suitable Euler direction \({\bar{z}}=(0,{\bar{w}},\mu _+{\bar{w}},{\bar{\sigma }},{\bar{q}})\in \Lambda \) such that \(z+t{\bar{z}}\in K_+'\) for \(\left| t\right| \) small enough. Indeed, by a change of basis we can restrict ourselves to the case that \(\mu _+w\otimes w-\sigma \) is diagonal. Denote the entries by \(\lambda _1\geqq \lambda _2 \geqq \cdots \geqq \lambda _n\), where \(\lambda _1 \leqq e\) and \(\lambda _n<e\). Let \(e_1,\ldots ,e_n\) denote the canonical basis of \({\mathbb {R}}^n\). We take \({\bar{w}}=e_n\) and

$$\begin{aligned} {\bar{\sigma }}=\mu _+e_n\otimes w+\mu _+w\otimes e_n-\alpha e_n\otimes e_n, \end{aligned}$$

where \(\alpha =2\mu _+w_n\) makes \({\bar{\sigma }}\) trace free. It follows that

$$\begin{aligned}&\mu _+(w+t{\bar{w}})\otimes (w+t{\bar{w}})-(\sigma +t{\bar{\sigma }})\\&\quad =\sum _{j=1}^n\lambda _je_j\otimes e_j+t(\mu _+e_n\otimes w+\mu _+w\otimes e_n-{\bar{\sigma }})+t^2\mu _+e_n\otimes e_n\\&\quad =\sum _{j=1}^{n-1}\lambda _je_j\otimes e_j + (\lambda _n+\alpha t+\mu _+t^2)e_n\otimes e_n. \end{aligned}$$

Clearly, \(\lambda _j\leqq e\), \(j=1,\ldots ,n-1\) and \(\lambda _n+\alpha t+\mu _+t^2\leqq e\) for all \(\left| t\right| \) small enough, since the inequaltiy holds strict for \(t=0\).

The same reasoning applies also to the case \(z\in K_-'{\setminus }K\).

Now let \(z\in {\overline{U}}_0\). If \(Q(z)<e\) or if \(T_+(z)=T_-(z)\) we take the Muskat direction \({\bar{z}}={\tilde{z}}(z)\). Because then \( T_\pm (z+t{\tilde{z}}(z))=T_\pm (z)\leqq e \), \(t\in (\mu _--\mu ,\mu _+-\mu )\) by Lemma 4.6 (iii). Moreover, a straightforward computation shows

$$\begin{aligned} Q(z)&=\frac{1}{n}{{\,\mathrm{tr}\,}}M(z)+\lambda _{\text {max}}(M(z)^\circ )\\&=\frac{\mu _+-\mu }{\mu _+-\mu _-}T_-(z)+\frac{\mu -\mu _-}{\mu _+-\mu _-}T_+(z) +\lambda _{\text {max}}(M(z)^\circ ) \end{aligned}$$

and thus by Lemma 4.6 (iii) we have

$$\begin{aligned} Q(z+t{\tilde{z}}(z))=Q(z)+t\frac{T_+(z)-T_-(z)}{\mu _+-\mu _-}. \end{aligned}$$

For \(\left| t\right| \left| T_+(z)-T_-(z)\right| \leqq (e-Q(z))(\mu _+-\mu _-)\) and \(\left| t\right| < {{\,\mathrm{dist}\,}}(\mu ,\left\{ \,\mu _-,\mu _+\,\right\} )\) we therefore conclude \(z+t{\tilde{z}}(z)\in {\overline{U}}_0\).

From now on we consider the remaining case \(Q(z)=e\) and \(T_+(z)\ne T_-(z)\). Note that then necessarily \(\lambda _{\text {min}}(M(z))<e\), because otherwise \(e=\lambda _{\text {max}}(M(z))=\lambda _{\text {min}}(M(z))\) yields \(M(z)^\circ =0\) and thus

$$\begin{aligned} e=Q(z)=\frac{\mu _+-\mu }{\mu _+-\mu _-}T_-(z)+\frac{\mu -\mu _-}{\mu _+-\mu _-}T_+(z). \end{aligned}$$

Since \(T_+(z)\leqq e\), \(T_-(z)\leqq 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 form \({\bar{z}}=(0,{\bar{w}},\mu _+{\bar{w}},{\bar{\sigma }},{\bar{q}})\), where \({\bar{w}}\in {\mathbb {R}}^n\) and \({\bar{\sigma }}\in {\mathcal {S}}_0^{n\times n}\) will be chosen later and \({\bar{q}}={\bar{q}}({\bar{w}},{\bar{\sigma }})\) by Lemma 4.6 (i). These Euler directions allow us to preserve \(T_-\) due to Lemma 4.6 (iv), that is, \(T_-(z+t{\bar{z}}_+)=T_-(z)\leqq e\) for all \(t\in {\mathbb {R}}\).

Now we need to guarantee that \(Q(z+t{\bar{z}})=Q(z)=e\) for small enough \(\left| t\right| \) and some choice of \({\bar{w}}\), \({\bar{\sigma }}\). As in the cases \(z\in K_\pm '{\setminus }K\) we can again assume that the matrix M(z) is diagonal with entries \(e=\lambda _1\geqq \lambda _2\geqq \cdots \geqq \lambda _n\) and \(\lambda _n<e\). As before we take \({\bar{w}}=e_n\), \({\bar{m}}=\mu _+e_n\) and the uniquely determined pair \(({\bar{\sigma }},\alpha )\in {\mathcal {S}}_0^{n\times n}\times {\mathbb {R}}\) satisfying

$$\begin{aligned} M(z+t{\bar{z}})&=M(z)+\alpha t e_n\otimes e_n+\frac{\mu _+(\mu _+-\mu _-)}{\mu -\mu _-}t^2e_n\otimes e_n. \end{aligned}$$

For small enough \(\left| t\right| \) we therefore conclude that this Euler perturbation does not affect the maximal eigenvalue, that is, \(Q(z+t{\bar{z}})=Q(z)=e\) for \(\left| t\right| \) small.

Furthermore, the last condition needed for \(z+t{\bar{z}}\in {\overline{U}}\) simply follows by the continuity of \(T_+\), that is, for all \(\left| t\right| \) small enough it holds that

$$\begin{aligned} T_+(z+t{\bar{z}})< T_-(z)\leqq e. \end{aligned}$$

\(\square \)

Now we have all ingredients for the proof of \(K^\Lambda =K^{co}={\overline{U}}\) at hand.

Proof of Proposition 4.2

Lemma 4.4 implies \(K^\Lambda \subset K^{co}\subset {\overline{U}}\), while Lemma 4.8 says that the \(\Lambda \)-extreme points of the up to the q-component compact set \({\overline{U}}\) are contained in K. The inclusion \({\overline{U}}\subset K^\Lambda \) follows by the Krein-Milman theorem for \(\Lambda \)-convex sets, cf. [23], Lemma 4.16. \(\square \)

4.3 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 \(\pi \) defined in (2.6). We have the following result:

Lemma 4.9

For any \(z\in U\) there exists \({\bar{z}}\in \Lambda \) such that we have

$$\begin{aligned} {[}z-{\bar{z}},z+{\bar{z}}]\subset U\ \text {and } |\pi ({\bar{z}})|\geqq \frac{1}{2N} d(\pi (z),\pi (K)), \end{aligned}$$

where \(N=\dim (Z)\) and d denotes the Euclidian distance on \(\pi (Z)\).

Proof

We proceed as in the proof of Lemma 4.7 from [17]. Since \(z\in U=\text {int}K^{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 \(\lambda _i\in (0,1)\), \(z_i\in K\), \(i\in {1,\ldots ,N+1}\), \(\sum _i \lambda _i=1\), such that

$$\begin{aligned} z=\sum _{i=1}^{N+1}\lambda _i z_i. \end{aligned}$$

We may also assume that the coefficients are ordered such that \(\lambda _1=\max _i{\lambda _i}\), then for any \(j>1\) we have

$$\begin{aligned} z\pm \frac{1}{2}\lambda _j(z_j-z_1)\in \text {int}K^{co}. \end{aligned}$$

Indeed, one may rewrite \(z\pm \frac{1}{2}\lambda _j(z_j-z_1)=\sum _i \kappa _i z_i\), where \(\kappa _1=\lambda _1\mp \frac{1}{2}\lambda _j\), \(\kappa _j=\lambda _j\pm \frac{1}{2}\lambda _j\) and \(\kappa _i=\lambda _i\) for \(i\not \in \{1,j\}\), such that these coefficients are in (0, 1).

Furthermore, since we have \(z-z_1=\sum _{i=2}^{N+1}\lambda _i(z_i-z_1)\), it follows that

$$\begin{aligned} |\pi (z)-\pi (z_1)|\leqq N \max _{i=2,\ldots ,N+1} \lambda _i |\pi (z_i)-\pi (z_1)|. \end{aligned}$$
(4.13)

Choose \(j>1\) such that \(\max _{i=2,\ldots ,N+1} \lambda _i |\pi (z_i)-\pi (z_1)|=\lambda _j|\pi (z_j)-\pi (z_1)|\), and let

$$\begin{aligned} {\bar{z}}=\frac{1}{2}\lambda _j(z_j-z_1). \end{aligned}$$

Then \([z-{\bar{z}},z+{\bar{z}}]\subset \text {int} K^{co}\) and

$$\begin{aligned} d(\pi (z),\pi (K))\leqq |\pi (z)-\pi (z_1)|\leqq 2N |\pi ({\bar{z}})|. \end{aligned}$$

To conclude the proof of the lemma, it would suffice to have \({\bar{z}}\in \Lambda \). 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 in \({\bar{z}}\). However, since the constraints in K, and respectively the inequalities in U do not involve the pressure, this can be done such that \(z\pm {\bar{z}}\in \text {int} K^{co}\) still remains valid. This concludes the proof. \(\square \)

4.4 Continuity of Constraints

We now go back to the \((x,t)\in {\mathscr {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({\mathscr {D}})\).

Lemma 4.10

Let \({\mathscr {U}}\subset {\mathscr {D}}\) be an open, bounded set and \(e:\Omega \times [0,T)\rightarrow {\mathbb {R}}_+\). If the map \((x,t)\mapsto e(x-1/2gt^2e_n,t)\in {\mathbb {R}}_+\) is continuous and bounded on \({\mathscr {U}}\), then it follows that the map \((x,t)\mapsto \pi (K_{(x,t)})\) is continuous and bounded on \({\mathscr {U}}\) with respect to the Hausdorff metric \(d_{{\mathcal {H}}}\).

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\subset {\mathbb {R}}^l\) for some \(l\in {\mathbb {N}}\) are compact sets and \(r>0\) such that

  • for any \(z\in A\) there exists \(z'\in B\cap B_r(z),\)

  • for any \(z\in B\) there exists \(z'\in A\cap B_r(z).\)

Then \(d_{{\mathcal {H}}}(A,B)< r\).

Proof

See [12]. \(\square \)

Proof of Lemma 4.10

Fix \(y=(x,t)\in {\mathscr {U}}\). For \(\varepsilon >0\) there exists \(\delta >0\) such that

$$\begin{aligned} \left| e(y)-e(y')\right|<\varepsilon \quad \text {and}\quad \left| \left( n\frac{e(y)}{\mu _\pm }\right) ^{1/2}-\left( n\frac{e(y')}{\mu _\pm }\right) ^{1/2} \right| <\varepsilon , \end{aligned}$$
(4.14)

for any \(y'\in B_\delta (y)\subset {\mathscr {U}}.\) Using Lemma 4.11 we will prove \(d_{{\mathcal {H}}}( \pi (K_y), \pi (K_{y'}))< c \varepsilon \) for any \(y'\in B_\delta (y)\subset {\mathscr {U}}\), with \(c>0\) depending only on \(\mu _+\), \(\mu _-\) and n.

Let

$$\begin{aligned} z=(\mu ,w,\mu w,\mu w\otimes w-e(y){{\,\mathrm{id}\,}},q)\in K_y, \end{aligned}$$

with \(\mu \in \{\mu _+,\mu _-\}\) and \(\mu |w|^2=ne(y)\). It follows that

$$\begin{aligned} w=\left( \frac{n}{\mu }e(y)\right) ^{1/2}b, \end{aligned}$$

for some \(b\in S^{n-1}\).

We define

$$\begin{aligned}z'=(\mu ,w',\mu w', \mu w'\otimes w'-e(y'){{\,\mathrm{id}\,}},q)\end{aligned}$$

by setting

$$\begin{aligned} w'=\left( \frac{n}{\mu }e(y')\right) ^{1/2}b. \end{aligned}$$

Note that \(z'\in K_{y'}\).

Furthermore, from (4.14) it follows that

$$\begin{aligned} |w-w'|<\varepsilon , \end{aligned}$$

from which one can conclude \(\left| z-z'\right| <c\varepsilon \) for some \(c=c(\mu _\pm ,n)>0\).

Due to the symmetry of this construction, one can similarly prove that for any \(z'\in K_{y'}\) there exists \(z\in K_y\) such that \(|z-z'|<c\varepsilon .\) The result then follows from Lemma 4.11.

The boundedness of \(\bigcup _{y\in {\mathscr {U}}}\pi (K_y)\) follows from Lemma 4.8 and the assumption that the function e is bounded. \(\square \)

5 Proof of Theorem 2.4

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

5.1 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:\Omega \times (0,T)\rightarrow {\mathbb {R}}_+\) with \((x,t)\mapsto e\big (x-\frac{1}{2}gt^2e_n,t\big )\) being continuous on an open set \({\mathscr {U}}\subset {\mathscr {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:{\mathscr {D}}\rightarrow 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 \({\mathscr {D}}\), \(\pi (z)\) is continuous in \({\mathscr {U}}\), \(z(y)\in K_y\) holds for almost every \(y\in {\mathscr {D}}{\setminus }{\mathscr {U}}\) and

$$\begin{aligned} z(y)\in U_y=\text {int}K_y^{co}\text { for any }y\in {\mathscr {U}}. \end{aligned}$$
(5.1)

We have the following convex integration result:

Theorem 5.2

Suppose that there exists a subsolution \(z_0\) in the sense of Definition 5.1. Then there exist infinitely many weak solutions \(z:{\mathscr {D}}\rightarrow Z\) of (3.5) which coincide with \(z_0\) almost everywhere in \({\mathscr {D}}{\setminus }{\mathscr {U}}\), satisfy \(z(y)\in K_y\) almost everywhere in \({\mathscr {D}}\), and for every open ball \(B\subset {\mathscr {U}}\) the solutions satisfy the mixing property

$$\begin{aligned} \int _B\mu _+-\mu (x,t)\,d(x,t)\cdot \int _B\mu (x,t)-\mu _-\,d(x,t)>0. \end{aligned}$$
(5.2)

Furthermore, among these weak solutions there exists a sequence \(\{z_k\}_{k\geqq 1}\) such that \(\pi (z_k)\) converges weakly to \(\pi (z_0)\) in \(L^2({\mathscr {U}};\pi (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 \(\pi \). 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.

Lemma 5.3

Suppose that there exists a subsolution z such that

$$\begin{aligned}\int _{\mathscr {U}} d(\pi (z(y)),\pi (K_y))^2\, \mathrm{d}y = \varepsilon >0.\end{aligned}$$

Then there exist \(\delta =\delta (\varepsilon )>0\) and a sequence of subsolutions \(\{z_k\}_{k\geqq 0}\) such that

  • \(z_k=z\) in \({\mathscr {D}}{\setminus }{\mathscr {U}}\), for any \(k\geqq 0\),

  • \(\int _{\mathscr {U}} |\pi (z_k(y)-z(y))|^2\, \mathrm{d}y \geqq \delta ,\) for any \(k\geqq 0\),

  • \(\pi (z_k)\rightharpoonup \pi (z)\) in \(L^2({\mathscr {U}};\pi (Z))\) as \(k\rightarrow +\infty \).

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\subset {\mathbb {R}}^n\) be a compact set. Then for any compact set \(C\subset \text {int}K^{co}\) there exists \(\varepsilon >0\) such that for any compact set \(K'\subset {\mathbb {R}}^n\) with \(d_{\mathcal {H}}(K,K')<\varepsilon \) we have \(C\subset \text {int}(K')^{co}\).

Proof of Lemma 5.3

Fix \(y\in {\mathscr {U}}\). From Lemma 4.9 it follows that there exists some \(C'>0\) independent of y and z, and some \({\bar{z}}(y)\in \Lambda \) such that

$$\begin{aligned} {[}z(y)-{\bar{z}}(y),z(y)+{\bar{z}}(y)]\subset U_y,\quad |\pi ({\bar{z}}(y))|\geqq C' d(\pi (z(y)),\pi (K_y)). \end{aligned}$$

Now Lemma 4.10, the continuity of \(\pi (z)\) and Lemma 5.4 applied to the projected sets imply that there exist \(r(y),R(y)>0\) such that

$$\begin{aligned} {[}z(y')-{\bar{z}}(y),z(y')+{\bar{z}}(y)]+\overline{B_{R(y)}(0)}\subset U_{y'},\\ d(\pi (z(y')),\pi (K_{y'})) \leqq 2 d(\pi (z(y)),\pi (K_y)), \end{aligned}$$

for any \(y'\in B_{r(y)}(y).\)

Using Lemma 4.1, we find a sequence \(\{z_{y,N}\}_{N\geqq 0}\subset C^\infty _c(B_1(0))\) solving (3.5) such that

  • \(z_{y,N}(y')\in [-{\bar{z}}(y),{\bar{z}}(y)]+\overline{B_{R(y)}(0)}\) for all \(y'\in B_1(0)\), \(N\geqq 0,\)

  • \(z_{y,N}\rightharpoonup 0\) in \(L^2\),

  • \(\int |\pi (z_{y,N})|^2 \, \mathrm{d}{\tilde{y}}\geqq C|\pi ({\bar{z}}(y))|^2\) for all \(N\geqq 0.\)

From here on the proof is the same as Step 2 of the proof of Lemma 2.4 from [12], using a standard covering argument, therefore the details are left to the reader. \(\square \)

Proof of Theorem 5.2

Let

$$\begin{aligned} X_0&=\left\{ z'\in L^2({\mathscr {D}};\pi (Z))\text { such that }z'=\pi (z)\text { for some subsolution } z \text { in the sense}\right. \\&\quad \left. \text {of Definition }5.1 \text { satisfying }z=z_0\text { on }{\mathscr {D}}{\setminus }{\mathscr {U}}\right\} , \end{aligned}$$

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(\cdot ,\cdot )\). 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'\in 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')=\int _{\mathscr {U}} |z'|^2\, \mathrm{d}y\) is a Baire-1 function on X. Indeed, setting

$$\begin{aligned} I_j(z')=\int _{\left\{ \,y\in {\mathscr {U}}: d(y,\partial {\mathscr {U}})>1/j\,\right\} } |(z'*\chi _j)(y)|^2\ \mathrm{d}y, \end{aligned}$$

where \(\chi _j\in C_c^\infty (B_{1/j}(0))\) is the standard mollifying sequence, one obtains that \(I_j\) is continuous on X and that \(I_j(z')\rightarrow I(z')\) as \(j\rightarrow +\infty .\)

It follows from the Baire category theorem that the set

$$\begin{aligned} Y=\{z'\in X:\ I\text { is continuous at }z'\} \end{aligned}$$

is residual in X. We claim that for any \(z'\in Y\) it follows that

$$\begin{aligned} J(z'):=\int _{\mathscr {U}} d(z'(y),\pi (K_y))^2\, \mathrm{d}y=0. \end{aligned}$$

Suppose the contrary. Then \(J(z')=\varepsilon >0\) for some \(z'\in Y\), and let \(z'_j\in 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\rightarrow z'\) strongly in \(L^2({\mathscr {U}};\pi (Z))\). Note that J is continuous with respect to the strong \(L^2\)-topology. Therefore we may assume that \(J(z'_j)>\varepsilon /2\) for all \(z'_j\).

Since \(z'_j\in X_0\), there exists some \(z_j:{\mathscr {D}}\rightarrow Z\) which is a subsolution in the sense of Definition 5.1 and such that \(z'_j=\pi (z_j).\) We may then apply Lemma 5.3 to deduce that there exists \(\delta =\delta (\varepsilon )>0\) and a subsolution \({\tilde{z}}_j\) such that \(\int _{\mathscr {U}} |\pi (z_j(y)-{\tilde{z}}_j(y))|^2\, \mathrm{d}y \geqq \delta \) and \(\pi (z_j-{\tilde{z}}_j)\rightharpoonup 0\) weakly in \(L^2\). Since \(z'_j=\pi (z_j)\rightarrow z'\) and \(z'\in Y\), we conclude as before \(\pi ({\tilde{z}}_j)\rightarrow z'\) strongly in \(L^2\) contradicting the fact that \(\pi ({\tilde{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 \({\mathscr {U}}\). The set

$$\begin{aligned} X_B^{\mu _\pm }=\left\{ \,z'\in X:\int _B\mu _\pm -\mu (x,t)\,d(x,t)=0\,\right\} \end{aligned}$$

is closed in X and has empty interior, since \(X_B^{\mu _\pm }\subset X{\setminus }X_0\). Therefore \(J^{-1}(0){\setminus }X_B^{\mu _\pm }\) is residual in X, as is \(J^{-1}(0){\setminus }\left( \bigcup _i\big (X_{B_i}^{\mu _+}\cup X_{B_i}^{\mu _-}\big )\right) \) for any countable union of balls \(B_i\subset {\mathscr {U}}\). By taking all balls \((B_i)_{i\in {\mathbb {N}}}\subset {\mathscr {U}}\) with rational centers and radii we can conclude the statement. \(\square \)

5.2 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 \(\Omega \subset {\mathbb {R}}^n\) and \(T>0\) we defined \({\mathscr {D}}=\left\{ \,(y,t)\in {\mathbb {R}}^n\times (0,T):y-\frac{1}{2}gt^2e_n\in \Omega \,\right\} \).

Now let \(z=(\mu ,w,m,\sigma ,q)\) be a weak solution of (3.5) with some suitable initial data. Defining again \(y:=x+\frac{1}{2}gt^2e_n\), as well as

$$\begin{aligned} \begin{aligned} \rho (x,t)&=\mu \left( y,t\right) ,\\ v(x,t)&=w\left( y,t\right) -g t e_n,\\ u(x,t)&=m\left( y,t\right) -\mu \left( y,t\right) g t e_n,\\ P(x,t)&=q\left( y,t\right) +gt\frac{1}{n}\left( gt\mu (y,t)-2m_n\left( y,t\right) \right) ,\\ S(x,t)&=\sigma \left( y,t\right) -gt \left( m\left( y,t\right) \otimes e_n+e_n\otimes m\left( y,t\right) \right) \\&\quad +g^2t^2\mu (y,t)e_n\otimes e_n -\left( gt\frac{1}{n}\left( gt\mu (y,t)-2m_n\left( y,t\right) \right) \right) {{\,\mathrm{id}\,}}, \end{aligned} \end{aligned}$$
(5.3)

one obtains through lenghty but straightforward calculations that \((\rho ,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:\Omega \times [0,T)\rightarrow {\mathbb {R}}_+\) the condition \(z(y,t)\in K_{(y,t)}\) for \(y=x+\frac{1}{2}gt^2\), \((x,t)\in \Omega \times (0,T)\) and with \(K_{(y,t)}\) defined in (3.6) translates to \((\rho ,v,u,S,P)(x,t)\in {{\mathcal {K}}}_{(x,t)}\) with \({{\mathcal {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)\in U_{(y,t)}\) translates to \((\rho ,v,u,S,P)(x,t)\in {{\mathcal {U}}}_{(x,t)}\) with \({{\mathcal {U}}}_{(x,t)}\) defined in (2.3),(2.4). Since the transformation is an affine bijection, we also see that \({{\mathcal {U}}}_{(x,t)}\) is the interior of the convex hull of \({{\mathcal {K}}}_{(x,t)}\).

We have now all pieces together to prove our main result.

Proof of Theorem 2.4

Let \(z=(\rho ,v,u,S,P):\Omega \times (0,T)\rightarrow Z\) be a subsolution (in the sense of Definition 2.3) of (1.1) associated with \(e:\Omega \times [0,T)\rightarrow {\mathbb {R}}_+\) bounded and initial data \((\rho _0,v_0)\in L^\infty (\Omega )\times L^2(\Omega ;{\mathbb {R}}^n)\) satisfying (1.2). We also define the transformed mixing zone

$$\begin{aligned} {\mathscr {U}}'=\left\{ \,(y,t)\in {\mathbb {R}}^n\times (0,T):\left( y-\frac{1}{2}gt^2e_n,t\right) \in {\mathscr {U}}\,\right\} . \end{aligned}$$

The inverse of the transformation (5.3) applied to z gives us a weak solution of (3.5) (in the sense of Definition 3.1) which we call \(z'=(\mu ,w,m,\sigma ,q):{\mathscr {D}}\rightarrow Z\). By the discussion of this section and Definition 2.3, \(\pi (z')\) is continuous on \({\mathscr {U}}'\), \(z'(y,t)\in U_{(y,t)}=\text {int} K_{(y,t)}^{co}\) for all \((y,t)\in {\mathscr {U}}'\) and \(z'(y,t)\in K_{(y,t)}\) for almost every \((y,t)\in {\mathscr {D}}{\setminus }{\mathscr {U}}'\).

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 \({\mathscr {U}}'\)). Theorem 5.2 therefore provides us with infinitely many solutions of our differential inclusion (3.5), (3.6) which outside of \({\mathscr {U}}'\) agree with \(z'\) and inside \({\mathscr {U}}'\) satisfy the mixing property (5.2), as well as with a sequence of solutions such that \(\pi (z'_k)\) converges \(L^2\)-weakly to \(\pi (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 \(\pi (z)\) is in addition of class \({{\mathcal {C}}}^0([0,T];L^2(\Omega ;\pi (Z)))\). As before let \(z'\) be the corresponding transformed subsolution defined on \({\mathscr {D}}\). Due to an improvement of the Tartar framework as in [7, 17] one can show that the induced sequence \(\{\pi (z'_k)\}_{k\in {\mathbb {N}}}\) not only converges weakly in \(L^2({\mathscr {D}})\) to \(\pi (z')\), but weakly on every time-slice \({\mathscr {D}}(t)\) uniformly in \(t\in [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

$$\begin{aligned} E_k(x,t):=\frac{n}{2} e(x,t)- gt e_n\cdot u_k(x,t)-\frac{1}{2}\rho _k (x,t)g^2t^2+\rho _k(x,t) gx_n \end{aligned}$$

satisfy

$$\begin{aligned} \int _\Omega E_k(x,t)\,\mathrm{d}x\rightarrow \int _\Omega E_{sub}(x,t)\,\mathrm{d}x \end{aligned}$$

uniformly in \(t\in [0,T]\) as \(k\rightarrow \infty \), 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

$$\begin{aligned} \varepsilon (t):=\int _\Omega E_{sub}(x,0)-E_{sub}(x,t)\,\mathrm{d}x>0,\quad t\in (0,T) \end{aligned}$$

goes to 0 as \(t\searrow 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

$$\begin{aligned} \left| \int _\Omega gt e_n\cdot (u(x,t)-u_k(x,t))+\left( \frac{1}{2}g^2t^2- gx_n\right) (\rho (x,t)-\rho _k (x,t))\, \mathrm{d}x \right| \leqq \varepsilon (t), \end{aligned}$$

for all \(t\in [0,T],\ k\geqq 0\). The statement follows. \(\square \)

6 Subsolutions

We now turn to the construction of Rayleigh–Taylor subsolutions. We start by observing that the relaxation inside the mixing zone \({\mathscr {U}}\subset \Omega \times (0,T)\) given in Definition 2.3 can be equivalently rewritten (in the spirit of [6]) as the system

$$\begin{aligned} \begin{aligned} \partial _t (\rho v+f)+\text {div } S + \nabla P&= - \rho g e_n,\\ \text {div } v&=0,\\ \partial _t \rho + \text {div }(\rho v+f)&=0, \end{aligned} \end{aligned}$$
(6.1)

where

$$\begin{aligned} f&:=\frac{\rho _+-\rho }{\rho _+-\rho _-}\sqrt{\frac{ne}{\rho _+}} (\rho -\rho _-)\xi +\frac{\rho -\rho _-}{\rho _+-\rho _-} \sqrt{\frac{ne}{\rho _-}}(\rho _+-\rho )\eta , \end{aligned}$$

for some functions \(\xi ,\eta :{\mathscr {U}}\rightarrow {\mathbb {R}}^n\) satisfying

$$\begin{aligned} \sqrt{ne}\left( \frac{\rho -\rho _-}{\sqrt{\rho _+}}\xi -\frac{\rho _+-\rho }{\sqrt{\rho _-}}\eta \right) =(\rho _+-\rho _-)(v+gte_n),\quad |\xi |< 1,\quad |\eta |< 1 \end{aligned}$$
(6.2)

in \({\mathscr {U}}\). The condition on \(\lambda _{\max }(A(z))\) from (2.4) with u replaced by \(\rho v+f\) is kept in accordance with Definition 2.3.

Indeed, in order to see this, given a subsolution \(z=(\rho ,v,u,S,P)\) it suffices to set

$$\begin{aligned} \xi :=\sqrt{\frac{\rho _+}{ne}}\frac{u-\rho _-v+(\rho -\rho _-)gte_n}{\rho -\rho _-},\quad \eta :=\sqrt{\frac{\rho _-}{ne}}\frac{u-\rho _+v+(\rho -\rho _+)gte_n}{\rho _+-\rho }. \end{aligned}$$
(6.3)

Conversely, given f, it suffices to set \(u:=\rho 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 \(\Omega \subset {\mathbb {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 \(\xi ,\eta \) 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_2e_2\), \(\xi =\xi _2e_2\), \(\eta =\eta _2e_2\) respectively. We further assume \(v\equiv 0\).

If we have chosen \(\xi \), \(\eta \), then condition (6.2) implies that e in the mixing zone is determined by

$$\begin{aligned} \sqrt{2e}=\frac{\sqrt{\rho _-\rho _+}(\rho _+-\rho _-)gt}{\sqrt{\rho _-} (\rho -\rho _-)\xi _2-\sqrt{\rho _+}(\rho _+-\rho )\eta _2}. \end{aligned}$$
(6.4)

Note also that under condition (6.2) the denominator will always be positive for \(t>0\). Outside the mixing zone we will have \(e=\frac{1}{2}\rho g^2t^2\) in accordance with (2.5).

The last equation in (6.1) then becomes

$$\begin{aligned} \partial _t \rho + gt \partial _{x_2}\left( \frac{(\rho _+-\rho )(\rho -\rho _-) (\sqrt{\rho _-}\xi _2+\sqrt{\rho _+}\eta _2)}{(\rho -\rho _-) \sqrt{\rho _-}\xi _2-(\rho _+-\rho )\sqrt{\rho _+}\eta _2} \right) =0. \end{aligned}$$
(6.5)

We recall (1.3) and hence that \(\rho _0\) only depends on the sign of \(x_2\). Using the change of coordinates \(\rho (x,t)=y(x,gt^2/2)\) and interpreting the \(\xi _2,\eta _2\) as functions of \(\rho \) only, one obtains equivalently

$$\begin{aligned} \partial _t y+ \partial _{x_2}G(y)=0, \end{aligned}$$
(6.6)

with

$$\begin{aligned} G(y)=\frac{(\rho _+-y)(y-\rho _-)(\sqrt{\rho _-}\xi _2(y) +\sqrt{\rho _+}\eta _2(y))}{(y-\rho _-)\sqrt{\rho _-}\xi _2(y)-(\rho _+-y) \sqrt{\rho _+}\eta _2(y)}. \end{aligned}$$

Now if \(G:[\rho _-,\rho _+]\rightarrow {\mathbb {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 \(\rho _0\) to obtain that

$$\begin{aligned} \begin{aligned} \rho (x_2,t)= \left\{ \begin{array}{lll} \rho _-,\text { when }x_2\leqq \frac{1}{2}gt^2 G'(\rho _-),\\ (G')^{-1}\left( \frac{2x_2}{gt^2}\right) ,\text { when }x_2\in \left( \frac{1}{2}gt^2 G'(\rho _-),\frac{1}{2}gt^2 G'(\rho _+)\right) ,\\ \rho _+,\text { when }x_2\geqq \frac{1}{2}gt^2 G'(\rho _+). \end{array} \right. \end{aligned} \end{aligned}$$
(6.7)

Observe that this already implies that the height of the mixing zone grows (up to a constant) like \(\frac{1}{2}gt^2\), more precisely we will have

$$\begin{aligned} {\mathscr {U}}=\left\{ \,(x,t)\in \Omega \times (0,T):\frac{1}{2}gt^2G' (\rho _-)<x_2<\frac{1}{2}gt^2 G'(\rho _+)\,\right\} . \end{aligned}$$
(6.8)

It is easy to check that if one is able to choose \(\xi _2,\eta _2\in (-1,1)\) such that G is indeed uniformly strictly convex and the above entropy solution exists, then defining

$$\begin{aligned} u_2(x_2,t)&:=gt G(\rho (x_2,t)),\quad S:=\frac{(\rho _++\rho _--\rho )u_2^2}{2(\rho _+-\rho )(\rho -\rho _-)}\begin{pmatrix} -1 &{} 0\\ 0 &{} 1 \end{pmatrix},\\ P(x_2,t)&:=S_1(x_2,t)-\int _{\frac{1}{2}gt^2 G'(\rho _-)}^{x_2}\partial _t u_2(x',t)-\rho (x',t)g\,\mathrm{d}x_2, \end{aligned}$$

with \(u_2\) and S extended by 0 outside \({\mathscr {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

$$\begin{aligned} e&>\frac{(\rho _++\rho _--\rho )u_2^2}{2(\rho _+-\rho )(\rho -\rho _-)}+gt u_2+\frac{1}{2}\rho g^2t^2\\&=\frac{\rho _+-\rho }{\rho _+-\rho _-}\frac{\rho _-}{2}\left( \frac{u_2}{\rho -\rho _+} +gt\right) ^2+\frac{\rho -\rho _-}{\rho _+-\rho _-}\frac{\rho _+}{2} \left( \frac{u_2}{\rho -\rho _-}+gt\right) ^2, \end{aligned}$$

which holds, since by our reformulation inequalities (2.3) are automatically satisfied for \(\xi _2,\eta _2\in (-1,1)\) and e defined in (6.4).

Therefore, all that remains to do in order to finish the construction of RT-subsolutions is to find \(\xi _2,\eta _2:(\rho _-,\rho _+)\rightarrow (-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

$$\begin{aligned} Q(\rho ):=(\rho -\rho _-)\sqrt{\rho _-}\xi _2(\rho )-(\rho _+-\rho ) \sqrt{\rho _+}\eta _2(\rho )>0, \end{aligned}$$

one has \(e(x_2,t)=g^2t^2{\tilde{e}}(\rho (x_2,t))\) with

$$\begin{aligned} {\tilde{e}}(\rho ):=\frac{1}{2}\frac{\rho _+\rho _-(\rho _+-\rho _-)^2}{Q(\rho )^2}. \end{aligned}$$

By the transformation \(x_2=\frac{1}{2}gt^2G'(\rho )\) the desired admissibility (2.8) in the strict sense is then equivalent to

$$\begin{aligned}&\int _{\rho _-}^{\rho _+}\left( {\tilde{e}}(\rho )-\frac{1}{2}\rho -G(\rho ) \right) G''(\rho )\,\mathrm{d}\rho \nonumber \\&\quad <\frac{1}{4}\int _{\rho _-}^{\rho _+} \big (\rho _0(G'(\rho ))-\rho \big )\big (G'(\rho )^2\big )'\,\mathrm{d}\rho . \end{aligned}$$
(6.9)

We further make the ansatz \({\tilde{e}}(\rho _\pm )=\frac{1}{2}\rho _\pm \), in other words that e is continuious across \(\partial {\mathscr {U}}\). Then partial integration shows that (6.9) is equivalent to

$$\begin{aligned} I_{\xi _2,\eta _2}:=\int _{\rho _-}^{\rho _+}\left( {\tilde{e}}'(\rho ) -\frac{3}{4}G'(\rho )\right) G'(\rho )\,\mathrm{d}\rho >0. \end{aligned}$$
(6.10)

Observe that the condition \({\tilde{e}}(\rho _\pm )=\frac{1}{2}\rho _\pm \) requires \(\xi _2(\rho _+)=1\), \(\eta _2(\rho _-)=-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 \(|\xi |=|\eta |=1\) throughout all of \([\rho _-,\rho _+]\), that is \(\xi _2\equiv -\eta _2\equiv 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.

Denote by \(Q_0\), \(G_0\), \({{\tilde{e}}}_0\) the functions associated with the choice \(\xi _2\equiv -\eta _2\equiv 1\), that is

$$\begin{aligned} Q_0(\rho )&=(\rho -\rho _-)\sqrt{\rho _-}+(\rho _+-\rho )\sqrt{\rho _+},\quad {\tilde{e}}_0(\rho )=\frac{1}{2}\frac{\rho _+\rho _-(\rho _+-\rho _-)^2}{Q_0(\rho )^2},\\ G_0(\rho )&=\frac{(\rho _+-\rho )(\rho -\rho _-)(\sqrt{\rho _-}-\sqrt{\rho _+})}{Q_0(\rho )} =-\frac{(\rho _+-\rho )(\rho -\rho _-)}{\rho _++\rho _-+\sqrt{\rho _+\rho _-}-\rho }. \end{aligned}$$

Through further calculations, one obtains that

$$\begin{aligned} G_0'(\rho )&=\frac{\sqrt{\rho _-\rho _+}(\rho _++\rho _-) +2\rho _-\rho _+}{(\rho _++\rho _-+\sqrt{\rho _+\rho _-}-\rho )^2}-1,\\ G_0''(\rho )&=\frac{2\sqrt{\rho _-\rho _+}(\rho _++\rho _-) +4\rho _-\rho _+}{(\rho _++\rho _-+\sqrt{\rho _+\rho _-}-\rho )^3}, \end{aligned}$$

in particular one sees that \(G_0\) is uniformly strictly convex on \([\rho _-,\rho _+]\). Furthermore, using the transformation \(\rho =\rho _++\rho _-+\sqrt{\rho _+\rho _-}-\rho _-s\) and the abbreviation \(r:=\frac{\sqrt{\rho _+}}{\sqrt{\rho _-}}\) we obtain

$$\begin{aligned} \frac{I_{1,-1}}{\rho _-}=\int _{1+r}^{r^2+r}\left( \frac{r^2(1+r)^2}{s^3} +\frac{3}{4}-\frac{3r(1+r)^2}{4s^2}\right) \left( \frac{r(1+r)^2}{s^2}-1\right) \,\mathrm{d}s=0. \end{aligned}$$

This means that with the choice \(\xi _2\equiv -\eta _2\equiv 1\) there holds equality in (2.8) for any \(t>0\).

We now turn to the perturbation. Let \(\varepsilon >0\) and consider

$$\begin{aligned} \xi _2(\rho ):=1+\varepsilon {{\bar{\xi }}}(\rho ),\quad \eta _2(\rho ):=-1+\varepsilon {{\bar{\eta }}}(\rho ), \end{aligned}$$
(6.11)

with functions \({{\bar{\xi }}},{{\bar{\eta }}}:[\rho _-,\rho _+]\rightarrow {\mathbb {R}}\) satisfying \({\bar{\xi }}<0\), \({\bar{\eta }}>0\) on \((\rho _-,\rho _+)\) and \({\bar{\xi }}(\rho _\pm )={\bar{\eta }}(\rho _\pm )=0\). Again, the last condition allows the function e defined via (6.4) to be continuous over the whole domain \(\Omega \times (0,T)\).

We will look for asymptotic expansions of the associated \(Q=Q_\varepsilon \), \(G=G_\varepsilon \), \({{\tilde{e}}}=\tilde{e}_\varepsilon \) with respect to \(\varepsilon >0\). It holds that

$$\begin{aligned} Q_\varepsilon (\rho )&=Q_0(\rho )+\varepsilon \left( (\rho -\rho _-) \sqrt{\rho _-}{\bar{\xi }}-(\rho _+-\rho )\sqrt{\rho _+}{{\bar{\eta }}} \right) =:Q_0(\rho )+\varepsilon {\bar{Q}}(\rho ),\\ {{\tilde{e}}}_\varepsilon (\rho )&=\frac{1}{2} \frac{\rho _+\rho _-(\rho _+-\rho _-)^2}{(Q_0(\rho ) +\varepsilon \bar{Q}(\rho ))^2}={{\tilde{e}}}_0(\rho ) -\varepsilon \rho _+\rho _-(\rho _+-\rho _-)^2\frac{{\bar{Q}}(\rho )}{Q_0(\rho )^3} +{\mathcal {O}}(\varepsilon ^2)\\&=:{{\tilde{e}}}_0(\rho )+\varepsilon {\bar{e}}(\rho )+{\mathcal {O}}(\varepsilon ^2),\\ G_\varepsilon (\rho )&=G_0(\rho )+\varepsilon \frac{(\rho _+-\rho ) (\rho -\rho _-)}{Q_0^2(\rho )}\sqrt{\rho _+\rho _-}(\rho _+-\rho _-) ({{\bar{\xi }}}+{{\bar{\eta }}})+{\mathcal {O}}(\varepsilon ^2)\\&=:G_0(\rho )+\varepsilon {\bar{G}}(\rho )+{\mathcal {O}}(\varepsilon ^2), \end{aligned}$$

while the expansion of \(I_\varepsilon :=I_{1+\varepsilon {\bar{\xi }},-1+\varepsilon {\bar{\eta }}}\) reads as

$$\begin{aligned} I_\varepsilon =\varepsilon \int _{\rho _-}^{\rho _+} \left( {\tilde{e}}_0'(\rho ){{\bar{G}}}'(\rho )+{\bar{e}}'(\rho ) G_0'(\rho )-\frac{3}{2} G_0'(\rho ){{\bar{G}}}'(\rho )\right) \, \mathrm{d}\rho +{\mathcal {O}}(\varepsilon ^2)=:\varepsilon {\bar{I}}+{\mathcal {O}}(\varepsilon ^2). \end{aligned}$$

Since \(G_0\) is uniformly convex on \([\rho _-,\rho _+]\), the perturbed function \(G_\varepsilon \) will also be uniformly convex for small enough \(\varepsilon >0\). Moreover, in order to have admissibility for \(\varepsilon >0\) small enough, it suffices to have \({{\bar{I}}}>0\).

By integration by parts we rewrite

$$\begin{aligned} {{\bar{I}}}&=-\int _{\rho _-}^{\rho _+}\left( {\tilde{e}}_0''(\rho )\bar{G}(\rho )+{\bar{e}}(\rho ) G_0''(\rho )-\frac{3}{2} G_0''(\rho )\bar{G}(\rho )\right) \, \mathrm{d}\rho \\ {}&=\int _{\rho _-}^{\rho _+}{{\bar{\xi }}}(\rho ) H_1(\rho ) \, \mathrm{d}\rho +\int _{\rho _-}^{\rho _+}{{\bar{\eta }}}(\rho ) H_2(\rho ) \, \mathrm{d}\rho , \end{aligned}$$

where

$$\begin{aligned} H_1(\rho )&=\frac{(\rho _+-\rho )(\rho -\rho _-)}{Q_0^2(\rho )} \sqrt{\rho _+\rho _-}(\rho _+-\rho _-)\left( \frac{3}{2}G''_0(\rho )-\tilde{e}''_0(\rho )\right) \\&\quad +\frac{\rho _+\rho _-(\rho _+-\rho _-)^2}{Q_0^3(\rho )} \sqrt{\rho _-}(\rho -\rho _-)G''_0(\rho ),\\ H_2(\rho )&=\frac{(\rho _+-\rho )(\rho -\rho _-)}{Q_0^2(\rho )} \sqrt{\rho _+\rho _-}(\rho _+-\rho _-)\left( \frac{3}{2}G''_0(\rho )-\tilde{e}''_0(\rho )\right) \\&\quad -\frac{\rho _+\rho _-(\rho _+-\rho _-)^2}{Q_0^3(\rho )} \sqrt{\rho _+}(\rho _+-\rho )G''_0(\rho ). \end{aligned}$$

It then follows that in order to have \({\bar{I}}>0\), it suffices to find \({{\bar{\rho }}}\in (\rho _-,\rho _+)\) such that either \(H_1({{\bar{\rho }}})<0\) or \(H_2({{\bar{\rho }}})>0\). Indeed, if \(H_1({{\bar{\rho }}})<0\), one may choose a smooth function \(\rho \mapsto {{\bar{\xi }}}(\rho )\) such that it is strictly negative on \((\rho _-,\rho _+)\), vanishes at the endpoints and concentrates at \({{\bar{\rho }}}\) sufficiently such that \(\int _{\rho _-}^{\rho _+}{{\bar{\xi }}}(\rho ) H_1(\rho ) \, \mathrm{d}\rho >0\). Then, regardless of the sign of \(H_2\), one may clearly choose a function \(\rho \mapsto {{\bar{\eta }}}={{\bar{\eta }}}(\rho )\) which is strictly positive on \((\rho _-,\rho _+)\), vanishes at the endpoints, and is small enough such that \({\bar{I}}>0\). The case \(H_2({{\bar{\rho }}})>0\) can be treated similarly.

Finally, to conclude the proof of Theorem 2.7, we will prove that in fact the first case \(H_1({{\bar{\rho }}})<0\) is not possible, while \(H_2({{\bar{\rho }}})>0\) is possible if and only if \(\sqrt{\frac{\rho _+}{\rho _-}}>\frac{4+2\sqrt{10}}{3}\).

Let us first prove the second statement. \(H_2({{\bar{\rho }}})>0\) is equivalent to

$$\begin{aligned} Q_0({{\bar{\rho }}})({{\bar{\rho }}}-\rho _-)\sqrt{\rho _-}\left( \frac{3}{2}G''_0({{\bar{\rho }}})-\tilde{e}''_0({{\bar{\rho }}})\right) -\rho _+\rho _-(\rho _+-\rho _-)G''_0({{\bar{\rho }}})>0. \end{aligned}$$

Plugging in the expressions for \(Q_0,G_0\) and \({\tilde{e}}_0\), one obtains that this is equivalent to

$$\begin{aligned} {{\bar{\rho }}}^2-(\rho _++2\rho _-){{\bar{\rho }}}+\frac{2}{3} \rho _+^{3/2}\rho _-^{1/2}+\frac{5}{3}\rho _+\rho _-+\rho _-^2<0. \end{aligned}$$

This is possible only if the discriminant with respect to \({{\bar{\rho }}}\) is strictly positive, which reads as

$$\begin{aligned} (\rho _++2\rho _-)^2-4\left( \frac{2}{3}\rho _+^{3/2}\rho _-^{1/2} +\frac{5}{3}\rho _+\rho _-+\rho _-^2\right) >0, \end{aligned}$$

or equivalently,

$$\begin{aligned} r^2\left( r^2-\frac{8}{3}r-\frac{8}{3}\right) =r^2\left( r -\frac{4-2\sqrt{10}}{3}\right) \left( r-\frac{4+2\sqrt{10}}{3}\right) >0, \end{aligned}$$

where we have denoted \(r:=\sqrt{\frac{\rho _+}{\rho _-}}>1\). The statement then follows by taking for instance \({\bar{\rho }}=\frac{\rho _++2\rho _-}{2}\in (\rho _-,\rho _+)\) due to \(\sqrt{\frac{\rho _+}{\rho _-}}>\frac{4+2\sqrt{10}}{3}\).

The case \(H_1({{\bar{\rho }}})<0\) being not possible is proven similarly, the same calculations yield the condition \(\frac{1}{r^2}-\frac{8}{3r}-\frac{8}{3}>0\), which is not possible for \(r>1\).

It remains to compute the precise growth rates of the mixing zone \({\mathscr {U}}\) given in (6.8). Observe that \(\partial _\xi G(\rho _\pm )=\partial _\eta G(\rho _\pm )=0\), such that \({\bar{\xi }}(\rho _\pm )={\bar{\eta }}(\rho _\pm )=0\) implies

$$\begin{aligned} G'_\varepsilon (\rho _\pm )&=G'_0(\rho _\pm )=\frac{\sqrt{\rho _\pm } -\sqrt{\rho _\mp }}{\sqrt{\rho _\mp }}. \end{aligned}$$

This concludes the proof of Theorem 2.7. \(\square \)

We would like to point out that the condition \(\sqrt{\frac{\rho _+}{\rho _-}}>\frac{4+2\sqrt{10}}{3}\) only enters in the admissibility of the subsolutions, more precisely it comes from our construction above for assuring \({{\bar{I}}}>0\). For an arbitrary ratio \(\frac{\rho _+}{\rho _-}>1\) the fact that in the unperturbed case \(I_{-1,1}=0\) shows that there exist infinitely many turbulently mixing solutions with the exact same growth rates \(c_\pm (t)\) violating the weak admissibility by an arbitrary small amount of energy.

Furthermore, we summarize the other ansatzes used during our construction and note that they can all be seen as not too restrictive for different reasons:

  • The independence of \(x_1\) can be interpreted as an averaging in the \(x_1\) direction.

  • \(v\equiv 0\) for the subsolution is in harmony with the vanishing initial velocity and the fact that the subsolution corresponds to an averaging of solutions.

  • \(\xi \) and \(\eta \) only depending on \(\rho \) allow us to find the density \(\rho \) as the unique entropy solution of a relatively simple conservation law, this generalizes the construction from [33, 34], where the unique viscosity solution of a Burgers equation was considered. In fact a similar conservation law also appeared in the relaxation of the two-phase porous media flow with different mobilities by Otto [29]. Our intuition behind choosing \(\xi \) and \(\eta \) to be perturbations of \(\pm e_2\) has been explained during the proof. Nonetheless, it would be interesting to see if other choices of \(\xi \) and \(\eta \) also lead to admissible subsolutions.

  • The continuity of e across \(\partial {\mathscr {U}}\) is not a huge jump from Definition 2.6, which combined with \(v\equiv 0\) already implied that \(e=\frac{1}{2}g^2t^2\rho _+\) in \(\{x_2>0\}\cap {\mathscr {D}}{\setminus }\overline{{\mathscr {U}}}\), respectively \(e=\frac{1}{2}g^2t^2\rho _-\) in \(\{x_2<0\}\cap {\mathscr {D}}{\setminus }\overline{{\mathscr {U}}}\), and therefore the continuity of e in each of the three pieces \(\{x_2<0\}\cap {\mathscr {D}}{\setminus }\overline{{\mathscr {U}}}\), \(\mathscr {U}\) and \(\{x_2>0\}\cap {\mathscr {D}}{\setminus }\overline{{\mathscr {U}}}\).

Finally, we would like to state further properties than those of the growth rates of the unperturbed “subsolution” associated with \(\xi _2\equiv 1\), \(\eta _2\equiv -1\) in an explicit way. Inversion of the derivative \(G_0':[\rho _-,\rho _+]\rightarrow \left[ -\frac{\sqrt{\rho _+}-\sqrt{\rho _-}}{\sqrt{\rho _+}}, \frac{\sqrt{\rho _+}-\sqrt{\rho _-}}{\sqrt{\rho _-}}\right] \) shows that the density profile, defined in (6.7), inside the mixing zone is given by

$$\begin{aligned} \rho \left( x_2,t\right) =\rho _++\sqrt{\rho _+\rho _-}+\rho _--\frac{(\sqrt{\rho _+} +\sqrt{\rho _-})\root 4 \of {\rho _+\rho _-}}{\sqrt{1+\frac{2x_2}{gt^2}}}, \end{aligned}$$

the relaxed momentum \(u_2(x_2,t)=gtG_0(\rho (x_2,t))\) and e defined in (6.4) inside \({\mathscr {U}}\) by

$$\begin{aligned} u_2(x_2,t)&=gt(\sqrt{\rho _+}+\sqrt{\rho _-}) \left( \frac{\root 4 \of {\rho _+\rho _-}}{\sqrt{1+\frac{2x_2}{gt^2}}} +\root 4 \of {\rho _+\rho _-}\sqrt{1+\frac{2x_2}{gt^2}}-\sqrt{\rho _+}-\sqrt{\rho _-}\right) \\ e(x_2,t)&=\frac{1}{2}g^2t^2\sqrt{\rho _-\rho _+}\left( 1+\frac{2x_2}{gt^2}\right) , \end{aligned}$$

from which an interested reader can obtain a formula of the associated energy density \(E_{sub}\) defined in (2.7). Here we would only like to state the conversion rate of total potential energy into total kinetic energy. Recall that the unperturbed “subsolution” satisfies (2.8) with equality. Hence the total kinetic energy at time \(t\geqq 0\) can be expressed as the difference in total potential energy, which is

$$\begin{aligned} \int _{\Omega }(\rho _0(x)-\rho (x,t))gx_2\,\mathrm{d}x&=\frac{g^3t^4}{8} \int _{\rho _-}^{\rho _+}(\rho _0(G_0'(\rho ))-\rho )\big (G'_0(\rho )^2\big )'\,\mathrm{d}\rho \\&=\frac{g^3t^4}{8}\int _{\rho _-}^{\rho _+}G'_0(\rho )^2\,\mathrm{d}\rho \\&=\frac{g^3t^4(\sqrt{\rho _+}+\sqrt{\rho _-})(\sqrt{\rho _+} -\sqrt{\rho _-})^3}{24\sqrt{\rho _+\rho _-}}. \end{aligned}$$

We conclude the paper by presenting a plot of the above density (blue) and momentum (red) profiles for the choice \(\rho _-=1/4\), \(\rho _+=4\), \(g=1\) at fixed time \(t=\left( \frac{2}{g(\sqrt{\rho _+}-\sqrt{\rho _-})}\right) ^{\frac{1}{2}}\). At this specific time the mixing zone extends from \(x_2=-\rho _+^{-1/2}=-1/2\) to \(x_2=\rho _-^{-1/2}=2\).

figure a