Approximation of a simple Navier-Stokes model by monotonic rearrangement

We consider a Navier-Stokes model for compressible fluids in one space dimension. We show that it can be approximated by a time-discrete scheme combining the discretization of a trivial stochastic differential equation and the application of a suitable monotonic rearrangement operator In addition, our result can be easily extended to a related Navier-Stokes-Poisson system.

where λ, ǫ are positive constants. We assume both the density field ρ(t, ·) > 0 and the velocity field v(t, x) to be 1-periodic in space, with unit mean for the density. We denote by X(t, a) = a + ξ(t, a) the "material" positions of the fluid parcels, so that ξ is oneperiodic in a ∈ R, ∂ t X(t, a) = v(t, X(t, a)) and ∂ a X(t, a)ρ(t, X(t, a)) = 1. Then, we show that X(t, a) = a + ξ(t, a) can be obtained as the unique limit, as h > 0 goes to zero and nh goes to t, of the time-discrete approximation X h,n (a) = a + ξ h,n (a) obtained from the recursion formulae X h,n (a) = [a + (1 + hλ)ξ h,n−1 (a) + hZ h,n−1 (a)) + √ 2ǫhN(a/h)]] ♯ Z h,n (a) = (1 − λh)Z h,n−1 (a) − λ 2 ξ h,n−1 (a)), (0. 2) where X h,0 , Z h,0 are suitably initialized, N is a 1-periodic function with average 0 and variance 1, and ♯ denotes the rearrangement operator in increasing order for functions depending on a ∈ R. Since the term N(a/h) can be interpreted as a deterministic approximation of a random variable W with zero mean and unit variance, we see that, in some sense, the resolution of the Navier-Stokes system has been reduced, through the crucial intervention of the monotonic rearrangement operator, to the discretization of a trivial stochastic differential equation, namely dξ = √ 2ǫdtW + (λξ + Z)dt, dZ = −(λξ + Z)λdt.
In addition, our result can be easily extended to a related Navier-Stokes-Poisson system.

