On the long time behavior of compressible fluid flows excited by random forcing

We are concerned with the long time behavior of the stochastic Navier--Stokes system for compressible fluids in dimension two and three. In this setting, the part of the phase space occupied by the solution depends sensitively on the choice of the initial state. Our main results are threefold. (i) The kinetic energy of a solution is universally and asymptotically bounded, independent of the initial datum. (ii) Time shifts of a solution with initially controlled energy are asymptotically compact and generate an entire solution defined for all $t\in R$. (iii) Every solution with initially controlled energy generates a stationary solution and even an ergodic stationary solution on the closure of the convex hull of its $\omega$--limit set.


Introduction
The role of random forcing incorporated in originally deterministic models is in many cases to substitute for the effect of external driving mechanism represented by inhomogeneous boundary conditions.This gives rise to a mathematically simpler model that should retain, however, the essential behavior of the original system at least in the long run.In the framework of continuum fluid mechanics, equations with random forcing should shed some light on more complex problems related to turbulence.In particular, the celebrated ergodic hypothesis asserts: Time averages along trajectories of the flow converge, for large enough times, to an ensemble average given by a certain probability measure.There is a common belief that such measure is in fact unique and completely characterizes the behavior of the fluid in the long run.This is supported by the pieces of evidence in the case of incompressible flows driven by the Navier-Stokes system with an additive stochastic forcing.More precisely, for the incompressible planar flow driven by a physically relevant very degenerate stochastic forcing, unique ergodicity was established by Hairer and Mattingly [22].The absence of a similar result in the physically relevant 3-D case is due to the existing gaps in the mathematical theory, in particular stability (uniqueness) of solutions with respect to the initial data.Nevertheless, it has been proved that noise has a beneficial impact when it comes to long time behavior and ergodicity.Da Prato and Debussche [8] obtained a unique ergodicity for 3-D stochastic incompressible Navier-Stokes equations with non-degenerate noise.The theory of Markov selections by Flandoli and Romito [18] provides an alternative approach which also allowed to prove ergodicity for every Markov solution, see Romito [29].The concept of (statistically) stationary solutions has been introduced both in the deterministic [19], [33] and stochastic framework [17].
The question of identifying a unique invariant measure for compressible fluid motions excited by random forces is substantially different from the incompressible setting.As a matter of fact, this possibility is naturally limited/excluded as there are certain invariant quantities, for instance the total mass, that persist under the action of random forcing.Accordingly, the part of the phase space occupied by the trajectories necessarily depends sensitively on the choice of the initial state.It is therefore desirable to show that there exist invariant measures/stationary solutions generated by solutions to the initial value problem supported on a suitably defined ω−limit set.
The main goal of the present paper is to investigate this question in the setting of the compressible Navier-Stokes system under stochastic perturbations.Roughly speaking, the result can be shown by means of the standard Krylov-Bogoliubov theory as long as the solutions of the stochastic problem enjoy the following properties: • Global existence.The problem admits a global-in-time solution -a random process ranging a suitable phase space -for any sample of initial data.• Global boundedness.The expected value of a suitable norm of global-in-time solutions is bounded independently of time.• Asymptotic compactness.The law of any global-in-time solution is tight in a suitable space of trajectories.The above outlined points also summarize the strategy of our proof.To be more precise, let ̺ = ̺(t, x) denote the mass density and u = u(t, x) the bulk velocity of a compressible viscous fluid occupying a bounded physical domain Q ⊂ R d , d = 2, 3.In this paper, we are concerned with the compressible Navier-Stokes system driven by a stochastic forcing: (1. 1) d̺ + div x (̺u) dt = 0, (1.2) d(̺u) + div x (̺u ⊗ u) dt + ∇ x p(̺) dt = div x S(∇ x u) dt + ̺g(̺, u)dt + ̺F(̺, u)dW, where we include a deterministic force g as well as a stochastic force driven by a Wiener process W .The problem is closed by imposing the no-slip boundary condition (1.4) u| ∂Q = 0.
We refer to [1,2] for the existing mathematical theory of the problem (1.1)-(1.4).The long-time behavior of global-in-time solutions of the deterministic problem was studied in [12], [27], [28], and the monograph [15].In particular, the problem may admit several (a continuum of) stationary solutions already in the deterministic setting, see [13].In [10], it was proved that every bounded solution to the deterministic system gives raise to a statistical stationary solution supported on its ω−limit set.The existence of stochastically stationary solutions to (1.1)-(1.4)with a given total mass was established in [3].
In the present paper, we go beyond the result of [3].In particular, we focus on a physically relevant hard sphere pressure-density equation of state (see Section 2.1) and show that solutions with initially controlled energy remain universally and asymptotically bounded in expectation, independently of the initial condition, see Theorem 2.6.In other words, all such solutions ultimately enter a bounded absorbing set.Furthermore, we establish an asymptotic compactness of solutions which in particular applies to sequences of time shifts of solutions and gives raise to an entire solution defined for all t ∈ R, see Theorem 2.7.Finally, we deduce that any solution with initially controlled energy generates a stationary solution, see Theorem 2.5.
Unlike in the deterministic setting (cf.[10]), we define an ω-limit set as a set of probability laws, not a set of trajectories.We then show that for any solution with initially controlled energy there is a stationary solution whose law belongs to the closure of the convex hull of its ω-limit set, see Corollary 6.1.Moreover, the method is constructive -the solution is obtained by a direct application of Krylov-Bogoliubov's method applied on the ω-limit set.Our approach is motivated by the pioneering work of Itô and Nisio [23] and the idea of Sell [31] replacing the natural phase space by the space of trajectories, see also Romito [30].Finally, we prove that there is an ergodic stationary solution on the closure of the convex hull of every such ω−limit set, see Theorem 6.3.In the setting of [3], the existence of the ω−limit set cannot even be proved and the procedure of the present paper cannot be repeated there.
As a byproduct of our strategy, we deduce two new results for the stochastic Navier-Stokes system (1.1)-(1.4)that are of independent interest themselves: • Bounded moments of the total energy.
is the energy of the fluid.The constants E ∞ (m), m = 1, 2, . . .are universal and independent of the initial condition.• Asymptotic compactness.The law of the time shifts of a fixed solution is tight in a suitable trajectory space.In the context of stochastic incompressible Navier-Stokes system, the analogous results are basically immediate due to the dissipative nature of the problem and the good compactness properties of solutions.However, the compressible Navier-Stokes system is a mixed parabolic-hyperbolic system with a very delicate structure.Its incompressible counterpart instead is a semilinear parabolic system, to some extend rather similar to the heat equation.As a consequence, the global boundedness as well as the asymptotic compactness both become substantially more difficult in the compressible case.The key difficulty is, on the one hand, the lack of energy dissipation stabilizing the system in the long run, on the other hand, the fact that the available energy and pressure estimates do not directly lead to strong convergence of the density necessary in order to pass to the limit.
The first issue is overcome by performing a higher order energy estimate together with a new dissipation balance estimate and proving the existence of a bounded absorbing set in expectation, see Section 3.2.The solution of the second issue leans on establishing the strong convergence of the approximate densities despite the fact that the initial conditions are lost in the limit process, see Section 4.4.This is a delicate issue and requires careful analysis of the oscillation "damping" in the renormalized equation of continuity.
The paper is organized as follows.In Section 2, we recall the basic definitions, formulate the hypotheses and state the main results.The global-in-time estimates are established in Section 3. Section 4 and Section 5 form the heart of the paper.In Section 4, we show tightness of the time shifts of global trajectories.This property is subsequently used in Section 5, where the Krylov-Bogoliubov method is applied to obtain a stationary solution.Section 6 is devoted to the study of the ergodic structure of the set of stationary solutions.A sketch of the proof of existence of global-in-time solutions is given in the Appendix.

