Digital Object Identifier (DOI) https://doi.org/10.1007/s00205-021-01672-1
A New Approach to the Rayleigh–Taylor Instability
Björn Gebhard , József J. Kolumbán & László Székelyhidi Jr.
Communicated byA. Bressan
Abstract
In this article we consider the inhomogeneous incompressible Euler equations describing two fluids with different constant densities under the influence of gravity as a differential inclusion. By considering the relaxation of the constitutive laws we formulate a general criterion for the existence of infinitely many weak solutions which reflect the turbulent mixing of the two fluids. Our criterion can be verified in the case that initially the fluids are at rest and separated by a flat interface with the heavier one being above the lighter one—the classical configuration giving rise to the Rayleigh–Taylor instability. We construct specific examples when the Atwood number is in the ultra high range, for which the zone in which the mixing occurs grows quadratically in time.
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].
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 724298-DIFFINCL).
The mathematical model (see for example Section 6.4 of [26]) is given by the inhomogeneous incompressible Euler equations
∂t(ρv)+div(ρv⊗v)+ ∇p= −ρgen, divv=0,
∂tρ+div(ρv)=0,
(1.1) which we consider on a bounded domain ⊂ Rn,n 2 and a time interval [0,T). Hereρ :× [0,T)→Rdenotes the fluid density,v :× [0,T)→Rn is the velocity field, respectively p : × [0,T)→ Ris the pressure,g > 0 is the gravitational constant andenis the n-th Euclidean coordinate vector. Compared to the homogenous density case,ρ ≡ 1, the solvability of the Cauchy problem of (1.1) for a general non-constant initial density distribution is more delicate even in the planar case; see Section 6.4 of [26]. Results concerning the local well- posedness have only been obtained under sufficiently strong regularity assumptions on the initial density; see [13–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–15,37].
More precisely, we consider (1.1) together with initial datav0 : → Rn, ρ0:→Rsatisfying
divv0=0 andρ0∈ {ρ−, ρ+} almost everywhere (1.2) with two fixed valuesρ+> ρ−>0. In fact our main focus lies on the flat unstable initial configuration
v0≡0 andρ0(x)=
ρ+whenxn>0,
ρ−whenxn0, (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 homo- geneous 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 solu- tions obtained this way as turbulent mixing. More precisely, we prove the existence of solutions with the following properties:
Forρ+ > ρ− >0 define the Atwood numberA= ρρ++−ρ+ρ−− and the quadratic functionsc±:R→R,
c+(t)= ρ++ρ−
2√ρ−(√ρ++ √ρ−)Agt2, c−(t)= ρ++ρ−
2√ρ+(√ρ++ √ρ−)Agt2. LetT >0 and=(0,1)×(−c−(T),c+(T))⊂R2.
Theorem 1.1.Let ρρ+
−
4+2√ 10 3
2
. The initial value problem(1.1),(1.3)has infinitely many weak admissible solutions(ρ, v)∈L∞(×(0,T);R×R2)with ρ ∈ {ρ−, ρ+}almost everywhere and such that
(i)ρ(x,t)=ρ+,v=0for x2c+(t), (ii)ρ(x,t)=ρ−,v=0for x2−c−(t),
(iii)for any open ball B contained in{(x,t)∈×(0,T):x2∈(−c−(t),c+(t))} it holds that
Bρ+−ρ(x,t)dxdt·
Bρ(x,t)−ρ−dxdt >0.
For the precise definition of weak admissible solutions we refer to Definitions2.1, 2.2.
We would like to point out that the infinitely many weak solutions differ only in their turbulent fine structure, while they all have a continuous coarse grained density profileρ¯ in common. The profile ρ(x,¯ t) = ¯ρ(x2,t)can be seen as an x1-average of the solutions, cf. Remark2.5(a) below and [7, Remark 1.2], and is found as the entropy solution of a conservation law
∂tρ¯+gt∂x2G(ρ)¯ =0,
which up to the factor t shows similarities to the conservation law appearing in Otto’s relaxation for the incompressible porous media equation [29]. Further de- tails and the explicit profileρ¯can be found in Section6.
The condition that the density ratioρ+/ρ−is larger than 4+2√
10 3
2
≈11.845, implies that the Atwood number is in the so-called (for example [5]) “ultra high”
range (0.845,1). The main obstruction in establishing a similar result to Theo- rem 1.1for a density ratio outside of this range comes from the fact that in our approach the weak admissibility condition on the solutions associated with the profileρ¯reduces to an algebraic inequality for the density ratio; see Section6for 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 Section2.
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αAgt2, but there is considerable disagreement about the value of the constantαand its possible dependence onA(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 ratesc±in Theorem1.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 insta- bility 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 C2and 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 Section2we 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 Section3we 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 dif- ferential inclusion. The differential inclusion fits in a modified version of the Tartar framework of convex integration, adapted from [17,33] to simultaneously handle the absence of the pressure from the set of constraints and the dependence of the set of constraints on(x,t)due to the prescribed energy density function.
In Section4we prove the ingredients of the topological framework, most im- portantly we calculate the-convex hull of the set of constraints, which forms the core of this paper.
In Section5we conclude the proof of our main convex integration result, while in Section6we construct appropriate subsolutions having the growth rates presented in Theorem1.1.
2. Statement of Results
Let⊂Rnbe a bounded domain andT >0. Our notion of solution to equation (1.1) on× [0,T)is as follows:
Definition 2.1.(Weak solutions) Let (ρ0, v0) ∈ L∞()×L2(;Rn)such that (1.2) holds almost everywhere in . We say that (ρ, v) ∈ L∞(×(0,T))× L2(×(0,T);Rn)is a weak solution to equation (1.1) with initial data(ρ0, v0) if for any test functions∈Cc∞(× [0,T);Rn),1∈Cc∞(× [0,T)),2 ∈ Cc∞(× [0,T))such thatis divergence-free, we have
T 0
[ρv·∂t+ ρv⊗v,∇ −gρn] dxdt +
ρ0(x)v0(x)·(x,0)dx=0, T
0
v· ∇1dxdt =0, T
0
[ρ∂t2+ρv· ∇2] dxdt+
ρ0(x)2(x,0)dx=0, and ifρ(x,t)∈ {ρ−, ρ+}for almost every(x,t)∈×(0,T).
Note that the definition ofvbeing 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 withv, 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 pis determined up to a function depending only on time, as in the case of the homogeneous Euler equations, see [36].
As in the homogeneous case, one can associate with a weak solution(ρ, v)an energy density functionE ∈L1(×(0,T))given by
E(x,t):= 1
2ρ(x,t)|v(x,t)|2+ρ(x,t)gxn. Furthermore, for smooth solutions of (1.1) one can show thatt →
E(x,t)dx is constant. For weak solutions this necessarily does not need to be true. As in the case of the homogeneous Euler equations or hyperbolic conservation laws, in order to not investigate physically irrelevant solutions we require our weak solutions to be admissible with respect to the initial energy.
Definition 2.2.(Admissible weak solutions) A weak solution(ρ, v)in the sense of Definition2.1is called admissible provided it satisfies the weak energy inequality
E(x,t)dx
E(x,0)dxfor almost everyt∈(0,T).
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
∂tu+divS+ ∇P= −ρgen,
∂tρ+divu=0, divv=0,
(2.1)
considered on×(0,T)and withz=(ρ, v,u,S,P)taking values in the space Z =R×Rn×Rn×S0n×n×R. HereS0n×ndenotes the space of symmetricn×n matrices with trace 0. We will also writeSn×nfor the space of symmetric matrices, id ∈Sn×nfor the identity andλmax(S), λmin(S)for the maximal, minimal resp., eigenvalue ofS∈Sn×n.
As usual, equations (2.1) will be complemented by a set of pointwise constraints.
Lete:×(0,T)→R+be a given function and define for(x,t)∈×(0,T) the sets
K(x,t) := {z∈ Z :ρ∈ {ρ−, ρ+}, u =ρv, ρv⊗v−S=
e(x,t)−2
nρgt en·v− 1 nρg2t2
id , (2.2) as well as the setsU(x,t)by requiring forz∈U(x,t)the following four inequalities to hold:
ρ−< ρ < ρ+, ρ+
n
|u−ρ−v+(ρ−ρ−)gt en|2
(ρ−ρ−)2 <e(x,t),
(2.3) ρ−
n
|u−ρ+v+(ρ−ρ+)gt en|2
(ρ−ρ+)2 <e(x,t), λmax(A(z)) <e(x,t)−2
ngt en·u−1
nρg2t2, (2.4) where
A(z)= ρρ−ρ+v⊗v−ρ−ρ+(u⊗v+v⊗u)+(ρ++ρ−−ρ)u⊗u
(ρ+−ρ)(ρ−ρ−) −S.
Note that by the definition ofK(x,t)in (2.2) and by recalling thatShas vanishing trace, every solution of (2.1) taking values inK(x,t)almost everywhere is a solution
to the inhomogeneous Euler equations (1.1) withρ ∈ {ρ−, ρ+}and associated energy
E =1
2ρ|v|2+ρgxn= n
2e(x,t)−ρgt en·v−1
2ρg2t2+ρgxn, (2.5) which is equivalent to saying that
1
2ρ|v+gt en|2= n 2e(x,t).
Conversely, if we have a solution (ρ, v,p)of (1.1) withρ ∈ {ρ−, ρ+}almost everywhere, we can introduce the variables u = ρv, S = ρv ⊗v−n1ρ|v|2id, P = p+ 1nρ|v|2to see thatz = (ρ, v,u,S,P)will satisfy system (2.1) while pointwise taking valuesz(x,t)∈K(x,t), whereK(x,t)is defined with respect to the function
e(x,t):= 1
nρ(x,t)|v(x,t)+gt en|2.
Since the pressure Pdoes not play a role in the set of constraintsK(x,t), it is convenient to consider the following projection: forz =(ρ, v,u,S,P) ∈ Z we denote
π(z)=(ρ, v,u,S)∈R×Rn×Rn×S0n×n. (2.6) Using the linear system (2.1) and the definition of U(x,t) we define relaxed solutions to (1.1) in the following way:
Definition 2.3.(Subsolutions) Lete:×[0,T)→R+be a bounded function. We say thatz=(ρ, v,u,S,P):×(0,T)→Z is a subsolution of (1.1) associated with eand the initial data (ρ0, v0) ∈ L∞()×L2(;Rn)satisfying (1.2) iff π(z)∈L∞(×(0,T);π(Z)),P is a distribution,zsolves (2.1) in the sense that vis weakly divergence-free (including the weak no-flux boundary condition),
T 0
[u·∂t+ S,∇ −gρn] dxdt+
ρ0(x)v0(x)·(x,0)dx=0, T
0
[ρ∂t+u· ∇] dxdt+
ρ0(x)(x,0)dx =0,
for any test functions∈Cc∞(× [0,T);Rn), div=0, ∈Cc∞(× [0,T)), and if there exists an open set U ⊂ ×(0,T), such that the maps(x,t) → π(z(x,t))and(x,t)→e(x,t)are continuous onU with
z(x,t)∈U(x,t), for all(x,t)∈U,
z(x,t)∈K(x,t), for almost every(x,t)∈×(0,T)\U.
We callU the mixing zone of z. Moreover, the subsolution is called admissible provided that
Esub(x,t):= n
2e(x,t)−gt en·u(x,t)−1
2ρ(x,t)g2t2+ρ(x,t)gxn (2.7) satisfies
Esub(x,t)dx
Esub(x,0)dxfor almost everyt ∈(0,T). (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 : × [0,T) → R+be bounded. If there exists a subsolution z associated with e in the sense of Definition2.3, then for the same initial data of the subsolution there exist infinitely many weak solutions in the sense of Definition2.1, which coincide almost everywhere on×(0,T)\U with z and whose total energy is given by E defined in (2.5). The solutions are turbulently mixing onU in the sense that for any open ball B⊂U it holds that
B
ρ+−ρ(x,t)d(x,t)·
B
ρ(x,t)−ρ−d(x,t) >0. (2.9) Among these weak solutions there exists a sequence{zk}k0such thatρk ρin L2(U). If in additionπ(z)is inC0([0,T];L2(;π(Z)))and satisfies(2.8)with strict inequality for every t ∈ (0,T], then infinitely many of the induced weak solutions are admissible in the sense of Definition2.2.
Remark 2.5.(a) The second to last two statements justify to callU the mixing zone and to interpret the subsolution densityρas a kind of coarse-grained or averaged density profile.
(b) The result carries over to the three- or higher-dimensional case by constructing suitable potentials analoguosly to [16], which is not done here, cf. Lemma4.1.
The other parts of the proof, for example the computation of the-convex hull in Section4.2, are carried out in arbitrary dimensions.
(c) We will see later that the open setU(x,t) is indeed the convex hull ofK(x,t). In particular we can conclude that weak limits of solutions are subsolutions in the following sense: Let(ρk, vk)k∈N be a sequence of essentially bounded weak solutions of (1.1) and define as beforeuk:=ρkvk,Sk:=ρkvk⊗vk−1nρk|vk|2id.
Assume thatzk :=(ρk, vk,uk,Sk) (ρ, v,∗ u,S)=:zinL∞(×(0,T);R× Rn×Rn×S0n×n). Assume further that there exists a continuous bounded function e∈ C0(×(0,T)), such thatek := n1ρk|vk+gt en|2→ einL∞(×(0,T)).
Thenzsupplemented by a possibly distributionalPis a weak solution of the linear system2.1with(z(x,t),P(x,t))∈ U(x,t)for almost every(x,t)∈×(0,T), whereU(x,t)is defined with respect to the functione.
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∩(Rn−1× {0})= ∅.
Definition 2.6.(Rayleigh–Taylor subsolution) We call a subsolution zof (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 everyt∈(0,T).
Theorem 2.7.Let n=2,=(0,1)×(−c−(T),c+(T)), where c−(t)= 1
2
1− ρ−
ρ+
gt2, c+(t)= 1 2
ρ+
ρ− −1
gt2.
Ifρ+>
4+2√ 10 3
2
ρ−, then there exists a RT-subsolution z which only depends on
x2
t2, and at time t >0the mixing zoneU(t):= {x∈:(x,t)∈U }associated with z is(0,1)×(−c−(t),c+(t)).
An explicit description of the subsolutions and further discussion can be found after the proof of Theorem2.7in Section6. Observe that by combining Theorems2.4 and2.7we arrive at the statement of Theorem1.1.
3. Reformulation as a Differential Inclusion
The proof of Theorem2.4will 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⊂Rnbe a bounded domain,T >0 and set D=
(y,t)∈Rn×(0,T):y−1
2gt2en∈ , such that fort∈(0,T)the slice is given by
D(t):=
y∈Rn:(y,t)∈D
=+1 2gt2en. Let(μ, w,q)be a weak solution of
∂t(μw)+div(μw⊗w)+ ∇q =0, divw=0,
∂tμ+div(μw)=0
(3.1)
onDfor some suitable initial data satisfying (1.2) and with weak boundary condi- tion
(w(y,t)−gt en)·νD(t)(y)=0 (3.2) fory∈∂D(t). More precisely, the notion of weak solution to (3.1), (3.2) is under- stood as in Definition2.1, except that nowg=0, in the momentum and continuity equation×(0,T),×[0,T)is replaced byD,D∪(×{0})resp., and the weak formulation of divw=0 including the weak boundary condition (3.2) becomes
Dw· ∇d(y,t)− T
0
∂D(t)(y,t)gt en·νD(t)(y)dS(y)dt
=0 for all ∈C∞(D). (3.3)
Then if we definey:=x+12gt2enand set ρ(x,t)=μ (y,t) ,
v(x,t)=w (y,t)−gt en, p(x,t)=q(y,t) ,
(3.4)
it is straightforward to check that(ρ, v)is a weak solution of (1.1) on×(0,T) with the same initial data(ρ0, v0)=(μ0, w0). 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
1
2μ(y,t)|w(y,t)|2. Let us write 1
2μ(y,t)|w(y,t)|2= n
2e(y−1/2gt2en,t)
for a functione:× [0,T)→R+. Then the total energyE(x,t)associated with the original system (1.1) is precisely given by (2.5).
We can now reformulate (3.1) as a differential inclusion by considering onD the system
∂tm+divσ+ ∇q=0, divw=0,
∂tμ+divm=0,
(3.5)
wherez:=(μ, w,m, σ,q)takes values inZ =R×Rn×Rn×S0n×n×R, together with the set of pointwise constraints
K(y,t)=
z∈Z :μ∈ {μ−, μ+}, m=μw, μw⊗w−σ=e
y−1 2gt2en,t
id ,
(3.6) where in analogy to the homogeneous Euler equationse : ×(0,T) → R+
is given and for the sake of consistency we have denoted μ± := ρ±. We will understand weak solutions of (3.5) in the following sense:
Definition 3.1.We say thatz:D→ Zis a weak solution of (3.5) with initial data π(z0)∈L2(;π(Z))iffπ(z)∈ L2(D;π(Z)),qis a distribution,wsatisfies (3.3) and one has
D[m·∂t+ σ,∇] dxdt+
μ0(x)w0(x)·(x,0)dx=0,
D[μ∂t+m· ∇] dxdt+
μ0(x)(x,0)dx=0,
for any∈C∞c (D∪(× {0});Rn), div=0, ∈Cc∞(D∪(× {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(ρ0, v0)∈ L∞()×L2(;Rn)be initial data satisfying(1.2), e ∈ L1(×(0,T);R+)be a prescribed function. If z = (μ, w,m, σ,q)is a weak solution of (3.5)in the sense of Definition3.1with initial dataμ(·,0)=ρ0, w(·,0) = v0 and if z(y,t) ∈ K(y,t) for almost every(y,t) ∈ D, then the pair (ρ, v)defined by(3.4)is a weak solution of (1.1)on×(0,T)with initial data (ρ0, v0). Moreover, the (possibly distributional) pressure is given by
p(x,t):=q(y,t)−1
nμ(y,t)|w(y,t)|2, y=x+1 2gt2en, 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 ⊂ Z such that for any ¯z ∈ , one can construct a localized plane wave associated with (3.5) oscillating in the direction ofz;¯
• calculating the-convex hull of K(x,t) (denoted by K(x,t)) and proving that one can perturb any element in its interior along sufficiently long-segments, provided that one is far enough fromK(x,t);
• deducing an appropriate set of subsolutions usingK(x,t)and proving that it is a bounded, nonempty subset ofL2(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 Theorem2.4 in Section5by 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):
=
⎧⎨
⎩z¯∈Z :ker
⎛
⎝σ¯ + ¯qid m¯
¯
mT μ¯
¯
wT 0
⎞
⎠= {0}, (μ,¯ m)¯ =0
⎫⎬
⎭.
It has the property that for z¯ ∈ there exists η ∈ Rn+1\{0} such that every z(x,t)= ¯zh((x,t)·η),h ∈C1(R)is a solution of (3.5). In Lemma4.1below we localize these solutions by constructing suitable potentials. Note that the condition (μ,¯ m¯)=0 serves to eliminate the degenerate case when the firstncomponents of ηvanish, that is when one is only allowed to oscillate in time.
Recall the projection operatorπ defined in (2.6).
Lemma 4.1.There exists C>0such that for any¯z∈, there exists a sequence zN∈Cc∞(B1(0);Z)
solving(3.5)and satisfying that (i)d(zN,[−¯z,z¯])→0uniformly, (ii)zN 0in L2,
(iii) |π(zN)|2dxdt C|π(¯z)|2.
Proof. We will only present the proof in the two-dimensional case, higher dimen- sions can be handled analogously to [16].
We start by observing that for any smooth functions ψ : R2+1 → R, φ : R2+1→S2×2, settingD(φ, ψ)=(μ, w,m, σ,q)with
μ=div divφ, w= ∇⊥ψ, m= −∂tdivφ, q = 1
2tr∂t tφ, σ =∂t tφ−qid, implies thatD(φ, ψ)solves (3.5).
LetS :R→ Rbe a smooth function,N 1 andz¯∈ with(μ,¯ m)¯ =0. It follows from the definition ofthat there exists
0=(ξ,c)∈ker
⎛
⎝σ¯ + ¯qid m¯
¯
mT μ¯
¯
wT 0
⎞
⎠. (4.1)
We then treat two cases.
Case 1:c=0
Note that in this case we also haveξ =0, sinceξ =0 would imply(μ,¯ m)¯ =0.
We then set
φN(x,t)= 1
c2(σ¯ + ¯qid) 1
N2S(N(ξ,c)·(x,t)), ψN(x,t)= | ¯w|sgn(ξ⊥· ¯w)
|ξ| 1
NS(N(ξ,c)·(x,t)), and we claim that
D(φN, ψN)= ¯z S(N(ξ,c)·(x,t)). (4.2)
Indeed, using (4.1), one has div divφN = 1
c2ξT(σ¯ + ¯qid)ξS(N(ξ,c)·(x,t))
= 1
c2ξT(−cm)S¯ (N(ξ,c)·(x,t))= ¯μS(N(ξ,c)·(x,t)),
∂tdivφN = 1
c(σ¯ + ¯qid)ξS(N(ξ,c)·(x,t))= − ¯m S(N(ξ,c)·(x,t)),
∂t tφN =c21
c2(σ¯ + ¯qid)S(N(ξ,c)·(x,t))=(¯σ + ¯qid)S(N(ξ,c)·(x,t)),
∇⊥ψN =ξ⊥| ¯w|sgn(ξ⊥· ¯w)
|ξ| S(N(ξ,c)·(x,t))= ¯wS(N(ξ,c)·(x,t)).
From here on, the localization is done in the standard fashion (for example as in [11,16]). We fix S(·) = −cos(·)and, forε > 0, consider χε ∈ Cc∞(B1(0)) satisfying |χε| 1 on B1(0),χε = 1 on B1−ε(0). It is then straightforward to check thatzN=D(χε(φN, ψ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ξ =0, so we may also assume without loss of generality that|ξ| =1. On the other hand, (4.1) implies that there exist constantsk1,k2,k3∈Rsuch that
¯
w=k1ξ⊥, m¯ =k2ξ⊥, σ¯ + ¯qid=k3ξ⊥⊗ξ⊥. (4.3) We set
φN(x,t)= ¯μid 1
N2S(Nξ ·x), ψN(x,t)= | ¯w|sgn(ξ⊥· ¯w)
|ξ| 1
NS(Nξ·x), from where with similar calculations as in Case 1, we obtain that
D(φN, ψN)=(μ,¯ w,¯ 0,0,0)S(Nξ ·x). (4.4) To handle the remaining terms(m,¯ σ ,¯ q), we introduce a different type of po-¯ tential, 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ω: R2+1→R2+1, definingW =curl(x,t)ωandD˜(ω)=(0,0,m, σ,q)by
m= −1
2∇⊥W3, σ+qid=
∂2W1 1
2(∂2W2−∂1W1)
1
2(∂2W2−∂1W1) −∂1W2
implies thatD(ω)˜ solves (3.5).
Now, if we considerωof the form ωN(x)=(a,b,a) 1
N2S(Nξ ·x),
for some constantsa,b∈R, withSas before, we obtain that ∂2W1 1
2(∂2W2−∂1W1)
1
2(∂2W2−∂1W1) −∂1W2
=aξ⊥⊗ξ⊥S(Nξ·x),
∇⊥W3=(ξ1b−ξ2a)ξ⊥S(Nξ·x).
Ifξ1=0, it follows from (4.3) that settinga=k3,b= −2k2ξ+1k3ξ2 gives us D(ω˜ N)=(0,0,m,¯ σ ,¯ q)S¯ (Nξ·x).
from where, using (4.4), we get
D(φN, ψN)+ ˜D(ωN)= ¯z S(Nξ·x).
The localization is then done as in Case 1, by consideringzN =D(χε(φN, ψN))+ D(χ˜ εωN).
Ifξ1=0, then choosinga =k3gives us that D˜(ωN)=
0,0,k3
2ξ2ξ⊥,σ ,¯ q¯
S(Nξ·x).
However, it is easy to see that for any smooth functionθ :R2+1 → R,D(θ)ˆ = (0,0,∇⊥θ,0,0)also solves (3.5). Therefore, we may consider the potential given by
θN(x)=
k2−ξ2
k3
2 1
NS(Nξ·x), we obtain that
∇⊥θN(x)=
k2−ξ2
k3
2
ξ⊥S(Nξ·x), and using (4.3), we get that
D(φN, ψN)+ ˜D(ωN)+ ˆD(θN)= ¯z S(Nξ·x).
One may then localize this potential by the usual means in order to conclude the
proof of the lemma.
4.2. The-Convex Hull
We now turn to the set of pointwise constraintsK(x,t),(x,t)∈ D defined in (3.6). The-convex hullK(x,t)is defined by saying thatz ∈ K(x,t)iff for all- convex functions f : Z → Rthere holds f(z)supz∈K(x,t) f(z), see [23] for more details. In our case it turns out that the-convex hull is nothing else but the usual convex hull, see Proposition4.2below.
For the computation of the hull we drop the(x,t)dependence of the setsK(x,t)
and consider a general set of pointwise constraints given by
K = {z=(μ, w,m, σ,q)∈ Z :μ∈ {μ−, μ+}, m=μw, μw⊗w−σ =eid}, (4.5) where 0< μ−< μ+,e∈R+are given constants.
DefineZ0:= {z∈ Z :μ∈(μ−, μ+)}andT+,T−,Q:Z0→R,M :Z0→ Sn×n,
M(z)= μμ−μ+w⊗w−μ−μ+(m⊗w+w⊗m)+(μ++μ−−μ)m⊗m
(μ+−μ)(μ−μ−) −σ,
Q(z)=λmax(M(z)), T±(z)= μ±
n
|m−μ∓w|2 (μ−μ∓)2 , as well as the open set
U = {z∈ Z :μ∈(μ−, μ+), T+(z) <e, T−(z) <e, Q(z) <e}. (4.6) Proposition 4.2.The-convex hull of K coincides with the convex hull of K and is given by U , that is, K=Kco=U .
Lemma4.4below shows that the closure ofUcan be written as U =K− ∪U0∪K+,
where U0=
z∈Z :μ∈(μ−, μ+), T+(z)e, T−(z)e, Q(z)e , K± = {z∈Z :μ=μ±, m=μ±w, λmax(μ±w⊗w−σ )e}.
Moreover, Lemma 4.8actually shows that K+ ,K− resp., is nothing but the- convex hull ofK+:=K∩ {μ=μ+},K−:=K ∩ {μ=μ−}resp..
Furthermore, notice that if one letsμ+−μ−→0, one recovers fromUexactly the convex hull of the constraints for the homogeneous Euler equations, cf. [17].
The proof of Proposition4.2relies on Lemmas4.4and4.8.
Lemma 4.3.The function Q is convex.
Proof. We write
Q(z)= sup
ξ∈Sn−1
ξTM(z)ξ = sup
ξ∈Sn−1
gξ(z)−ξTσ ξ ,
where for every fixedξ ∈Sn−1the functiongξ :Z0→Ris given by gξ(z)=ξTM(z)ξ+ξTσξ
= μμ−μ+(w·ξ)2−2μ−μ+(m·ξ)(w·ξ)+(μ++μ−−μ)(m·ξ)2
(μ+−μ)(μ−μ−) .
We will show that everygξis convex. As a consequenceQis convex as a supremum of convex functions. In order to do this let us complementξ ∈Sn−1to a orthonormal basis(ξ, v2, . . . , vn)ofRn. Expressingwandmwith respect to this basis one sees that it is enough to show that the functiong:(μ−, μ+)×R2→R,
g(μ,x)= μμ−μ+x12−2μ−μ+x1x2+(μ++μ−−μ)x22 (μ+−μ)(μ−μ−)
is convex. We writeg(μ,x)=xTA(μ)xwith
A(μ):= 1
(μ+−μ)(μ−μ−)
μμ−μ+ −μ−μ+
−μ−μ+ μ++μ−−μ
.
Let us fix (μ,x) ∈ (μ−, μ+)×R2 and observe that A(μ) is positive definite becauseμμ−μ+>0 and
det[(μ+−μ)(μ−μ−)A(μ)] =μ−μ+(μ+−μ)(μ−μ−) >0.
Thus the restricted functiong(μ,·)is convex, or equivalentlyD2g(μ,x)[0,y]20 for all y ∈ R2. It therefore remains to show that D2g(μ,x)[1,y]2 0 for all y∈R2. By the positive definiteness ofA(μ)we obtain
D2g(μ,x)[1,y]2=xTA(μ)x+4yTA(μ)x+2yTA(μ)y
=2
y+A(μ)−1A(μ)xT
A(μ)
y+A(μ)−1A(μ)x +xTA(μ)x−2xTA(μ)A(μ)−1A(μ)x
xT
A(μ)−2A(μ)A(μ)−1A(μ) x.
Now we claim that in factA(μ)=2A(μ)A(μ)−1A(μ), which finishes the proof.
Indeed, differentiation of the identity (μ+−μ)(μ−μ−)A(μ)=
μμ−μ+ −μ−μ+
−μ−μ+μ++μ−−μ
shows that
(μ+−μ)(μ−μ−)A(μ)=(2μ−μ−−μ+)A(μ)+C, (4.7) (μ+−μ)2(μ−μ−)2A(μ)=2((μ+−μ)(μ−μ−)+(2μ−μ−−μ+)2)A(μ) +2(2μ−μ−−μ+)C, (4.8) where
C :=
μ−μ+ 0
0 −1
. Moreover, in a straightforward way one can check that
(μ+−μ)(μ−μ−)
A(μ)C−12
+(μ−+μ+−2μ)A(μ)C−1=idR2,
which implies that
C A(μ)−1C =(μ+−μ)(μ−μ−)A(μ)+(μ−+μ+−2μ)C. (4.9) Now (4.7)–(4.9) imply the identityA(μ)=2A(μ)A(μ)−1A(μ). Lemma 4.4.The set U is convex and U =K− ∪U0∪K+. In particular K ⊂U . Proof. Forμ ∈ (μ−, μ+) the two conditions T+(z) < e, T−(z) < ecan be rewritten as
|m−μ−w|<c+(μ−μ−),
|m−μ+w|<c−(μ+−μ), (4.10) wherec±=
μne±
1/2
.Using the basic triangle inequality one can check that the two conditions in (4.10) define a convex set. By Lemma4.3we already know that Qis a convex function. Hence we have shown thatUis convex.
Now we turn to the characterization ofU. ClearlyU0 ⊂ U. Let us show that K+ ⊂ U. The inclusion K− ⊂ U can be obtained in the same way. Let z∗∈K+ . Take anyz∈ Kwithμ=μ−and some sequence(μj)j∈N⊂(μ−, μ+) converging toμ+. Define
zj = μ+−μj
μ+−μ−z+ μj −μ− μ+−μ−z∗.
Clearlyzj →z∗as j → ∞. Sincez∗∈K+ andz∈ K−a short calculation shows T+(zj)= μ+
n |w∗|2= 1
ntr(μ+w∗⊗w∗−σ∗)λmax(μ+w∗⊗w∗−σ∗)e. Similarly we obtainT−(zj)=e. In a third, slightly longer computation we plugzj
into M, sort with respect to the termsw⊗w,w⊗w∗,w∗⊗w,w∗⊗w∗and find
M(zj)= μ+−μj
μ+−μ−
μ−w⊗w−σ
+ μj−μ− μ+−μ−
μ+w∗⊗w∗−σ∗
= μ+−μj
μ+−μ−eid+μj−μ− μ+−μ−
μ+w∗⊗w∗−σ∗ .
We conclude Q(zj)=λmax(M(zj)) e. Hence everyzj and therefore also the limitz∗is contained inU. So far we knowK− ∪U0∪K+ ⊂U.
For the other inclusion consider(zj)j∈N⊂U,zj →z∗. The interesting case of course isμ∗∈/(μ−, μ+), sayμ∗=μ+. By (4.10) we directly see thatm∗=μ+w∗. Moreover, rewriting
M(z)=μ−m−μw
μ−μ− ⊗m−μ+w
μ+−μ +m−μ−w
μ−μ− ⊗m−σ and looking at (4.10) yields
jlim→∞M(zj)=μ+w∗⊗w∗−σ∗.
Thusλmax(M(zj)) <e, j ∈Nimpliesz∗∈ K+. The caseμ∗=μ−can again be treated by obvious adaptations. ConsequentlyU=K− ∪U0∪K+ .
