Macole Sabat

Aymeric Vié
email: aymeric.vie@centralesupelec.fr

Adam Larat

Marc Massot

Statistical description of turbulent particle-laden flows in the very dilute regime using the Anisotropic Gaussian Moment Method

Keywords: Disperse phase, Eulerian models, Anisotropic Gaussian, Particle Trajectory Crossing, Turbulence

. In order to validate this model further, we investigate here 3D turbulent flows that are more representative of the mixing processes, which occurs in realistic applications. The chosen configuration is a 3D statistically-stationary Homogeneous Isotropic Turbulence (HIT) loaded with particles in a very dilute regime. The analysis focuses on the description of the first three lowest order moments of the particulate flow: the number density, the Eulerian velocity and the internal energy. A thorough numerical study on a large range of particle inertia allows us to show that the AG closure extends the ability of the Eulerian models to correctly reproduce the particle dynamics up to a Stokes number based on the Eulerian turbulence macro-scale equal to one, but also highlights the necessity of high-order numerical schemes to reach mesh convergence, especially for the number density field.

Introduction

In several industrial systems such as aeronautical combustors or fluidized beds, turbulent particle-or droplet-laden flows play a central role. Simulating such flows is still challenging because of the multi-scale character of the interactions between the gaseous carrier phase and the disperse phase. Several methods are available in the literature, ranging from the simulation at the particle scale [START_REF] Subramaniam | Lagrangian-Eulerian methods for multiphase flows[END_REF] to the ensemble-averaged simulation of a full system [START_REF] Capecelatro | Strongly coupled fluid-particle flows in vertical channels[END_REF].

The present work focuses on the prediction of the dynamics of particles subject to drag force with a turbulent carrier phase in the very dilute regime. In this context, the particles do not affect the carrier phase and the inter-particle collisions are neglected. A point-particle approximation is considered, for which particle-scale flow fluctuations are dissipated fast enough to be neglected. Under these assumptions, the statistics of the disperse phase, i.e. its Number Density Function (NDF), can be determined by a Generalized Population Balance Equation (GPBE). As stated in [START_REF] Février | Partitioning of particle velocities in gas-solid turbulent flow into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF], a key aspect of this statistical description is the particle velocity dispersion, called Random Uncorrelated Motion (RUM) due to the Particles Trajectory Crossing (PTC). This velocity dispersion is typical for inertial particles that are weakly correlated to the local fluid velocity field.

Solving the PBE directly is not feasible in realistic configurations due to the high dimensionality of the phase space. To circumvent this issue, moment methods can be used. This efficient strategy significantly reduces the computational cost by transporting moments of the NDF that only lie in the physical space.

However, for a given moment set to be solved, unsolved higher order moments are required in order to close the moment equations. The closure must properly describe the velocity dispersion in order to reproduce the physics of the flow.

In the literature, several closures can be found in the particulate flow community, but also in the rarefied gases community. They can be separated into two families depending of the type of approximations: the Algebraic-Closure-Based Moment Methods (ACBMM, [START_REF] Kaufmann | Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence[END_REF][START_REF] Masi | Algebraic-Closure-Based Moment Method for unsteady Eulerian simulations of non-isothermal particle-laden turbulent flows at moderate Stokes numbers in dilute regime[END_REF][START_REF] Masi | Development of an algebraic-closure-based moment method for unsteady Eulerian simulations of particle-laden turbulent flows in very dilute regime[END_REF][START_REF] Le Lostec | Numerical description of dilute particle-laden flows by a quadrature-based moment method, Proceedings of the Summer Program[END_REF][START_REF] Sakiz | Numerical experiments and modelling of nonequilibrium effects in dilute granular flows[END_REF]), which make assumptions on the structure and evolution of the moments, and the Kinetic-Based Moment Methods (KBMM, [START_REF] Groth | Towards physically-realizable and hyperbolic moment closures for kinetic theory[END_REF][START_REF] Lam | Numerical prediction of three-dimensional non-equilibrium gaseous flows using the gaussian moment closure[END_REF][START_REF] Laurent | A hierarchy of Eulerian models for trajectory crossing in particle-laden turbulent flows over a wide range of stokes numbers[END_REF][START_REF] Mcdonald | Towards realizable hyperbolic moment closures for viscous heat-conducting gas flows based on a maximum-entropy distribution[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF][START_REF] Chalons | Multivariate gaussian extended quadrature method of moments for turbulent disperse multiphase flow[END_REF][START_REF] Kong | Eulereuler anisotropic gaussian mesoscale simulation of homogeneous clusterinduced gas-particle turbulence[END_REF] ), which make assumptions directly on the shape of the NDF itself.

By essence, the two families share the objective of controlling moments up to a given order. However, ACBMM always require fewer moments to be solved because, when KBMM add an equation for any additional moments to be controlled, ACBMM can rely on an algebraic relationship. In [START_REF] Vié | On the direct numerical simulation of moderate-stokes-number turbulent particulate flows using Algebraic-Closure-Based and Kinetic-Based Moment Methods[END_REF], a first comparison between the 2 -ΦEASM model [START_REF] Masi | Development of an algebraic-closure-based moment method for unsteady Eulerian simulations of particle-laden turbulent flows in very dilute regime[END_REF] and the Anisotropic Gaussian (AG) model [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] has been performed, showing a comparable accuracy in a 2D turbulent configuration. Both methods were targeting to control the second order moments, but where AG has to solve the full set of second order moments (6 in 3D), 2 -ΦEASM model uses only one moment (the trace of the pressure matrix), leading to a reduced computational effort. However, the design of numerical methods is more difficult for ACBMM because of the presence of extra nonlinear terms preventing from a complete characterization of the mathematical structure. On the other hand, KBMM lead to hyperbolic or weakly-hyperbolic systems of equations [START_REF] Chalons | Beyond pressureless gas dynamics: quadrature-based velocity moment models[END_REF], for which well-designed solvers exist in the literature (see [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical introduction[END_REF][START_REF] Larat | A stable, robust and high order accurate numerical method for Eulerian simulation of spray and particle transport on unstructured meshes[END_REF] and references therein for example).

Even if the AG model has already been validated in 2D and compared with a member of the ACBMM, it still has to be validated in 3D turbulent configurations in order to evaluate its ability to tackle realistic industrial configurations. This is the aim of the present paper.

In the following, the AG model is first presented along with the MonoKinetic model (MK): a model only accurate for small inertia particles [START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF]. Numerical methods are also briefly mentioned, while a more detailed description is provided in Appendix B. Next, we aim at showing that the second-order AG model can be viewed as an extension of the first-order MK model. To this purpose, a 3D turbulent configuration is presented and the results obtained with the AG and MK models are compared to a Lagrangian statistically-converged reference. The numerical study is performed in an as extensive manner as possible, by sweeping a large range of the disperse regimes and looking at numerous representative quantities to assess the strengths and weaknesses of the AG Eulerian model. By this thorough numerical study, we assess the ability of the AG model to recover the right statistics for a specified range of particle inertia.

Statistical modelling of particle-laden flows

The present work relies on the following assumptions:

• Point-particle approximation [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF]: particles are considered to be much smaller than any of the carrier phase scales.

• Negligible mass loading and volume fraction: the particles load does not affect the carrier phase. The two phases are coupled only one-way: the carrier phase feeds the disperse phase with energy, but the particles retroaction is neglected.

• Very dilute regime: the particles mean free path is much greater than any dimension of the problem. The Knudsen number is considered to be infinite and thus, inter-particles collisions are neglected.

• Large density ratio between particle and gas phases: the drag force is the only force acting on the particles.

• Fixed-size monodisperse spherical particles: there is no need to track the distribution in size.

• Stokes drag: particules Reynolds number does not affect their relaxation time toward the underlying carrier velocity field.

In this framework, the NDF, f (t, x, c), of the disperse phase is described by the following GPBE:

∂ t f + ∂ x • ( cf ) + ∂ c • u g -c τ d f = 0, ( 1
)
where c is the velocity of the particles in the phase space, u g is the gas velocity and τ d is the characteristic relaxation time of the particles.

Eulerian method of moments

After multiplying equation ( 1) by c i 1 c j 2 c k 3 and integrating over the phase space, one gets a system of moment equations. In a three-dimensional physical space, the general equation on the velocity moments is:

∂ t M i,j,k + ∂ x •      M i+1,j,k M i,j+1,k M i,j,k+1      = 1 τ d      u g •      iM i-1,j,k jM i,j-1,k kM i,j,k-1      -(i + j + k)M i,j,k      , (2)
where the general (i + j + k) th order moment in velocity is:

M i,j,k = c i 1 c j 2 c k 3 f dc 1 dc 2 dc 3 . (3)
This system is not closed: actually for every set of moments of order N which contains the moments of order (i + j + k) ≤ N , moments of order N + 1 are needed to describe the higher order fluxes in physical space: M i,j,k where

i + j + k = N + 1 .
In order to close the moment system (2), one needs to provide:

• the moment set of order N :

{M i,j,k , i + j + k ≤ N } ,
• a closure relationship which allows to model the unknown fluxes. In the Kinetic Based Moment Methods (KBMM), this is done by keeping the link between the kinetic level and the macroscopic level, at the price of the reconstruction of the NDF, which is generally supposed to have a given shape.

The choice of the presumed shape of the NDF is of paramount importance for a KBMM. It can be driven by the following guidelines:

• a one-to-one mapping between the moments and the parameters of the NDF is essential to ensure the well-posedness of the system of equations;

• an appropriate shape can also guarantee that the final moment system is hyperbolic and has an entropic structure, usually relying on hydrodynamic equilibrium assumption as in the Levermore hierarchy [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF];

• the shape can also be chosen based on the physics of the problem in sight, which can for instance suggest an equilibrium distribution (such as a Maxwellian distribution in the case of collision-driven particulate flows);

• peculiar attention must be given to the border of the moment space (zero pressure, zero density), to avoid numerical issues when dealing with transition between regions at the border and regions inside the moment space.

Monokinetic model

When no PTC occur, the disperse phase has a unique deterministic velocity at each time-space point (t, x). Consequently, the velocity distribution is monokinetic:

f (t, x, c) = n(t, x)δ( c -u(t, x)), (4)
where n(t, x) is the total number density of the disperse phase, and u(t, x)

the Eulerian velocity. In this specific case, the system of moments closes at first order, since the pressure tensor, its second order centered-moment, is null:

P = ( c -u(t, x)) ⊗ ( c -u(t, x))f (t, x, c) d c = 0.
Therefore, we get the following weakly-hyperbolic system, identical to the Pressureless Gas Dynamics (PGD) system [START_REF] Bouchut | On zero pressure gas dynamics[END_REF][START_REF] Marble | Dynamics of dusty gases[END_REF]:

         ∂ t n + ∂ x • (n u) = 0, ∂ t (n u) + ∂ x • (n u ⊗ u) = n( u g -u) τ d . (5)
The MK model correctly reproduces the formation of depletion and accumulation regions obtained for small Stokes numbers and dedicated numerical schemes can be designed to properly account for such hypercompressible dynamics [START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF]. However, when PTC occur for moderate to high inertial particles, the velocity dispersion cannot be captured by the single Dirac δ-function used for the velocity distribution. Instead, at each PTC location, unphysical δ-shocks are generated, and these are highly constraining from the numerical point of view. Another constraint of this model is that it does not conserve the total kinetic energy of the disperse phase [START_REF] Bouchut | On zero pressure gas dynamics[END_REF][START_REF] Brenier | Sticky particles and scalar conservation laws[END_REF].

Anisotropic Gaussian model

To reproduce the velocity dispersion induced by the PTC, we consider the AG closure, which is a second order model of the KBMM family:

f (t, x, c) = n(t, x)N ( c -u(t, x), Σ(t, x)), (6)
where N is a joint Gaussian distribution centered in u and of covariance matrix

Σ = (σ ij ) in the space of dimension N d : N ( c -u, Σ) = |Σ| -1/2 (2π) N d /2 exp - 1 2 ( c -u) T Σ -1 ( c -u) . (7)
This closure was inspired by the work of [START_REF] Levermore | The Gaussian moment closure for gas dynamics[END_REF] in the context of out-of-equilibrium rarefied gas dynamics and was called the 10-moment closure, and is the first member of the Levermore hierarchy [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF]. In the rarefied gases context, the good mathematical properties of this model at the kinetic and at the moment levels was highlighted in [START_REF] Holway | New statistical models for kinetic theory: methods of construction[END_REF][START_REF] Andries | The Gaussian-BGK model of Boltzmann equation with small Prandtl number[END_REF] and a numerical approximation was proposed by [START_REF] Brown | Numerical solution of a 10-moment model for non-equilibrium gas dynamics[END_REF] and [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF].

The resulting system of moments reads:

                       ∂ t n + ∂ x • (n u) = 0, ∂ t (n u) + ∂ x • (n u ⊗ u + P) = n( u g -u) τ d , ∂ t (nE) + ∂ x • ((nE + P) ∨ u) = n( u g ∨ u -2E) τ d , (8)
where ∨ denotes the symmetric tensor outer product 1 , the total energy tensor E is given by E = 1 2 u ⊗ u + P 2n and the anisotropic pressure tensor is given by the equation of state P = nΣ. The resulting system is hyperbolic and admits entropies [START_REF] Berthon | Numerical approximations of the 10-moment Gaussian closure[END_REF]. This model was introduced for the first time for the modeling of the disperse phase in the work of [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF].

Realizable numerical schemes

An important requirement for the moment methods is to preserve the realizability of the moment set, i.e. the fact that the moment set is the set of moments of a positive NDF. Because the closure of the moment set is obtained through assumptions on the NDF itself, the derived systems of moments preserve the realizability of the moment set at the continuous level. When it comes to solve the moment system using a given time and space discretization, introducing numerical errors, it is not straightforward to maintain the realizability of the moments at the discrete level, especially since we want to cover the whole moment space, including its frontier (zero density, zero pressure, ...).

In the present work, the two systems of interest are solved using numerical methods expressly developed to ensure realizability and stability, as well as high accuracy: a dedicated second order MUSCL-HLL scheme.

For the AG model, the MUSCL-HLL scheme developed in [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] has been extended to the 3D case. It is of second-order in time and space, also preserves realizability of the moment set and behaves well at moment space boundaries.

More details are given in Appendix B.

It is important to note that for the Monokinetic approach, a second-order Kinetic Scheme [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] can also be used. The realizability properties of this latter are inherited from a linear and conservative reconstruction of the central moments n and u, and an exact evolution in time of the reconstructed NDF at the kinetic level.

Numerical setup: 3D Forced Homogeneous Isotropic Turbulence

The configuration of interest is based on a 3D Forced Homogeneous Isotropic Turbulence (FHIT) for the carrier phase, solved in a cubic domain, with periodic boundary conditions. This classical academic test case, auto-similar in terms of Reynolds number, is essential for the validation of the disperse phase models, since it represents one of the canonical samples of the physics faced in industrial applications. These simulations are carried out at a low Reynolds number.

However, let us emphasize here that the evaluation and the application of the Anisotropic Gaussian model is not limited to low turbulence flows.

All the values of the variables presented in this work are dimensionless:

u * = u u ref , t * = t t ref , L * = L L ref , (9)
and the normalizing Reynolds number is Re

0 = u ref L ref /ν = 10 3 .

Velocity field and general properties

The FHIT is solved using the ASPHODELE code [START_REF] Reveillon | Effects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF], which solves the Navier-Stokes equations in the low-Mach number limit.

