From Particle Fluctuations to Macroscopic Evolution Equations: The Case of Exclusion Dynamics ? ?? Shenglin Huang1[0000−0002−7450−6966] , Santiago Aguado-Montero2 , Nicolas Dirr3[0000−0003−3634−7367] , Johannes Zimmer4[0000−0002−7606−2422] , and Celia Reina1[0000−0001−9302−3401] 1 Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia PA 19104, USA {slhuang, creina}@seas.upenn.edu 2 Escuela Técnica Superior de Ingenierı́a, Universidad de Sevilla, Sevilla 41092, Spain saguado@us.es 3 School of Mathematics, Cardiff University, Cardiff CF24 4AG, UK DirrNP@cardiff.ac.uk 4 Fakultät für Mathematik, Technische Universität München, 85748 Garching, Germany jz@ma.tum.de Abstract. We consider the symmetric simple exclusion process in one dimension as an example of a stochastic particle process exhibiting anoma- lous diffusion; this process is so slow that the mean-square displacement of a tagged particle is sub-linear. This implies that the standard mean- square displacement method (Einstein relation) cannot be applied to determine the diffusivity for the associated macroscopic equation de- scribing the hydrodynamic limit. We demonstrate that a recent approach developed by the authors and collaborators, based on the covariation of the fluctuations, is applicable to this process and does not only deliver the transport coefficient, but the entire evolution operator associated with the formulation of the macroscopic equation as an entropic gradi- ent flow. Furthermore, the approach relies on fluctuations of the density field (macroscopic fluctuations) as opposed to particle level data. This data could in principle be obtained from experimental observations. Keywords: Anomalous diffusion · Fluctuation-dissipation relation · Fluc- tuating hydrodynamics. 1 Introduction Anomalous diffusion is characterised by a nonlinear dependence of the mean squared displacement (MSD) as a function of time, in contrast with the classical ? Supported by the University Research Foundation Grant from the University of Pennsylvania (SAM and CR) and a Royal Society Wolfson Research Merit Award (JZ). ?? Copyright c 2020 for this paper by its authors. Use permitted under Creative Com- mons License Attribution 4.0 International (CC BY 4.0). Particle fluctuations to macroscopic evolution for exclusion dynamics 141 diffusion of Brownian particles, where MSD ∝ t (at large times). Examples of anomalous diffusion include sublinear diffusion observed in biological applica- tions, such as diffusion in the cell nuclei and membranes; other examples occur for instance in disordered media, astronomy, fluid mechanics and network the- ory [10]. A classic paper giving a broad overview on this topic is [5]. Anomalous diffusion poses numerous theoretical difficulties (see, for exam- ple, [5]). However, the challenges are also of very practical nature. Simulation of particle processes can quickly become prohibitively expensive. One thus seeks macroscopic descriptions in form of partial differential equations (PDEs). A com- mon procedure is to assume phenomenologically the structure of the govern- ing equation and then determine the transport coefficients appearing in it from particle simulations. Prominent methods to determine diffusion coefficients are Green-Kubo relations (based on time-correlation functions) and Einstein rela- tion (based on the mean square displacement) [7]. However, these methods do not provide information on the structure of the PDE. Furthermore, the Einstein relation has as underlying assumption that the MSD scales linearly with time, MSD ∝ t, and is thus not directly applicable to anomalous diffusion processes. In this contribution, we show that an alternative approach developed by some of the authors and their collaborators [6, 9] delivers the full structure of the PDE (and parameters therein) for an example of anomalous diffusion, namely for the symmetric simple exclusion process. More specifically, this approach delivers the operator of the macroscopic evolution equation (in discretised form), written as a gradient flow of the entropy functional, without making further phenomeno- logical assumptions on the macroscopic evolution. We now explain the symmetric simple exclusion process (SSEP) in more detail. In this process, particles sit on a lattice and jump stochastically to one of its neighbouring sites with a constant rate, which can be assumed to be one without loss of generality. The central feature of the process is that such a jump is only possible if the destination site is empty, which results in every site being occupied at most by one particle. If the target site is already occupied, then the particle cannot jump and remains at its current location. Mathematically, the jump rate from a site X to a neighbouring site X̃ can be expressed as gX→X̃ (η) := 21d η(X)(1 − η(X̃)), where d is the dimension of the space (here, d = 1), and η(X) is zero (one) if site X is empty (occupied); see [8, Section 2.2] for further details. In one space dimension, the resulting process is slow – as in a single-lane motorway, one slow particle is enough to impede the passage of all the other particles behind, effectively blocking them. Thus (albeit the hydrodynamic limit is a linear diffusion equation), the individual particles themselves are not following a Brownian motion on a microscopical level [2]. They are slowed down by the exclusion condition, which results in a sublinear power law for the mean square displacement, namely MSD ∝ t1/2 . The methodology of [6] has been shown to be applicable in this context, though, to determine the governing transport coefficient only. In this contribu- tion, we show that the method can be extended using the approach of [9] to 142 Sh. Huang et al. determine the full macroscopic evolution numerically from particle simulations, assuming the governing entropy is known. The key requirement for the application of the strategy outlined in [9] is that the macroscopic evolution is a gradient flow, and that fluctuations around such limit are Gaussian. This means that the density evolution has the following thermodynamic structure in the limit of infinitely many particles δS δS ∂t ρ = K(ρ) (ρ) =: Kρ (ρ), (1) δρ δρ where Kρ is a symmetric operator, and δS δρ is the variational derivative of the entropy functional, while for large but finite number of particles N , the evolution shall be formally described by the stochastic partial differential equation of the form r δS(ρ) Cp ∂t ρ = K(ρ) + 2K(ρ)Ẇt,x , (2) δρ N with Ẇt,x a space-time white noise, and C a constant such that  = C/N is the individual lattice site volume. Equation (2) highlights a very important relation between fluctuations and dissipation: the same operator K appears in the deter- ministic (macroscopic) part and in the fluctuations. This important connection precisely lies at the core of the approach in [9]. For the symmetric simple exclusion process, the macroscopic (hydrodynamic) limit under parabolic scaling is known to be the linear diffusion equation ∂t ρ = ∆ρ [8, Section 2.2]. Yet, the linear nature disguises nonlinear features that are apparent when written in the thermodynamic form, given in Eq. (1). Indeed, the macroscopic equation is a gradient flow of the entropy (further details in Section 2), with S being the mixing entropy, i.e., Z S(ρ) = − [ρ log ρ + (1 − ρ) log(1 − ρ)] dx (3) in dimensionless form, and K the operator K(ρ)ξ := − div(ρ(1 − ρ)∇ξ), i.e., δS δS ∂t ρ = ∆ρ = − div(ρ(1 − ρ)∇ (ρ)) = Kρ (ρ). (4) δρ δρ 2 Gradient Flow Structure of the SSEP In this section, we explain why the formulation (4) is a gradient flow. We include this calculation as the setting does not seem to be described in the literature, though readers with a focus on the computational aspects can safely skip this section. We show that (4) is the steepest descent in a weighted version of the standard L2 Wasserstein geometry. The argument is sketched adapting ideas by Benamou and Brenier [3]; see also [1]. In this context, a gradient flow of a functional S is defined as the evolution Z δS (∂t ρ, s2 )K−1 = s2 dx, (5) δρ Particle fluctuations to macroscopic evolution for exclusion dynamics 143 where the geometry of the gradient evolution is defined by an inner product, here symbolically denoted K−1 ; Eq. (5) is a weak formulation in the sense that it has to be satisifed for all suitable test functions s2 . We now define the inner product in (5). For the SSEP, the evolution takes place on the space of probability measures with bounded Lebesgue density 0 < ρ < 1 (almost surely). For such measures, we define formally an inner product on the tangent space, Z ρ (s1 , s2 )K−1 := ∇Ψ1 ∇Ψ2 dx, (6) 1−ρ where sj = − div(ρ∇Ψj ) for j = 1, 2. (The tangent space can be interpreted as time derivatives satisfying the continuity equation.) The last equation also defines the class of test functions in Eq. (5). We now compute the expressions of the gradient flow evolution for the SSEP. On one hand, we find with undetermined Ψ1 and Ψ2 and ∂t ρ = − div(ρ∇Ψ1 ) and s2 = − div(ρ∇Ψ2 ) (7) for the expression on the left of (5) Z 1 (∂t ρ, s2 )K−1 = ρ∇Ψ1 ∇Ψ2 dx. (8) 1−ρ On the other hand, the expression of the right of (5) becomes with the mixing entropy (3) Z Z   δS ρ s2 dx = − log div(ρ∇Ψ2 ) dx δρ 1−ρ Z   Z ρ 1 = ρ∇ log ∇Ψ2 dx = ∇ρ ∇Ψ2 dx. (9) 1−ρ 1−ρ As, in the absence of topological obstructions, any gradient vector field ∇f can be written in the form 1 ∇f = ∇Ψ2 (10) 1−ρ   1 (by solving div 1−ρ ∇Ψ2 = ∆f ), the combination of (8) and (9) yields Z Z ρ∇Ψ1 ∇f dx = ∇ρ∇f dx. With the first continuity equation in (8) and integration by parts, this becomes Z Z ∂t ρf dx = ∇ρ∇f dx, which is the weak form of the hydrodynamic limit (4) associated with the SSEP. 144 Sh. Huang et al. For the SSEP, we call Eq. (6) the thermodynamic metric. It is a weighted version of the classic Wasserstein metric. We notice that there are different ways to obtain the linear diffusion equa- tion (4), ∂t ρ = ∆ρ as gradient flow. For example, we can use theR standard Wasserstein metric and the standard (Boltzmann) entropy S = − ρ log ρ dx (see, for example, [1]), or the thermodynamic metric (6) of the SSEP and the mixing entropy (3). This shows an advantage of the entropic gradient flow for- mulation: here the differences between different proceses become visible, even if they have the same macroscopic limit. For example, Brownian particles define the classic (Boltzmann) entropy gradient flow, and the SSEP defines the thermo- dynamic flow of the mixing entropy. Even if the hydrodynamic limit of the two processes is the same, their fluctuations are different, as evident from Eq. (2). The entropy-metric pair also records these differences. 3 Computational Framework 3.1 Discretisation of the Dissipative Operator The partial differential equation (4) and the operator Kρ that the PDE en- codes are infinite dimensional objects. Thus, their full characterisation from particle level data either requires infinite amount of data, which is of course unattainable from numerical simulations, or, alternatively, optimal fits from a finite library of infinite dimensional operators. To avoid the phenomenology of the latter approach, we here use particle data to physically infer a discretised version of the dissipative operator (i.e., now a finite dimensional object), as pre- viously proposed by the authors in [9]. Specifically, we consider a finite element approximation scheme of the density field ρ and the thermodynamic driving force F := δS/δρ of the form, X X ρ(x, t) ≈ ρa (t)γa (x), and F (x, t) ≈ Fa (t)γa (x), (11) a a where {γa (x)} is a basisPof functions with local support. ThesePsatisfy the so- called partition of unity a γa (x) = 1, linear field reproduction a γa (x)xa = x and Kronecker-Delta property γa (xb ) = δab (the latter is particularly convenient when numerically solving the PDE with Dirichlet boundary conditions). For the case of SSEP, linear finite element shape functions provide sufficient regularity to approximate the above fields and are here chosen for their simplicity. For a regular mesh with nodal spacing ∆x in one dimension, such shape functions can be written as 1 γa (x) = max (|x − xa | , 0) . (12) ∆x Substituting the approximation given in Eq. (11) into the evolution equation Eq. (4), and further testing it with a shape function γb , one obtains the weak form of the evolution equations, X ∂ρa X hγa , γb i = hKρ γa , γb i Fa , for all b, (13) a ∂t a Particle fluctuations to macroscopic evolution for exclusion dynamics 145 where the bracket h·, ·i denotes the integral inner product in the L2 space. The resulting matrix hKρ γa , γb i represents the discretised dissipative operator, which we will seek to estimate from particle data, with the strategy outlined in Sec- tion 3.2. We note that for processes satisfying conservation of mass, as the SSEP of interest here, the discretised operator entries should satisfy the constraint X hKρ γa , γb i = 0. (14) b P Equivalently, this constraint can be written as b hKρ γb , γa i = 0 for symmetric operators, as the one associated to SSEP or any other purely dissipative processes with the structure of Eq. (1). Despite the finite-dimensional nature of hKρ γa , γb i, this operator is, a priori, a function of the entire density profile, or, for the already assumed finite ele- ment approximation, a function of all the nodal density values. Consequently, its domain is extremely vast and dependent on the specific simulation domain, rendering the sought-after approach computationally intractable and of limited generality. Yet, many physical processes, such as SSEP are governed by local operators. This property, combined with the fact that the chosen basis functions have local support, implies that hKρ γa , γb i is only non-zero for the cases where node a is equal or near to node b (concretely, a = {b − 1, b, b + 1}), and that only the local density profile is needed to determine the discretised operator. This important observation allows us to approximate the density profile by means of a Taylor expansion, and rewrite Eq. (13) as X ∂ρa X D E hγa , γb i ≈ K(ρa +∇ ρ| (x−xa )+...) γa , γb Fa . (15) a ∂t a a∈{b−1,b,b+1} For SSEP, we will truncate the Taylor expansion after the liner term, rendering hKρ γa , γb i a function of only two variables: ρa and ∇ ρ|a . From a practical per- spective, this implies that only linear profiles, encompassing different densities and gradients, are to be simulated to characterise the discretised operator from the associated particle data. Once the discrete operator has been estimated, this may then be used to simulate the evolution of arbitrary initial profiles, i.e., not necessarily linear ones. For this, we will use a forward Euler time discretisation scheme, and rewrite Eq. (15) as ρi+1 − ρi M = Ki Fi . (16) ∆t where the superscript ”i” denotes the time step ti , and all quantities have been written in the form of vectors and matrices. Specifically, ρi = ρia N ×1 is the γ vector of densities at the nodal positions, ∇ρi = ∇ρ|ia Nγ ×1 is the vector of    density gradients, and Fi = Fi ρi = Fai ρi N ×1 is the vector of thermody- γ namic driving forces, where Nγ is the number of basis functions. Furthermore, the matrix forms of M (independent of time) and the dissipative operator Ki read 146 Sh. Huang et al. 2/3 1/6 0 · · · 1/6   .   1/6 2/3 1/6 . . . ..     M = (hγa , γb i)Nγ ×Nγ = ∆x  0 1/6 . . . . . . 0  ,   (17)    . .  .. . . . . . . . . 1/6   1/6 · · · 0 1/6 2/3 Nγ ×Nγ Ki = Ki ρi , ∇ρi = Kba i   Nγ ×Nγ  i i i  K11 K12 0 ··· K1N γ   Ki Ki Ki .. ..   21 22 23 . .   =  i .. ..   0 K32 . . 0    . .. .. ..   . . . . i  . KN  γ −1,Nγ  i i i KN γ 1 · · · 0 K N γ ,N γ −1 K N N γ γ (18) Nγ ×Nγ  0 i p i  K (X1 ) K (X2 ) 0 ··· K m (XN i γ )   K m (X i ) K 0 (X i ) K p (X i ) .. ..   1 2 3 . .   =  . . . .  , 0 K m (X2i ) . . 0     .. .. .. .. p i    . . . . K (XNγ )   p i m i 0 i K (X1 ) ··· 0 K (XNγ −1 ) K (XNγ ) Nγ ×Nγ where K 0 , K p and K m represent the three operator entries hKρ γa , γb i = Kba with a = b − 1, b, b + 1, respectively, Xai = (ρia , ∇ρ|ia ), and the periodic boundary conditions have been considered. 3.2 Fluctuation-Dissipation Statement to Compute the Discretised Operator from Density Fluctuations For any given density profile, hKρ γa , γb i can be computed from the covariation of the rescaled local density fluctuation as [9] 1 hKρ γa , γb i = lim E [(Yγa (t0 + h) − Yγa (t0 )) · (Yγb (t0 + h) − Yγb (t0 ))] , (19) h&0 2h where t0 is an arbitrary initial time that does not affect the result for the op- erator entries as long as the system has reached a local equilibrium, E [·] repre- sents the expectation, and√ Yγ are local rescaled density fluctuations, defined as Yγ = lim→0 hρ − ρ, γi /  with ρ = E [ρ ]. In practice, the particle system is of course simulated with a finite number of sites, and thus the limit as  → 0 in Yγ is numerically approximated. Similarly, the right hand side of Eq. (19) is Particle fluctuations to macroscopic evolution for exclusion dynamics 147 computed for a finite time step h, that is sufficiently small to capture the limit numerically, yet, sufficiently large for the particle system to exhibit some particle jumps. Finally, the expectation in this same equation will be approximated as the ensemble average over many but finite realisations of macroscopically identical density profiles. 3.3 Polynomial Regression of the Discretised Operator Equation (19) will enable the calculation of the discrete operator D E Kba = K(ρa +∇ ρ| (x−xa )) γa , γb a for a finite set of linear profiles, resulting in a discrete space Vdiscr of ρa and ∇ ρ|a . Yet, its use in Eq. (15) to simulate the evolution of arbitrary density profiles will require evaluating Kba at different density and density gradient values, for which a suitable interpolation scheme is therefore needed. Towards this goal, we recall that the three operator entries Kba with a ∈ {b − 1, b, b + 1} should satisfy the constraint given by Eq. (14) in order for the evolution to be mass preserving. Although such identity may not be exactly satisfied for the particle-inferred operator in Vdiscr , we here perform a polynomial regression on such data that, by construction, exactly satisfies the conservation constraint. More specifically, denoting by φj (X) (j = 1, 2, . . . , n) the fitting basis, with X = (ρ, ∇ρ), we approximate the three non-zero entries of the discretised operator as n X K0 = a0j φj (X), (20) j=1 Xn Km = am j φj (X), (21) j=1 n X p 0 m a0j + am  K = −K − K =− j φj (X), (22) j=1 where K 0 , K m and K p have been defined after Eq. (18). For simplicity, we use polynomials up to order-k as the basis functions, i.e., φj (X) ∈ {1, ρ, ∇ρ, ρ2 , ρ∇ρ, (∇ρ)2 , ..., ρk , ..., (∇ρ)k }. Specifically, we will use k = 4 to fit the discrete operator data for SSEP. Towards the goal of finding the optimal fitting parameter a0j and am j , we denote the data points from the numerical simulation as {Xi0 , Ki0 } with i = 1, ..., N0 for the diagonal entry Kbb , {Xim , Kim } with i = 1, ..., Nm for the sub- diagonal entry Kb,b−1 , and {Xip , Kip } with i = 1, ..., Np for the super-diagonal entry Kb,b+1 , where N0 , Nm , Np are the number of data points for the three entries, respectively. We then define the loss function L as the sum of the least 148 Sh. Huang et al. square error for all three entries, i.e.,  2  2 N0 n Nm X n 1 X X 1 X L=  φj (Xi0 )a0j − Ki0  +  φj (Xim )am m j − Ki 2 i=1 j=1 2 i=1 j=1  2 Np n 1 X X φj (Xip ) a0j + am − Kip  . (23)  + − j 2 i=1 j=1 The minimum of the loss function can be analytically computed by setting its partial derivatives with respect to a0 and am to zero, i.e., ∂L = Φ0T Φ0 a0 − K0 + ΦpT Φp a0 + am + Kp = 0,     (24) ∂a0 ∂L = ΦmT (Φm am − Km ) + ΦpT Φp a0 + am + Kp = 0.    m (25) ∂a We remark that the above equations use the matrix/vector forms for the fit- ting parameters a0 = (a0j )n×1 , am = (am 0 j )n×1 , the operator entries K = p 0 m m p (Kj )N0 ×1 , K = (Kj )Nm ×1 , K = (Kj )Np ×1 , and the basis functions Φ0 = (φj (Xi0 ))N0 ×n , Φm = (φj (Xim ))Nm ×n , Φp = (φj (Xip ))Np ×n . Denoting  0T 0 Φ Φ + ΦpT Φp ΦpT Φp  Φ= , ΦpT Φp ΦmT Φm + ΦpT Φp  0   0T 0 (26) Φ K − ΦpT Kp  a a= , and b = , am ΦmT Km − ΦpT Kp Equations (24)–(25) can be equivalently written in compact form as Φa = b, (27) delivering a simple linear system of equations from which the fitting parameters a can be obtained. Alternatively, a constrained least-square solver available in some commercial packages could be used to find the optimal fitting parameters. 4 Numerical Results 4.1 Discretised Operator from Particle Fluctuations In order to compute the operator entries numerically from density fluctuation using Eq. (19), we use the lattice kinetic Monte-Carlo (KMC) method [4], imple- mented in C++, to simulate the particle jumps for SSEP over a unit interval with Nsites = 5000 sites (i.e.,  = 1/5000). We perform simulations with 43 different initial linear (or piecewise linear) density profiles, which results in the discrete space Vdiscr shown in Fig. 1, with data points within the range of ρ ∈ [0.05, 0.95] and ∇ρ ∈ [−3, 3]. Particle fluctuations to macroscopic evolution for exclusion dynamics 149 Fig. 1. The probed discrete data space Vdiscr of (ρ, ∇ρ) used to evaluate the discretised operator. For each initial profile, we perform R = 105 realisations using the strategy of [6] aimed at partially alleviating the computational burden associated to the equilibration time (this strategy requires additional parameters, which are here chosen as R1 = 50, R2 = 2000, equilibration time tprep − tini = 4 × 10−6 and randomisation time t0 −tprep = 4×10−9 ). The actual time interval over which the covariation in density fluctuations are computed is h = 4 × 10−10 . Additionally, Nγ = 50 shape functions are used to discretise the macroscopic domain. The results for the three operator entries hKρ γa , γb i with a ∈ {b − 1, b, b + 1} are shown in Figs. 2(a)–(c) (blue circles) together with the analytical surfaces, obtained by performing a Taylor expansion on xa as mb 1 ∂2m hKρ γb , γb i = hm∇γb , ∇γb i = 2 + ∆x + O(∆x2 ), (28) ∆x 3 ∂x2 b mb−1 1 ∂m 1 ∂2m hKρ γb−1 , γb i = hm∇γb−1 , ∇γb i = − − − ∆x ∆x 2 ∂x b−1 6 ∂x2 b−1 + O(∆x2 ), (29) mb+1 1 ∂m 1 ∂2m hKρ γb+1 , γb i = hm∇γb+1 , ∇γb i = − + − ∆x ∆x 2 ∂x b+1 6 ∂x2 b+1 + O(∆x2 ), (30) with ∂m ∂2m 2 ma = ρa (1 − ρa ), = ∇ρ|a (1 − 2ρa ) , and = −2 (∇ρ|a ) + O(∇ρ2 ). ∂x a ∂x2 a (31) The diagonal entry Kbb is symmetric with respect to ρ = 0.5 and is indepen- dent with ∇ρ, while the two off-diagonal entries are non-symmetric and have 150 Sh. Huang et al. Fig. D 2. (a–c) NumericalE results (blue circles) for the discretised operator entries K(ρa +∇ ρ| (x−xa )) γa , γb with a ∈ {b − 1, b, b + 1} as a function of density ρ and a density gradient ∇ρ (at node a) in Vdiscr , shown in Fig. 1. The analytical predictions (smooth surfaces) based on Eqs. 28–30 are jointly shown. (d–f) Corresponding relative errors (in %), denoted as errKba with a ∈ {b − 1, b, b + 1}, for the numerical results, between the data and analytical results (blue circles) together with the relative errors for the polynomial regression of data (smooth surface with non-zero relative error). an opposite dependence due to the mass conservation constraint. The relative error between the data points and the analytical predictions is quantified in Figs. 2(d)–(f), where the relative error of the fit, following Section 3.3, is also shown. As it may be observed in the figure, the fourth order polynomial regres- sion with mass conservation strictly enforced significantly decreases the error: while the standard deviation of the relative errors of the original data points for three entries Kbb , Kb,b+1 and Kb,b−1 are 1.4%, 2.3% and 2.5%, respectively, these get reduced to 0.27%, 0.44% and 0.31%, respectively, for the polynomial fit. Similarly, the maximum relative errors decrease from 7.9%, 12.7% and 21.6%, to 2.6%, 5.7% and 3.8%. We remark that the relative errors are larger near ρ = 0 and ρ = 1, where the analytical values of operator entries are close to 0, and the relative errors thus become singular. Overall, the numerical strategy outlined in this section delivers a high accuracy for the discretised dissipative operator governing the simple exclusion process. 4.2 Macroscopic Simulation We now test the capability of the particle-inferred operator to predict the density evolution for an arbitrary initial profile. Concretely, we consider an initial density ρ(x, 0) = 0.5 − 0.3 cos(4πx) over the unit interval x ∈ [0, 1], discretised with Nγ = 50 shape functions (i.e., ∆x = 1/Nγ = 0.02). This profile is evolved Particle fluctuations to macroscopic evolution for exclusion dynamics 151 using the forward Euler scheme described in Eq. (16) with ∆t = ∆x2 /1000, the discrete operator obtained in Section 4.1 and the analytical driving force, given by the entropy in Eq. (3). Periodic boundary conditions are considered, and the system is evolved till almost reaching an equilibrium state with a flat density profile at the average initial density. Snapshots of the evolution are represented in Fig. 3, together with the an- alytic PDE (with identical spatio-temporal discretisation, and discretised op- erator given by Eqs. (28)–(31)), and the average over multiple realisations of long-time KMC simulations. Due to the large computational cost of these KMC simulations, these are performed using Nsites = 1000 sites and averaging over R = 200 realisations (we recall that the discretised operator was inferred from particle simulations using Nsites = 5000). The numerical results depict an excel- lent agreement between the three curves, making them almost indistinguishable in the figure. Fig. 3. Snapshots of the density evolution for the symmetric simple exclusion process at times t = 0 (a), 0.002 (b), 0.004 (c), 0.006 (d), 0.01(e), and 0.02 (f). Such evolution is computed with three different strategies: PDE with analytic operator, PDE with particle-inferred operator, and long-time KMC simulations. The present contribution therefore shows that the computational strategy of [9] is able to extract the dissipative operator for the symmetric simple ex- clusion process, an anomalous diffusion process, and hence, it can be used to predict the non-equilibrium macroscopic evolution of arbitrary initial profiles. References 1. Adams, S., Dirr, N., Peletier, M., Zimmer, J.: Large deviations and gradient flows. Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 371(2005), 20120341 152 Sh. Huang et al. (2013). https://doi.org/10.1098/rsta.2012.0341 2. Arratia, R.: The motion of a tagged particle in the simple sym- metric exclusion system on Z. Ann. Probab. 11(2), 362–373 (1983). https://doi.org/10.1214/aop/1176993602 3. Benamou, J.D., Brenier, Y.: A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numer. Math. 84(3), 375–393 (2000). https://doi.org/10.1007/s002110050002 4. Bortz, A.B., Kalos, M.H., Lebowitz, J.L.: A new algorithm for Monte Carlo simulation of Ising spin systems. J. Computational Phys. 17(1), 10–18 (1975). https://doi.org/10.1016/0021-9991(75)90060-1 5. Bouchaud, J.P., Georges, A.: Anomalous diffusion in disordered media: statisti- cal mechanisms, models and physical applications. Phys. Rep. 195(4–5), 127–293 (1990). https://doi.org/10.1016/0370-1573(90)90099-N 6. Embacher, P., Dirr, N., Zimmer, J., Reina, C.: Computing diffusivities from par- ticle models out of equilibrium. Proc. Roy. Soc. London Ser. A 474(2212) (2018). https://doi.org/10.1098/rspa.2017.0694 7. Frenkel, D., Smit, B.: Understanding molecular simulation: From algorithms to applications. 2nd edn. Academic Press (2001) 8. Kipnis, C., Landim, C.: Scaling limits of interacting particle systems, Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sci- ences], vol. 320. Springer-Verlag, Berlin (1999) 9. Li, X., Dirr, N., Embacher, P., Zimmer, J., Reina, C.: Harnessing fluctuations to discover dissipative evolution equations. J. Mech. Phys. Solids 131, 240–251 (2019). https://doi.org/10.1016/j.jmps.2019.05.017 10. Oliveira, F.A., Ferreira, R.M.S., Lapas, L.C., Vainstein, M.H.: Anomalous diffusion: A basic mechanism for the evolution of inhomogeneous systems. Front. Phys. 7, 18 (2019). https://doi.org/10.3389/fphy.2019.00018 11. Öttinger, H.C.: Beyond equilibrium thermodynamics. Wiley Online Library (2005)