2.1.1 Scalar conservation laws

The simplest examples of conservation laws are scalar conservation laws, that is, $N=1$ in (2.1). In this case, every convex function ${\it\eta}$ gives rise to an entropy pair with corresponding entropy flux $q(u):=\int ^{u}{\it\eta}^{\prime }({\it\xi})f^{\prime }({\it\xi})\,\text{d}{\it\xi}$. This rich family of entropy pairs was used by Kružkov (Reference Kružkov1970) to obtain the following existence, uniqueness and stability of solutions for scalar conservation laws.

Theorem 2.4 (Kružkov Reference Kružkov1970).

Let $u_{0}\in L^{1}(\mathbb{R}^{d})\cap L^{\infty }(\mathbb{R}^{d})$. Then there exists a unique entropy solution $u\in L^{\infty }(\mathbb{R}^{d}\times \mathbb{R}_{+})\cap C(\mathbb{R}_{+},~L^{1}(\mathbb{R}^{d}))$ of the scalar conservation law (2.1). Furthermore, if $u(\cdot ,t)=\mathbf{S}_{t}u_{0}$ denotes this unique solution, then the corresponding solution solution operator

$$\begin{eqnarray}\mathbf{S}_{t}:L^{1}(\mathbb{R}^{d})\cap L^{\infty }(\mathbb{R}^{d})\rightarrow L^{1}(\mathbb{R}^{d})\cap L^{\infty }(\mathbb{R}^{d})\end{eqnarray}$$

satisfies

(2.4)$$\begin{eqnarray}\Vert \mathbf{S}_{t}u_{0}-\mathbf{S}_{t}v_{0}\Vert _{L^{1}(\mathbb{R}^{d})}\leq \Vert u_{0}-v_{0}\Vert _{L^{1}(\mathbb{R}^{d})}\quad \text{for all}~t>0.\end{eqnarray}$$

Thus, the entropy solution operator for a scalar conservation law forms an $L^{1}$ contractive semi-group.

The existence of entropy solutions for scalar conservation laws can be established as the limit of vanishing viscosity approximations or approximations generated by monotone finite volume schemes (Godlewski and Raviart Reference Godlewski and Raviart1991). The key aspect in this context is propagation of compactness: if $\{u^{{\it\varepsilon}}\}$ is a family of uniformly bounded approximate solutions with a limit $\lim _{n\rightarrow \infty }u^{{\it\varepsilon}_{n}}=u_{\infty }$, then (2.4) implies their bounded variation compactness, and one can pass to the strong limit in (2.2) and (2.3), $\lim _{n\rightarrow \infty }f(u^{{\it\varepsilon}_{n}})=f(u_{\infty })$, to recover that the limit $u_{\infty }$ is the desired entropy solution.

2.1.2 One-dimensional systems

Next, we consider one-dimensional systems of conservation laws (2.1) with $d=1$,

(2.5)$$\begin{eqnarray}\begin{array}{@{}rcll@{}}\partial _{t}u+\partial _{x}f(u)\; & =\; & 0,\qquad & (x,t)\in \mathbb{R}\times \mathbb{R}_{+},\\ u(x,0)\; & =\; & u_{0}(x)\qquad & x\in \mathbb{R}.\end{array}\end{eqnarray}$$

A key role is played by the flux Jacobian $A(u):=\partial _{u}f(u)$. The system of conservation laws (2.5) is called hyperbolic if the eigenvalues ${\it\lambda}_{1},\ldots ,{\it\lambda}_{N}$ of the Jacobian $A$ are real and the corresponding eigenvectors $r_{1},\ldots ,r_{N}\in \mathbb{R}^{N}$ are linearly independent. The system is called strictly hyperbolic if, in addition, the eigenvalues of $A$ are distinct. We say that the $i$th wave family is genuinely non-linear if ${\rm\nabla}_{u}{\it\lambda}_{i}(u)\cdot r_{i}(u)>0$ for all $u\in \mathbb{R}^{N}$, and it is linearly degenerate if ${\rm\nabla}{\it\lambda}_{i}(u)\cdot r_{i}(u)\equiv 0$ for all $u$.

The first (global) existence result for one-dimensional systems was obtained by Lax (Reference Lax1957) in the special case of Riemann problems, where the initial data consist of two sufficiently close constant states separated by a jump discontinuity. The solution in this case consists of constant states, separated by discontinuous shocks, contact discontinuities or self-similar continuous rarefaction waves.

The fundamental result concerning global existence of entropy solutions subject to general initial data with small total variation (TV) is due to Glimm.

Theorem 2.5 (Glimm Reference Glimm1965).

Assume that the system (2.5) is strictly hyperbolic and that each wave family is either genuinely non-linear or linearly degenerate. Then there exists a constant ${\it\varepsilon}>0$ such that for every initial condition $u_{0}\in L^{1}(\mathbb{R},\mathbb{R}^{N})$ with

(2.6)$$\begin{eqnarray}\text{TV}(u_{0})\leq {\it\varepsilon},\end{eqnarray}$$

the corresponding Cauchy problem (2.5) admits an entropy solution for all $T>0$.

Glimm’s proof relied on constructing approximations using the Glimm scheme (or the random choice method) and obtaining uniform bounds on the total variation of the approximating sequence. Alternative proofs using the wave-front tracking method and the convergence of vanishing viscosity approximations have been obtained by Dafermos (Reference Dafermos2010) and Bianchini and Bressan (Reference Bianchini and Bressan2005), respectively. One advantage of the constructions using either the wave-front tracking method or the vanishing viscosity approximation is the fact that the corresponding operator forms an $L^{1}$ stable semi-group.

Theorem 2.6. There exists a domain ${\mathcal{U}}\subset L^{1}(\mathbb{R},\mathbb{R}^{N})$ containing all functions $u$ with sufficiently small total variation, and a constant $L>0$ such that for all initial data $u_{0}\in {\mathcal{U}}$, approximate weak solutions constructed by front tracking or vanishing viscosity methods converge to a unique entropy solution $u=\mathbf{S}_{t}u_{0}$ of (2.5), with the map $\mathbf{S}_{t}:{\mathcal{U}}\rightarrow {\mathcal{U}}$ forming a semi-group that satisfies

(2.7)$$\begin{eqnarray}\Vert \mathbf{S}_{t}u_{0}-\mathbf{S}_{t}v_{0}\Vert _{L^{1}(\mathbb{R})}\leq L\Vert u_{0}-v_{0}\Vert _{L^{1}(\mathbb{R})}.\end{eqnarray}$$

For the general case of entropy solutions, the above result has been extended by Bressan, Crasta and Piccoli (Reference Bressan, Crasta and Piccoli2000) to obtain the following global uniqueness result.

Theorem 2.7. Let $u\in C([0,T],\,{\mathcal{U}})$ be an entropy solution to (2.5) which satisfies the following tame oscillation condition: there exist constants $C$, ${\it\lambda}>0$ such that, for every $x\in \mathbb{R}$ and every $t,h>0$, we have

(2.8)$$\begin{eqnarray}|u(x,t+h)-u(x,t)|\leq C\cdot \text{TV}_{[x-{\it\lambda}h,x+{\it\lambda}h]}(u(\cdot ,t)).\end{eqnarray}$$

Then $u(\cdot ,t)=\mathbf{S}_{t}u_{0}$. In particular, the entropy solution is unique.

The assumption of initial data with sufficiently small total variation seems to be an essential requirement for the existence of general one-dimensional systems of conservation laws. This assumption can be relaxed, however, for global existence results of special classes of such systems, for example, $2\times 2$ systems ($N=2$ in (2.5)) or systems in the Temple class (Bressan Reference Bressan2000). As in the scalar case, the common approach to establishing existence in these cases is that of bounded variation compactness arguments and their refinement (in particular, compensated compactness arguments: see Tartar Reference Tartar1979), which allow the passage of a strong limit in the non-linear flux $\lim _{n\rightarrow \infty }f(u^{{\it\varepsilon}_{n}})=f(u_{\infty })$.

2.1.3 Multi-dimensional systems

Compared to the scalar case or to one-dimensional systems, there are very few global existence results for multi-dimensional systems. Local (in time) existence and uniqueness results for a system of conservation laws (2.1) with a strictly convex entropy were obtained by Kato (Reference Kato1975). Global existence results for entropy solutions have only been obtained for very special systems; see Benzoni-Gavage and Serre (Reference Benzoni-Gavage and Serre2007) for a comprehensive review. Notable exceptions are the striking results of Chen and Feldman (Reference Chen and Feldman2010) addressing the von Neumann conjecture concerning the global existence of steady two-dimensional Euler equations. However, there are no similar global existence results for physically realistic systems of multi-dimensional time-dependent conservation laws. Indeed, except for those special cases, there does not seem to exist an invariant regularity space which propagates compactness for general classes of multi-dimensional systems of conservation laws ($d\geq 2$), and the question of existence of weak solutions for such systems remains wide open.

2.1.4 Non-uniqueness

The above-mentioned results cover the global existence and uniqueness of entropy solutions for scalar conservation laws, and for one-dimensional systems (with small initial variation). However, there are no global existence results for general classes of multi-dimensional systems, and there are severe difficulties with the question of uniqueness.

Adapting the original construction of De Lellis and Székelyhidi Jr (Reference De Lellis and Székelyhidi2009), a recent paper by Chiodaroli, De Lellis and Kreml (Reference Chiodaroli, De Lellis and Kreml2015) proves the non-uniqueness of solutions of the two-dimensional isentropic Euler equations. This is a reduced system of the compressible Euler equations (1.1), where the density ${\it\rho}$ and velocity $v=(v_{1},\,v_{2})^{\top }$ are governed by the system of conservation law (2.1) with

(2.9)$$\begin{eqnarray}u=\left(\begin{array}{@{}c@{}}{\it\rho}\\ {\it\rho}v_{1}\\ {\it\rho}v_{2}\end{array}\right)\!,\quad f^{j}(u)=\left(\begin{array}{@{}c@{}}{\it\rho}v_{j}\\ {\it\rho}v_{1}v_{j}+p{\it\delta}_{1j}\\ {\it\rho}v_{2}v_{j}+p{\it\delta}_{2j}\end{array}\right)\!,\quad j=1,2.\end{eqnarray}$$

Here, the closure for the pressure is given by the ${\it\gamma}$-law $p({\it\rho})={\it\kappa}{\it\rho}^{{\it\gamma}}$. The system is equipped with an entropy pair $(E,\,v(E+{\it\gamma}p))$, where the total energy

$$\begin{eqnarray}E=\frac{1}{2}\biggl({\it\rho}|v|^{2}+\frac{p}{{\it\gamma}-1}\biggr)\end{eqnarray}$$

serves as the convex entropy.

Theorem 2.8 (De Lellis and Székelyhidi Jr Reference De Lellis and Székelyhidi2009, Chiodaroli et al. Reference Chiodaroli, De Lellis and Kreml2015). There exist initial data $u_{0}\in L^{1}(\mathbb{R}^{2},\mathbb{R}^{3})$ for which the two-dimensional isentropic Euler equations (2.9) have infinitely many entropy solutions.

In principle, this construction can be extended to the full compressible Euler system (1.1).

2.2.1 Classical well-posedness results

One of the earliest theoretical results about the incompressible Euler equations pertains to the following result of Lichtenstein (Reference Lichtenstein1925) on the short-time existence of classical solutions,

Theorem 2.11 (Lichtenstein Reference Lichtenstein1925).

Consider the divergence-free initial data $u_{0}\in C^{1,{\it\alpha}}(D)$ with $0<{\it\alpha}<1$. Then there exists a time $T^{\ast }>0$ depending on $u_{0}$, a velocity field

