The Dean–Kawasaki Equation and the Structure of Density Fluctuations in Systems of Diffusing Particles

The Dean–Kawasaki equation—a strongly singular SPDE—is a basic equation of fluctuating hydrodynamics; it has been proposed in the physics literature to describe the fluctuations of the density of N independent diffusing particles in the regime of large particle numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\gg 1$$\end{document}N≫1. The singular nature of the Dean–Kawasaki equation presents a substantial challenge for both its analysis and its rigorous mathematical justification. Besides being non-renormalisable by the theory of regularity structures by Hairer et al., it has recently been shown to not even admit nontrivial martingale solutions. In the present work, we give a rigorous and fully quantitative justification of the Dean–Kawasaki equation by considering the natural regularisation provided by standard numerical discretisations: We show that structure-preserving discretisations of the Dean–Kawasaki equation may approximate the density fluctuations of N non-interacting diffusing particles to arbitrary order in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^{-1}$$\end{document}N-1 (in suitable weak metrics). In other words, the Dean–Kawasaki equation may be interpreted as a “recipe” for accurate and efficient numerical simulations of the density fluctuations of independent diffusing particles.


Introduction
The Dean-Kawasaki equation -with ξ denoting a vector-valued space-time white noise -has been proposed by Dean [7] and Kawasaki [26] as a model for density fluctuations in a system of N indistinguishable particles undergoing diffusion in the regime of large particle numbers N 1. Its mathematical analysis is complicated by its highly singular nature: A scaling argument shows that (1) is not renormalisable by the theory of regularity structures by Hairer et al., even in one spatial dimension d = 1.
Recently, Konarovskyi, Lehmann, and von Renesse [29] have obtained a striking rigidity result for the Dean-Kawasaki equation (1): They show that the only martingale solutions to (1) are of the form of an empirical measure for N independent Brownian motions where the {w k } N k=1 denote the N independent Brownian motions. In particular, no solution exists for non-integer values of N . This result may be viewed as casting doubts on the mathematical meaningfulness of the Dean-Kawasaki equation: It amounts to stating that the Dean-Kawasaki equation is just a mathematically complex way of rewriting the diffusion of N particles. In turn, this naturally raises the question what the advantages of the Dean-Kawasaki equation (1) might be over the particle formulation of diffusion (2) from the point of view of physics.
In the present work, we provide a rigorous justification for the Dean-Kawasaki equation: We show that standard numerical discretisations of the Dean-Kawasaki equation (1) -such as finite difference or finite element schemes -provide accurate descriptions of the density fluctuations in a system of N diffusing particles if measured in suitably weak metrics. Roughly speaking, we show that, under certain conditions, the solutions ρ h to the discretised Dean-Kawasaki equation achieve the approximation quality where j ∈ N is arbitrary, h is the spatial discretisation parameter, p + 1 is the order of convergence of the numerical scheme in the Sobolev space H −1 , and d weak,2j−1 is a suitable weak metric of negative Sobolev type. In particular, the accuracy is of arbitrarily large order in N −1/2 and hence only limited by the numerical discretisation error and the negative part ρ − h . In addition, we show that E ρ − h decays exponentially fast inroughly speaking -(hN 1/d ) 1/2 , demonstrating that the term becomes quickly negligible in the scaling regime (where we have dropped logarithmic corrections in N and contributions on the final time horizon for the sake of this introductory exposition). In a nutshell, the bound (3) implies that the Dean-Kawasaki equation can be used as a "recipe" for accurate simulations of density fluctuations in systems of diffusing particles. Note that our scaling regime (4) is not an actual restriction in the context of numerical simulations: It ensures that the average number of particles per cell is strictly larger than one. In fact, in the opposite regime h ≤ N −1/d , the direct simulation of particles would become less expensive than the approximation by the Dean-Kawasaki equation, as the numerical effort for the Dean-Kawasaki equation is strictly larger than h −d .
While the Dean-Kawasaki equation correctly describes the fluctuations around the mean-field limit to arbitrarily large order in N −1/2 , the well-known linearised description of fluctuations given by the solutionρ to is limited to the approximation quality Here,ρ denotes the mean-field limit given as the solution to    ∂ tρ = 1 2 ∆ρ, ρ(·, 0) = ρ 0 .
In fact, the linearised descriptionρ only captures the leading-order fluctuation correction to the mean-field limit correctly and hence carries a relative error of order N −1/2 with respect to the fluctuation scaling. We provide numerical evidence of such difference between the two models, and we also numerically verify convergence rates for suitable discretisations of the Dean-Kawasaki model (1).
1.1. Related Literature. The theory of fluctuating hydrodynamics describes fluctuations in interacting particle systems in the regime of large particle numbers using suitable SPDEs; see e. g. [37]. In its form (1), the Dean-Kawasaki equation describes noninteracting particles, with similar equations being available for weakly interacting particles undergoing overdamped Langevin dynamics. In the recent contribution [12], the authors also address quantitative fluctuation bounds in the non-interacting particle case, but by means of a suitable approximated SPDE model rather than a numerical scheme. While their setting grants several well-posedness results (non-negativity of the solution, comparison principles, entropy estimates) and allows to consider initial particle profiles which are more general than those treated here, their relative fluctuation error is however limited to N −1/(d/2+1) log N , while the rate of fluctuations in (3) is -in suitable metrics -arbitrarily high. For a more general particle setting, the SPDE of fluctuating hydrodynamics for the zero range process given by ∂ t ρ = ∆(Φ(ρ)) + ∇ · Φ(ρ)ξ (6) has been addressed in [11], and linked it to the large deviation principle for such process in a suitable thermodynamic setting. A corresponding well-posedness result for truncated (low spatial frequency) noise and regularised nonlinearity has been obtained in [18], see also [19]. In [20], the construction of random dynamical systems for conservative SPDE is discussed, together with well-posedness theory of invariant measures and mixing of the related Markov process. In [17], a large deviation principle for regularised variants of (6) is shown; in a suitable limit, the rate functional of such large deviations principle and the corresponding one of the interacting particle system are shown to approach each other. The paper [10], written independently of -and simultaneously to -the present manuscript, gives a rigorous justification of the fluctuating hydrodynamics SPDE associated with the symmetric simple exclusion process While in contrast to our work the authors in [10] only consider the continuum SPDE, they regularise it by truncating the noise for small spatial wavelengths. In a certain sense, this introduces a regularisation in the same spirit as our numerical approach; however, their truncation criterion is somewhat more restrictive than our condition h N −1/d . While they face a more challenging problem with the more complex noise intensity factor ρ(1 − ρ) (whose square is not a linear function of the density ρ) and also prove convergence results for the rate functions for large deviation principles, they only establish a leading-order description of fluctuations in the low deviations regime. In other words, in contrast to our present work, they do not show superiority of fluctuating hydrodynamics to a linearised approach on fluctuations for the "bulk" of the probability distribution. For recent numerical approaches to fluctuating hydrodynamics, we refer the reader e. g. to [36,9,14,3,2,24,13,33,1] (in particular, [2] contains the extension of the current work to the case of weakly interacting particles). Note that the small prefactor of the noise term in the Dean-Kawasaki equation (1) enables the use of certain higher-order timestepping schemes [22], a fact that we also use in our numerical simulations.
Concerning further mathematical results on Dean-Kawasaki models, the link between fluctuating hydrodynamics and Wasserstein geometry has long been understood, and extensively studied in several works, see for instance [25,38,28,31,32,30,8].
Dean-Kawasaki type models including the effects of inertia have been derived and analysed by the first author, Shardlow, and Zimmer [4,5,6].
The fluctuation-dissipation relation -implicitly contained for instance in the Dean-Kawasaki equation -may be used to recover macroscopic diffusion properties from fluctuations in finite particle number simulations, see for instance [16,33]. Outside of the realm of physics, the concept of fluctuating hydrodynamics has also been applied to systems of interacting agents, see e. g. [24,13,27].
Finally, conservative stochastic PDEs have recently been shown to give optimal convergence rates in the mean-field limit approximation of stochastic interacting particle systems, such as those encountered in the stochastic gradient descent methods for overparametrised, shallow neural networks [21]. Remark 1. Given the nature of the metric d weak,2j−1 in (3), it is natural to ask whether or not the high-order fluctuation error of (3) could be formally derived from suitable a priori estimates of negative Sobolev type for the continuous Dean-Kawasaki model (1). An a purely formal level, testing the mild formulation of (1) -where G is the heat semigroup kernel -with trigonometric functions, and performing elementary computations, one arrives at the a priori estimate which is valid in the regime j > d/2. Despite its formal validity -which, however, relies on the non-trivial negativity requirement for the density ρ -the inequality (8) does not give any information beyond the leading order N −1 , and therefore is too weak an estimate to justify the high-order fluctuation error bound in (3).