Mathematical framework and main results
2.1.Pressure-density equation of state.The uniform bounds on the total energy require strong control of the fluid density.To this end, we consider the physically relevant hard sphere pressure-density equation of state.Specifically, there is a limit density ̺ > 0 such that for some exponents γ > 1, β > 3. (2.1) The restriction on γ and β are technical and can be possibly relaxed.The essential feature of (2.1) is the singularity of the pressure at ̺ yielding the (deterministic) bound ̺ ≤ ̺.This hypothesis is relevant for any real fluid, see e.g.Carnahan and Starling [6], Kolafa et al. [24].
2.2.Driving force.The deterministic driving force is given by g Let T ≥ 0 and let (Ω, F, (F t ) t≥−T , P) be a complete probability space with a complete right-continuous filtration (F t ) t≥−T .The stochastic process W is a cylindrical (F t )-Wiener process on a separable Hilbert space U normalized so that W (0) = 0.It is formally given by the expansion W (t) = ∞ k=1 e k W k (t) where (W k ) k∈N is a sequence of mutually independent real-valued Wiener processes relative to (F t ) t≥−T normalized so that W k (0) = 0, and (e k ) k∈N is a complete orthonormal system in U. Accordingly, the diffusion coefficient F is defined as a superposition operator F(̺, u) :