Introduction
Optimal transport theory has been succesfully used to treat many different types of parabolic equations as gradient flows of various "entropy functionals" for various "transportation metrics", the canonical example being the regular scalar heat equation viewed by Jordan, Kinderlehrer and Otto [21] as the gradient flow of the Boltzmann entropy for the quadratic Monge-Kantorovich MK2 (frequently nammed Wasserstein) metric. There has been a huge literature on this subject in the last 20 years (study of the heat equation in a very general framework, porous-medium equations, thin-film flow equations, chemotaxis models, etc.., see, as a sample, [1,23,2,18]). Therefore, it is natural to try to treat the Navier-Stokes or the Euler equations of fluid mechanics in a similar way. Let us mention some attempts such as the introduction of a JKO scheme for the Euler equation by Gangbo and Westdickenberg [17] or the related treatment of the Navier-Stokes equation by Gigli and Mosconi [20]. However, we are probably still very far to be able to solve the initial value problem for the Euler or the Navier-Stokes equations by optimal transport tools. Nevertheless, there is more hope for fluid models in one space dimension. Indeed, in one dimension, optimal transport according to the MK2 metric is rather simple, as well known, since it can be entirely rephrased in terms of monotonic rearrangement and pseudo-inverse functions. Very simple models, such as the inviscid Burgers equation (and more generally multidimensional scalar conservation laws) or one-dimensional pressure-less Euler equations, have been successfully written, using pseudo-inverse functions and monotonic rearrangements, directly as subdifferential equations with nice contraction properties in all L p spaces [3,10,12]. In the present paper, we show that similar ideas can be extended to some Navier-Stokes models for compressible fluids in one space dimension, leaving widely open the the much more challenging case of 2D and 3D Navier-Stokes equations. More precisely, let us consider the (over-simplified) fluid mechanics model of Navier-Stokes type, without energy equation, in one space variable: x), respectively denote the density and the velocity of the fluid at time t > 0 and point x, on the real line, and p, µ denote the pressure and the viscosity of the fluid that are supposed to be just given functions of the density. (The genuine Navier-Stokes equations would include an additional equation for the temperature field, on which both p and µ would depend.) Our approach is restricted to the (rather unphysical) case: [As a matter of fact, if we neglect temperature effects, as we have done, it would be physically more consistent to assume the viscosity to be a constant and the pressure to be proportional to the density, but our method would not apply then.] The ǫ parameter will be interpreted subsequently as a level of noise in the formulation we will propose for the Navier-Stokes system. In particular, the "pressure-less viscosity-less" case p = µ = 0 will correspond to a zero level of noise. The λ parameter just scales the pressure with respect to the viscosity.
Then, we find (after standard calculations with a multiple use of the chain rule) or, equivalently, Our key observation is that this system looks like a very mild modification of X(t, a) = a + ξ(t, a), which is nothing but than the "material" version of the linear heat equation (written in "JKO style" [21]) through the same change of variables (1.5). This suggests the following discrete scheme, based on the simulation of the heat equation by random walk, to approximate the Navier-Sokes equation written in material coordinates (1.7). Given a uniform time step h > 0 and a uniform grid on R with mesh size 1/M. we approximate X(t, a), at each discrete time t = nh, n = 0, 1, 2, · · ·, by a non-decreasing sequence k → X n,k corresponding to each subinterval (k − 1)/M < a < k/M, for k ∈ Z. Because of spatial periodicity, it is convenient to introduce which is 1-periodic in a. Accordingly, we require X n,k+M = 1 + X n,k (which makes the scheme effective for computation) and introduce ξ n,k = X n,k − k/M. We first perform a predictor stepX n+1,k =ξ n+1,k + k/M ξ n+1,k = (1 + hλ)ξ n,k + hZ n,k + √ 2ǫhN n,k 11) where N n,k is a sequence, M-periodic in k, of M independent random numbers, with expectation 0 and unit variance. So far, this can be seen as a discretization of the rather trivial stochastic differential equation where W is a normalized white noise. However, we observe that there is no reason for k →X n+1,k to be increasing. So, we perform a corrector step just by sorting this sequence in increasing order, which is enough to define the next sequence k → X n+1,k . As discussed below, we even simplify the scheme by using a deterministic simulation of the white noise, for instance by substituting for the random number N n,k the value +1 if k is even, and −1 if k is odd, independently of n. This simplification is possible, because the corrector step generates a large amount of mixing. Let us point out that the sorting step is crucial to lift the rather uninteresting stochastic ordinary differential system (1.12) to the level of the more intriguing Navier-Stokes equations! It is now very easy to define the natural continuous version of this scheme with respect to the space variable a, while keeping t discrete. To do that we use the concept of rearrangement of functions in increasing order, which is just the continuous version of the sorting operator for sequences of real numbers. However, we have to be a little bit careful to handle the required space periodicity. Let us consider a given function Y ∈ I + L 2 (T), where I denotes the identity map on R and L 2 (T) is the space of all locally square integrable 1-periodic functions on R. Then, its rearrangement in increasing order is the unique function Y ♯ ∈ I + L 2 (T) such that (1.13) for all 1-periodic continuous function G, G ∈ C 0 (T). Then, the semi-discrete version of the scheme just reads (with obvious notation): X n+1 (a) =ξ n+1 (a) + a, ξ n+1 (a) = (1 + hλ)ξ n (a) + hZ n (a) + √ 2ǫhW Z n+1 (a) = (1 − λh)Z n (a) − hλ 2 ξ n (a), where W is a normalized white noise. As a matter of fact, we rather use a deterministic approximation of the white noise term, namely where N is a fixed 1-periodic bounded function with zero mean and unit variance and 1/r is a sufficiently large integer. This prevents us to rely on any stochastic analysis. Let us point out, before stating our convergence result, that, due to its very definition, the semi-discrete scheme immediately becomes a fully discrete scheme, under the following conditions. We first approximate the initial data (X 0 (a), Z 0 (a)) and the identity map a by piecewise constant functions in a, on a uniform grid, with constant values in each interval ](k−1)/M, k/M[, for k ∈ Z. Observe that the uniformity of the grid is crucial, because the rearrangement operator in increasing order preserves piecewise constant functions only on a uniform grid. We also crucially need the "noise function" N(a/r) to be itself piecewise constant on the same grid. Let us consider for instance the simplest noise function N(a) = −1 for 0 < a < 1/2 and N(a) = 1 for 1/2 < a < 1. Since L = 1/r is supposed to be an integer, we see that N(a/r) will be piecewise constant on the grid of mesh width 1/M as soon as we choose M to be a multiple of 2L. If we choose, with economy, Then the time-discrete scheme, without modification, coincides with the fully-discrete scheme discrete (1.11) where we substitute for the random number N n,k the value +1 if k is even, and −1 if k is odd. Let us also observe that, by construction, the fully discrete scheme can be trivially coded (with the help of a fast sorting algorithm, taking into account periodicity), and the resulting complexity, to reach the solution at a finite fixed time t, is of order O(h −1 r −1 log(1/r)) on a sequential machine. It is now time to state our main result.
Theorem 1.1. Let us consider the time-discrete scheme (1.14,1.15), where N is the 1periodic binary function with value −1 for 0 < a < 1/2 and +1 for 1/2 < a < 1, and 1/r an integer large enough so that r = O(h). Let us consider initial conditions X 0 = I + ξ 0 , Z 0 , where I is the identity map and ξ 0 and Z 0 are smooth 1-periodic functions, with inf X ′ 0 (a) > 0. Let us linearly interpolate in time the (X h,n , Z h,n ) and denote the result by (X h , Z h ). Then, as h → 0, the entire family converge to a solution (X, Z) of the material formulation of the Navier-Stokes system (1.7), under assumption (1.4).
Symbolically, to describe this limit, we suggest the following continuous formulation: where W stands for a normalized white noise and ♯ denotes the rearrangement operator in increasing order.
Our method of proof will involve three main steps: i) The Navier-Stokes system, written in material coordinates, is first interpreted as a modified gradient flow with delay terms for a suitable convex functional, in a slightly more general situation. As a byproduct, we establish the global L 2 well-posedness of the system, in a general setting of non smooth initial conditions (including vacuum, i.e. zones where the density field vanishes). ii) In the special case (1.4), we show that the previous system can be further reduced to a non-autonomous viscous scalar conservation laws, for which existence of global smooth solutions can be routinely established. iii) Finally, we observe that our time-discrete scheme is an adaptation of the "transportcollapse-kinetic" method [4,5,19,6,8] to the discretization of a suitable coupled scalar parabolic equation (3.37) and prove the desired convergence theorem.
Let us conclude this introduction by mentioning that our approach can be adapted to some interesting cases when a self-consistent forcing term is added to the NS equations. This is in particular true for the Navier-Stokes-Poisson (NSP) system where β is a constant, φ = φ(t, x) can be interpreted ad a gravitational potential if β > 0, or as an electrostatic potential if β < 0, with a "neutralizing" background of unit mean. In material coordinates, this system reads For this Navier-Stokes-Poisson system, under assumption (1.4), we just have to modify the noisy differential system in a straightforward way where F is given by (1.19). So we may immediately include the Navier-Stokes-Poisson case, just by adding a smooth forcing term, F (y) with bounded derivative in y and F (0) = 0, to the Navier-Stokes equations, as in (1.18), without necessarily assuming that F has form (1.19). This is what we will do through the paper.

