Uniform approximation of 2D Navier-Stokes equations with vorticity creation by stochastic interacting particle systems

We consider a stochastic interacting particle system in a bounded domain with reflecting boundary, including creation of new particles on the boundary prescribed by a given source term. We show that such particle system approximates 2D Navier–Stokes equations in vorticity form and impermeable boundary, the creation of particles modeling vorticity creation at the boundary. Kernel smoothing, more specifically smoothing by means of the Neumann heat semigroup on the space domain, allows to establish uniform convergence of regularized empirical measures to (weak solutions of) Navier–Stokes equations.

where the choice of Dirichlet conditions for the Green function (−∆ Dir ) −1 makes it so that the vorticity dynamics encodes the impermeability condition u| ∂D • n = 0 at the boundary.The tangential part u ⊥ | ∂D • n = 0 of the no-slip boundary condition, on the other hand, is harder to express in terms of ω only, giving rise to a non-local condition on vorticity (see [3,Section 1.1] for an explicit expression).Physically speaking, the condition u ⊥ | ∂D • n = 0 is expected to force production of vorticity, especially in the proximity of the boundary, so that the tangential velocity at ∂D induced by the bulk of the fluid be compensated [5].This fact becomes an issue in the study of approximation methods for Navier-Stokes equations in vorticity form, which is a natural approach in the 2-dimensional setting.
Previous approaches to approximating vorticity creation at the boundary for 2-dimensional Navier-Stokes equations [8,5,12,11,29,36] have considered (unphysical) Neumann boundary conditions for the vorticity dynamics, (NNS) where g ∈ L 2 ([0, T ], L 2 (∂D)) and ν > 0. The latter system acts as a natural proxy for the (physical) no-slip condition if one seeks to mathematically describe vorticity production in the boundary layer by means of the source term of Neumann boundary condition or with an explicit addition of vorticity sheets, recovering the correct PDE dynamics (1.1) in a macroscopic or iterative limit.The derivation of (1.1) on a generic domain D by means of an approximation method explicitly describing the effect of the boundary as a singular source of vorticity remains a mostly open problem.The present contribution aims to lay the groundwork for the construction of an approximation method for (1.1) based on interacting particle systems uniformly converging to the PDE system, that is, seeking to provide at the same time a good microscopic description of vorticity creation, and a good approximation scheme.In order to do so, we generalize to the bounded domain setting the analytic and probabilistic tools introduced in [17,20] for the uniform approximation of PDEs by (stochastic) interacting particle systems in full-space.
We consider a system of N particles evolving according to the stochastic differential equations with reflecting boundary (1.2) each particle starting its evolution at time t i ∈ [0, T ] at x i (t i ), either located in the interior of D if t i = 0 or on ∂D if t i > 0. The system is closely related to the one considered in [29], the main difference being the generation mechanism.
Let us briefly and informally describe the various parts of the dynamics (precise definitions are all deferred to section 3).Generation of particles, that is initialization of each of the SDEs (1.2), happens either in the interior D at time 0 (so to model the initial datum of the limiting PDE) or at the boundary ∂D at later times (so to model the source term g of the Neumann b.c.).Concretely, we will consider a grid spanning ({0} × D)∪((0, T ]× ∂D) indexed by i, and generate the particle x i at grid point i with an intensity ξ i determined by a local average of the boundary datum around x i .The mesh of the grid (both in time and space) will be 1/n, thus n ∈ N, n → ∞ is a convenient parameter to rule the macroscopic limit: for instance the global number of particles will have order N = N (n) ≃ n 2 .Weights (or intensities) ω i of particles will correspond to local averages of the initial datum ω 0 or the boundary datum g around the starting point of the trajectory of x i .
Interaction of particles is given by a regularized version of the singular kernel K, obtained by applying the Neumann heat semigroup, (1.3) K n (x, y) = D K(x, z)p ε (z, y)dz, ε = ε(n), x, y ∈ D, with p ε the Neumann heat kernel on D and ε → 0 with a slow enough rate in terms of n.Kernel smoothing by means of the heat semigroup P t , compared for instance to smoothing by convolution with bump functions, has many advantages in our setting: it is well suited for regularizing functions keeping their support in D, P t preserves the L 1 norm of non-negative functions, and precise estimates on p t and its derivatives are available.The interaction term of (1.2) also includes a cutoff, (1.4) (which is Lipschitz continuous, and its Lipschitz constant is uniformly bounded in M > 0).In principle, the presence of the cutoff makes it so that the macroscopic behavior of the particle system is described by a different PDE, that is coinciding with (NNS) only when the velocity field in the latter is bounded.As revealed by suitable a priori estimates, this is in fact the case, so the cutoff F is just a technical device that will naturally disappear in the macroscopic limit taking M > 0 large enough.Regularization of interactions is important to ease uniform convergence of empirical measures and well-posedness of the SDE system.Brownian noise acting independently on each particle models viscosity, and boundary terms −dk i are defined so that ∂D acts on particles as a reflecting boundary, in the standard setting of [42,31].In terms of boundary conditions for the macroscopic limit, this would produce the null Neumann condition n • ∇ω = 0, so an additional effect must be included to model the source g.
Our main result concerns the uniform convergence (in space) of kernel-smoothed empirical measures of the particle system under consideration, (1.5) S n (t) = i∈A n t ω i δ xi(t) , ω n (t, x) = P ε S n (t), where, P ε stands for the action of Neumann heat semigroup on measures, and A n t collects the indices of particles generated up to time t.The generation of new particles at the boundary produces discontinuities in time at the level of ω n (t, x), therefore our limit theorem is going to be set in the space of C( D)-valued càdlàg functions, D([0, T ], C( D)), endowed with the usual Skorohod topology.Informally, our main result (Theorem 3.4) is the following: the regularized empirical measures ω n (t, x) converge as n → ∞, almost surely in the topology of D([0, T ], C( D)), to the unique weak solution of (NNS) with Lipschitz continuous initial datum ω 0 and bounded Neumann boundary datum g.As detailed below, we actually obtain a stronger convergence in Sobolev norms.We will in fact assume that g ≥ 0 in the exposition of our argument, and remove this further assumption in the later section 6, see Theorem 6.1.
The proof or Theorem 3.4 consists in a compactness argument, therefore the technical core of this paper resides in uniform estimates on the particle system.The main additional difficulty compared to the works [17,20] inspiring our technique is of course the presence of the boundary, which complicates many analytic operations.For instance, the absence of translation invariance in the models under consideration prevents derivatives and heat semigroups from commuting, so one has to resort instead to suitable gradient estimates on the heat kernel in order to exploit the regularizing properties of heat semigroup.Uniform convergence at time 0 and at the boundary also requires a careful control due to particle generation, but we are able to deal with rather natural assumptions on ω 0 , g (as compared with the obtained convergence), and to avoid further technical hypothesis on initial data such as [17,Assumption 1.1.3](a restricting requirement on the relation between the scale of kernel smoothing and sample size).
Beside the improved notion of convergence to the PDE we are able to obtain, the main point of our work is the approach to particle approximation for (NNS): our approach is almost completely functional analytic and based on semigroup theory, in contrast with results in the perhaps more classical framework of propagation of chaos such as [29,21].The motivation for our approach is that our technique circumvents the need for fine controls on particle dynamics.We regard this as a step forward in the scope of attempting to build a similar particle approximation of the physical no-slip Navier-Stokes equations (1.1), since in that case particle creation has to be linked with the velocity at the boundary induced by the bulk of the fluid, making single particle dynamics very hard to treat (as already commented in [29]).
The particle system under analysis can be thought of as a (regularized) stochastic version of the well-known point vortex system.However, there are some relevant differences concerning in particular the effect of the boundary on the vortex dynamics, see Remark 3.3 below.Let us only briefly discuss here some qualitative aspects for the sake of a comparison with the related literature.Due to the regularization of interaction at the microscopic level our approximants share some features with the systems of moderately interacting particles (see [16,17] for a comparison), but our focus is on the approximation procedure rather than on the scaling of ε = ε(n).If we were to omit the regularization of the interaction kernel and the cutoff F (and neglecting for a moment the complications of boundary effects), the scaling of weights ω i would make our particle system (1.2) a mean-field rescaled version of the stochastic point vortex model.While the mean field limit of deterministic point vortices system towards 2d Euler's equation in the well-posedness class is a recent result [39], our (stochastic) result (towards 2d Navier-Stokes) is of course closer to the ideas of [33,35].Concerning the relation with propagation of chaos theory, besides the above references closely related to our arguments, we refer for completeness to the work of Jabin-Wang [28] on stochastic particle systems with general interactions (on full space), including stochastic vortex dynamics, for which we also mention [10] discussing more general noise.The study of propagation of chaos in stochastic particle systems in bounded domain dates back at least to [43], we refer to [9] for a discussion and recent results.