The initial turbulence field is generated using an isotropic random mode generation and inverse Fourier transform in a cubic domain of size 5 3 using a 128 3 mesh. The energy spectrum function introduced in [START_REF] Pope | Turbulent flows[END_REF] is used as the model spectrum:

E m (κ) = C ε 2/3 κ -5/3 f L (κL)f η (κη), (10)
with:

f L (κL) = κL [(κL) 2 + c L ] 1/2 5/3+p0 , (11)
f η (κη) = exp -β[((κη) 4 + c 4 η ) 1/4 -c η ] , (12)
where

η = ν 3 ε 1/4
is the Kolmogorov length scale given as a function of the carrier phase viscosity ν and the mean dissipation rate of turbulent kinetic energy ε . In addition, the length scale L = TKE The turbulent characteristics of the statistically-stationary HIT, including the mean turbulent kinetic energy ( TKE ), the dissipation ε , the Kolmogorov time and length scales (time τ K and length η), the Lagrangian integral time scale τ L , the Eulerian eddy turn-over time τ eddy , the longitudinal integral length scale L 11 , and the turbulent Reynolds number Re = u g,rms L 11 /ν, where u g,rms is the gas velocity fluctuations, are given in Table 1. been reported in the literature, see for instance [START_REF] Zaichik | Statistical models for predicting pair dispersion and particle clustering in isotropic turbulence and their applications[END_REF], and is of great interest for modelling purpose, as it could suggest that even low-Stokes number flows require the modelling of the second order moments.

Eulerian and Lagrangian simulations for the particles

Here we give details about the simulations that have been performed. For the initialisation, 100M lagrangian particles are randomly loaded in the domain following a homogeneous distribution. Eulerian simulations are also initially

Results and discussions

The results are presented in a sequential manner, comparing the AG and MK model to a Lagrangian reference. Eulerian quantities for the Lagrangian data are obtained using a box average based on the computational mesh. We first investigate the ability to accurately reproduce the number density field. Then, results on the Eulerian velocity are presented. Finally, the ability to predict the energy partition in the particles dynamics is evaluated.

Segregation

In order to get a quantitative estimate of the level of accumulation and of depletion in the number density field, we investigate the segregation, which represents the spatial correlation of the number density field at a given length scale [36]. This quantity is of paramount importance in the domain of combustion since it strongly affects the evolution of the mean vapour mixture fraction [START_REF] Reveillon | Effects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF].

The segregation G ∆ pp is given by the following equation:

G ∆ pp = n 2 ∆ n ∆ 2 , ( 13
)
where • is the averaging operator over the whole domain, and n ∆ is the number density field projected onto a mesh of grid size ∆, using a box filter. In the following, ∆ is chosen to be equal to the largest mesh spacing used for the Eulerian simulation, i.e. ∆ = L 64 . The time evolution of the segregation is plotted in Fig. 3 for three Stokes numbers and several mesh refinements.

For St K = 1, both Eulerian models underestimate the segregation compared to the Lagrangian reference. The MK model gives higher segregation than the AG model, and is therefore closer to the Lagrangian segregation. When refining the mesh, both Eulerian models yield higher segregation and seem to "converge" to the Lagrangian reference. However for the 512 3 mesh, the segregation of MK slightly overestimates the Lagrangian one for t < 20, showing that this model is not converging to the Lagrangian result. This can be explained by the fact that the occurence of local strain rate higher than the average value used to compute τ k can lead to PTC that cannot be handled correctly by the MK model.

In the case of moderately inertial particles St K = 3, the segregation is overestimated by the MK model and diverges from the Lagrangian segregation curve, when refining the mesh. On the other hand, the segregation level for the AG model approaches to the Lagrangian one, when refining the mesh.

For inertial particles at St K = 8, the segregation is still largely overestimated by the MK model, whereas the AG model is capturing the overall time evolution, and is getting closer to the Lagrangian reference as we refine the mesh. However, it should be noted that for long time t = 80, the most refined solution of the AG model gives a slightly higher segregation than the Lagrangian reference.

From both high Stokes numbers, we clearly see the two main sources of "diffusion" of the Eulerian solutions, the numerical error and the model:

• the numerical error is present in both MK and AG simulations, and exhibits a first-order convergence in terms of segregation because the numerical schemes are necessarily switching to first order in accumulation zones and close to vacuum. The high sensitivity of Eulerian simulations to numerical error has already been evidenced in the literature [START_REF] Riber | Modeling turbulent two-phase flows using large-eddy simulation[END_REF][START_REF] Kaufmann | Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence[END_REF][START_REF] Dombard | Direct Numerical Simulation of non-isothermal dilute sprays using the Mesoscopic Eulerian Formalism[END_REF][START_REF] Sierra | Modeling the dispersion and evaporation of sprays in aeronautical combustion chambers[END_REF];