Reduction to a gradient flow with delay
In this section, we consider the Navier-Stokes model (written in material coordinates) (1.18), in the case when the viscosity coefficient µ and the pressure p are linked to each other through a given smooth strictly convex function τ > 0 → θ(τ ), in the following way: for some constants ǫ, λ > 0. We recover the case (1.4) discussed in the introduction for the special choice: Another interesting example is θ(τ ) = τ (log τ − 1), for which p = −λǫ log τ , and µ = ǫ (i.e. constant viscosity). Under assumption (2.21), the NS equation written in material coordinates, with forcing F , (1.18) becomes (2.23) ∂ 2 tt X(t, a) = λǫ ∂ a (θ ′ (∂ a X(t, a))) + ǫ ∂ a (θ"(∂ a X(t, a))∂ 2 ta X(t, a)) + F (X(t, a) − a) which can be immediately blown up as a first order system: The resulting system is a mild modification of a gradient flow for a convex functional. More precisely, let us introduced the Hilbert space H = L 2 (T) and the closed convex cone K of all nondecreasing functions in I + H. otherwise.
[Notice that whenever X = I + ξ ∈ K, its derivative X ′ (a) = 1 + ξ ′ (a) can be seen as a locally bounded nonnegative Borel measure on T. Then Θ[ξ] can be more precisely defined as "a convex function of measures" by Legendre duality, as done, for example, in [15]: where f describes the set of continuous function on T and θ * is the Legendre-Fenchel transform of θ:  )).] So, the first order system (2.25) can be viewed, in the Hilbert space H, as a sub-differential inclusion coupled to a linear ODE in Z, namely: (2.28) By integrating out Z from the second equation, we may view the system as a gradient flow with delay terms (which are harmless since we have assumed F (y) to have bounded derivative in y and F (0) = 0, as in the Navier-Stokes-Poisson case (1.19)): Such formulations are very useful to define global "generalized" solutions for a large set of initial conditions, with well-posedness in L 2 , just by using, in a routine way, the general theory of maximal monotone operators [13]. In addition, in this framework, the limit ǫ → 0 is trouble-free. Indeed, in the sense of maximal monotone operator theory, the limit of ǫΘ[ξ] is just the indicator function 1 K [X] of the cone K, with value 0 whenever X = I + ξ belongs to K and +∞ otherwise. Thus, we may conclude this section by asserting: Theorem 2.1. Let (ξ 0 , Z 0 ) given in the Hilbert space H = L 2 (T), with X 0 = I + ξ 0 in the closed convex cone K of all non-decreasing functions in I + H. Then, there is a unique global solution t → (ξ(t), Z(t)) ∈ H ×H to (2.28), depending continuously in H on both t and the initial condition (ξ(0) = ξ 0 , Z(0) = Z 0 ). This solution can be viewed as a "generalized" solution to the Navier-Stokes system (1.18), written in material coordinates, with pressure and viscosity laws given by (2.21), and, in particular for (1.4). In addition, as ǫ → 0, these solutions converge to those of: is the indicator function of the cone K, with value 0 whenever X = I + ξ belongs to K and +∞ otherwise.
Remark 1. Formulations of type (2.30) have already been introduced in the absence of viscosity for the pressure-less Euler equations and some of their variants. Let us quote the model of order-preserving vibrating strings, the Chaplygin gas, the pressureless Euler-Poisson system and also multi-dimensional scalar conservation laws etc ... ([9], [10], [22], [11], for example). In all of these earlier works, λ has never been taken into account. It is unclear to us, whether or not λ plays (or should play) a role in these zero viscosity limit equations. We leave that for further investigations.
Remark 2. Subdifferential formulation (2.28), or (2.29), leads to stability of solutions with respect to their initial condition not only in L 2 but also in all L q spaces. Indeed, let us consider, on a fixed finite time interval [0, T ], two solutions (ξ, Z) and (ξ,Z). Using formulation (2.29), we get where c depends only on q, T, λ and the Lipschitz constant of F , but not on ǫ. This easily leads to the L q stability estimate, uniformly on ǫ, with (another) c depending only on q, T, λ, F , but not on ǫ.