The coefficients
uniformly in x ∈ Q.Finally, we define the auxiliary space U 0 ⊃ U via endowed with the norm Note that the embedding U ֒→ U 0 is Hilbert-Schmidt.Moreover, trajectories of W are P-a.s. in C loc ([−T, ∞); U 0 ).(c) the density 0 ≤ ̺ ≤ ̺ belongs to the space C weak,loc ([−T, ∞); L q (Q)) P-a.s. for any 1 ≤ q < ∞ and is (F t )-adapted; (d) the momentum ̺u belongs to the space C weak,loc ([−T, ∞); L2 (Q, R d )) P-a.s. and is (F t )-adapted; (e) the velocity u belongs to L 2 loc ([−T, ∞); W1,2 0 (Q, R d )) P-a.s. and is (F t )-adapted 1 ; (f) the total energy for all −T ≤ τ 1 < τ 2 , ψ ∈ C 1 (Q), P-a.s.; (i) the momentum equation holds for all φ ∈ C 1 c ((−T, ∞)), φ ≥ 0, P-a.s.We recall that the total energy is a lower semi-continuous function defined for any t ∈ [−T, ∞) P-a.s. for any dissipative martingale solution (̺, u).In addition, it follows from the energy inequality (2.6) that the limit (2.7) is well defined and F t −measurable for any t ≥ −T .Moreover, the function E is càdlàg in [0, ∞), thus progressively measurable, and (2.8) In the remaining part of the paper, we use E to denote the total energy keeping in mind (2.8).

Stationary solutions.
The concept of stationary solution is motivated by the approach of Itô and Nisio [23].We first introduce the space of trajectories.Despite the fact that the (weakly) time continuous quantities are the conservative variables ̺, m = ̺u, it is more convenient to consider the standard variables where C loc,0 denotes the space of continuous functions vanishing at 0. We define the trajectory space Next, we introduce the time shift operator, It is easy to check that the time shift

Definition 2.3 (Stationary solution)
. We say that an entire solution

Main results.
Having collected all the necessary material we are ready to state our main results.
Then there is a sequence T n → ∞ and a stationary solution The proof of Theorem 2.5 leans on two auxiliary results that are of independent interest.Theorem 2.6 (Ultimate boundedness).Under the hypotheses of Theorem 2.5, there exists a universal constant E ∞ (m) such that for any dissipative martingale solution for all t > 0.
In the following result, we employ the notation E n for the càdlàg version of the energy associated to (̺ n , ̺ n u n ) and defined through (2.7).
Then there exists a subsequence (not relabeled) such that where (̺, u, W ) is an entire solution of the same problem defined on a certain probability space (Ω, F, P) with a filtration (F t ) t∈R .
The energy momenta estimates claimed in Theorem 2.6 are new in the context of the stochastic problem and depend essentially on the properties of the hard-sphere pressure equation of state in (2.1).The crucial point is to control the density (pressure) in terms of the dissipation term With these estimates at hand, the proof of Theorem 2.7 follows the steps of the proof of existence.There is, however, a key difference, namely, the initial data represented by the value of the time shifts at the time −T n are "lost" in the limit process.In particular, the crucial ingredient of the existence proof -compactness of the initial data -is no longer available.Instead, the deterministic argument on propagation of density oscillations proved originally in [14] must be adapted to the stochastic framework.