Functional Analytic Setup
In what follows, dx denotes the Lebesgue measure restricted to the domain D, and dσ(x) the (1-dimensional) volume measure on the smooth boundary ∂D.In order to lighten notation, we will simply write L p = L p (D, dx), avoiding ambiguity by explicitly indicating the underlying space for any different L p (X).
We denote by B(X), B b (X) respectivly the spaces real-valued Borel and bounded Borel functions on a topological space X.
Landau O's and o's have their usual meaning, subscripts indicating eventual dependence on parameters.The symbol C will denote a positive constant, possibly differing in any of its occurrence even in the same formula, depending only on eventual subscripts.
The following statement summarizes standard results, we refer to monographs [27,44,45] for a complete discussion.Proposition 2.1.It holds: • (Sobolev embeddings) if 1 < p ≤ q < ∞, and α ≥ β ≥ 0, embeddings are compact if the latter inequality is strict; then the trace of ∇u on ∂D (defined as usual as the extension of the restriction to ∂D for smooth functions) is of class L p (∂D).Therefore, if n : ∂D → R 2 is the inward pointing unit normal vector field of ∂D, the following definition of Sobolev space with Neumann boundary conditions is meaningful: We refer to [45,Sec. 4.3.3]for a more general discussion of boundary conditions for Sobolev spaces, and to [45, Sec.2.9.3] for further regularity of the mentioned trace operator in terms of Besov spaces on ∂D.
Consider the heat equation with Neumann boundary conditions (2.1) x ∈ D, with f a bounded measurable function on D. It is well-known that the PDE problem is well-posed, we denote its solution by P t f (x) = u t (x), t ≥ 0, x ∈ D; moreover P t f ∈ C ∞ (D) for all t > 0. Operators P t , t > 0 form a Markov semigroup (in the sense of [2]), associated to the reflected Brownian motion in D (see subsection 3.2 below).We now recall classical facts concerning the action of the semigroup P t on L p (D) and Sobolev spaces: we refer to [38,Section 7.3] for a detailed discussion.For all p ∈ [1, ∞), P t extends to an analytic semigroup of contractions P t = e −tAp : L p → L p with infinitesimal generator , coinciding with the Laplace operator ν∆ on C ∞ c (D).The spectrum of A p is included in [0, +∞).A first consequence of the fact that P t is an analytic semigroup is ultracontractivity on L p spaces: we refer to [1, Section 7.3.2]for the following result.
Proposition 2.2.If 1 ≤ p < q ≤ ∞, then it holds (2.2) P t L p →L q ≤ C T,ν,p,q t 1/p−1/q , t ∈ [0, T ].Another consequence of analiticity of P t is that we can introduce fractional powers of I + A p as operators on L p by means of where the integral converges in the uniform operator topology and defines injective operators (cf.[38,Section 2.6]).We set (I + A p ) 0 = I, (I + A p ) α = ((I + A p ) α ) −1 for α > 0, and denote by H α,p = D((I + A p ) α/2 ) the domain of fractional powers for all α ∈ R. Proposition 2.3.Let p ∈ (1, ∞).It holds: • (group property) for α, β ∈ R and u ∈ H 2γ,p , γ = α ∨ β ∨ (α + β), • (identifications of the domain 1 ) H 2α,p can be identified with the interpolation space [L p , D(A p )] α for all α ∈ [0, 1]; moreover it holds for all t > 0 and operators P t and (I + A p ) α commute on H 2α,p ; moreover The action of P t on f ∈ L 1 is given by where the heat kernel p t (x, y) ≥ 0, t > 0, x, y ∈ D, is a smooth function, p t ∈ C ∞ ( D2 ), t > 0, satisfying (2.1) with the initial condition replaced by u t (x) → δ y as t → 0. The Neumann heat kernel can be controlled with the free heat kernel (on the whole plane) up to a multiplicative constant at the exponent, that is: As a consequence, for all t > 0 and p ∈ (1, ∞), As reported in [46,Theorem 3.2], the latter implies the following: 1 These identifications and Sobolev embeddings on D are closely related to the ultracontractivity of Pt, we refer again to [1] for a discussion.
We will also employ a gradient estimate for the heat semigroup acting on distribution spaces.
Lemma 2.6.For all t > 0, α > 0, p ∈ (1, ∞), Proof.Observe first that, by definition of H −α,p and density of H α,p in L p , it holds for all s > 0 the last step using (2.3).For 0 < s < t we can combine (2.7) and the estimate of above in order to obtain (2.9) from which the thesis follows minimizing the right-hand side with respect to the parameter s.
The following statement collects technical passages involving heat semigroups that will appear often in our computations, since they concern (duality) relations between operator norms of semigroups and uniform estimates of their kernels over the space domain D.
Proof.The first two claims follow from (2.4).For y ∈ D consider a sequence y n ∈ D such that y n → y; as for (2.10), the following chain of equalities holds due to the regularity of f and p ε : This implies (2.10) thanks to the fact that (I + A p ′ ) α/2 and the heat semigroup commute.The exchange between limit and integral is allowed thanks to (2.4).The proof of (2.11) follows from a similar argument.Indeed, Equation 2.11 follows from the regularity of P ε (I + A p ′ ) α/2 P t and the fact that (I + A p ′ ) α/2 and the heat semigroup commute.After these preliminaries, (2.12) and (2.13) are easy to prove by duality.Indeed, for each y ∈ D, Considering the supremum of both sides for y ∈ D (2.12) follows immediately due to the fact that (I + A p ′ ) α/2 P t+ε ∈ C( D). (2.13) is similar.Indeed Equation 2.13 then follows considering the supremum of both sides for y ∈ D.
We conclude the paragraph recalling the regularizing property of Biot-Savart kernel, following from the one of Green's function for Dirichlet boundaries, for which we refer to [27,Section 3].
Lemma 2.8.The linear operator K Moreover, for all the above extensions K[ω] is divergence-less (in the sense of distributions), and its normal trace on ∂D (when defined) vanishes.

