{'fraction_non_alphanumeric': 0.04104386606738441, 'fraction_numerical': 0.012490853977213338, 'mean_word_length': 4.162140287769784, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fick and Fokker--Planck diffusion law in inhomogeneous media

We discuss diffusion of particles in a spatially inhomogeneous medium. From the microscopic viewpoint we consider independent particles randomly evolving on a lattice. We show that the reversibility condition has a discrete geometric interpretation in terms of weights associated to un--oriented edges and vertices. We consider the hydrodynamic diffusive scaling that gives, as a macroscopic evolution equation, the Fokker--Planck equation corresponding to the evolution of the probability distribution of a reversible spatially inhomogeneous diffusion process. The geometric macroscopic counterpart of reversibility is encoded into a tensor metrics and a positive function. The Fick's law with inhomogeneous diffusion matrix is obtained in the case when the spatial inhomogeneity is associated exclusively with the edge weights. We discuss also some related properties of the systems like a non-homogeneous Einstein relation and the possibility of uphill diffusion.


Introduction
The modelling of the diffusion of a physical quantity encoded by a density field ρ(x, t) is usually constructed by assuming a continuity equation expressed in terms of the flux vectorial field J(x, t) and a relation between the flux and the density field. The most popular choice is the Fick's law (see [16] for a very introductory discussion) where the positive function D is called diffusion coefficient. In general D = D(ρ, x). When there is a dependence on ρ we obtain a nonlinear equation. For spatially homogeneous systems D does not depend on x.
Let us for simplicity consider the cases of a diffusion coefficient that does not depend on ρ. In many experimental situations [1, 3, 5-8, 13, 19, 21, 22, 24, 25] one should consider a not constant diffusion coefficient D(x). In this cases it is not clear if Fick's law is the correct equation expressing the connection between the density and the flux fields. A different choice is the Fokker-Planck diffusion law (see the books [15,17] for an introduction to the Fokker-Planck equation) J = −∇(Dρ) (1.3) which adds to the standard Fick's law a drift with velocity −∇D, see Section 5.1. In correspondence of these two different assumptions one finds two possible equations for the diffusion problem ∂ρ ∂t = ∇ · (D∇ρ) (1.4) and ∂ρ ∂t = ∆(Dρ) (1.5) which will be respectively called the Fick and the Fokker-Planck diffusion equation; note that they reduce to the same equation if D is constant. These two equations can be studied in Λ × [0, T ] with Λ ⊂ R d and T > 0 with D ∈ C 2 (Λ) and with initial condition ρ(x, 0) = ρ 0 (x) ∈ C 2 (Λ). Possible boundary conditions are Dirichlet or Neumann conditions on ∂Λ. In case Λ is a parallelepiped, it is possible to consider periodic boundary conditions.
In the applied science literature there are many situations in which the two different points of view are assumed. We just mention the paper [27] where the Fick's law is used to study the transport of nutrients in cartilaginous tissues and the paper [22] where it is discussed an experiment in which a not uniform stationary density profile is produced starting from a uniform distribution of particles flowing inside a medium with not constant diffusion coefficient obtained by adding gelatine to water. This experimental observation is obviously in contrast with the Fick's law prediction.
The fact is that, as clearly explained in [25,26], the question "what is the right generalization of the Fick's law to inhomogeneous systems" is too naive. A more detailed knowledge of the microscopic system is necessary to model correctly the macroscopic behavior. In [25] the authors, in particular, discuss a convincing and simple example based on two systems in which a closed box contains a very dilute gas moving through a dense mesh of iron wool. Model one: the iron wool density is uniform and the box experiences a fixed temperature gradient so that the typical particle speed varies continuously throughout the box. Model two: the temperature is uniform, but the iron wool density varies continuously in the box. The systems are designed so that the effective diffusion coefficient, which can be defined as the ratio between the square of the mean free path and the mean free time, is the same function of the space coordinates in the two systems. The authors remark that, since the temperature is uniform in box two and not uniform in box one they expect a stationary uniform particle density distribution in box two and not uniform in box one; indeed, they also deduce Fokker-Planck behavior for the first model and Fick for the second.
Our work is very much in the spirit of [25,26], indeed, we assume the microscopic point of view and prove that two different models behave in the hydrodynamic limit [14,18] respectively according to the Fick and the Fokker-Planck diffusion law. In our modelling particles move in a discrete space and jump from one site to another following an edge. We find the Fick's behavior if the inhomogeneity is associated with edges and the Fokker-Planck one if inhomogeneity is associated with sites.
Our modelling provides a deep physical interpretation of the phenomenon, indeed, it suggests that the Fokker-Planck's law is associated with locally isotropic inhomogeneities, whereas inhomogeneity accompanied to anisotropy results into Fick's behavior. More precisely, suppose that in a small interval of time the number of particles leaving a site of the system is equally distributed among the edges intersecting that site, then the macroscopic behavior is Fokker-Planck. On the contrary, suppose that the number of particles leaving a site are not equally distributed among the edges intersecting that site, but assume also that if two sites connected by an edge are occupied by the same number of particles then the number of particles moving along the bond in the two directions is equal. In such a case the macroscopic behavior is Fick. The second assumption assures that there is no preferred direction along an edge, in particular it rules out the possibility to have external fields acting on the system. We note, finally, that our results are coherent with the simple example discussed in [25]. Consider a small portion of volume in the box one, the number of particles exiting the volume depends on its location due to velocity gradient. But, since the wool mesh is uniformly distributed, particles move with the same speed in all directions, so that the system is locally isotropic and this, accordingly to our results, implies the Fokker-Planck behavior. On the other hand, in box two the non-uniformity of the iron wool distribution breaks the local isotropy and this is why the Fick's behavior is found.
As we mentioned above the main goal of the paper is the derivation of the Fick and Fokker-Planck diffusion laws starting from a microscopic model in which the spatial inhomogeneity is differently implemented. The paper contains also a final section in which we discuss some relevant phenomena connected with inhomogeneous diffusion. In particular, we note that coupling a Fick channel with a Fokker-Planck one with suitable boundary conditions gives rise to the phenomenon of uphill currents, in the sense that the current will flow in the standard downhill direction in the Fick channel, namely, from the higher density end to the lower density one, whereas it will flow uphill in the Fokker-Planck channel. Moreover, in the same section we discuss the validity of an inhomogeneous Einstein relation.
The paper is organized as follows. In Section 2 we introduce the microscopic model and discuss some elementary properties connected to invariant measures. In Section 3 we first introduce the basic notions which are needed to state our main result on the scaling limit which is, indeed, stated in Section 3.3 and proven in Sections 3.4-3.8. Some heuristics and numerical simulations are given in Section 4. Finally, in Section 5 we report some additional remarks as the above mentioned uphill current and Einstein relation.

Models
We discuss here the microscopic structure of our inhomogeneous media.

Preliminaries
At microscopic level we have a graph with vertices V , and directed edges E. The corresponding set of unordered edges is denoted by F . A generic directed edge is denoted by (x, y) ∈ E while an undirected one by {x, y} ∈ F . We consider always finite graphs such that if {x, y} ∈ F then both (x, y) and (y, x) belong to E.
Two vertices x, y ∈ V are said to be neighbors if and only if {x, y} ∈ F . We assume that the graph is connected, namely, for any pair of vertices x, y ∈ V there exists a sequence of unordered edges e 1 , . . . , e n ∈ F such that x ∈ e 1 , y ∈ e n , and e m ∩ e m+1 = ∅ for m = 1, . . . , n − 1. For any x ∈ V we let C(x) ⊂ V be the set of vertices that are neighbors of x. The directed graph (V, E) is called strongly connected if for any pair of vertices x, y ∈ V there exists a directed path going from x to y. We assume that our graphs are always strongly connected.

Random walks and particle systems
We consider one particle performing a Random Walk on the graph (V, E) with rates r(x, y) > 0 when (x, y) ∈ E. We say that the random walk is reversible if and only if there exists a probability measure µ(x) on V such that the detailed balance condition is satisfied. This condition can be satisfied only if {x, y} ∈ F implies that both (x, y) and (y, x) belong to E. We stress again that this will be always true. If the condition (2.6) is satisfied then µ is invariant for the dynamics. This means that if the walker is distributed initially like µ its distribution does not change with time.
The inhomogeneous random walk (IRW) is the Markov jump process on the graph with transition rate from x to y given by where α : V → R + and Q : F → R + are arbitrary functions. We stress that Q is a function on un-ordered edges so that Q({x, y}) = Q({y, x}). To avoid irreducibility problems we assume that such functions are strictly positive. Sometimes we shall consider two particular cases in which the inhomogeneity is associated exclusively either with sites or bonds. The site inhomogeneous random walk (SIRW) is the IRW with Q(e) = 1 for any e ∈ F and the edge inhomogeneous random walk (EIRW) is the IRW with α(x) = 1 for any x ∈ V .
We can pass from the case of one single particle to that of M independent and indistinguishable particles letting η(x) be the number of particles at site x ∈ V and considering η(x)r(x, y) as the rate at which one particle jumps from site x to site y ∈ C(x). More formally, a configuration of particles is an element of the set Ω is the number of particles at x ∈ V and it is usually called the occupation variable at x. If x, y ∈ V and η ∈ Ω such that η(x) ≥ 1, we denote by η x,y the configuration obtained by η letting one particle jump from x to y. This means that, η x,y (x) = η(x) − 1 and η x,y (y) = η(y) + 1 while all the remaining occupation variables remain the same. The stochastic evolution is encoded by the generator , Ω). This is the space of the maps η · : [0, t] → Ω that are right continuous and have limit from the left. We endow this space by the Skorokhod topology [4]. In the following we will denote by P ν the probability measure on D([0, t], Ω) determined by the Markovian stochastic evolution given by (2.8) when the particles are distributed at time 0 according to the measure ν. The corresponding expected value will be denoted by E ν . The probability and the expected value with respect to a probability measure ν on Ω will be instead denoted respectively by E ν and P ν (or simply ν).