Main results and summary
The methodology of this paper can be applied to several standard numerical discretisations of the Dean-Kawasaki model (1), including finite difference and finite element schemes. In the interest of brevity, we only focus on a finite difference discretisation: The corresponding results in the finite element case are given in Appendix B. Specifically, on the periodic domain T d := [−π, π) d , we denote the uniform square grid with spacing h by G h,d , the discrete inner product of L 2 (G h,d ) by (·, ·) h , the interpolating operator on G h,d by I h , and define the distance d −j [X, Y ] between two R M -valued random variables as Our first main result reads as follows.
Theorem 2 (Accuracy of description of fluctuations by the finite-difference discretised Dean-Kawasaki model of order p + 1 ∈ N). Assume the validity of Assumption FD1 (discretised differential operators), Assumption FD2 (Brownian particle system), Assumption FD3 (scaling assumptions), and Assumption FD4 (discretised mean-field limit), all given below. In particular, assume that the mean-field limit ρ h in (20) Then, for any j ∈ N, the discrete Dean-Kawasaki model FD-DK captures the fluctuations of the empirical measure µ N in the sense that, for any T = (T 1 , . . . , =: Err neg + Err num + Err f luct,rel holds for any ϕ = (ϕ 1 , . . . , ϕ M ) ∈ [W p+Θ+j+1,∞ (T d )] M such that ϕ m L 2 = 1, ∀m = 1, . . . , M and´T d ϕ k ϕ l dx = 0 whenever T k = T l . Finally, we have the a posteriori bound where we have set We make some observations in order to better illustrate the meaning of Theorem 2: • The quantities (ρ h (T m ), I h ϕ m ) h , and µ N Tm , ϕ m are rescaled with the factor N 1/2 , as the natural order of density fluctuations is N −1/2 . In other words, our main error estimate basically provides an estimate for the relative error in the fluctuations.
• The distances d −j [X, Y ] correspond to negative Sobolev norm differences of the probability measures on R M given by the laws of X and Y . In particular, it holds where W 1 is the 1-Wasserstein distance. • The above estimates contain three types of error terms. The term Err neg quantifies the a priori lack of knowledge concerning non-negativity of the solution ρ h ; the term Err num encodes the numerical precision of the scheme; finally, the term Err f luct,rel bounds the relative error in the fluctuations.
• The order of differentiation required for the functions ϕ should be thought of as the sum of p + 2 + Θ (accounts for the requirements of the spatial discretisation, discussed below) and j − 1 (necessary due to an induction argument over j).
If one is only interested in moment bounds (i.e., in a polynomial ψ) then the following estimate with no relative error in the fluctuations can be produced.
Theorem 3 (Estimates on the error for stochastic moments). In the same setting of Theorem 2, fix times T = (T 1 , . . . , Then the difference of moments between ρ h and the empirical density µ N (2) reads with constants C, C 1 , . . . , C 4 > 0 independent of j, h, N , T , and where we have the bound where E(N, h) has been defined in (11).
2.1. Structure of the paper. Section 3 lays out the finite difference discretisation of the Dean-Kawasaki model. Subsection 3.1 (respectively, Subsection 3.2) lays out the necessary notation (respectively, the relevant assumptions and definitions) related to the model. Subsection 3.3 -which has an informal flavour -brings forward some of the main ideas used in the paper. This section lays the ground for Subsection 3.4 (respectively, Subsection 3.6), which contains preparatory results for the proofs of Theorem 2 (respectively, Theorem 3). The proof of Theorem 2 (respectively, Theorem 3) is finalised in Subsection 3.5 (respectively, Subsection 3.7). Technical details are deferred to Subsection 3.8 (bounds for all moments of ρ h , and exponentially decaying bound for the negative part ρ − h ), and Appendix A (deterministic finite difference arguments and relevant Itô calculus). The statements of results for finite element schemes are given in Appendix B. Finally, Section 4 contains numerical simulations associated with Theorem 2, using a first-order finite difference discretisation (i.e., p = 1) in the one-dimensional case d = 1.

Analysis for finite difference discretisations
3.1. Notation. Domain and function spaces. Let N d ≤ 3, and let T d := [−π, π) d . Let h := 2π/L, for some L ∈ 2N, be the discretisation parameter of the periodic square grid We always work with periodic functions (defined either on T d or G h,d ). From now on, this fact will be implicitly assumed and no longer stated. In particular, we abbreviate C β = C β (T d ) and W r,p = W r,p (T d ). We use bold characters to denote vector fields. For and admits an orthonormal basis {e m x, } (x, )∈(G h,d ,{1,...,m}) , whose elements are defined as ] m as the function agreeing with φ on G h,d . When there is no ambiguity, we simply write φ instead of I h φ.
Discrete differential operators. We use the notation ∂ h,x to denote a finite difference operator approximating the partial derivative ∂ x . We denote by ∇ h := [∂ h,x1 , . . . , ∂ h,x d ] the associated finite difference gradient operator. Furthermore, for each , we define the discrete second partial derivative D 2 h,x as the operator for which the standard integration by parts formula holds, where D h,x is some (possibly different) finite difference operator approximating the partial derivative ∂ x . We abbreviate ∇ D,h := [D h,x1 , . . . , D h,x d ]. As a result of (13), the discrete operators D 2 h,x are symmetric (in the sense of finite difference operators). We abbreviate h,x to indicate the discrete Laplace operator. Specific details on ∇ h and ∆ h will be provided in the following subsection.
Remark 4. The operators ∇ h and ∇ D,h (both providing an approximation of the continuous gradient ∇) may be different, and have different uses in our discretised Dean-Kawasaki model (Definition FD-DK below). The operator ∇ h is deployed in the noise, while the operator ∇ D,h in the integration by parts formula (13).
For reasons which will become clear in Subsection 3.3 (see Block 3 therein), we set the notation for suitable continuous and discrete backwards heat flows. Specifically, for a sufficiently regular function ϕ and a final time T , we denote by φ t the solution the continuous backwards heat equation with final datum φ T = ϕ. Analogously, we denote by φ t h the solution to the discrete backwards heat equation with final datum φ T h = I h ϕ. In the following, we also use the alternative notation , to stress that P z (ϕ) (respectively, P z h (I h ϕ)) is the result of evolving a backwards heat equation (respectively, a discrete backwards heat equation) starting from ϕ (respectively, from I h ϕ) for a timespan z.
For y ∈ R, we define y + := max{y; 0} and y − := − min{y; 0}. In addition, as usual, we use the letter C to denote a generic constant, whose value may change from line to line in the computations.

Assumptions and discretised Dean-Kawasaki model.
Assumption FD1 (Discrete differential operators). Let p ∈ N be fixed. We make the following assumptions on the discrete operators ∂ h,x and D 2 h,x : • the discrete operators ∂ h,x and D 2 h,x are finite difference operators of order p + 1. Explicitly, this means that for any φ ∈ C p+2 (T d ); • The operators ∂ h,x and D 2 h,x commute. Assumption FD2 (Brownian particle system and initial datum of Dean-Kawasaki dynamics). Let p be as in Assumption FD1. We assume to have N ∈ N independent d-dimensional Brownian motions {w k } N k=1 moving in T d . Moreover: • the initial positions {w k (0)} N k=1 are deterministic; • there exists a deterministic function ρ 0,h ∈ L 2 (G h,d ) (which will serve as the initial datum of the discretised Dean-Kawasaki dynamics in Definition FD-DK below), satisfying the following properties: there exist h-independent constants ρ min and ρ max such that the empirical density of the initial configuration µ N 0 := N −1 N k=1 δ w k (0) approximates ρ 0,h with accuracy p + 1, in the sense that the inequality holds for each function η ∈ C p+1 .
Assumption FD3 (Scaling of relevant parameters). We assume the scaling for some T > 0, and where ρ min and ρ max have been introduced in Assumption FD2. This scaling will be needed to produce an exponentially decaying estimate associated with ρ − h , see (66) below.
Assumption FD4 (Mean-field limit). The solution to the discrete heat equation is such that ρ min ≤ ρ h ≤ ρ max (where ρ min and ρ max have been introduced in Assumption FD2) for all times up to T (where T has have been introduced in Assumption FD3).
We can now state the precise definition of our finite difference Dean-Kawasaki model.
Definition FD-DK (Finite difference Dean-Kawasaki model of order p + 1). Assume the validity of Assumptions FD1-FD4. We say that the L 2 (G h,d )-valued process ρ h solves a finite difference Dean-Kawasaki model of order p + 1 if it solves the system of stochastic differential equations where {β (y, ) } (y, )∈(G h,d ,{1,...,d}) are standard independent Brownian motions, and where F ρ ∈ L 2 (G h,d ) is defined as Remark 5. If (20) admits a discrete maximum principle, then Assumption FD4 is satisfied for any T > 0 and any non-negative datum ρ 0,h . For example, the discrete maximum principle applies for the second-order symmetrical discrete Laplace operator where y ∼ x indicates that y and x are adjacent grid points.
Remark 6. One may also omit the contribution (1+T ) in the scaling (19), at the expense of obtaining results with a worse dependency on the final time T . We are not interested in optimising time dependencies in this work, and we simply include the term 1 + T in order to get cleaner final results.
3.3. Key ideas behind the proofs of the main results. The proofs of Theorems 2 and 3 are of inductive type. In order to simplify their exposition, it is useful to first list a skeleton of the main building blocks.
Block 1. Discrete Dean-Kawasaki model: cross-variation analysis. At their core, both proofs use basic Itô calculus to describe the time evolution of suitable nonlinear functionals ψ of the quantities and of their expected values, where φ h and φ are suitable test functions. The quantities in (24) are linear functionals of ρ h and µ N , respectively. What is crucial, is that the cross-variation of the processes (24) are -up to a small error -also linear functionals of ρ h and µ N . The argument for µ N is straightforward, and we can thus defer it to the proofs themselves. As for ρ h , we use Definition FD-DK to write for two different test functions φ i,h , i ∈ {1, 2}. Using the Itô formula and the Parseval identity in [L 2 (G h,d )] d , one finds that the cross-variation of the stochastic noise of (25) is The first term in (29) is indeed a linear functional of ρ h . The second term (which we will show to be negligible for suitable scaling regimes, see Subsection 3.8) takes into account the a priori lack of knowledge concerning the non-negativity of solutions to the discrete Dean-Kawasaki model (21). We also stress that the validity of the computations above is independent of the order of the finite difference scheme (i.e., p). Expression (28) crucially preserves the cross-variation structure associated with the continuous Dean-Kawasaki (1) for nonnegative densities. More precisely, formally testing where the last inequality if justified by the representation ξ = s∈Z d e sβs , where {e s } s∈Z d is an orthonormal basis of [L 2 (T d )] d and {β s } s∈Z d are independent Brownian motions.
The noise cross-variation is then obtained using the Itô formula and the Parseval idendity and thus the cross-variations (31) and (28) are (modulo positive part ρ + h ) structurally identical.
Block 2. Numerical error. There are two contributions to the numerical error, namely: -the difference of initial data µ N 0 and ρ h (0), and -the difference in the evolution of test functions (say, φ and φ h ), and both are proportional to h p+1 . While the first contribution has the correct bound by Assumption FD2, the second contribution needs to be estimated: The main difficulty is that the interpolation of the test function arising from the cross-variation of the second quantity in (24) , the cross-variation of the first quantity in (24)). We therefore need to show the bound in order not to lose h-regularity in consecutive steps of our inductive proofs (more details in Block 5 below). The necessary tools for this task are contained in Subsection A.1.
Block 3. Deterministic dynamics of the test functions. As we are interested only in the analysis of the fluctuations for the Dean-Kawasaki model, it is convenient to choose the deterministic functions ψ, φ, φ h in such a way that as many drift terms as possible in relevant Itô differentials vanish. This is the reason behind the choice of the backwards heat equation (14) (respectively, (15)) for φ (respectively, for φ h ), which directly compensates the diffusive nature of the particle system (respectively, of the Dean-Kawasaki model). In practice, this is reflected in the useful equalities (which follow from Lemma 15) for φ, φ h as in (14), (15). The discussion for ψ in the case of Theorem 2 is conceptually analogous, but technically more involved, and is devolved to the proof itself. As for Theorem 3, ψ is chosen to be static, therefore this discussion does not apply. We expand these considerations in Appendix A.2.
Block 4. Stretched exponential bounds for centred moments of the particle system and the Dean-Kawasaki solution. This block associates the scaling regime of Assumption FD3 to the validity of the moment bounds max where T 1 , . . . , T m ∈ [0, T ], and Θ was introduced in (10). The difference in the norms of the test functions stems from a difference in underlying mathematical arguments (depending on the circumstance, we will either use the maximum principle or the Sobolev embedding Theorem). The necessary tools for this point are contained in Subsection A.2.
Block 5. Inductive argument. Block 1 essentially states that computing cross-variations of discrete Dean-Kawasaki models yields linear functionals (24), as well as negligible corrections related to the negative part ρ − h . Taking Block 2 also into account, this leads to the following crucial observation.
The Itô correction term in the Itô differential of smooth enough nonlinear functions ψ applied to (24) and their expected values is a sum of: -negligible terms featuring ρ − h and the numerical error, and -yet another (possibly different) nonlinear functionψ applied to (24) and their expected values.
This property allows to set up both proofs using an induction argument whose inductive step is the change in nonlinear function (from ψ toψ): the residual terms (featuring ρ − h and the numerical error) are estimated at each step, and are not fed to the next step.