2.2.
Càdlàg Functions and Aldous' Criterion.For T > 0, we denote by D([0, T ], S) the space of càdlàg functions on [0, T ] taking values in a complete metric space (S, d).We will always endow D([0, T ], S) with the Skorohod metric: we refer to [37, Section 2.1] for a definition of the latter, which we will not be using directly since we can rely on the following tightness criterion, [37,Theorem 3.2].Proposition 2.9 (Aldous' Criterion).Consider a sequence of filtered probability spaces (Ω n , F n , P n ) n∈N on each of which it is defined a càdlàg adapted process (X n t ) t∈[0,T ] taking values in a complete separable metric space (S, d).The laws of processes (X n t ) t∈[0,T ] are tight on D([0, T ], S) if the following two conditions are satisfied: (1) for every t in a dense subset of [0, T ] the laws of X n t are tight on S; (2) for all ε, δ > 0 there exists r 0 > 0 and n 0 ∈ N such that, for any sequence (τ n ) n∈N with τ n a F n -stopping time, it holds

2.3.
Well-Posedness of Navier-Stokes Equations with Neumann Boundary.In order to identify the limit dynamics of our interacting particle system with solutions of (NNS), we will need a uniqueness result and a priori estimates.We set the discussion of the limit PDE in the space of càdlàg functions since convergence of the particle system will take place in that space.
Remark 2.11.The latter is a good definition: thanks to Lemma 2.8, integrability in space of the nonlinear term ∇φ • K[ω]ω follows from Hölder's inequality, since both ∇φ and K[ω] are L q for all q > 2.
Proposition 2.12.Given ω 0 ∈ L 2 and g ∈ L 2 ([0, T ], L 2 ), there exists a unique weak solution of (NNS) in the sense of Definition 2.10 with ω(0) = ω 0 .Moreover, the unique solution belongs to C([0, T ]; The same statement holds for the cutoff PDE (c-NNS) for which in fact, thanks to the bounded nonlinearity, much stronger well-posendenss results can be obtained.
Proof.Consider two weak solutions ω, ω ∈ D([0, T ]; L 2 (D)), let w t = ω t − ωt .With standard passages (we refer for instance to [18,Section 3]) one can extend the weak formulation of the PDE to time-dependent tests (this makes use of the right continuity hypothesis on ω), then, taking φ t (s, x) = P t−s ψ(x), ψ ∈ H 2,2 transforms the PDE into the variation-of-constants form A similar formula holds for ωt .Therefore w t satisfies By Sobolev embedding, for all q > 2, hence, by Holder inequality, again for all q > 2, We can thus control nonlinear terms by means of ultracontractivity (2.2), Hölder inequality and Lemma 2.8: The uniqueness statement now follows from Grönwall's Lemma.
The last statement of the Proposition follows from existence of a weak solution belonging to ).This can be proved with a standard Galerkin approximation together with the usual energy estimate (see the beginning of the proof of Proposition 2.13 below), we refer to [29, Theorem 2.6] for details.
Proposition 2.13.Any weak solution of (NNS) in the sense of Definition 2.10 in As a consequence, there exists M 0 = M 0 (ω 0 , g) such that for all M ≥ M 0 the unique weak solution of (NNS) coincides with the one of the cutoff PDE (c-NNS) with the same initial and boundary data.
In order to establish this a priori estimate we first recall a convenient representation for boundary terms appearing in the weak formulation of (NNS).
We refer to [44,Chapter 4] for the (standard) proof.It follows that Proof of Proposition 2.13.For the sake of completeness we sketch the classical argument for the a priori estimate in the L 2 setting.By Proposition 2.12 we know that the weak solution actually belongs to C([0, T ]; The second term on the right-hand side vanishes, and the third one can be controlled thanks to the properties of the trace of functions in H 1 (D): Therefore, we obtain 1 2 from which by Grönwall's lemma we conclude that Arguing as in the proof of Proposition 2.12, extending the weak formulation of the PDE to time-dependent tests and taking φ t (s, x) = P t−s ψ(x), ψ ∈ H 2,p ′ , standard passages allow to rewrite the PDE into the variation-of-constants form (2.17) We control the nonlinear term by means of ultracontractivity (2.2), Hölder inequality and Lemma 2.8: for 2p p+2 < q < 2 < p it holds the last step following from Hölder inequality and the fact that t 0 (t−s) 2/q−2/p ds < ∞.As for the boundary term, by Proposition 2.3, Lemma 2.14 and (2.3), for δ > 0, where the integral in parentheses is finite if we choose 1 2p > δ.By Lemma 2.8, and combining the estimates above, we deduce that , from which the first statement of the Proposition follows.
As for the second statement, it suffices to take , where C T,p,ν is the constant implied in the previous inequality, making it so that the nonlinear term of solutions of (NNS) and (c-NNS) coincide by definition of the cutoff F .2.4.Stochastic Integrals and Burkholder-Davis-Gundy Inequalities.Let V be a separable Hilbert space, and consider a cylindrical Wiener process W on V with identity covariance operator; for our purposes V will in fact always be finitedimensional.Let (Ω, F , P) be the standard filtered probability space on which W is defined (and adapted), and consider a progressively measurable stochastic process (ψ t ) t∈[0,T ] taking values into HS(V, L 2 ), the space of Hilbert-Schmidt operators between V and L 2 = L 2 (D).Let us also assume that for all t ∈ [0, T ], P-almost surely, We refer to the classical monograph [13] for the definition of the stochastic integral t 0 ψ s dW s (a continuous L 2 -valued local martingale) and its basic properties.
The following collects (two versions of) the Burkholder-Davis-Gundy inequality, on which we will rely to bound uniformly in time terms due to the Brownian noise in the particle dynamics.Proposition 2.15.Under the assumptions of above, for any p > 0 it holds Moreover if (S t ) t≥0 is an analytic semigroup on L2 , for any p > 2 it holds The proof can be found in [41], to which we refer for a detailed and more general discussion on this kind of inequalities.

Definition of the Model and Main Results
Since we are considering a finite, fixed time horizon T > 0, for the sake of lightening notation we assume from now, without loss of generality, T = 1.
3.1.Generation of Particles.Let us introduce a grid of mesh 1/n spanning ({0}× D)∪((0, 1]× ∂D).Let γ : S 1 ≃ (0, 1] → ∂D be a diffeomorphism parametrizing the boundary of D; for n ≥ 1 we set We will write to enumerate the points of the grid L n , N (n) being the cardinality of L n .From now on we simply write N = N (n) implying the dependence on n.We now introduce a partition 2 of ({0} × D) ∪ ((0, 1] × ∂D) whose elements are centered at points ζ i : Notice that N (n) ≃ n 2 and the area of the set of indices relative to particles created before time t.
From our assumptions on the mesh of the grid, ω 0 and g it follows that

Diffusion Processes and Reflecting
Boundaries.Let (Ω, F , P) a complete, filtered probability space satisfying the standard hypothesis, on which it is defined a sequence (B i ) i∈N of independent F -Brownian motions.We have the following (probabilistically strong) well-posedness result: Proposition 3.1.For all n ≥ 1, on (Ω, F , P), there exists a unique continuous DN -valued adapted process x n (t) = (x n 1 (t), . . ., x n N (t)) t∈[0,T ] and a continuous R 2×N -valued adapted process k n (t) = (k n 1 (t), . . ., k n N (t)) t∈[0,T ] with bounded variation trajectories such that: for i = 1, . . ., N , t ∈ [t n i , 1], while for t < t n i we impose k n i (t) = 0 and x n i (t) = ζ n i .The proof can be straightforwardly adapted from the one given in [31,43] for systems of SDEs with regular coefficients and reflecting boundaries: the only difference is generation of new particles at the boundary, which is taken care of applying the well-posedness result on each time interval of length 1/n during which particles are not generated.We also refer to [29,Section 4] for details on a SDE system (closely related to ours) including boundary generation at random times.Itô's formula for the process x n (t) takes the following form (for which we refer again to [31]): Notice that the hypothesis ∇φ • n = 0 on ∂D makes it so that reflection terms −k n i do not appear in the Itô formula.In our applications, this assumption will be verified thanks to our choice of regularizing the empirical measure with the heat kernel under Neumann bouondary conditions.
Remark 3.3.The classical vortex dynamics in a bounded domain is a system of singular ODEs of the form ẋi = including self-interaction terms induced from the boundary effects; this is necessary for the system to satisfy in a weak form (as in [40]) the 2-dimensional Euler equations.We refer to [34,Chapter 4] for a general introduction to the topic, to [15,22,23,24] on the issue of well-posedness of the dynamics and vortex collisions, and to [32,25,26] on the statistical mechanics point of view.Self-interactions diverge logarithmically at ∂D, and this prevents us to include them in our model because the Brownian part of the dynamics (which we must include to model viscosity) might drive particles onto the boundary causing blowup of the dynamics at finite time.Nevertheless, under the Mean-Field scaling of particle intensities we are considering, self-interactions should be negligible in macroscopic limit.Indeed (heuristically) a self-interaction term ω i ∇ ⊥ γ(x i , x i ) in our model would scale as n −2 , while the nonlinear and noise terms are of order 1 as n → ∞.Hence, self-interactions due to the boundary appear to be irrelevant for the purpose of our discussion.
, ω 0 , g ≥ 0, and let the related notation introduced above prevail.There exists M > 0 (only depending on ω 0 , g) such that for every η ∈ 2 p , α , as n → ∞, the kernel-smoothed empirical measure ) to the unique weak solution (given by Proposition 2.12) of (NNS) with initial datum ω 0 and boundary source g.By Morrey's inequality, the stated convergence implies that on D([0, 1], C( D)).