Global in time estimates
Our goal is to show the momentum estimates claimed in Theorem 2.6.

3.1.
A Gronwall-type estimate for BV-functions.Let us start with a standard version of Gonwalls' lemma.Lemma 3.1.Let F ∈ BV loc (0, ∞) be such that for any 0 < τ 1 ≤ τ 2 with some D > 0. Then for all t > 0 Proof.we first note that (3.1) is equivalent to It follows that the function is non-increasing on (0, ∞); whence 3.2.Higher energy moments.Let us introduce the so-called Bogovskii operator B enjoying the following properties: see Galdi [20, Chapter 3], Geißert, Heck, and Hieber [21].Here n denotes the outer unit normal on ∂Q.
In what follows, we neglect the deterministic forcing ̺g(̺, ̺u) as its treatment does not present any additional difficulties.Recalling our convention (2.8) we deduce from the energy inequality (2.6) for any 0 ≤ τ 1 < τ 2 P-a.s.
Next, we use the quantity for any 0 ≤ τ 1 < τ 2 P-a.s., where We combine the energy inequality (3.3) with the pressure estimates (3.4) to obtain the total dissipation balance that is a crucial tool in the subsequent analysis.First, consider β(̺) = ̺ in (3.4) to gain P-a.s. where Next, we recall the uniform bounds on the density and hypothesis (2.9) implying for some deterministic constant δ > 0. Setting Consequently, it follows from (3.5) P-a.s.
where c is a deterministic constant.
Next, we are going to estimate the second, third and fourth term on the right-hand side of (3.7) by dissipation.By virtue of the Korn-Poincaré inequality, and the uniform bound (3.6) we obtain (3.9) Similarly, the Korn-Sobolev inequality yields (3.10) Of course, the estimate holds for arbitrary exponent q if d = 2. Consequently, we have P-a.s.
using properties of the Bogovskii operator, see (3.2), and boundedness of ̺.Plugging this into (3.7)we conclude that P-a.s.
with deterministic constants c 1 , c 2 .Multiplying (3.11) by a sufficiently small (deterministic) constant ε > 0 and adding the resulting expression to the energy inequality (3.3) we obtain the dissipation balance P-a.s.
where C is a deterministic constant, and where we have used with a deterministic constant c, cf.(2.3).
As the next step we aim at deriving higher moment estimates for the process in the spirit of Itô's formula applied to (3.12).For τ 1 > 0 arbitrary, we consider the continuous process In accordance with the dissipation balance (3.12), the process In order to apply Itô's formula to the process D = D M +D C it suffices to regularise D M in time.To avoid problems with progressive measurability we introduce the backward regularization of a function F = F (t) given by where χ κ is a standard family of regularizing kernels.Applying the standard Itô formula to for any Φ ∈ C 2 .If, in addition, Φ ′ ≥ 0, we have M is non-increasing.All other terms in (3.14) converge to their counterparts as κ → 0 and we obtain Φ (D) for all Φ ∈ C 2 , Φ ′ ≥ 0. We clearly have using (3.9) as well as for an arbitrary κ ∈ (0, 1) using (2.3), (3.6) and continuity of B. Plugging these estimates into (3.15) and applying expectations yields for all integers m > 0 with some positive constants C m , D m .Here the passage from τ 2 − to τ 2 + follows from the fact that D is a sum of a non-increasing function and a continuous one, specifically, D(τ 2 +) ≤ D(τ 2 −).Also note that we approximated the mapping D → |D| m by a sequence of smooth functions Φ with bounded derivatives.
Thus Lemma 3.1 together with (3.17) give rise to the uniform bound In view of the bound which holds true for all τ ≥ 0, we deduce that (2.10) follows.Consequently, the first claim of Theorem 2.6 holds and the proof is complete.
As a consequence, we may also control the supremum over time inside expectation.
Corollary 3.2.Let T > 0. Under the assumptions of Theorem 2.6 it holds where the function T → c T > 0 is locally bounded on [0, ∞).
Proof.We consider (3.15) and take first supremum over time and then expectation.In order to estimate the stochastic integral, we apply Burkholder-Davis-Gundy's inequality, (3.16) and Young's inequality to obtain The first term on the above right hand side can be absorbed into the left hand side of the estimate.The second term on the right hand side is controlled in view of (3.17) by the initial value.Altogether, we deduce T , which yields the claim by using (3.18).