Invariant measures
Let us first discuss the case of one single particle. We claim that the class of all the reversible random walks on the graph G indeed coincides with the class of IRW.
Proof. Consider first a random walk with rates (2.7) and consider the probability measure µ(x) = 1/(α(x)Z). Then the detailed balance condition (2.6) holds and the random walk is then reversible and the invariant measure is µ. Conversely consider a random walk for which (2.6) holds. Define then Q({x, y}) := µ(x)r(x, y) = µ(y)r(y, x) and α(x) = µ −1 (x). Then with this choice of the weights formula (2.7) holds and we have therefore an IRW.
For the many particle system, the dynamic conserves the total number of particles and consequently if there are not sources there will be a family of invariant measures depending on the number of particles. On each subset Ω M the dynamics is irreducible and there will be a corresponding unique invariant measure. This is the canonical invariant measure with M particles ν M defined by ν M (η) = 0 if η ∈ Ω M and otherwise By the multinomial theorem, the normalization constant is It is easy to prove that the canonical measure is (2.10) by showing that it satisfies the detailed balance condition for a system of independent IRW where we recall definition (2.9) for the rates c x,y . We note that the average number of particles at site x ∈ V under the stationary measure ν M is (2.13) where Ω x M −k denotes the set {η ∈ N V \{x} , y∈V \{x} η(y) = M − k}. Hence, using the expression of the partition function for M − k particles on V \ {x}, one has yielding (2.13) after some straightforward algebra. An alternative way of looking at this is by labeling the particles. Since the particles are independent, if we distribute initially the particles independently they will be independent at any later time. In particular considering very long times the particles will be independent in the stationary state. Calling X i ∈ V the position of the particle with label i in the stationary state we have that the variables X i are independent and each of them has distribution coinciding with the invariant measure of one single walker described in Lemma 2.1. We have therefore that is exactly the right hand side of (2.13).
It will be more convenient to work with the grand canonical invariant measures that are obtained as special convex combinations of the canonical ones. The family of grand canonical invariant measures is parameterized by a parameter related to the averaged density. Given a function λ(·) : V → R we define an associated inhomogeneous product Poisson measure (2.14) When λ(·) = λ is a constant function we call simply µ λ the corresponding homogeneous product measure. The measure (2.14) satisfies a detailed balance condition similar to (2.12) provided λ(x) = cα −1 (x) for an arbitrary constant c. We obtain in this way a family of grand canonical invariant measures depending on the free parameter c. We note that the average number of particles at site x ∈ V under the measure µ λ(·) is E µ λ(·) [η(x)] = λ(x). We have therefore for the grand canonical stationary measures E µ cα −1 (·) (η(x)) = cα −1 (x). The canonical measures are obtained by the grand canonical ones conditioning on the total number of particles. More precisely we have and the conditioning is independent from the parameter c of the grand canonical measure.