Reduction to a scalar parabolic equation with mild coupling
The sub-differential formulation (2.28) of the Navier-Stokes equations in material coordinates (1.18), that we have introduced in the previous section, is very powerful since it provides L 2 (and, more generally L q ) well-posedness for very general data (including vacuum). However, it is not clear that the corresponding solutions are smooth, for appropriate initial conditions. At least in the special case (2.22) where θ(τ ) = − log τ , this can be done by reducing (1.18) to the scalar parabolic equation where Z is coupled to u in a suitable way. For this purpose, let us introduce the inverse function u of X = I + ξ (3.33) u(t, X(t, a)) = a.
At the moment, let us assume X(t, a) = a + ξ(t, a) to be smooth with ∂ a X > 0 (which will be justified a posteriori in the case (1.4)). We have: u(t, X(t, a)) = a, (∂ x u)(t, X(t, a))∂ a X(t, a) = 1,

Using (3.34), we get
Thus, we get for (u, Z) the following system: ∂ t Z(t, a) + λZ(t, a) = F (X − a) − λ 2 (X(t, a) − a), u(t, X(t, a)) = a.. By integrating out Z, we get, more explicitly, Next, under assumption (2.21,2.22), i.e. θ(τ ) = − log τ , we find  Observe that the coupling of Z with u, through X is very mild, so that the theory for this system differs very little from the standard theory of viscous conservation laws, such as (3.38) ∂ t u + ∂ x (f (u)) = ǫ∂ 2 xx u. So, without entering in details, we claim that for u 0 given in the class C of all orientationpreserving diffeomorphisms of R/Z and Z 0 smooth 1-peridoic function, then there is a unique smooth solution (u(t, ·) ∈ C, , Z(t, ·)) to the coupled scalar parabolic equation (3.37). L 1 stability with respect to initial conditions. In addition, thanks to the the subdifferential formulation (2.28), we get an L 1 stability result for every pair (u, Z) and (ũ,Z) of solutions to the coupled scalar parabolic equation (3.37). on every fixed finite time interval [0, T ]: with c depending only on T, λ, F , but not on ǫ. This immediately follows from the L q stability estimate (2.31) for the subdifferential formulation (2.28), in the special case q = 1.