3.4.
The key step for the accuracy estimate for fluctuations in Theorem 2. For use in the next proposition, we define the two function spaces L q pow,r ,L q pow,r as Furthermore, we emphasise that we use the shorthand notations i.e., we implicitly multiply vectors in an element-wise fashion respectively evaluate vectorial functions by a vector of (time) parameters in an element-wise way.
Theorem 2 will be seen to be an easy consequence of the following crucial proposition and an inductive argument.
Proposition 7. Let µ N t denote the empirical measure of N independent Brownian particles as defined in (2).
Let ρ h be a solution to the Dean-Kawasaki equation discretised using finite differences on a uniform grid (21). Suppose furthermore that Assumption FD1 (details of operators ∆ h and ∇ h ), Assumption FD2 (initial condition on Brownian particle system), Assumption FD3 (scaling assumptions), and Assumption FD4 (positivity-preserving properties of mean-field limit) hold.
Let M , p ∈ N, q ∈ N, and r ∈ N 0 . Let ψ : Then there exist test functionsψ t kl ,φ t kl , ψ 0 , and φ 0 as well asT kl ∈ R M +1 such that Furthermore, Err num and Err neg are subject to the estimate where E(N, h) is defined in (11).
Under the additional assumption that ϕ k L 2 = 1 and´T d ϕ k dx = 0 for all k, that we have the additional bounds and as well as The proof is split into four steps. In Step 1, we provide deterministic estimates of suitable backwards diffusive equations of relevance, as well as basic stochastic estimates associated with the Dean-Kawasaki dynamics FD-DK.
Step 4 bounds the residual terms Err num , Err neg in (34b).
Proof of Proposition 7.
Step 1: definitions and elementary estimates. Let φ t m satisfy the backwards heat equation (14) subject to φ Tm m := ϕ i . Define the function ψ t : R M → R by setting ψ T := ψ and by evolving ψ t backward in time using the backward diffusion equation The purpose of the definitions of φ t m and ψ t will become clear in Step 2 and 3 below. Note that these definitions entail (where for simplicity we have assumed that the eigenvalues of Λ are nondegenerate; otherwise, we replace the formula by its natural analogue) with This implies and thus Arguing similarly, we deduce and therefore Using the estimate (49), under the additional assumptions on the ϕ k stated above we infer whenever T k > t. This in particular implies (35e). A similar argument yields Now fix η ∈ W 1+Θ . Let η h satisfy the discrete backwards heat equation (15) subject to η T h := I h η. We observe that the moment estimate holds for any j ∈ N. To see this, we use (21) and deduce that for any t > 0 Doob's martingale inequality, the moment bound (65), and the estimate sup (which is, depending on Θ, a consequence of either the discrete maximum principle or the Sobolev embedding theorem) yield (41). It is also straightforward to notice that Furthermore, we write where P · and P · h have been introduced in Subsection 3.1. Term T 1 is bounded using (73), while T 2 is settled using (18) from Assumption FD2. Altogether, this leads to Step 2: proof of (34a). Using Itô's formula and the fact that , plugging in the equation (14) satisfied by φ t , and taking the expected value, we obtain Integrating in t, recalling that φ T = ϕ, and plugging in the equation (36) satisfied by ψ t , we obtain We then defineψ t kl : Moreover, we setT kl := (T 1 , . . . , T M , min{T k , T l }). With these definitions, and in view of (which is a consequence of the maximum principle) and the definition ofφ t kl . Likewise, the estimate (35b) is immediate by the definition ofψ t kl , the estimate (38), and the definition of the norms · L q pow,r . Finally, from (39) and the definition ofψ t kl we deduce (35e).
Step 3: proof of (34b). Using Itô's formula and (21), we infer Using the fact that −∂ t φ h,k = χ t≤T k 1 2 ∆ h φ h,k and taking the expected value, we obtain Using the cross-variation identity (28), we get Switching to integral notation, using (36) as well as φ T h = I h ϕ, and adding zero, we obtain Adding zero once more and using the fact that ρ h (0) = E[ρ h (0)] (which is a consequence of Assumption FD2), we arrive at where we have set as well as and Err num,2 := 1 Using the definitions (47a) and (47b) and setting Err num := Err num,1 + Err num,2 , this yields the representation (34b).
Step 4: estimates for Err neg and Err num,i . We begin with Err neg ; it is easily seen to be bounded by This entails (35d). Furthermore, the analogue of (39) for the second derivative, and the time integrability of the singularity {min m:Tm≥t (T m − t)} −1/2 entail (35g).
Suppose that all ϕ m have vanishing average and are normalized in the sense ϕ m L 2 (T d ) = 1; suppose furthermore that whenever T m = Tm, the corresponding ϕ m and ϕm are orthogonal to each other in L 2 (T d ). Define Denoting the pseudo-inverse of the (possibly degenerate) nonnegative symmetric matrix Λ t defined in (37) by Λ −1 t , we have the estimate Proof. To simplify notation, we define T 0 := 1 2 T 1 . Estimating the matrix in (37), writing ϕ k (x) := n∈Z d a k,n exp(−in · x), and using the fact that −∂ t φ t k = 1 2 ∆φ t k , we get for Using the fact that φ T k k = ϕ k , that φ t k L 2 (T d ) ≤ 1 for all t, that the ϕ k have vanishing average, and our assumption on the orthogonality of the ϕ k with the same T k , we deduce Note that (Λ t ) kl = 0 whenever T k < t or T l < t. This concludes our proof.
3.5. Proof of Theorem 2. For finite difference discretization schemes, Theorem 2 is an easy consequence of Proposition 7.
Proof of Theorem 2 in the finite difference case. Taking the difference of (34b) and (34a) and using (35f) and (35g), we see that Proposition 7 implies The inequality (40) implieŝ In case j = 1, (51) entails the desired bound by the estimate onψ t kl upon replacing ψ in (50) by its convolution with a mollifier on the scale N −1/2 , which we denote by η N −1/2 . This is a straightforward result of the convolutional inequalities pow,0 . For j > 1, taking the difference of (34b) and (34a), using the bounds (35a), (35c), (35d), and iterating this estimate j − 1 times (i. e. using in each step again (34b) and (34a) to estimate the terms of the form only bounding the terms in ψt using (52) in the last step), we deduce We use estimate (53) in (51) to bound the terms of the form Therefore, estimate (53) turns into . This, together with the fact that m (1/2)T1 is controlled by ρ min , proves Theorem 2 in the case of finite difference discretisations.
3.6. Recursive step for Theorem 3. In Theorem 2, one is forced to distinguish between the different final times T 1 , . . . , T M due to the singular nature of the evolution equation for ψ (36). In contrast, ψ is static in Theorem 3: therefore, its proof can be detailed in the (notationally much more convenient) case of equal final times T 1 = · · · = T m = T without losing in generality.
We first introduce some handy notation. For t ≤ T , we abbreviate solves the backwards heat equation (14) (respectively, the backwards discrete heat equation (15)) with datum ϕ (respectively, I h ϕ) at time T . Given a multi-index j = (j 1 , . . . , j M ) such that |j| 1 = j ∈ N and a set of smooth test functions ϕ = (ϕ 1 , . . . , ϕ M ), we abbreviate and we set In order to show Theorem 3, we first provide a series of preliminary results.
Lemma 9 (First moments). The first moments of the Dean-Kawasaki model in Definition FD-DK agree with those of the Brownian particle system. Namely, for ϕ ∈ C 1 , we have Proof. This follows promptly from Lemma 15, as neither S N (I h ϕ, T, t) nor T N (ϕ, T, t) admits drift.
where P · h is the solution operator for the discrete backwards heat equation, see Subsection 3.1. On the other hand Lemma 15 also implies where P · is the solution operator for the backwards heat equation, see Subsection 3.1. We get A 1 = B 1 since the first (centred) moments agree (see Lemma 9). Furthermore, (66) and (19) grant The bounds (75) and (19) promptly givê We decompose A 2 − B 2 as follows where we have also used that ρ h (0) is deterministic. The term C 1 is bounded using (18) applied to the function η := P t (∇φ t 1 · ∇φ t 2 ), while the term C 2 is dealt with using (73) with choice ϕ := ∇φ t 1 · ∇φ t 2 . All together, we obtain the bound The proof is complete.
Proposition 11 (Recursive formula for higher moments). Let Θ be as in (10). Fix ϕ = (ϕ 1 , . . . , ϕ K ) ∈ C p+3+Θ M , a vector j = (j 1 , . . . , j M ) such that |j| 1 = j. For each pair i, j ∈ {1, . . . , M }, let j ij be as defined in Lemma 15. Let E(N, h) be as defined in (11). Assume the validity of Assumptions FD1, FD2, FD4, FD3. We recall the abbreviation for the difference of moments (see (54), (55)) Then we have the recursive formula Proof. We use Lemma 15 to deduce In analogy to the notation of Lemma 10, we define Let P · (respectively, P · h ) be the solution operator for the backwards heat equation (respectively, for the discrete backwards heat equation), see Subsection 3.1. We then proceed above as On the other hand It is straightforward to notice that A 1 − B 1 can be settled using the estimates for the moments of order j − 1, as (for each pair k, l) the exponent vector j is decreased by two units to j kl , while the additional test function ∇φ t k · ∇φ t l is picked up. The bound for A 3 relies on the Cauchy-Schwartz inequality, Corollary 17, (19), and (66). It reads The term A 4 may be bounded as followŝ The difference A 2 − B 2 is rewritten as where equality (60) is valid because the term is deterministic. The term T 1 is dealt with using the estimates of order j − 2 (as, for each k, l, the exponent vector is decreased by two units to j kl ). The term T 2 is settled with the same arguments as for term C 2 in (57), with the additional use of the Hölder inequality and of (92). We obtain Putting together all the estimates and integrating in time gives (59).
Remark 12. The finite-difference error in (59) accounts for two different errors: • the difference between the initial conditions ρ h,0 and the empirical density µ N 0 , as well as the difference between the solutions to continuous and discrete backwards heat equations. This is captured in the term A 2 − B 2 for the second order moment, and in the term T 2 for higher moments.
• the difference between I h (∇φ t k ·∇φ t l ) and ∇ h φ t k,h ·∇ h φ t l,h . As anticipated in Subsection 3.3, Block 3, the high-order accuracy of the difference between the solutions to continuous and discrete backwards heat equations relies on the discrete final datum to be the interpolant of the continuous final datum.
Step 1: Interpreting (59). The recursive relation (59) may be visualised in the following way: i) Each moment of order j produces residuals Err neg and Err num . ii) Each moment of order j is linked recursively to a collection of moments of order j − 1 (A j−1 recursion ) and a collection of moments of order j − 2 (A j−2 recursion ).
iii) The overall bound for D(j, ϕ, T ) is given by summing all the residuals for all moments found by exhausting the recursive relation. More specifically, it holds where R K is the sum of all residuals associated with the moments explored after exactly K steps. Therefore, we only need to suitably control R K for K = 0, . . . , j−2.
In order to do this, we need the following auxiliary bound.
Step 2: Auxiliary bound. At every step of the recursive relation, the sets of test functions which are fed into the lower order terms A j−1 recursion and A j−2 recursion are modifications of the current set of test functions, specifically: • in the case of A j−1 recursion , one instance for each of two functions ϕ k , ϕ l is replaced by the product ∇ϕ k · ∇ϕ l ; • in the case of A j−2 recursion , one instance for each of two functions ϕ k , ϕ l is removed from the set of test functions, and a pre-factor ϕ k C 1+Θ ϕ l C 1+Θ is gained.
It is thus natural to define the object where r is a given way of exhausting the recursive relation for K steps (i.e., a sequence of K moves dictating whether moments of type A j−1 recursion or A j−2 recursion are explored at each step), where ψ K,r is the set of test functions after K steps with sequence r, where j K,r is the corresponding set of powers, and where Y K,r is the overall pre-factor cumulated from all the moments of type A j−2 recursion for the sequence r. For each γ ∈ N 0 , we have the bound which is justified by the following observations: • The number of occurrences of the original functions ϕ (i.e., j) is preserved, regardless of the path r. This is straightforward to verify by direct inspection of how the recursive terms A j−1 recursion and A j−2 recursion handle the test functions. • The factor j 2K provides a bound on the product of the number of individual addends making up the functions {ψ K,r,m } m and of the number of individual addends making up the functions of type ψK ,r,m (whereK < K) found in the term Y K,r . This is a simple consequence of the fact that, whenever a step of type A j−1 recursion is performed, such product can be multiplied by at most K · K = K 2 (i.e., by the product of the maximum lengths of the addends making up the two functions φ k and φ l which give rise to the new test function ∇φ k · ∇φ l ). When a step of type A j−2 recursion is performed, such product does not increase.
• The factor M m=1 ϕ m jm C max{γ;1+Θ}+K takes into account the evaluation of the norms for all functions (both {ψ K,r,m } m and those associated with Y K,r ) by using the most restrictive exponent between 1 + Θ (needed in any step of type A j−2 recursion ) and γ (which is the exponent we are interested in), and adding K (to reflect the unitary increment of differentiation entailed by each step of type A j−1 recursion ). • The term j (max{γ;1+Θ}+1) is associated with the pre-factor of the inequality applied with ≤ j (j is the maximum number of factors in the addends of type i=1 f i making up any function ψ K,r,m and any function associated with Y K,r ), and with β = max{γ; 1 + Θ}. The overall pre-factor j j(max{γ;1+Θ}+1) results from multiplying j (max{γ;1+Θ}+1) by itself j times (j being an upper bound for the total number of functions ψ K,r,m together with all functions associated with Y K,r ).
Crucially, (61) only depends on K and j, and not on the specific path r.
Step 3: Bounding R K . The quantity 2 K j 4(K+1) = 2 K × j 4K × j 4 is a bound for both the number of residuals of type Err neg and Err num associated with the moments explored after exactly K steps: Such a quantity is the product of 2 K (accounting for the recursive splitting of (59) into two families of moments of lower order), of j 4K (accounting for a bound of the pre-factor M k,l=1 (j k − δ kl )j l /2 multypling each of the two families of moments), and of j 4 (accounting for a bound of the pre-factor M k,l=1 (j k − δ kl )j l /2 multypling the residual terms). Using (59) and (61), we obtain Step 4: Concluding the argument. Since D(j, ϕ, T ) ≤ , which -up to trivial rescaling in N 1/2 -is precisely (12). Proposition 13. Let the assumptions and notation of the finite difference case of Theorem 2 be in place; in particular, let ρ h be a solution to the Dean-Kawasaki equation discretised using finite elements in the sense of (21). Assuming in addition the scaling