Microscopic and macroscopic observables
In order to perform the scaling limits we need to introduce a general framework and some observables. We will give a microscopic and a macroscopic description of the system. The macroscopic domain Λ is in general a bounded domain of R d , but to avoid dealing with boundary conditions we consider the d dimensional torus [0, 1] d with periodic boundary conditions. The discretization of the macroscopic domain is Λ N := (Z/N ) d ∩ Λ that will be the set of vertices denoted before as V , with edges between nearest neighbors sites. We call respectively E N and F N the oriented and the un-oriented edges of the graph. We denote by L N the generator of the process (2.8) when the underlying graph is (Λ N , E N ). In general, a lower index N is used to denote the fact that the graph that we are considering is the lattice Λ N with the corresponding edges.
A discrete vector field φ is a map φ : . In this case we write φ = ∇f . We use the same notation for the discrete and continuous gradient and divergence since they are one a discretized version of the other. To understand if the symbol means the discrete or the continuous operator we have to observe on which object it is acting.
Given a smooth function f : Λ → R, its discretized version f N on the lattice Λ N is defined by f N (x) = f (x), x ∈ Λ N (with abuse of notation we drop sometimes the index N ). Given a smooth vector field ψ : Λ → R d a natural discretization is obtained for example considering the line integral (3. 16) We have that ψ N is a discrete vector field. We will use repeatedly the following integration by parts formula that can be easily checked. Consider a function f : Λ N → R and a discrete vector field φ N we have We have also the following relationship between sums over ordered edges and unordered ones. Given two discrete vector fields φ N , ψ N we have Note that the right hand side in (3.18) is not ambiguously written since the term to be summed is symmetric in the exchange of x with y. Consider a collection of smooth weight functions Q = (Q 1 , . . . , Q d ) : Λ → (R + ) d . We consider a corresponding discretized version as a weight function Q N taking values on R + and defined on the un-oriented edges by where i in (3.19) has to be fixed in such a way that y = x ± e i where e i is the vector of modulus N −1 and directed as the i coordinate axis. Note that this discretization is very different with respect to (3.16) since in that case ψ N (x, y) is of order 1/N while in this case Q N ({x, y}) is of order one. The general situation that we imagine is that the weights on the edges are the discretization Q N of positive smooth weight functions while the weights on the vertices are the discretization α N of a positive smooth function.
There is a natural mathematical object to be introduced in order to describe the scaling limit of the models. This is the empirical measure π N (η) that is a positive measure on Λ, with finite total mass, i.e. an element of M + (Λ), associated to a configuration of particles η and defined by where δ x is the delta measure. According to this definition, given a continuous function We endow M + (Λ) with the weak topology. We say that a sequence of configurations η (for each N we have a configuration of particles on Λ N , for simplicity of notation the dependence on N is understood) is associated to a density profile . This means that for any continuous function Likewise a sequence of probability measures µ N on the configurations of particles N Λ N is said to be associated with a density profile ρ if for any continuous function f and for any > 0 we have

Large deviations and free energy
We discuss firstly the scaling limit for the empirical measure when the particles are distributed according to a grand canonical invariant measure. We perform the computation for a generic continuous function λ(·) recalling that the grand canonical invariant measure is obtained setting λ(·) = cα −1 (·) for a suitable c. Since the measure is of product type we can discuss this problem following classic strategies and obtaining not only the scaling limit but also the corresponding large deviations asymptotic [18,23]. In this case it is indeed possible to compute exactly the scaled cumulant generating function. Let f be a continuous function; we can compute Since the invariant measure is product, (3.22) can be developed as The last equality follows by the fact that we have in the previous step the corresponding Riemann sums. According to general results on large deviations [23] the corresponding large deviations rate functional, on M + (Λ) endowed with the weak convergence, is given by This gives a rate functional V that is +∞ if the positive measure ρ is not absolutely continuous and when ρ = ρ(x) dx we have where f (ρ) = ρ log ρ is the density of free energy for a system of independent particles.
Here and hereafter with call with the same name an absolutely continuous measure and the corresponding density. The form of the rate functional (3.25) has a structure similar to the one corresponding to a spatially homogeneous system. The only difference is that in (3.25) λ(x) has to be substituted by a constant corresponding to the typical density. Recall instead that λ(x) = cα −1 (x) for the inhomogeneous grand canonical measure.
The functional (3.25) plays the role of a thermodynamic potential and its probabilistic interpretation is that roughly we have where ∼ means closeness in the weak topology and means asymptotic logarithmic equivalence (see [23] for a precise statement). In particular, since V (ρ) = 0 if and only if ρ(x) = λ(x), from (3.26) we can deduce the scaling limit of the empirical measure when the particles are distributed according to the invariant measure. We have indeed that a.e..