• for a given mesh, comparing AG and MK shows the "diffusive" effect of solving for an additional internal energy: the AG model always giving a lower segregation than the MK model. Arrows indicate increasing mesh sizes for Eulerian simulations.

Distribution functions

To further investigate the number density field, we look at its distribution in the physical domain for each Stokes number. At the critical Stokes number, St K = 1, both Eulerian results approach the number density PDF of the Lagrangian result (Fig. 4) with mesh refinement. The number density PDF of the MK-512 3 result is however closer to the Lagrangian PDF than the most refined AG one. This kind of differences has also been evidenced on the literature, see [START_REF] Kaufmann | Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence[END_REF]. Nonetheless, the MK result might diverge from the Lagrangian one with even more refined meshes, since the PTC might occurs because of the smallest time scales of the turbulence.

For a moderately inertial disperse phase (St K = 3), see Fig. 5, the number density PDF of the AG result approaches the Lagrangian one when refining the mesh, except for the underestimation of the vacuum zone. The MK number density distribution is not converging to the Lagrangian one when refining the mesh and gives a higher number of mesh cells where the number density is zero or very small due to the singularities created by this model at each PTC. At an even higher Stokes number (St K = 8), the distribution of the number density is plotted in Fig. 6, and the AG model is shown to properly capture the number density distribution compared to the Lagrangian reference.

Eulerian velocity

The second moment of interest is the Eulerian velocity u. Here, we directly investigate the distribution function of its component in the x-direction, see Fig. 7. For St K = 1, we see again that the MK and AG models are close to the Lagrangian reference, as only a negligible amount of PTC occur and the main part of the particle dynamics is contained into the Eulerian velocity. For St K = 3, the AG result matches the Lagrangian one except near zero velocity.

For the MK results the velocity in the x-direction is less concentrated near the zero velocity and more spread toward the maximum and minimum values (Fig. 7). The velocity PDF of the MK result is not affected significantly by mesh refinement. For the AG result, the distribution level of the AG velocity approaches the Lagrangian one with mesh refinement, this can be clearly seen near zero-velocity. In fact, in this case the PTC now affect the motion of the disperse phase because St K > 1. Finally, for St K = 8, the AG model is again close to the Lagrangian reference, and the velocity distribution is not impacted by mesh refinement.

Total and internal energies

As a final validation, we investigate the mean total and internal energies:

MIE = n tr (Σ) 2 n , ( 14
) MTE = n i u 2 i + tr (Σ) 2 n . ( 15
)
The study of the evolution of the mean internal energy (MIE) is essential, since it is through the velocity variance or the internal agitation energy that the AG model is capable of reproducing statistically the PTC. Indeed, at each crossing an amount of the kinetic energy n 2 i u 2 i is transferred into internal energy n 2 tr (Σ). Then, the MIE is an evaluation of the occurrence of PTC. This mechanism of transfer of energy between the Eulerian velocity and the MIE has been clearly depicted in [START_REF] Février | Partitioning of particle velocities in gas-solid turbulent flow into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF], when looking at the equation of the internal energy:

∂ t (n tr (Σ)) + ∂ x • (n tr (Σ) u) = - 2ntr (Σ) τ d -nΣ : ∂ x u. ( 16
)
where : is the contraction product of two tensors. In this equation the last term corresponds to the transfer of energy from the Eulerian velocity to the internal energy. Therefore, at crossing events, energy is transferred from the eulerian velocity to the internal energy.

Next, since the mean total energy (MTE) is the total amount of energy available in the whole domain, the knowledge of both MIE and MTE allows to evaluate the transfer of mesoscopic momentum between the kinetic (correlated motion) and the internal (uncorrelated motion) energies. In theory, this energy transfer can be quantified, see [START_REF] Février | Partitioning of particle velocities in gas-solid turbulent flow into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF].

In Fig. 8, the statistics are plotted for St K = 1. In Fig. 8(b) and 8(c), the MTE is underestimated by the MK model, because of the destruction of the kinetic energy during PTC. On the contrary, the AG model performs better and better with mesh refinement. Its MIE is higher compared to the Lagrangian reference but converges toward the right level when the numerical dissipation decreases, thanks to the mesh refinement. In the case of moderately inertial particle St K = 3 (Fig. 9), the departure of the MK from the Lagrangian result is obvious. For the AG model, refinement improves the quality of the results for the MTE, but MIE seems to converge to a higher value than the Lagrangian reference.

For an even more inertial disperse phase (St K = 8, Fig. 10), the AG model has got the right long-time behaviour, but the MTE is now underestimated.

PDF of the trace of the pressure tensor

To accurately evaluate the degree of PTC captured by the AG model, the distribution functions of the trace of the pressure tensor for the Lagrangian and AG simulations on different meshes are plotted in Fig. 11. For St K = 1, the distribution functions of the AG results converge to the Lagrangian one when refining the mesh. Most of the particle eulerian field has zero or very small pressure, since, as explained previously, the probability of having local Stokes numbers higher than the critical one is non-negligible.

For moderately inertial particles, St K = 3, the form of the PDF of the trace of the pressure tensor of the AG result matches the Lagrangian one at high pressures. However, the number of low-pressure particles is underestimated by the AG results, especially on the coarsest mesh. When refining the mesh, the difference between the AG and Lagrangian distribution functions of the internal energy decreases. Actually, the internal energy includes two contributions:

• one is purely related to the model and is the reason why AG is able to capture PTC;

• an other one is due to numerical error; because the numerical scheme is "mixing" states of different velocities and energies, it necessarily leads to the creation of "numerical" PTC, which contribute to the overall internal energy.

This explains the convergence of the PDF of the internal energy using the AG model when refining the mesh, since by doing so, the second contribution is decreased and the physical behaviour becomes predominant; whereas for the coarsest mesh, the dominant effect is the numerical diffusion.

For the highest Stokes number, St K = 8, the PDF of the trace of the pressure tensor for the AG result is still close to the Lagrangian one. for which these PTC will contribute to an additional subgrid scale pressure [START_REF] Chalons | On the Eulerian Large Eddy Simulation of disperse phase flows: an asymptotic preserving scheme for small stokes number flows[END_REF].

density field using MK is close to the Lagrangian one, even if some small details are slightly different. The result on the 512 3 grid captures the right structures, the correct vacuum zones and even the thin high concentration filaments and very small difference appears between the 256 3 and 512 3 meshes.

For the velocity, the results are not highly affected by mesh refinement, as shown in Fig. 15(c). For the coarsest mesh the right motion is already captured, except near the vacuum regions. Same conclusions can be drawn from the energy fields. can reach such a regime using the segregation scalar field in such a case. In these conditions, it is obvious that the asymptotic regime has not been reached. This is even clearer when looking at the maximum of the extrapolated error (A.3) in time, see table A.3. When the convergence error for the mean total energy is globally lower than two percent, it is higher than 10% when looking at the segregation of the disperse field and reacher 60% for the most segregated Eq. B.5 shows a non-linear coupling between the reconstruction of each primitive variables. Therefore, ensuring realizability constraint on a specific variable may affect the other variables. In the following, we propose a sequential strategy that limits the need of a non-linear optimization to find the best choice:

1. First the density field n i is limited using a classical minmod limiter that ensures its positivity. To ensure the positivity of the determinant of the covariance matrix Σ, an additional procedure is used (not detailed here), in the same way as in [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] for the 2D case: if the computed slopes generate a negative determinant, an additional limitation is applied on the slopes, knowing that for zero slopes, realizability is ensured. where Σ * = Σ j -∆x 2 12 αD 2 u . To ensure the positivity of the determinant of Σ, a non-linear optimisation is still required, as in [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF].

Appendix B.3. Flux evaluation

The numerical flux evaluation is performed by means of a classical HLL approximate Riemann solver [START_REF] Harten | On upstream differencing and Godunovtype schemes for hyperbolic conservation laws[END_REF]: This standard flux can be shown to be realizability-preserving in the sense that, if S denotes the convex realizability constraints, then W n i-1 , W n i , W n i+1 ∈ S at time step n implies that the explicit first order Euler update

F HLL (W L , W R ) = 1 2 (F (W L ) + F (W R )) - 1
W n+1 i = W n i - ∆t ∆x F HLL W n i , W n i+1 -F HLL W n i-1 , W n i
also belongs to S under the classical CFL condition on the ratio ∆t ∆x for first order finite volume schemes [START_REF] Zhang | On positivity preserving high order Discontinuous Galerkin schemes for compressible euler equations on rectangular meshes[END_REF].

Then, by means of a strong stability preserving (SSP) 2-step Runge-Kutta method, we reach second order in time while preserving the realizability of the update.

3 / 2 ε

 32 characterizes the large eddies, where TKE is the mean turbulent kinetic energy. The initial parameters of the spectrum are C = 1.5, β = 5.2 and p 0 = 4[START_REF] Pope | Turbulent flows[END_REF], c L = 0.019, TKE = 0.00667, c η = 0.051, L = 4 and η = 0.027.

3. 2 . 1 2

 21 Characteristic time scale of the fluid The Stokes number based on the Kolmogorov time scale St K is classically used to characterize the occurrence of PTC, since it relies on the scales of the turbulence of the shortest characteristic time. This choice, although very common in DNS, can be tempered because the Kolmogorov time scale τ K = ν/ ε represents an average time scale. Indeed, ε = ν < 2s ij s ij > is obtained as an averaging over the spatial domain, where s ij = may exist in the turbulence. To illustrate this possibility, the distribution of the characteristic time scales of the fluid τ local = 1/ √ s ij s ij is plotted in Fig. 1. It shows the existence of short time scales that may generate PTC, even if St K is below the critical value of 1. Such a behaviour has already

Figure 1 :

 1 Figure 1: Distribution function of the local time scale of the gas velocity field for the FHIT at t = 40 in black. The red line identifies the Kolmogorov time scale τ k .