Exponentially decaying estimate for E ρ −
for any B ≥ 1. In particular, we can deduce for any j ≥ 1, as well as Proof. We split the proof into several steps.
Step 1: energy estimates for test functions. In order to evaluate ρ h (x 0 , T ) at a given point x 0 , we choose φ h (·, T ) ∈ L 2 (G h,d ) as the function satisfying (φ h (·, T ), η h ) h = η h (x 0 ) for all η h ∈ L 2 (G h,d ) and evolve φ h in time by the backward heat equation By the standard energy estimate for the discrete heat equation we get Step 2: exponentially decaying bounds for |ρ h − E[ρ h ]|(x 0 ) for chosen point x 0 . Using (96), (97), and (67), we obtain by the Itô formula for any positive integer j In particular, Integrating in time up to a stopping time T s and taking the expected value, we obtain Choosing T s for arbitrary but fixed B ≥ 1 as we get using ρ max ≥ ρ min and the assumption |E Using Doob's martingale inequality, we deduce for nonnegative even integers j E sup Raising both sides to the power j/2 and using Chebyshev's inequality, we get after optimizing in j P sup In particular, we deduce by the definition of φ h (·, T ) Step 3: extending the estimate to finitely many time points in [0, T ∧ T s ]. Applying the previous estimate for all x 0 ∈ G h,d (there are ∝ h −d of such points), and for all times h β , 2h β , 3h β , . . ., for some β > 0 to be chosen, we obtain Step 4: extending the estimate to all times in [0, T ]. It only remains to pass from the discrete times ih β to all times t and to remove the restriction to times t ≤ T S . Let e k ∈ L 2 (G h,d ) be nodal function satisfying e k (x j ) = δ kj . Then the differential entails, using in a second step also Doob's maximal inequality and abbreviating W(ρ + h , e k ) : Using the triangle inequality for the first term on the right-hand side and a (straightforward but rather pessimistic) estimate on the quadratic variation of W, we obtain By absorption, the triangle inequality, the fact that Using Young's inequality and absorbing as well as using the fact that for ih β ≤ T S we have |ρ h | ≤ (B + 1)ρ max , we obtain For β ≥ 6d + 8 and for all h ≤ c(ρ min , ρ max ), we obtain Step 5: obtaining (63). Overall, from (69) and (70) we conclude P sup
Step 7: obtaining (66). We use the Hölder inequality and the lower bound E[ρ h ] ≥ ρ min and obtain (66) via the estimate

