{'fraction_non_alphanumeric': 0.07244795025714805, 'fraction_numerical': 0.021117560483537565, 'mean_word_length': 3.5124990672337884, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hydrodynamic Limits and Propagation of Chaos for Interacting Random Walks in Domains

A new non-conservative stochastic reaction-diffusion system in which two families of random walks in two adjacent domains interact near the interface is introduced and studied in this paper. Such a system can be used to model the transport of positive and negative charges in a solar cell or the population dynamics of two segregated species under competition. We show that in the macroscopic limit, the particle densities converge to the solution of a coupled nonlinear heat equations. For this, we first prove that propagation of chaos holds by establishing the uniqueness of a new BBGKY hierarchy. A local central limit theorem for reflected diffusions in bounded Lipschitz domains is also established as a crucial tool.


Introduction
The original motivation of this project is to study the transports of positive and negative charges in solar cells. We model a solar cell by a domain in R d that is divided into two disjoint sub-domains D + and D − by an interface I, a (d − 1)-dimensional hypersurface, which can be possibly disconnected. D + and D − represent the hybrid medium that confine the positive and the negative charges, respectively. At microscopic level, positive and negative charges are initially modeled by N independent reflected Brownian motion (RBM) with drift on D + and on D − respectively. (In this paper, they are actually modeled by N independent random walks on lattices inside D + and D − that serve as discrete approximation of RBM with drifts.) These random motions model the transport of positive (respectively negative) charges under an electric potential (See Figure 1). Figure 1: I is the interface of D + and D − These two types of particles annihilate each other at a certain rate when they come close to each other near the interface I. This interaction models the annihilation, trapping, recombination and separation phenomena of the charges. The interaction distance is of microscopic order ε where N ε d is comparable to 1, and the intensity of annihilation per pair is of order λ/ε where λ ≥ 0 is a given parameter. This means that, intuitively and roughly speaking, according to a random time clock which runs with a speed proportional to the number of pairs (one particle of each type) of distance ε, we annihilate a pair (picked uniformly among those pairs of distance less than ε) with an exponential rate of parameter λ/ε. The above scaling guarantees that in the limit, a nontrivial proportion of particles are annihilated in any open time interval. We investigate the scaling limit of the empirical distribution of positive and negative charges; that is, the hydrodynamic limit of this interacting diffusion system. We show that in the macroscopic level, the empirical distribution converges to a deterministic measure whose density satisfies a system of partial differential equations that has non-linear interaction at the interface.
The study of hydrodynamic limits of particle systems with stochastic dynamics is of fundamental importance in many areas. This study dates back to the sixth Hilbert problem in year 1900, which concerns the mathematical treatment of the axioms of physics, and to Boltzmann's work on principles of mechanics. Proving hydrodynamic limits corresponds to establishing the law of large number for the empirical measure of some attributes (such as position, genetic type, spin type, etc.) of the individuals in the systems. It contributes to our better understanding of the asymptotic behavior of many phenomena, such as chemical reactions [31], population dynamics [20,32,34], super-conductivity [40], quantum dynamics [21], fluid dynamics [24], etc. It reveals fascinating connections between the microscopic stochastic systems and deterministic partial differential equations that describe the macroscopic pictures. It also provides approximations via stochastic models to some partial differential equations that are hard or impossible to solve directly.
Since the work of Boltzmann and Hilbert, there have been many different lines of research on stochastic particle systems. Various models were constructed and different techniques were developed to establish hydrodynamic limits. Among those techniques, the entropy method [25] and the relative entropy method [40] are considered to be general methods. Unfortunately these methods do not seem to work for our model due to the singular interaction near the interface.
Many models studied in literature are conservative, for example exclusion processes [27,28] and Fleming-Viot type systems [8,9]. Reaction-diffusion systems (R-D system) constitute a class of models that are typically non-conservative. These are systems which have hydrodynamic limits of the form ∂u ∂t = 1 2 ∆u + R(u) (a reaction-diffusion equation, or R-D equation in short), where R(u) is a function in u which is thought of as the reaction term. R-D systems arise from many different contexts and have been studied by many authors. For instance, for the case R(u) is a polynomial in u, these systems contain the Schlögl's model and were studied in [18,19,29,30] on a cube with Neumann boundary conditions, and in [4,5] on a periodic lattice. Recently, perturbations of the voter models which contain the Lotka-Volterra systems are considered in [14]. In addition to results on hydrodynamic limit, [14] also established general conditions for the existence of non-trivial stationary measures and for extinction of the particles.
Our model is a non-conservative stochastic particle system which consists of two types of particles. In [9], Burdzy and Quastel studied an annihilating-branching system of two types of particles, for which the total number of particles of each type remains constant over the time. Its hydrodynamic limit is described by a linear heat equation with zero average temperature. In contrast, besides being non-conservative, our model gives rise to a system of nonlinear differential equations that seems to be new. Moreover, the interaction between two types of particles is singular near the interface of the two media, which gives rise to a boundary integral term in the hydrodynamic limit. The approach of this paper provides some new tools that are potentially useful for the study of other non-equilibrium systems.
We now give some more details on the discrete approximation of the spatial motions in our modeling. We approximate D ± by square lattices D ε ± of side length ε, and then approximate reflected diffusions on D ± by continuous time random walks (CTRWs) on D ε ± . The rigorous formulation of the particle system is captured by the operator L ε in (2.8).
Let X ± i (t) be the position of the particle with index i in D ± at time t. We prescribe each particle a mass 1/N and consider the normalized empirical measures  Here 1 y (dx) stands for the Dirac measure concentrated at the point y, while α ∼ t if and only if the particle X + α is alive at time t, and β ∼ t if and only if the particle X − β is alive at time t. For fixed positive integer N and t > 0, X N,± t is a random measure on D ± . We want to study the asymptotic behavior, when N → ∞ (or equivalently ε → 0), of the evolution in time t of the pair (X N,+ t , X N,− t ).

Main results
Our first main result (Theorem 2.20) implies the following. Suppose each particle in D ± is approximating a RBM with gradient drift 1 2 ∇(log ρ ± ), where ρ ± ∈ C 1 (D ± ) is strictly positive. Then under appropriate assumptions on the initial configuration (X N,+ 0 , X N,− 0 ), the normalized empirical measure (X N,+ t (dx), X N,− t (dy)) converges in distribution to a deterministic measure (u + (t, x)ρ + (x)dx, u − (t, y)ρ − (y)dy) for all t > 0, where (u + , u − ) is the solution of the following coupled heat equations: and where n ± is the inward unit normal vector field on ∂D ± of D ± and 1 I is the indicator function on I.
Note that ρ ± = 1 corresponds to the particular case when there is no drift.
The above result tells us that for any fixed time t > 0, the probability distribution of a randomly picked particle in D ε ± at time t is close to c ± (t)u ± (t, x) when N is large, where c ± (t) = ( D± u ± (t)) −1 is a normalizing constant. In fact, the above convergence holds at the level of the path space. That is, the full trajectory (and hence the joint law at different times) of the particle profile converges to the deterministic scaling limit described by (1.1) and (1.2), not only its distribution at a given time. Question: How about the limiting joint distribution of more than one particles ?
Our second main result (Theorem 2.22) answers this question. It asserts that propagation of chaos holds true for our system; that is, when the number of particles tends to infinity, their positions appear to be independent of each other. More precisely, suppose n and m unlabeled particles in D ε + and D ε − , respectively, are chosen uniformly among the living particles at time t. Then, as N → ∞, the probability joint density function for their positions converges to c (n,m) (t)

Key ideas
A key step in our proof of propagation of chaos (Theorem 2.22) is Theorem 3.7. The latter establishes uniqueness of solution for the infinite system of equations satisfied by the correlation functions of the particles in the limit N → ∞. Such infinite system of equations is sometimes called BBGKY-hierarchy in statistical physics. Our BBGKY hierarchy involves boundary terms on the interface, which is new to the literature. Our proof of uniqueness involves a representation and manipulations of the hierarchy in terms of trees. This technique is related to but different from that in [21] which used Feynman diagrams. It is potentially useful in the study of other stochastic models involving coupled differential equations.
To establish hydrodynamic limit result (Theorem 2.20), we employ the classical tightness plus finite dimensional distribution approach. Tightness of (X N,+ , X N,− ) in the Skorokhod space is proved in Theorem 4.4. This together with the propagation of chaos result (Theorem 2.22) establishes the hydrodynamic limit of the interacting random walks.
Two new tools for discrete approximation of random walks in domains are developed in this article. Namely, the local central limit theorem (local CLT) for reflected random walk on bounded Lipschitz domains (Theorem 2.12) and the 'discrete surface measure' (Lemma 2.4). We believe these tools are potentially useful in many discrete schemes which involve reflected Brownian motions.
Weak convergence of simple random walk on D ε ± to RBM has been established for general bounded domains in [6] and [7]. However, we need more for our model; namely a local convergence result which guarantees that the convergence rate is uniform up to the boundary. For this, we establish the local CLT. We further generalize the weak convergence result and the local limit theorem to deal with RBMs with gradient drift. There are two reasons for us to consider gradient drift. First, it is physically natural to assume the particles are subject to an electric potential. Second, the maximal extension theorem, [13, Theorem 6.6.9], which is a crucial technical tool used in [6] and [7], has established only in symmetric setting. The proof of the local CLT is based on a 'discrete relative isoperimetric inequality' (Theorem 5.5) which leads to the Poincaré inequality and the Nash inequality. The crucial point is that these two inequalities are uniform in ε (scaling of lattice size) and is invariant under the dilation of the domain D → aD.
The paper is organized as follows. In section 2, we introduce the stochastic model and some preliminary facts that will be used later. We then prove the existence and uniqueness of solution for the coupled PDE. The main results, Theorem 2.20 and Theorem 2.22, will be rigorously formulated. We also mention various extensions of our main results in Remark 2.25. Section 3 and section 4 contains the proof of Theorem 2.22 and Theorem 2.20 respectively. Section 5 is devoted to the proofs of the discrete relative isoperimetric inequality and the local CLT.

Notations and Preliminaries
For the reader's convenience, we list our notations here: filtration induced by the process (X t ), i.e. F X t = σ(X s , s ≤ t) 1 x indicator function at x or the Dirac measure at x (depending on the context) L −→ weak convergence of random variables (or processes) We also use the following abbreviations: given functions which are strictly positive, λ > 0 is a fixed parameter.

Interacting random walks in domains
In this subsection, we describe the interacting random walk model. We start with some key ingredients needed in discrete approximation.

Discrete approximation of surface measure
To capture the boundary behavior of our processes near the interface I in the discrete scheme, we need a discrete approximation of the surface measure σ on I. The construction of I ε and σ ε in the following lemma is a key to our approximation scheme. For us, N := {0, 1, 2, · · · } denotes the set of non-negative integers.
Then there exist finite subsets I (j) = I εj of I and functions σ (j) = σ εj : such that (a) and (b) below hold simultaneously: where #A denotes the number of elements in the finite set A, B(x, s) = {y ∈ R d : |y − x| < s} is the ball with radius s centered at x, and C is a constant that depends only on D.
(b) For any equi-continuous and uniformly bounded family F ⊂ C(I), Proof We can always split I into small pieces. The point is to guarantee that each piece is not too small, so that σ (j) /ε d−1 ≥ C and that (2.1) holds. Since I is H d−1 -rectifiable, we have for R ∈ (0, 1], where C does not depend on R. Since I is closed, it is regular with dimension d − 1 in the terminology of section 1 of [16]. Hence by [15] or section 2 of [16], we can build "dyadic cubes" for I. More precisely, there exists a family of partitions {∆ j } j∈Z of I into "cubes" Q such that for all Q ∈ ∆ j and r > 0. Here the constant C is independent of j, Q, or r. Note that H d−1 is the surface measure σ of ∂D and that property (iii) tells us that the cubes have relatively small boundary. In particular, (iii) implies σ(∂Q ∩ I) = 0 for all cube Q.
. We pick one point z i ) and the Lipschitz property of ∂D. It remains to check (2.2).
Remark 2.5. (2.2) implies that we have the weak convergence z∈I (j) σ (j) δ z → σ I on the space M + (I) of positive finite measure Borel measures on I. Here δ z is the dirac delta measure at z, and σ I is the surface measure restricted to I. (2.1) is a control on the number of points locally in I εj . We call I ε the 'ε-point approximation' of I and σ ε the 'discrete surface measure' associated to I ε .
Remark 2.6. The above lemma remains true if I is the finite union of disjoint closed connected and H d−1 -rectifiable subsets of ∂D. This enables us to deal with disconnected interface I.