Uniform Bounds
The proof of Theorem 3.4 essentially relies on uniformly bounding (in n) the approximating process ω n (t, x) in terms of strong norms, allowing the application of Aldous' tightness criterion, Proposition 2.9, and to pass to the limit the dynamics obtaining (NNS).The way we exploit the regularizing effect of the Brownian noise driving particles, that is the parabolic nature of the corresponding PDE dynamics, consists in formulating the evolution problem in Duhamel form (that is the variation-of-constants form, or mild formulation), in complete analogy with the a priori estimates for the limiting PDE in the previous section.
4.1.Duhamel Formulation of Empirical Measure Dynamics.We adopt, for the remainder of the Section, the notation of subsection 3.2.Let thus (x n 1 , . . ., x n N ) be the DN -valued stochastic process on [0, 1] defined by Proposition 3.1 as the unique strong solution of the approximating particle system, and A direct application of the Itô formula in Corollary 3.
Let us now considered the kernel-smoothed empirical measure which, in sight of (4.1), must satisfy for all t ∈ [0, 1] and x ∈ D As in the proof of Proposition 2.13, we can derive from the latter formulation of the dynamics the variation-of-constants form with standard passages.For all t ∈ [0, 1] and x ∈ D it holds: 4.2.Preliminary Estimates on Particle Creation and Diffusion.We begin by estimating separately the terms due to particle generation at t = 0 and at the boundary at later t > 0.
Proof.We have We bound the term I ω0 : (4.4) Concerning R, recall that, for i such that t n i = 0, We then have, by (2.13), .
Thanks to the gradient bound Lemma 2.6, we get Putting together the bounds (4.4) and (4.5), we conclude that By the assumption ε = ε(n) n −2/(3+α) , we get the desired bound.
As for the previous lemma, the idea of the proof is that, for n large and ε small, and the right-hand side is bounded if, say, g is bounded.
Proof.We have We bound the term I g first: by (2.10), we get t 0 ∂D P t−s p ε (•, y)g s (y)dsdσ(y) ds.
Using the trace theorem and the contractivity bound (2.3), we get, for some fixed δ > 0 such that α/2 + 1/(2p ′ ) + δ/2 < 1, Hence we can split R as To bound the term R 1 , we start with an observation: setting ∈ D for all ξ since D is convex), since by definition of heat semigroup it holds ∂ t P t = −A p P t , we can write Applying the latter to R 1 we obtain: Concerning R 11 , recall that |s − t n i | ≤ 1/(2n) for every s in Q n i .Hence by (2.12) we have g s (y)dsdσ(y)dξ.
Since H 2,p ′ is embedded into C( D) for p ′ > 1, we have Note that, if t n i + 1/(2n) < t, then t − [t n i , s](ξ) ≥ (t − s)/2 for every s in Q n i .Thanks to the contractivity bound (2.3), we get Hence by (2.13) we have g s (y)dsdσ(y)dξ.
Thanks to the gradient bound Lemma 2.6, we get We turn now to the bound on R 21 .By (2.13) we have g s (y)dsdσ(y).
Notice that i: Therefore, by the embedding H Finally, we turn to the bound on R 22 .Similarly to the case of R 21 , we notice that for each i ∈ A 0 (t), Hence, proceeding as for the term R 21 , we get Putting together the bounds (4.6), (4.7), (4.8), (4.9), (4.10), we conclude that .
We also estimate in a dedicated Lemma the martingale terms appearing in the mild formulation of particle dynamics.
Proof.By Sobolev embedding H 1−2/p,2 ֒→ L p and Burkholder-Davis-Gundy inequality (2.19) it holds The inner integrand in the right-hand side can be controlled uniformly with respect to the particles' positions: by Lemma 2.6 Notice that the denominator is integrable in ds for all δ > 0. We now take into account the fact that ) uniformly in i, and that N (n) ≃ n 2 .Hence, the left-hand side of (4.11) is bounded from above by Taking ε(n) as in the hypothesis, the last quantity is uniformly bounded in n.