Numerical examples
In this section, we give numerical examples that illustrate that the Dean-Kawasaki equation correctly captures the fluctuations of diffusing non-interacting particles 1 . We limit our attention to the case d = 1.
To compute the motion of N Brownian particles, we perform a direct simulation based on the transition probabilities; this is feasible as our numerical experiments only concern empirical measures µ N t at two different times T 1 and T 2 (see below). Our discretisation of the Dean-Kawasaki equation is obtained as follows: • For the spatial discretisation of the Dean-Kawasaki equation (1), we use the finite difference scheme from Definition FD-DK with order p = 1. • To discretise the spatially semi-discrete equation in time, we use the (two-step) BDF2 scheme (see, e.g., [22]). The first timestep is performed using an explicit treatment for the noise and a mixed implicit-explicit Euler scheme for the deterministic diffusion. 1 The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.
• Overall, our discrete scheme for the Dean-Kawasaki equation (1) reads for the first timestep and for the (m + 1)-th timestep, m ≥ 1, where (β y ) are independent Brownian motions. • We place the initial positions {w k (0)} N k=1 of the Brownian particles only at grid points of G h,1 . Consequently, we define the initial condition ρ h (0) by requiring that the equality µ N 0 , η = (ρ h (0), I h η) h holds for any test function η. This way, we avoid any error caused by deviating initial conditions. • As we are primarily interested in scaling in h and N , we make the following choices: we set the time-step ∆t := 0.001, which, according to our numerical convergence tests, is small enough for the spatial discretisation error to dominate over the time error, and we keep the discretisation parameter h above or equal to the threshold 2π · 2 −7 ≈ 0.05, so that the finite difference error dominates over the error associated with the negative part of ρ h .
Using a Monte-Carlo approach with M 1 realizations, we next computed the centered stochastic moments for test functions ϕ 1 , ϕ 2 , times T 1 , T 2 , and integer exponents j 1 , j 2 specified below. We then compared these stochastic moments to the corresponding centered stochastic moments of the empirical density µ N M Brownian j1,j2 the latter being also computed by a Monte Carlo approximation with M realizations. We have performed various simulations in order to assess the convergence of the moments with respect to h, N , and to compare the discretisations to the linearised Dean- i. e. our relative error is basically independent of the particle number N and only depends on the grid size h. . For the same choice of initial data ρ 0 (x), test functions ϕ i (x), and times T i , and two different choices of N , we investigate the difference of performance between the time-discretised Dean-Kawasaki is the natural counterpart to M DK j1,j2 , in Figure 6 as a function of the discretisation parameter h. . We observe that the two models show the same behaviour for the second moment associated with the exponents (j 1 , j 2 ) = (2, 0). This is expected, as both models share the same quadratic variation structure of the noise (more explicitly, one can readapt Lemma 10 to the linearised case). On the contrary, the nonlinear model visibly outperforms the linearised model for the higher moment associated with (j 1 , j 2 ) = (2, 1). The reason for this is that one can not readapt Proposition 11 to the linearised case, as doing so would result in lower order moments comprising both the Dean-Kawasaki solution and its mean-field limit, thus breaking the very recursive structure of the Proposition.