Reflected diffusion and random walk approximation
We now describe the motion of each underlying particle. First we fix a bounded Lipschitz domain D ⊂ R d and any ε > 0. Without loss of generality, we assume that the origin 0 ∈ D. Let εZ d be the union of all closed line segments joining nearest neighbors in εZ d , and (D ε ) * the connected component of D ∩ εZ d that contains the point 0. Set D ε = (D ε ) * ∩ εZ d . We can view D ε as the vertices of a graph whose edges coming from (D ε ) * . We also denote the graph-boundary Since D is Lipschitz, (W 1,2 (D), E) is a regular Dirichlet form on L 2 (D; ρ) and so there is a ρ-symmetric diffusion X associated with it (cf. for example [12]).
Definition 2.7. We call X the (I d×d , ρ)-reflected diffusion, where I d×d is the d × d identity matrix. When ρ = 1, X is called the reflected Brownian motion (RBM) in D. Hence a (I d×d , ρ)-reflected diffusion is a RBM inD with drift 1 2 ∇ log ρ.
The L 2 -infinitesimal generator of X is Moreover, X has the Skorokhod representation: where n is the inward unit normal of ∂D and L is the positive continuous additive function (PCAF) of X whose Revuz measure is 1 2 σ (c.f. [12]). We call L the boundary local time of X. Next, we define X ε to be a continuous time random walk (CTRW) on D ε with exponential waiting time of parameter d ε 2 and one step transition probabilities where {µ xy : x, y ∈ D ε } are symmetric weights (conductances) to be constructed in two steps as follows: First, for every x ∈ D ε \ ∂D ε and i = 1, 2, · · · , d, define Clearly, µ xy = µ yx for all x, y ∈ D ε \ ∂D ε . Note that since ρ is C 1 and strictly positive on D, when ε is sufficiently small, µ x,x+ε ei and µ x,x−ε ei are strictly positive for every x ∈ D ε \ ∂D ε and i = 1, 2, · · · , d. Second, we define , if x, y ∈ ∂D ε are adjacent in D ε Now µ xy = µ yx for all x, y ∈ D ε . A heuristic reason of the above construction can be found in [23].
We call X ε the ε-approximation of X. Clearly, X ε is symmetric with respect to the measure m ε defined by Since ρ ∈ C 1 (D), there exists a constant C > 0 such that A special but important case is when ρ ≡ 1. In this case, X is simply the reflected Brownian motion on D, and X ε is a simple random walk on the graph D ε . It is proved in [6] that X ε converges weakly to the reflected Brownian motion X as ε → 0.