$$\begin{eqnarray}u\in C_{b}([0,T^{\ast }];C^{1,{\it\alpha}}(D))\cap C^{1}([0,T^{\ast }];C^{0,{\it\alpha}}(D))\end{eqnarray}$$

and pressure $p\in C_{b}([0,T^{\ast }];C^{1,{\it\alpha}}(D))$, such that the pair $(u,p))$ solves the incompressible Euler equations (1.3) in the time interval $[0,T^{\ast }]$, in the classical sense.

Note that the above short-time existence result is valid for both two and three space dimensions. The question of whether one can extend this short-time existence for all times depends on the uniform bound

$$\begin{eqnarray}\int _{0}^{T}\Vert {\rm\nabla}u(\cdot ,t)\Vert _{L^{\infty }}\,\text{d}t<\infty .\end{eqnarray}$$

In fact, certain linear combinations of the entries of ${\rm\nabla}u$ will suffice for global regularity. The following celebrated result of Beale, Kato and Majda (Reference Beale, Kato and Majda1984) provides a sufficient condition for the interval of existence, which reflects the prominent role of the vorticity.

Theorem 2.12 (Beale et al. Reference Beale, Kato and Majda1984).

Let

$$\begin{eqnarray}u\in C([0,T);H^{s}(D)^{3})\cap C^{1}([0,T];H^{s-1}(D)^{3}),\end{eqnarray}$$

with $s>5/2$, be a smooth solution of the three-dimensional Euler equations (1.3) in the time interval $[0,T)$, and assume that

(2.14)$$\begin{eqnarray}\int _{0}^{T}\Vert \text{curl}(u(\cdot ,t))\Vert _{L^{\infty }(D)}\,\text{d}t<\infty .\end{eqnarray}$$

Then $u$ can be uniquely extended to a smooth solution in the interval $[0,T+{\it\delta})$ for some ${\it\delta}>0$.

Thus, it is enough to control the vorticity ${\it\omega}=\text{curl}(u)$, rather than the full gradient ${\rm\nabla}u$, in order to obtain global solutions. In fact, an even more refined result (Constantin, Fefferman and Majda Reference Constantin, Fefferman and Majda1996) shows that control of variations in the direction of the vorticity ${\it\omega}$ suffices to extend the interval of existence of regular solutions. Obtaining an $L^{\infty }$ bound on the three-dimensional vorticity, however, appears to be very difficult, due to possible vortex stretching driven by the right-hand side of the vorticity equation (2.13). In fact, the question of whether smooth solutions of the incompressible Euler equations blow up in finite time is still open. Numerical results (see Hou Reference Hou2008 and references therein) do suggest finite-time blow-up, particularly for problems with no-penetration boundary conditions.

The vorticity equation becomes considerably simpler in the case of two space dimensions, precisely because the vortex stretching term on the right-hand side of (2.13) vanishes. Consequently, one obtains the following global existence and uniqueness results for two-dimensional weak solutions.

For the mere existence of two-dimensional weak solutions, one can relax the assumption on initial regularity by considering larger classes of initial vorticities, ${\it\omega}_{0}\in X$, which are compact in $H^{-1}(D)$, including $X=L^{p}(D)$ and the larger Morrey spaces $X=M^{p}(D)$ with $p>1$ or Lorentz space $X=L^{(1,2)}(D)$ (Lopes Filho, Nussenzveig and Tadmor Reference Lopes Filho, Nussenzveig and Tadmor2000). This brings us to the ‘largest’ borderline space of Radon measures $X={\mathcal{M}}(D)$ which, among others, realizes initial vortex sheets (DiPerna and Majda Reference DiPerna and Majda1987b). In this context, we recall the celebrated existence result of Delort.

The proof is based on mollifying the initial data, resulting in the generation of a sequence of approximate (smooth) solutions to the Euler equations. The resulting vorticity will be of a definite sign as it satisfies a maximum principle. The strong compactness of the approximating sequence is based on localized $L^{1}$ control of the vorticity and uses the fact that the vorticity is of definite sign in a crucial manner; see also Majda (Reference Majda1993) and Schochet (Reference Schochet1995). The question of uniqueness of solutions in the Delort class remains open.

2.2.2 Non-uniqueness

A recent series of remarkable results by De Lellis and Székelyhidi prove the non-uniqueness of admissible weak solutions for the incompressible Euler equations.

The results of De Lellis and Székelyhidi significantly extend the constructions of Scheffer (Reference Scheffer1993) and Shnirelman (Reference Shnirelman2000) that led to the realization of compactly supported (in space-time) non-trivial solutions of the Euler equations. This set of infinitely many admissible weak solutions is based on an intricate construction that adds very high-frequency oscillations in an iterative procedure to construct a convergent sequence of approximate solutions of the Euler equations. The class of initial data, called wild initial data, is very rich and forms a dense subset of $L^{2}$. In particular, vortex sheets are admissible wild initial data in this sense, although the resulting solutions no longer lie in the Delort class. Recent results by Buckmaster et al. (Reference Buckmaster, De Lellis, Isett and Székelyhidi2015) have further refined the original construction of De Lellis and Székelyhidi Jr (Reference De Lellis and Székelyhidi2009), and demonstrate the existence of infinitely many admissible Hölder-continuous weak solutions of the Euler equations.

3.1.1 Probability measures

We let ${\mathcal{M}}(\mathbb{R}^{N})$ denote the set of finite, signed Radon measures on $\mathbb{R}^{N}$ (Folland Reference Folland1999, Chapter 7). Let $C_{0}(\mathbb{R}^{N})$ be the space of continuous real-valued functions on $\mathbb{R}^{N}$ which vanish at infinity, equipped with the supremum norm. By the Riesz–Kakutani theorem, ${\mathcal{M}}(\mathbb{R}^{N})$ can be identified with the dual space of $C_{0}(\mathbb{R}^{N})$ via the pairing

$$\begin{eqnarray}\langle {\it\mu},\,g\rangle =\int _{\mathbb{R}^{N}}g({\it\xi})\,\text{d}{\it\mu}({\it\xi}),\quad {\it\mu}\in {\mathcal{M}}(\mathbb{R}^{N}),~g\in C_{0}(\mathbb{R}^{N})\end{eqnarray}$$

(Folland Reference Folland1999, Section 7.3). We do not distinguish between these two equivalent definitions of ${\mathcal{M}}$. By a slight abuse of notation, we shall sometimes write

$$\begin{eqnarray}\langle {\it\mu},\,g({\it\xi})\rangle =\int _{\mathbb{R}^{N}}g({\it\xi})\,\text{d}{\it\mu}({\it\xi}).\end{eqnarray}$$

For instance, the first moment of ${\it\mu}$ is

$$\begin{eqnarray}\langle {\it\mu},\,{\it\xi}\rangle =\int _{\mathbb{R}^{N}}{\it\xi}\,\text{d}{\it\mu}({\it\xi}).\end{eqnarray}$$

The duality between $C_{0}(\mathbb{R}^{N})$ and ${\mathcal{M}}(\mathbb{R}^{N})$ gives rise to a topology on ${\mathcal{M}}(\mathbb{R}^{N})$, that of weak* convergence. A sequence ${\it\mu}^{n}\in {\mathcal{M}}(\mathbb{R}^{N})$ converges weak* to ${\it\mu}\in {\mathcal{M}}(\mathbb{R}^{N})$ provided $\langle {\it\mu}^{n},\,g\rangle \rightarrow \langle {\it\mu},\,g\rangle$ for all $g\in C_{0}(\mathbb{R}^{N})$. (This is also called vague convergence: see Folland Reference Folland1999.)

The set of probability measures on $\mathbb{R}^{N}$ is the subset

$$\begin{eqnarray}{\mathcal{P}}(\mathbb{R}^{N}):=\{{\it\mu}\in {\mathcal{M}}(\mathbb{R}^{N}):{\it\mu}\geq 0,~{\it\mu}(\mathbb{R}^{N})=1\}.\end{eqnarray}$$

3.1.2 Young measures

A Young measure from $D\subset \mathbb{R}^{k}$ to $\mathbb{R}^{N}$ is a function which maps $z\in D$ to a probability measure on $\mathbb{R}^{N}$. More precisely, a Young measure is a weak*-measurable map ${\it\nu}:D\rightarrow {\mathcal{P}}(\mathbb{R}^{N})$, meaning that:

$$\begin{eqnarray}\text{The mapping }z\mapsto \langle {\it\nu}_{z},\,g\rangle \text{is Borel-measurable for every }g\in C_{0}(\mathbb{R}^{N}).\end{eqnarray}$$

The set of all Young measures from $D$ into $\mathbb{R}^{N}$ is denoted by $\mathbf{Y}(D,\mathbb{R}^{N})$.

In the particular context of our discussion above, Young measures provide a complete description for bounded sequences of approximate solutions $u_{n}=u^{{\it\varepsilon}_{n}}(x,t)$ with weak* limits $u_{n}(x,t)\rightharpoonup u_{\infty }(x,t)$. We ask what can be said about the corresponding weak* limits $f(u_{n})\rightharpoonup f_{\infty }$ for all continuous functions $f\in C_{0}(\mathbb{R}^{N})$. Fix $(x,t)\in D\times \mathbb{R}_{+}$. Since $f_{\infty }(x,t)$ depends linearly and positively on $f\in C_{0}$, it follows that there exists a Young measure ${\it\nu}_{x,t}$ such that $f_{\infty }=\langle {\it\nu}_{x,t},f({\it\xi})\rangle$.

Example 3.1. The most familiar example of Young measures are the so-called atomic Young measures associated with Borel-measurable functions $u:D\rightarrow \mathbb{R}^{N}$:

(3.1)$$\begin{eqnarray}{\it\nu}_{z}={\it\delta}_{u(z)}.\end{eqnarray}$$

The strong convergence $u_{n}\rightarrow u_{\infty }$ is characterized by such atomic measures:

$$\begin{eqnarray}f(u_{n}(x,t))\rightarrow \langle {\it\delta}_{u_{\infty }(x,t)},f\rangle .\end{eqnarray}$$

This is in sharp contrast to the persistence of oscillations in weak convergence.

Example 3.2. Let $u:D\rightarrow \mathbb{R}^{N}$ be a bounded function on the $2{\it\pi}$-periodic torus $D=\mathbb{T}^{k}$, and consider $u_{n}(x)=u(nx)$. Then

$$\begin{eqnarray}u_{n}\rightharpoonup u_{\infty }=\frac{1}{(2{\it\pi})^{k}}\int _{\mathbb{T}^{k}}u(x)\,\text{d}x,\end{eqnarray}$$

and hence the corresponding weak* limits $f_{\infty }$ differ from $f(u_{\infty })$,

$$\begin{eqnarray}f(u_{n})\rightharpoonup f_{\infty }=\frac{1}{(2{\it\pi})^{k}}\int _{\mathbb{T}^{k}}f(u(x))\,\text{d}x,\end{eqnarray}$$

and are described by the Young measure

$$\begin{eqnarray}\langle {\it\nu},f\rangle =(2{\it\pi})^{-k}\int _{\mathbb{T}^{k}}f(u(x))\,\text{d}x.\end{eqnarray}$$

The Young measure ${\it\nu}\in \mathbf{Y}(D,\mathbb{R}^{N})$ is uniformly bounded if there is a compact set $K\subset \mathbb{R}^{N}$ such that $\text{supp}{\it\nu}_{z}\subset K$ for almost every $z\in D$. Note that if ${\it\nu}$ is atomic, ${\it\nu}={\it\delta}_{u}$, then ${\it\nu}$ is uniformly bounded if and only if $\Vert u\Vert _{L^{\infty }(D)}<\infty$.