Dynamic scaling limit
We deduce in this section the diffusive scaling limit of many independent IRW's on the lattice Λ N . This means that we consider a system of particles defined by the rates (2.9). This system has a diffusive behavior and this means that we have to multiply by N 2 the rates of jump that corresponds to accelerate by the same scale factor the time.
Recall that we consider the situation where the weights on the lattice are inherited by discretization of C 2 inhomogeneities. In particular we fix some C 2 and strictly positive weights Q = (Q 1 , . . . , Q d ) and a C 2 and strictly positive function α. The parameters of the models are fixed discretizing these functions as discussed before.
The proof of our result follows the general strategy outlined in [18] for gradient reversible models with the simplifying feature that we have independent particles. We give an outline of the proof underlying the modifications that we have to do in order to keep into account the spatial inhomogeneity of the models.
Given ν N and µ N two sequences of probability measures on the configuration of particles Ω and such that ν N is absolutely continuous with respect to µ N we introduce their relative entropy defined by A key mathematical object to understand the hydrodynamic behavior of the system is the instantaneous current. This is a discrete vector field depending on configurations of particles and representing the rate at which particles cross the bonds. If c x,y (η) is the rate at which one particle jumps from x to y in the configuration η we have that the corresponding instantaneous current is given by (3.28) For each fixed configuration η this is a discrete vector field. The intuitive interpretation of the instantaneous current is the rate at which particles cross the bond (x, y). Let N x,y (t) be the number of particles that jumped from site x to site y up to time t in the stochastic evolution. The current flown across the bond (x, y) up to time t is defined as This is again a discrete vector field. It is important to point out however that (3.29) depends on the whole trajectory on the time window [0, t] of the system of particles while instead the instantaneous current (3.28) depends just on a configuration of particles η. The importance of the instantaneous current is based on the key observation (see for example [28] Section II 2.3) that is a martingale. Recalling (2.9) we have that the instantaneous current is given by Recall also that to get a non-trivial scaling limit we will accelerate the process by a factor of N 2 so that the instantaneous current (3.31) will be multiplied by N 2 . Our result is the following.
Theorem 3.1. Consider a collection of IRW's associated to the discretization of C 2 smooth and strictly positive weights α and Q. Consider ρ 0 an element of L 1 (Λ, dx). Let ν N be a sequence of probability measures on the configuration of particles Ω associated to the profile ρ 0 in the sense of (3.21) and such that there exists a positive constant K and a constant λ such that When the rates in (2.8) are multiplied by N 2 we have that for any t, for any continuous function f and for any > 0 where ρ(x, t) is the unique weak solution of the equation Proof. The proof is organized into different steps.

Preliminaries
First of all we recall some basic facts about martingales and Markov processes (see for example [18] Appendix 1 Section 5). Consider a function g(s, η) that for each configuration η is C 2 in the time variable s. We have that is a martingale. Moreover we have that is a martingale too. The N 2 factor is due to the rescaling of the time of the process. Since B 0 = M 0 = 0 we have mean zero martingales.
As an example consider the discrete continuity equation for the process that is This is true for any trajectory of the process. Using (3.30) we obtain that is a martingale. A direct computation shows that so that (3.37) is a martingale of the form (3.35) with g(η) = η(x). We recall that in (3.38) the lower index N on the generator simply stress the fact that the underlying graph is the lattice Λ N . Consider a smooth test function f (s, x) : R + × Λ → R and the associated martingale The martingale (3.39) is a martingale of the form (3.35) corresponding to the function The corresponding martingale of the form (3.36) is given by The second term (without he minus sign) on the right hand side of (3.40) is called the quadratic variation of the martingale M f . A direct computation gives This is obtained by the following elementary facts and simple algebraic manipulations. If We have also Finally when {x, y} ∈ F N we have ) .
Since f, α, Q are C 2 , using (3.41), we have that for a suitable constant C. This is a key estimate in our computations that is similar to the estimate that holds in the homogeneous case. This fact allows to extend the results in the homogeneous case to the non-homogeneous one.
With a discrete integration by parts (3.17) the third term on the right hand side of (3.39) (without the minus sign) becomes Using the expression (3.31) of the rates and performing another discrete integration by parts, the second term in (3.43) becomes Inside squared parenthesis in the above formula we have a discrete operator acting on the test function f and not depending on configurations of particles. We need to understand which is the corresponding continuous differential operator. Since our rates are obtained by discretizing smooth functions we obtain with a Taylor expansion of Q that the term inside the squared parenthesis in (3.44) can be written, up to a term O(1/N ), as .
Recall that e i is the vector associated to the i Cartesian axis and having modulus 1/N . The expression inside the squared parenthesis in (3.45) is then equal to up to a infinitesimal term uniform over x, where the divergence ad gradient operators are the continuous ones. We obtain, therefore, that where the residual term R N (t) can be bounded by We used the fact that the dynamics is conservative and we have Since the initial configuration is associated to an integrable profile ρ 0 , selecting as a test function in the definition (3.21) (with ρ replaced by ρ 0 and µ N by ν N ) a function constantly equal to 1, we deduce The general strategy of our proof is the following. Let us call P N ∈ M 1 D([0, t]; M + (Λ)) the probability measure corresponding to the distribution of (π N (η s )) s∈[0,t] ∈ D([0, t]; M + (Λ)). We write shortly P N = P ν N ·π −1 N that means that for any measurable set We will first prove that the sequence of probability measures P N is relatively compact. By Prohorov Theorem this is equivalent to prove that P N is tight. Then we will prove that any possible limiting measure P * of any possible converging subsequence extracted from P N is concentrated on elements of D([0, t]; M + (Λ)) that are absolutely continuous for each s ∈ [0, t] and that satisfy a suitable weak formulation of the equation (3.34). As a final step we prove uniqueness of the weak solution to (3.34). This implies that the whole sequence P N converges weakly to P * = δ ρ id , where we call ρ id the unique weak solution to (3.34). The convergence (3.21) follows by the weak convergence of P N and the fact that ρ id is an element of D([0, t]; M + (Λ)) that is weakly continuous in the time variable.