4.3.
Estimates for Approximating Solutions.The following is the core estimate of our argument, providing uniform boundedness in strong norms for the approximating processes.
Proposition 4.4.Let p > 2, 2 p < α < 1, and ε = ε(n) n −1/2 .For all q ≥ 2 it holds (4.12) sup Proof.In order to lighten notation we denote by v n t (x) := F D K n (x, y)dS n s (y) the vector field acting on a single, smoothened particle.The mild formulation (4.3) allows to bound, by Minkowski inequality, The initial and boundary creation terms I 1 and the martingale term I 3 are uniformly bounded (both in n and in t ∈ [0, 1]) respectively thanks to Lemma 4.1, Lemma 4.2 and Lemma 4.3.We thus focus on the nonlinear interaction term I 2 , where the regularizing effect of the heat semigroup is essential.Using the gradient estimate (2.5) we thus bound, for f ∈ L p ′ , i∈A n (s) By duality, this implies that i∈A n (s) 1/q ds with the constant C ε decreasing as ε → 0, in particular uniformly bounded in n.
Altogether, we arrive to the integral inequality 1/q ds which, by Grönwall's lemma (as α < 1), implies In order to complete the proof, we apply again the mild formulation of the dynamics (4.3) to estimate Once again, the initial and boundary creation terms J 1 and the martingale term J 3 are uniformly bounded in n respectively thanks to Lemma 4.1, Lemma 4.2 and Lemma 4.3.As for J 2 , since we can repeat the computations performed to control I 2 obtaining where we now apply Hölder inequality, , in which we can choose a small enough r > 0 (depending on α) so that the first factor on the right-hand side is finite and q/r ′ < 1, hence which is uniformly bounded in n by (4.13).
It is worth noticing that the gradient estimates for P t are crucial in controlling the nonlinear term of the dynamics, but the regularizing effect of P t was neglected in estimating the initial empirical measure in Lemma 4.1 and the martingale terms in Lemma 4.3 (due to the application of BDG inequality).The singularity of ∇p ε appearing in the initial empirical and in the stochastic integrals is not improved by the heat semigroup and produces a restriction on the asymptotic behavior of ε = ε(n).
Thanks to the good control provided by Proposition 4.4, we can estimate time increments -which we need for the equicontinuity part of our compactness argumentin a much weaker norm, thus allowing to exploit the weak formulation of the dynamics, easier to deal with compared to the variation-of-constants form but producing bounds in weaker norms.
Proof.Substituting the weak formulation of particle dynamics (4.2), we estimate by Minkowski inequality: We now proceed estimating each term separately.
Estimates on generation term I 1 .
The first factor on the right-hand side of the latter is controlled by so we are left to control the second factor.For f ∈ L p ′ , y ∈ ∂D and a sequence y k ∈ D converging to y, it holds the exchange between the limit and the integral being allowed thanks to the fact that, by Proposition 2.4, p ε (x, y) ≤ C ε e −|x−y| 2 /(cε) ∀x, y ∈ D. We can thus estimate: p 1 follows by Morrey's inequality.
Estimates on the diffusion term I 2 follow from a straightforward application of Hölder inequality, Estimates on the nonlinear term I 3 .We have: ds , in which we control, somewhat analogously to Proposition 4.4, Therefore I 3 is bounded, up to some constant independent of n, by Estimates on the martingale term I 4 .By Sobolev embedding and Burkholder-Davis-Gundy inequality, Proposition 2.15, uniformly in i, and replacing the expectation involving particles's positions with a supremum over D, we obtain , which we combine with a consequence of Lemma 2.6, and the assumption ε(n) n −1/2 , concluding that I 4 ν,p,q r q/2 .