Random walks with interaction
Fix ε = ε j = 2 −j (j ∈ N) and N = 2 jd such that N ε d = 1. Assume there are N "+" particles in D ε + and N "−" particles in D ε − at t = 0. Each particle moves as an independent CTRW X ε,± (see the previous subsection) in its respective domain D ε ± . Let I ε be the finite subset of I defined in Lemma 2.4. For each z ∈ I ε , pick an z + ∈ D ε + and an z − ∈ D ε − which are closest to z (See Figure 2). A pair of particles of opposite charges at (z + , z − ) is being killed with a certain rate to be explained. Note that for ε small enough, we have sup z∈I ε |z ± − z| ≤ 2M ε, where M is the Lipschitz constant of I.
The state space of the particle system is the collection of configurations The state of the particle system at time t is a random element η ε t = (η ε,+ t , η ε,− t ) ∈ E ε . Here η ε,± t (x) stands for the number of "± particles at x ∈ D ε ± at time t. We omit ε and N for convenience when there is no ambiguity. For example, we write η t and m(x) in place of η ε t and m ε (x) respectively. The function ξ such that ξ(x) = 1 and ξ(y) = 0 for y = x is denoted as 1 x . Definition 2.8. η t is defined to be the unique strong Markov process which has the generator L = L ε given by where L ε 0 is the generator of two families of independent random walks in D ε + and D ε − , respectively, with no annihilation between them, namely and K ε is the operator corresponding to annihilation between particles of opposite signs at the interface I ε , namely where p ± xy is the one-step transition probabilities for the CTRW X ε,± on D ε ± (without any interaction) and with σ ε and I ε being constructed by Lemma 2.4.
The expression for K ε comes from the underlying assumptions of the model: First, the term η + (z + )η − (z − ) is combinatorial in nature. Since there are η + (z + )η − (z − ) pairs of particles at position (z + , z − ), the chance of killing is proportional to the number of ways of selecting a pair of particles near the interface. Second, each pair of particles near I disappears at rate (λ/ε) Ψ ε (z) where λ is a parameter. Intuitively, in the limit, the amount of annihilation in a neighborhood of a point is proportional to the surface area of the interface I in that neighborhood. The scaling 1/ε is suggested by the observation that there are about 1/ε "layers" starting from the interface I, so that the chance for a particle to arrive near I is of order ε. Ψ ε (z) is comparable to 1 and can be viewed as a normalizing constant with respect to the lattice. This choice (2.11) is justified in the proof of Theorem 3.5.
On other hand, suppose g ∈ B b ([0, T ] × ∂D). Then for t ∈ [0, T ] and x ∈ D, g(s, y)p(s, x, y)σ(dy)ds. (2.15) Now let X ε be the ε-approximation of X with symmetrizing measure m ε . The transition density p ε of X ε with respect to the measure m ε is Clearly, p ε is strictly positive and is symmetric in x and y.
We will prove in Section 5 that p ε enjoys two-sided Gaussian bound and is jointly Hölder continuous uniform in ε ∈ (0, ε 0 ) for some ε 0 > 0, and that p ε converges to p uniformly on compact subsets of (0, ∞) × D × D. In rigorous terms, we have the following three results. The important point is that the constants involved in these results do not depend on . Theorem 2.9. (Gaussian upper bound) There exist C k = C k (d, D, ρ, T ) > 0, k = 1, 2, and ε 0 = ε 0 (d, D, ρ, T ) ∈ (0, 1] such that for every ε ∈ (0, ε 0 ) and x, y ∈ D ε ,
Proof Fix θ ∈ [0, T ]. By the Gaussian upper bound in Theorem 2.9, we have y∈∂D ε p ε (θ, x, y) Here C i are all constants which depend only on d, D and T . Note that in the second last line, we used the fact that |{∂D ε ∩ B(x, s)}| ≤ C((s/ε) d−1 ∨ 1), which follows from Lemma 2.4. The proof is now complete.
In general, we use " ± " for quantities related to X ± . For example, A ± , (P ± t ) t≥0 and p ± (t, x, y) denote the generator, semigroup and transition density of the reflected diffusion X ± in D ± with respect to ρ ± (x)dx. In addition, we use "ε" for quantities related to the CTRWs in the discrete domains D ε ± . For example, p ε,± (t, x, y) denotes the transition density of the CTRW X ε,± on D ε ± with respect to the measure m ± ε . We also denote p ε (t, ( r, s), ( r , s )) : By applying Lemma 2.13 to p ε,± (t, x, y), then by the boundedness of Φ ε in (2.11), Theorem 2.12, Lemma 2.4 and LDCT, we have the following approximation for the local time of X ± on I.