4. 1 .

 1 Number density field 4.1.1. Qualitative visualization The first level of validation consists in a qualitative comparison of the number density fields for three different Stokes numbers, see Fig. 2.For St K = 1, the spatial structure to be captured is challenging from the numerical point of view: we have to deal with very high concentration in thin structures close to large depletion zones. In fact, the closer to the critical Stokes number, the greater the probability to generate δ-shocks. At this particular Stokes number, AG and MK models provide the same global structure of the number density field, in agreement with the Lagrangian reference, despite the fact that Lagrangian solution presents thinner structures.For St K = 3, the underlying physics is expected to lead to significant PTC In this case, the MK model is not expected to capture the right dynamics of the particle-laden flow. This is clearly seen in Fig.2, where particles accumulations are highly overestimated by the MK model compared to the Lagrangian and AG results. On the other hand, the AG model still captures the global structures, their sizes and does not overestimate the accumulations. These conclusions still hold for St K = 8, for which the MK model still overestimates the level of particles accumulation. The Lagrangian reference shows a lower level of accumulation compared to St K = 3, and the AG model follows the same trend.

Figure 2 :

 2 Figure 2: Number density field at t = 40 on the 128 3 mesh for St K = 1 (top), St K = 3 (center) and St K = 8 (bottom), for the Lagrangian reference (left), the Mono-kinetic model (center) and the Anisotropic Gaussian model (right).

K = 8 Figure 3 :

 83 Figure 3: Time evolution of the segregation for the Lagrangian result (black solid line) and Eulerian results using 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (grey) meshes, for St K = 1, St K = 3 and St K = 8, for the MK model (dashed lines) and the AG model (dot-dashed lines).

Figure 4 :

 4 Figure 4: Number density PDF for the Lagrangian (black histogram), AG (dot-dashed lines)and MK (dashed lines) results on 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (gray) meshes for St K = 1. Arrows indicate increasing mesh sizes for Eulerian simulations.

Figure 5 :

 5 Figure 5: Number density PDF for the Lagrangian (black histogram), AG (dot-dashed lines) and MK (dashed lines) results on 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (gray) meshes for St K = 3. Arrows indicate increasing mesh sizes for Eulerian simulations.

Figure 6 :

 6 Figure 6: Number density PDF for the Lagrangian (black histogram), AG (dot-dashed lines) and MK (dashed lines) results on 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (gray) meshes for St K = 8. Arrows indicate increasing mesh sizes for Eulerian simulations.

K = 8 Figure 7 :

 87 Figure 7: PDF of the x-velocity for the Lagrangian (black solid line), AG (dot-dashed line) and MK (dashed line) results on 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (grey) meshes for St K = 1, St K = 3 and St K = 8.

Figure 8 :

 8 Figure 8: Statistics evolution with time for Lagrangian results (black solid line) and Eulerian results using 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (grey) meshes with projection on the same 64 3 mesh for St K = 1, for the MK model (dashed lines) and the AG model (dot-dashed lines). Arrows indicate increasing mesh sizes for Eulerian simulations.

Figure 9 :

 9 Figure 9: Statistics evolution with time for Lagrangian results (black solid line) and Eulerian results using 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (grey) meshes with projection on the same 64 3 mesh for St K = 3, for the MK model (dashed lines) and the AG model (dot-dashed lines). Arrows indicate increasing mesh sizes for Eulerian simulations.

Figure 10 :

 10 Figure 10: Statistics evolution with time for Lagrangian results (black solid line) and Eulerian results using 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (grey) meshes with projection on the same 64 3 mesh for St K = 8, for the MK model (dashed lines) and the AG model (dotdashed lines). Arrows indicate increasing mesh sizes for Eulerian simulations.

K = 8 Figure 11 :

 811 Figure 11: PDF of the pressure for the Lagrangian (black histogram) and AG results on 64 3 (green), 128 3 (blue), 256 3 (red) and 512 3 (grey) meshes for St K = 1, St K = 3 and St K = 8. Arrows indicate increasing mesh sizes for Eulerian simulations.

4. 4 .Figure 12 :Figure 13 :

 41213 Figure 12: Segregation (a), Total energy (b) and internal energy (c) as a function of the Stokes number for Lagrangian result and Anisotropic Gaussian results on the 128 3 and 512 3 meshes with projection on the same 64 3 mesh at t f . The vertical lines are: St K = 1 in black, St L = 1 in blue and St E = 1 in red.

Figure A. 15 :

 15 Figure A.15: From left to right: Lagrangian, MK on 64 3 , 128 3 , 256 3 and 512 3 meshes at t = 40 for St = 0.5

K = 8 Figure A. 16 :

 816 Figure A.16: Time evolution of the MTE for the Lagrangian result (black solid line), the AG results (dot-dashed lines) with 128 3 (grey), 256 3 (red) and 512 3 (blue) cells, and the estimate of the segregation at convergence (dashed line with crosses), for three Stokes number.

2 uiD 2 viDD

 22 ui D vi D ui D wi D ui D vi D vi D wi D ui D wi D vi D wi D σ11,i D σ12,i D σ13,i D σ12,i D σ22,i D σ32,i D σ13,i D σ23,i D σ33,i

2 .- Dn i ni ∆x 6 ,

 26 Next, the slopes for velocity are also limited to ensure the constraint of positive internal energy. This limitation reads:D u,i = 1 2 (sign (u i+1 -u i ) + sign (u i -u i-1 )) min   |u i+1 -u i | ∆x 1 |u i -u i-1 |

3 .

 3 Finally, the components of the covariance matrix Σ are limited in order to ensure the preservation of the positivity of its diagonal terms:diag (D Σ,i ) = 1 2 sign diag Σ * i+1 -diag (Σ * i ) + sign diag (Σ * i ) -diag Σ *

2 1 2

 21 |λ min | (W * -W L ) -|λ max | (W R -W * ) , (B.9)where the intermediate state W * isW * = λ min W L -λ max W R λ min -λ max -F (W L ) -F (W R ) λ min -λ max , (B.10)R and L stand respectively for the right and left side of the interface and λ min and λ max are respectively the slowest and fastest characteristic waves of the states on the left and right side of the interface.

Table 1 :

 1 Turbulence properties of the HIT.

	TKE	ε	τ K	η	τ L	τ eddy	L 11	Re
	0.00667 0.000635 1.2545 0.03535 4.7905 9.7864 0.65217 43.46

Table 2

 2

	for St K = 8.0, for a computation up to t = 80, on 512 intel Xeon E5-2690V3
	processors.
	Simul.	Lag		MK		AG
	Mesh	64 3 128 3 256 3 512 3 64 3 128 3 256 3	512 3
	Time 1h37m 33s	3m	32m 7h6m 54s	5m	57m 15h19m

Table 2 :

 2 Simulation time (wall time) for the particle solver using Lagrangian or Eulerian solvers. Simulations are performed on 512 Intel Xeon E5-2690V3 on the supercomputer Oc-

cigen (CINES, France).

The symmetric tensor outer product acts on a symmetric k-tensor and a symmetric ltensor by symmetrizing their (k + l)-order outer product tensor, see[START_REF] Levermore | The Gaussian moment closure for gas dynamics[END_REF].

For large enough discretization steps, we may be in a non-asymptotic regime and the ordre of the method can change with the step size as analyzed in[START_REF] Descombes | Analysis of operator splitting in the non-asymptotic regime for nonlinear reaction-diffusion equations. Application to the dynamics of premixed flames[END_REF].

Acknowledgements

The authors acknowledge SAFRAN-SNECMA for the financial support of the PhD of M. Sabat. This work was granted access to the HPC resources of CINES under the allocation x2016026172 made by GENCI (Grand Equipement National de Calcul Intensif). The authors also acknowledge IDRIS for the computation resources on Ada. The help of M. Boileau, J.C. Brändle de Motta, O. Thomine and L. Freret on Asphodele/Muses3D codes is gratefully acknowledged.

Appendix A. Numerical convergence of the Eulerian simulations

In order to compare the Eulerian models in the context of a refinement process, the results for MK and AG have been presented on the 64 3 , 128 3 , 256 3 and 512 3 meshes in the result section. Hence, we aim at separating the modeling influence from the numerical contributions on the capacity of our strategy in predicting the proper physical behavior of the flow. Ideally, we would be interested in the comparison of the models at numerical convergence, which is hard to reach.

In this Appendix, we wish to go one step further in the convergence analysis, and look at a qualitative assessment of the convergence of the different considered fields, as well as at a grid convergence index and regression tool in order to estimate a converged solution. meshes and if we suppose that we have reached the asymptotic convergence regime of the numerical scheme, we may postulate that the variable ε evolves with the mesh refinement as:

where ε * is the theoretical value of the quantity ε at convergence, C is a constant, h k is the characteristic size of mesh k, and p is an estimate of the order of convergence. If we are in the asymptotic regime 2 , these two constants, C, p are fixed and from the three equations (A.1) we can obtain the values of the constants as well as an estimate of the scalar field at convergence:

We thus get an estimate of the finest result in the asymptotic regime using:

Such an analysis relies heavily on the fact that we have reached the asymptotic regime for the numerical scheme. If not, it could lead to inaccurate orders and meaningless converged solutions.

In the following, this analysis is performed for the evolution of the mean total energy plotted and of the segregation in Figs. Appendix B. A Realizable Second-Order Dedicated Finite Volume

Scheme

Here, the realizable second-order MUSCL/HLL scheme [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] discussed in Sec. 2.4 is detailed. It is based on the MUSCL strategy [START_REF] Van Leer | Towards the ultimate conservative difference scheme V. A second order sequel to Godunov's method[END_REF], using a linear conservative reconstruction of the primitive variables, a realizable HLL scheme, and a

Strong-Stability-Preserving Runge-Kutta method. Multi-dimensional problems are solved by a dimensional splitting strategy, so we focus on x-direction numerical building block. For more details, one may refer to the work of [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF].

Appendix B.1. Realizable and conservative reconstruction

Realizability constraints on the moment set are expressed more easily on the primitive variables. Therefore, we apply the linear reconstruction on these variables, and we impose moment conservation in each cell to ensure realizability of the numerical scheme [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF].

It starts by the evaluation of local slopes by means of the neighboring mean values V:

