From Particle Fluctuations to Macroscopic Evolution Equations: The Case of Exclusion Dynamics ? ??

Abstract

We consider the symmetric simple exclusion process in one dimension as an example of a stochastic particle process exhibiting anomalous di usion; this process is so slow that the mean-square displacement of a tagged particle is sub-linear. This implies that the standard meansquare displacement method (Einstein relation) cannot be applied to determine the di usivity for the associated macroscopic equation describing the hydrodynamic limit. We demonstrate that a recent approach developed by the authors and collaborators, based on the covariation of the uctuations, is applicable to this process and does not only deliver the transport coe cient, but the entire evolution operator associated with the formulation of the macroscopic equation as an entropic gradient ow. Furthermore, the approach relies on uctuations of the density eld (macroscopic uctuations) as opposed to particle level data. This data could in principle be obtained from experimental observations.

Keywords

• Anomalous di usion Fluctuation-dissipation relation Fluctuating hydrodynamics

-

Anomalous di usion is characterised by a nonlinear dependence of the mean squared displacement (MSD) as a function of time, in contrast with the classical di usion of Brownian particles, where MSD / t (at large times). Examples of anomalous di usion include sublinear di usion observed in biological applications, such as di usion in the cell nuclei and membranes; other examples occur for instance in disordered media, astronomy, uid mechanics and network theory [10]. A classic paper giving a broad overview on this topic is [5].

