A Hamiltonian interacting particle system for compressible flow

The decomposition of the energy of a compressible fluid parcel into slow (deterministic) and fast (stochastic) components is interpreted as a stochastic Hamiltonian interacting particle system (HIPS). It is shown that the McKean-Vlasov equation associated to the mean field limit yields the barotropic Navier-Stokes equation with density dependent viscosity. Capillary forces can also be treated by this approach. Due to the Hamiltonian structure the mean field system satisfies a Kelvin circulation theorem along stochastic Lagrangian paths.


1.
A. The barotropic Navier-Stokes equations. Consider a compressible barotropic fluid in an n-dimensional domain with periodic boundary conditions. The velocity, u = u(t, x), and density, ρ = ρ(t, x), are a time-dependent vector field and function, respectively, defined on the torus M = R n /Z n .
For background regarding the barotropic Navier-Stokes equations with viscosities which depend linearly on the density we refer to [43] and references therein.
The capillary tensor (1.5) appears in this form also in [3,Equ. (4)] and goes back to Korteweg [37]. Analytic aspects of the barotropic system with capillary forces (1.1)-(1.2) are treated in [10]. Navier-Stokes equations with more general third-order spatial derivative terms are discussed in [36].
1.B. Description of results. This paper is concerned with a mean field representation of solutions to the compressible Navier-Stokes system (1.1)- (1.2), and this mean field is derived from a stochastic Hamiltonian interacting particle system (Hamiltomian IPS or HIPS). The HIPS picture follows from a combination of a many particle approach to fluid dynamics and a decomposition of the energy into slow (deterministic) and fast (stochastic) components. The interpretation as a slow-fast decomposition is consistent with the multi-time formulation of [15], where it is shown that advection along stochastic transport fields in the Eulerian representation of stochastic fluid dynamics can be obtained by homogenization.
To describe the many particle approach, consider a tiny Eulerian volume (fixed in space and time) ∆V which is divided into a very large number N of equal subvolumes. It is assumed that the continuum hypothesis holds in each of the infinitesimal subvolumes ∆V α . Thus there is a mass density ρ α for each grid index α. Since the subvolumes are equal, so are the initial conditions for ρ α , and ρ α | 0 = ρ| 0 which is the initial condition for the overall mass density in ∆V . Now, the fluid parcels in all of the subvolumes interact because the deterministic component of the energy of the blob of fluid in ∆V depends on the total momentum and the overall mass density. ('Momentum' shall always refer to momentum per unit volume, i.e. its dimension is density times velocity.) On the other hand, molecular diffusion is modeled as a set of N independent (multidimensional) Brownian motions such that all N individual parcels undergo their own stochastic process. Since molecules are incompressible these processes are set up as stochastic perturbations along divergence free vector fields.
Hence the fluid parcel in ∆V consists of N identical subparcels, and each subparcel follows the flow of the ensemble of subparcels (is dragged along or advected) but also undergoes its own diffusion. Heuristically, this means that the interaction in the IPS is due to the deterministic part of the motion where each infinitesimal subparcel follows the mean flow of the ensemble.
The corresponding total Hamiltonian (3.22) for the fluid in M is then obtained by integrating the energies over the infinitesimal subvolumes ∆V α and summing over all indices α. Since this Hamiltonian describes a system consisting of N subparcels, it is a function which is the N-fold product of the phase space of compressible fluid mechanics. (The semidirect product notation is explained in Section 2.C.) The canonical symplectic structure on the phase space then yields a Hamiltonian system of Stratonovich SDEs by adapting the construction of [39] to the infinite dimensional setting. This system of N interacting SDEs will be called the 'HIPS equations of motion'. In Section 3.A it is explained how the Hamiltonian H N is a sum of terms involving: translational kinetic energy of the particle ensemble, equilibrium internal energy associated to the hydrostatic pressure, equilibrium internal energy associated to capillary forces, non-equilibrium internal energy due to expansion/compression along the flow, and stochastic energy associated to molecular bombardment. The HIPS equations of motion are derived in Section 3.B, see (3.25)- (3.27). Under the assumption that the mean field limit exists, as N → ∞, the stochastic mean field equations are obtained in Section 3.C.
The mean field SDEs (3.32)-(3.33) represent the Eulerian description of the motion of a fluid parcel associated to a subdivision ∆V α for very large N. The expected fluid flow is then obtained by averaging over momenta and mass densities of all the smaller fluid parcels. But these are just the mean fields. It therefore remains to calculate the equations of motions for the latter. These equations are deterministic and given in (3.34)- (3.35).
In Section 3.D it is shown that, if the stochastic perturbation corresponds to a Brownian motion for each Fourier mode (i.e., is a cylindrical Wiener process in the space of solenoidal vector fields), then the momentum and mass density mean fields solve the compressible Navier-Stokes system (1.1)-(1.2). This is the content of Theorem 3.4. Moreover, since the HIPS is invariant under the particle relabeling symmetry (invariance with respect to the group of diffeomorphisms) Noether's Theorem implies that a Kelvin Circulation Theorem holds along stochastic Lagrangian paths (Proposition 4.1).
The density dependence of the viscosity, as manifest by the ρ-factor in (1.3), is a consequence of the form of the Hamiltonian (3.22). The compressible Navier-Stokes equations are often expressed with respect to a density-independent viscosity, but it is not clear how to realize this independence in the HIPS framework.
Mean field representations of solutions to the incompressible Navier-Stokes equation have been previously obtained by [13,34] by considering a Weber functional along stochastic Lagrangian paths. In [29,30] the above described Eulerian (H)IPS formulation has been used to also obtain a stochastic mean field representation for solutions to the incompressible Navier-Stokes equation.
To the best of my knowledge, Theorem 3.4 is the first stochastic representation of solutions to the compressible Navier-Stokes equations.
Related approaches to fluid dynamics from the perspective of stochastic variational principles include [54,14,18], which characterize solutions to the incompressible Navier-Stokes equations by variational principles for stochastic Lagrangian paths, and generalizations in [4,11].
1.C. Applications: NatCat models and Solvency II capital requirements. Uses of the Navier-Stokes equation range from semi-conductor engineering to astrophysics and it is not the goal of this section to attempt a review of these topics. Rather, I want to briefly describe an application where the interaction between academia and industry is perhaps not very well established. This concerns models of natural perils (NatCat models) that are used in the insurance industry to calculate risk capital requirements. These risk capital requirements are a determining factor for the solvency of a given company. In the EU the relevant regulatory framework is called Solvency II ( [55]). There exist different NatCat models, which are in general proprietary, for natural disasters such as earthquakes, flooding, tropical cyclones, and extratropical cyclones (European winter storms). Storm models are often based on numerical weather prediction (NWP) systems, and therefore inherit all their advantages and flaws.
However, it is not claimed that the stochastic HIPS formulation of this paper is appropriate to generate stochastic NatCat storm scenarios. Thus the material in this section is only intended as additional background information concerning possible applications of stochastic fluid mechanics, and it is logically independent from the main result, Theorem 3.4.
Numerical weather prediction (NWP) and climate modeling. NWP systems and climate models are important for a number of apparent reasons such as daily weather forecasts or climate change quantification.
Geophysical flows are modeled by the compressible Navier-Stokes equations (without capillary term, i.e. C = 0 in (1.1)). These equations are deterministic. However, any implementation of these equations needs to introduce temporal and spatial discretizations. Physical processes which occur below these chosen grid scales ('subgrid phenomena') cannot be accounted for by any numerical model of the deterministic flow equations. Since the subgrid processes are inherently unknown and uncertain it seems reasonable to model these by a stochastic dynamics point of view. Indeed, such a position has been adopted quite early by Kraichnan [38]. See also [46] for a modern version and fluid equations with stochastic force terms.
Reviews regarding NWP systems and climate models are contained in [8,9,47]. These also address the need for stochastic parameterizations of unknown subgrid processes.
Recent advances in the stochastic modeling of geophysical flows include the 'location uncertainty' approach of Mémin, Resseguier and collaborators [44,49,50,51] as-well as the 'stochastic advection by Lie transport (SALT)' theory of Holm and collaborators [32,15,16,22,2]. The SALT approach is based on the observation that subgrid phenomena represent unknown physical processes and should therefore be derived from a stochastic variational principle. As a consequence, these models preserve circulation along stochastic Lagrangian paths.
NatCat models and solvency capital requirement (SCR). With the implementation of the Solvency II regulatory regime ([55]) per 1. January 2016, applied insurance mathematics has become a surprisingly diverse and multi-disciplinary subject. The relevant tools extend beyond classical actuarial science to, e.g., modeling of local general accounting principles ( [31,26]), no-arbitrage principles and stochastic interest rate models ( [53,24]), and stochastic fluid dynamics. The significance of the last point in this (certainly incomplete) list is briefly explained below.
One of the basic quantitative principles of Solvency II can be summarized as follows: Insurance and reinsurance undertakings are required to quantify all relevant risk factors over a one-year horizon and derive a corresponding loss distribution. Now, the own funds (i.e. excess of assets over liabilities) have to cover the 99.5 percentile of this distribution ('survival of the 200 year event'). This 99.5 percentile corresponds to the so-called solvency capital requirement (SCR). The SCR can be calculated from a prescribed standard formula or a company specific internal model. Medium and large sized companies generally use internal models. If a company has chosen the internal model approach and covers windstorm risks, then the corresponding loss distribution for a one-year period has to be derived. To do so the following three-step procedure, or a variant of it, is used ( [27]): (1) Hazard module: physical model. A large set of windstorm scenarios is generated. These are the so-called stochastic scenarios of the NatCat model. (2) Vulnerability module. This step quantifies the vulnerability (i.e., damage done) of the insured structures for each of the stochastic scenarios.
(3) Financial loss module. Finally, structural damage is converted to loss by taking into account contract specifics (e.g., sum insured) for each damaged structure and, possibly, reinsurance.
Since the relevant time period for NatCat models is one year they operate on a new temporal scale when compared to short term weather forecast models and medium or long term climate models. The most advanced NatCat models are based on a coupling of a global circulation model (GCM) and a NWP system. Since these models are proprietary it is not possible to cite a suitable model documentation. However, [12] contains a review of the basic method, which is still quite up to date. In particular, it is explained how (deterministic) NWP systems are used to generate a set of scenario events by statistically sampling the initial conditions. "Using NWP technology, a large set of potential future storms is generated by taking data sets comprising the initial pressure fields of historical storms, perturbing them both temporally and spatially, and moving them forward in time through the application of a set of partial differential equations governing fluid flow. The resulting event set is rigorously tested to ensure that it provides an appropriate representation of the entire spectrum of potential storm experience -not just events of average probability, but also the extreme events that make up the tail of the loss distribution." ( [12]) This point will be taken up again in Section 5. A recent reanalysis of European winterstorm events, which also highlights the potential financial risks for the insurance sector, has been carried out by [33]. The L 2 scalar product ., . on g is defined by for ξ, η ∈ g, where dx is the standard volume element in M, and ., . is the Euclidean inner product. Via R g this can be extended to a right invariant Riemannian metric on Diff(M). See [7,25,40,45] for further background.