Contrary to the classical MUSCL reconstruction [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF], the reconstruction of the primitive variables V = (n, u, Σ) T needs to be corrected in order to ensure the conservation of the conserved variables W = (n, nu, nE) T :

1 ∆x

where index i refers to the cell i. These latest constraints imply the following corrected primitive mean values:


Narjes Jomaa

Paolo Torrini

David Nowak

Gilles Grimaud

Samuel Hym Cristal

Proof-Oriented Design of a Separation Kernel with Minimal Trusted Computing Base

Keywords: protokernel, memory isolation, formal proof, Coq

The development of provably secure OS kernels represents a fundamental step in the creation of safe and secure systems. To this aim, we propose the notion of protokernel and an implementation -the Pip protokernel -as a separation kernel whose trusted computing base is reduced to its bare bones, essentially providing separation of tasks in memory, on top of which non-influence can be proved. This proof-oriented design allows us to formally prove separation properties on a concrete executable model very close to its automatically extracted C implementation. Our design is shown to be realistic as it can execute isolated instances of a realtime embedded system that has moreover been modified to isolate its own processes through the Pip services.

Introduction

The development of provably secure kernels addresses the trusted computing base (TCB) that has privileged access to the hardware, thus playing a fundamental role in the development of secure systems [START_REF] Alves-Foss | The MILS architecture for highassurance embedded systems[END_REF]. Making such development manageable and cost-effective remains a major challenge in formal methods even after major breakthroughs [KEH + 09, GSC + 16]. Memory management and isolation are at the core of the functionalities of separation kernels [START_REF] Rushby | Design and verification of secure systems[END_REF][START_REF] Zhao | High-assurance separation kernels: a survey on formal methods[END_REF], designed to ensure that users can access resources and communicate with each other only according to a given security policy.

The centrality of memory management and access control provides our motivation for introducing the notion of protokernel, as a minimal separation kernel that in fact provides only services related to virtual memory management and context switching between applications. In this paper, we present the development and verification of the Pip protokernel [START_REF]The Pip Development Team[END_REF], relying on a comparatively lightweight approach, yet providing an efficient, real-world application. We use Coq, a proof assistant based on the Calculus of Constructions [START_REF] Bertot | Interactive Theorem Proving and Program Development[END_REF], to develop the monadic, executable model of the services on top of a hardware abstraction layer, and we prove security properties based on a low-level memory isolation criterion, using Hoare logic [START_REF]The Pip Development Team[END_REF]. The model is then automatically translated to a small fragment of C, using an Haskell-implemented translator [START_REF] Hym | Digger[END_REF].

Pip has been developed for industrial application as part of a European project1 . Its development relies on a proof-oriented design, narrowing the TCB needed to implement the kernel services in a machine independent way on top of the hardware abstraction layer. This design,

INTRODUCTION

minimalistic yet efficient, is the result of a collaborative effort with our colleagues from the operating system community (cf. [BGI18, YGG + 18, BJY + 18] for benchmarks and more details) to support on-demand secure isolation [The17a]. Crucially, Pip supports efficient memory management and a dynamic notion of memory separation based on a hierarchical TCB architecture, allowing for isolated partitions to be created and removed after initialisation.

Contributions

In this paper, we report the formal verification viewpoint of the development of the Pip protokernel. We present: 1) a Coq model that includes the relied-upon hardware components (memory and MMU), an executable specification of the memory manager supporting hierarchical TCB, and a highly modular API consisting of ten services; 2) a low-level specification of memory isolation; 3) an abstract information flow model which allows us to prove a non-influence result; 4) a verification methodology based on a Hoare logic on top of an enhanced state monad, designed to carry out well-structured proofs of memory isolation on the executable model; 5) the application of this methodology to the verification of three of the API services.

Related work

The development of secure kernels has been a long-standing goal in the formal methods community [START_REF] Zhao | High-assurance separation kernels: a survey on formal methods[END_REF]. Memory separation in its basic form can be characterised in terms of access control. Stronger notions of information flow separation have been introduced with the noninterference property [START_REF] Goguen | Security policies and security models[END_REF], ensuring that events are not leaked [START_REF] Rushby | Noninterference, transitiviy, and channel-control security policies[END_REF], and the non-leakage property [vO04], ensuring that data are not leaked. Early work on access control policies for the UCLA security kernel [START_REF] Walker | Specification and verification of the UCLA Unix security kernel[END_REF] focused on the verification of an abstract model, stressing the much higher cost of verifying executable code. The notion of separation kernel, originally introduced by Rushby [START_REF] Rushby | Design and verification of secure systems[END_REF] to ensure information flow separation, is at the foundations of security-oriented kernels based on the MILS architecture [START_REF] Alves-Foss | The MILS architecture for highassurance embedded systems[END_REF] and of safety-oriented ones [START_REF] Zhao | Reasoning about information flow security of separation kernels with channel-based communication[END_REF].

SeL4 [KEH + 09] is a microkernel of the L4 family formally verified in Isabelle-HOL. It provides virtual memory management, thread scheduling, inter-process communication and access control. High-level verification is carried out using Hoare logic on a abstract model. An executable specification that refines the abstract model is obtained by automated translation from a Haskell prototype. The manually written C implementation is proved to refine the executable specification using an automatic translation of C to Isabelle based on the semantics in [START_REF] Tuch | Types, Bytes, and Separation Logic[END_REF]. The refinement proofs are functional correctness ones, thus quite big, ensuring that the abstract model provides a complete specification of the C behaviour. Security properties are proved relying on the abstract model and functional correctness [MMB + 13].

PROSPER is a separation kernel that supports dataflow functionalities of the ARMv7 processor [DGK + 13], including inter-partition communication, scheduling, and in an extended version [START_REF] Schwarz | Formal verification of secure user mode device execution with DMA[END_REF] direct memory access. Security properties have been proved in HOL4. Bisimulation is used to carry higher-level proofs, based on an abstract model including communication channels, over to a concrete one based on an intermediate language that represents the ARM instruction set [START_REF] Fox | A trustworthy monadic formalization of the ARMv7 instruction set[END_REF].

CertiKOS [GSC + 16] provides a certified programming method to formally develop kernels in Coq, targeting an extension of CompCert Clight and assembly code [START_REF] Leroy | Mechanized semantics for the Clight subset of the C language[END_REF]. The specification approach is based on certified abstraction layers, corresponding to Clight external functions, enhanced with a notion of abstract state, to allow for functional correctness specifications in the code. CertiKOS supports fine-grained layering of programming and verification, relying on a notion of contextual refinement [GKR + 15], and it has also been used to prove noninterference results [START_REF] Costanzo | End-to-end verification of information-flow security for C and assembly programs[END_REF].

In [START_REF] Baumann | Proving memory separation in a microkernel by code level verification[END_REF] the PikeOS memory manager is verified at the C source code level using an automated static verifier, proving that initialisation-time partitioning of memory is preserved at runtime. In [START_REF] Barthe | Formally verifying isolation and availability in an idealized model of virtualization[END_REF] a Xen-based paravirtualization-style hypervisor is abstractly modelled in Coq and information flow separation is proved. In [START_REF] Zhao | Reasoning about information flow security of separation kernels with channel-based communication[END_REF] Isabelle-HOL is used to model abstractly the ARINC-653 channel-based inter-partition communication for safety-oriented separation kernels. In [BCG + 16] an abstract model of the Nova kernel is formalised in Coq, data security properties are proved, and program extraction is used to generate testing code. In [START_REF] Sanan | Separation kernel verification: the XtratuM case study[END_REF] an abstract security model is formalised in Isabelle-HOL to verify properties of the XtratuM kernel using refinements.

Pip

Our system relies on a notion of memory separation, on top of which noninterference can be proved for isolated partitions. Inter-partition communication and scheduling are not part of Pip, and therefore communication policies are not part of the kernel TCB. On the other hand, Pip allows for partition reconfiguration at runtime, and tree-like partitioning can naturally support any user TCB hierarchy. The security proof we provide does not rely on any specific architecture, and has been the driving factor in the development. The verification of memory separation is comparatively lightweight, and this has been achieved by avoiding to go through a full functional correctness proof. From the development point of view, Pip relies on a separation of concerns between modelling and verification on one side and translation to C on the other. Nonetheless, unlike work based on abstract models, Pip provides an executable model of its services, written in a shallow embedding of a small fragment of C. Indeed, our Coq development is based on a direct formalisation of the executable model, relying on consistency properties to capture the abstract requirement in a bottom-up way, following the approach already used in [JNGH18] for a simpler system with a flat memory model. The translation of the service layer to C source code is fully automated, and its verification, on which the correctness of the model does not depend, is under way as an independent workload (see [START_REF] Torrini | Formalising executable specifications of low-level systems[END_REF] for more discussion of this point).

Outline

In Section 2 we present the design of Pip. In Section 3 we present the Coq model, including hardware components and partitioning manager. In Section 4 we present the security properties and the abstract information flow model. In Section 5 we present the isolation proofs and the underlying methodology. Conclusions and further work are discussed in Section 6.

Proof-oriented design

Pip is a kind of security kernel conveniently described as a protokernel [The17b, BJY + 18, BIG17], designed to minimize the TCB. The goal of such minimization is twofold: to reduce the attack surface, and to increase the feasibility of a formal proof. Pip only provides the hardware management that suffices to enforce security in the form of memory separation: in this sense, we speak of proof-oriented design, emphasising the feedback of theorem-proving in the iterative development of the system. Like exokernels [START_REF] Engler | Exokernel: An Operating System Architecture for Application-level Resource Management[END_REF], Pip does not offer any higher-level abstraction: in particular, it offers neither threads, nor processes, nor a file-system, nor a network stack. Pip provides a runtime environment that is a kind of virtual machine and that we call a partition. On top of that we can implement any high-level abstraction. Each partition allows configuring a virtual memory management unit (MMU) and managing virtual interrupt request (IRQ) via dedicated Pip services that form the Pip API. Using these configuration services, it is possible to port on Pip a conventional kernel (e.g. Linux) as well as an embedded system (e.g. FreeRTOS) or even a hypervisor (e.g. Xen). Little effort is needed to adapt one of them to run as a partition (see [BJY + 18] for details, including more than satisfactory benchmarks), relying on the paravirtualization approach [BDF + 03].