Proof of the Main Result
The proof of Theorem 3.4 proceeds as follows.We combine the estimates of the previous section with Aldous' Lemma in order to obtain tightness on the space of càdlàg functions taking values in Sobolev spaces.In order to verify that the limiting dynamics coincides with the Navier-Stokes equations (NNS) we need almost sure convergence in such a space, which we obtain by changing the underlying probability space by Skorohod theorem.Convergence in the original probability space is recovered by uniqueness of the deterministic limit.
In what follows, according to the hypothesis of Theorem 3.4, we tacitly assume that p > 2 and 2 p < α < 1 are fixed, as well as ω 0 ∈ Lip( D) and g ∈ B b ([0, 1] × ∂D).We denote by ω the unique solution of (NNS) given by Proposition 2.12 with initial datum ω 0 and Neumann source g.

5.
1. Compactness Argument.Let (Ω, F , P) be a complete, filtered probability space satisfying the standard hypothesis, on which it is defined a sequence (B i ) i∈N of independent F -Brownian motions.For all n ≥ 1, on this same probability space we can consider the well-posed dynamics of Proposition 3.1, since the latter is a probabilistically strong existence and uniqueness result.In other words, we can consider for all n the dynamics of particles (x n 1 (t), . . ., x n N (t)) t∈[0,1] and the one of the regularized empirical measure ω n t = i∈A n (t) ω n i p ε (•, x n i (t)) as stochastic processed defined on (Ω, F , P) and F -adapted.
We denote by L n the law of ω n on D([0, 1], H α,p ); we can actually consider any parameters α, p since samples of the process ω n are smooth in the space variable, but time dependence is at best càdlàg due to the creation of new particles.
We need first to show tightness of the laws (L n ) n∈N : Lemma 5.1.For all 2 p < η < α, the sequence of laws (L n ) n∈N is tight on D([0, T ]; H η,p ).
We then need to show that the limiting dynamics is actually a solution of (NNS).Lemma 5.2.For all 2 p < η < α, any weak limit point of the sequence (L n ) n∈N of measures on D([0, T ]; H η,p ) is concentrated on the unique solution ω of (NNS) given by Proposition 2.12; in other words there exists a unique limit point, δ ω .
These two Lemmas imply Theorem 3.4.
Proof of Theorem 3.4.By Lemma 5.1 and Lemma 5.2, every subsequence L n(k) admits a sub-subsequence which converges to the unique limit point δ ω , where ω is the deterministic solution of (NNS).Then, for example by [7,Theorem 2.6], the whole sequence L n converges weakly to δ ω , and then also in probability as the limit point is a Dirac delta (see e.g. the argument in [7, page 27]).The proof is complete.5.2.Proof of Lemma 5.1.Thanks to the uniform estimates of section 4, we can straightforwardly apply Aldous' criterion, Proposition 2.9, to the processes ω n t on D([0, T ]; H η,p ).
Condition (1), that is tightness of laws at fixed time in H η,p , actually holds for all t ∈ [0, 1].This is a trivial consequence of Proposition 4.4, Markov's inequality and the fact that H η ′ ,p ֒→ H η,p with η ′ > η is a compact embedding.
Condition (2) of Proposition 2.9 is the harder one.Let us show that (2.14) holds: given a sequence (τ n ) n of F -stopping times we estimate, for q > 2, , where the second step is the interpolation inequality between H −2,p and H α,p , and the third is Hölder inequality with exponents 1 < u, u ′ < ∞.A simple computation reveals that, since q > 2, one can choose u so that the exponents of the norms in the expected values are strictly larger than 2, therefore we can now apply the uniform estimates of Proposition 4.4 and Proposition 4.5, obtaining > δ ν,M,α,p,q δ −q r uq(α−η) It is now possible, given δ > 0, to choose r 0 > 0 small enough and n 0 large enough (taking then r ≤ r 0 and n ≥ n 0 ) so that the right-hand side of the latter inequality is arbitrarily small, thus satisfying Aldous' condition and implying tightness of (L n ) n∈N .
5.3.Proof of Lemma 5.2.We consider a weak limit point L in D([0, T ]; H η,p ); for simplicity, we still denote by L n the convergent subsequence.
Consider the joint law of ω n and the sequence of Brownian motions (B i ) i∈N0 on the product space (R 2 ) N × D([0, T ]; H η,p ) (which we regard as a product metric space, considering a distance on (R 2 ) N that makes it separable).In particular, the relation between Brownian motions and ω n is encoded in the weak formulation of the dynamics: P-a.s. for each t ∈ [0, 1] and φ ∈ H s+2,p ′ with sp ′ > 2 (and 1 < p ′ < 2 since p > 2).By Skorohod coupling theorem on (R 2 ) N × D([0, T ]; H η,p ), [30,Theorem 3.30], there exists a complete probability space ( Ω, A, P) with random elements (ω n ) n∈N , ( Bi,n ) i∈N having the same joint laws of (ω n ) n∈N , (B i ) i∈N , and ω with law L = L , such that ωn converges P-almost surely to ω. Moreover we can take the filtration F n,0 = ( F n,0 t ) t generated by (ω n , ( Bi,n ) i ) and the P-null sets on (Ω, A) and we define the filtration F n by F n t = ∩ s>t F n,0 t ; by a standard argument, see e.g. the proof of [4, Proposition 2.5, part 1], we get that (ω n and ( Bi,n ) i are progressively measurable with respect to F n and) ( Bi,n ) i is still a cylindrical Brownian motion with respect to F n .Arguing for example as in [6,Section 4.3.4](see also [14,Proposition 4.1] or [19,Section 2.4.4]),we can show that (5.1) is satisfied by ωn , ( Bi,n ) i=1,...,N on the filtered probability space ( Ω, F n , P).
If we now show that L = δ ω , so that also L = δ ω , we obtain the first statement of the Lemma.In proving such claim, we drop tilde signs from ωn , ( Bi,n ) i to lighten notation 3 .The aim is to pass to the limit (5.1), we do so term by term.We then need to prove that ω n 0 , φ → ω 0 , φ , that is, ω0 = ω 0 ; we do so by showing that as n → ∞, the following quantities vanish: 3 The tilde on our limit point ω is retained since a priori it can differ from the solution ω of (NNS).
We treat I 1 exploiting the fact that φ is Lipschitz continuous, As for I 2 , using the regularity of ω 0 we assumed by hypothesis, Let us define I n (t) = (t n j − 1/(2n), t n j + 1/2n] with t n j such that t ∈ I n (t).We control J 1 exploiting again the fact that φ is Lipschitz, Once again, the second term J 2 is the easier one: . We now proceed to show that To do so, we add and subtract the same quantity three times, apply the triangular inequality and estimate differences.For the sake of shorter formulas and clearer passages we will write Notice that T 1 and T 5 coincide with the terms of the difference in (5.2), therefore the latter vanishes if The first difference is the harder one, since it involves the pointwise difference between the empirical measure and its regularization: From there, since F is Lipschitz continuous, and applying the pointwise estimate on heat kernel (2.4), The other differences now are controlled by the almost sure convergence assumption and functional analytic estimates already repeatedly applied.Bounding L p ′ , it is clear that differences on the right-hand side converge P-almost surely to 0 as n → ∞.