2.B.
Derivatives. The adjoint with respect to ., . to the Lie derivative L, given by Y j e i with respect to the standard basis e i , i = 1, . . . , n. The notation ad(X)Y = [X, Y ] = −L X Y and ad(X) ⊥ = −L ⊤ X will be used.
The variational derivative of a functional F : g → R will be denoted by δF/δX, that is for X, Y ∈ g. The semi-direct product structure is defined by the right action The corresponding phase space is trivialized with respect to this right multiplication as The variables µ ∈ g * and ρ ∈ F (M) * represent momentum and mass density, respectively. Making use of the Euclidean volume form dx, the duals can be identified as g * = Ω 1 (M) and F (M) * = F (M). Details on Hamiltonian mechanics on semi-direct products are given in [41], where also the case of compressible ideal fluids is treated.
2.D. Stochastic dynamics. Let (Ω, F , (F t ) t∈[0,T ] , P ) be a filtered probability space satisfying the usual assumptions as specified in [48]. In the following, all stochastic processes shall be understood to be adapted to this filtration. The symbol (2.9) δ t will be used to denote the Stratonovich differential to distinguish it from the variational derivative δ. The Ito differential does not appear in this paper. The exterior differential is d.
2.E. Brownian motion in g 0 . Let . . , k ⊥ n−1 denote a choice of pairwise orthogonal vectors in R n such that |k ⊥ i | = |k| and k ⊥ i , k = 0 for all i = 1, . . . , n − 1. Consider the following system of vectors in g 0 ( [17,19]): where e j ∈ R n is the standard basis and s is the Sobolev index from Section 2.A. By slight abuse of notation we identify these vectors with their corresponding right invariant vector fields on Diff(M) 0 .
Further, in the context of the ζ r vectors we shall make use of the multi-index notation r = (k, i, a) where k ∈ Z + n and a = 0, 1, 2 such that Thus by a sum over ζ r we shall mean a sum over these multi-indices, and this notation will be used throughout the rest of the paper. It can be shown (see [19,Appendix] for details) that the ζ r form an orthogonal system of basis vectors in g 0 , such that where W r t are independent copies of Brownian motion in R. Then W defines (a version of ) Brownian motion (i.e., cylindrical Wiener process) in g 0 .