Comparison with linearised
We have chosen a relatively low number of particles N a particular couple of test functions (with ϕ 2 approximately matching the quadratic variation associated with the second test function after some time, i.e., ϕ 2 ≈ ∇|ϕ 1 (T /4)| 2 , thus giving non-trivial correlation between ϕ 1 and ϕ 2 ) in order to make the difference between the two models more pronounced. Such difference is not completely clear cut though, as one can see for the lowest values of h in the bottom figure. This behaviour is likely caused by: • the reduced accuracy of the BDF2 integration method for low h; • in the case of Figure 6 (Bottom), an accuracy saturation.   N = 4096 (Bottom)). We observe that the discretised Dean-Kawasaki model outperforms -to a good extentthe linearised version for the moment associated with (j 1 , j 2 ) = (2, 1).
Proof. We recall the relation N L = 2π/h and definition It is easy to use the continuous and discrete backwards heat equations (14)- (15) to deduce that the Fourier coefficients of φ t i and φ t i,h , i ∈ {1, 2}, arê for some functional P (h, ξ). As the discrete Laplacian ∆ h is a (p + 1)-th order approximation of the true Laplacian with order p + 1, it is easy to see that Furthermore, since ∆ h is a symmetric finite difference operator, it is easy to see that P (h, d) is nonnegative. This fact, together with the convexity of the exponential function (which in turn implies the monotonicity of the ratio (e x − e y )/(x − y) in either one of the two variables, provided the other one is kept fixed) gives and therefore We can deduce that the discrete Fourier expansion of I h φ t i from (76) is where we have also used the fact that Z d = I h + LZ d . We deduce Since I h ϕ i = φ T i,h , we carry on in and write We estimate where the final step is justified bŷ which is valid since 2(p + 1) > d, as d ≤ 3 and p ≥ 1. We continue in (81) as and (73) is proved.
We adapt the arguments carried out so far and write After simple algebraic rearrangements in (83), we get The term T 1 is estimated using (79), giving T 1 ≤ Ch 2(p+1) ϕ i 2 C p+2 . The term T 2 is estimated by relying on (16), giving T 2 ≤ Ch 2(p+1) ϕ i 2 C p+1 . As for T 3 , we rely on the fact that ∇ϕ i (ξ + Lz) = −iξφ i (ξ + Lz) and write and (74) is proved. To prove (75), we write Now let φ be the solution to (15) with final datum φ T = ϕ ∈ C 1+Θ . If (15) admits a discrete maximum principle, then If (15) We focus on T 4 . It is easy to see that The estimate for T 5 is even more straightforward, and it reads and (75) is proved.
A.2. Stretched exponential moment bounds for the Dean-Kawasaki solution and the particle system. We compute the Itô differential of the quantities in (54). (F ρ (t)e d h,y, , ∇ h φ t i,h ) h dβ y, ∇φ t m (w k (t)) · dw k (t) ∇φ t k (w r (t)) · ∇φ t l (w r (t)) dt. (86b) Proof. All differentials in this proof are with respect to the variable t. We prove (86a) in three steps.
Lemma 16. Let φ solve the heat equation (14), with final datum φ T = ϕ and t ≤ T . Let ρ h be as given in Definition FD-DK. For any 2 ≤ j ∈ N, we have (F ρ (t)e d h,y, , ∇ h φ t h )dβ (y, ) Taking the expected value, we obtain .
We prove (89) by induction on j. The case j = 2 is easily settled. Now take j > 2 and assume the validity of (89) for j − 1. We use (91)  where N −1/2 W(ρ + h , e i ) is a real-valued martingale with quadratic variation given by Remark 18. Unlike in the finite-difference case, we can only provide an explicit representation of the Dean-Kawasaki noise in the case p = 1. This is due to the fact noise is nonlinear, and only preserves piece-wise constant functions (these being gradients of test functions in X p h , p = 1). We are not aware of any finite-dimensional representation of the martingale term N −1/2 W(ρ + h , e i ) in (96) in the case p > 1.
We now present the finite element counterparts of Theorems 2-3.   where E(N, h) has been defined in (11).