Statistical solutions of the incompressible Euler equations

We propose and study the framework of dissipative statistical solutions for the incompressible Euler equations. Statistical solutions are time-parameterized probability measures on the space of square-integrable functions, whose time-evolution is determined from the underlying Euler equations. We prove partial well-posedness results for dissipative statistical solutions and propose a Monte Carlo type algorithm, based on spectral viscosity spatial discretizations, to approximate them. Under verifiable hypotheses on the computations, we prove that the approximations converge to a statistical solution in a suitable topology. In particular, multi-point statistical quantities of interest converge on increasing resolution. We present several numerical experiments to illustrate the theory.


Introduction
Many interesting incompressible fluid flows are characterized by high to very high values of the Reynolds number. Taking the infinite Reynolds number limit of the underlying Navier-Stokes equations (formally), results in the Incompressible Euler equations governing the motion of an ideal i.e. inviscid and incompressible fluid, given by      ∂ t u + u · ∇u + ∇p = 0, div(u) = 0, u| t=0 = u. (1) Here, the velocity field is denoted by u ∈ R d (for d = 2, 3), and the pressure is denoted by p ∈ R + . The pressure acts as a Lagrange multiplier to enforce the divergence-free constraint [36]. The equations need to be supplemented with suitable boundary conditions. Throughout this work, we shall assume periodic boundary conditions, and take as our domain D, the d-dimensional torus D = T d , d ∈ {2, 3}.  (1) in three space dimensions, even with sufficiently smooth initial dataū, is not yet resolved. Moreover in two space dimensions, where one can prove well-posedness of classical solutions as long as the initial data is sufficiently regular, many interesting initial data of interest, such as vortex sheets, do not possess this regularity. Hence, it is imperative to consider the so-called weak solutions of (1), defined as, Definition 1.1. A vector field u ∈ L ∞ ([0, T ); L 2 (D; R d )) is a weak solution of the incompressible Euler equations with initial data u ∈ L 2 (D; R d ), if T 0ˆD u · ∂ t φ + (u ⊗ u) : ∇φ dx dt = −ˆD u · φ(x, 0) dx, (2) for all test vector fields, φ ∈ C ∞ c ([0, T ) × D; R d ), div(φ) = 0, and D u · ∇ψ dx = 0, for all test functions ψ ∈ C ∞ (D).
It is customary to require additional admissibility criteria in order to recover uniqueness of weak solutions. A natural criterion in this context is given by the so-called dissipative or admissible weak solutions, i.e. weak solutions u such that u(t) L 2 ≤ u L 2 . Although the global existence of admissible weak solutions in three space dimensions is open, one can prove global existence of admissible weak solutions in two dimensions with very general initial data. The most general result of this kind is the celebrated work of Delort [9], where he proved global existence of weak solutions of (1) as long as the initial data was in the so-called Delort class, i.e. with initial vorticity being the sum of a signed Radon measure (in space) and a p-integrable function. On the other hand, it is well-known that admissible weak solutions, in both two and three space dimensions are not unique [40,41,8] and references therein.

Numerical schemes.
A variety of numerical methods have been proposed in order to approximate the incompressible Euler equations (1). These include the spectral and spectral viscosity methods [46], that are particularly suitable for periodic boundary conditions. On bounded domains, efficient methods such as the finite difference projection method [7,3] and discontinuous Galerkin (DG) finite element method [20] have been proposed. Alternatives include numerically approximating the Euler equations in the vorticity-stream function formulation. Methods such as the Lagrangian vortex blob [31], point vortex [39] and central finite difference schemes [26] have been proposed in this context.
Rigorous convergence results for numerical approximations of the incompressible Euler equations (1) are mostly restricted to two space dimensions and to sufficiently regular initial data, see [2] for spectral viscosity methods, [36] for vortex methods and [20] for DG methods. Notable exceptions include [33] where the central schemes of [26] were shown to converge to weak solutions of the two-dimensional Euler equations as long as the initial vorticity was in L p , for p > 1. Similarly in [31,32], the authors showed convergence of vortex point and vortex blob methods in two space dimensions as long as the initial vorticity was a signed Radon measure. Finally in a recent paper [23], the authors showed convergence of a spectral viscosity method for initial data in Delort class.
Furthermore, careful numerical experiments, for instance those presented in [22,15] have shown that numerical schemes may not necessarily converge, even in two space dimensions, for rough initial data. Even when the numerical approximations converge (such as in [23]), the inherent instabilities of the Euler equations lead to a very slow convergence rate (see also Figure 8 (Left)) and render the computation of weak solutions of (1) prohibitively expensive.
1.3. Measure-valued and Statistical solutions. Given the lack of well-posedness results for weak solutions and the lack of convergent numerical approximations, there is considerable scope for the design of alternative solution frameworks for (1). One such framework is that of measure-valued solutions [10], where the sought for solutions are no longer functions but space-time parametrized probability measures on state space. The global existence of measure-valued solutions, even in three space dimensions, was shown in [10]. A convergent numerical method (of the spectral viscosity) type and an efficient algorithm to compute measure-valued solutions was proposed in [22]. However, measure-valued solutions are generically non-unique. This holds true even for the much simpler case of the one-dimensional Burgers' equation [42]. In [12], the authors implicated the lack of information about multi-point (spatial) correlations in the non-uniqueness of measure-valued solutions. Moreover, they also proposed a framework of statistical solutions as an attempt to recover uniqueness.
In the formulation of [12], statistical solutions are time-parameterized probability measures on L p , for 1 ≤ p < ∞, that are consistent with the underlying PDE in a weak sense. They were shown to be equivalent to a family of correlation measures, with the k-th member of this family is a Young measure representing correlations (or joint probabilities) of the solution at k distinct spatial points. Thus, one can interpret statistical solutions as measure-valued solutions, augmented with information about all possible multi-point correlations. A priori, statistical solutions contain much more information than a measure-valued solution. Moreover, statistical solutions encode statistical (ensemble averaged) properties of the solutions of the underlying PDE. Thus, statistical solutions provide a suitable framework for uncertainty quantification (UQ) [12,1]. This is particularly relevant for the incompressible Euler equations as it is well-known that the flow of fluids, at very high-Reynolds numbers, can be turbulent and only averaged (or statistical) properties can be inferred from measurements [19].
Statistical solutions for scalar conservation laws were considered in [12], wherein well-posedness was shown under an entropy condition. In [13,14], a Monte Carlo algorithm, based on the ensemble averaging algorithm of [11], was proposed and analyzed for scalar conservation laws and multi-dimensional hyperbolic systems of conservation laws, respectively. Independent notions of statistical solutions of the incompressible Navier-Stokes equations have been proposed in [18] and in [48]. The zero-viscosity limit of these statistical solutions of Navier-Stokes equations are considered in a forthcoming paper [16]. 1.4. Aims and scope of the current article. The main goal of this article is to propose a suitable notion of statistical solutions for the incompressible Euler equations (1) and to analyze and approximate these statistical solutions. To this end, we will • propose a notion of dissipative statistical solutions, i.e. a time-parameterized probability measure on L p (D) that is consistent with (1) in a suitable sense and prove well-posedness in some special cases, namely local in time wellposedness and global well-posedness in two space dimensions for sufficiently regular initial data, • propose a numerical algorithm, based on ensemble averaging and a spectral viscosity spatial discretization, to approximate statistical solutions of (1) and prove that the approximations converge in an appropriate topology to a statistical solution, under reasonable and verifiable hypotheses on the numerical method, • perform and present extensive numerical experiments to verify the theory and to illustrate interesting properties of statistical solutions.
The rest of the paper is organized as follows: in section 2, we present timeparameterized probability measures on L 2 (D; R d ) and characterize convergence in a suitable topology on this space of measures. In section 3, we define statistical solutions of (1) and present well-posedness results. The numerical approximation of statistical solutions and its convergence is presented in section 4 and numerical experiments are presented in section 5.

2.
Time-parameterized probability measures on L 2 (D; R d ) As mentioned in the introduction, statistical solutions are time-parameterized probability measures on L p . For incompressible Euler equations, the energy bound enforces that p = 2. In this section, we will describe time-parameterized probability measures, characterize them and describe a suitable topology on them.
2.1. Notation. For spatial dimensions d = 2, 3, we denote the spatial domain D = T d , i.e. d-dimensional torus [−π, π] d , with endpoints identified, the co-domain (or phase space) U = R d , so that a solution of (1) is given as a vector field u : D × [0, T ) → U . On occasion x → u(x, t) will also be identified with a 2π-periodic function on R d in the obvious way.
As we are interested in time-dependent vector fields u(x, t), we will fix a final time T > 0 and consider u(x, t) as a mapping u : [0, T ) → L 2 x = L 2 (D; U ). In general, for a Banach space Y , we will denote by L q t (Y ) the space of measurable Of particular interest in the context of the incompressible Euler equations are spaces We will say that ω(r) is a modulus of continuity if ω : [0, ∞) → [0, ∞) is a map such that lim r→0 ω(r) = 0.
On multiple occasions in this paper, we will need to mollify a given function x → u(x). Let us therefore fix a standard mollifier ρ(x) ∈ C ∞ (R d ), supported in a ball of radius 1 around the origin, i.e. ρ ≥ 0 is a function, such that