HIPS for compressible flow
This section is concerned with the HIPS equations of motion and the mean field limit. Background on mean field theory can be found in [52,1,20,42,35].

3.
A. The Hamiltonian. Consider a barotropic fluid in M = R n /Z n . At a (macroscopic) position x ∈ M consider a tiny volume element ∆V x . Suppose ∆V x is further divided into N infinitesimal volume elements, ∆V α x , labeled by α = 1, . . . , N, which are all assumed to have identical dimensions. In each volume element there is a mass density ρ α = ρ α (x). Thus it is assumed that the continuum description of the fluid also holds at the level of the subdivsion. Let the fluid element in volume ∆V α x have velocity v α (x). Thus the energy of the particle ensemble in the total infinitesimal volume is determined by the velocities, v α , and densities, ρ α , in the subdivisions. To arrive at the total energy (3.22) we consider five contributions: (1) Translational kinetic energy; (2) Equilibrium internal energy U which gives rise to the hydrostatic pressure p; (3) Equilibrium capillary energy; (4) Non-equilibrium expansion/compression energy; (5) Stochastic energy due to molecular bombardment; Translational kinetic energy. The total velocity at x, v(x), is the weighted average .
The total momentum, µ(x), is therefore and the translational kinetic energy K N (x) of the particle system is Equilibrium internal energy. Notice that the overall mass density in ∆V x is expressed, in this subdivision picture, as ρ(x) = ρ α (x)/N. Consequently, the mass contained in ∆V x is The barotropicity assumption implies that the specific equilibrium internal energy in ∆V depends on the overall mass density ρ α /N. Then the equilibrium internal energy in ∆V is Capillary energy. Let κ ≥ 0 be a constant. The capillary energy is defined as which depends, again, on the overall mass density. This form coincides with the capillary contribution to the Helmholtz free energy [3, Equ. (5)].
Non-equilibrium expansion energy. Let ν ≥ 0 be a constant. Assume temporarily that ∆V x is a box which is aligned along Cartesian coordinates e 1 , e 2 , e 3 and that the flow has only velocity components pointing in the direction of the e 1 -axis. Let L denote the set of labels α such that the corresponding subvolumes ∆V α constitute the left wall of ∆V x while R denotes those which correspond to the right wall (viewed along e 1 ) of the volume. Now, if (3.17) α∈L ρ α v α α∈L ρ α − α∈R ρ α v α α∈R ρ α is greater than 0, then particles moving into ∆V are faster than those moving out, and the corresponding energy difference should contribute to the (non-equilibrium) internal energy of the system. Since this expression is proportional to minus the divergence of the barycentric velocity, we propose a non-equilibrium internal energy expression Stochastic energy. The energy due to molecular bombardment along a vector field ξ r in ∆V α x ≈ dx is See [39,28,32,29,30]. This stochastic perturbation corresponds to individual molecules imparting their velocities, namely ξ r , on the macroscopic fluid element in ∆V α . Since individual molecules are incompressible the ξ r are assumed to be divergence free.
Total energy. The above subdivision formulation implies that the total configuration space is of the form Π N α=1 (Diff(∆V α ) F (∆V α ))). But now the ∆V α , which are all identical by assumption, are identified with the infinitesimal element dx in M. Thus each index α corresponds to a copy of M, and this can be done since the for the state variables, velocity and density, it does not make a difference whether these are regarded at the ∆V or at the ∆V α level. Therefore, letting the position x ∈ M range over the full domain, the total configuration space is Π N α=1 (Diff(M) F (M)) = (Diff(M) F (M)) N . Let us switch from velocity and density to momentum (density), µ α = ρ α v α , and density, ρ α , as state variables. The phase space is thus T * (Diff(M) F (M)) N = (T * (Diff(M) F (M))) N . We use the Euclidean metric to identify each copy in the phase space, which is the regular dual, as where the last identification follows from right-multiplication in the semi-direct product group, see Section 2.C. The resulting Hamiltonian of the IPS will therefore be a function (3.20) H N : the total Hamiltonian is the sum of (3.14), (3.15), (3.16), (3.18) and (3.19), and given in semi-martingale notation as where W r,α are pairwise independent Brownian motions such that [W r , W s ] t = δ r,s t, the solenoidal vector fields ξ r ∈ X 0 (M) are fixed and ε ≥ 0 is a constant. The Hamiltonian is right-invariant by construction as it does not depend on (Φ α , f α ) ∈ Diff(M) F (M).
3.B. HIPS equations of motion. The phase space (3.20) is an N-fold direct product of a tangent bundle identified with its dual. It carries therefore the corresponding direct product canonical symplectic form. Since the Hamiltonian H N does not depend on (Φ α , f α ) N α=1 ∈ (Diff(M) F (M)) N we can pass via Lie-Poisson reduction to the phase space (X (M) F (M)) N .
The Hamiltonian IPS equations of motion follow therefore from the variational derivatives (again, in the semi-martingale notation) by using the N-fold product of the semi-direct product structure ( [41]). Here the abbreviations µ N := β µ β N , ρ N := β ρ β N and u N := µ N ρ N are used. Remark 3.3. The quantity u N is the specific momentum (viewed as a vector field) of the ensemble average. But it is (for non-constant density) not equal to the empirical average of velocities, that is The stochastic Hamilton equation associated to (3.22) for (Φ α , µ α , ρ α ) are: Here ρ α is viewed as a density whence L is the Lie derivative of a density, not of a function. The momentum variable is identified, via the Euclidean metric, as an element whence the transpose Lie derivative L ⊤ is used instead of L * . The diamond notation in (3.26) is defined by f ⋄ ρ = ρ∇f and this term arises because of the semi-direct product structure.