5.3.5.
Removing the cutoff.The argument detailed so far implies that the limit ω is a weak solution, thus the unique weak solution of the cutoff PDE (c-NNS).However, by Proposition 2.13 (second statement) it is now possible to choose M large enough so that K and thus ω = ω is actually the unique weak solution of the PDE (NNS) without cutoff.This concludes the proof of the first part of Lemma 5.2.

General Initial and Boundary data
In this section we provide a brief description of changes to the above arguments required to deal with general boundary data ω 0 ∈ Lip( D) and g ∈ B b ([0, 1] × ∂D), dropping the non-negativity assumption.The idea is simply to divide boundary data into positive and negative parts, and to consider particle and PDE dynamics for the two parts, properly taking into account their interactions.
6.1.Splitting of PDE and particle dynamics.Let us set )) ∨ 0. The limit dynamics (NNS) is then equivalent to a coupled system of PDEs, (6.1) The notion of weak solution to the latter system is completely analogous to the one of Definition 2.10.Moreover, a couple (ω + , ω − ) is a weak solution of the system above if and only if the couple (ω = ω + − ω − , ω + ) is a weak solution of (6.2) Well-posedness of (6.2) follows easily from Proposition 2.12.Therefore, also (6.1) is well posed in D([0, T ]; L 2 (D) 2 ).As for the approximating dynamics, coherently with notation of subsection 3.1, we define From our assumptions on the boundary data and on the mesh of the grid, it follows that ω n,± i 1 n 2 uniformly in i = 1, . . ., N , for all n.We introduce two mutually independent sequences of independent F -Brownian motions (B i,± ) i∈N and particles with positions x n,± i (t) satisfying the following evolution equations: for t > t n i , x n,± where k n,± i are continuous, adapted, R 2 -valued processes with bounded variation trajectories, whereas for t ≤ t n i , x The following statement extends Theorem 3.4 to the general case of ω 0 , g not necessarily non-negative.
In order to obtain the required compactness we need to show that ω n,− t q H α,p < C M,q,p,ν,α , E ω n,+ τn+r − ω n,+ τn q H −2,p + E ω n,− τn+r − ω n,− τn q H −2,p M,q,p,ν,α r q/2 + 1 n q . (6.7) The proof of these inequalities follows the same strategy described in section 4. Indeed, linear terms, stochastic integrals and generation terms can be treated exactly as in the case of positive data.In fact, nonlinear terms present no additional difficulties and the crucial estimate for the nonlinear term of Proposition 4.4 is easily adapted.For f ∈ L p ′ such that f L p ′ = 1 the following chain of inequalities holds: Then the proof of (6.6) follows exactly as in Proposition 4.4.Similarly, in the proof of (6.7) the only differences concern the nonlinear terms, and they can be handled with analogous simple modifications.
Once uniform bounds are recovered, one can proceed to replicate the tightness argument of subsection 5.1.Passing to the limit dynamics is similar to that explained in subsection 5.3, once again we only outline the (slightly) different treatment required by nonlinear terms.
n,± i (t) = ζ n i and k n,± i (t) = 0. Well-posedness of this particle system follows from the same arguments proving Proposition 3.1.Similarly to section 4 we introduce the empirical measures S ), x ∈ D.

6. 2 .
Uniform Estimates for split dynamics.The regularized empirical measures satisfy the integral equations and the mild formulations below: ) ∇ y p ε (x, y)dS n,) ∇ y p ε (x, y)dS n,± s
The expression A ≃ a,b B indicates that B is both an upper and lower bound for A up to multiplicative constants depending only on possible subscripts a, b.Expressions A a,b B or A a,b B indicate respectively an upper and lower bound in the same sense.