Asymptotic compactness
Our next goal is to prove Theorem 2.7.
, where the implicit constant as well as c 2T > 0 is universal, i.e. independent of the solution [̺ n , u n , W n ].Recall that [̺ n , u n , W n ] solves the system on [−T n , ∞) and that T n → ∞.Hence, employing Theorem 2.6, in particular (2.10), and the uniform bound (2.11), we conclude that if n is sufficiently large so that −T + T n > M for some M > 0 large enough, then which proves (4.1) for n ≥ n(T ) large enough.On the other hand, for n < n(T ) and T n > T the left hand side of (4.1) is bounded using Corollary 3.2 by , where the last inequality follows from the fact that T → c T is locally bounded on [0, ∞) and −T + T n ≤ M.
The bound (4.1) for the case T n < T follows from (2.11) and Remark 2.4.
where we have set In view of hypothesis (2.1), the pressure potential satisfies Consequently, if ω > 0 is chosen small enough, all integrals on the right-hand side of (4.2) except I n are controlled by the energy bounds (4.1).
As for I n , the smoothing properties of the operator B specified in (3.2) can be used to control I n by the energy as long as the quantity 2 which is possible as β > 3. Passing to expectations in (4.2) we may therefore infer that (4.4) Note that, in view of hypothesis (2.1),

4.3.
Limit process.The energy estimate (4.1) and the pressure estimate (4.4) are exactly the same as those obtained in the existence theory.Following the stochastic compactness arguments of [1,Chapter 4] or rather [1,Chapter 7] which also gives the necessary additional details regarding the trajectory space T , we may use Jakubowski-Skorokhod's representation theorem and find a new sequence of random variables ̺n , ũn , with associated cylindrical Wiener processes Wn defined on the standard probability space (Ω, F, P) = ([0, 1], B, dy) such that (up to a subsequence) for any 1 ≤ q < ∞.In addition, there exists a process [̺, u, W ] such that for any T > 0 P−a.s.In particular,