Tightness
The first step consists in proving that the sequence of probability measures P N is relatively compact. As it is discussed in [18] chapters 4 and 5, we need to prove relative compactness of the marginals for any fixed time and in addition we need to have a control concerning oscillations in time.
Since the total mass is preserved by the dynamics to prove the relative compactness of any marginal it is enough to prove it for the initial condition. Since Λ is compact we need just to control the total mass. In particular we need to prove This is obtained by the same argument used for (3.47).
To control oscillations we use the Aldous criterion (see [18] chapter 4 Proposition 1.6). By the arguments again in [18] chapter 4 Section 2, we need to prove that for any δ > 0 and for any C 2 test function f . In the above formula τ is varying among all the stopping times bounded by t while θ is a real number varying in [0, γ]. We use (3.39) for a function f that does not depend on time and we obtain that (3.49) is true if we have (3.50) The integrand in the upper condition above can be manipulated up to the form (3.44) that according to the subsequent computations can be written up to negligible terms as By the regularity of the functions involved, the integral in the upper condition in (3.50) is bounded by where the inequality follows by the fact that the dynamics is conservative and C is a suitable constant. Here and hereafter we denote by the same letter C a generic constant that may depend just on the weight and the test functions. The values of the constants in different equations may be different. Since we have (3.48) and θ is going to zero we deduce easily the upper condition in (3.50), with an argument like the one for (3.47).
For the lower condition in (3.50) we use Chebysev inequality and get Since τ is a bounded stopping time then M t (τ + θ) − M f (τ ) is again a martingale (with time parameter θ) and having quadratic variation τ +θ τ Γ f (s)ds (see [18]). We have therefore that the right hand side of (3.51) is equal to (3.52) Using (3.42) and the conservative property of the dynamics the last term above is bounded by If we prove that the expected value in the above formula is bounded then, recalling that θ ≤ γ, γ → 0 and N → +∞, we proved also the lower condition in (3.50). This fact does not follow by the fact that ν N is associated to an integrable profile. At this point it is relevant the entropy condition. Recall the basic entropy inequality (see for example [18] appendix 1 Section 8). Given two probability measures µ and ν and a function f we have where β is an arbitrary parameter. We apply this inequality considering ν = ν N , µ = µ λ N , β = N d and finally f (η) = Λ dπ N (η). We obtain where we used the hypothesis on the relative entropy of the initial condition and the explicit form of the generating function of a Poisson distribution. We proved therefore the validity also of the lower condition in (3.50) and we proved therefore (3.49). The proof of tightness is concluded.

Absolute continuity
First of all we observe that the bound on the relative entropy for the initial distribution is still valid with respect to a slowly varying product of exponentials µ λ(·) N . This is obtained using again the entropy inequality (3.54) with ν = ν N , µ = µ λ N , β = 1 and f = log Since we have product measures we can perform explicitly the computations obtaining Since λ(·) is continuous and strictly positive the first term on the right hand side is a Riemann sum and converges while the second one is bounded by assumption.
Considering λ(·) = cα −1 (·) we have that µ λ(·) N is invariant for the dynamics and we have therefore (see [18] appendix 1 Section 9) that H ν N (t)|µ λ(·) N is decreasing in time where ν N (t) is the distribution of particles at time t. This means that for any t ≥ 0 we have H ν N (t)|µ λ(·) N ≤ N d C for a suitable constant C. This is the basic fact on which it is based the argument in [18] Section 1. In particular Lemma 1.6 there, should be rewritten considering in this case I 0 coinciding with the large deviations rate functional V in (3.25).
We deduce that any possible limit point P * of any subsequence in P N is concentrated on elements of D([0, t], M + ) that are of the form ρ(x, s)dx for any s ∈ [0, t] and ρ(x, s) ∈ L 1 (Λ).

Characterization of limit points
Since the sequence of probability measures P N is relatively compact we can extract a converging subsequence. For simplicity of notation we call again P N this converging subsequence and P * its limit point.
Let us consider the martingale (3.39). By the Chebysev and the Doob inequality we have Recalling the bounds (3.42) and (3.55) we have that the right hand side of (3.56) is bounded by 4Ct 2 N d for a suitable constant C and this is converging to zero when N → +∞.
Let us call (3.57) First we recall that by (3.46) we have that is uniformly negligible in probability according to (3.47).
Second we observe that the map that associate to any π(s) ∈ D([0, t], M + (Λ)) the number is a continuous function in the Skorokhod topology of D([0, t], M + (Λ)).
Since by assumption we have that the subsequence P N is weakly converging to P * , by Portmanteau Theorem we have for any > 0 (3.58) By estimates (3.47) and (3.56) the right hand side in (3.58) is zero and this happens for any > 0. We obtain therefore that for any limiting measure P * we have

Uniqueness
In the above steps we proved that any possible limit point P * of a converging subsequence in P N gives full measure to elements π ∈ D([0, t], M + ) such that: π(0) = ρ 0 (x)dx (this follows by the assumption on the initial condition), for any s ∈ [0, t] π(s) ∈ M + is absolutely continuous π(s) = π(s, x)dx and with total finite mass given by Λ ρ 0 (x)dx (this follows by the conservative nature of the dynamics and the initial condition), and finally for any test function f that is C 1 in time and C 2 in space we have Let us now show that there is a unique π(s, x)dx with π(s) ∈ L 1 (Λ) satisfying (3.59). If π 1 , π 2 are two solutions, from (3.59) we readily obtain for π = π 1 − π 2 Here g ∈ C 1 ([0, t] × R d ) is Λ-periodic, as well as all other functions, and vanishes near s = t.
The existence of f in the class above follows from classical results ( [20] chapter 4 Section 5). Then we get in fact t 0 ds Λ g(s)dπ(s) = 0 , for all g as above, yielding therefore π = 0. We conclude therefore that any possible limiting measure P * needs necessarily to be δ ρ id , Since any possible converging subsequence is converging to the same limiting measure we have that the whole sequence P N is converging to δ ρ id .
Any weak solution of the hydrodynamic equation is an element of D([0, t], M + (Λ)) that it is indeed weakly continuous in t i.e. it is an element of C([0, t], M + (Λ)). Indeed by (3.59) we have for any where the constant C depends on the weights, on the function f and on the total mass. The same estimate for any continuous function can be deduced by approximations. The map that associates to any π ∈ D([0, t], M + (Λ)) the real number Λ f dπ(s), for a given time s ∈ [0, t] and a continuous function f , is in general not continuous. We have however that P * is concentrated on weakly continuous paths so that the discontinuity points of this map have P * probability zero and by Portmanteau Theorem we deduce that Λ f dπ N (s) weakly converges to the constant random variable Λ f (x)ρ(x, s)dx where ρ(x, s) is the solution of (3.34). Since weak convergence to a constant random variable implies convergence in probability we deduce (3.33).