Hydrodynamic limit: system of nonlinear PDEs
In this subsection, we provide suitable notion of solutions for the coupled PDE (1.1) and (1.2), and then prove the existence and uniqueness of the solution. Throughout this subsection, D is a bounded Lipschitz are the semigroup and the transition density of X, with respect to the measure ρ(x)dx.
Observe that (1.1) is a second order parabolic equation for u + with Robin boundary condition, and (1.2) is a similar equation for u − . This leads us to consider the following Robin boundary problem, where (2.23) By Itó's formula and the Skorokhod representation for X, we see that a classical solution of (2.23), should it exists, has the probabilistic representation where L is the boundary local time of X.
Definition 2.15. u defined by (2.24) is called a probabilistic solution of (2.23).
First we show that the function u defined by (2.24) is continuous.
Proof Observe that for any r ∈ [0, t], (2.25) By Markov property, the first term is Since X has the strong Feller property (see [1]) and u is bounded, The second term of (2.25) converges to 0 uniformly on (0, T ]×D, as r → 0. This is because its absolute value is bounded by Hence, u is continuous in x ∈ D. By a similar calculation as in (2.25), we have Hence, Both terms go to 0 uniformly in x ∈ D as a goes to 0. (In fact, the first term goes to 0 uniformly since the semigroup P t is strongly continuous on C(D). For the second term, a 0 ∂D p(s, x, z) σ(dz)ds ≤ 2C 1 √ a + C 2 a also goes to 0 uniformly in x.) Hence u is continuous in t ∈ (0, T ] uniformly in x ∈ D. Therefore, u ∈ C((0, T ] × D). If ϕ ∈ C(D), we can extend the above argument to show that u ∈ C([0, T ] × D).
Remark 2.17. In fact, one can allow g to be unbounded and show that the conclusion of Lemma 2.16 remains true if gσ satisfies a Kato class condition: Then that satisfies the following integral equation: . Moreover, by Markov property and (2.15) we have Hence u satisfies the integral equation. It remains to prove uniqueness. Supposeũ ∈ C([0, T ] × D) also satisfies the integral equation.
By a Gronwall type argument and (2.14), we can show that w = 0. More precisely, let ψ(s) = sup x∈D |w(s, x)|. Then Combining the above two inequalities, we have Integrating both sides with respect to t on the interval [0, t 0 ], we have and repeat the argument to obtainψ = 0 on [0, Now we come to our coupled equation.
Then there is a unique element (u + , u − ) ∈ Λ T which satisfies the coupled integral equation (2.28) where L ± is the boundary local time of X ± on the interface I.
Functions (u + , u − ) satisfying equation (2.28) will be called a weak solution of (1.1) and (1.2), as it can be shown that they are weakly differentiable and solve the equations in the distributional sense.
A similar result holds for S − u 1 − S − u 2 . Hence, for some γ < 1 when T is small enough.
Hence there is a T 0 > 0 such that S : Λ T0 → Λ T0 is a contraction map. By Banach fixed point theorem, there is a unique element (u , v ) ∈ Λ T0 such that (u , v ) = S(u , v ). By Proposition 2.18, (u , v ) is the unique weak solution to the coupled PDE on [0, T 0 ].
Repeat the above argument, with u 0 (·) replaced by u (T 0 , ·), and v 0 (·) replaced by v (T 0 , ·). We see are the same, we can extend the solution of the coupled PDE uniquely to [T 0 , 2T 0 ]. Iterating the argument, we have for any T > 0, the coupled PDE has a unique weak solution in Λ T . Invoke Proposition 2.18 once more, we obtain the desired implicit probabilistic representation (2.29).
Finally, by using Markov property as in the proof of Proposition 2.18, we see that (2.29) and (2.28) are equivalent.

Main results (rigorous statements)
In this paper, we always assume the scaling N ε d = 1 holds, so that the interacting random walk model is parameterized by a single parameter N which is the initial number of particles in each of D ε + and D ε − . More precisely, for each fixed N , we set ε = N −1/d and let (η ε t ) t≥0 be a Markov process having generator L ε defined in (2.8) and having initial distribution satisfies x∈D ε where M ≤1 (E) denotes the space of non-negative Borel measures on E with mass at most 1. M ≤1 (E) is a closed subset of M + (E), where the latter denotes the space of finite non-negative Borel measures on E equipped with the following metric: where {φ k : k ≥ 1} is any countable dense subset of C(E). The topology induced by this metric is equivalent to the weak topology (i.e. µ n − µ → 0 if and only if µ n , f → µ, f for all f ∈ C(E)). Under this metric, M + (D) is a complete separable metric space, hence so are E and the Skorokhod space D([0, T ], E) (see e.g. Theorem 3.5.6 of [22]). Here is our first main result.
Theorem 2.20. (Hydrodynamic Limit) Suppose Assumption 2.3 holds and the sequence of initial configurations η ε 0 satisfies the following conditions: has unit mass for all N .) Then for any T > 0, as ε → 0 along the sequence ε j = 2 −j , we have Theorem 2.20 gives the limiting probability distribution of one particle randomly picked in D ε ± at time t. This is the 1-particle distribution in the terminology of statistical physics.
Question: What is the limiting joint distribution of more than one particles?
Before stating the answer, we need to introduce a standard tool in the study of stochastic particle systems: the notion of correlation functions 1 . Recall that the state space of η ε = (η ε t ) t≥0 is E ε defined in (2.7). We denote by the set of configurations with n and m particles in D ε + and D ε − respectively. We then define A : E ε ×E ε → R in such a way that whenever ξ ∈ Ω ε n,m , where n → A n k is the Poisson polynomial 2 of order k, namely A n 0 := 1 and A n k := n(n − 1) · · · (n − k + 1) for k ≥ 1 (in particular, A n k = 0 for k > n). Note that A n k is the number of permutations of k objects chosen from n distinct objects. So A(ξ, η) is the total number of possible site to site pairings between labeled particles having configuration ξ with a subset of labeled particles having configuration η. An alternative representation of (2.31) will be given in (3.5). Convention: Definition 2.21. Let P η is the law of a process with generator L ε and initial distribution η satisfying when ξ = ( i δ ri , j δ sj ). By convention, we also have γ ε (( r, s), t) := γ ε (ξ, t). Note that γ ε depends on the initial configuration of η.
Intuitively, suppose we randomly pick n and m living particles in D + and D − respectively at time t, then ( r, s) → γ ε,(n,m) (( r, s), t) is the joint probability density function for their positions, up to a normalizing constant. Therefore, it is natural that γ ε,(n,m) defined by (2.33) is called the (n, m)-particle correlation function.
The next is our second main result, propagation of chaos, for our system. It says that when the number of particles tends to infinity, they appears to be independent of each other. Mathematically, the correlation function factors out in the limit N → ∞.
To investigate the intensity of killing near the interface, we define J N,± ∈ D([0, ∞), M + (D ± )) by 1 , which measures the number of encounters of the two types of particles near I. An immediately corollary of Theorem 2.22 is the following, which is what we need to identify the limit of (X N,+ , X N,− ).
Proof We only need to apply Theorem 2.22 for the cases (n, m) = (1, 1) and (n, m) = (1, 0). By definition, Using (2.6) and Lemma 2.4, we get the first two equations via Theorem 2.22. Using (2.6) and the assumption that ρ ± ∈ C(D ± ), we have the last two equations again by Theorem 2.22.
Remark 2.24. (Conditions on η 0 ) The two conditions for the initial configuration η 0 in Theorem 2.20 are mild and natural. They are satisfied, for example, when each particle has the same random initial distribution . Condition (ii) guarantees that, asymptotically, there is no "blow up" of number of particles at any site. More precisely, this technical condition is imposed so that we have The above can be easily checked by comparing with the process η that has no annihilation (i.e. η has generator L ε 0 ). Alternatively, we can use the comparison result (3.11) to prove (2.37). Remark 2.25. (Generalization) We can generalize our results in a number of ways by the same argument. For example, the initial number of particles in D + and D − can be different, the condition N ε d = 1 can be relaxed to lim N →∞ N ε d → 1 where ε depends on N . The annihilation constant λ can be replaced by a space and time dependent function λ(t, x) ∈ C([0, ∞) × I). The diffusion coefficients in D + and D − can be different. The condition "X N,± 0 has mass one for all N " can be replaced by "the mass of X N,± 0 is uniformly bounded in N ". More generally, the same method can be extended to deal with similar models with more than two types of particles.
The remaining part of this paper is devoted to the proof of Theorem 2.20 and Theorem 2.22. We first prove Theorem 2.22 because the proof of Theorem 2.20 relies on Theorem 2.22.

Duality
The starting point of our analysis is the discrete integral equation for γ ε in Lemma 3.2. At the heart of its proof is the dual relation in Lemma 3.1, which says that the two independent processes ξ 0 = (ξ 0 t ) t≥0 and η 0 = (η 0 t ) t≥0 of independent ransom walks with no interaction are dual to each other with respect to the function A(ξ,η) αε(ξ) , where ξ, η ∈ E ε . Such kind of dual formula for the whole grid Z d appeared in [5] and in Chapter 15 of [11].
∈ Ω ε n,m for any (n, m). We first focus on D + in Step 1 and Step 2 below.