Just as for the space ${\mathcal{M}}(\mathbb{R}^{N})$, there is a natural topology on $\mathbf{Y}(D,\mathbb{R}^{N})$: a sequence ${\it\nu}^{n}\in \mathbf{Y}(D,\mathbb{R}^{N})$ converges weak* to ${\it\nu}\in \mathbf{Y}(D,\mathbb{R}^{N})$ if $\langle {\it\nu}^{n},\,g\rangle \overset{\ast }{\rightharpoonup }\langle {\it\nu},\,g\rangle$ in $L^{\infty }(D)$ for all $g\in C_{0}(\mathbb{R}^{N})$, that is,

$$\begin{eqnarray}\int _{D}{\it\varphi}(z)\langle {\it\nu}_{z}^{n},\,g\rangle \,\text{d}z\rightarrow \int _{D}{\it\varphi}(z)\langle {\it\nu}_{z},\,g\rangle \,\text{d}z\quad \text{for all}~{\it\varphi}\in L^{1}(D).\end{eqnarray}$$

The fundamental theorem of Young measures was first introduced by Tartar for $L^{\infty }$ bounded sequences (Tartar Reference Tartar1979) and then generalized by Schonbek (Reference Schonbek1982) and Ball (Reference Ball, Rascle, Serre and Slemrod1989) for sequences of measurable functions. A further generalization has been presented by Fjordholm et al. (Reference Fjordholm, Käppeli, Mishra and Tadmor2016a): every sequence ${\it\nu}^{n}\in \mathbf{Y}(D,\mathbb{R}^{N})$ which does not ‘leak mass at infinity’ (condition (3.2)) has a weak* convergent subsequence.

Theorem 3.3. Let ${\it\nu}^{n}\in \mathbf{Y}(D,\mathbb{R}^{N})$ for $n\in \mathbb{N}$ be a sequence of Young measures. Then there exists a subsequence ${\it\nu}^{m}$ which converges weak* to a non-negative measure-valued function ${\it\nu}:D\rightarrow {\mathcal{M}}_{+}(\mathbb{R}^{N})$ in the sense that

1. (i) $\langle {\it\nu}_{z}^{m},\,g\rangle \overset{\ast }{\rightharpoonup }\langle {\it\nu},\,g\rangle$ in $L^{\infty }(D)$ for all  $g\in C_{0}(\mathbb{R}^{N})$,

and moreover satisfies

1. (ii) $\Vert {\it\nu}_{z}\Vert _{{\mathcal{M}}(\mathbb{R}^{N})}\leq 1$ for a.e. $z\in D$, and

1. (iii) if $K\subset \mathbb{R}^{N}$ is closed and $\text{supp}{\it\nu}_{z}^{n}\subset K$ for a.e. $z\in D$ and $n$ large, then $\text{supp}{\it\nu}_{z}\subset K$ for a.e. $z\in D$.

Suppose further that there exists a non-negative function ${\it\kappa}\in C(\mathbb{R}^{N})$ with $\lim _{|{\it\xi}|\rightarrow \infty }{\it\kappa}({\it\xi})=\infty$ such that

(3.2)$$\begin{eqnarray}\sup _{n}\int _{D}\langle {\it\nu}_{z}^{n},\,{\it\kappa}\rangle \,\text{d}z<\infty .\end{eqnarray}$$

Then

1. (iv) $\Vert {\it\nu}_{z}\Vert _{{\mathcal{M}}(\mathbb{R}^{N})}=1$ for a.e. $z\in D$,

whence ${\it\nu}\in \mathbf{Y}(D,\mathbb{R}^{N})$.

The proof of this theorem is given in Appendix A.1 of Fjordholm et al. (Reference Fjordholm, Käppeli, Mishra and Tadmor2016a).

3.3.1 Concentration and generalized Young measures

We begin with the following corollary of Theorem 3.3 (see Remark 2 of Ball Reference Ball, Rascle, Serre and Slemrod1989). Let $N\in \mathbb{N}$, fix a bounded $D\subset \mathbb{R}^{k}$ and let $\{u^{{\it\varepsilon}}\}_{{\it\varepsilon}>0}$ be an $L^{p}$ bounded sequence $u^{{\it\varepsilon}}\in L^{p}(D,\mathbb{R}^{N})$ for some $p>1$. Then, there is a subsequence $u^{{\it\varepsilon}_{n}}$ and a Young measure ${\it\nu}\in \mathbf{Y}(D,\mathbb{R}^{N})$ such that

$$\begin{eqnarray}g(u^{{\it\varepsilon}_{n}})\overset{\ast }{\rightharpoonup }\langle {\it\nu},\,g\rangle \quad \text{in }~~{\mathcal{M}}(D)\text{ as }n\rightarrow \infty\end{eqnarray}$$

for all $g$ with growth strictly slower than $|u|^{p}$, namely,

$$\begin{eqnarray}\int _{D}{\it\varphi}(x)g(u^{{\it\varepsilon}_{n}}(x))\,\text{d}x=\int _{D}{\it\varphi}(x)\langle {\it\nu}_{x},\,g\rangle \,\text{d}x\end{eqnarray}$$

for all $g$ of the form $g(u)=(1+|u|^{p})g_{0}(u)$ with $g_{0}\in C_{0}(\mathbb{R}^{N})$. The above result does not address the question of convergence of $g(u^{{\it\varepsilon}})$ for $g$ with growth equal to $|u|^{p}$,

(3.6)$$\begin{eqnarray}\exists g_{b}\in C_{b}(\mathbb{R}^{N})\quad \text{such that}\quad g(u)=(1+|u|^{p})g_{b}(u)\quad \text{for all}~u\in \mathbb{R}^{N}.\end{eqnarray}$$

What can go wrong in such cases is that $u^{{\it\varepsilon}}$ can concentrate mass at a rate of $|u^{{\it\varepsilon}}|^{p}$ as ${\it\varepsilon}\rightarrow 0$, so that there exist sets $K_{{\it\varepsilon}}\subset D$ of Lebesgue measure $|K_{{\it\varepsilon}}|\rightarrow 0$ with

(3.7a)$$\begin{eqnarray}\int _{K_{{\it\varepsilon}}}|u^{{\it\varepsilon}}(x)|^{p}\,\text{d}x\geq c\quad \text{for all}~{\it\varepsilon}>0\end{eqnarray}$$

for some $c>0$, while for any $q<p$ we have

(3.7b)$$\begin{eqnarray}\displaystyle \int _{K_{{\it\varepsilon}}}|u^{{\it\varepsilon}}(x)|^{q}\,\text{d}x\leq \biggl(\int _{K_{{\it\varepsilon}}}|u^{{\it\varepsilon}}(x)|^{p}\,\text{d}x\biggr)^{q/p}|K_{{\it\varepsilon}}|^{p/(p-q)}\leq C|K_{{\it\varepsilon}}|^{p/(p-q)}\rightarrow 0 & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

as ${\it\varepsilon}\rightarrow 0$. Thus, the concentration of mass is not detected by functions behaving like $|u|^{q}$ for $q<p$. This scenario may be encountered in the setting of systems of conservation laws (2.1) (and the incompressible Euler equations (1.3)), with a generic $L^{2}$ bound $\Vert u(\cdot ,t)\Vert _{L^{2}(\mathbb{R}^{d})}\leq C$. We can view $\{u(\cdot ,t)\}_{t\geq 0}$ as a sequence of functions indexed by $t$, such that $u(\cdot ,t)$ might exhibit concentration of mass at a rate of $|u|^{2}$ as $t\rightarrow T$ for some $T>0$, in the sense that there are spatial subsets $K_{t}\subset \mathbb{R}^{d}$ of Lebesgue measure $|K_{t}|\rightarrow 0$ as $t\rightarrow T$ such that (cf. (3.7a))

$$\begin{eqnarray}\int _{K_{t}}|u(x,t)|^{2}\,\text{d}x\geq c>0\quad \text{for all}~t<T.\end{eqnarray}$$

Example 3.8. Let $T=1$ and

$$\begin{eqnarray}u(x,t):=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{1}{\sqrt{1-t}}\quad & \text{if }|x|<1-t,\\ 0\quad & \text{else},\end{array}\right.\quad \text{for }0\leq t\leq 1.\end{eqnarray}$$

Then, for $t<1$,

$$\begin{eqnarray}\int _{\mathbb{R}}|u(x,t)|^{2}\,\text{d}x=\int _{-(1-t)}^{1-t}|u(x,t)|^{2}\,\text{d}x\equiv 2.\end{eqnarray}$$

Hence, if $K_{t}:=[-(1-t),1-t]$ then (3.7a) holds with $p=2$, $c=2$, while