Pip can be better described as protokernel rather than hypervisor because it does not include multiplexing management, thus reducing the kernel TCB to its bare bone. Indeed, in a sense Pip runs a single primary partition. This partition, created at boot time, represents the whole and only real machine. Pip assigns to it the entire address space (except for the memory used by Pip itself) and the whole CPU time (except for the CPU time spent by Pip to handle real IRQs).

Nonetheless, Pip supports a hierarchical partitioning model, allowing each partition to create its own subpartitions and split part of its own address space between them. This recursively results in the creation of a partition tree whose root is the primary one (henceforth the root partition). The management of the partition tree takes place according to a recursive scheme. Only one partition can be executing at each time in a monocore system, but Pip does neither multiplex interrupt requests nor share CPU time between distinct partitions. Pip only deals with software interrupts that are system calls to its own API. Other software interrupts are forwarded to the parent of the caller. On the other hand, all the hardware interrupts are forwarded to the root. In this way, Pip delegates multiplexing (i.e. determining the control flow between partitions) to the root, which can implement any possible multiplexing policy, recursively allowing each partition to share CPU time between its children.

Executing multiplexing in user mode allows us to cut drastically the size of the kernel TCB. On the other hand, the hierarchical architecture allows each partition to rely on the environment managed by its parent and to split its environment between its children. This naturally induces a hierarchical characterisation of the TCB associated with each partition, which can be expressed in terms of a simple transitive policy: each partition only relies on its ancestors and the kernel. The primary purpose of this design is to allow for a cost-effective verification that this policy is enforced, by means of a memory isolation proof.

The memory isolation property ensured by Pip is based on hardware functionalities that control physical memory accesses. Partitions can not use physical addresses: in order to access information they can only use virtual addresses that will be converted to physical ones using the Memory Management Unit. Each partition has an MMU configuration, and on every user access, the MMU is invoked to decide on whether authorising or rejecting the access. These decisions are based on the configuration pages, called page tables or indirections, written by the kernel so as to allow translation of any given virtual address to the corresponding physical address. In order to guarantee security, page tables should never be accessible by any entity except the kernel. Pip satisfies this requirement by construction.

Partitioning model

The physical memory of a computer is split into fixed-size chunks called pages which can be allocated, accessed and used to store information. Whenever a partition is created, Pip associates it with a page which we call its partition descriptor (PD), denoted by an identifier referred to as PDI.

The memory is split between primary kernel pages, solely accessible by the kernel, and user pages, initially assigned to the root partition and accessible by it. When a partition creates a child in the course of its execution, it allocates some of its accessible pages exclusively for that child; in particular, the executing partition hands over to the kernel the sole accessibility of some of the newly allocated pages, to be used for the configuration of the newly created partition; the remaining of those pages, already assigned to the parent are also assigned to the child. In line with the hierarchical model, the executing partition can always delete a child (and with it its descendants). When this happens, the parent regains accessibility to the configuration pages of the descendants. Notice that in general, the pages a partition can access are a subset of those assigned to it. We refer to the pages assigned to a partition together with those used for its configuration as the pages allocated for that partition.

The API of Pip is constituted of ten services that can be called by the executing partition, eight of which are for managing its memory space. createPartition adds a new child to the caller. deletePartition deletes a child and gives all the pages allocated for it back to the caller. addVaddr lends physical memory to a child. removeVaddr removes a page from a child. mappedInChild returns the child to which a given page is assigned. prepare gives pages to the kernel to manage a child's configuration. pageCount computes the number of pages required to configure a child. collect gives pages lent to the kernel back to the calling partition. The remaining two services are for handling IRQs. dispatch notifies a partition about a given interrupt. resume restores the context of a previously interrupted partition.

Fig. 1 serves to illustrate partition management in Pip: the available memory makes up the root partition P root . User code running in any partition can use the memory management services exposed by Pip to create child partitions, here P 1 and P 2 , lending some of its accessible pages to each of them. For instance, code running in P 1 can create P 1.1 using available memory in P 1 .

While the parent partition can read and write in the memory accessible by its children (we call this property vertical sharing), it cannot access anymore to the pages given to Pip (we obviously want to prevent a partition from messing up Pip data structures, we call this property kernel isolation). Sibling partitions (i.e. partitions that have the same parent) cannot access each other's memory (we call this property horizontal isolation). If siblings want to communicate, they cannot do it through shared memory. Instead, their parent can use the services of Pip to implement a communication protocol, for instance by flipping a memory page between the siblings.

For example, as result of horizontal isolation, the partition P 1 in Figure 1 Figure 1: An example of partition tree partition P 2 . Thus, code running in the partition P 1 cannot access the memory in the partition P 2 . Moreover, since partitions P 1.1 and P 1.2 are included in the partition P 1 , they are also disjoint from the partition P 2 . However their parent partition P 1 can access their memory and the memory of their descendant partitions P 1.1.1 , . . . , P 1.2.2 as result of vertical sharing.

The executable specification

The development of Pip is based on a layered model formalised in Coq, as shown in Fig. 2. At the top level, the service layer provides a source-code level executable specification of architectureindependent system management functions. The service layer consists of the kernel service API that allows for context switching and partition tree management (PTM) which is the core of the Pip engine. The service layer is implemented as an algorithmic model, using the shallow embedding of a C-like fragment based on a monadic encoding, which can be automatically translated to C code. The service layer is built on top of a hardware abstraction layer (HAL), including memory abstraction (MAL) and interrupt abstraction (IAL). The MAL itself can be structured into two layers: hardware memory abstraction (HMAL) and monadic memory abstraction (MMAL).

The HMAL provides an abstract model of the physical memory and the MMU. The IAL is also structured into two layers: hardware interrupts abstraction layer (HIAL) and monadic interrupts abstraction layer (MIAL) required to perform context switching. The MMAL and MIAL provide high-level but executable specifications of low-level, architecture-dependent functions which have been manually implemented in C and assembly. The executable model that can be actually run in Coq thus includes the algorithmic model of the service layer as well as the MMAL and MIAL.

Hardware memory abstraction layer

Memory size is determined in our model by two architecture-dependent parameters, the positive integers pageSize for memory page size and nmbPages for the number of pages. Pages are pointed to by page identifiers of type page, which are natural numbers less than nmbPages, modelled in Coq using dependent records and defined as {p :> nat; Hp: p < nmbPages}.

We use 0 as default value for the identifier that points to the null page. The physical address of a memory cell consists of a page identifier and a position in that page, given as an offset value called index. It is modelled by type paddr and defined as (page * index). Also indices are typed as a subset of positive integers index and defined as {i:> nat; Hi: i< pageSize}.

In general, the physical memory state of a computer can be modelled as an association list which maps physical addresses to values. However, our model is completely abstract with respect to the content of the user space. We only need to model the partition tree, and thus we only provide the content of the kernel-owned pages, from which the tree structure can be computed. The part of the hardware state that is relevant to Pip consists of the partition descriptor of the currently executing partition and of the parts of the physical memory where Pip stores its own data (essentially related to the configuration and management of the partition tree), as expressed by the record type state which is defined as: Here Pentry stands for physical entry, Ventry for virtual entry, and vaddr for virtual address. We prove that memory management is well-typed with respect to these value types, as an invariant consistency property of our model. Physical entries (PTEs) make it possible to associate page addresses with information concerning whether each page is assigned to the current partition (present) and whether it is owned by it (accessible). A physical entry is modelled by

SH 1 SH 2 L MMU PD kernel pages

Partition tree management

The configuration of each partition (other than the root) is stored in the memory allocated for it by its parent partition, thus assigned to the parent, though in the part of it that is only accessible to the kernel, and it is defined by four entities. The first and principal one is the MMU configuration, the other ones include two tree-like structures and a list-like one. Each partition descriptor (PD) is a page that contains the addresses of four configuration entities, as illustrated by Fig. 3. The PD together with its entities forms what we call a configuration node.

The MMU configuration of each partition has a hierarchical structure that can be abstractly described as a graph where nodes are pages, each being either a page table or a terminal page, and each terminal page being either a page assigned to the configured partition or the null page. Each entry in a page table points to another page, thus providing the arcs in the graph, while providing accessibility information. The MMU structure has two fundamental properties which have been verified as invariant consistency properties of our model. The first property ensures that the one and only physical entry marked as not present corresponds to the null page -in fact, the MMU configuration only contains pages that are assigned to the parent partition. The second property ensures that the subgraph formed by non-null pages has the structure of a tree, which we call the MMU tree. The tree depth is fixed by the architecture-dependent natural parameter levelNum. Virtual addresses are modelled by type vaddr as lists of indices of length levelNum+1 and defined as:

Record vaddr := {va :> list index; Hva: length va = levelNum + 1}.

Each virtual address is translated either to the null address or to a physical address in a leaf, by interpreting each index in the list as offset in the page table at the corresponding level in the MMU tree. The last index in the list is the offset of the corresponding physical address, which may be accessible or not to the partition, depending on whether the page is owned or not, as specified by the last PTE. Crucially, in Pip only the kernel can access directly physical addresses. Virtual addresses are the objects provided by the MMU to user applications as indirect references to the physical memory. Indirect access relies on the mmu_translate function, which models the hardware mapping of virtual addresses to physical ones, and which has been proved to respect the access policy.

Intuitively, non-null virtual addresses correspond to maximal branches in the MMU tree (plus the final offset), and they include those pointing to the descendants in the partition tree. In fact, the content of a virtual address (i.e. the value of the physical address it gets translated to) can be the PDI of another partition (indeed, this is the only kind of content Pip cares for). In this way configuration nodes can be linked together. The resulting structure is a graph which we prove to be a tree, and thus to represent the partition tree, as an invariant consistency property of our model.

The partition tree management makes use of two auxiliary configuration entities, called shadows, that mirror the corresponding MMU tree. The first shadow is used to find out which pages are assigned to children, an information that is needed to ensure horizontal isolation. It uses the type Ventry of virtual entries which is defined using the following record:

Record Ventry : Type:={pd: bool; va: vaddr}.