Heuristics and numerics
In this section we discuss an heuristic argument which explains the hydrodynamic limits stated in Section 2. Moreover, we shall illustrate numerically the behavior of the SIRW and EIRW stochastic models for many particles in connection with the Fokker-Planck and Fick diffusion equations. In this section, for notation convenience, we shall not use the set Λ N as above, but we will directly work on the graph V = {0, 1, . . . , N }.

Heuristics for the hydrodynamic limit
Consider the SIRW process on V with periodic boundary conditions for M indistinguishable and independent particles. We show that in the limit N → ∞ the evolution of the Markov process density profile converges to that of the Fokker-Planck diffusion problem provided the diffusive scaling is considered. Let a < b be two reals and set z x = a + (b − a)x/N so that z x ∈ [a, b]. Consider a positive function D ∈ C 2 ([a, b]) and set α(x) = D(z x ) for x ∈ V . Denote by η x (t) the particle profile at time t, informally speaking, η x (t) is the average number of particles occupying the site x at time t. The change of the number of particles at site x in a small interval ∆t can be computed as This equality can be rewritten as Thus, if time is rescaled as t/N 2 → t (diffusive scaling), then in the limit N → ∞ the particle density profile η x (t)/(1/N ) will tend to a function ρ(z, t) solving the equation We consider the EIRW process on V = {0, 1, . . . , N } with periodic conditions for M indistinguishable and independent particles. and we use the same notation introduced above in the SIRW process case. We let Q({x, x + 1}) = D((z x + z x+1 )/2) be the rate associated with the edge {x, x + 1} for x ∈ V , where {N, N + 1} is identified with {N, 0}. The change of the number of particles at site x in a small interval ∆t can be computed as and, hence, This equality can be rewritten as Thus, if time is rescaled as t/N 2 → t (diffusive scaling), then in the limit N → ∞ the particle density profile η x (t)/(1/N ) will tend to a function ρ(z, t) solving the equation which is the Fick diffusion equation.

Numerical solution of the diffusion equations
We discuss some numerical results for the periodic boundary condition Fick and Fokker-Planck diffusion problem on [0, 1]×[0, 1] with the following choices of the diffusion coefficient: and The numerical solution was found using the NDSolve routine in Mathematica. The initial condition is u 0 (z) = 6z(1 − z) in all simulations. We did not use a constant profile as initial condition, since that would have been the stationary solution of the Fick diffusion process so that no dynamics would have been observed.
Note that in the case (4.61), which mimics a discontinuous diffusion coefficient, the Fick diffusion problem has a constant profile as stationary solution, whereas the Fokker-Planck problem tends to profile rapidly varying in correspondence of the diffusion coefficient "discontinuities".
The stationary solutions of the Fick and Fokker-Planck equations can be derived explicitly. In the Fokker-Planck case we have that at stationarity (Du) must be constant. But, for mass conservation, it must indeed be equal to zero, so that at stationarity u(z) = c/D(z) where the constant c is such that where, we recall, u 0 denotes the initial condition. In the Fick case we have that at stationarity Du must be constant. But, for mass conservation, it must indeed be equal to zero, so that the stationarity solution is the constant 1 0 u 0 (z) dz.

SIRW process and Fokker-Planck equation
We now compare the evolution of the SIRW process introduced in Section 2 to that of the Fokker-Planck diffusion equation on [0, 1] × [0, 1]. The stationary profile can be discussed explicitly, indeed, in Section 2.3 we have stated that at stationarity the average number of particles at site where, we recall, M is the total number of particles. Comparing (4.62) and (4.63) we have that, for N large, b ≈ c/N . Hence, for N large the stationary particle density profiles (b/α(x))/(1/N ) of the SIRW process is a very good approximation of the Fokker-Planck stationary solution c/D(z).
For the time dependent results we simulate the stochastic model as follows: we let z x = x/N and recall α(x) = D(z x ) for x ∈ V . Recalling η x (t) is the number of particles at site x and time t, we extract an exponential random time τ with parameter N x=0 2α(x)n x (t) and set the time equal to t + τ . We associate the probability 2α(y)n y (t)/ N x=0 2α(x)n x (t) to each site y ∈ V and select at random a site according to such a distribution. We move a particle from the selected site to one of the two adjacent sites with probability 1/2. To compute the stationary particle profile we let the system evolve for 10 3 full sweeps (in one sweep M particles are moved). Then, we average the value of the number of particles occupying each site of the lattice by considering one configuration each 10 sweeps. The numerical experiment is stopped after about 10 5 more sweeps.
In Figure 4.4 we compare the stationary solution of the Fokker-Planck diffusion processes with the stationary particle profile of the Random Walk. The stationary particle profile is divided times the spacing 1/N to get the stationary particle density profile and is divided times M since the Fokker-Planck diffusion equation has been solved with an initial state having total mass equal to one. The match is perfect.
In Figure 4.5 we compare the evolution of the Fokker-Planck diffusion processes with the Random Walk particle profile. As for the stationary state, the Random Walk particle profile has been divided times the spacing 1/N to get the particle density profile and divided times M since the Fokker-Planck diffusion equation has been solved with an initial state having total mass equal to one. Moreover, the time measured in the stochastic evolution has been divided times N 2 . Averages have been computed by considering 50 independent realizations of the process and averaging the particle distribution at equal times. The match is striking.

EIRW process and Fick diffusion equation
We now compare the evolution of the EIRW process introduced in Section 2 to that of the Fick diffusion equation on [0, 1] × [0, 1]. In this case the stationary state is trivial, indeed, we compute the stationary particle distribution profile as outlined for the SIRW case and we find that it is constant with very high precision.
In Figure 4.6 we compare the evolution of the Fick diffusion processes with the Random Walk particle profile. As for the stationary state, the Random Walk particle profile has been divided times the spacing 1/N to get the particle density profile and divided times M since the Fick diffusion equation has been solved with an initial state having total mass equal to one. Moreover, the time measured in the stochastic evolution has been divided times N 2 . Averages have been computed by considering 50 independent realizations of the process and averaging the particle distribution at equal times. The match is striking.