$$\begin{eqnarray}\int _{\mathbb{R}}|u(x,t)|^{q}\,\text{d}x=2(1-t)^{1-q/2}\rightarrow \left\{\begin{array}{@{}ll@{}}0\quad & \text{if }q<2,\\ \infty \quad & \text{if }q>2,\end{array}\right.\quad \text{as }t\rightarrow 1.\end{eqnarray}$$

Since $u(x,t)\rightarrow 0$ as $t\rightarrow 1$ for every $x\in \mathbb{R}$, the Young measure associated with the sequence $\{u(\cdot ,t)\}_{t<1}$ is simply ${\it\nu}_{x}\equiv {\it\delta}_{0}$, but clearly $\int _{\mathbb{R}}|u(x,t)|^{2}\,\text{d}x$ does not converge to $\int _{\mathbb{R}}\langle {\it\nu}_{x},\,|{\it\xi}|^{2}\rangle \,\text{d}x=0$.

Given the above discussion, we need a framework to describe the limits of functions (or, more generally, Young measures) that exhibit both oscillations as well as concentration. General frameworks which provide a detailed description of both oscillations and concentrations include Tartar’s $H$-measures (Tartar Reference Tartar1990) and Gérard’s microlocal defect measures (Gérard Reference Gérard1991). Here we use the notion of generalized Young measures, introduced by DiPerna and Majda (Reference DiPerna and Majda1987a) precisely to address the question of concentration in the context of the incompressible Euler equations. We follow the general framework of Alibert and Bouchitté (Reference Alibert and Bouchitté1997).

Definition 3.9 (Alibert and Bouchitté Reference Alibert and Bouchitté1997, Székelyhidi Jr and Wiedemann Reference Székelyhidi and Wiedemann2012). Let $D$ be an open set in Euclidean space, and let $S^{N-1}$ denote the unit sphere in $\mathbb{R}^{N}$. A generalized Young measure on $D$ is a triple $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ consisting of a Young measure ${\it\nu}:D\rightarrow {\mathcal{P}}(\mathbb{R}^{N})$, a measure ${\it\lambda}\in {\mathcal{M}}_{+}(D)$ which is singular with respect to Lebesgue measure, and a weak*-measurable map ${\it\nu}^{\infty }:D\rightarrow {\mathcal{P}}(S^{N-1})$ defined for ${\it\lambda}$-a.e. $z\in D$. The objects ${\it\nu}$, ${\it\lambda}$ and ${\it\nu}^{\infty }$ are called the oscillation measure, the concentration measure and the concentration-angle measure, respectively.

3.3.2 Convergence of $L^{1}$ bounded sequence of functions

For the sake of notational simplicity we will consider sequences of functions which are bounded in $L^{1}(D)$, and study functions $g$ with at most linear growth, that is, (3.6) with $p=1$. The general setting is considered in Section 3.3.3. Here and below, $D$ denotes an arbitrary open set in Euclidean space, and ${\mathcal{L}}$ will denote Lebesgue measure on $D$.

We let $B:=\{z\in \mathbb{R}^{N}:|z|<1\}$ denote the open unit ball in $\mathbb{R}^{N}$, and $\bar{B}$ its closure. Let $C(E)$ be the space of continuous, real-valued functions on $E$. We can identify $C(D\times \mathbb{R}^{N})$ with $C(D\times B)$ via

(3.8)$$\begin{eqnarray}\widetilde{g}(x,z):=g\biggl(x,\frac{z}{1-|z|}\biggr)(1-|z|),\quad (x,z)\in D\times B.\end{eqnarray}$$

Note that if $g$ has linear growth (i.e., it satisfies (3.6) with $p=1$) then $\widetilde{g}\in C_{b}(D\times B)$. We let $C_{1}$ denote the set of those $g\in C(D\times B)$ with linear growth for which $\widetilde{g}$ can be extended to a continuous function on $D\times \bar{B}$:

$$\begin{eqnarray}C_{1}:=\{g\in C(D\times \mathbb{R}^{N}):\exists ~{\rm\Phi}\in C(D\times \bar{B})\text{ such that }{\rm\Phi}\big|_{D\times B}=\widetilde{g}\}.\end{eqnarray}$$

By a slight abuse of notation, we denote this unique extension by $\widetilde{g}\in C(D\times \bar{B})$, for $g\in C_{1}$. The recession function of $g$ is defined by $g^{\infty }:=\widetilde{g}\big|_{D\times S^{N-1}}$, the restriction of $\widetilde{g}$ to the unit sphere $S^{N-1}=\partial B$. Thus, for $g\in C_{1}$,

$$\begin{eqnarray}\widetilde{g}(x,z)=\left\{\begin{array}{@{}ll@{}}\displaystyle g\biggl(x,\frac{z}{1-|z|}\biggr)(1-|z|)\quad & \text{if }z\in B,\\ g^{\infty }(x,z)\quad & \text{if }x\in S^{N-1}.\end{array}\right.\end{eqnarray}$$

Conversely, for ${\rm\Phi}\in C(D\times \bar{B})$, we define the inverse operation

(3.9)$$\begin{eqnarray}\widehat{{\rm\Phi}}(x,u):={\rm\Phi}\biggl(x,\frac{u}{1+|u|}\biggr)(1+|u|),\quad (x,u)\in D\times \mathbb{R}^{N}.\end{eqnarray}$$

Then $\widehat{{\rm\Phi}}\in C_{1}$, and in particular, $\widehat{\widetilde{g}}=g$ for all $g\in C_{1}$. Clearly, $C_{0}(D\times \mathbb{R}^{N})\subset C_{1}$, and if $g\in C_{0}(D\times \mathbb{R}^{N})$ then $g^{\infty }\equiv 0$.

Theorem 3.10 (Alibert and Bouchitté Reference Alibert and Bouchitté1997).

Let $\{u^{n}\}_{n}\!\subset \!L^{1}(D,\mathbb{R}^{N})$ be a sequence of functions such that

(3.10)$$\begin{eqnarray}\int _{D}|u^{n}(x)|\,\text{d}x\leq C\quad \text{for all}~n.\end{eqnarray}$$

Then there exists a subsequence, still denoted $\{u^{n}\}_{n}$, and a generalized Young measure $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ such that

(3.11)$$\begin{eqnarray}g(u^{n}){\mathcal{L}}\overset{\ast }{\rightharpoonup }\langle {\it\nu},\,g\rangle {\mathcal{L}}+\langle {\it\nu}^{\infty },\,g^{\infty }\rangle {\it\lambda}\quad \text{in }{\mathcal{M}}(D)\end{eqnarray}$$

for all $g\in C_{1}$, that is,

$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{D}{\it\varphi}(x)g(x,u^{n})\,\text{d}x\rightarrow \nonumber\\ \displaystyle & & \displaystyle \qquad \int _{D}{\it\varphi}(x)\int _{\mathbb{R}^{N}}g(x,u)\,\text{d}{\it\nu}_{x}(u)\,\text{d}x+\int _{D}{\it\varphi}(x)\int _{S^{N-1}}g^{\infty }(x,z)\,\text{d}{\it\nu}_{x}^{\infty }(z)\,\text{d}{\it\lambda}(x)\nonumber\end{eqnarray}$$

for all ${\it\varphi}\in C_{0}(D)$.

3.3.3 Convergence of $L^{p}$ bounded functions

Although the framework described in the previous section only deals with sequences of Young measures (or, as a particular case, functions) that are $L^{1}$ bounded, this theory may be easily extended to $L^{p}$ bounds as follows.

Fix $p\geq 1$. The $L^{p}$ recession function $g^{\infty }$ of a function $g\in C(D\times \mathbb{R}^{N})$ is defined as

(3.12)$$\begin{eqnarray}g^{\infty }(x,z):=\lim _{s\rightarrow \infty }\frac{g(x,sz)}{s^{p}},\quad z\in S^{N-1}\end{eqnarray}$$

(provided the limit exists). We define $C_{p}$ as

$$\begin{eqnarray}C_{p}:=\{g\in C(D\times \mathbb{R}^{N}):g(x,u)=g_{1}(x,u)|u|^{p-1}\text{ for some }g_{1}\in C_{1}\}.\end{eqnarray}$$

Corollary 3.11. Let $\{u^{n}\}_{n}$ be a sequence of functions such that

(3.13)$$\begin{eqnarray}\int _{D}|u^{n}|^{p}\,\text{d}x\leq C\quad \text{for all}~n\end{eqnarray}$$

for some $p\geq 1$. Then there exists a subsequence, still denoted $\{u^{n}\}_{n}$, and a generalized Young measure $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ such that

(3.14)$$\begin{eqnarray}g(u^{n}){\mathcal{L}}\overset{\ast }{\rightharpoonup }\langle {\it\nu},\,g\rangle {\mathcal{L}}+\langle {\it\nu}^{\infty },\,g^{\infty }\rangle {\it\lambda}\quad \text{in }{\mathcal{M}}(D)\end{eqnarray}$$

for all $g\in C_{p}$.

Proof. Define $\widetilde{u}^{n}:=u^{n}|u^{n}|^{p-1}$. Then $\widetilde{u}^{n}$ satisfies the integrability condition (3.10). Let $(\widetilde{{\it\nu}},\widetilde{{\it\lambda}},\widetilde{{\it\nu}^{\infty }})$ be the generalized Young measure associated with $u^{n}$. For a fixed $g\in C_{p}$, let $\widetilde{g}(x,\widetilde{u}):=g(x,\widetilde{u}|\widetilde{u}|^{1/p-1})$, and define $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ as

$$\begin{eqnarray}\langle {\it\nu}_{x},\,g\rangle :=\int _{\mathbb{R}^{N}}\widetilde{g}(x,\widetilde{{\it\xi}})\,\text{d}\widetilde{{\it\nu}}_{x}(\widetilde{{\it\xi}}),\quad {\it\nu}^{\infty }:=\widetilde{{\it\nu}}^{\infty },\quad {\it\lambda}:=\widetilde{{\it\lambda}}\end{eqnarray}$$

for $g\in C_{0}(\mathbb{R}^{N})$. Clearly, both ${\it\nu}$ and ${\it\nu}^{\infty }$ are Young measures. We conclude that

$$\begin{eqnarray}g(x,u^{n})=\widetilde{g}(x,\widetilde{u}^{n})\overset{\ast }{\rightharpoonup }\langle \widetilde{{\it\nu}},\,\widetilde{g}\rangle {\mathcal{L}}+\langle \widetilde{{\it\nu}}^{\infty },\,\widetilde{g}^{\infty }\rangle \widetilde{{\it\lambda}}=\langle {\it\nu},\,g\rangle {\mathcal{L}}+\langle {\it\nu}^{\infty },\,g^{\infty }\rangle {\it\lambda}.\square\end{eqnarray}$$

Theorem 3.10 and Corollary 3.11 can be easily generalized to sequences of ‘$L^{p}$ bounded’ Young measures ${\it\nu}^{n}:D\rightarrow {\mathcal{P}}(\mathbb{R}^{N})$, such that

$$\begin{eqnarray}\langle {\it\nu}_{x}^{n},\,|{\it\xi}|^{p}\rangle =\int _{\mathbb{ R}^{N}}|{\it\xi}|^{p}\,\text{d}{\it\nu}_{x}^{n}({\it\xi})<\infty .\end{eqnarray}$$

Theorem 3.12. Let $p\geq 1$ and let $\{{\it\nu}^{n}\}_{n}$ be a sequence of Young measures such that

(3.15)$$\begin{eqnarray}\int _{D}\langle {\it\nu}_{x}^{n},\,|{\it\xi}|^{p}\rangle \,\text{d}x\leq C.\end{eqnarray}$$

Then there exists a subsequence, still denoted $\{{\it\nu}^{n}\}_{n}$, and a generalized Young measure $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ such that

$$\begin{eqnarray}\langle {\it\nu}^{n},\,g\rangle {\mathcal{L}}\overset{\ast }{\rightharpoonup }\langle {\it\nu},\,g\rangle {\mathcal{L}}+\langle {\it\nu}^{\infty },\,g^{\infty }\rangle m\quad \text{in }{\mathcal{M}}(D)\end{eqnarray}$$

for all $g\in C_{p}$, that is,

$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{D}{\it\varphi}(x)\int _{\mathbb{R}^{N}}g(x,{\it\xi})\,\text{d}{\it\nu}_{x}^{n}({\it\xi})\,\text{d}x\rightarrow \nonumber\\ \displaystyle & & \displaystyle \qquad \int _{D}{\it\varphi}(x)\int _{\mathbb{R}^{N}}g(x,{\it\xi})\,\text{d}{\it\nu}_{x}({\it\xi})\,\text{d}x+\int _{D}{\it\varphi}(x)\int _{S^{N-1}}g^{\infty }(x,z)\,\text{d}{\it\nu}_{x}^{\infty }(z)\,\text{d}{\it\lambda}(x)\nonumber\end{eqnarray}$$

for all ${\it\varphi}\in C_{0}(D)$.

The proof of this theorem is presented in the Appendix.

3.3.4 Generalized measure-valued solutions of incompressible Euler equations

Equipped with above notation, we define generalized measure-valued solutions of the incompressible Euler equations as follows.

Definition 3.13. A family of generalized Young measures $({\it\nu}_{t},{\it\lambda}_{t},{\it\nu}_{t}^{\infty })$ on $D$, parametrized by $t\geq 0$, is a measure-valued solution of the incompressible Euler equations (1.3) with initial data $u_{0}$ if it satisfies

(3.16a)$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{\mathbb{R}_{+}}\int _{D}\langle {\it\nu}_{x,t},\,{\it\xi}\rangle \partial _{t}{\it\varphi}(x,t)+\langle {\it\nu}_{x,t},\,{\it\xi}\otimes {\it\xi}\rangle :{\rm\nabla}_{x}{\it\varphi}(x,t)\,\text{d}x\,\text{d}t\nonumber\\ \displaystyle & & \displaystyle \qquad +\int _{\mathbb{R}_{+}}\int _{D}\langle {\it\nu}_{x,t}^{\infty },\,{\it\theta}\otimes {\it\theta}\rangle :{\rm\nabla}_{x}{\it\varphi}(x,t)\,\text{d}{\it\lambda}_{t}(x)\,\text{d}t+\int _{D}u_{0}(x){\it\varphi}(x,0)\,\text{d}x=0\end{eqnarray}$$

for all ${\it\varphi}\in C_{c}^{\infty }(D\times \mathbb{R}_{+};\,\mathbb{R}^{N})$ with ${\rm\nabla}_{x}\cdot {\it\varphi}=0$, and

(3.16b)$$\begin{eqnarray}\displaystyle \int _{0}^{\infty }\int _{D}\langle {\it\nu}_{x,t},\,{\it\xi}\rangle \cdot {\rm\nabla}_{x}{\it\psi}(x,t)\,\text{d}x\,\text{d}t=0 & & \displaystyle\end{eqnarray}$$

for all ${\it\psi}\in C_{c}^{\infty }(D\times [0,\infty ))$.

Note that if the Young measure is atomic, that is, ${\it\nu}_{x,t}={\it\delta}_{u(x,t)}$ and ${\it\lambda}\equiv 0$, then the definition of measure-valued solutions reduces to the standard notion of weak solutions (2.11).

The global existence of measure-valued solutions was shown by DiPerna and Majda (Reference DiPerna and Majda1987a).

Thus, measure-valued solutions are a framework of solutions for which global existence is guaranteed. Furthermore, they do characterize the (physically relevant) zero-viscosity limit of the Navier–Stokes equations, thus approximating very high Reynolds number flows. On the other hand, the question of uniqueness of measure-valued solutions is more delicate. At the very least, we require the following notion of admissibility.

Definition 3.15. A measure-valued solution $({\it\nu}_{t},{\it\lambda}_{t},{\it\nu}_{t}^{\infty })$ with initial data $u_{0}\in L^{2}(D,\mathbb{R}^{2})$ is called admissible if

(3.17)$$\begin{eqnarray}\int _{D}\langle {\it\nu}_{x,t},|{\it\xi}|^{2}\rangle \,\text{d}x+{\it\lambda}_{t}(D)\leq \int _{D}|u_{0}(x)|^{2}\,\text{d}x\end{eqnarray}$$

for almost all $t\in [0,\infty )$.

Remark 3.16. The $L^{2}$ recession function of the function $g({\it\xi})=|{\it\xi}|^{2}$ is $g^{\infty }\equiv 1$ (see Section 3.3.3), so the ‘squared $L^{2}$ norm’ of $({\it\nu}_{t},{\it\lambda}_{t},{\it\nu}_{t}^{\infty })$ is

$$\begin{eqnarray}\int _{D}\langle {\it\nu}_{x,t},\,g\rangle \,\text{d}x+\int _{D}\langle {\it\nu}_{x,t}^{\infty },\,g^{\infty }\rangle \,\text{d}{\it\lambda}_{t}(x)=\int _{D}\langle {\it\nu}_{x,t},\,|{\it\xi}|^{2}\rangle \,\text{d}x+{\it\lambda}_{t}(D).\end{eqnarray}$$

Hence the left-hand side of (3.17).

Given the fact that weak solutions are also measure-valued solutions (in particular, atomic MV solutions), Theorem 2.15 implies that there are also infinitely many admissible measure-valued solutions of the incompressible Euler equations. However, one can prove the following MV-strong uniqueness result.

Theorem 3.17 (Brenier, De Lellis and Székelyhidi Jr Reference Brenier, De Lellis and Székelyhidi2011). Let $u\in C([0,T];L^{2}(D,\mathbb{R}^{d}))$ (with $d=2$ or 3) be a weak solution of (1.3) with

$$\begin{eqnarray}\int _{0}^{T}\Vert {\rm\nabla}u+{\rm\nabla}u^{\top }\Vert _{L^{\infty }}\,\text{d}t<\infty ,\end{eqnarray}$$

and let $({\it\nu}_{t},{\it\lambda}_{t},{\it\nu}_{t}^{\infty })$ be an admissible measure-valued solution with initial data ${\it\nu}_{0}={\it\delta}_{u_{0}(x)}$ and $u_{0}\in L^{2}(D)$. Then ${\it\nu}_{x,t}={\it\delta}_{u(x,t)}$ and ${\it\lambda}=0$, that is, $v$ is the unique admissible MV solution.

Thus, admissible measure-valued solutions of the Euler equations coincide with classical solutions when they exist. The question of uniqueness (stability) of admissible measure-valued solutions is further investigated in Section 9.

5.5.1 The ELW scheme

The entropy-stable Lax–Wendroff (ELW) schemes introduced in Fjordholm (Reference Fjordholm2013, Section 4.2), are finite difference schemes of the form (5.1a) with a numerical flux of the form

(5.17)$$\begin{eqnarray}F_{i+1/2}:=\widetilde{F}_{i+1/2}^{(k)}-d_{i+1/2}|\unicode[STIX]{x27E6}v\unicode[STIX]{x27E7}_{i+1/2}|^{k-1}\unicode[STIX]{x27E6}v\unicode[STIX]{x27E7}_{i+1/2},\quad d_{i+1/2}>0.\end{eqnarray}$$

Here, the jump $\unicode[STIX]{x27E6}v\unicode[STIX]{x27E7}_{i+1/2}:=v_{i+1}-v_{i}$ is expressed in terms of the entropy variables, $v:={\it\eta}^{\prime }(u)$, and $\widetilde{F}_{i+1/2}^{(k)}$ is a $k$th-order accurate entropy conservative numerical flux which is characterized by the compatibility requirement with the corresponding entropy flux $q(u)$, introduced in Tadmor (Reference Tadmor1987):

(5.18)$$\begin{eqnarray}\displaystyle \langle v_{i+1}-v_{i},\widetilde{F}_{i+1/2}^{(k)}\rangle ={\it\zeta}_{q}(v_{i+1})-{\it\zeta}_{q}(v_{i}),\quad {\it\zeta}_{q}(v):=\langle v,f(\!u(v)\rangle -q(u(v)). & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The schemes based on such fluxes $\widetilde{F}_{i+1/2}^{(k)}$, that is, (5.17) and (5.18) with $d_{i+1/2}\equiv 0$, are entropy conservative in the sense of satisfying the equalities (5.6) and (5.15). This results in the desired entropy inequalities (5.6) and (5.15) by adding a judicious amount of $k$th-order numerical viscosity $d_{i+1/2}|\unicode[STIX]{x27E6}v\unicode[STIX]{x27E7}_{i+1/2}|^{k-1}\unicode[STIX]{x27E6}v\unicode[STIX]{x27E7}_{i+1/2}$ parametrized with positive viscosity amplitudes $d_{i+1/2}>0$.

A general class of such second-order entropy conservative fluxes was constructed in Tadmor (Reference Tadmor2003), and a particularly ‘affordable’ version for the compressible Euler equations was given in Ismail and Roe (Reference Ismail and Roe2009) in terms of the rescaled variables $z:=({\it\rho}/p)^{1/2}(1,v_{1},v_{2},p)^{\top }$:

$$\begin{eqnarray}\displaystyle \widetilde{F}_{i+1/2}^{(2),1} & = & \displaystyle (\overline{z}_{2})_{i+1/2}(z_{4})_{i+1/2}^{\ln },\quad \overline{z}_{i+1/2}:=\frac{1}{2}(z_{i}+z_{i+1}),~z_{i+1/2}^{\ln }:=\frac{\unicode[STIX]{x27E6}z_{i+1/2}\unicode[STIX]{x27E7}}{\unicode[STIX]{x27E6}\log (z)_{i+1/2}\unicode[STIX]{x27E7}},\nonumber\\ \displaystyle \widetilde{F}_{i+1/2}^{(2),2} & = & \displaystyle \frac{(\overline{z}_{4})_{i+1/2}}{(\overline{z}_{1})_{i+1/2}}+\frac{(\overline{z}_{2})_{i+1/2}}{(\overline{z}_{1})_{i+1/2}}\widetilde{F}_{i+1/2}^{(2),1},\nonumber\\ \displaystyle \widetilde{F}_{i+1/2}^{(2),3} & = & \displaystyle \frac{(\overline{z}_{3})_{i+1/2}}{(\overline{z}_{1})_{i+1/2}}\widetilde{F}_{i+1/2}^{(2),1},\nonumber\\ \displaystyle \widetilde{F}_{i+1/2}^{(2),4} & = & \displaystyle \frac{1}{2(\overline{z}_{1})_{i+1/2}}\biggl(\!\frac{{\it\gamma}+1}{{\it\gamma}-1}\frac{1}{(z_{1})_{i+1/2}^{\ln }}\widetilde{F}_{i+1/2}^{(2),1}+(\overline{z}_{2})_{i+1/2}\widetilde{F}_{i+1/2}^{(2),2}+(\overline{z}_{3})_{i+1/2}\widetilde{F}_{i+1/2}^{(2),3}\!\biggr).\nonumber\end{eqnarray}$$

A general recipe to convert such second-order methods to arbitrarily high-order entropy conservative fluxes, $\widetilde{F}_{i+1/2}^{(2)}{\rightsquigarrow}\widetilde{F}_{i+1/2}^{(k)}$, was derived in LeFloch, Mercier and Rohde (Reference LeFloch, Mercier and Rohde2002) and utilized in Fjordholm et al. (Reference Fjordholm, Mishra and Tadmor2012) and Fjordholm (Reference Fjordholm2013). The resulting ELW schemes are readily seen to be (i) (formally) $k$th-order accurate and (ii) entropy-stable in the sense of satisfying the entropy consistency (5.6) and (5.15), and (iii) they can be shown to satisfy the weak BV bounds (5.5) and (5.14) and converge strongly for scalar conservation laws (Fjordholm Reference Fjordholm2013, Propositions 4.2 and 4.3).

5.5.2 TeCNO finite difference schemes

The TeCNO schemes, introduced in Fjordholm et al. (Reference Fjordholm, Mishra and Tadmor2012) and Fjordholm (Reference Fjordholm2013), are finite difference schemes of the form (5.1a) with numerical flux

(5.19)$$\begin{eqnarray}F_{i+1/2}:=\widetilde{F}_{i+1/2}^{(k)}-\frac{1}{2}D_{i+1/2}(v_{i+1/2}^{+}-v_{i+1/2}^{-}).\end{eqnarray}$$

Here, $\widetilde{F}_{i+1/2}^{(k)}$ is the same entropy conservative flux satisfying (5.18), $D_{i+1/2}$ is a positive definite (viscosity) matrix, and $v_{i+1/2}^{\pm }$ are the cell interface values obtained by the following $k$th-order accurate ENO reconstruction (Harten et al. Reference Harten, Engquist, Osher and Chakravarty1987). Specifically, starting with the given piecewise constant approximate solution, $\sum _{i}v_{i}1_{{\mathcal{C}}_{i}}(x)$, one constructs a piecewise polynomial approximation $v^{{\rm\Delta}}(x)$:

$$\begin{eqnarray}\mathop{\sum }_{i}v_{i}1_{{\mathcal{C}}_{i}}(x)~~{\rightsquigarrow}~~v^{{\rm\Delta}}(x):=\mathop{\sum }_{i}p_{i}(x)1_{{\mathcal{C}}_{i}}(x),\end{eqnarray}$$

where $\{p_{i}(x)\}$ are a judicious choice of $k$th-order accurate essentially non-oscillatory (ENO) polynomials introduced in Harten et al. (Reference Harten, Engquist, Osher and Chakravarty1987). We then set the ENO reconstructed pointvalues across the cell interface, $v_{i+1/2}^{\pm }=v^{{\rm\Delta}}(x_{i+1/2}\pm )$, namely,

$$\begin{eqnarray}v_{i+1/2}^{+}=p_{i+1}(x_{i+1/2}),\quad v_{i+1/2}^{-}=p_{i}(x_{i+1/2}).\end{eqnarray}$$

The multi-dimensional (Cartesian) version was designed in Fjordholm et al. (Reference Fjordholm, Mishra and Tadmor2012) (see also Fjordholm Reference Fjordholm2013). The class of TeCNO schemes was shown by Fjordholm et al. (Reference Fjordholm, Mishra and Tadmor2012) and Fjordholm (Reference Fjordholm2013) to have the following properties.

• ∙ They are (formally) $k$th-order accurate.

• ∙ They are entropy-stable. Using the sign property $\text{sign}(v_{i+1/2}^{+}-v_{i+1/2}^{-})=\text{sign}(v_{i+1}-v_{i})$ proved in Fjordholm, Mishra and Tadmor (Reference Fjordholm, Mishra and Tadmor2013), it follows that the TeCNO schemes are entropy-stable in the sense of satisfying the discrete entropy inequality (5.6), and the corresponding two-dimensional inequality (5.15) (Fjordholm et al. Reference Fjordholm, Mishra and Tadmor2012).

• ∙ They have weakly bounded variation when $k=1$ and $k=2$, that is, they satisfy a bound of the form (5.5) in the one-dimensional case and (5.14) in two dimensions; see Fjordholm (Reference Fjordholm2013, Theorem 6.6) and Fjordholm (Reference Fjordholm2013, Section 3.2) for the general multi-dimensional case. The weak total variation bound for $k\geq 3$ depends on a conjectured result for the ENO reconstruction method which remains to be proved; see Fjordholm (Reference Fjordholm2013, Section 5.5).

In light of these properties, the approximate measure-valued solutions generated by the TeCNO scheme converge to an EMV solution of (2.1).

5.5.3 Shock-capturing space-time discontinuous Galerkin schemes

Finite difference and finite volume schemes of the type (5.1a) are particularly suitable for Cartesian grids. For problems which involve spatial domains with complex geometries that necessitate the use of unstructured grids (triangles, tetrahedra), an alternative discretization procedure is the space-time discontinuous finite element procedure of Johnson and Szepessy (Reference Johnson and Szepessy1987), Jaffre, Johnson and Szepessy (Reference Jaffre, Johnson and Szepessy1995), Barth (Reference Barth, Kroner, Ohlberger and Rohde1999) and Hiltebrand and Mishra (Reference Hiltebrand and Mishra2014). In this discontinuous Galerkin (DG) procedure, entropy-stable numerical fluxes expressed in terms of the entropy variables, such as those in (5.19), are used at cell interfaces. Hiltebrand and Mishra (Reference Hiltebrand and Mishra2014) showed that the space-time DG method augmented with a shock-capturing streamline diffusion term satisfies a suitable version of the weak BV bound (5.5) and is entropy-stable in the sense of satisfying the discrete entropy inequality (5.15). Arguing along the lines of Hiltebrand and Mishra (Reference Hiltebrand and Mishra2014, Theorem 4.1), it follows that this DG method converges to a dissipative measure-valued solution of (2.1). We remark that the space-time DG methods are fully discrete, in contrast to semi-discrete finite difference schemes such as (5.1a).

6.1.1 The numerical method $\mathbf{S}_{t}^{N}$

Spectral methods approximate the Euler equations (1.3) in Fourier space (Gottlieb et al. Reference Gottlieb, Hussaini, Orszag and Voigt1984). Let $\mathbb{T}^{d}$ denote the $d$-dimensional Torus. If $(u,p)$ is a solution of the Euler equation (1.3) with periodic boundary conditions, then (formally) it also satisfies the convective form of the equations

(6.1)$$\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}r@{}}\partial _{t}u+u\cdot {\rm\nabla}u+{\rm\nabla}p=0\\ {\rm\nabla}_{x}\cdot u=0\end{array}\right.\quad \text{in }\mathbb{T}^{d}\times \mathbb{R}_{+}. & & \displaystyle\end{eqnarray}$$

Consider the spatial Fourier expansion

$$\begin{eqnarray}u(x,t)=\mathop{\sum }_{k\in \mathbb{Z}^{d}}\widehat{u}_{k}(t)\,\text{e}^{\text{i}k\cdot x}\end{eqnarray}$$

with coefficients given by

$$\begin{eqnarray}\displaystyle \widehat{u}_{k}(t)=\frac{1}{(2{\it\pi})^{d}}\int _{\mathbb{T}^{d}}u(x,t)\,\text{e}^{-\text{i}k\cdot x}\,\text{d}x,\quad k\in \mathbb{Z}^{d}. & & \displaystyle \nonumber\end{eqnarray}$$

If $u$ is a solution of (6.1), the above expression yields

$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\widehat{u}_{k} & = & \displaystyle \frac{1}{(2{\it\pi})^{d}}\int _{\mathbb{T}^{d}}u_{t}\,\text{e}^{-\text{i}k\cdot x}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle -\frac{1}{(2{\it\pi})^{d}}\int _{\mathbb{T}^{d}}(u\cdot {\rm\nabla}u+{\rm\nabla}p)\,\text{e}^{-\text{i}k\cdot x}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle -\frac{\text{i}}{(2{\it\pi})^{d}}\mathop{\sum }_{\substack{ \ell ,m}}(\widehat{u}_{\ell }\cdot m)\widehat{u}_{m}\int _{\mathbb{T}^{d}}\,\text{e}^{\text{i}(\ell +m-k)\cdot x}\,\text{d}x-\frac{\text{i}k}{(2{\it\pi})^{d}}\int _{\mathbb{T}^{d}}p\,\text{e}^{-\text{i}k\cdot x}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle -\text{i}\mathop{\sum }_{\substack{ \ell ,m \\ \ell +m-k=0}}(\widehat{u}_{\ell }\cdot m)\widehat{v}_{m}-\text{i}k\widehat{p}_{k}.\nonumber\end{eqnarray}$$

We note that

$$\begin{eqnarray}\text{div}(u)=\text{i}\mathop{\sum }_{k}(\widehat{u}_{k}\cdot k)\,\text{e}^{\text{i}k\cdot x}=0\end{eqnarray}$$

is equivalent to $\widehat{u}_{k}\bot k$ for all $k$. Using $m=k-\ell$ and $\widehat{u}_{\ell }\bot \ell$ for all terms in the summation, we can thus rewrite the last equation in the form

(6.2)$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\widehat{u}_{k}=-\text{i}\mathop{\sum }_{\substack{ \ell ,m \\ \ell +m-k=0}}(\widehat{u}_{\ell }\cdot k)\widehat{u}_{m}-\text{i}k\widehat{p}_{k}.\end{eqnarray}$$

This is the Fourier space version of the Euler equations (1.3). It becomes evident that the pressure term $-\text{i}k\widehat{p}_{k}$, which is parallel to $k$, serves as the orthogonal $L^{2}$ projection of the non-linear term

$$\begin{eqnarray}-\text{i}\mathop{\sum }_{\substack{ \ell ,m \\ \ell +m-k=0}}(\widehat{u}_{\ell }\cdot k)\widehat{u}_{m}\end{eqnarray}$$

to the orthogonal complement of $k$, thus keeping $u$ divergence-free.

For the coefficient $\widehat{u}_{k}$ with $k=0$, equation (6.2) yields $(\text{d}/\text{d}t)\widehat{u}_{0}=0$. This corresponds to conservation of momentum. Using the Galilean invariance of the Euler equations, we can assume without loss of generality that

$$\begin{eqnarray}\widehat{u}_{0}=\frac{1}{(2{\it\pi})^{d}}\int _{\mathbb{T}^{d}}u\,\text{d}x=0.\end{eqnarray}$$

To obtain a discretized approximation to system (6.2), we restrict our attention to the Fourier modes below some threshold $N$. We thus consider divergence-free fields of the form

$$\begin{eqnarray}u(x,t)=\mathop{\sum }_{|k|\leq N}u_{k}(t)\,\text{e}^{\text{i}k\cdot x},\end{eqnarray}$$

and use the $N$-term Leray projection

$$\begin{eqnarray}\mathbb{P}_{N}\biggl(\mathop{\sum }_{k\in \mathbb{Z}^{2}}\widehat{w}_{k}\,\text{e}^{\text{i}k\cdot x}\biggr):=\mathop{\sum }_{|k|\leq N}(\widehat{w}_{k}-k\widehat{q}_{k})\,\text{e}^{\text{i}k\cdot x},\quad \widehat{q}_{k}:=\frac{\widehat{w}_{k}\cdot k}{|k|^{2}}\end{eqnarray}$$

to ensure the divergence-free constraint $(\widehat{w}_{k}-k\widehat{q}_{k})\bot k$. We can also add a small amount of numerical viscosity to ensure stability of the resulting scheme. This idea results in the following scheme. For initial data $u_{0}(x)$, we obtain an approximate solution $u_{N}(x,t)\approx u(x,t)$ by solving the finite-dimensional problem

(6.3)$$\begin{eqnarray}\begin{array}{@{}r@{}}\partial _{t}u_{N}+\mathbb{P}_{N}\,(u_{N}\cdot {\rm\nabla}u_{N})={\it\varepsilon}_{N}{\rm\Delta}u_{N},\\ u_{N}(x,0)=\mathbb{P}_{N}u_{0}(x).\end{array}\end{eqnarray}$$

In this scheme, the numerical diffusion coefficient ${\it\varepsilon}_{N}>0$ depends on $N$ and goes to zero as $N\rightarrow \infty$.

A refined version of this basic scheme was introduced by Tadmor (Reference Tadmor1989). In this version we choose an integer $m=m(N)\sim {\it\theta}N$ for some ${\it\theta}\in (0,1)$ which serves as a threshold between small and large Fourier modes. We apply a viscous regularization only to the large Fourier modes. With a judicious choice of ${\it\varepsilon}_{N}$ and $m(N)$, the resulting method can be shown to be spectrally accurate (see Tadmor Reference Tadmor1989, Bardos and Tadmor Reference Bardos and Tadmor2015, and Section 3.4 of Tadmor Reference Tadmor2012). We obtain the corresponding spectral viscosity approximation

(6.4)$$\begin{eqnarray}\begin{array}{@{}rcl@{}}\partial _{t}u_{N}+\mathbb{P}_{N}\,(u_{N}\cdot {\rm\nabla}u_{N})\; & =\; & {\it\varepsilon}_{N}{\rm\nabla}\cdot (\mathbb{Q}_{N}{\rm\nabla}u_{N}),\\ u_{N}(x,0)\; & =\; & \mathbb{P}_{N}u_{0}(x).\end{array}\end{eqnarray}$$

Here $\mathbb{Q}_{N}$ denotes a smooth cut-off onto the higher modes, $|k|\geq {\it\theta}N$, for some ${\it\theta}<1$:

$$\begin{eqnarray}\mathbb{Q}_{N}w(x):=\mathop{\sum }_{{\it\theta}N\leq |k|\leq N}{\it\sigma}\biggl(\frac{|k|}{N}\biggr)\widehat{w}_{k}\,\text{e}^{\text{i}k\cdot x},\quad {\it\sigma}(t)=(t-{\it\theta})_{+}^{2p}.\end{eqnarray}$$

The fully viscous spectral method (6.3) corresponds to the special choice of ${\it\theta}=0$ and $p=0$. The application of spectral viscosity for simulations of three-dimensional Euler equations can be found in Karamanos and Karniadakis (Reference Karamanos and Karniadakis2000), Pasquetti (Reference Pasquetti2006) and Avrin and Xiao (Reference Avrin and Xiao2014).

Lanthaler and Mishra (Reference Lanthaler and Mishra2015) derived the following stability estimate.

Lemma 6.1 (stability).

If $u_{N}$ is the solution of the semi-discrete system (6.4), then

(6.5)$$\begin{eqnarray}\displaystyle \frac{1}{2}\Vert u_{N}(t)\Vert _{L^{2}}^{2}+{\it\varepsilon}_{N}\int _{0}^{t}\Vert \mathbb{Q}_{N}{\rm\nabla}u_{N}(s)\Vert _{L^{2}}^{2}\,\text{d}s=\frac{1}{2}\Vert \mathbb{P}_{N}u_{0}\Vert _{L^{2}}^{2}\leq \frac{1}{2}\Vert u_{0}\Vert _{L^{2}}^{2}. & & \displaystyle\end{eqnarray}$$

In particular, we have $\Vert u_{N}(t)\Vert _{L^{2}}\leq \Vert u_{0}\Vert _{L^{2}}$, independent of the choice of $N,m$ and ${\it\varepsilon}_{N}$.

The following consistency estimate was also proved in Lanthaler and Mishra (Reference Lanthaler and Mishra2015).

Lemma 6.2 (consistency).

If $u_{N}$ is the solution of the semi-discrete system (6.4), then for all divergence-free test functions ${\it\varphi}$ in $C_{c}^{\infty }([0,T)\times \mathbb{T}^{d})$,

(6.6)$$\begin{eqnarray}\lim _{N\rightarrow \infty }\int _{\mathbb{R}_{+}}\int _{\mathbb{T}^{d}}\partial _{t}{\it\varphi}\cdot u_{N}+{\rm\nabla}{\it\varphi}:u_{N}\otimes u_{N}\,\text{d}x\,\text{d}t=0.\end{eqnarray}$$

Furthermore, we have uniform Lipschitz continuity

$$\begin{eqnarray}u_{N}\in \text{Lip}([0,T];H^{-d/2-1}(\mathbb{T}^{d})).\end{eqnarray}$$

6.1.2 Approximation of measure-valued solutions

Next, we apply the numerical method described above to the FKMT algorithm. Relabelling ${\rm\Delta}{\rightsquigarrow}N$, we let $\mathbf{S}_{t}^{N}u_{0}:=u_{N}$ denote the approximate spectral viscosity solution computed by (6.4), and apply this to Step 2 of Algorithm 4.1. The following convergence theorem for Algorithm 4.1 was proved in the recent paper by Lanthaler and Mishra (Reference Lanthaler and Mishra2015).

Theorem 6.3. Consider the approximate measure-valued solution ${\it\nu}^{N}$ computed by Algorithm 4.1 using the spectral viscosity method (6.4) as the numerical solution operator $\mathbf{S}_{t}^{N}$. Assume that the (kinetic) energy of the initial Young measure ${\it\sigma}$ is finite, that is,

$$\begin{eqnarray}\int _{\mathbb{T}^{d}}\langle {\it\sigma}_{x},|{\it\xi}|^{2}\rangle \,\text{d}x<\infty .\end{eqnarray}$$

Then ${\it\nu}^{N}$ converges weak* (up to a subsequence) to an admissible measure-valued solution $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ of the incompressible Euler equations (1.3) with initial data ${\it\sigma}$.

The following theorem proves the analogous result for Algorithm 4.3 for atomic initial data.

Theorem 6.4. Let $u_{0}\in L^{2}(\mathbb{T}^{d})$. Let ${\it\nu}^{N,{\it\varepsilon}}$ be the family of approximate measure-valued solutions constructed by Algorithm 4.3, with the spectral viscosity method (6.4) defining the numerical solution operator $\mathbf{S}_{t}^{N}$. Then there exists a subsequence $N_{n}\rightarrow \infty$, ${\it\varepsilon}_{n}\rightarrow 0$ such that

$$\begin{eqnarray}{\it\nu}^{N_{n},{\it\varepsilon}_{n}}\overset{\ast }{\rightharpoonup }({\it\nu},{\it\lambda},{\it\nu}^{\infty })\quad \text{as }n\rightarrow \infty\end{eqnarray}$$

for an admissible measure-valued solution $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ of the incompressible Euler equations (1.3) with atomic initial data $u_{0}$.

It was also shown by Lanthaler and Mishra (Reference Lanthaler and Mishra2015) that the Monte Carlo version (using Algorithm 4.4) also converges to an admissible measure-valued solution of the incompressible Euler equations (1.3).

6.2.1 The numerical method $\mathbf{S}_{t}^{h}$

For simplicity of the exposition, we restrict ourselves to the two-dimensional version of the incompressible Euler equations (1.3) and discretize the domain into Cartesian cells with mesh size ${\rm\Delta}:=({\rm\Delta}x,{\rm\Delta}y)$. These two-dimensional cells,

$$\begin{eqnarray}{\mathcal{C}}_{i,j}:=[x_{i-1/2},x_{i+1/2})\times [y_{j-1/2},y_{j+1/2}),\end{eqnarray}$$

are assumed to have similar mesh sizes, that is, ${\rm\Delta}x=x_{i+1/2}-x_{i-1/2}$ and ${\rm\Delta}y=y_{j+1/2}-y_{j-1/2}$, such that ${\rm\Delta}y=ch,{\rm\Delta}x=h$ for some $c,h>0$. Let

$$\begin{eqnarray}(x_{i},y_{j})=\biggl(\frac{x_{i-1/2}+x_{i+1/2}}{2},\frac{y_{j-1/2}+y_{j+1/2}}{2}\biggr)\end{eqnarray}$$

denote the midpoint of the cells. For a scalar function ${\it\varphi}$ that is represented on the grid with its point values at $c$, we let $\text{grad}^{h}$ denote a cell-centred central, second-order discrete gradient operator defined as

$$\begin{eqnarray}\text{grad}^{h}{\it\varphi}_{i,j}=\left(\begin{array}{@{}c@{}}\frac{{\it\varphi}_{i+1,j}-{\it\varphi}_{i-1,j}}{{\rm\Delta}x}\\ \frac{{\it\varphi}_{i,j+1}-{\it\varphi}_{i,j-1}}{{\rm\Delta}y}\end{array}\right)\!.\end{eqnarray}$$

We introduce the notation ${\rm\Delta}^{h}$ for a standard five-point discrete Laplace operator:

$$\begin{eqnarray}{\rm\Delta}^{h}{\it\varphi}_{i,j}=\frac{{\it\varphi}_{i+1,j}-2{\it\varphi}_{i,j}+{\it\varphi}_{i-1,j}}{{\rm\Delta}x}+\frac{{\it\varphi}_{i,j+1}-2{\it\varphi}_{i,j}+{\it\varphi}_{i,j-1}}{{\rm\Delta}y}.\end{eqnarray}$$

For a function $u=(u_{1},u_{2})$ (the alias for the two-dimensional velocity field), represented on the grid by its point values at $(x_{i},y_{j})$, we define the discrete divergence operator $\text{div}^{h}$ as a cell-centred central divergence operator:

$$\begin{eqnarray}\text{div}^{h}u_{i,j}=\frac{(u_{1})_{i+1,j}-(u_{1})_{i-1,j}}{{\rm\Delta}x}+\frac{(u_{2})_{i,j+1}-(u_{2})_{i,j-1}}{{\rm\Delta}y}.\end{eqnarray}$$

Equipped with the above notation, we define the following fully discrete fractional step predictor–corrector finite difference scheme for approximating the velocity field $u_{i,j}\approx u(x_{i},y_{j})$:

(6.7a)$$\begin{eqnarray}\displaystyle \frac{u^{\ast ,n+1}-u^{n}}{{\rm\Delta}t}+{\mathcal{C}}(u^{n},\bar{u}^{n+1/2}) & = & \displaystyle {\mathcal{D}}\bar{u}^{n+1/2},\end{eqnarray}$$

(6.7b)$$\begin{eqnarray}\displaystyle \frac{u^{n+1}-u^{n}}{{\rm\Delta}t} & = & \displaystyle {\mathcal{P}}\biggl(\frac{u^{\ast ,n+1}-u^{n}}{{\rm\Delta}t}\biggr),\end{eqnarray}$$

(6.7c)$$\begin{eqnarray}\displaystyle \text{grad}^{h}p^{n+1/2} & = & \displaystyle {\mathcal{Q}}\biggl(\frac{u^{n+1}-u^{n}}{{\rm\Delta}t}\biggr).\end{eqnarray}$$

(i) Time-averaged velocity. We use a time-averaged velocity $\bar{u}^{n+1/2}:={\it\theta}u^{n}+(1-{\it\theta})u^{\ast ,n+1},{\it\theta}\in [0,1/2)$ in order to prove some essential stability properties. This time-averaged velocity is used in both the convective and numerical diffusion terms. It makes the convective step (6.7a) implicit.

(ii) Numerical diffusion operator ${\mathcal{D}}$. We use Lax–Wendroff-type diffusion operators of the form

$$\begin{eqnarray}\displaystyle {\mathcal{D}}u & = & \displaystyle \frac{\text{d}_{i+1/2,j}|\unicode[STIX]{x27E6}u_{1}\unicode[STIX]{x27E7}_{i+1/2,j}^{x}|^{k-1}\unicode[STIX]{x27E6}u_{1}\unicode[STIX]{x27E7}_{i+1/2,j}^{x}+d_{i-1/2,j}|\unicode[STIX]{x27E6}u_{1}\unicode[STIX]{x27E7}_{i-1/2,j}^{x}|^{k-1}\unicode[STIX]{x27E6}u_{1}\unicode[STIX]{x27E7}_{i-1/2,j}^{x}}{{\rm\Delta}x}\nonumber\\ \displaystyle & & \displaystyle +\frac{\text{d}_{i,j+1/2}|\unicode[STIX]{x27E6}u_{2}\unicode[STIX]{x27E7}_{i,j+1/2}^{y}|^{k-1}\unicode[STIX]{x27E6}u_{2}\unicode[STIX]{x27E7}_{i,j+1/2}^{y}+d_{i,j-1/2}|\unicode[STIX]{x27E6}u_{2}\unicode[STIX]{x27E7}_{i,j-1/2}^{y}|^{k-1}\unicode[STIX]{x27E6}u_{2}\unicode[STIX]{x27E7}_{i,j-1/2}^{y}}{{\rm\Delta}y},\nonumber\end{eqnarray}$$

expressed in terms of the respective jumps $\unicode[STIX]{x27E6}a\unicode[STIX]{x27E7}^{x}:=a_{i+1,j}-a_{i,j}$ and $\unicode[STIX]{x27E6}a\unicode[STIX]{x27E7}^{y}:=a_{i,j+1}-a_{i,j}$, where $d_{{\it\alpha},{\it\beta}}$ are positive viscosity amplitudes and $k$ is an integer ($k=2$ being a common choice).

(iii) Projection operator ${\mathcal{P}}$. We let ${\mathcal{P}}$ denote a discrete projection operator onto a discretely divergence-free vector space (i.e., such that $\text{div}^{h}u=0$) and ${\mathcal{Q}}=\text{I}-{\mathcal{P}}$ (the gradient, or irrotational, part of the vector field, where $\text{I}$ is the identity map). The discrete projection is implemented by setting ${\mathcal{P}}u:=u-\text{grad}^{h}({\rm\Delta}_{2}^{h})^{-1}\text{div}^{h}u$, where ${\rm\Delta}_{2}^{h}=\text{div}^{h}\text{grad}^{h}$ is a standard cell-centred five-point Laplacian acting on cell-centres in checkerboard formation:

$$\begin{eqnarray}{\rm\Delta}_{2}^{h}{\it\varphi}_{i,j}=\frac{{\it\varphi}_{i+2,j}-2{\it\varphi}_{i,j}+{\it\varphi}_{i-2,j}}{4{\rm\Delta}x^{2}}+\frac{{\it\varphi}_{i,j+2}-2{\it\varphi}_{i,j}+{\it\varphi}_{i,j-2}}{4{\rm\Delta}y^{2}},\end{eqnarray}$$

for a grid function ${\it\varphi}$.

(iv) The non-linear (convective) term ${\mathcal{C}}$. The non-linear convective term is evaluated at two different discrete quantities, namely the ‘lagged’ velocity $u^{n}$ and the time-averaged $\bar{u}^{n+1/2}$. Moreover, we use an energy conservative (second-order, consistent) flux term:

$$\begin{eqnarray}\displaystyle {\mathcal{C}}(u,v)_{i,j} & = & \displaystyle \frac{F^{x}(u_{i,j},~u_{i+1,j},~v_{i,j},~v_{i+1,j})-F^{x}(u_{i-1,j},~u_{i,j},~v_{i-1,j},~v_{i,j})}{{\rm\Delta}x}\nonumber\\ \displaystyle & & \displaystyle +\frac{F^{y}(u_{i,j},~u_{i,j+1},~v_{i,j},~v_{i,j+1})-F^{y}(u_{i,j-1},~u_{i,j},~v_{i,j-1},~v_{i,j})}{{\rm\Delta}y}.\nonumber\end{eqnarray}$$

We employ the following numerical fluxes:

$$\begin{eqnarray}\displaystyle F^{x}(u^{-},u^{+},v^{-},v^{+}) & = & \displaystyle \frac{1}{4}\left(\begin{array}{@{}c@{}}(u_{1}^{+}+u_{1}^{-})(v_{1}^{+}+v_{1}^{-})\\ (u_{1}^{+}+u_{1}^{-})(v_{2}^{+}+v_{2}^{-})\end{array}\right)\!,\quad \nonumber\\ \displaystyle F^{y}(u^{-},u^{+},v^{-},v^{+}) & = & \displaystyle \frac{1}{4}\left(\begin{array}{@{}c@{}}(u_{2}^{+}+u_{2}^{-})(v_{1}^{+}+v_{1}^{-})\\ (u_{2}^{+}+u_{2}^{-})(v_{2}^{+}+v_{2}^{-})\end{array}\right)\!.\nonumber\end{eqnarray}$$

6.2.2 Approximation of measure-valued solutions

We apply the same Algorithm 4.1, but this time with the projection evolution $\mathbf{S}_{t}^{h}$ replacing $\mathbf{S}_{t}^{{\rm\Delta}}$ in Step 2. The following convergence result is proved in Leonardi and Mishra (Reference Leonardi and Mishra2016).

Theorem 6.5. Consider the approximate measure-valued solution ${\it\nu}^{h}$ generated by Algorithm 4.1 using the finite difference projection method (6.7) as the numerical solution operator $\mathbf{S}_{t}^{h}$. Assume that the (kinetic) energy of the initial Young measure ${\it\sigma}$ is finite, that is,

$$\begin{eqnarray}\int _{\mathbb{T}^{n}}\langle {\it\sigma}_{x},|{\it\xi}|^{2}\rangle \,\text{d}x<\infty .\end{eqnarray}$$

Then ${\it\nu}^{h}$ converges weak* (up to a subsequence) to an admissible measure-valued solution $({\it\nu},{\it\lambda},{\it\nu}^{\infty })$ of the incompressible Euler equations (1.3) with initial data ${\it\sigma}$.

The proof is based on several properties of the scheme (6.7), such as the equivalence between the conservative and convective forms of the scheme, $L^{2}$ stability bounds and suitable weak BV estimates. We refer the reader to Leonardi and Mishra (Reference Leonardi and Mishra2016) for further details. Similarly, one can prove convergence of Algorithm 4.3 for computing admissible measure-valued solutions of (1.3) with atomic initial data.

7.1.1 Single realization

As a first test, we fix ${\it\omega}\in {\rm\Omega}$ and compute approximate solutions using the second-order accurate Monte Carlo limiter and HLLC flux-based finite volume scheme of Käppeli et al. (Reference Käppeli, Whitehouse, Scheidegger, Pen and Liebendörfer2011). The computed density, at time $T=4$, on a sequence of successively refined grids ranging from $128^{2}$ to $1024^{2}$, is shown in Figure 1.2. As stated in the Introduction, the solution is quite complex at this time as the leading shock wave has exited the domain but has re-entered from the corners on account of the periodic boundary conditions. Furthermore, this re-entry shock wave interacts with and strongly perturbs the interface, forming a very complex region of small-scale eddy-like structures. As seen in Figure 1.2, there seems to be no convergence as the mesh is refined. This lack of convergence is quantified in Figure 1.3, where we present differences in $L^{1}$ for successive mesh resolutions (1.2), and we see that the approximate solutions for a single sample do not even form a Cauchy sequence, let alone converge.

Figure 7.1.

The mean density for the Richtmyer–Meshkov problem with initial data (7.2) for different grid resolutions at time $t=4$. All results are obtained with $400$ Monte Carlo samples.

7.1.2 Ensemble of solutions

Next, we apply Algorithm 4.3 using the Monte Carlo approximation of Algorithm 4.4, with the perturbed initial data (7.2). We use $M=400$ Monte Carlo samples and compute up to $t=4$ using grid resolutions from $128^{2}$ up to $1024^{2}$ grid points.

In Figure 7.1, for each grid resolution ${\rm\Delta}x$ we plot the mean of the density variable:

(7.4)$$\begin{eqnarray}\bar{{\it\rho}}^{{\rm\Delta}x}(x,t):=\frac{1}{M}\mathop{\sum }_{k=1}^{M}{\it\rho}^{{\rm\Delta}x,k}({\it\omega};x,t),\end{eqnarray}$$

where ${\it\rho}^{{\rm\Delta}x,k}({\it\omega})$ is the mass density of each individual Monte Carlo sample. We observe that small-scale features are averaged out in the mean, and only large-scale structures, such as the strong re-entrant shocks (recall the periodic boundary conditions) and mixing regions, are retained through the averaging process. The figure indicates that, unlike the individual samples shown in Figure 1.2, the mean converges as the mesh is refined. This convergence is quantified in Figure 7.2(a), where we plot the difference in the mean density for successive resolutions:

(7.5)$$\begin{eqnarray}\Vert \bar{{\it\rho}}^{{\rm\Delta}x}(\cdot ,t)-\bar{{\it\rho}}^{{\rm\Delta}x/2}(\cdot ,t)\Vert _{L^{1}}.\end{eqnarray}$$

The figure indicates that this quantity goes to zero, so the approximate means form a Cauchy sequence and hence converge.

Figure 7.2.

Cauchy rates (7.5) for the mean (a) and variance (b) versus grid resolutions ($x$-axis) at time $t=4$ for the Richtmyer–Meshkov problem (7.2). All results are obtained with $400$ Monte Carlo samples.

Figure 7.3 shows the variance of the mass density:

$$\begin{eqnarray}\text{Var}({\it\rho}^{{\rm\Delta}x})(x,t):=\frac{1}{M}\mathop{\sum }_{k=1}^{M}({\it\rho}^{{\rm\Delta}x,k}({\it\omega};x,t)-\bar{{\it\rho}}^{{\rm\Delta}x}({\it\omega};x,t))^{2}.\end{eqnarray}$$

Figure 7.3.

Variance of the density with initial data (7.2) for different grid resolutions at time $t=4$. All results are obtained with $400$ Monte Carlo samples.

As in Figure 7.1, the variances of the approximate Young measures seem to converge as the mesh is refined. This is again quantified in Figure 7.2(b), where the $L^{1}$ differences of the variances at successive mesh resolutions are plotted. Note from Figure 7.3 that the variance is concentrated at the shocks, and even more so in the mixing layer around the original interface.

7.2.1 Convergence as ${\it\varepsilon}\rightarrow 0$

Here we demonstrate the convergence of Algorithm 4.3. First, observe that the perturbed initial data (7.7) converge (in any reasonable metric) to the initial data (7.6) as ${\it\varepsilon}\rightarrow 0$. Following Algorithm 4.3, we wish to study the limit behaviour of approximate solutions ${\it\nu}^{{\rm\Delta}x,{\it\varepsilon}}$ as ${\it\varepsilon}\rightarrow 0$. To this end, we compute approximate solutions using the TeCNO2 scheme (see Fjordholm et al. Reference Fjordholm, Mishra and Tadmor2012) at a very fine mesh resolution of $1024^{2}$ points for different values of ${\it\varepsilon}$.

Figure 7.6.

Approximate sample means of the density for the Kelvin–Helmholtz problem (7.6) at time $t=2$ and different values of perturbation parameter ${\it\varepsilon}$. All the computations are on a grid of $1024^{2}$ mesh points and $400$ Monte Carlo samples.

Results for a single sample at time $t=2$ for different values of ${\it\varepsilon}$ are presented in Figure 7.4. The figures indicate that there is no convergence as ${\it\varepsilon}\rightarrow 0$. The spread of the mixing region seems to remain large even when the perturbation parameter is reduced. This lack of convergence is further quantified in Figure 7.5, where we plot the $L^{1}$ difference of the approximate density for successively reduced values of ${\it\varepsilon}$. This difference remains large even when ${\it\varepsilon}$ is reduced by an order of magnitude.

Next, we apply Algorithm 4.3 using the Monte Carlo approximation of Algorithm 4.4, with the perturbed initial data (7.7). We use $M=400$ Monte Carlo samples and compute up to $t=2$ on a fixed grid resolution of $1024^{2}$ points, using different values of the perturbation parameter ${\it\varepsilon}$. The mean of the mass density (7.4) is plotted in Figure 7.6. The figure clearly shows convergence as ${\it\varepsilon}\rightarrow 0$, to a limit different from the steady-state solution (7.6). This convergence of the mean with respect to decaying ${\it\varepsilon}$ is quantified in Figure 7.7(a), where we compute the $L^{1}$ difference (7.5) of the mean for successive values of ${\it\varepsilon}$. We observe that the approximations of the mean form a Cauchy sequence, and hence converge.

Figure 7.7.

Cauchy rates for the sample mean and the sample variance of the density ($y$-axis) for the Kelvin–Helmholtz problem (7.6) for different values of ${\it\varepsilon}$ ($x$-axis). All the computations are on a grid of $1024^{2}$ mesh points and $400$ Monte Carlo samples.

Similar computations of the variance for different values of ${\it\varepsilon}$ are presented in Figures 7.8 and 7.7(b). These figures clearly show convergence of variance as ${\it\varepsilon}\rightarrow 0$. Moreover, Figure 7.8 clearly indicates that in the ${\it\varepsilon}\rightarrow 0$ limit, the limit variance is non-zero. This strongly suggests that dissipative measure-valued solution can be non-atomic, even for atomic initial data. These results are consistent with the claims of Theorem 5.5.

Figure 7.8.

Approximate sample variances of the density for the Kelvin–Helmholtz instability at time $t=2$ and different values of perturbation parameter ${\it\varepsilon}$. All the computations are on a grid of $1024^{2}$ mesh points and $400$ Monte Carlo samples.

To further demonstrate the non-atomicity of the resulting measure-valued solution, we have plotted the probability density functions (PDFs), approximated by empirical histograms, for the mass density at the points $x=(0.5,0.7)$ and $x=(0.5,0.8)$ in Figure 7.9. We see that the initial unit mass centred at ${\it\rho}=2$ (${\it\rho}=1$, respectively) at $t=0$ is smeared out over time, and at $t=2$ the mass has spread out over a range of values of ${\it\rho}$ between 1 and 2.

Figure 7.10 shows the same quantities, but for a fixed time $t=2$, over a series of meshes. Although a certain amount of noise seems to persist on the finer meshes – most likely due to the low number of Monte Carlo samples – the probability density functions seem to converge with mesh refinement.

Figure 7.9.

Time evolution of the approximate PDF for density ${\it\rho}$ at the points $x=(0.5,0.7)$ (a–e) and $x=(0.5,0.8)$ (f–j) on a grid of $1024^{2}$ mesh points.

Figure 7.10.

The approximate PDF for density ${\it\rho}$ at the points $x=(0.5,0.7)$ (a–d) and $x=(0.5,0.8)$ (e–h) on a series of meshes.