Stochastic Kelvin Circulation Theorem
In [21] it is shown that stochastic Euler-Poincaré fluid equations are characterized by preserving circulation along Lagrangian paths. Since (3.32)-(3.33) are obtained as a mean field limit of a Hamiltonian IPS, and can be viewed of as mean field generalization of the stochastic fluid system in [21], there should be a Kelvin Circulation Theorem: Proposition 4.1. Let C be a smooth closed loop which is transported by the Lagrangian flow Φ t , defined through the mean field limit (3.29) and characterized by Let µ t and ρ t be solutions of (3.32) and (3.33). Then where ♭ is the Euclidean isomorphism to one-forms (since µ is treated as a vector field).

5.
A. HIPS approach to the compressible Navier-Stokes equation. The mean field system (3.32)-(3.33) is derived from the interacting particle point of view under the basic assumption that the equations of motion follow from stochastic Hamiltonian mechanics. Therefore, circulation is preserved along stochastic Lagrangian paths. If the perturbation fields ξ r run over the orthogonal system ζ r , defined in Section 2.E, such that the stochastic perturbation is given by a cylindrical Wiener process, then the mean fields E[µ] and E[ρ] solve the compressible Navier-Stokes equations. (Theorem 3.4 and Proposition 4.1.) While the HIPS formulation relies on a system of N interacting SDEs, the mean field equations (3.32)-(3.33) is a single SDE system for momentum and mass density. In contrast to statistical mechanics, this mean field formulation is obtained without any closure assumptions. (However, in this paper the existence of the mean field limit is not proved but assumed.) In the mean field limit (3.29), the HIPS evolution equation (3.25) becomes which is of similar form as the LA SALT advection field [22,2], where u L t is a stochastic velocity field. However, there are a few crucial differences: while (5.50) is the starting point for LA SALT theory, the HIPS formulation is based on the ensemble Hamiltonian (3.22) and the Lie transport along (5.49) in the mean field equations (3.32)-(3.33) is a consequence of the Hamiltonian structure of the IPS (and the ensuing passage to the mean field limit). Moreover, unless the density is constant, it is not clear how to identify the drift in (5.49) with the expectation of a velocity. Thus, both, the starting points and the advection fields are different. However, the perturbation fields ξ r can be interpreted in the same manner. 5.B. NatCat modeling of windstorm events. NatCat models used to calculate the solvency capital requirement (SCR as defined in [55]) for storm risks rely on NWP systems. These NWP systems are deterministic and, to arrive at a set of 'stochastic' NatCat scenarios, the initial conditions are statistically sampled. As discussed in the Introduction, all such numerical schemes suffer from subgrid phenomena, and for geophysical flow models a well-established means for treating these deficiencies is by stochastic fluid mechanics ( [2,9,32,44,49]). Since SCR calculation is concerned with predicting extreme events in the 99.5 percentile, and not only average storm patters, it seems reasonable to expect that also NatCat models would benefit from a stochastic dynamics approach.
However, it is not claimed that the stochastic HIPS formulation of this paper is appropriate to generate stochastic NatCat storm scenarios.