Putting this into (3.3), we have
where E is the expectation corresponding the probability measure under which each coordinate processes of { Y +,ε t ; t ≥ 0} are independent CTRWs with Y +,ε 0 = r and are independent of (η 0 t ) t≥0 .
Denote by P η 0 the law of the unlabeled process (η 0 t ) t≥0 starting from η ∈ Ω ε N,M and has generator L ε 0 . Since all processes on D ε + are independent of those on D ε − , we can proceed as in step 2 (via (3.5)) to obtain , s), and its components are mutually independent CTRWs on D ε ± , respectively. The proof is now complete by taking θ = 0. We now formulate the discrete integral equations that we need. Recall the definition of K ε from (2.10) and the definition of γ ε (( r, s), t) from (2.33).
Proof Starting from (3.1), we can obtain Lemma 3.2 by 'integration by parts' as follows.
Let P (ξ 0 ) and P (η 0 ) be the laws of ξ 0 and η 0 respectively. (3.1) is equivalent to saying that for any ξ and η, we have Taking w = t − s, we see that (3.8) is in turn equivalent to where P (ξ 0 ) t and P (η 0 ) t are the transition semigroup of ξ 0 and η 0 , respectively, and they act on the ξ and η variables in A(ξ,η) αε(ξ) , respectively. Therefore, with L ε 0 acting on the η variable, we have Recall that η t is the configuration process of our interacting system with generator L ε 0 + K ε (see Definition 2.8). Fix ξ and consider the function (s, η) → F s (η) := F (ξ) s (η). We have is a F η s -martingale for s ∈ [0, t]. By (3.10) and the fact that E η [M t ] = E η [M 0 ] = 0, where P η is the law of (η t ) t≥0 starting from η, we have for all ξ and η. This is equivalent to the stated equation in the lemma.
It is clear that K ε A(ξ, η) ≤ 0. Hence, as an immediate consequence of Lemma 3.2, we have the following comparison result: for all t > 0 and ( r, s)
Observe also that for x ∈ D ε + and y ∈ D ε − , Similarly, From the above calculations and (3.12), we see that for ξ ∈ Ξ, which gives the desired formula.

Uniform bound and equi-continuity
We extend to define γ ε,(n,m) (·, t) continuously on D n + × D m − while preserving the supremum and the infinmum in each small ε-cube. We can accomplish this by the interpolation described in [2] or [39], or by a sequence of harmonic extensions along simplexes with increasing dimensions (described in [23]). Recall that the definition of γ ε,(n,m) (·, t) depends on the initial configuration η 0 of the interacting random walks (see Definition 2.21), which has the normalization (2.32). Proof We first prove uniform boundedness. By (3.11) and the Gaussian upper bound in Theorem 2.9, we have γ ε (( r, s), t) ≤ ( r , s ) γ ε (( r , s ), 0) p ε (t, ( r, s), ( r , s )) m( r , s ) Since the initial distribution η 0 = (η + 0 , η − 0 ) has the property that x∈D ε Thus there exist ε 0 = ε 0 (d, D, ρ) and C = C(d, D, ρ) > 0 such that for all t ∈ (0, ∞), (n, m) ∈ N × N and ε ∈ (0, ε 0 ), (3.17) We next show that both terms on the right hand side of (3.7) are equi-continuous. Recall that we can rewrite the equation (3.7) as where we have used (3.16) in the last line. By the uniform Hölder continuity of p ε (t, ( r, s), ( r , s )) (see Theorem 5.12 below) and the fact that p ± (t, x, y) ∈ C((0, ∞) × D ± × D ± ), we see that {F ε } is equicontinuous at (( r, s), t). For the second term, note that where and In the remaining, we will show that G ε is equi-continuous. We first deal with H (1) in (3.18).
As in (3.16), we have On other hand, using (3.11), the Chapman Kolmogorov equation and assumption (ii) for η 0 , in this order, we have sup θ1, θ2>0 for large enough N , where C > 0 is a constant.
From Theorem 3.4 and a diagonal selection argument, it follows that for any sequence ε k → 0 there is a subsequence along which γ ε converges on D n + × D m − × (0, T ), uniformly on the compacts, to some γ (n,m) ∈ C(D n + × D m − × (0, T )), for every (n, m) ∈ N × N. Our goal is to show that We will achieve this by first showing that both Γ = {γ (n,m) } and Φ