Approximation by a "transport-collapse-kinetic" method with noise
The transport-collapse-kinetic method is a time-discrete scheme for hyperbolic and viscous conservation laws introduced and studied in [4,19,6,8]. It was suggested to the author by an earlier algorithm proposed by A. Chorin [16] for reaction-diffusion equations. This scheme turns out to be very well suited for the formulations (3.37) and (2.28) we have obtained in the previous sections for the Navier-Stokes model (1.18). Let us now describe the time-discrete scheme: Given a uniform time step h > 0, we denote by (X h,n (a) = a + ξ h,n (a), Z h,n (a)) the approximate solution at time t = nh, for n = 0, 1, 2, · · ·, defined in two steps as follows. We first define Z h,n+1 andX h,n+1 = a +ξ h,n+1 (a) bŷ where N is a 1-periodic function with average 0 and variance 1, and 1/r is a large integer that may depend on h > 0. [Let us recall that this first step has been designed to approximate the noisy ordinary differential system (1.12), namely: where W is a normalized white noise, with a deterministic approximation to W .] Next, we crucially introduce a rearrangement step. Namely, we rearrangeX h,n+1 (a) = a +ξ h,n+1 (a) in increasing order with respect to a, and obtain X h,n+1 (a) = a + ξ h,n+1 (a). In other words, we set (4.42) X h,n+1 =X ♯ h,n+1 , where we denote by Y → Y ♯ the rearrangement in increasing order, defined by (1.13), for functions differing from the identity map I by a L q 1-periodic function (i.e. belonging to I + L q (T)). Of course, the resulting scheme contains the scheme (1.14) discussed in the introduction, up to the addition of the (harmless) term F , that we suppose, as previously mentioned, smooth with bounded derivative. Notice that the scheme is almost translation invariant in the variable a. Indeed, the rearrangement step is clearly translation invariant, and, in the predictor step, the only non-autonomous term is the "noise" term which is proportional to N(a/r). Thus the scheme is invariant by any translation a → a + k/r, with k ∈ Z. This observation will be very important for our analysis. We also observe that Z h is entirely "slaved" by ξ h and depends linearly on it through (4.40), so we are only concerned about monitoring ξ h . (Of course we use that F is smooth and Lipschitz with F (0) = 0.) For the subsequent analysis, it is also simpler to restrict ourself to the case when N is the binary function with value −1 for 0 < a < 1/2 and +1 for 1/2 < a < 1. Indeed, in that case, remarkably enough, (ξ 0 , Z 0 ) = (0, 0), is a fixed point of the scheme, as F = 0. [This can be easily checked graphically: we see the exact compensation between the addition of the "noise" term √ 2ǫhN(a/r) and the rearrangement step.] Let us now prove the following convergence result (from which our main result Theorem 1.1 immediately follows as the particular case when the force term F is absent): Theorem 4.1. Let us consider the time-discrete (4.40,4.42), where N is the 1-periodic binary function with value −1 for 0 < a < 1/2 and +1 for 1/2 < a < 1 and 1/r an integer large enough so that r = O(h). Let us consider initial conditions X 0 = I + ξ 0 , Z 0 , where I is the identity map and ξ 0 and Z 0 are smooth 1-periodic functions, with inf X ′ 0 (a) > 0. Let us linearly interpolate in time the output of the semi-discrete scheme (X h,n , Z h,n ) and denote the result by (X h , Z h ). Then, as h → 0, the entire family converge to (X, Z) solution of the material formulation of the Navier-Stokes system (1.18), under assumption (1.4).