The flag pd indicates if the associated physical address is the PDI of a child partition. The second shadow is used to associate each partition descriptor to the virtual address it has in its parent partition. The fourth auxiliary entity is a linked list used to remember which pages have been lent to the kernel. The information in the second shadow and in the linked list is needed to process efficiently the deletion of a partition.

The monad

In order to implement the Pip services imperatively in Coq as sequential programs, we rely on a monadic approach. Monads [Mog91, W.92] allow threading effects through computations, and can be used to interpret semantically imperative languages in purely functional ones. In our case, side effects include the state, corresponding to the hardware state, and a notion of exception corresponding to the possibility of undefined behaviours that we actually prove to never happen. For example, any attempt by Pip to access a physical address that is not defined in the memory state would result in an undefined behaviour.

We define our monad (Low Level Interface monad) as an abstract datatype LLI A that wraps together hardware state and undefined behaviours, following [Mog91, W.92].

Definition LLI (A : Type) : Type := state ! result (A * state).

Here result is an inductive type with two constructors: the first one corresponds to a result of type A, and the second one to an undefined behaviour. We write perform x := m in e for the monadic binding operation bind m (fun x ) e), which can also be written m ;; e when x is not used in e, to represent sequential composition of the stateful actions m and e. The state monad ensures we can define the usual stateful functions, e.g. get and put to read and update the state. We use the monadic functions to define the specific ones that form the MAL, matching the corresponding architecture-dependent implementations in C and assembly. Crucially, the Pip services are implemented so to access the state only through the MMAL functions.

Monadic memory abstraction layer

The MMAL specifies a C and assembly library for the architecture-dependent part of Pip, which consists of kernel atomic operations, such as reading and writing values in physical memory, and performing simple computations and comparisons on them. This includes physical memory access operations needed to configure partitions, virtual memory activation and a set of auxiliary operations. It mainly consists of bitwise operations of constant computational complexity. The functions defined in the MMAL can be organized along four categories. The first one is HAL_write, allowing the kernel to write values into physical memory through physical addresses and including the primitives writeVirtual, writePhysical, writeVirEntry, writePhyEntry, writeAccessible, writePresent, writePDflag, writeIndex. For example, the writeVirtual primitive has the following definition: It stores a virtual address at the physical address given by a page and an index. The monadic operator modify is used to modify the association list that represents the memory. This requires that there is a value of the right type at the given physical address. If the condition is not met, the result is an undefined behaviour. The third category (named HAL_op) allows the kernel to perform simple operations, including e.g. successor and equality comparison, on each value type. The last category (named HAL_const) is about accessing global constants. Notice that since Coq types are more discriminating than C ones, different MMAL functions may be actually implemented as the same C function.

Linguistic aspects

The program extraction capabilities of Coq made it possible to extract Haskell code which we used to test our model. However, such code does not meet the runtime requirements of a realistic kernel. The cost of the runtime environment needed by functional programs, and particularly garbage collection, stand in the way of a reliable runtime behaviour. Therefore, we rely on a translation to C. The monadic code we use in our executable specification corresponds to a quite simple imperative language: essentially, a first-order sequential language with call-by-value, primitive recursion and mutable references. From a low-level point of view, all the datastructures used by Pip, including trees, can be treated as linked lists. The abstract memory model in Coq matches closely this low-level characterisation: Pip datastructures are represented using an encoding of linked lists, given in terms of lists and access to the monadic state, rather than by using specific inductive datatypes, and then by proving invariant consistency properties as appropriate to ensure that the low-level representation is correct.

In fact, our monadic code corresponds to a shallow embedding of the denotational semantics of a language, formally described in [START_REF] Torrini | Formalising executable specifications of low-level systems[END_REF], which matches closely a comparatively small fragment of C (we use no arrays, neither structures nor unions, no loop instructions, no pointer arithmetics). We use a Haskell-implemented tool called Digger [START_REF] Hym | Digger[END_REF] to generate automatically C code for the service layer from the abstract syntax of the monadic Coq code. The code generated by Digger from the executable model can then be compiled (we currently use GCC) together with the manually implemented MAL functions.

Relying on a low-level modelling of data structures has a drawback in the development of recursive functions based on such structures: we cannot rely on structural recursion as Coq would provide for the corresponding inductive type. This can be a problem already at the level of function definition, as Coq requires that we ensure termination. Nonetheless, we can easily deal with this issue, relying on the fact that Pip functions only require recursion bounded by parameters of the hardware architecture such as the size of a memory page. However, this means we need to prove that fuel suffices, for each top-level function call.

Security properties

Security in Pip relies on the MMU functionalities that control user access to physical memory, as user access to physical memory is only allowed through the MMU. Proved consistency properties ensure, for instance, that the MMU is configured correctly, consistently with the partition tree, and that memory access through the MMU is consistent with the access policy encoded in the configuration pages. Such properties are needed to prove the actual security properties of Pip. We now introduce the main properties that define the Pip partitioning model, provably preserved by the API services: horizontal isolation, kernel isolation, and vertical sharing.

In a tree, two nodes are unrelated whenever neither of them is a descendant of the other one. With respect to the partition tree, horizontal isolation means that unrelated partitions cannot access each other's memory. This is the case in particular for sibling partitions.

Definition 1 (HI) Horizontal isolation holds for a state s whenever for all partition descriptors m 0 , m 1 , m 2 2 partition tree(s), with m 1 , m 2 2 children(m 0 , s) and pdi(m 1 ) 6 = pdi(m 2 ), it holds that

allocated pages(m 1 , s) \ allocated pages(m 2 , s) = / 0.
where partition tree is the list of partition identifiers. We use consistency properties to ensure that partitions are organized into a tree-like structure.

Kernel isolation means that no partition can access the pages owned by the kernel.

Definition 2 (KI) Kernel isolation holds for a state s whenever for each partition descriptor m 2 partition tree(s), it holds that accessible pages(m, s) \ kernel owned pages(s) = / 0.

Vertical sharing means that all the pages allocated for a partition are included in the pages assigned to its ancestors, and therefore, owing to the access policy implemented by the MMU, a partition has read and write access on all the memory accessible by its descendants.

Definition 3 (VS) Vertical sharing holds for a state s whenever for all partitions m 0 , m 1 2 partition tree(s), with m 1 2 children(m 0 , s), it holds that allocated pages(m 1 , s) ⇢ assigned pages(m 0 , s).

Relying on these invariants we can prove formally a non-influence property for isolated partitions [START_REF]The Pip Development Team[END_REF] (i.e. excluding inter-partition communication), relying on an abstract information flow model. We call P-machine a state machine at the level of physical memory, where transition steps are specified by the type pstep (ls : list page) : state ! state . Here ls is the list of pages which can be accessed in the execution of the physical action associated with the step. We can define state equivalence with respect to a list of pages (where select returns the value stored at a location).

s 1 ⇠ ps s 2 := 8(p : page) (i : index), p 2 ps ! select s 1 p i = select s 2 p i
We write s 1 ⇠ s 2 for s 1 ⇠ ps s 2 when ps are all the pages, and s 1 ⇠ p s 2 for s 1 ⇠ [p] s 2 . We can specify the access requirement intended by pstep in terms of two conditions. First, the read condition states that the step depends only on the values read at locations in ps (i.e. the physical action depends only on its accessible locations).

readCond (p : pstep) : Type := 8(ps : list page) (s 1 s 2 : state), s 1 ⇠ ps s 2 ! p ps s 1 ⇠ ps p ps s 2

Second, the write condition states that the only locations possibly affected by the step are those in ps (i.e. the physical action only affects its accessible locations).

writeCond (p : pstep) : Type := 8(ps : list page) (s : state) (p : page), ¬ p 2 ps ! s ⇠ p p ps s

Any state machine can be modelled as a P-machine. Runs can be modelled as action sequences.

We can now introduce V-machines as virtual-level state machines which lift our information flow specification at the level of virtual memory. We can define the notion of virtual step by lifting a physical step, relying on the MMU translation function that represents virtual memory access in the HAL. Here m is the executing partition and vs is the list of virtual addresses accessed in the execution of the virtual action associated with the step.

vstep (p : page) (p : pstep) (vs : list vaddr) (s : state) : state := let ps := map (mmu translate p) vs in p ps s

The fact that the MMU is configured to return physical addresses in the memory owned by the executing partition, together with the access requirement on the P-machine steps, makes it straightforward in Coq to prove non-influence (a notion defined in [vO04]), with respect to a policy that excludes communication between separated partitions. We write @m to denote all the pages that are in a partition with address m. We write s 1 ⇠ m!ps s 2 whenever s 1 ⇠ ps s 2 , the executing PDI is m in both states, and s 1 ⇠ @m s 2 . We prove stepwise non-leakage, defined following [vO04].

s 1 ⇠ a!as s 2 ! 8pt vs, vstep a pt vs s 1 ⇠ as vstep a pt vs s 2 This property can be easily extended to action sequences. Similarly we can prove a noninterferencerelated property (called local-respect in [START_REF] Rushby | Noninterference, transitiviy, and channel-control security policies[END_REF]).

¬ a 2 @a x ! 8pt vs, vstep a x pt vs s ⇠ a s

This, together with non-leakage suffices to prove classical noninterference and non-influence [START_REF] Rushby | Noninterference, transitiviy, and channel-control security policies[END_REF][START_REF] Oheimb | Information Flow Control Revisited: Noninfluence = Noninterference + Nonleakage[END_REF]. Our proof shows how the low-level invariants (HI, KI and VS) can be used to prove a higher-level result. The extent of the result presented here is comparatively limited, as we consider a definition of virtual step which is very strict, ruling out the mechanism we need for inter-partition communication (i.e. interrupts). On the other hand, this model could be extended to deal with given user-level communication policies, by a reformulation of the write condition that takes interrupts into account.

The verification approach

Hoare logic on top of a state monad