2.2.
Weak convergence of probability measures. Given a topological space X, we denote by P(X) the space of all probability measures on the Borel σ-algebra of X. Given a sequence µ n ∈ P(X) (n ∈ N) and µ ∈ P(X), we say that µ n converges weakly to µ, written µ n µ, if We shall call a family of probability measures {µ ∆ } ∆>0 ⊂ P(X) defined on a separable Banach space X tight, provided that for any > 0 there exists a compact subset K ⊂ X such that It is a classical result due to Prokhorov (see e.g. theorems 8.6.7, 8.6.8 of the monograph [4]) that a family µ ∆ ∈ P(X), with X a separable Banach space, is tight if and only if µ ∆ is relatively compact under the weak topology.
We recall that weak convergence on a Banach space X is metrized by the p-Wasserstein metric W p : where the infimum is taken over all transfer plans π ∈ Π defined as, for all measurable F, G. More precisely, let µ ∆ be a sequence of probability measures on X. If there exists a M > 0, such that µ ∆ ∈ P(X) are uniformly concentrated on the ball B M (0) = {u ∈ X | u X < M } of radius M and centered at the origin in X, then In the present work, we shall only consider p = 1 and p = 2.
We also note that for p = 1, the following duality formula holds: where the supremum is taken over all 1-Lipschitz continuous functions Φ : X → R.
Next, we characterize the compactness properties of families of probability measures µ ∆ ∈ P(L 2 x ). The following quantity, the so-called structure function, is of particular interest in the present context: We first state the following uniform approximation principle, whose proof is an exercise in general topology.
Lemma 2.1. Let (X, d) be a complete metric space with metric d. Let K ⊂ X. If for any > 0, there exists a mapping i : K → X, such that the image i (K) ⊂ X is precompact (i.e. has compact closure), and d(x, i (x)) < for all x ∈ K, then K is precompact.
We now state our main result characterizing certain compact subsets of P(L 2 x ).
Then the following statements are equivalent: x has compact closure (with respect to the weak topology), (2) There exists a modulus of continuity ω, such that we have a uniform bound on the structure function: Proof. We shall only use the implication (2) =⇒ (1) in the present work. We prove this result here. For a proof of the converse, (1) =⇒ (2), we refer instead to [21].
Step 1: Fix a mollifier ρ for > 0. Define i : P(L 2 x ) → P(L 2 x ) by mollification against ρ , Note that for any fixed > 0, the image i (B M (0)) is precompact in L 2 x . Let us denote by i # : P(L 2 x ) → P(L 2 x ) the push-forward mapping associated with the mollification map i . Note that, since all of the probability measures µ ∈ F are supported on B M (0) ⊂ L 2 x , the corresponding push-forward measures µ :=i #µ have support in the (compact) set i (B M (0)) ⊂ L 2 x . In particular, this implies that, for fixed > 0, the family F := {µ ; µ ∈ F} is tight. By Prokhorov's theorem, this implies that F is precompact in the weak topology.
Step 2: We claim that there exists a constant C > 0 (independent of F and > 0), such that for all µ ∈ F. Here we recall that ω(r) is the uniform modulus of continuity which exists by assumption on F.
To see (7), we note that given any Lipschitz continuous function F : The last integral can be bounded using where C 2 = ρ L ∞ depends only on the (fixed) choice of mollifier ρ. Thus, It follows from the uniform modulus of continuity assumption on F that Taking the supremum over all such F on the left and recalling the duality formula for the 1-Wasserstein distance W 1 , eq. (6), we recover the claimed estimate (7).
Step 3: To conclude the proof of this theorem, we note that by Steps 1 and 2, F is uniformly approximated by precompact sets in the sense of lemma 2.1 (with X = P(L 2 x ) under the 1-Wasserstein distance, K = F and i # : K → X the pushforward induced by mollification), and hence is itself precompact.
2.3. Time parameterized probability measures. As mentioned before, statistical solutions are time-parameterized probability measures. We define them below, x ) and with the property that Denoting by δ 0 the Dirac measure concentrated on 0 ∈ L 2 x , the above condition can equivalently be written asˆT This leads us to define a natural metric on L 1 ([0, T ); P) by We then have the following proposition, whose proof is presented in appendix A.
Our objective is to find a suitable topology on L 1 t (P) and characterize compactness in this topology. It would be natural to try and extend the compactness theorem 2.2 to time-parameterized probability measures and find a suitable version of the weak topology. This necessitates formalizing some notion of time-continuity or time-regularity of underlying functions.
To this end, fix a (time-independent) divergence-free test function φ ∈ C ∞ c (D; R d ). Formally, solutions of the incompressible Euler equations (1) satisfy for s, t ∈ [0, T ), x , in terms of the initial data u 0 . If L > 0 is large enough such that by Sobolev embedding where the constant C > 0 depends only on the initial data. Taking the supremum over all φ ∈ H L x with φ H L x ≤ 1, it follows, at least formally, that (10) Given these considerations, it is natural to assume that statistical solutions of the Euler equations satisfy some version of this time continuity. A formalization is provided in the following definition, Definition 2.5. A weak- * measurable, time-parametrized probability measure t → µ t ∈ P(L 2 x ) is called time-regular, if there exists a constant L > 0, and a mapping s, t → π s,t ∈ P(L 2 x × L 2 x ), such that for almost all s, t ∈ [0, T ): • The measure π s,t is a transport plan from µ s to µ t , • There exists a constant C > 0, such that π s,t satisfies the following regularity condition A family {µ ∆ t } ∆>0 of time-parametrized probability measures is uniformly time-regular, provided that each µ ∆ t is time-regular, and the constants L, C > 0 above can be chosen independently of ∆ > 0.
Remark 2.6. Note that if µ t is of the form with t → u(t) weak solutions of the incompressible Euler equations satisfying (10), then we can define suitable transfer plans The time-regularity property follows from the estimate (10) for the u i .
We now show that a family µ ∆ t , ∆ > 0, of uniformly time-regular probability measures is relatively compact, provided that they satisfy a time-averaged version of the second property of theorem 2.2. To this end, we define the time-averaged structure function of (t → µ t ) ∈ L 1 t (P) (weak- * measurable) as the following quantity: We first state the following lemma: Lemma 2.7. There exists C > 0, such that if µ t ∈ L 1 t (P) has structure function S 2 r (µ t , T ), bounded by a modulus of continuity ω(r), i.e. S 2 r (µ t , T ) 2 ≤ ω(r), then for any > 0, we havê Proof. Follows directly from estimate (8).
The main result of the present section is the following Theorem 2.8. Let µ ∆ t ∈ L 1 t (P) be a family of uniformly time-regular probability measures, for ∆ > 0, for which there exists M > 0, such that µ ∆ t (B M (0)) = 1 for all ∆ > 0, a.e. t ∈ [0, T ). Here B M (0) := { u L 2 x < M }. If there exists a modulus of continuity ω(r) such that The idea behind theorem 2.8 is to use the spatial regularity of the sequence to show that the weak time-regularity assumption of definition 2.5 implies a similar time-regularity with respect to a stronger spatial norm, where H −L is replaced by L 2 . The details of this argument are provided in appendix B.
Let us also remark that a limit µ ∆ t → µ t of a uniformly time-regular sequence µ ∆ t is itself time-regular: Proposition 2.9. Let µ ∆ t ∈ L 1 t (P) be a family of uniformly time-regular probability measures, for ∆ > 0. And such that there exists , then µ t is time-regular (with the same time-regularity constants C, L > 0 as for the family µ ∆ t ).
Proof. To prove proposition 2.9 and according to the definition of time-regularity definition 2.5, we need to show that there exists a map (s, t) → π s,t ∈ P(L 2 x × L 2 x ), which defines a transport plan from µ s to µ t , and satisfies the H −L -estimate of definition 2.5. A proof of existence of such π s,t will be achieved by viewing each π ∆ s,t as an element of P(H −L x × H −L x ), and observing that the π ∆ s,t are concentrated on the compact subset for almost all s, t. Let us now choose a sequence ∆ k → 0. By assumption, for each k ∈ N, there exists a subset I k ⊂ [0, T ) of full Lebesgue measure, such that π ∆ k s,t is concentrated on for all s, t ∈ I k , and satisfieŝ Let I = ∞ k=1 I k . Then I ⊂ [0, T ) is of full Lebesgue measure, and π ∆ k s,t is concen- for all s, t ∈ I and (12) is satisfied for all s, t ∈ I, uniformly for all k ∈ N.
Fix now s, t ∈ I. Since the π ∆ k s,t are concentrated on x ) (and hence relatively compact under weak convergence by Prokhorov's theorem). We can pass to a weakly convergent subsequence π Then, by definition of weak convergence, we have for any When passing to the limit in the last step, we have used the convergence µ ∆j s µ s as elements of P(L 2 x ), as well as the fact that if F ∈ C b (H −L x ), then the restriction x ) that proj 1 #π s,t = µ s . A similar argument shows that also proj 2 #π s,t = µ t . In particular, this also implies that π s,t is in fact concentrated on Then the complement of the Cartesian In the above estimate, the equality on the second line follows from the fact that proj 1 #π s,t = µ s and proj 2 #π s,t = µ t , and the last equality follows from the fact that µ s , µ t are supported on X = L 2 Since s, t ∈ I were arbitrary and as I ⊂ [0, T ) is of full Lebesgue measure, we have thus shown that for almost all s, t ∈ [0, T ), there exists a transfer plan π s,t from µ s to µ t satisfying the assumptions of definition 2.5. This shows that the limit µ t is itself time-regular.

2.4.
Time-dependent Correlation measures and their compactness. It has been shown in [12] that there is a one-to-one correspondence between probability measures on L 2 x and so-called correlation measures. Correlation measures are defined as infinite hierarchies of Young measures, taking into account spatial correlations, or more precisely, Definition 2.10. A correlation measure is a collection ν = (ν 1 , ν 2 , . . . ) of maps ν k : D k → P(U k ) satisfying the following properties: (1) Weak- * measurability: Each map ν k : D k → P(U k ) is weak- * -measurable, in the sense that the map x → ν k x , f from x ∈ D k into R is Borel measurable for all f ∈ C 0 (U k ) and k ∈ N. In other words, ν k is a Young measure from D k to U k .
It has been shown in [12], that if µ ∈ P(L 2 x ), then we can associate to it a unique correlation measure ν, with the interpretation that for A 1 , . . . , A k ⊂ U : . More precisely, we have the following theorem [12]: Theorem 2.11. For every correlation measure ν ∈ L 2 (D, U ) there exists a unique probability measure µ ∈ P(L 2 (D; U )) satisfying for all g ∈ C 0 (D k ×U k ) and k ∈ N (where u(x) denotes the vector (u(x 1 ), . . . , u(x k ))). Conversely, for every probability measure µ ∈ P(L 2 (D; U )) with finite moment (15), there exists a unique correlation measure ν ∈ L 2 (D, U ) satisfying (16). The relation (16) is also valid for any measurable g : Moreover, the moments . . ⊗ ξ k , uniquely determine the correlation measure ν and hence the underlying probability measure µ.
The following result now follows from theorem 2.8: there exists a subsequence ∆ j → 0 (j ∈ N), and a time-parametrized probability measure µ t ∈ L 1 t (P), such that Furthermore, denoting by ν t = (ν 1 t , ν 2 t , . . . ) the correlation measure corresponding to the limit µ, we have t,x,y , |ξ 1 − ξ 2 | 2 dy dx dt ≤ ω(r), • We define admissible observables, in terms of test functions g ∈ C([0, T ) × D k × U k ), which satisfy the following bounds, Then, these admissible observables converge strongly in L 1 t,x , in the sense that Remark 2.13. We note that the uniform modulus of continuity estimate in theorem 2.12 can equivalently be expressed aŝ We now come to the proof of theorem 2.12.
Proof of theorem 2.12. The compactness property follows easily from theorem 2.8, above. We only provide a proof of strong convergence of the observables here. We drop the subscript j in the following, and assume that ∆ → 0 is a sequence such that µ ∆ t → µ t in L 1 t (P). We recall that by the assumption of this theorem, there exists M > 0, such that µ ∆ t (B M (0)) = 1, for all ∆. Fix now g ∈ C([0, T )×D k ×U k ) as in the statement of the theorem. Denote By the assumed bound on x and > 0, let u denote the mollification of u. For fixed > 0, and x , and so is u → g(t, x, u (x)). Here we recall that in our notation, Indeed, using the assumed Lipschitz bound on ξ → g(t, x, ξ), we find for u, v ∈ B M (0) from (18), (19): Using the Hölder estimate |u (x)| ≤ ρ L 2 x u L 2 x for the point-values of the mollification, it follows that The previous observation allows us to define and similarly U (t, x) for µ t . Then, we have for We can thus integrate over these variables to find x ≤ M }, by assumption. On the other hand, from Hölder's inequality, we havê on the support of µ ∆ t . Combining these estimates, we conclude that In the last estimate, we have used Jensen's inequality. Finally, we recall that by (8), there exists a constant C > 0, such that ˆT We have thus shown that where C = C(g, k, M, T ) depends on the observable g(t, x, ξ), the number of pointcorrelations k, the uniform a priori L 2 x -bound M > 0, and the maximal time T > 0, but is independent of ∆, and φ.
Taking the supremum over all φ The same inequality holds also true for U (x, t), following the same argument but dropping the subscript ∆ in the above estimates. By lemma 2.7, and the uniform (in ∆) upper bound on the structure function, it now follows that there exists an absolute constant C > 0, such that From the convergence µ ∆ t → µ t , it furthermore follows that for any > 0 fixed: Indeed, this follows from the fact that where F x, (u) has been defined in (20) above, and the following estimate: By (21), we have F x, Lip < ∞, so that last term converges to 0 as ∆ → 0, as follows from the fact that µ ∆ t → µ t in L 1 (P). Thus, combining (22), (23) and (24), we find for any fixed > 0. Finally, noting that the left-hand side is independent of , we may let → 0 to show that