Miscellany
In this section we collect some interesting remarks on the behavior of the system that we have studied above.

Einstein relation
A very general modelization of the presence of an external field is obtained perturbing the rates as follows. Let φ : Λ → R d be a smooth vector field that acts on particles. The action of the field is encoded in the perturbed transition rates that are defined as that correspond to movements of the particles with an associate negative work of the field are decreased.
Let us first discuss the influence of an external field in the case of spatially homogeneous models [2]. The hydrodynamic scaling limit of diffusive particle systems under the action of a weakly asymmetric external field is associated to equations of the form (5.65) The symmetric and positive definite matrix D is the diffusion matrix while the symmetric and positive definite matrix M is the mobility matrix. For independent particles we have that the diffusion matrix coincides with the identity matrix D = I while instead M = ρI.
In the homogeneous case a relevant thermodynamic relationship is the so called Einstein relation between the diffusion matrix and the mobility given by that says that the two matrices D and M are proportional and the proportionality factor is the second derivative of the density of free energy f (that is f (ρ) = ρ log ρ in the independent particles case as discussed after (3.25)).
It is very natural to interpret this equation introducing the space dependent diffusion matrix D(x, ρ) = α(x)Q(x) and the space dependent mobility matrix M(x, ρ) = α(x)ρ(x)Q(x). Note that they satisfy the Einstein relation for each x ∈ Λ. Indeed recalling that the density of free energy is f (ρ) = ρ log ρ for independent particles we have With this identification we have that the inhomogeneity determines space dependent diffusion and mobility matrices. The form of these matrices depend both on the weights on the edges and on the weights on the vertices. The spatial inhomogeneity of the material generates however also an external field that depends just on the site inhomogeneity. This external field is exactly −(1/2)∇ log α. We show that this interpretation is correct. This is done switching on a weak external field and showing that the hydrodynamic equation is modified with the appearance of a term proportional to the mobility matrix M(x, ρ) like in the homogeneous case (5.65). In presence of an external field the rates are modified according to (5.64) and correspondingly the instantaneous current becomes Recall that the values of φ N are infinitesimal (3.16) so that we have The instantaneous current is therefore Substituting (5.68) to j η in the second term in (3.43) and ignoring negligible terms we obtain the extra factor With computations similar to the ones in the proof of Theorem 3.1 we have that the term inside squared parenthesis in the above formulas coincides up to uniform infinitesimal terms with 2Q(x)φ(x) · ∇f (x, s) .
This means that the hydrodynamic equation in presence of a weak external field becomes We deduce that M(x, ρ) plays the role of the mobility matrix and we obtain a version of the Einstein relation in the non-homogeneous framework.

Alternative proof
Since we are considering a system of independent particles we can obtain an alternative proof under some special initial conditions. In particular we consider the case when the initial condition is obtained with identical particles distributed independently. Note that Theorem 3.1 covers much more general initial conditions. In this special case, the collective behavior of the occupation variables can be deduced by the scaling behavior of one single particle. We could however not find a specific reference for the scaling limit of one single IRW. The following is a sketch of the general argument that can be used once the scaling limit of one single IRW is established.
Consider the initial condition ρ 0 in the hydrodynamic equation (3.34) and define the corresponding probability measureρ 0 (y) = ρ 0 (y)/ Λ ρ 0 (x)dx. We consider a sequence of probability measures p N on Λ N such that where the convergence is the weak one.
A simple generalization of the law of large numbers says the following. Suppose that for each natural number N we have a random variable Y N 1 taking values on a Polish space A and such that the law of Y N 1 is converging weakly to γ ∈ M 1 (A) when N diverges. We called M 1 (A) the set of probability measures on A with the Borel sigma algebra. For each N let us consider (Y N i ) i∈N be a collection of i.i.d. random variables each of them having the same distribution of Y N 1 . Then we have that where the convergence is the weak one in probability (indeed even a.e.). More precisely the above statement means that for any continuous and bounded function f : We consider at time zero Λ ρ 0 (x)dx N d particles independently distributed and each of them distributed on Λ N according to p N . Let X N i (0) be the random position in Λ N at time 0 of the particle number i. We consider D([0, t]; Λ) the Skorokhod space of trajectories. The trajectory of the particle number i is denoted by X N i (·) := (X N i (s)) s∈ [0,t] . This is a random variable taking values on D([0, t]; Λ). We consider each particle evolving with an IRW with rates of jump accelerated by a factor of N 2 .
We assume in this argument that the law of the trajectory of one single particle converges to the law Pρ 0 of a diffusion process (see next Section 5.3 ) with initial distributionρ 0 and Kolmogorov evolution equation for the distribution given by the hydrodynamic equation (3.34) (with initial conditionρ 0 ). This is an assumption because we could not find a precise reference for this result.
We have therefore the convergence (5.69) that in this specific case implies that a.e., and therefore in probability, we have Consider a continuous and bounded function f : Λ → R and the functional F : D([0, t], Λ) → R defined by F (X(·)) := f (X(s)) where s ∈ [0, t] is a fixed time. The functional F is not continuous with respect to the Skorokhod topology. We have however that under the probability measure Pρ 0 the set of discontinuous points of this functional has probability zero. This is because the probability Pρ 0 is concentrated on continuous paths. We can therefore deduce by Portmanteau Theorem that a.e., and therefore in probability, we have the convergence whereρ(s) is the solution of (3.34) with initial conditionρ 0 . We can deduce the hydrodynamic behavior of the model observing that