Anomalous di usion poses numerous theoretical di culties (see, for example, [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 di erential equations (PDEs). A common procedure is to assume phenomenologically the structure of the governing equation and then determine the transport coe cients appearing in it from particle simulations. Prominent methods to determine di usion coe cients are Green-Kubo relations (based on time-correlation functions) and Einstein relation (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 di usion 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 di usion, namely for the symmetric simple exclusion process. More speci cally, this approach delivers the operator of the macroscopic evolution equation (in discretised form), written as a gradient ow of the entropy functional, without making further phenomenological 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, e ectively blocking them. Thus (albeit the hydrodynamic limit is a linear di usion 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 coe cient only. In this contribution, we show that the method can be extended using the approach of [9] to 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 ow, and that uctuations around such limit are Gaussian. This means that the density evolution has the following thermodynamic structure in the limit of in nitely many particles

S ( ); where K is a symmetric operator, and S is the variational derivative of the entropy functional, while for large but nite number of particles N , the evolution shall be formally described by the stochastic partial di erential equation of the form

S( ) + r C p2K( )W_ t;x;

N with W_ 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 uctuations and dissipation: the same operator K appears in the deterministic (macroscopic) part and in the uctuations. 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 di usion 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 ow of the entropy (further details in Section 2), with S being the mixing entropy, i.e.,

S( ) =

Z [ log + ( 1 ) log( 1 )] dx in dimensionless form, and K the operator K( ) := div( ( 1 )r ), i.e., 2

Gradient Flow Structure of the SSEP In this section, we explain why the formulation ( 4 ) is a gradient ow. 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 ow of a functional S is de ned as the evolution (@t ; s2)K 1 =

Z where the geometry of the gradient evolution is de ned 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 de ne 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 de ne formally an inner product on the tangent space,

Z

1 (s1; s2)K 1 := r 1r 2 dx; where sj = div( r j ) for j = 1; 2. (The tangent space can be interpreted as time derivatives satisfying the continuity equation.) The last equation also de nes the class of test functions in Eq. ( 5 ).

We now compute the expressions of the gradient ow evolution for the SSEP. On one hand, we nd with undetermined 1 and 2 and ( 6 ) ( 7 ) ( 8 ) ( 10 ) for the expression on the left of ( 5 ) div( r 1) and s2 =

div( r 2) rf = 1

r rf dx; r log 1 r 2 dx =

Z r 1 1 r 2 dx: ( 9 ) As, in the absence of topological obstructions, any gradient vector eld rf can be written in the form (by solving div 1 1 r 2 =

f ), the combination of ( 8 ) and ( 9 ) yields With the rst continuity equation in ( 8 ) and integration by parts, this becomes which is the weak form of the hydrodynamic limit ( 4 ) associated with the SSEP.

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 di erent ways to obtain the linear di usion equation ( 4 ), @t = as gradient ow. For example, we can use the standard Wasserstein metric and the standard (Boltzmann) entropy S = R 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 ow formulation: here the di erences between di erent proceses become visible, even if they have the same macroscopic limit. For example, Brownian particles de ne the classic (Boltzmann) entropy gradient ow, and the SSEP de nes the thermodynamic ow of the mixing entropy. Even if the hydrodynamic limit of the two processes is the same, their uctuations are di erent, as evident from Eq. ( 2 ). The entropy-metric pair also records these di erences. 3 3.1

Computational Framework

Discretisation of the Dissipative Operator

The partial di erential equation ( 4 ) and the operator K that the PDE encodes are in nite dimensional objects. Thus, their full characterisation from particle level data either requires in nite amount of data, which is of course unattainable from numerical simulations, or, alternatively, optimal ts from a nite library of in nite 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 nite dimensional object), as previously proposed by the authors in [9]. Speci cally, we consider a nite element approximation scheme of the density eld and the thermodynamic driving force F := S= of the form, (x; t)

X a(t) a(x); and

F (x; t) a

X Fa(t) a(x); a ( 11 ) where f a(x)g is a basis of functions with local support. These satisfy the socalled partition of unity Pa a(x) = 1, linear eld reproduction Pa 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 nite element shape functions provide su cient regularity to approximate the above elds 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 (jx xaj ; 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 a; bi Fa; for all b; (13) where the bracket h ; i denotes the integral inner product in the L2 space. The resulting matrix hK a; bi represents the discretised dissipative operator, which we will seek to estimate from particle data, with the strategy outlined in Section 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 (14) Equivalently, this constraint can be written as Pb hK b; ai = 0 for symmetric operators, as the one associated to SSEP or any other purely dissipative processes with the structure of Eq. ( 1 ).

Despite the nite-dimensional nature of hK a; bi, this operator is, a priori, a function of the entire density pro le, or, for the already assumed nite element approximation, a function of all the nodal density values. Consequently, its domain is extremely vast and dependent on the speci c 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; bi is only non-zero for the cases where node a is equal or near to node b (concretely, a = fb 1; b; b + 1g), and that only the local density pro le is needed to determine the discretised operator. This important observation allows us to approximate the density pro le by means of a Taylor expansion, and rewrite Eq. (13) as

X a

X a2fb 1;b;b+1g

D

E K( a+r ja(x xa)+:::) a; b Fa: (15) For SSEP, we will truncate the Taylor expansion after the liner term, rendering hK a; bi a function of only two variables: a and r ja. From a practical perspective, this implies that only linear pro les, encompassing di erent 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 pro les, i.e., not necessarily linear ones. For this, we will use a forward Euler time discretisation scheme, and rewrite Eq. (15) as

M i+1 i

= KiFi: t where the superscript "i" denotes the time step ti, and all quantities have been written in the form of vectors and matrices. Speci cally, i = ia N 1 is the vector of densities at the nodal positions, r i = r jia N 1 is the vector of density gradients, and Fi = Fi i = Fai i N 1 is the vector of thermodynamic 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 (16)

M = (h a; bi)N

N = ;

(17) where K0, Kp and Km represent the three operator entries hK a; bi = Kba with a = b 1; b; b + 1, respectively, Xai = ( ia; r jia), and the periodic boundary conditions have been considered. For any given density pro le, hK a; bi can be computed from the covariation of the rescaled local density uctuation as [9] hK a; bi = lim h&0 2h 1

E [(Y a (t0 + h)

Y a (t0)) (Y b (t0 + h)

Y b (t0))] ; (19) where t0 is an arbitrary initial time that does not a ect the result for the operator entries as long as the system has reached a local equilibrium, E [ ] represents the expectation, and Y are local rescaled density uctuations, de ned as Y = lim !0 h ; i =p with = E [ ]. In practice, the particle system is of course simulated with a nite number of sites, and thus the limit as ! 0 in Y is numerically approximated. Similarly, the right hand side of Eq. (19) is computed for a nite time step h, that is su ciently small to capture the limit numerically, yet, su ciently 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 nite realisations of macroscopically identical density pro les. 3.3

Polynomial Regression of the Discretised Operator Equation (19) will enable the calculation of the discrete operator Kba =

D

K( a+r ja(x xa)) a; b

E for a nite set of linear pro les, resulting in a discrete space Vdiscr of a and r ja. Yet, its use in Eq. (15) to simulate the evolution of arbitrary density pro les will require evaluating Kba at di erent 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 2 fb 1; b; b + 1g should satisfy the constraint given by Eq. (14) in order for the evolution to be mass preserving. Although such identity may not be exactly satis ed for the particle-inferred operator in Vdiscr, we here perform a polynomial regression on such data that, by construction, exactly satis es the conservation constraint. More speci cally, denoting by j (X) (j = 1; 2; : : : ; n) the tting basis, with X = ( ; r ), we approximate the three non-zero entries of the discretised operator as K0 = Km = n X aj0 j (X); j=1 n X am

j j (X); j=1 Kp =

K0

Km = n X aj0 + ajm j=1 j (X); (20) (21) (22) where K0; Km and Kp have been de ned after Eq. (18). For simplicity, we use polynomials up to order-k as the basis functions, i.e., j (X) 2 f1; ; r ; 2; r ; (r )2; :::; k; :::; (r )kg. Speci cally, we will use k = 4 to t the discrete operator data for SSEP.

Towards the goal of nding the optimal tting parameter aj0 and ajm, we denote the data points from the numerical simulation as fX0; K0 i i g with i = 1; :::; N0 for the diagonal entry Kbb, fXim; Km i g with i = 1; :::; Nm for the subdiagonal entry Kb;b 1, and fXip; Kipg 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 de ne the loss function L as the sum of the least square error for all three entries, i.e.,

L = 2 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., 0a0

K0 + pT p a0 + am

+ Kp = 0; mT ( mam

Km) + pT p a0 + am + Kp = 0: We remark that the above equations use the matrix/vector forms for the tting parameters a0 = (aj0)n 1, am = (ajm)n 1, the operator entries K0 = (Kj0)N0 1, Km = (Kjm)Nm 1, Kp = (Kjp)Np 1, and the basis functions 0 = ( j (Xi0))N0 n, m = ( j (Xim))Nm n, p = ( j (Xip))Np n. Denoting = a0 am 0T 0 + pT p

pT p a = ; and b = mT 0TK0 mTKm pT p m + pT p pTKp pTKp ; ; Equations (24){(25) can be equivalently written in compact form as a = b; delivering a simple linear system of equations from which the tting parameters a can be obtained. Alternatively, a constrained least-square solver available in some commercial packages could be used to nd the optimal tting parameters. 4 4.1

Numerical Results

Discretised Operator from Particle Fluctuations In order to compute the operator entries numerically from density uctuation using Eq. (19), we use the lattice kinetic Monte-Carlo (KMC) method [4], implemented 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 di erent initial linear (or piecewise linear) density pro les, which results in the discrete space Vdiscr shown in Fig. 1, with data points within the range of 2 [0:05; 0:95] and r 2 [ 3; 3].

For each initial pro le, 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 uctuations 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; bi with a 2 fb 1; b; b + 1g are shown in Figs. 2(a){(c) (blue circles) together with the analytical surfaces, obtained by performing a Taylor expansion on xa as hK hK b 1; bi = hmr b 1; r bi = b+1; bi = hmr b+1; r bi = + O( x2); + O( x2); with ma = a(1 a); hK x + O( x2); (31) The diagonal entry Kbb is symmetric with respect to = 0:5 and is independent with r , while the two o -diagonal entries are non-symmetric and have mb 1 x 2 (r ja)2 + O(r 2):

K( a+r ja(x xa)) a; bE with a 2 fb 1; b; b + 1g as a function of density and density gradient r (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 2 fb 1; b; b + 1g, 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 quanti ed in Figs. 2(d){(f), where the relative error of the t, following Section 3.3, is also shown. As it may be observed in the gure, the fourth order polynomial regression with mass conservation strictly enforced signi cantly 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 t. 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 pro le. Concretely, we consider an initial density (x; 0) = 0:5 0:3 cos(4 x) over the unit interval x 2 [0; 1], discretised with N = 50 shape functions (i.e., x = 1=N = 0:02). This pro le is evolved 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 at density pro le at the average initial density.

Snapshots of the evolution are represented in Fig. 3, together with the analytic PDE (with identical spatio-temporal discretisation, and discretised operator 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 excellent agreement between the three curves, making them almost indistinguishable in the gure.

The present contribution therefore shows that the computational strategy of [9] is able to extract the dissipative operator for the symmetric simple exclusion process, an anomalous di usion process, and hence, it can be used to predict the non-equilibrium macroscopic evolution of arbitrary initial pro les.

References

1. Adams , S. , Dirr , N. , Peletier , M. , Zimmer , J.: Large deviations and gradient ows . Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci . 371 ( 2005 ), 20120341 ( 2013 ). https://doi.org/10.1098/rsta. 2012 .0341

2. Arratia , R.: The motion of a tagged particle in the simple symmetric 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 uid 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 di usion in disordered media: statistical 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 di usivities from particle 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 Sciences] , vol. 320 . Springer-Verlag, Berlin ( 1999 )

9. Li , X. , Dirr , N. , Embacher , P. , Zimmer , J. , Reina , C. : Harnessing uctuations 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 di usion: A basic mechanism for the evolution of inhomogeneous systems . Front. Phys. 7 , 18 ( 2019 ). https://doi.org/10.3389/fphy. 2019 .00018

11. Ottinger, H.C.: Beyond equilibrium thermodynamics . Wiley Online Library ( 2005 )