Dissipative Statistical solutions and their well-posedness
Given the discussion on time-parameterized probability measures in the last section, we can now define statistical solutions of (1) as, is a statistical solution of the incompressible Euler equations with initial dataμ, if t → µ t is time-regular, and the associated correlation measure ν t satisfies: Hereν is the correlation measure corresponding to the initial dataμ. We denote F (ξ) := ξ ⊗ ξ and the contraction in the second term is more explicitly given by The above PDEs specify the time-evolution of the moments (17) for all k and by theorem 2.11, determine the evolution of the probability measure µ t . Remark 3.2. As ν 1 above is a standard Young measure, it is straightforward to observe that the corresponding identity for the evolution of ν 1 corresponds to the definition of measure-valued solution of (1) in the sense of [10] under the further assumption that there is no concentration. Hence, one can think of statistical solutions as measure-valued solutions coupled with information about all possible multi-point correlations.
We first show that the second property of definition 3.1 is equivalent to the requirement that µ t be supported on divergence-free vector fields for almost all t.
x ), with associated correlation measure ν. Then µ is concentrated on divergence-free vector fields if, and only if, Proof. Let ψ ∈ C ∞ c (D). Then we have the following identitŷ Therefore, if µ is concentrated on the divergence-free vector fields, then D u · ∇ψ dx = 0, µ-almost surely, and hence, from equation (25), we obtain To prove the converse, let us assume that relation (26) We note that for any > 0, we have where the last equality follows from (25) and the assumption (26). Letting → 0, we conclude that µ [F n (u) > 0] = 0 for all n ∈ N. Equivalently, we have µ ˆD u · ∇ψ n = 0 = 0, for all n ∈ N.
In analogy with work of [12,14] on entropy statistical solutions for hyperbolic systems of conservation laws, we define , div(φ) = 0, and all i = 1, . . . , M . And, in addition, we have for almost every t ∈ [0, T ): 3.1. Existence and uniqueness of dissipative solutions. In this section, we show based on purely topological arguments, that if the set of C 1 -regular initial data admitting classical solutions of (1), over a given time-interval [0, T ) is dense in L 2 x , then there exists a (topologically) generic set G ⊂ L 2 x , containing these regular initial data, with the following property: For any initial dataμ ∈ P(L 2 x ) which is concentrated on this generic set G ⊂ L 2 x , we have existence and uniqueness in the class of dissipative statistical solutions. By a "generic" set G, we denote a set whose complement E = L 2 x \ G is a countable union of nowhere dense sets (implying that E is a meagre set in the topological sense). We say thatμ is concentrated on G, if µ(G) = 1.
The construction of such a generic G under the above mentioned assumption has first been carried out in [30]. Let us first review the construction of G. We let Finally, we let G = n∈N G n .
Remark 3.5. If there exists a dense set of initial data v ∈ C, then G is generic in the topological sense (more precisely a G δ set), being the countable intersection of the dense open sets G n . By the Baire category theorem, the set G is non-empty and dense in this case. In particular, this would hold true if there is no finite-time blowup for sufficiently smooth classical solutions of the incompressible Euler equations (e.g. for C 1,α initial data v possessing a Hölder continuous derivative), which is an established fact in two space dimensions but an open question in three space dimensions.
We can now state the main theorem of the present section: Theorem 3.6. Define the generic set G as above (cp. equation (27)). Ifμ ∈ P(L 2 x ) is an initial datum such thatμ(G) = 1 and there exists M > 0 such that µ(B M (0)) = 1, then there exists a unique dissipative statistical solution µ t of the incompressible Euler equations with initial dataμ.
The proof of theorem 3.6 is somewhat technical and is left to appendix C.
Remark 3.7. Even though a topologically generic set G is large in some sense (e.g. in the sense that a countable intersection of generic sets is again generic, and in particular non-empty), it might be very small from a different point of view. For example, there exist generic subsets of R n , which have Lebesgue measure zero. It is therefore a priori not clear whether there are any "interesting"μ, such thatμ(G) = 1, beyond thoseμ which are concentrated on smooth initial data. Nevertheless, theorem 3.6 implies suitable weak-strong uniqueness results as we show below.
We can derive the following corollaries from theorem 3.6.
Corollary 3.8 (Short-time existence and uniqueness). If m ≥ d/2 + 2, and if there exists a C > 0, such that µ ∈ P(L 2 x ) is concentrated on Proof. Classical short-time existence results for the Euler equations show [36] that there exists T * > 0, such that for initial data Since H m x →C 1 , this implies that µ is concentrated on u ∈ C. In particular, we conclude that µ(G) = 1, and the result now follows from theorem 3.6. Corollary 3.9 (Weak-Strong uniqueness in 2d). Let d = 2, and let α ∈ (0, 1). If µ is concentrated on C 1,α (D; U ) and if there exists M > 0, such that µ(B M (0)) = 1, then there exists a dissipative statistical solution µ t with initial data µ. Furthermore, µ t is unique in the class of dissipative statistical solutions with initial data µ.
Proof. Again, we observe that for any u ∈ C 1,α , there exists a unique solution u(t) ∈ C 1,α . Hence, we have u ∈ C for all such u. In particular, it follows that µ is concentrated on G. The claim follows from theorem 3.6.