Reversible diffusions
As we observed in the previous section, in the case of independent particles the hydrodynamic equation describing the collective behavior of several particles is linear and coincides with the equation of the evolution of the probability distribution of one single particle. Since the scaling limit of one single particle is a diffusion process and since our discrete models are reversible it is natural to compare the class of hydrodynamic equations that we obtained with the possible Fokker Plank equations associated to reversible diffusions. At the microscopic level we obtained that the reversibility condition has a geometric interpretation. We have indeed that the models are reversible if and only if the rates are chosen according to some weights associated to the edges and the vertices of the graph (see Lemma 2.1). In the case of continuous diffusion process we have a similar geometric characterization of reversibility, indeed reversible diffusions can be parameterized by a positive function and a symmetric and positive definite matrix, that can be interpreted as the metric tensor. These are the continuous counterparts of the discrete weights on the graph.
We refer to [15,17] for the basic facts about diffusion processes. For simplicity we consider the processes on R d instead that on the torus. Consider a diffusion process of the form where Note that while in the equation (5.71) appears the matrix B, the evolution of the probability distribution depends just on the symmetric matrix BB T . The condition of reversibility (see [15,17]) is that the vector We have therefore that all the reversible diffusion processes can be parameterized in terms of the function ψ and the symmetric and positive definite matrix BB T . This is because you can fix arbitrarily these two objects and then A is completely determined by (5.73). If we use instead the positive function α related to ψ by ψ = log α and the symmetric positive definite matrix Q(x) = BB T (x)α −1 (x)/2 we have that the Fokker Plank equation (5.72) is given by that is exactly of the type of our hydrodynamic equation (3.34). It is important to note however that in (3.34) the matrix Q has to be diagonal while instead this is not the case in (5.74). As we will discuss in the next section this is due to the special lattice that we are considering in Theorem 3.1. We can obtain non diagonal matrices considering different lattices.

Different lattices
Here we show that we obtained just equations with diagonal matrices Q since we are considering a squared lattice. We briefly discuss how to handle different situations obtaining non diagonal matrices Q. From the proof of Theorem (3.1) we known that the basic computation to identify the limiting equation is to approximate up to uniformly infinitesimal corrections the term inside square parenthesis in (3.44) that is The generalized framework that we consider now is a lattice having vertices coinciding again with Λ N but having more edges than the usual square lattice. This corresponds to allowing more possible jumps to the particles. The graph on which the particles are evolving is obtained as follows. We start with Z d with more edges with respect to the usual ones that are connecting just the minimal distance vertices. The collection of directed edges exiting form any vertex x ∈ Z d are of the form (x, x +ṽ i ) whereṽ i for i = 1, . . . , k is a collection of vectors such that x +ṽ i ∈ Z d . Since we are always requiring that an un-oriented edge can be crossed on both directions then k has to be necessarily an even number and for any vector v i there should be a corresponding label j such thatṽ j = −ṽ i so that both (x, x +ṽ i ) and (x +ṽ i , x) are elements of the directed edges E N . The lattice that we consider is obtained scaling by a factor of N −1 this lattice. In particular we call v i := N −1ṽi .
We have therefore that on each lattice site x ∈ Λ N there are k different edges incident that correspond to k possible jumps of one particle from x to x + v j , j = 1, . . . , k. In the case of the square lattice we had k = 2d and each v j is equal to ±e i for some i. Note that we have now |C(x)| = k. More general frameworks are of course possible but for simplicity we restrict to this generalization.
We need to give weights to the vertices and the edges of the lattice suitably discretizing smooth objects. The weights on the vertices are associated as before computing a smooth function α on the corresponding point. For the edges we need to generalize the construction done before.
We consider a smooth metrics Q(x) that is a symmetric and positive definite d×d matrix depending in a regular way (C 2 for example) on the continuous variable x ∈ Λ. We associate the weight to an edge of the form {x, x + v i } as (5.76) The appearance of theṽ i vectors above is due to the fact that we have |v i | ∼ 1/N (since the vectors without tilde are comparable with the mesh of the lattice) and we want that the weights to be associated to the edges are not infinitesimal in N but are of order one. The last equality in (5.76) is just the definition of a shorthand for the weights. With a suitable Taylor expansion we get that (5.75) coincides up to uniformly infinitesimal terms with where H(x) is the Hessian matrix at x of the function f having elements (H(x)) l,m = ∂ x l ∂ xm f (x). Recall that k is an even number ad if v i is the vector associated to a possible jump then also −v i is a vector associated to a possible jump. Due to this, we have that the leading term in the product in (5.77) that is is identically zero. This is because we can pair the edges exiting from x in such a way that if the label i is paired to the label j then v i = −v j and consequently Q i (x) = Q j (x). Of the remaining three terms obtained when we develop the product in (5.77) we have that one is infinitesimal. The two relevant ones that survive are The above expression coincides up to uniform infinitesimal terms with where the matrix Q is defined as With the same arguments of the proof of Theorem 3.1, but using this expansion, we can prove that the limiting equation is again of the form (3.34) but the matrix Q is given by (5.79) that in general is non-diagonal.

Uphill currents
A current is said to move "uphill" when particles migrate up the gradient, namely towards regions of higher concentration, thus violating the basic tenets of Fick's law of diffusion.
The onset of such uphill currents can be traced back to the action of an external field, to the presence of mutual interactions in a multi-component system or, for single-component systems, to a phase transition, and was recently investigated in a variety of lattice gas models, cf. Refs [5,[9][10][11][12]. We look, here, at the case where two inhomogeneous diffusion processes take place in two intervals of length L > 0, for two concentration functions v, u, being connected by conditions of equality of concentration and of flux at the two endpoints. The latter is meant in the sense that the outflux of v equals the influx of u. However, v solves Fick's equation, while u solves a Fokker-Planck type equation. The diffusivities are assumed to be piecewise constant.
We consider the stationary case, see also Refs. for a more general discussion about the observation of uphill currents.
Thus the problem is, in a distributional formulation, We assume here D 1 = D 2 ; see also Remark 1. We refer to the following weak formulation of this problem: find v ∈ H 1 (0, L), u ∈ L ∞ (0, L) such that Du ∈ H 1 (0, L) and Note that more generally one should write e.g., (5.94) as (D 1 u) x (a−) = (D 2 u) x (a+) .
but this is not relevant under our assumption of piecewise constant D. A similar remark applies to (5.84), (5.85), which indeed are valid in a pointwise sense. Clearly problem (5.80)-(5.85) is invariant for multiplication by a constant, and always has the null solution. Therefore for the sake of precision we'll impose also the following normalization condition v(0) = 1 .