Asymptotic compactness of densities.
To finish the proof of Theorem 2.7, it remains to show that [̺, u, W ] is an entire solution.This can be done similarly to [1,Section 4.5] as soon as we are able to show strong convergence of the density sequence {̺ n } ∞ n=1 P−a.s.This is a delicate issue as we have no information on compactness of "initial data".
First observe that (4.5) yields the equation of continuity for the limit functions, namely we may use the standard regularization technique of DiPerna and Lions [9] to deduce the renormalized version of (4.6), To be able to let n → ∞ in (4.8), we must extend the convergence stated in (4.5) to nonlinear functions of (̺, u, ∇ x u).This is possible as the Skorokhod argument can be extended to any composition (cf.[1, Proposition 4.5.5]) as long as (4.9) uniformly for n → ∞ for some r > 1.Consequently, we may assume, in addition to (4.5) that (4.10) for any T > 0 P-a.s. as soon as (4.9) holds.In particular, we may let n → ∞ in (4.8) to obtain s. Subtracting (4.7) from (4.11) and using spatially homogeneous test functions ϕ yields an ODE for the oscillation defect The existence theory for the compressible Navier-Stokes system leans on Lions identity see Lions [26].Validity of (4.13) has been extended to sequences of (approximate) solutions in [1,Section 4.5].Plugging (4.13) in (4.12) yields P-a.s.
whence the defect D is a non-increasing function of time.This immediately yields the desired conclusion D ≡ 0 as soon as we know that D(t 0 ) = 0 for some t 0 , which is for instance the case for solutions of the initial-value problem emanating from a compact sequence of initial data.
On the other hand, as the pressure satisfies (2.1), Consequently, we deduce from (4.
≤ 0 for a.a.t ∈ R for some θ > 0. Since 0 ≤ D ≤ D for any t ∈ R, we obtain the desired conclusion D = 0 yielding strong L 1 −convergence of {̺ n } ∞ n=1 P-a.s.We have proved Theorem 2.7.

Construction of stationary solutions
The goal of this section is to prove Theorem 2.5.Therefore, let [̺, u, W ] be a dissipative martingale solution on [0, ∞) defined on some stochastic basis (Ω, F, (F t ) t≥0 , P) and satisfying (2.9).We define the probability measures (5.1) More precisely, the time average is defined as a narrow limit of Riemann sums, i.e. for every F ∈ BC(T ) we have for a sequence of equidistant partitions As explained in Remark 2.4, also in (5.1) we tacitly regard functions defined on time intervals [−t, ∞) as trajectories on R by extending them to s ≤ −t by the value at −t.
The proof of Theorem 2.5 now proceeds in two main steps.First, we prove tightness of the above measures and apply Prokhorov's theorem in order to obtain a narrowly converging subsequence.Second, in view of Theorem 2.7 we identify the limit measure as a law of a stationary solution.
Proposition 5.1.The family of measures {ν S ; S > 0} is tight in T .
Proof.Choose a bounded interval [−T, T ].In order to prove tightness of ν S , we first prove tightness of the laws of the time shifts In view of (2.9), we may apply the considerations of Section 4.1 applied to the time shifts S t [̺, u, W ] to deduce The important point is that the bound depends on the length of the time interval but not on the time shift.As a consequence, the time shifts u(• + t) are tight on equipped with the weak topology.Moreover, from the continuity equation we get for all t ≥ 0 This implies tightness of ̺(• + t) on C weak,loc (R; L q (Q)), 1 ≤ q < ∞, using also the boundedness of ̺, .
This already implies tightness of the projection of the time averaged measures ν S to the first two components.Indeed, if ε > 0 is given and K ε is the associated compact set in C weak,loc (R; L q (Q)), such that The argument for the projection to u is the same.
In addition, we recall that the shift chosen on the noise is S t W = W (• + t) − W (t). Consequently, every S t W is a Wiener process with S t W (0) = 0.That means that all S t W have the same law which is tight on C loc,0 (R, U 0 ).Altoghether, the claim follows.
As the next step, we observe that limits of the ergodic averages are invariant under various shifts.This in particular implies shift invariance of any accumulation point of the time averages ν S as we will see below.
Lemma 5.2.It holds Proof.Let G ∈ BC(T ).Using the continuity of the time shifts t → S t on T , we have Using boundedness of G, the above has the same narrow asymptotic limit as which finishes he proof.
As a consequence, we observe that if the narrow limit of in P(T ) as n → ∞ exists for some τ = τ 0 ∈ R then it exists for all τ ∈ R and is independent of the choice of τ .In view of Proposition 5.1 and Lemma 5.2 together with Jakubowski-Skorokhod's theorem, there exists a sequence S n → ∞ and ν ∈ P(T ) so that ν 0,Sn → ν narrowly in P(T ) as well as ν τ,Sn → ν narrowly for all τ ∈ R. Accordingly, the limit measure ν is shift invariant in the sense that for every G ∈ BC(T ) and every τ ∈ R we have To conclude the proof of Theorem 2.5, it remains to show that ν is a law of an entire solution to (1.1)- (1.4) in the sense of Definition 2.2.We begin with an auxiliary proposition.
By [1, Lemma 2.1.35],W is a cylindrical Wiener process with respect to its canonical filtration and W (0) = 0.As the next step, we want to strengthen this and show that W is non-anticipative with respect to the joint filtration (F t ) t≥−T , which in view of [1,Corollary 2.1.36]implies that it is a cylindrical Wiener process with respect to (F t ) t≥−T , as required in Definition 2.1.To this end, we observe that for any F : T → R bounded Borel we have This is due to the fact that the time average is defined as a narrow limit of Riemann sums and the extension to bounded Borel functions follows by the dominated convergence theorem.
We know that for every t ≥ 0 the joint canonical filtration generated by S t [̺, u, W ] is non-anticipative with respect to S t W in the sense that for every s ≥ −T , every bounded continuous functions h 1 : T | (−T,s] → R and h 2 : U 0 → R and every τ ≥ 0. Therefore, the integrand in (5.4) can be written as a composition of a bounded continuous function F : T → R with S t [̺, u, W ], where F does not depend on t ≥ 0. Using this function in (5.3) we obtain where we slightly abused the notation: each expected value E possibly refers to a different probability measure as the processes [̺, ũ, W ] and [̺, u, W ] can be defined on different probability spaces.Thus, W is non-anticipative with respect to the joint filtration where G : R → [0, ∞) is continuous, bounded, G > 0 on R\{0} and G(0) = 0. Plugging this into (5.3)we deduce that Since [̺, u, W ] is a solution on (−T, ∞), it follows that S t [̺, u, W ] is a solution on (−T −t, ∞) and, in particular, the continuity equation holds on (−T, ∞).Consequently, the integrand on the right-hand side of (5.5) vanishes for all t ≥ 0. This implies that the continuity equation is also satisfied by [̺, ũ, W ] on (−T, ∞).The same argument applies to the renormalized continuity equation (2.4).
For the momentum equation (2.5) as well as for the energy inequality (2.6) we need to proceed differently since the corresponding stochastic integrals are generally not defined as functions on the space of trajectories T .Recall that the momentum equation (2.5) is solved on (−T, ∞) by S t [̺, u, W ] for all t ≥ 0. Furthermore, as it was showed for instance in the proof [1, Theorem 2.9.1], stochastic Itô integrals of the form T S G(̺, u) dW can be written as a composition H(̺, u, W ) : Ω → R where H : T → R is a measurable function which is universal in the sense that it depends on S, T but is independent of the process (̺, u, W ) provided W is a cylindrical Wiener process with respect to some filtration and G(̺, u) is stochastically integrable with respect to W .As a consequence, also for the momentum equation, there is a bounded Borel function F : T → R such that F (̺, u, W ) = 0 P-a.s. if and only if [̺, u, W ] satisfies (2.5).This function can now be used in (5.5) to deduce that [̺, ũ, W ] satisfies (2.5) on (−T, ∞).
A similar argument can be applied for the energy inequality as well.More precisely, we put all the terms in the energy inequality on the left hand side and write the left hand side as a composition H(̺, u, W ) for some Borel function H : T → R. Then we define F = G • H where G : R → R is bounded and continuous such that G(z) = z + for z ≤ 1. Applying (5.5) we finally conclude that the energy inequality (2.6) is satisfied by [̺, ũ, W ]. The remaining points of Definition 2.1 are immediate and hence the proof is complete.
We recall that the probability measure ν was obtained as a narrow limit of the time averages ν τ,Sn for any τ ∈ R and a sequence S n → ∞.Since ν is shift invariant, any process with law ν is stationary as required in Definition 2.3.To conclude the proof of Theorem 2.5, it remains to prove that ν is a law of an entire solution to (1.1)-(1.4) in the sense of Definition 2.2.
Recall that [̺, u, W ] from the statement of Theorem 2.5 solves the system on [0, ∞).
Hence in particular the mth moment of the energy of S t+τm [̺, u, W ] at time s = −τ m is controlled this way.Since the law of [̺ m , ũm , Wm ] is given by ν τm,S n(m) −τm , it follows which yields (2.11).More precisely, this follows by applying (5.3) to a bounded truncation of the energy and then passing to the limit.Thus, by Theorem 2.7, [̺, ũ, W ] is an entire solution to (1.1)- (1.4).This completes the proof of Theorem 2.5.

Ergodic structure
In this section we study the ergodic structure of the system (1.1)-(1.4).In particular, we show that each dissipative martingale solution as in Theorem 2.5 gives raise to an ergodic stationary solution on the closure of its limit set.This is the best we can say at the moment, as we generally do not expect stationary solutions to (1.1)-(1.4) to be unique.As a matter of fact, already the deterministic counterpart of (1.1)-(1.4)may admit infinitely many equilibrium states for a given total mass.
For a dissipative martingale solution [̺, u, W ] satisfying (2.9), we define the ω-limit set as a subset of P(T ) given by Ξ[̺, u, W ] = L T [r, w, B]; there exists T n → ∞ so that S Tn [̺, u, W ] → [r, w, B] in law in T .
In addition, according to Theorem 2.6 and Theorem 2.7, the ω−limit set Ξ[̺, u, W ] is a non-empty set of laws of globally bounded entire solutions, which is shift invariant and compact.Moreover, we observe that similarly to the proof of Proposition 5.3 (see also [1, Theorem 2.9.1]), it holds that every process [r, w, B] having the law of an entire solution is an entire solution itself.
Let co(Ξ[̺, u, W ]) denote the closure of the convex hull of Ξ[̺, u, W ] with respect to the narrow convergence of probability measures.Theorem 2.5 then implies the following.Corollary 6.1.For every dissipative martingale solution [̺, u, W ] as in Theorem 2.5, there exists a stationary solution whose law belongs to co(Ξ[̺, u, W ]).
Proof.As discussed above, the ω−limit set Ξ[̺, u, W ] consists of laws of globally bounded entire solutions.Let [r, w, B] be a process whose law belongs to Ξ[̺, u, W ].
for some deterministic constants ̺ min , ̺ max , and (Ω, F, (F t ) t≥0 , P α ), ̺ α , u α , W ) to (1.1)-(1.2) in the sense of Definition 2.1 with the pressure p α follows from [1, Thm.4.0.2.] (see also [4,Thm. 2.4]).Although only a pressure of the form a̺ γ with a > 0 is treated in [1], it is clear that the same arguments apply for any monotone pressure function which behaves asymptotically as ̺ γ .Also note that [1] deals only with periodic boundary conditions.However, as demonstrated in [5] for the full Navier-Stokes-Fourier system the approach also applies to the case of bounded domains with Dirichlet boundary conditions for the velocity field (see also [32]).From the energy inequality (2.6) we obtain uniformly in α for all T > 0 and all n = 1, . . ., r 0 .Next we aim at establishing uniform bounds for the pressure.As in (3.4) we obtain for any T > 0 with A α = 1 |Q| Q (β(̺ α ) − β ′ (̺ α )̺ α ) div x u α dx.After taking expectations the stochastic integral vanishes whereas all the other terms can be estimated as in the deterministic case based on (A.2) for the choice β(̺) = (̺ − ̺) −ω with ω > 0 small.Recalling the arguments from Section 4.2, in particular (4.3), we can bound all terms on the right-hand side and obtain consequently for all T > 0 using also (A.2) for n = 1.
With estimates (A.2) and (A.4) at hand one can apply the stochastic compactness method based on the Jakubowski-Skorokhod representation theorem exactly as in [1,Chapter 4.4].Also, we can pass to the limit in all terms in the equations and the energy inequality apart from the pressure arguing as in [1,Chapter 4.4].Note that, since (A.4) implies higher integrability of the pressure, the method from [25] applies directly as explained in [16,Sec. 3.6].This is in fact a purely deterministic argument and the