Numerical approximation of statistical solutions
In this section, we will propose an algorithm for computing statistical solutions of the incompressible Euler equations (1). As mentioned before, this algorithm is very similar to the one proposed in [14] for computing statistical solutions of hyperbolic systems of conservation laws, which in turn was inspired by the ensemble averaging algorithms of [11,22] for computing measure-valued solutions. This algorithm requires a spatio-temporal discretization and a Monte Carlo sampling of the underlying probability space. We propose to use a spectral viscosity spatial discretization which is described below. 4.1. Spectral hyper-viscosity scheme. We write u ∆ (x, t) = |k|∞≤N u ∆ k (t)e ik·x , where now and in the following we shall consistently denote ∆ = 1/N , and we denote |k| ∞ := max i=1,...,d |k i |. We consider the following spectral viscosity approximation [46,47] of the incompressible Euler equations, Here P N is the spatial Fourier projection operator, mapping an arbitrary function f (x, t) onto the first N Fourier modes: and we assume 0 ≤ Q k ≤ 1.
The idea behind the SV method is that dissipation is only applied on the upper part of the spectrum, i.e. for |k| > m N , thus preserving the formal spectral accuracy of the method, while at the same time enabling us to enforce a sufficient amount of energy dissipation on the small scale Fourier modes which is needed to stabilize the method. The additional hyperviscosity parameter s ≥ 1 in (28) can be chosen larger to enforce more numerical dissipation on the high Fourier modes, thus allowing a larger part of the Fourier spectrum to remain free of numerical diffusion, while still ensuring stability of the resulting numerical scheme.
Note that Q N is defined via its Fourier coefficients Q k , which, given an additional parameter θ > 0, are assumed to satisfy the following constraints: In [47], the parameters m N , N , θ are chosen such that Multiplying the evolution equation (28) by u ∆ and integrating by parts, we obtain the following energy balance, 4.2. Monte Carlo algorithm. Following [14], the computation of statistical solutions of (1) requires combining the spectral viscosity scheme (28) with the following Monte Carlo sampling, denotes the solution operator, defined by the scheme (28). • The approximate statistical solution µ ∆ t is given by the so-called empirical measure, We remark that in practice, the samplesū i for 1 ≤ i ≤ m are random realizations with respect to a certain underlying probability space.
Remark 4.2. The Monte Carlo algorithm 4.1, when restricted only to the computation of the first correlation marginal ν 1 , reduces to the ensemble averaging algorithm proposed in [22] for computing measure-valued solutions of the incompressible Euler equations.

4.3.
Convergence to Statistical Solutions. In this section, we will investigate the convergence of the empirical measure µ ∆ t (33), generated by the Monte Carlo algorithm 4.1, to a statistical solution of (1). To this end, we seek to apply the convergence theorem 2.8 to these approximations. We start by verifying the temporal regularity of the empirical measures in the following lemma, Lemma 4.3. There exists L ∈ N, such that if u ∆ is obtained from the spectral hyper-viscosity method (28), with ∆ = 1/N and initial data u ∈ L 2 x , then Proof. From the definition of the spectral hyperviscosity scheme (28), we obtain, . Clearly, the first error term E ∆ 1 can be estimated by x ), where we have used that N ∼ N −2s+1 ≤ N −1 in the last step (assuming s ≥ 1). On the other hand, to estimate the second term, let φ ∈ C ∞ (D; U ) be a given vector field. Then x . Choosing > 0 sufficiently large, we have by the Sobolev embedding theorem H x →L ∞ x . Hence, we can further estimate for some absolute constant C > 0: Thus, we have shown that x . Let now L := max(2s, + 2). Then, combining the above two estimates, and noting that H x , we now conclude that there exists an absolute constant C > 0, such that For the time-regularity estimate, we write By an argument analogous to the above, we then find that x )|t − s|. This concludes our proof.
From lemma 4.3, it is now easy to see that if µ ∆ t is generated by the Monte-Carlo algorithm 4.1, i.e.
with u ∆ i (t) computed by the spectral hyper-viscosity scheme (28), then the transport plan defined by satisfies the properties required by the definition of time-regularity, definition 2.5. This provides the required temporal regularity required by theorem 2.12. Next, we turn our attention to the spatial regularity bounds of theorem 2.8. In particular, we need to obtain uniform estimates on the structure function (11). We start with the following simple observation, Lemma 4.4. For any r ≥ 0, we have This implies that for |h| ≤ r, and τ ∈ [0, 1]: Furthermore, the upper bound |f (t)| ≤ 2 is obvious, so that |e ik·h − 1| 2 ≤ 4 min(|k| 2 r 2 , 1).
Averaging over h ∈ B r (0), it follows that The next result is an estimate on the structure function (11) at the grid scale ∆.
Lemma 4.5. If µ ∆ t is an approximate statistical solution obtained from the spectral hyper-viscosity method with ∆ = 1/N , and initial dataμ for which there exists , for some absolute constant C > 0. The same estimate is also true for r ≤ ∆, i.e. we have Proof. Our goal is to estimate the structure function S 2 r (µ ∆ t , T ) based on the energy estimate (32). By construction, we have and u ∆ i is obtained by the spectral hyper-viscosity method (28) with initial data u i . By Plancherel's theorem, we obtain where the last estimate is shown in lemma 4.4 below. Let us now fix some α with θ ≤ α < 1. We split the summation over modes |k| ≤ 2N α and |k| > 2N α : Recall that for all |k| ≥ 2N α ≥ 2m N we have as a consequence of (30), (31). We thus find where we have used the prescribed scaling N ∼ N −2s+1 in the last step. Combining both estimates, and noting that for all i, u i L 2 ≤ M by assumption, we conclude that with an implied constant that depends on s, T , but is independent of the initial data, and N . If we now choose r = ∆ = 1/N , so that ∆ is a length at the grid-scale, then we find We finally would like to determine a constant α, which minimizes this last expression. This is achieved when the term in brackets is of order 1, i.e. with α = 2s−1 2s . With this choice of α, we have 2α − 2 = −1 s , and Let us finally sum over i = 1, . . . , m. Thus, if µ ∆ t is an approximate statistical solution obtained from the spectral vanishing hyperviscosity method of order s, and initial data µ, with ∆ = N −1 , then it satisfies the following bound on the structure function at the grid scale: Inspection of (35) shows that also As in [14] section 4.2, we have uniform estimates on the structure function at (or below) the grid scale. Large scale features are in any case independent of the resolution ∆. However, we lack any information on the intermediate scales, in between the two. To close this information gap, we follow [14] and make an assumption on scaling of the structure function (11) . Then the approximate statistical solutions µ ∆ t converge (up to a subsequence still denoted by ∆), as ∆ → 0, to some µ t ∈ L 1 t (P).
Proof. By lemma 4.5, there exists a constant C > 0, such that for all r ≤ ∆. If r > ∆, then by the assumed approximate scaling property, we write r = ∆, with > 1, and obtain for some constant CC > 0. The convergence now follows from theorem 2.12.
Remark 4.7. The scaling assumption (37) can be interpreted as a weaker version of the scaling assumptions of Kolmogorov (see hypothesis H2, equation 6.3, page 75 of [19]) that was instrumental in the K41 theory for homogeneous, isotropic turbulence. We also note that the inequalities in (37) can accommodate intermittency in the form of deviations for the standard Kolmogorov determination of the exponent 1/3 for the structure function (11). 4.4. Decay of energy spectrum. In this section, we will provide an alternative criterion to ensure convergence of probability measures with respect to the metric (9).
This criterion is motivated from well-known experimental and theoretical concepts in the study of turbulent flows and is based on the energy spectrum E(u; K) (K ∈ N 0 ) associated to a vector field u, defined as Note that the kinetic energy is obtained as a sum Given a probability measure µ ∈ P(L 2 x ), let us similarly define: so that E(δ u ; K) = E(u; K), for u ∈ L 2 x . Finally, we denote by E T (µ t ; K) the time-integrated energy spectrum It is an experimentally observed fact [19] that the typical energy spectrum of turbulent flows with a sufficiently strong dissipation mechanism at small scales, typically takes a shape similar to the one shown in Figure 1: Visible are three parts of the energy spectrum. The left-most part (small K) corresponds to large-scale features for the flow, the middle part (intermediate K) is referred to as the inertial range, while the right-most part (large K) may be referred to as the dissipation range. The appearance of these three parts is heuristically explained as follows. Starting from initial data (with a sufficiently fast decay of the energy spectrum) initially fixes the large-scale features of the flow. Due to the non-linear nature of the evolution equation, these large-scale features decay to smaller scales, corresponding to energy cascading from small values of K to larger values of K. While a satisfactory mathematical treatment of the precise nature of this energy cascade remains an outstanding challenge, there is evidence by physical reasoning and as well as from numerical and real-world experiments that typically the energy spectrum resulting from this cascade process satisfies at least an upper bound of the form E(K) K −γ , for some fixed γ that is associated with the non-linearity. In the presence of a dissipative mechanism acting on small scale features of the flow, this "free" energy cascade to larger values of K due to the non-linearity is finally interrupted by the dissipation. Thus, energy is dissipated at dissipative scales.
From this heuristic point of view, we would expect the large-scale features to depend mostly on the initial data, while the decay of the energy spectrum at the largest values of K can be controlled in a numerical approximation scheme by a suitable choice of the numerical dissipation. On the other hand, there is no a priori information on the decay of the spectrum in the intermediate, inertial range. Hence, we make the following, rather natural, assumption, Assumption 4.8. There exist β > 0 and constant C > 0 such that the computed energy spectra with algorithm 4.1 scale as, Under this assumption on the energy spectrum, we have the following convergence theorem, Theorem 4.9. If µ ∆ t is obtained by the spectral viscosity method through algorithm 4.1, and if the energy spectra E T (µ ∆ t ; K) satisfy the inertial range assumption 4.8 with β > 1/2, then there exists a subsequence (not relabeled) ∆ → 0 and a time-parametrized probability measure µ t , such that µ ∆ t → µ t in L 1 t (P).
Proof. From Plancherel's identity and lemma 4.4, we have Hence, based on the assumption 4.8, we now obtain the estimate, Therefore, the scaling assumption on the average energy spectrum leads to the uniform diagonal continuity: From theorem 2.12, we obtain compactness of the sequence µ ∆ t .  Figure 1 (B), a convenient way to check the scaling assumption 4.8 in practice is to consider the compensated energy spectrum, which is defined as K γ E(K), where γ is the (proposed) scaling exponent in the inertial range. proposition 4.9 says that if there exists γ > 1, such that the compensated energy spectrum K γ E T (µ ∆ t ; K) is uniformly bounded by a constant, and independently of ∆, then {µ ∆ t | ∆ > 0} is compact in L 1 (P).
Remark 4.11. If d = 3 and p = 2, then Kolmogorov's theory states that for fully developed turbulence S 2 r ∼ r 1/3 . Based on our estimate, this requires β = 5 6 . So that the (expected) energy spectrum is E(K) ∼ K −2β ∼ K −5/3 . Such an assumed scaling is consistent with many real, as well as numerical, experiments reported in the literature, and is sufficient for compactness in the space of probability measures L 1 t (P) (cp. proposition 4.9).