Proof.
Sup-norm and Lipschitz estimates. The rearrangement operator is a non-expansive map on all spaces I + L q (T), in particular for q = ∞. So we can compare in sup-norm any fixed solution (ξ h,n , Z h,n )(a) of the time-discrete scheme to: i) (ξ h,n , Z h,n )(a + kr), its space translation by kr, for any integer k ∈ Z, which is also solution of the scheme (due to translation invariance, as already mentioned); ii) the fixed point solution already (ξ, Z) = (0, 0), obtained in the case F = 0, Let us perform the first comparison. We easily get: ||(ξ h,n , Z h,n ) − (ξ h,n , Z h,n )(· + kr)|| ∞ ≤ (1 + hc) n ||(ξ 0 , Z 0 ) − (ξ 0 , Z 0 )(· + kr)|| ∞ ≤ |k|rLip(ξ 0 , Z 0 ) exp(nhc), (4.43) where Lip denotes Lipschitz constants and c is a constant depending only on λ and the Lipschitz constant of F . [Indeed, we first get, for the predictor step, where ǫ and the "noise" term do not play any role, since the translation is an integer multiple of r. Next, using the contraction property of the rearrangement operator, we get which is enough to get (4.43.] As a matter of fact, we can improve estimate (4.43) using the fact that X h,n (a) = a+ξ h,n (a) is increasing in a. Indeed, we get, for all ω ∈ [kr −r, kr], ξ h,n (a + ω) − ξ h,n (a) = X h,n (a + ω) − X h,n (a) ≤ X h,n (a + kr) − X h,n (a). Thus, Since Z h,n is slaved by X n,h , we have obtained ||(ξ h,n , Z h,n ) − (ξ h,n , Z h,n )(· + ω)|| ∞ ≤ (|ω| + r)Lip(ξ 0 , Z 0 ) exp(nhc), (4.44) for all ω ∈ R, where c is a constant depending only on λ and the Lipschitz constant of F and not on ǫ. This is our first key estimate for the scheme. Next, let us compare the solution to the fixed point of the scheme (ξ = 0, Z = 0), already mentioned and obtained for F = 0. For this fixed point, the predictor step provides the values ( √ 2ǫhN(a/r), 0) Thus, we get (because the rearrangement operator is non-expansive) (by definition of the predictor step, using that F is Lipschitz and F (0) = 0), where c is a constant depending only on λ and the Lipschitz constant of F . So, we have obtained our second important bound where c depends only on λ and the Lipschitz constant of F , and not on ǫ.
Consistency estimates. From now on, we fix a time interval [0, T ] and denote by c any constant depending only on the data T, X 0 = I + ξ 0 , Z 0 , N, F, ǫ, λ. Let G be a smooth 1-periodic test function. Let us fix a time step n. Since G is 1-periodic, because of the rearrangement step, we have By Taylor expansion in h, we get for the right-hand side of (4.47). +G"(X h,n (a))ǫh}da + O(h 3/2 ), There are two terms to deal with: Since 1/r is assumed to be an integer, by 1-periodicity of G, X h,n − I and M, we first find I 3 = 0. Then, using that X h,n (a) is non-decreasing in a, we get where c denotes a generic constant depending on the data and the test function G. Thanks to the sup norm estimate (4.45), we finally get |I 1 | ≤ cr √ h. Thus I 1 can be absorbed in the error term O(h 3/2 ), just by assuming r = O(h). Let us now consider I 2 = I 5 + I 6 + I 7 , where  g(X h,n (a))(λξ h,n (a) + Z h,n (a))da, ∀g ∈ C 0 (T). (4.49) Notice that both are uniformly bounded in (h, n) as Borel measures, since X h,n and X h,n are uniformly bounded, as nh ≤ T , for all finite T . We obtain, as a substitute for (4.48), in the sense of distributions on T.
Compactness. We first perform a linear interpolation in time of the measures (ρ h,n , Q h,n ) which generates a continuous piecewise linear function of time t → (ρ h (t), Q h (t)) valued in the dual space C 0 (T) ′ ( i.e. the space of all bounded Borel measures on T, by Riesz identification theorem). For each fixed T > 0, this family is bounded in Lip([0, T ], C 2 (T) ′ ), because of (4.50), and therefore compact in C 0 ([0, T ], C 0 (T) ′ ). Let us consider a convergent subsequence, for a sequence of time steps h m going to zero, and denote the limit (ρ, Q).
Since Z h is "slaved" by ξ h and depends "linearly" on it through (4.40), we are just concerned by the compactness of ξ h . Because of the sup norm and Lipschitz estimates (4.45,4.44), we already know that (ξ h (t)) is valued in a fixed compact set of L 2 (T). Next, we use in a crucial way the following lemma which comes from rearrangement theory (or, alternately, from "optimal transport theory") [7,23]:  Using this proposition, we deduce from (4.48) that (ξ h ) is uniformly equi-continuous from [0, T ] to H, for all 0 < T < +∞. We conclude, by the Arzela-Ascoli theorem that the family (ξ h , Z h ) is relatively compact in C 0 (R + , H), remembering that Z h is linearly slaved by ξ h .
Consistency. We now have enough compactness to pass to the limit: up to the extraction of a subsequence of h → 0, we have at least a limit (ξ, Z) in C 0 (R, H) for the family (ξ h , Z h ). Because of the sup-norm and Lipschitz estimates (4.45,4.44), we also know that (ξ(t, a), Z(t, a)) are bounded and Lipschitz continuous in a, uniformly in t for all bounded time interval t ∈ [0, T ]. Concerning Z h , passing to the limit in the second equation of (4.40) is straightforward (using that F is smooth and Lipschitz). We get: Next, from (4.48), we deduce that for all smooth 1-periodic function G d dt {G ′ (a + ξ(t, a))(λξ(t, a) + Z(t, a)) + ǫG"(a + ξ(t, a))}da. (4.53) Uniqueness of the limit. Let us consider any limit (ξ, Z) of the scheme, on a finite interval of time [0, T ]. We know that it belongs to L ∞ ([0, T ], Lip(T)). Thus, for each t, X(t, a) = a + ξ(t, a) has an inverse function u(t, x) which is increasing and satisfies u(t, x + 1) = 1 + u(t, x). In addition, there is a constant α > 0 such that ∂ x u(t, x) ≥ α. (However, at this stage, u(t, x) may have jumps in x). Its derivative ρ(t, x) = ∂ x u(t, x) ≥ α is, for each t, a probability measure on T, bounded away from zero. We may also introduce, for each time t, the real valued bounded (uniformly in t) measure Q(t, x) on T, defined by g(X(t, a))(λξ(t, a) + Z(t, a))da, ∀g ∈ C 0 (T), (4.54) which is absolutely continuous with respect to ρ(t, x) with a bounded density (since X and Z are bounded). Thus, (ρ, Q) belongs to L ∞ ([0, T ], C 0 (T) ′ ), where C 0 (T) ′ is the space of all bounded (Borel) measures on T. The consistency relation (4.53), written in terms of ρ and Q, exactly means ∂ t ρ + ∂ x Q = ǫ∂ xx ρ, (4.55) in the sense of distribution. Their Fourier coefficients F ρ(t, k), F Q(t, k) are uniformly bounded in t ∈ [0, T ] and k ∈ Z, and satisfy (∂ t + 4π 2 ǫk 2 )F ρ(t, k) = 2iπkF Q(t, k). This implies that ρ belongs to L ∞ ([0, T ], H s (T)) (where H s (T) denotes the standard Sobolev space of 1-periodic function with derivatives in L 2 up to order s, for s ∈ R) for all s < 1/2, and, therefore, by Sobolev embedding theorem to all L ∞ ([0, T ], L q (T)) for 1 ≤ q < +∞. Since the density of Q with respect to ρ is bounded, we deduce that Q ∈ L ∞ ([0, T ], L q (T)) for all 1 ≤ q < +∞. Since ρ = ∂ x u, we also deduce u ∈ L ∞ ([0, T ], H s (T)) for all s < 3/2, and, thus u ∈ L ∞ ([0, T ], C 0,σ (T)) for all σ < 1. Thus, from the definition of Q (4.54), we can now write ρ(t, x) = ∂ x u(t, x), Q(t, x) = ((x − u(t, x))λ + Z(t, u(t, x))∂ x u(t, x) (4.57) and deduce that Q, just like ρ, belongs to L ∞ ([0, T ], H s (T)) for all s < 1/2. [For that, instead of using Fourier coefficients, we estimate translations of Q in L 2 . Using that u is bounded and Z Lipschitz, we get I = Since ρ belongs to L ∞ ([0, T ], H s (T)) for all s < 1/2, we immediately get that I 1 ≤ c s |ω| 2s , for all s < 1/2. Next, we use that u belongs to L ∞ ([0, T ], H s (T)) for all s < 3/2. This implies by Sobolev embedding theorem, that u and Z(u) certainly belong to L ∞ ([0, T ], W 1,4 (T)), while ρ belongs to L ∞ ([0, T ], L 4 (T)). This implies that I 2 ≤ c|ω| 2 . So we see that I ≤ c s |ω| 2s for all s < 1/2, which, indeed, means that Q belongs to L ∞ ([0, T ], H s (T)) for all s < 1/2.] So we can now bootstrap the regularity of ρ and Q, using (4.56). We first get that ρ, belongs to L ∞ ([0, T ], H s (T)) for all s < 3/2 and u to the this space for all s < 5/2. By Sobolev embedding theorem, this means that u also belong to L ∞ ([0, T ], C 1,σ (T)) for all σ < 1. We already know that ∂ x u(t, x) ≥ α > 0. Since X(t, a) is the inverse function of u, it follows that ξ(t, a) = X(t, a) − a also belongs to L ∞ ([0, T ], C 1,σ (T)) for all σ < 1. This is also true for Z, which is slaved by ξ through (4.52). At this point, the boostrap argument becomes pure routine and we conclude that u, ξ, Z are smooth, just as the initial conditions. From (4.55,4.57,4.52), we see they are the unique solution in classical sense of the system u(t, X(t, a)) = a, X(t, a) = a + ξ(t, a), ∂ t u(t, x) + [Z(t, u(t, x)) + λ(x − u(t, x))]∂ x u(t, x) = ǫ∂ 2 xx u, ∂ t Z(t, a) + λZ(t, a) = F (ξ(t, a)) − λ 2 ξ(t, a). This concludes the proof of our theorem since this system is just the formulation (3.37) of the Navier-Stokes system (1.18).