Limiting hierarchy
Note that D n + × D m − is a bounded Lipschitz domain in R (n+m)d , and that the boundary ∂(D n We define the function ρ = ρ (n,m) : where p ± is the transition density of the reflected diffusion X ± on D ± with respect to the measure ρ ± (x)dx. We now characterize the subsequential limits of {γ ε } ε>0 : Theorem 3.5. Let η ε 0 be a sequence of initial configurations that satisfy (2.32) with ε = N −1/d , and the conditions of Theorem 2.20; that is, their corresponding empirical measures X N,± 0 converges weakly to u ± 0 (z)dz in M ≤1 (D ± ) for some u ± 0 ∈ C(D ± ) and Denote by Γ ε = {γ ε, (n,m) ; t ≥ 0, n, m ∈ N} the correlation functions for the interacting random walks with initial configuration η ε 0 . Let Γ = {γ (n,m) t ; t ≥ 0, n, m ∈ N} be any subsequential limit (as ε → 0) of Γ ε = {γ ε, (n,m) ; t ≥ 0, n, m ∈ N}. Then the following infinite system of hierarchical equations holds: Observe that the coordinate processes of X (n,m) consist of n independent copies of reflected diffusions in D + and m independent copies of reflected diffusions in D − .
(iii) It is easy to check by using (ii) and Proposition 2.18 that is a solution of (3. γ ε (( a, b), 0) p ε (t, ( r, s), ( a, b)) m( a, b) Fix any (n, m) ∈ N × N, t > 0 and ( r, s) ∈ (D ε + ) n × (D ε − ) m . By a simple counting argument and condition (3.24) for η ε 0 , we see that the first term in (3.27) equals which converges to E ( r, s) [Φ (n,m) (X (n,m) (t))] by Theorem 2.12 and assumption (i) for the initial distributions. Here P ( r, s) is the probability measure for X (n,m) starting at ( r, s).
For the third term in (3.27), we split the integrand with respect to dθ into three terms corresponding to (3.13), (3.14) and (3.15) respectively. The term corresponding to (3.13) equals By Theorem 2.12 and Lemma 2.4, and the convergence is uniform for r i ∈ D ε + . Therefore, by applying Theorem 2.12 again, the term corresponding to (3.13) We repeat the same argument for the term corresponding to (3.14). Moreover, note that the term corresponding to (3.15) will not contribute to the limit as ε → 0, by the same argument we used for the second term in (3.27). Therefore, the integrand of the second term in (3.27) The integral for θ ∈ (0, t) in the third term in (3.27) then converges to the desired quantity, by (3.21) and LDCT. The proof is complete.
In view of Remark 3.6(iii), the proof of Theorem 2.22 (Propagation of Chaos) will be complete once we establish the uniqueness of the solution of the limiting hierarchy (3.25). This will be accomplished in Theorem 3.7 in the next subsection.

Uniqueness of infinite hierarchy
Uniqueness of BBGKY hierarchy is an important issue in statistical physics. For example, it is a key step in the derivation of the cubic non-linear Schrödinger equation from the quantum dynamics of many body systems obtained in [21]. Our BBGKY hierarchy (3.25) is new to the literature and the proof of its uniqueness involves a representation and manipulations of the hierarchy in terms of trees. The technique is related but different from that in [21], which used the Feynman diagrams.
Note that, by Theorem 3.5, γ The remaining of this subsection is devoted to give a proof of this theorem. Convention in this subsection: β = {β (n,m) } will always denote the functions stated in Theorem 3.7. For notational simplicity, we will also assume λ = 2 and ρ ± = 1. The proof for the general case is the same. we will also drop T from the notation β (n,m) t (T,n,m) . It is convenient to rewrite the infinite hierarchy (3.25) in a more compact form as Here σ (n,m)

Duhamel tree expansion
We now describe the infinite hierarchy in detail. It is natural and illustrative to represent the infinite hierarchy in terms of a tree structure, with the 'root' at the top and the 'leaves' at the bottom. Fix two positive integers n and m. We construct a sequence of finite trees {T (n,m) N : N = 0, 1, 2, · · · } recursively as follows.

T
For mnemonic reason, we use the same notation T (n,m) N to denote the collection of paths that connect the root to a leaf in T (n,m) N . By induction, the total number of paths (or the total number of leaves) is where V θi (for i = 1, 2, · · · , N ) is defined by V +i = V + i and V −j = V − j . For example, when (n, m) = (2, 5), N = 3 and the path is θ = (+ 1 , − 6 , − 5 ). Then

Telescoping via Chapman-Kolmogorov equation
By (3.31), the right hand side of (3.32) is a sum of (n + m)(n + m + 1) · · · (n + m + N − 1) terms of multiple integrals. We will apply the bound β (p,q) t (p,q) ≤ C p+q to each term, and then simplify the integrand using Chapman-Kolmogorov equation.
We demonstrate this for the twelve terms for the case (n, m, N ) = (1, 2, 2). The twelve terms on the right hand side of (3.32) are . (3.33) The first four terms came from the leftmost leaf of the previous level, we group them together to obtain, for (x, y 1 , Note the telescoping effect upon using the Chapman-Kolmogorov equation for the middle two terms in the last equality above gives rise to t − t 3 rather than t 2 − t 3 . We apply (2.14) to obtain where C ± = C(D ± , T ) are positive constants. Repeat the above argument for the remaining eight terms of (3.33), we obtain The key is to visualize the twelve terms on the right as 12 paths of T (1,2) 2 with the edges relabeled. We denote this relabeled tree by S , we only need to record the indexes (p, q) and the ± sign. For example, the first four terms can be represented by and hence it appears twice. In S (1,2) 2 , these four paths are formed by a + 1 edge followed by four edges with labels In general, we obtain S  2. At level k ≥ 2, consider the set Λ + := {+ 1 , · · · , + 1 , + 2 , + 3 , · · · , + k } in which we have n copies of + 1 (hence there are n + k − 1 elements in Λ + , in which n of them are + 1 ). Similarly, let k−1 , let R ξ be the labels of (the edges of) the path from the root to ξ in S (n,m) k−1 , counting with multiplicity. Finally, the collection of new labels of the edges below ξ, denoted by L ξ , is chosen in such a way that Λ + ∪ Λ − = R ξ ∪ L ξ (counting multiplicity).
For mnemonic reason, we use the same notation S (n,m) N to denote the collection of paths from the root to the leaves of S (n,m) N . Any such path is represented by the ordered (new) labels of the edges. We now 'forget' the sign of the edges and only record the integer labels. For example, the path (− 1 , + 1 , − 2 , + 3 ) is replaced by (1, 1, 2, 3).
Using the hypothesis β (n,m) t (n,m) ≤ C n+m (of Theorem 3.7) and applying Chapman-Kolmogorov equation to (3.32), and then applying (2.14), we obtain the following lemma by the same argument that we used to obtain (3.34).
Our goal is to show that I (n,m) N (t) ≤ (C t) N/2 for some C = C(n, m) > 0. This will imply β (n,m) t (n,m) = 0 for t > 0 small enough. Clearly we have .
Unfortunately, this crude bound is asymptotically larger than (Ct) N/2 for any C > 0.

Comparison with a 'dominating' tree
Note that This is obtained by comparing the labels in S (1,2) 3 with a 'dominating' labeling, in which the labels of the edges below every leaf of S (1,2) 2 is {3, 3, 2, 2, 1} (the ± sign is discarded). This trick enables us to group the terms at each level.
We can now group the terms in each level k as a sum of k terms to obtain In the last term, we have used the observation that when k > n + m, the smallest element inL ξ is at least 2 and so the sum stops before reaching 1/ √ t − t k+1 . From this and the simple estimates like we have derived the following Lemma 3.9. (3.38)

Estimating J N
Our goal in this section is show that J N (t) ≤ (C t) N/2 for some C > 0. Our proof relies on the following recursion formula pointed out to us by David Speyer: We assume (3.39) for now and use it to establish the following lemma. The proof of (3.39) will be given immediately after it. Moreover, Proof J N (t) = J N · t N/2 is obvious from (3.39) after a change of variable. Let M N be the collection of functions f : {2, 3, · · · , N + 1} → {1, 2, · · · , N } satisfying f (i) < i. We can rewrite J N (t 1 ) as This is a sum of N ! terms. When we put t 1 = 1, the smallest term is Hence we have the lower bound 2 N ≤ J N . Unfortunately, the largest term is exactly . which grows faster than C N for any C ∈ (0, ∞). Hence for the upper bound, we will employ the recursion formula (3.39). We apply the homogeneity to the right hand side of (3.39) to obtain The integrals are now simple one dimensional and can be evaluated: where W is the Lambert W function. By Lagrange inversion Theorem (see Theorem 5.4.2 of [37]), W (z) = ∞ k=1 (−k) k−1 z k /k! (for |z| < 1/e). Hence by comparing coefficients in the series expansion of φ(x), we have K N = (N +1) N π N/2 (N +1)! as desired.
Remark 3.11. By Stirling's formula, (N +1) b(N ) = 1). Hence J N ≤ C N for some C > 0. Monte Carlo simulations suggests that J N ∼ π N . The recursion (3.43) also makes it clear that J N 's are all in Q[π] (polynomials in π with rational coefficients) and makes it easy to compute them recursively. This is because √ π Γ((nj +1)/2) Γ((nj +2)/2) is rational if n j is odd and is a rational multiple of π if n j is even. For example, J 1 = 2, J 2 = 2 + π and J 3 = 4 + 10π 3 .
We now turn to the proof the recursion formula (3.39) which is restated in the following lemma.
Proof The proof is based on standard combinatorial methods for working on sums over planar rooted trees (see [37]).
Step 1: Summing over labeled trees. Recall (3.42). There are N ! elements in M N . We can visualize each of them as a rooted tree with vertex set {1, 2, · · · , N + 1} and a directed edge from i to f (i) for each i. For example, the 6 elements of M 3 can be represented by The trees are drawn so that arrows point upwards and the children of a given vertex are listed from left to right. Note that the second and forth tree of the list are the same up to relabeling the vertices. The idea is to group terms in (3.42) like this together. First, we rewrite (3.42) in terms of trees. Let D N be the set of 'decreasing trees', which are trees whose vertices are labeled by {1, 2, · · · , N + 1} and such that i < j whenever there is an edge i ← j. Then (3.45) Step 2: Summing over unlabeled trees. A planar tree is a rooted unlabeled tree where, for each vertex, the children of that vertex are ordered. We draw a planar tree so that its children are ordered from left to right. Here are the 5 planar rooted trees on 3+1 vertices: Let T k be the set of planar rooted trees with k vertices. In general, there are (2N )! N !(N +1)! (the Catalan number) elements in T N +1 , see exercise 6.19 in [37]. We now group all the integrals in (3.45) with the same planar tree. For example, two different labeled trees (the second and forth in our list of labeled trees) both give the same unlabeled planar tree (which is the second in the above list). We redraw this unlabeled planar tree T 0 below and attach letters {a, b, c, d} to T 0 for later use.
The integrands corresponding to the second and forth labeled trees are and They are the same to once we relabel the variables by the vertices of T 0 . That is, (1, 2, 3, 4) → (a, b, c, d) for the first term and (1, 2, 3, 4) → (a, b, d, c) for the second.
We now go back and keep track of the bounds of integration. In the first integral, they are t a ≥ t b ≥ t c ≥ t d and, in the second integral, they are t a ≥ t b ≥ t d ≥ t c . We can group these together as which is the same as t a ≥ t b ≥ t c and t b ≥ t d .
In general, the inequality constraints we have are of two types. First, whenever we have an edge u ← v, we get the inequality t u ≥ t v . Second, if v and w are children of u with v to the left of w, then t v ≥ t w . Let P (T, t 1 ) be the polytope cut out by these inequalities where t 1 is the variable at the root. We have proved

46)
Step 3: Grouping terms for which the root has degree k. We abbreviate if T has more than one vertex (otherwise ω(T, t 1 ) := 1). Then (3.46) translates into Fix an integer k and let T be a tree whose root has degree k. Removing the root leaves behind k children, denoted in chronical order by t 2 , · · · t k+1 , and k planar subtrees T j having t j as its root for 2 ≤ j ≤ k + 1. Then Hence, group together the terms where the root has degree k, we have Here the summation conditions include that k+1 j=2 |T j | = N and each different ordering of (T 2 , T 3 , · · · , T k+1 ) are considered to be different, where |T j | is the number of vertices in T j . This abbreviation applies whenever T2, ··· , T k+1 appears.
On other hand, we have by applying (3.47) Step 4: Identifying the integrals. It remains to show that (3.48) is equal to (3.49). Let S k denote the space of permutations of {2, 3, · · · , k + 1}. Then the right hand side of (3.49) is equal to which is J N (t) by (3.48). This completes the proof of the lemma.

Proof of uniqueness
Proof of Theorem 3.7. By Lemma 3.8, Lemma 3.9 and Lemma 3.10, we have for all t ∈ [0, T ] and N ∈ N. This implies that there is a constant τ > 0 so that β (n,m) t (n,m) = 0 for t ≤ τ and for all (n, m) ∈ N × N. Note thatβ t := β τ +t also satisfies the hierarchy (3.30), and that β 0 = 0. Using the hypothesis β

Hydrodynamic Limits
This section is devoted to the proof of Theorem 2.20. Throughout this section, we assume the conditions of Theorem 2.20 hold.

Constructing martingales
Since η t = η ε t has a finite state space, we know that for all bounded function F : R + × E ε → R that is smooth in the first coordinate with sup (s,x) ∂F ∂s (s, x) < C < ∞, we have two F η t -martingales below: and where L = L ε is the generator defined in (2.8). See Lemma 5.1 (p.330) of [27] or Proposition 4.1.7 of [22] for a proof. We will use this fact to construct some important martingales in Lemma 4.1 below.
In general, suppose X = (X t ) t≥0 is a CTRW in a finite state space E, whose one step transition probability is p xy and mean holding time at x is h(x). Its infinitesimal generator of X is the discrete operator The formal adjoint A * of A is defined by It can be easily checked that We denote by A ± ε the generator of the CTRW X ±,ε on D ε ± , respectively, and by A * ,± ε the corresponding formal adjoint. In our case, h(x) = h ε (x) = ε 2 /d for all x. We can check that if f ∈ C 2 (D ± ), then is an F η t -martingale for t ≥ 0, where J N,+ is the measure-valued process defined by (2.35). Moreover, if φ ∈ C 1 (D + ), then there is a constant C > 0 independent of N so that for every T > 0, Similar statements hold for X N,− .
Proof The lemma follows by applying (4.1) and (4.2) to the function We spell out the details here for completeness. Observe that f (η t ) = φ, X N,+ t . Fix x 0 ∈ D ε + , and define η + x0 to be the function from E ε to R which maps η to η + (x 0 ). Then by the definition of L = L ε in (2.8), Similarly, for y 0 ∈ D ε − , we have Hence, by linearity of L, (4.7) and then (4.3), we have and M (t) is an F η t -martingale by (4.1). Next, we compute E[ M t ]. Note that where L(η a η b ) can be computed explicitly using (2.8). Hence from (4.2), we can check that After taking expectation for g(η r ), the first term is the last display is of order at most 1/N since φ ∈ C 1 (D + ), while the second term inside the bracket is of order at most 1/N , uniformly in r ∈ [0, t], by (2.37). Hence E[M 2 (t)] ≤ C N for some C = C(φ, d, D ± , λ). Doob's maximal inequality then gives (4.6).
Remark 4.2. From the second term of (4.9), we see that if the parameter of the killing time is of order λ/ε, then we need N ε d to be comparable to 1.

Tightness
The following simple observation is useful for proving tightness when the transition kernel of the process has a singularity at t = 0. It says that we can break down the analysis of the fluctuation of functionals of a process on [0, T ] into two parts. One part is near t = 0, and the other is away from t = 0 where we have a bound for a higher moment. Its proof, which is based on the Prohorov's theorem, is simple and is omitted (detail can be found in [23]).   Proof We write X ± in place of X N,± for convenience. By Stone-Weierstrass Theorem, C 2 (D ± ) is dense in C(D ± ) in uniform topology. It suffices to check that {( X + , φ + , X − , φ − )} is relatively compact in D([0, T ], R 2 ) for all φ ± ∈ C 2 (D ± ) (cf. Proposition 1.7 (p.54) of [27]) for this weak tightness criterion). By Prohorov's theorem (see Theorem 1. (1) For all t ∈ [0, T ] and ε 0 > 0, there exists a compact set K(t, ε 0 ) ⊂ R 2 such that (2) For all ε 0 > 0, We first check (1). Since φ ± is bounded on D ± and | X + t , 1 | ≤ 1 for all t ∈ [0, ∞), we can always take . To verify (2), since |(x 1 , y 1 ) − (x 2 , y 2 )| R 2 ≤ |x 1 − x 2 | + |y 1 − y 2 |, we only need to focus on X + . By Lemma 4.1, So we only need to verify (2) with φ, X + t − φ, X + s replaced by each of the 3 terms on RHS of the above equation (4.10).
For the first term of (4.10), we apply Lemma 4.3 for the case q = 2 and Y N (r) = A + ε φ, X + r . Since φ ∈ C 2 (D + ), we have for some constant C(φ) which only depends only on φ. Using Lemma 2.13, we have For the second term of (4.10), by ( We have proved that (2) is satisfied. Hence {(X N,+ , X N,− )} is relatively compact. Using (2) and the metric of E, we can check that any subsequential limit L * of the laws of {(X N,+ , X N,− )} concentrates on C([0, ∞), E). For example, The reason is that the two processes can jump at different times (t = 1 and t = 1 + 1 n ) that become identified in the limit (only one jump at t = 1); this can be avoided if one of the two processes is C-tight (i.e. has only continuous limiting values), which is satisfied in our case since X N,+ and X N,− turns out to be both C-tight.

Identifying the limit
Suppose (X ∞,+ , X ∞,− ) is a subsequential limit of (X N,+ , X N,− ), say the convergence is along the subsequence {N }. By the Skorokhod representation Theorem, the continuity of the limit in t and [  where G is the collection of open subsets U ⊂ D such that |U | ≤ r|D| and ∂U ∩ D is a manifold of class C 0,1 (i.e. each point x ∈ ∂U ∩ D has a neighborhood which can be represented by the graph of a Lipschitz function). Moreover, S(r, D) = S(r, aD) for all a > 0.
In this subsection, we establish a discrete analogue for the relative isoperimetric inequality (Theorem 5.5).
First, we study the scaled graph aD ε and gather some basic properties of a continuous time random walk on a finite set.

CTRW on scaled graph aD ε
Recall the m -symmetric CTRW X ε on D ε defined in the Subsection 2.1.2. The Dirichlet form (E ( ) , l 2 (m )) of X ε in l 2 (m ) is given by where µ xy = µ D ε x,y are the conductance on the graph D ε defined in the Subsection 2.1.2. The stationary measure π = π D ε of X ε is given by We now consider the scaled graph aD ε = (aD) aε , which is an approximation to the bounded Lipschitz domain aD by square lattice aεZ d . Clearly the degrees of vertices are given by v aD ε (ax) = v D ε (x). Define the function ρ (aD) on aD by ρ (aD) (ax) := ρ(x). Then define the CTRW X aD ε using ρ (aD) as we have done for X ε using ρ.
The mean holding time of X aD ε is (aε) 2 /d. Clearly, the symmetrizing measure m aD ε and the stationary probability measure π aD ε have the scaling property m aD ε (ax) = a d m D ε (x) and π aD ε (ax) = π D ε (x). Let p aε aD be the transition density of X aD ε with respect to the symmetrizing measure m aD ε . Then a d p aε aD (a 2 t, ax, ay) = p ε D (t, x, y) for every t > 0, ε > 0, a > 0 and x, y ∈ D ε .
We will simply write m and π for the symmetrizing measure and the stationary probability measure when there is no ambiguity for the underlying graph.

An extension lemma
Following the notation of [36], we let G be a finite set, K(x, y) be a Markov kernel on G and π the stationary measure of K. Note that a Markov chain on a finite set induces a natural graph structure as follows. Let Q(e) := 1 2 (K(x, y)π(x) + K(y, x)π(y)) for any e = (x, y) ∈ G × G. Define the set of directed edges E := {e = (x, y) ∈ G × G : Q(e) > 0}.
We use the following 2 different notions for the "boundary" of A ⊂ G: Observe that each edge in ∂ e A is counted twice. Set (K(x, y)π(x) + K(y, x)π(y)). We call 1/S π (r, G) an isoperimetric constant of the chain (K, π). It provides rich information about the geometric properties of G and the behavior of the chain (cf. [36]).
In our case, G = aD ε , π(x) = m(x) m(aD ε ) and K(x, y) = p xy in aD ε , where p x,y is the one-step transition probabilities of X aD ε defined in subsection 2.1.2. For a = 1 and A ⊂ D ε , we have In this notation, we have Figure 6 for an illustration.
It is easy to check that A is grid-connected if and only if for every x, y ∈ A, there exists The following is a key lemma which allows us to derive the relative isoperimetric inequality for the discrete setting from that in the continuous setting, and hence leads us to Theorem 5.5. such that if ε ∈ (0, ε 1 ), then for all grid-connected A ⊂ D ε with π srw (A) ≤ r, we can find a connected open subset U ⊂ D which contains A and satisfies: Proof Since the proof for each r ∈ (0, 1) is the same, we just give a proof for the case r = 1/2. Figure 6 for an illustration.) Note that we cannot simply take U = (W 1 ) o because (d) may fail, for example when ∆A contributes too much to ∂U ∩ D, i.e., when (W 1 ) ∆ := ∂W 1 ∩ x∈∆A ∂U x is large. However, ∆A ⊂∂D ε is close to ∂D and so we can fill in the gaps between ∆A and ∂D to eliminate those contributions. In this process, we may create some extra pieces for ∂U ∩ D, but we will show that those pieces are small enough. Following this observation, we will eventually take U = (W 1 ∪ W 2 ) o where W 2 ⊂ D h for some small enough h > 0.
Points in ∆ 1 A are marked in solid black in Figure 7. For x ∈ ∆ 1 A, consider the following closed cube centered at x: Let Θ x be the union of all connected components of T x ∩ D whose closure intersects U x and define Step 2: (Fill in the gaps between ∂D and (W 1 ) ∆ near ∆ 2 . See Figure 8) Note that ∪ x∈∆1A ∂U x does not contribute to ∂(W 1 ∪ W 3 ) ∩ D. Let W 4 be the union of all connected components of D \ (W 1 ∪ W 3 ) whose closure intersects U x for some x ∈ ∆ 2 A. • (W 1 ) ∇ = ∅, where (W 1 ) ∆2 is defined analogously as (W 1 ) ∆ , with ∆ replaced by ∆ 2 , Now σ(U ∂ ) ≤ σ((W 1 ) ∂ ) ≤ |∂A| 2dε d−1 . Moreover, each x ∈ ∂A is adjacent to at most 3 d − 1 points in ∆ 1 A ∪ ∆A, and for each x ∈ ∆ 1 A, there are at most | B(10Rε)| ≤ (20R + 1) d cubes in T x . So we have Hence (d) is satisfied.
To complete the proof, it suffices to show that W 4 ⊂ D (10R)ε . This is equivalent to show that any curve in D \ W 1 ∪ W 3 starting from any point in (W 1 ) ∆2 must lie in D (10R)ε .
Henceγ(t * ) ∈ ∂U x for some x ∈ ∆ 1 A. We therefore have γ(s j ) ∈ Θ x . The proof is now complete.

Discrete relative isoperimetric inequality
Let π srw be the stationary measure of the simple random walk (SRW) on the graph under consideration and recall Definition 5.2.
Proof We can also assume that A is grid-connected. This is because We thus have the first inequality since 4d|∂A| ≥ |∂ e A| ≥ 2|∂A|. The second inequality follows from the first since q(e) = (ε) d 2d m(D ε ) .

Proof of local CLT
The following weak convergence result for RBM with drift is a natural generalization of [6, Theorem 3.3].
Theorem 5.13. Let D ⊂ R d be a bounded domain whose boundary ∂D has zero Lebesque measure. Suppose D also satisfies: C 1 (D) is dense in W 1,2 (D).
(ii) (X 2 −k , P x k ) converges weakly to (X, P x ) in the Skorokhod space D([0, T ], D) whenever x k converges to x ∈ D.
Proof For (i), the proof follows from a direct modification of the proof of [6,Theorem 3.3]. Recall the definition of the one-step transition probabilities p xy , defined in the paragraph that contains (2.4) and (2.5). Observe that, since ρ ∈ C 1 (D), approximations using Taylor The process X ε has a Lévy system (N ε (x, dy), t), where for x ∈ D ε , N ε (x, dy) = d ε 2 z: z↔x p xz δ {z} (dy).
Following the same calculations as in the proof of [6, Theorem 3.3], while noting that [13, Theorem 6.6.9] (in place of [6, Theorem 1.1]) can be applied to handle general symmetric reflected diffusions as in our present case, we get part (i). Part (ii) follows from part (i) by a localization argument (cf. [7, Remark 3.7]).
We can now present the proof of the local CLT.
Proof of Theorem 2.12. For each ε > 0 and t > 0, we extend p ε (t, ·, ·) to D × D in such a way that p ε is nonnegative and continuous on (0, ∞) × D × D, and that both the maximum and the minimum values are preserved on each cell in the grid εZ d . This can be done in many ways, say by the interpolation described in [2], or a sequence of harmonic extensions along the simplexes (described in [23]).
On other hand, by part (ii) of Theorem 5.13, if the original sequence {ε k } is a subsequence of {2 −k }, then q = p. More precisely, the weak convergence implies that for all t > 0, D φ(y)p(t, x, y)dy = D φ(y)q(t, x, y)dy for all φ ∈ C c (D) and x ∈ D.
Then by the continuity of both p and q in the second coordinate, we have q = p on (0, ∞) × D × D. Since p(t, ·, ·) and q(t, ·, ·) are continuous on D × D (cf. [1]), we obtain p = q on (0, ∞) × D × D. In conclusion, we have p ε converges to p locally uniformly through the sequence {ε n = 2 −n ; n ≥ 1}.