4.5.
Lax-Wendroff type theorem. We have used a compactness argument to show that under some reasonable hypotheses on the approximations, numerical solutions computed by the spectral hyper-viscosity converge to a limiting timeparametrized probability measure. In this section, we show that such a limit necessarily is a statistical solution of the incompressible Euler equations in the sense of definition 3.1.
Theorem 4.12 (Lax-Wendroff type theorem). Let µ ∆ t be computed by the spectral hyper-viscosity scheme with initial dataμ, and assume µ ∆ t → µ t in L 1 t (P), as ∆ → 0. Then µ t is a statistical solution of the incompressible Euler equations with initial dataμ.
) be given solenoidal test functions. Set φ := φ 1 ⊗ · · · ⊗ φ k and denote ν k = ν k x1,...,x k ,t . Let u ∆ be obtained from the spectral method, with initial dataū. Let us denote (u, φ) :=´D u · φ dx. Then, as a consequence of lemma 4.3, we can write d dt where there exists L > 0 independent of ∆ and the initial dataū, such that the error term ). Taking the product over i = 1, . . . , k, we find where F (u) := u⊗u. Recognizing the special structure of the empirical measure µ ∆ t (33) as a convex combination, denoting by ν k,∆ = ν k,∆ x1,...,x k ,t the k-point correlation measure corresponding to µ ∆ t , we obtain from the above identity that, The right-hand side can be bounded bŷ

which, by lemma 4.3 is further bounded by
Note that ifμ is supported on B M (0) ⊂ L 2 x , then it follows that µ ∆ t is supported on B M (0), as well. This is a consequence of the a priori L 2 -bound (32). Hence the error term in equation (40) is in this case bounded by C∆, where C = C(φ, k, M, T ) is a constant independent of ∆.
Let us also note that the terms on the left-hand side of (40) converge strongly in L 1 t,x as ∆ → 0. Indeed, it is not difficult to see that all terms on the left -hand side, e.g. g(t, x, ξ) := (ξ 1 ⊗ · · · ⊗ F (ξ i ) ⊗ · · · ⊗ ξ k ) : ∇ xi φ(x, t), are admissible observables in the sense of (18). For such observables, the L 1 t,xconvergence of ν k,∆ t,x , g(t, x, ξ) → ν k t,x , g(t, x, ξ) , as ∆ → 0, has been established in theorem 2.12. The same holds true for the other two terms on the left-hand side.
Passing to the limit µ ∆ t → µ t , it thus follows that The fact that µ t is concentrated on incompressible vector fields follows immediately from the corresponding property of the approximations µ ∆ t (cp. lemma 3.3). Furthermore, from proposition 2.9, it also follows that the limit µ t is time-regular. This finishes the proof that µ t is a statistical solution of the incompressible Euler equations with initial dataμ. Remark 4.13. It is straightforward to show that if µ ∆ t are generated from the spectral hyper-viscosity scheme (28), and if they satisfy the assumptions of theorem 4.12, the limit µ t is in fact a dissipative statistical solution in the sense of definition 3.4.

Numerical Experiments
In this section, we will present a suite of numerical experiments to demonstrate the effectiveness of the Monte Carlo algorithm 4.1 in computing statistical solutions of the incompressible Euler equations. 5.1. Implementation. For our numerical experiments, we use the implementation of the spectral hyper-viscosity scheme (28) provided by the SPHINX code, which was first presented in [24]. In SPHINX, the non-linear advection term is implemented using the O(N 2 log N )-costly fast Fourier-transform. Aliasing is avoided via the use of a padded grid, as e.g. described in [24,23]. The spectral scheme is implemented based on the primitive variable formulation.
For the numerical experiments reported below, we use a spectral viscosity operator of order s = 1 (cp. equation (28)), with N = /N , = 1/20 unless otherwise stated. The Fourier multiplier Q N is chosen with Fourier coefficients otherwise.
corresponding to m N = √ N . For each sample, SPHINX solves the following system of ODEs for the Fourier coefficientsû k (t), |k| ∞ ≤ N of u(x, t) = |k|∞≤N u k (t)e ikx : for |k| ∞ ≤ N , k = 0. Here u k denote the Fourier coefficients of the initial data. The multiplication by the matrix (I −(k ⊗ k)/|k| 2 ) implements the Leray projection onto divergence-free vector fields. The dot product with ik · (. . .) corresponds to taking the divergence in Fourier space. We note that the 0-th Fourier component of u is constant in time, reflecting conservation of momentum. In SPHINX, this component is set equal to u 0 ≡ 0. The above system of ODEs (41) is integrated in time using an adaptive explicit third-order Runge-Kutta method. Although the theory of section 4 is valid for both two and three space dimensions and the SPHINX code is available for both cases, we restrict our focus to two space dimensions in this section, on account of affordable computational costs.

5.2.
Flat Vortex Sheet. Vortex sheets occur in many models in physics and are an important test bed for numerical experiments for the Euler equations [22] and references therein. We first consider a randomly perturbed version of the flat vortex sheet that corresponds to the following initial data also considered in [22], 5.2.1. Initial data. Given a smoothing parameter ρ > 0, and a parameter δ ≥ 0 (measuring the size of the random perturbation of the interface), this vortex sheet initial data is of the form (42) u ρ,δ (x) = P(U ρ (x 1 , x 2 + σ δ (x 1 ))), where P denotes the Leray projection, U ρ (x) = (U ρ 1 (x), U ρ 2 (x)) is the following smoothened flat vortex sheet initial data: and σ δ (x) is a random function, which for a given (random) choice of parameters α 1 , . . . , α q ∈ (0, δ), β 1 , . . . , β q ∈ [0, 2π), is defined by We will also consider the discontinuous case of initial data that are obtained in the limit ρ → 0 resulting in For our simulations, we fix q = 10 modes for the perturbations. The coefficients α k are drawn independently, uniformly in (0, 1), and then multiplied by δ. The coefficients β k are i.i.d., with a uniform distribution on [0, 2π). The initial data for the statistical solutionμ δ ρ ∈ P(L 2 x ) is defined as the law of these random perturbations. It depends on the two parameters ρ ≥ 0, δ ≥ 0. While ρ controls the smoothness of the initial data, δ measures the amplitude of the perturbation. We fix δ = 0.025 in the following and consider different values of ρ. Note that the choice ρ = 0 corresponds to an initial measure supported on discontinuous flows with a very sharp transition (see figure 2 (A,B) for realizations (samples) of this initial data). In figure 2 (C,D), we present the initial mean and variance that correspond to the random variations of the initial interface location.  Clearly when ρ > 0, the corresponding initial data for every sample is smooth. Consequently, smooth solutions of (1) are well-posed and the spectral viscosity method converges to this solution as N → ∞ [2]. However for ρ = 0, which corresponds to the case of a discontinuous vortex sheet, there are no well-posedness results even for weak solutions, as the vorticity corresponding to the initial datum (for each sample) is a sign changing measure and does not belong to the Delort Class. In [22], the authors had presented multiple numerical experiments to illustrate the approximate solutions, computed with a spectral viscosity method, may not converge (or converge too slowly to be of practical interest) for individual samples (see figures 5 and 6 of [22]). Hence, it would be interesting to study if approximate statistical solutions, generated by algorithm 4.1 converge in this case.