Our verification targets the properties KI, HI and VS, defined in Section 4, as well as the consistency properties (which we will denote jointly by C in the following), with the goal of ensuring that each system call to an API service preserves all of them. Thus our proofs consist mainly of reasoning about invariants. We do this by using a Hoare logic which we define on top of the LLI monad. We use the syntax {{ P }} m {{ Q }} for our Hoare triples, defined as follows. Here the unary predicate P gives the precondition and the binary predicate Q the postcondition of running the computation m. If performing m yields an undefined behavior, the triple does not hold, thus ensuring a basic functional correctness property: the services of Pip never lead to an undefined behavior. We created a library of general Hoare logic rules to reason about monadic programs. We also created a specific library of Hoare triples (in general, either weakest preconditions or strongest postconditions) for the MAL primitives that sequentially form each of the API service. For example, the following has the weakest precondition of readVirtual as precondition. Our use of Hoare logic relies on the assumption that each system call is atomic: it is required that interrupts are blocked during a system call. This is one of the reasons the system calls were designed to be elementary. Moreover, our Hoare logic deals only with the global state, therefore it is required that in multicore architectures at most one core is executing a system call at any time. Concerning the hardware, we naturally assume that the physical memory behaves correctly, as specified (i.e. essentially it is not volatile), and that so does the MMU component. We also assume that the non-algorithmic MMAL is correctly implemented in C and assembly, and that the system was booted in an isolated and consistent state (as our proof is essentially a preservation one).

The code of Pip has a sequential character. Thus, verification can proceed backward, relying on weakest precondition triples, or forward relying on strongest postcondition ones. The two approaches are actually equivalent in our case, since we can ensure termination and our triples rule out undefined behaviours. In our concrete proofs we found it more convenient to move forward, following the actual execution flow.

The separation proof for addVaddr

In the case of AddVaddr, the top level invariant is the following:

Lemma addVaddrInvariant (src dst child: vaddr): {{fun s ) HI s ^KI s ^VS s ^C s}} addVaddr src dst child {{fun _ s ) HI s ^KI s ^VS s ^C s }}.

The AddVaddr service is called by the executing partition p 0 to assign one of its pages (further down denoted by srcp) to one of its children (child). Here src is the virtual address of srcp in the address space of p 0 , whereas dst is the virtual address that gets associated with srcp in the address space of child. The first part of the service consists in making some checks, before making any update to the state. These checks include consistency ones needed to avoid undefined behaviours, and security ones on each of the given parameters. If any of the security checks fails, the service will abort leaving the state unchanged. To ensure that child is actually a child of p 0 , AddVaddr starts by checking if the PDI of child is included in the first shadow of p 0 . Then it makes sure that srcp is already assigned to p 0 and accessible by it, checking the control bits present and accessible of the corresponding table entry in the MMU tree of p 0 . It checks that srcp is not already assigned to any child using the flag pd stored in the first shadow of p 0 . Finally it verifies that dst is not already associated to a physical page by going through the MMU tree of child. This prevents overwriting data. These properties gets propagated through the sequential execution until the last part, where state updates take place. In fact, it is part of our design of the Pip services to carry out state updates in a way that maximises invariant propagation without compromising efficiency. In the case of AddVaddr state updates are confined to the last three instructions, for which the following triple holds: writeVirtual modifies the second shadow of child, by storing src at dst (thus ensuring that the page can be efficiently located in the parent address space when it needs to get revoked).

writeVirEntry starts the process of assigning the page associated with src to the partition associated with child, by storing the address child at src in the first shadow of p 0 . writePhyEntry concludes the process, by associating the physical page srcp to dst in the MMU tree of child. As a result of relying on comparatively flat representations (e.g. our trees are represented using lists), a significant part of the proof consists of ensuring that consistency properties are preserved through sequential updates. This is particularly the case for the consistency of the first shadow, a data structure which has been specifically put in place to support isolation, by ensuring that a page cannot be assigned to distinct children. Our policy of ordering steps to maximise propagation pays with respect to several properties: for example accessibleChildPageIsAccessibleIntoParent,

CONCLUSIONS AND FURTHER WORK

requiring that each page assigned to a partition has a back-pointer in its second shadow, is preserved by the fact of performing writeVirtual first. However, there are some propagated properties that get temporarily broken and need to be patched up with weaker ones. For example, the property notShared srcp s, stating that srcp is not marked as shared according to the kernel information stored at p 0 , is no more valid after executing writeVirEntry. Luckily, a weaker property that here we denote H holds. It states that srcp is not assigned to any child according to MMU configuration, and suffices to prove the final triple: In particular, proving HI s (as by definition 1) involves showing that srcp is assigned to child and not assigned to any distinct sibling partition. Proving KI s (as by definition 2) involves showing that srcp is not a kernel page. Finally proving VS s (as by definition 3) involves showing that srcp is assigned to the parent partition of child.

Overview of the development and the verification

The top level executable specification includes about 1300 lines of Gallina. We estimate 3 months for kernel conception and 9 months for its development. We have currently verified three services: createPartition, addVaddr and mappedInChild. We started from the proof of createPartition which introduced the majority of the consistency properties resulting in about 60K lines of proof as detailed in table 1. The other two services were proved relying largely on lemmas already proved for createPartition and resulting in about 20K additional lines of proof. The current proof has required about one person-year of verification time. Along the proof, we made several changes to the code of the kernel. Most of time it was about fixing bugs such as adding some missing checks or reordering instructions to simplify the proof.

Invariants

Lines of proof duration createPartition (⇡ 300loc)

⇡ 60000 ⇡ 10 months createPartition + addVaddr (⇡ 110loc) ⇡ 78000 ⇡ 2 monthes createPartition + addVaddr + mappedInChild (⇡ 40loc) ⇡ 78300 ⇡ 4 hours ensured by the kernel is memory separation. Pip was designed with special attention to increasing the feasibility of its verification. In this sense, not only the minimisation of the TCB, but also the modular design of the API, the structuring of the configuration information, and the modelling economy of the monadic encoding have played a role in allowing for comparatively simple proofs.

Applications

In order to show that our proof-oriented design is realistic, the real-time embedded system FreeR-TOS and Linux 4.10.4 have been ported to Pip [YGG + 18, BJY + 18]. The porting consists essentially in removing privileged instructions and other operations, replacing them with system calls to Pip, thus allowing for isolating tasks by running them in isolated partitions. Pip has also been extended to support different multi-core configurations [START_REF] Bergougnoux | Porting the Pip proto-kernel's model to multi-core environments[END_REF], and this has been possible without changing its current implementation. The modification of existing proofs on the mono-core version to account this extension is still under study.

Translations to C and DEC

A verified translation to C of the service layer is currently being implemented, following one of the work plans discussed in [START_REF] Torrini | Formalising executable specifications of low-level systems[END_REF], relying on a small intermediate language called DEC [TN17] which has been formalised as a deep embedding to allow for a definition in Coq of a translation to CompCert C [START_REF] Leroy | Mechanized semantics for the Clight subset of the C language[END_REF]. The Digger translator [START_REF] Hym | Digger[END_REF], which works on the AST of the shallow embedding and is currently used to obtain uncertified C code, can also be used to translate automatically the service layer to DEC. Ongoing work is focusing on proving the adequacy of the deep embedding with respect to the shallow one. In [START_REF] Torrini | Formalising executable specifications of low-level systems[END_REF] the deep embedding has also been used to experiment with syntax-driven semi-automation of Hoare logic proofs.

Ongoing and future work

Concerning the formal track, the verification of the remaining services is currently under way, and so is the verification of the translation to C. Future work could involve, on one hand, extending the information flow model to deal with user-level interpartition communication, and on the other hand, extending verification to the architecture-dependent part.

Figure 2 :

 2 Figure 2: The design of Pip

Figure 3 :

 3 Figure 3: Partition tree configuration

  Inductive result (A : Type) : Type := |val: A ! result A |undef: nat ! state ! result A. We define the monadic operations ret and bind, which can be easily proved to satisfy the monadic laws [Mog91, W.92]. Definition ret : A ! LLI A := fun a s ) val (a, s). Definition bind : LLI A ! (A ! LLI B) ! LLI B := fun m f s ) match m s with |val (a, s') ) f a s' |undef a s' ) undef a s' end.

  Definition writeVirtual (addr:page)(idx:index)(va:vaddr): LLI unit:= modify (fun s ) {| currentPartition := s.(currentPartition); memory := add paddr idx (VA va) s.(memory) beqPage beqIndex|}).

  {{fun s ) HI s^KI s ^VS s ^C s ^isChild child s îsPresent srcp s ^isAccessible srcp s Ĥ srcp s ^isEmpty dst s}} writePhyEntry MMUtableDst idx srcp {{fun _ s ) HI s ^KI s ^VS s ^C s}}.

  Direct memory access is represented by a lookup function named select. The value datatype sums up the types of values that can be stored by Pip.

	Record state: Type:={ currentPartition: page;
	memory: list(paddr * value)}.

Inductive value : Type:= |PE: Pentry ! value |VE: Ventry ! value |PP: page ! value |VA: vaddr ! value |I: index ! value.

  The second category is HAL_read. It is about reading values from physical memory and includes the following primitives: readVirtual, readPhysical, readVirEntry, readPhyEntry, readAccessible, readPresent, readPDflag, readIndex. For example, readVirtual allows reading a virtual address from a physical location.

	match lookup p idx s.(memory) beqPage beqIndex with
	|Some (VA v) ) ret v |Some _ ) undefined 3 |None ) undefined 2 end.
	Definition readVirtual (p: page) (idx: index): LLI vaddr:=
	perform s := get in

  HI s ^KI s ^VS s ^C s}}.

	{{fun s ) HI s ^KI s ^VS s ^C s ^isChild child s srcp s ^isAccessible srcp s notShared îsPresent srcp s ^isEmpty dst s}} writeVirtual shadow2TableDst dst src;;
	writeVirEntry shodow1TableSrc dst child;;
	writePhyEntry MMUtableDst dst srcp
	{{fun _ s )

Table 1 :

 1 The overview of verification 6 Conclusions and further work In this paper we have presented the formalization of the protokernel Pip which is written in Coq and automatically translated to C code. It supports a hierarchical partitioning model that only requires managing partition virtual spaces and context switching. The main security property

This work was funded by the European Celtic-Plus Project ODSI C2014/2-12.

PROOF-ORIENTED DESIGN