5.2.2.
Structure functions and Compensated Energy spectra. The convergence theorem 4.6, based on the compactness theorem 2.8, provides us with verifiable criteria to check convergence of algorithm 4.1. In particular, we need to check certain decay conditions on the structure function (11) for small correlation lengths. To this end, we consider the following instantaneous version of the structure function (11), , Note that the above is a formal definition and it can be made rigorous in terms of the time-dependent correlation measures. It is much simpler to compute the instantaneous quantity (44) than the time-averaged version (11). Our objective is to check whether the structure function (11), or rather its instantaneous version (44), decays (uniformly in resolution ∆) as r → 0. Such a decay would automatically imply convergence of the approximations to a statistical solutions by theorems 2.8 and 4.12. Clearly if ρ > 0 in (42), the spectral viscosity method converges to the unique classical solution as ∆ → 0. Moreover, a straightforward calculation shows that the structure function (44) should scale as, (45) S 2,∆ r,t (µ t ) ≈ r, ∀∆, t. This is indeed verified from figure 3 (A) where we plot the structure function (44) at t = 0.4 and ρ = 0.1 for different values of the mesh parameter. We see from this figure that S 2,∆ r,0.4 (µ t ) ≈ r 0.9 , at fine resolutions, which is very close to the expected value of 1 for the scaling exponent of the structure function.
On the other hand, for ρ = 0, corresponding to the discontinuous flat vortex sheet, the lack of smoothness inhibits us from inferring a particular form of decay of (11) (or (44)) a priori.
At least initially, calculations in [21] imply that S 2,∆ r,0 (μ) ∼ r 1 2 , for the discontinuous flat vortex sheet. Surprisingly, we find from figure 3 (B) that at fine resolutions, S 2,∆ r,0.4 (µ t ) ≈ r 0.52 , which agrees with the decay of the structure function of the initial data. Although we do not present the results there, we observe that the structure function (44) scales as r θt , with θ t ≥ 0.5 for all t. This implies an uniform decay of the structure function (11) and convergence of the approximations to a statistical solution of the Euler equations (1), even for this case of discontinuous vortex sheet data. Note that the computed structure functions (44) in figure 3 clearly satisfy the approximate scaling hypothesis (37) and thus imply convergence through theorem 4.6.  An alternative criterion for convergence of statistical solutions is provided by the energy spectrum decay in the inertial range (38). To check whether this criterion is satisfied, we follow theorem 4.9 and compute the following instantaneous compensated energy spectrum, (46) C ∆ γ,t (µ t ; K) := K γ E(µ t , K). Following the arguments in the proof of theorem 4.9, we can relate the decay of the instantaneous energy spectrum to the corresponding decay of the structure function (44) by a direct analogue of (39).
For ρ = 0.1 in (42), we plot the compensated energy spectrum C ∆ 3,0.4 (µ t ; K) for all K and at time t = 0.4, with compensating factor γ = 3 in figure 4 (A). Note that this choice of γ is consistent with a decay exponent of 1 for the structure function in (39), i.e. S 2 r (µ ∆ t ; T ) r. We observe from this figure that as expected for this case, the compensated energy spectrum is clearly bounded and in fact, decays faster than the expected rate for the entire range of wave numbers.
On the other hand, we plot the compensated energy spectrum C ∆ 2,0.4 (µ t ; K) (46) for the discontinuous flat vortex sheet case, i.e. ρ = 0 in (42), in figure 4 (B). In this case, we expect from the structure function computations (see figure 3(B)) that the instantaneous structure function decays with an exponent of ≈ 0.5. From (39), we see that this corresponds to the choice of γ = 2 as the exponent of compensation in (46). Moreover, in figure 4 (B), we also plot the line corresponding to wave number m N ≈ √ N , which for the spectral viscosity method (28) represents the wave number after which the spectral viscosity is activated and hence, demarcates the separation between inertial and dissipation ranges. We observe from figure 4 (B) that the compensated energy spectrum is clearly uniformly bounded (in terms of the resolution ∆) for the whole inertial range and for all resolutions ∆ barring the coarsest resolution, and decays fast in the dissipation range, although there is a slight kink upwards at the very end of the dissipation range, almost at the grid scale. This might be attributed to numerical errors, which are dominant at this range. Translating these results to the energy spectrum, we see that the spectrum decays as K −2 in the inertial range uniformly with respect to resolution. Hence, according to theorems 4.9 and 4.12, the sequence of approximations will converge to a statistical solution of (1).

Convergence in Wasserstein Metrics.
Given the computational results on the structure function and the compensated energy spectra, results in section 4 clearly imply convergence of the approximations µ ∆ t , generated by the Monte Carlo algorithm 4.1 to a statistical solution of the incompressible Euler equations. Moreover from the discussion in section 2, we should observe with respect to the following Cauchy rates: Unfortunately, the calculation of the Wasserstein distance between probability measured defined on high-dimensional (or indeed ∞-dimensional) spaces is a highly non-trivial issue, which we cannot tackle with present computational resources.
On the other hand, one can compute finite-dimensional marginals of (47) by utilizing the complete characterization of L 1 t (P) in terms of correlation measures as given in theorem 2.11. Following [14] (theorem 5.7), one can prove that, Here, k ≥ 1 and ν ∆,k t,x is the k-th correlation marginal corresponding to the approximate statistical solution µ ∆ t . Note that we consider instantaneous versions of the Wasserstein metric (47) for reasons of computational convenience. Figure 5. The Wasserstein distances between correlation ) dx for k = 1, 2, 3, at time t = 0.4 with respect to resolution We remark that computing the Wasserstein distances W 1 (ν ∆,k t,x , ν ∆/2,k t,x ) for small k is much more tractable. We have computed these Wasserstein distances using the algorithm of [5] (as implemented in [17]) and the corresponding results for k = 1, 2, 3, at time t = 0.4 for the discontinuous flat vortex sheet, i.e. ρ = 0 in (42) are presented in figure 5. As seen from this figure, we observe a clear convergence of these Wasserstein distances (in the Cauchy sense as in (48)) for the one-point, two-point and three-point correlation measures, albeit at a slow rate for the second and third correlation marginals. This, together with the results on the structure function and compensated energy spectra, provides considerable evidence that the approximate statistical solutions, generated by algorithm 4.1, converge to a statistical solution of (1). Moreover, given theorem 2.12, results shown in figure  5 establish convergence with respect to any admissible observable in the sense of (18), corresponding of one-point, two-point and three-point statistical quantities of interest. These include mean, variance, structure functions, energy spectra as well as three-point correlation functions.

Sinusoidal vortex sheet.
In this section, we will consider a random perturbation of the so-called sinusoidal vortex sheet, i.e. the initial vorticity is concentrated on a sine curve. This test case was extensively studied in a recent paper [23] in the context of numerical approximation of weak solutions (in Delort class) of the two-dimensional incompressible Euler equations. 5.3.1. Initial Data. We fix a sinusoidally perturbed vortex sheet, where the initial vorticity is a Borel measure of the form such that´T 2 ω 0 dx = 0, and up to a constant, ω 0 is uniformly distributed along a curve Γ, which is defined as the graph: We chose d = 0.2 for our simulations.  The numerical initial data is obtained from the mollification of this initial data with a parameter ρ > 0. As a mollifier, we consider the third-order B-spline ψ(r) := 80 7π (r + 1) 3 + − 4(r + 1/2) 3 Next, we define ψ ρ (x) = ρ −2 ψ(|x|/ρ). The numerical approximation of the perturbed vortex sheet is now defined by setting where ρ determines the thickness (smoothness) of the approximate vortex sheet. The convolution at x = (x 1 , x 2 ) ∈ T 2 is evaluated via numerical quadrature: with ξ i = x 1 + iρ/Q equidistant quadrature points in x 1 , and g(ξ) the function whose graph is Γ, i.e. g(ξ) = d sin(2πξ). We choose Q = 400 quadrature points. We denote by U ρ (x 1 , x 2 ) the velocity field such that div(U ρ ) = 0 and curl(U ρ ) = ω ρ 0 . Similar to the case of the flat vortex sheet, we carry out random perturbations of the sinusoidal vortex sheet as follows: Here, P : L 2 x → L 2 x denotes the Leray projection onto divergence-free vector fields, and we again fix a random function σ δ (x), depending on a parameter δ ≥ 0 and a choice of (random) coefficients α 1 , . . . , α q ∈ (0, δ), β 1 , . . . , β q ∈ [0, 2π). For our simulations, we fix q = 10 modes for the perturbations. In practice, the coefficients α k are first drawn independently, uniformly in (0, 1), and then multiplied by δ. The coefficients β k are i.i.d., with a uniform distribution on [0, 2π). The initial data for the statistical solutionμ δ ρ ∈ P(L 2 x ) is defined as the law of these random perturbations. It depends on the two parameters ρ ≥ 0, δ ≥ 0. While ρ controls the smoothness of the initial data, δ measures the amplitude of the perturbation. We fix δ = 0.003125 in the following and vary ρ as a function of the grid size N . To approximate vortex sheet initial data, we must scale ρ = ρ(N ) with N , such that ρ → 0 as N → ∞. We use ρ = 5/N for our simulations. The additional diffusion parameter of the spectral viscosity scheme is set to = 0.01. With this choice of parameters, we will drop the sub-and superscripts and denote the initial data at a given resolution simply by µ ∈ P(L 2 x ).    Therefore, by the results of [23], the approximate solutions generated by the spectral viscosity method (28) will converge, on increasing resolution, to a weak solution of (1). However, as noted in [23], this convergence can be very slow as the flow breaks down into smaller and smaller vortices. In fact, this phenomenon is also seen from figure 7 (Top row), where we plot the horizontal component of velocity u 1 at time t = 1.2 and different resolutions. At this time, the initial vortex sheet has rolled over and broken down into a succession of small vortices, whose location and amplitude are different for different resolutions. This very slow convergence is also displayed in figure 8 (A), where we plot the Cauchy rates u ∆ (t) − u ∆/2 (t) L 2 x , with u ∆ denoting the approximate solution computed with the spectral viscosity method (28), for three different times t = 0, 0.6, 1.2. As seen from this figure, the rate of convergence decreases very rapidly and at time t = 1.2, it appears as if there is no convergence on mesh refinement.

5.3.3.
Structure functions and Compensated energy spectra. Given this apparent non-convergence of individual samples, it is pertinent to investigate if computing the statistics will be more convergent. To this end, we consider the initial data to be the initial probability measureμ. The mean and variance (of the horizontal component u 1 ) are plotted in figure 6 (B,C). From this figure, we observe that the initial probability measure is concentrated on very small perturbations of the underlying sinusoidal vertex sheet, as reflected in the initial variance.
In order to investigate the convergence of approximations to the statistical solution, generated by the algorithm 4.1, we follow the template of the previous numerical experiment and compute the (instantaneous) structure function (44) and the compensated energy spectrum (46) in figure 9. From this figure, we observe that the structure function at time t = 1.2 scales with an exponent of ≈ 0.7 at the finest resolutions. From (39), this implies roughly a γ = 2.4 in the scaling of the energy spectrum (46). A better fit to the scaling of the energy spectrum is found with γ = 2.2. We plot the compensated energy spectrum with the latter value of γ in figure 9 (B). From this figure, we see that for the inertial range, the energy spectrum clearly decays (faster than) a rate of 2.2. Thus, the assumptions of theorems 2.8, 4.9 are satisfied and the approximations will converge to a statistical solution of (1).
(a) Instant. structure function (44) (b) Compensated energy spectrum (46) with γ = 2.2 Figure 9. The structure function and compensated energy spectrum for the sinusoidal vortex sheet at time T = 1.2

Convergence of observables and Wasserstein Distances.
Given the results on the computed structure functions and energy spectra, the approximations will converge. But is this convergence at a better rate than that of single samples? To investigate this issue, we consider two different sets of computations. First, we compute the mean and the variance of the velocity field at different resolutions and plot them (for the horizontal velocity at time t = 1.2) in figure 7 (Middle and Bottom rows). Clearly, the one-point statistics appear much more convergent than the single sample results. The mean flow consists of a coherent set of large vortices, which is in stark contrast to the large number of vortices formed in the single sample simulations. Moreover, we also plot Cauchy rates for the mean and the variance, corresponding to the norm u 2 1 + u 2 2 at time t = 1.2, and different resolutions in figure 8 (B). Again, we observe that these one-point statistics converge at a significantly faster rate than the single sample. These results indicate that one can expect significantly better convergence of approximations for statistics than for individual realizations of fluid flows, even if the initial probability measure is a small perturbation of the underlying deterministic data and further reinforces the results of [11,22,14] in this direction. )dx for the one-and two-point correlation marginals for the sinusoidal vortex sheet at time t = 1.2 Finally, we plot the Wasserstein distances (48) for k = 1, 2, corresponding to the one-and two-point correlation marginals, at time t = 1.2, in figure 10. The results clearly show convergence in these metrics at a significantly faster rate than for individual samples and indicate possible convergence in the metric (9) on probability measures on L 2 .

5.4.
Fractional Brownian Motion. The study of the evolution of initial ensembles corresponding to (fractional) Brownian motion stems from [43,45], where the authors model interesting aspects of Burgers turbulence by evolving Brownian motion initial data for the (scalar) Burgers' equation, see [13] for a more recent numerical study. Similarly in [14], the authors consider the compressible Euler equations with (fractional) Brownian motion initial data. Following these articles, we will consider the two-dimensional Euler equations (1) with initial data corresponding to fractional Brownian motion, i.e. the following initial data: where B H 1 and B H 2 are two independent two-dimensional fractional Brownian motions with the Hurst index H ∈ (0, 1). Standard Brownian motion corresponds to a Hurst index of H = 1/2. The initial probability measureμ is the law of the above random field. To generate fractional Brownian motion, we use the random midpoint displacement method originally introduced by Lévy [27] for Brownian motion, and later adapted for fractional Brownian motion, see [14] section 6.6.1.
Considering fractional Brownian motion initial data (49) is a significant deviation from the vortex sheet initial data in the following respects, • For the vortex sheet initial data, the initial measureμ ∈ P(L 2 x ) was concentrated on a 20-dimensional subset of L 2 x (corresponding to the choice of 20 free parameters α k , β k ). On the other hand, in the limit of infinite resolution (∆ → 0), the fractional Brownian motion initial data corresponds to a measure concentrated on an infinite dimensional subset of L 2 x . • For any 0 < H < 1, and for any sample ω ∈ Ω, the initial vorticity for (49) is not a Radon measure. Consequently, the initial data does not belong to the Delort class and there are no existence results for the corresponding samples. Hence, fractional Brownian motion does not fall within the ambit of any of the available well-posedness theories for two-dimensional Euler equations. • The Hurst index H in (49) controls the regularity (and also roughness) of the initial data (pathwise). Roughly speaking, each sample is Hölder continuous with exponent H. Hence, we can consider a very wide range of scenarios in terms of roughness of the initial data by varying the Hurstindex H, see figure 11 for realizations of the horizontal velocity field for three different Hurst indices. In particular, one can observe from this figure that lowering the value of H leads to oscillations of both higher amplitude and frequency in the initial velocity field.
5.4.1. Structure functions and Compensated energy spectra. In order to verify convergence of the approximations, generated by algorithm 4.1, for the fractional Brownian motion initial data (49), we will check if the computed structure functions (44) decay uniformly with respect to resolution, on decreasing correlation lengths. In figure 12 figure 12, the compensated energy spectra remain bounded up to the final time t = T , independent of the spectral resolution. Hence, the energy spectrum decays at least at the rate of K −γ for increasing wave number K, in the inertial range. Consequently, we can readily apply proposition 4.9 and conclude that the approximations, generated by the algorithm 4.1, converge to a statistical solution, for all three values of the Hurst index H in (49).

5.4.2.
Convergence in Wasserstein distance. Next, we seek to verify convergence of observables (statistical quantities of interest). To this ends, we follow the previous section and compute the Wasserstein distances´D k W 1 (ν ∆,k t,x , ν ∆/2,k t,x )dx, corresponding to the k-point correlation marginals for the three different Hurst indices of H = 0.75, 0.5, 0.15. In figure 13, these metrics are computed at time T = 1, for k = 1, 2, corresponding to one-point and two-point statistical quantities of interest. As observed from the figure, the approximations clearly converge in this metric for both one-and two-point statistics, at rates which are independent of the underlying initial Hurst index. The two-point correlation marginals appear to converge at a slower rate than the one-point Young measures. These results validate convergence of all one-and two-point statistical quantities of interest. Taken together with the results for the structure function, compensated energy spectra and theorem 2.8, they strongly suggest convergence in metric d T (9) on L 1 T (P). )dx for k = 1 (Top Row) and k = 2 (Bottom Row) for Fractional Brownian motion initial data with three different Hurst indices at time t = 1.

5.5.
Stability of the computed statistical solution. The afore-presented numerical experiments clearly validate the convergence theory developed in this article. Therefore, algorithm 4.1 provides a practical way to compute statistical solutions. We use this algorithm to study the sensitivity/stability of the computed statistical solution in the context of the discontinuous flat vortex sheet initial data (ρ = 0 in (42)) to the following perturbations, • Type of the initial perturbation, with respect to the underlying probability distribution. • Type of underlying numerical method.
• Amplitude of the initial perturbation. We start by changing the type of the underlying initial probability measure in the flat vortex sheet initial data. Instead of choosing the i.i.d. random variables α k ∼ U[−1, 1] according to a uniform distribution as before in (42), we choose them from a normal distribution α k ∼ N (0, 1/3) with mean 0 and variance 1/3. The mean and variance are chosen so the distribution's mean and variance are consistent with those of U[−1, 1].
In order to compare the corresponding approximate statistical solutions, we compute the Wasserstein distance´D k W 1 (ν ∆,k t,x ,ν ∆/2,k t,x )dx, at different resolutions ∆.
x )dx withν the correlation marginal corresponding to a Normally distributed initial measure in (42). Right Column: Distancé )dx with ν * the correlation marginal computed with the finite difference projection method.
Here,μ ∆ is the statistical solution computed with normally distributed initial data andν ∆ is the corresponding correlation measure. We set the initial perturbation amplitude δ = 0.05 in (42), t = 0.4 and k = 1, 2 and plot the computed Wasserstein distances in figure 14 (Left column). As seen from this figure, the two approximate solutions clearly converge to each other in this metric as the resolution is increased. This indicates that the computed statistical solutions are stable with respect to the variation of underlying initial probability measures.
Next, we consider if the computed statistical solution depends on the underlying numerical method. To this end, we compute the statistical solution for the discontinuous flat vortex sheet initial data (42) with algorithm 4.1 but replace the spectral viscosity method with a finite difference projection method [7,3,24]. The convergence of this method to a statistical solution is considered in a forthcoming paper [38]. We denote the computed statistical solutions with this method as µ * ,∆ t and the corresponding correlation measures as ν * ,∆ and compute the Wasserstein distance´D k W 1 (ν ∆,k t,x , ν * ,∆/2,k t,x )dx at time t = 0.4 and k = 1, 2 and plot the results in figure 14 (Right Column). From this figure, we readily conclude that the statistical solutions computed with the finite difference projection method converges to that computed with the spectral viscosity method on increasing resolution. This strongly suggests the stability of the computed statistical solutions to the underlying (convergent) numerical method.  Left). Clearly, the computed structure functions are very close and decay with decreasing correlation length with approximately the same exponent. Thus, we can appeal to theorem 2.8 and claim that the approximate statistical solutions converge. This is further reinforced when we compute the instantaneous compensated energy spectrum (39) at time T = 0.4 and with γ = 2. We observe from figure 15 (Top Row, Right) that the computed energy spectra with different values of δ are very close and decay at the approximately the same rate. Hence, from theorem 4.9, the underlying approximate statistical solutions will converge.
After establishing convergence of statistical solutions when the amplitude of perturbations in the initial datum (42) is decreased, it is natural to ask what these solutions converge to. One possibility is the initial discontinuous flat vortex sheet itself. After all, setting ρ, δ → 0 in (42), it is easy to see that the discontinuous flat vortex sheet is a stationary solution of the incompressible Euler equations. To check if this stationary solution is indeed the zero perturbation limit of the computed statistical solution, we compute the Wasserstein distances´D W 1 (ν δ,1 t,x , δū (x) ) dx and x,y , δū (x)⊗ū(y) ) dx dy for different values of the perturbation parameter δ with t = 0.4 and k = 1, 2. Here,ū is the stationary solution corresponding to the flat vortex sheet initial data, i.e. ρ, δ → 0. We display these distances in figure 15 (Bottom Row). We observe from this figure that there does not seem to be any perceptible evidence of convergence to the stationary solution. This is consistent with the findings in [22], where the authors had observed a non-atomic measure-valued solution as the limit of the perturbations to the flat vortex sheet.
In order to further identify the limit of the computed statistical solutions as δ → 0, we compute a reference solution µ ref t , by setting δ = 0.05 32 in (42) and N = m = 1024 as the spectral resolution and Monte Carlo samples. We compute the Wasserstein distances´D k W 1 (ν δ,k t,x , ν ref,k t,x ) dx for different values of the perturbation parameter δ with t = 0.4 and k = 1, 2. Here, ν ref is the correlation measure, corresponding to the reference statistical solution. The distances, plotted in figures 15 (Bottom Row), do decrease as δ is decreased, albeit at a slow rate. However, this decrease is clearly visible when compared to the lack of convergence to the stationary solution in the same figures. Thus, we have established that the computed statistical solutions are stable with respect to the amplitude of perturbations and moreover, converge to a probability measure that is different from the one concentrated on the stationary solution.

Discussion
We consider the incompressible Euler equations (1) in this paper. The existence of classical (or weak) solutions is an outstanding open question in three space dimensions. Although weak solutions are known to exist in two space dimensions, even for very rough initial data, they may not be unique. Similarly, numerical experiments reveal that standard numerical methods may not converge, or converge very slowly, to weak solutions on increasing resolution.
Given these inadequacies of traditional notions of solutions, it is imperative to find solution concepts for (1) that are well-posed and amenable to efficient numerical approximation. In this context, we consider the solution framework of statistical solutions in this paper. Statistical solutions are time-parametrized probability measures on L 2 (D; R d ). Given the characterization of probability measures on L p spaces in [12], these measures are equivalent to so-called correlation measures, i.e. Young measures on tensor-products of the underlying domain and phase space that represent multi-point spatial correlations. Furthermore, we require statistical solutions to satisfy an infinite number of PDEs (see definition 3.1) for the moments of the underlying correlation measure. Hence, a statistical solution can be interpreted as a measure-valued solution (in the sense of [10]), augmented with information about the evolution of all possible multi-point spatial correlations.
Our aim in this paper was to study the well-posedness and efficient numerical approximation of statistical solutions. To this end, first, we had to characterize convergence on a weak topology on the space L 1 t (P(L 2 (D; R d ))), under an assumption of time-regularity on the underlying measures. Convergence in this topology amounted to convergence of a very large class of observables (or statistical quantities of interest). We then proposed a notion of dissipative statistical solutions and also proved well-posedness results for them in a generic sense, namely when the initial measure is concentrated on functions sufficiently near initial data for which smooth solutions exist. This led to short-time well-posedness if the initial probability measure is concentrated on smooth functions. In two space dimensions, we proved global well-posedness for statistical solutions when the initial data is concentrated on smooth functions. Moreover, we also proved a suitable variant of weak-strong uniqueness.
Our main contribution in this paper is the proposal of an algorithm 4.1 to approximate statistical solutions of the Euler equations. This Monte Carlo type algorithm is a variant of the algorithms proposed recently in [11,22,14] and is based on an underlying spectral hyper-viscosity spatial discretization. Under verifiable hypotheses, we prove that the approximations converge in our proposed topology to a statistical solution. These hypotheses either rely on a suitable scaling (or uniform decay) for the structure function, or equivalently, on finding an inertial range (of wave numbers) on which the energy spectrum decays (uniformly in resolution). These hypotheses are very common in the extensive literature on turbulence (see [19] and references therein). A key novelty in this article is the rigorous proof of the fact that easily verifiable conditions on the structure functions or energy spectrum imply a rather strong form of convergence for (multi-point) statistical quantities of interest. For instance, we observe a surprising fact that a bound on the compensated energy spectrum (46) implies that k-point statistics of interest, even for large k, converge. The convergence results also provide a conditional global existence result for statistical solutions in both two and three space dimensions.
We present results of several numerical experiments for the two-dimensional Euler equations. From the numerical experiments, we observe that • Our convergence theory is validated by all the numerical experiments. The assumptions on the structure functions and energy spectra appear to be very clearly fulfilled in practice. Moreover, the computed solutions converge to a statistical solution in suitable Wasserstein metrics on multi-point correlation marginals. In particular, all admissible observables of interest such as mean, variance, higher moments, structure functions, spectra, multipoint correlation functions, converge on increasing resolution and sample augmentation. • In clear contrast to the deterministic case where computed solutions may converge very slowly even if one can prove convergence of the underlying numerical method (see [23] and figure 8), statistical quantities of interest seem to be better behaved and converge faster. • For our numerical examples, we observe convergence of approximations even when the initial data was quite rough such as when the initial vorticity may not have definite sign (as in the flat vortex sheet) or may not even be a Radon measure (as in the fractional Brownian motion with any Hurst index H ∈ (0, 1)). For such initial data, the samples are not in the Delort class and the convergence (and existence) theory for two-dimensional Euler equations is no longer valid. On the other hand, we find neat convergence to a statistical solution. • The computed statistical solutions where observed to be stable with respect to amplitude and type of perturbations of initial data. Moreover, we observed that two very different numerical methods converge to the same statistical solution for a fixed initial data. This is in contrast to the deterministic case, where the computed solutions can differ significantly [24].
Based on the above discussion, we conclude that statistical solutions are a promising solution framework for the incompressible Euler equations. In particular, there is some scope for proving well-posedness results within this class, possibly with further admissibility criteria. Moreover, numerical approximation of statistical solutions is feasible with ensemble averaging algorithms. Statistical solutions can be a suitable framework for uncertainty quantification and Bayesian inversion for the Euler equations and to encode and explain numerous computational and experimental results for turbulent fluid flows.
There are several limitations of the current paper, which provide directions for future work. At the theoretical level, we seek to either relax the criteria on scaling of structure functions or prove it. This will pave the way for global existence results. Similarly, the weak-strong uniqueness results of this paper could be improved.
In terms of numerical approximation, the main issue with the Monte Carlo type algorithm 4.1 is the slow convergence (in terms of number of samples). This necessitates a very high computational cost, particularly in three space dimensions. We plan to consider efficient variants such as multi-level Monte Carlo [13,37,25], Quasi-Monte Carlo [34] and deep learning algorithms [35], for computing statistical solutions of the incompressible Euler equations in three space dimensions, in forthcoming papers.
Proof. We note that for any τ ∈ [0, T 1 ], we have is precompact for any T 1 < T . Combining these results (e.g. setting T 1 = T /2), we finally conclude that Appendix C. Proof of existence and uniqueness, Thm. 3.6 We will follow the notation used in section 3.1. Before coming to the proof of theorem 3.6, we observe that In particular, G n is covered by the open balls B r(v)/n (v). Since L 2 x is a separable metric space, we can choose a countable subcovering, i.e. we can find a dense set of initial data {u i } i∈N ⊂ C in such a way that Furthermore, choosing such a countable subset {u (n) i } i∈N for each G n , n ∈ N, and considering finally the union i } i∈N , we can in fact find a single countable subset of C (not depending on n), such that (54) holds for all G n , n ∈ N.
After this preliminary observation, we now come to the proof of theorem 3.6.
Proof of theorem 3.6. Letμ ∈ P(L 2 x ) be given, such thatμ(G) = 1. The idea of the proof is to first choose a countable set {u i } i∈N satisfying (54) for all n ∈ N, and to observe thatμ can be well approximated by atomic measures of the form ρ = ∞ i=1 α i δ ui . We then use the regularity of the solution t → u i (t) of the Euler equations with initial data u i to show that a dissipative solution µ t exists and that it is unique.
Construction of suitable approximants: For any n ∈ N, let us first construct a suitable approximantρ (n) ≈μ, of the form We define the coefficients α and inductively Note that S (n) i can be thought of as the support ofμ close to u i , and Σ (n) i keeps track of the set that has already been assigned to u j for some j < i. By this construction, we have and furthermore A dissipative statistical solution with initial dataρ (n) is given by In the following, we will show that (ρ (n) t ) n=1,2,... is Cauchy for any t ∈ [0, T ), implying the existence of a limit ρ (n) t µ t . Since the definition of a dissipative statistical solutions is linear in the probability measure, it is not hard to see that such a limit µ t must itself be a dissipative statistical solution. Furthermore, we will show that the limit µ t hasμ as initial data. Finally, we will show that this µ t is in fact unique in the class of dissipative statistical solutions.
The sequence ρ (n) t is Cauchy: Let m, n ∈ N be arbitrary. Note that implies that upon defining S we can define a transport plan π t ∈ P(L 2 where we set α 0 = i>N α i , and u 0 (t) ≡ 0. Here, N is chosen sufficiently large to ensure that where M ≥ 1 is chosen such thatμ(B c M ) = 0. Such M exists by assumption onμ. Note in particular, that for this choice of ρ t , we have Define nowμ i ∈ Λ(α,μ) bȳ x dµ i (u). Repeating the argument used in the proof of weak-strong uniqueness for measurevalued solutions [6], it follows from the fact that u i (t) is a strong solution, that by assumption on the L 2 -boundedness of µ.