Hydrodynamic and hydrostatic limit for a generalized contact process with mixed boundary conditions

We consider an interacting particle system which models the sterile insect technique. It is the superposition of a generalized contact process with exchanges of particles on a finite cylinder with open boundaries (see Kuoch et al., 2017). We show that when the system is in contact with reservoirs at different slowdown rates, the hydrodynamic limit is a set of coupled non linear reaction-diffusion equations with mixed boundary conditions. We also prove the hydrostatic limit when the macroscopic equations exhibit a unique attractor.


Introduction
In this paper, we consider the interacting particle system introduced in [19] to model the sterile insect technique.This technique was developed, among others, by E. Knipling (see [18]) to eradicate New World screw worms in the 1950s, a serious pest for warm blooded animals.The method is still used today, for instance in France, to protect crops from the very invasive Mediterranean flies, and it is also being tested to fight mosquitoes which transmit dengue in countries like Panama or Brazil.We refer to [7] for a detailed list of trials and programs regarding that method.The sterile insect technique works as follows: male insects are sterilized in captivity using gamma rays.They are then released in the wild population, where females mate only once, giving rise to no offspring if they mate with a sterile male.When enough sterile individuals are released, the wild population eventually becomes extinct.From a mathematical perspective, the sterile insect technique has mainly been modelled in a deterministic way through the study of partial differential equations (see [1]).
This technique was studied from a probabilistic perspective in [19] and [20] using interacting particle systems.In [19], a phase transition result was proved at the microscopic level.Recently, another probabilistic model was studied in [16], also at the microscopic level.In [20], the study is carried out at the macroscopic level (hydrodynamic limit) in both finite and infinite volume with reservoirs, in order to account for the migration/immigration mechanism.
Here, we aim at studying the hydrodynamic and hydrostatic limits of that interacting particle system under the effect of slow reservoirs in any dimension d ≥ 1.The slow-down mechanism models the fact that beyond the boundary through which insects arrive into the system or leave it, there are very few insects (the exterior of the system might be a territory which is much less favorable to the development of these insects).
In the perspective of interacting particle systems, the sterile insect technique is modelled as follows: individuals evolve on a d-dimensional finite set N , where N ≥ 1 and T d−1 N = (Z/N Z) d−1 .The evolution of the population is described by a continuous time Markov process (η N t ) t≥0 with state space E B N where E is a countable set.In this model, the gender does not come into account so we refer to sterile individuals rather than sterile males.The quantity of interest here is not the number of insects per site but the types of insects present at a given site.Precisely, E = {0, 1, 2, 3} and for x in B N , 0 if there are no insects in x, 1 if there are only wild insects in x, 2 if there are only sterile insects in x, 3 if there is a combination of wild and sterile insects in x.
The dynamics of the Markov process is the superposition of three Markovian jump processes: (i) An exchange dynamics which models the fact that insects move in an isotropic way within the bulk B N and which is parameterised by a diffusivity constant D > 0. Precisely, for a configuration η and x, y two sites in B N , the states of sites x and y in η are exchanged at rate D.
(ii) A birth and death dynamics within B N which models births of individuals due to the mating of a wild individual with wild or sterile insects, as well as deaths of individuals.This is the dynamics introduced in [19] that was referred to as a contact process with random slowdowns (CPRS).It is parameterised by a release rate r > 0 and growth rates λ 1 , λ 2 > 0. Sterile individuals are injected on a site at rate r independently of everything else.The rate at which wild individuals give birth (to wild individuals) on neighbouring sites is λ 1 at sites in state 1, and λ 2 at sites in state 3. Sterile individuals do not give birth.We take λ 2 < λ 1 to reflect the fact that fertility is reduced at sites in state 3. Deaths for each type of insects occur independently and at rate 1.
Note that without the presence of sterile insects, the CPRS would be a basic contact process (as defined for instance in [23]) with parameter λ 1 , and the presence of sterile insects can be interpreted as a random decrease of the fertility rate due to the presence of sites containing sterile and wild individuals.In [19], the microscopic study of the contact process with random slowdowns in dimension d ≥ 1 leads to the following phase transition result: for certain values of λ 1 and λ 2 , when r is large enough, the healthy population almost surely becomes extinct and survives otherwise.In this paper, the CPRS will be called generalized contact process.In [20], the hydrodynamic limit of the superposition of the three dynamics above, where the first and the third one are accelerated in the diffusive scaling N 2 , and where θ = θ r = 0, is proven to be a system of non linear reaction-diffusion equations with Dirichlet boundary conditions in any dimension.
In this paper, we first prove the finite volume hydrodynamic limit of this particle system for any values of θ , θ r ≥ 0 and in any dimension.The hydrodynamic equation obtained has mixed boundary conditions which depend on the values of θ , resp.θ r .Precisely, for θ ∈ [0, 1), resp.θ r ∈ [0, 1), we get a Dirichlet type boundary condition at the left-hand side, resp.right-hand side of the system.For θ = 1, resp.θ r = 1, we get a Robin type boundary condition at the left-hand side, resp.right-hand side of the system.For θ > 1, resp.θ r > 1, we get a Neumann type boundary condition at the left-hand side, resp.right-hand side of the system.
We then prove the finite volume hydrostatic limit of the particle system for a specific class of parameters regarding the dynamics.Within that class of parameters, the sequence of invariant measures of the interacting particle system is associated to a profile which is the stationary solution of the hydrodynamic equation with corresponding mixed boundary conditions.
Our paper is, up to our knowledge, the first one regarding the effect of mixed reservoirs in and out of equilibrium (hydrodynamic and hydrostatic limit) for a multi species process in a bounded d-dimensional cylinder.Note that all our results can be extended to the d-dimensional hypercube, [−1, 1] d , following the method in [21].We believe that the analysis for general domains would require more effort, in particular, regarding the choice of a suitable discretization of the underlying macroscopic space.The discretization issue has been addressed for some conservative interacting particle systems evolving on a bounded Lipschitz domain (see [5] and references therein).For domains such as manifolds, we refer to [25] for the symmetric simple exclusion process with no reservoirs.Both papers [5] and [25] rely on duality techniques.The effect of reservoirs on a one dimensional conservative system has been widely studied in finite volume (see for instance [6], [10]).Much is now known both at the microscopic and macroscopic level.Recently, the effect of slow reservoirs has aroused considerable interest for the symmetric simple exclusion process in one dimension (see for instance [2], [12], [13], [14] and references therein).In [11], the authors proved a hydrostatic principle for a boundary driven gradient symmetric exclusion process using the fact that the stationary profile is a global attractor for the hydrodynamic equation.This method inspired our proof for the hydrostatic limit.However, the fact that we obtained coupled equations for the hydrodynamic limit, and that we work in any dimension make the analysis more subtle and require general analytical tools.
This paper is organized as follows.In Section 2 we introduce the notations and state our results.The proof of the hydrodynamic limit for each of these regime is established in Section 3 via the Entropy Method.Among other things, as we work in arbitrary dimension, some care must be taken to perform the replacement lemma, and also to define and characterize the solution of the hydrodynamic limit at the boundary, through the use of the Trace Operator (see subsection 3.4).The difficulties due to different boundary slowing exponents are purely analytical, and appear when proving uniqueness of the hydrodynamic equation with mixed boundary conditions.The proof of the hydrostatic limit, established in Section 4 relies on the use of a well chosen change of coordinates for the coupled equations.Under this change of coordinates (inspired by some simulations see Appendix B), a comparison principle holds.It allows us to find a unique attractor when some conditions on the parameters are satisfied.Outside that class of parameters, although uniqueness of the invariant measure holds, we do not even know whether there is uniqueness of the stationary solution of the hydrodynamic equation and simulations show that for Neumann type boundary conditions there are several stationary profiles.However, we believe that a more general hydrostatic principle in the spirit of the one proved in [22] is valid.

The microscopic model
The dynamics of our interacting particle system is given by three generators, one for the exchange dynamics, one for the generalised contact dynamics and one for the boundary dynamics.In order to explicit each one of those generators, let us give a few notations.Let N ∈ N and d ≥ 1.For p ≥ 1, we write T p N , resp.T p the discrete, resp.continuous torus (Z/N Z) p , (R/Z) p .Denote by the boundary, resp.right-hand side boundary, resp.left-hand side boundary of the bulk.Denote B = (−1, 1) The microscopic state space is denoted Ω N := {0, 1, 2, 3} B N and its elements, called configurations, are denoted η.Therefore, for x ∈ B N , η(x) ∈ {0, 1, 2, 3}.To describe the dynamics of our model, we will use the correspondence introduced in [20] between the state space Ω N and Σ N := ({0, 1} × {0, 1}) B N where the correspondence between an element (ξ, ω) ∈ Σ N and η ∈ Ω N is given as follows: for x ∈ B N , In other words, ω ∈ {0, 1} represents the presence of sterile insects, and ξ ∈ {0, 1} that of wild ones on a given site, i.e., (ξ(x), ω(x)) = (0, 0) if x is in state 0, (1, 0) if it is in state 1, (0, 1) if it is in state 2 and (1, 1) if it is in state 3. Also, in order to describe the evolution of the density of sites in state 1, 2, 3, resp 0, we define for x in B N and a configuration η ∈ Ω N with associated configuration (ξ, ω) Finally, we also express the correspondence (2.1) by the following application from Σ N to Ω N : η = η(ξ, ω), where, for any x ∈ B N , η(x) = 2ω(x) + ξ(x). (2.3) • Generator for the exchange mechanism: it corresponds to the usual stirring mechanism where each site has an exponential clock with rate D and independent from all the other clocks, where D is a fixed positive parameter.When the clock rings, a neighbouring site is chosen uniformly at random and the states of both sites are exchanged.The action of the generator on functions f : Σ N → R is therefore given by: where (e • Generator for the generalised contact process in the bulk: following the description of the generalised contact process in the introduction, let us give the rates of the generalised contact in the bulk.For , where x ∼ y means that x and y are neighbouring sites in B N .
Births and arrivals of sterile individuals at x happen at the following rates: and 0 → 2 at rate r, 1 → 3 at rate r, Deaths of individuals at x happen at the following rates: 1 → 0 at rate 1, 2 → 0 at rate 1, 3 → 1 at rate 1, 3 → 2 at rate 1, (2.6) Therefore, using the correspondences (2.1) and (2.2), the generator L N = L N,λ1,λ2,r of the generalised contact process acts as follows on functions f : Σ N → R: where for x ∈ B N , where, for ζ ∈ {0, 1} B N , σ x ζ is the configuration obtained from ζ by flipping the configuration at x, i.e.
• Generator for the boundary dynamics: the generator of the dynamics at the boundary is parametrized by θ = (θ , θ r ) with θ , θ r ≥ 0 and a positive function b and, the restriction of g to Γ is equal to b.The dynamics at the boundary can then be described as follows: a site x ∈ Γ − N , resp.x ∈ Γ + N flips from state i ∈ {0, 1, 2, 3} to state j ∈ {0, 1, 2, 3} \ {i} at rate N −θ b j (x/N ), resp.N −θr b j (x/N ) .In order to express the generator of the boundary dynamics, we make use of η i = η i (ξ, ω) for i ∈ {0, 1, 2, 3} which is the configuration in {0, 1} B N obtained from (ξ, ω) ∈ Σ N according to (2.2).For f : Σ N → R, the boundary generator acts on f as follows: where b 0 (x/N ) := 1 − 3 i=1 b i (x/N ) and with σ i,x (ξ, ω) := σ i,x η(ξ, ω), the configuration in Σ N associated to σ i,x η, where Fix a time horizon T > 0 and denote {(ξ N t , ω N t ), t ∈ [0, T ]} the Markov process associated to the generator Let D Σ N ([0, T ]) be the path space of càdlàg trajectories with values in Σ N .Given a measure µ N on Σ N , denote by P µ N the probability measure on D Σ N ([0, T ]) induced by µ N and (ξ t , ω t ) t∈[0,T ] and denote E µ N the expectation with respect to P µ N .
Invariant measures for the exchange and boundary dynamics: (2.12) Denote ν N α the Bernoulli compound product measure on B N with parameter α: for (ξ, ω) ∈ Σ N , where α 0 = 1 − α 1 − α 2 − α 3 and where Z α,N is the normalizing constant Note that ν N α is such that for every 1 ≤ i ≤ 3 and x ∈ B N , ).One can verify the following statements: • Consider α a smooth profile satisfying (2.12) and ∀x ∈ Γ, α(x) = b(x). (2.14) Then, ν N α is an invariant and reversible measure for the boundary dynamics so for any f : • Consider α a constant profile.Then ν N α is an invariant and reversible measure for the exchange dynamics so for any f : Σ N → R, For any θ ∈ (R + ) 2 , at fixed N , the dynamics defined by (2.11) is irreducible and the state space is finite.Therefore, the dynamics has a unique invariant probability measure that in the sequel we denote µ ss N ( θ).
Useful (in)equalities: For any A, B > 0, (2.17) For any a, b, A and N ∈ N, For any sequences of positive numbers (a

The macroscopic equations
Let us first introduce a few notations.We will write functions with values in R with normal letters (for instance G) and the ones with values in R 3 with letters with a hat (for instance G) .For n, m ∈ N, denote C n,m ([0, T ] × B) the space of functions that are n times differentiable in time and m times differentiable in space, C n,m 0 , resp.C n,m 0,− , resp.C n,m 0,+ the ones in C n,m ([0, T ] × B) which are zero on Γ, resp.Γ − , resp.Γ + .Denote C ∞ k (B) the space of smooth functions with compact support in B, C m (B) the space of functions that are m times differentiable in space with C m 0 (B), resp.C m 0,− (B), resp.C m 0,+ (B) those which are zero on Γ, resp.
Γ − , resp.Γ + , and C(B) the space of continuous functions on B .For θ = (θ , θ r ) in (R + ) 2 , we will use the following notation to denote these functional spaces: (2.20) Let , be the L 2 inner product and , µ the inner product with respect to a measure µ.For f = (f 1 , f 2 , f 3 ) and g = (g 1 , g 2 , g 3 ) in L 2 (B) 3 , f , g = 3 i=1 f i , g i .Denote (e 1 , • • • , e d ) the canonical basis of Z d .Introduce the Sobolev space H 1 (B) which we recall to be the set of functions g ∈ L 2 (B) such that for any 1 ≤ k ≤ d, there exists an element denoted where ∂ e k ϕ is the usual partial derivative.The H 1 (B) norm is then defined as follows: We will write g 2 2 instead of g 2 L 2 (B) when no confusion is possible.Introduce In order to define the value of an element G in H 1 (B) at the boundary, we need to introduce the notion of trace of functions on such Sobolev spaces.The trace operator in the Sobolev space H 1 (B) can be defined as a bounded linear operator, Tr : H 1 (B) → L 2 (Γ) such that Tr extends the classical trace, that is Tr(G) = G |Γ , for any G ∈ H 1 (B) ∩ C(B).We refer to [8, Part II Section 5] for a detailed survey on the trace operator.
To lighten notations, for a function G depending on time and space we will often write G s instead of G(s, .).Finally, for θ ∈ (R + ) 2 , introduce the following linear functional on L 2 [0, T ], H 1 (B) parameterised by a test function G in C θ : for t ∈ [0, T ], where with The hydrodynamic equation is a system of coupled reaction diffusion equations with mixed boundary conditions depending on θ.If θ , resp.θ r belongs to [0, 1), the boundary conditions are of Dirichlet type on Γ − , resp.Γ + .If θ = 1, resp θ r = 1, they are of Robin type on Γ − , resp.Γ + .If θ > 1, resp.θ r > 1, they are of Neumann type on Γ − , resp.Γ + .We will focus on the cases where θ ∈ [0, 1), θ r = 1 resp.θ > 1, θ r = 1 corresponding to a Dirichlet boundary condition on Γ − and a Robin boundary condition on Γ + , resp. a Neumann boundary condition on Γ − and a Robin boundary condition on Γ + .All the other cases can be adapted from those ones (see Table 1).Definition 1.Let γ : B → R 3 be a continuous function.
is a weak solution of the Dirichlet + Robin mixed boundary problem for any function G ∈ C θ , for any t ∈ [0, T ], where n 1 (r) is the outward normal unit vector to the boundary surface Γ and dS(r) is an element of surface on Γ.And, ρ(0, .)= γ(.)almost surely. (2.26) if ρ satisfies conditions (2.24) and (2.26) as well as the following condition (2.28): for any G ∈ C θ , for any t ∈ [0, T ], (2.28) Remark 1.In (2.25), the integral over Γ − corresponds to the Dirichlet boundary condition.In (2.28) the integral over Γ − comes from an integration by parts of the terms involved in the bulk.Both in (2.25) and (2.28) the first integral over Γ + comes from an integration by parts of the terms involved in the bulk and the second integral over Γ + corresponds to the Robin boundary condition.

Hydrodynamic and hydrostatic results
Let us state the main results proved in this paper.The first one (Theorem 1) establishes the hydrodynamic limit of the dynamics defined above and the second one (Theorem 2) establishes its hydrostatic limit.Before stating Theorem 1, let us first define the empirical measure π N (ξ, ω) = π N associated to a given configuration (ξ, ω).Recall how in (2.2), we built where δ x/N is the point mass at The empirical measure is therefore the triplet of empirical measures associated to the density of sites in state 1, resp.2, resp.3. Denote M the set of positive measures on B with total mass bounded by 3. The process ( π N t ) t≥0 = ( π N (ξ t , ω t )) t≥0 , is a Markov process with state space M 3 and its trajectories belong to D([0, T ], M 3 ), the path space of càdlàg time trajectories with values in M 3 .We endow the path space with the Skorohod topology (we refer to [3] for a detailed presentation on the Skorohod topology).For θ ∈ (R + ) 2 and µ N a probability measure on Σ N , denote Q θ N = P µ N ( π N ) −1 the law of the process ( π N (ξ t , ω t )) t≥0 when (ξ 0 , ω 0 ) ∼ µ N and where (ξ t , ω t ) t≥0 evolves according to the dynamics given by (2.11), with parameter θ for the boundary reservoirs.The hydrodynamic result states as follows: Theorem 1. (Hydrodynamic limit).For any sequence of initial probability measure (µ N ) N ≥1 on Σ N , the sequence of probability measures (Q θ N ) N ≥1 is weakly relatively compact and all its converging subsequences converge to some limit Q θ, * concentrated on the set of weak solutions of the hydrodynamic equation that are in L 2 ([0, T ]; H 1 (B)), in the sense of Definition 1. Furthermore, if there is an initial continuous profile γ : B → [0, 1] 3 such that for any δ > 0 and any converges to the Dirac mass Q θ concentrated on the unique weak solution ρ of the boundary value problem associated to θ and with initial condition γ.Therefore, for any t ∈ [0, T ], δ > 0 and any function We prove Theorem 1 in Section 3.
Intuitively, the sterile insect technique is more effective when r is large and λ 2 is small.This result has been made precise at the microscopic level for the generalised contact process, in [19,Theorem 2.5], where the author proved a phase transition result: when λ 1 > λ 2 are properly tuned, there is a critical value r c > 0 below which the wild population survives with strictly positive probability and above which the wild population dies out almost surely.
To establish the hydrostatic limit, depending on the parameters at the boundary, we will need one of the following sets of conditions to be satisfied: where δ 1 is the smallest eigenvalue of the Laplacian with Dirichlet boundary conditions (see (3.68)).These conditions will appear technically in the proof of Theorem 2. However, note that at fixed D and λ 1 , this corresponds to having r large and λ 2 small enough, which is consistent with the effectiveness of the sterile insect technique.
Recall that µ ss N ( θ) denotes the sequence of unique invariant measures for the irreducible dynamics defined by (2.11).The hydrostatic result states as follows.Theorem 2. (Hydrostatic limit).Suppose that conditions (H 1 ) hold.There exists a unique stationary solution of (2.23) that we denote ρ D,R , and a unique stationary solution of (2.27) that we denote ρ Ne,R .Furthermore, the following statements hold.
Remark 2. For all the other mixed boundary regimes corresponding to other values of θ, the hydrostatic principle states in the same way, replacing ρ D,R i or ρ Ne,R i by the stationary solution of the associated hydrodynamic equation.In the cases where only Dirichlet and Robin boundary conditions are involved, we can slightly weaken the conditions (H 1 ) by using conditions (H 2 ) or (H 3 ) instead.Precisely: in the (D;D), (D;R), (R;D) regimes, the hydrostatic principle holds under conditions (H 2 ) and in the (N e ; N e ) regime, it holds under conditions (H 3 ).
The proof of Theorem 2 is done in Section 4. It essentially relies on an intermediate result stated in Theorem 6 regarding the convergence of solutions of the hydrodynamic equation towards the unique stationary state.This result is non standard as it involves a system of coupled equations and we prove it in the second Subsection of Section 4.

Proof of the hydrodynamic limit
As said before, we focus on the cases where θ ∈ [0, 1), θ r = 1 and θ > 1, θ r = 1.We follow the entropy method introduced by Guo, Papanicolaou and Varadhan in [15] to prove the hydrodynamic limit.First, we prove tightness of the sequence of measures (Q θ N ) N ≥1 .Then, we show that any limit point of ( Finally, we prove uniqueness of the solution of the hydrodynamic equations at fixed initial data.We do not give details for the standard steps but rather, insist on the specific difficulties arising in our case, namely, the d-dimensional replacement lemmas (subsections 3.2.3 and 3.2.4) and the uniqueness of the solution of the hydrodynamic equation (sections 3.5).

The martingale property and tightness
Recall from (2.11) the definition of the total generator L N .By Dynkin's formula (see [17, is a martingale with respect to the natural filtration F t = σ(η s , s ≤ t) and with quadratic variation given by: We then have that is also a martingale whose quadratic variation is known.In order to develop the integral terms in (3.1), introduce the discrete second derivative in the direction e k (for 1 ≤ k ≤ d) in the bulk, the discrete Laplacian, and the discrete gradient in the direction e 1 at the boundary: Computations yield where we used that with The second and third lines in (3.3) correspond to the computation of the time integral associated to N 2 L N , the fourth line in (3.3) corresponds to the time integral associated to L N and the last term, to the integral associated to N 2 L b, θ,N .
For i ∈ {1, 2, 3}, a computation of the quadratic variation of the martingale M N i,t ( G) shows that its expectation vanishes as N ↑ ∞.Therefore, by Doob's inequality, for every δ > 0, lim sup Proposition 1.The sequence of probability measures We refer to [17,Chapter 4] for details regarding the proof of tightness of a sequence of probability measures.It is enough to show that for every H in a dense subset of C(B) for the L 2 norm, for every 1 , so that H vanishes at the boundary.To prove that, we use the martingale and its quadratic variation introduced in (3.1) and (3.2), and show that lim sup δ→0 lim sup and lim sup δ→0 lim sup We get (3.7) using the triangular inequality and (3.5).To prove (3.8), we show that there is a constant C depending only on H such that for every r ∈ [0, T ], For that, we use the decomposition of L N and the fact that H vanishes at the boundary as well as explicit computations and the fact that the f i 's are uniformly bounded in N .

Replacement Lemmas
In order to characterize the limit points of a sequence (Q θ N ) N ≥1 , we need to close the equation (3.3).That means that we want to show that each term of the martingale converges to a term that appears in the weak formulation of the solution of the hydrodynamic equation, and that the martingale converges to zero.For that, we perform a replacement lemma in the bulk and one at the boundary.The replacement lemma in the bulk (Proposition 2) is exactly the same as in [20,Lemma 4.2] and we refer to that article for a detailed proof.Here we focus on the replacement lemmas at the boundary and more specifically on the left-hand side boundary (the same statements hold on the right-hand side).There are two replacement lemmas: one for θ ∈ [0, 1) whose formulation coincides with the replacement lemma at the boundary in [20,Proposition 4.3] (corresponding to a Dirichlet condition), and one for θ r ≥ 1, whose formulation involves particle densities over small macroscopic boxes.

Dirichlet forms
Let us recall the expressions introduced in [20, Section 5] of the Dirichlet forms associated to each dynamics.For that, recall the correspondences (2.1) and (2.2).For f : Σ N → R and µ a measure on Σ N , and In the proofs of the Replacement lemmas, we will widely make use of the following inequalities.
Lemma 1. (i) Consider α a smooth profile which satisfies (2.12).There is a constant C 1 > 0 such that for any density function f : Σ N → R with respect to the measure ν N α , (3.10) (ii) Consider α a smooth which satisfies (2.12) including constants.There is a constant C 2 > 0 such that for any density function f : Σ N → R with respect to the measure ν N α , (iii) Consider α a smooth profile which satisfies (2.12) and (2.14), then for any density function f : Σ N → R with respect to the measure ν N α , We refer to [20], Lemma 6.1 for the proof.The authors use some change of variable formulas in the same spirit as those given in (A.1) and (A.2).For point (iii), they use an alternative expression of the boundary generator expressed in terms of (ξ, ω).

Replacement lemma in the bulk.
Let us first introduce a few notations.Given a smooth profile α, and a function φ : where y − x = max{|y i − x i |, 1 ≤ i ≤ d}, and denote η i (x) the average of η in Λ x , that is, (3.14) Introduce the vector and for ε > 0, In the sequel, we will write εN instead of εN .The replacement lemma in the bulk stated and proved in [20, Lemma 4.2] is the following: and for any function φ : 3 Replacement lemma at the left-hand side boundary for θ ∈ [0, 1).
Here we fix θ in [0, 1) and prove the replacement lemma at the left-hand side boundary.It essentially states that when performing the macroscopic limit N → ∞, we can replace η i (x) by b i (x/N ).For θ r ∈ [0, 1), the replacement lemma at the right-hand side boundary is exactly the same.Recall that this result has been proved for θ = θ r = 0 in [20, Section 6] and we generalize it here to the case where the left-hand side (or right-hand side) parameter θ is allowed to vary in [0, 1).
Proposition 3.For any sequence of measures (µ N ) N ≥0 on Σ N , for any G ∈ C 1,2 ([0, T ] × B) and any i ∈ {1, 2, 3}, for any t ∈ [0, T ], for all δ > 0, lim sup Note that the replacement lemma at the right-hand side boundary for θ r ∈ [0, 1) states as above, with the sum in x carrying over Γ + N rather than Γ − N .Proof.Fix an i ∈ {1, 2, 3}.It is enough to show that lim sup Consider α a smooth profile satisfying conditions (2.12) and (2.14).For a > 0, We used, in the first inequality, that the Radon-Nikodym derivative of µ N with respect to ν N α is bounded by exp(K 0 N d ) with K 0 a constant, and Chebychev's inequality in the second line.Therefore, It is enough to show that the last term is uniformly bounded in a and N and then, take a → ∞ with N .Since e |x| ≤ e x + e −x , using inequality (2.19), we show that the last term in (3.16) without the absolute values, is uniformly bounded in a and N .Applying Feynman-Kac's inequality (see [17,Appendix A.1]) with we get that where the supremum is taken over densities with respect to ν N α .Note that for and, for j = i, by the change of variable presented in (A.2), Therefore, where we used (2.18) in the last line replacing A by AN , with A > 0. Summing (3.18) over Γ − N and multiplying by a N d−1 yields, where the second term in the first inequality comes from Cauchy-Schwarz's inequality, the fact that f is a density, the change of variable formula (3.17) and the fact that each coordinate of b is bounded by 1.Therefore, using (3.10), (3.11) and (3.12) to bound L N √ f , √ f ν N α and the fact that a Dirichlet form is positive we are left with (3.20) and taking a → ∞, the result follows.

3.2.4
Replacement lemma at the left-hand side boundary for θ ≥ 1.
For θ ≥ 1, the replacement lemma at the boundary involves particle densities over small macroscopic boxes.Again, the same replacement lemma holds at the right-hand side boundary for θ r ≥ 1.In fact, we will see in the proof that the lemma holds for any positive value of θ , resp.θ r regardless of whether θ resp.θ r ≥ 1.
Here, as we are working in arbitrary dimension, some care must be taken in the proof when adapting the argument used for instance in [17,Chapter 5].Proof.For a vector First, consider the expression in the expectation without absolute value and the time integral, and rewrite it for any s ∈ [0, t] as where we recall that Λ εN x is defined in (3.13).For since G is twice differentiable in space, a Taylor expansion allows us to bound the first term on the right hand side of (3.23) by dεC G (N ) where C G (N ) is uniformly bounded in N by a constant C G , depending only on G. Now, rewrite the last term in (3.23) as follows: with C G (N ), uniformly bounded by C G , a constant that only depends on G. Therefore, we are left to prove that lim sup ε→0 lim sup Consider α a smooth profile satisfying conditions (2.12) and (2.14) .By the entropy inequality (see [17, Appendix 1]), for any A > 0, As B N is finite, there is a constant K 0 > 0 such that H(µ N |ν N α ) ≤ K 0 N d so the first term in (3.25) is bounded by K 0 /A.Let us show that the second term tends to zero when N → ∞ and ε → 0 and then take A arbitrarily big.Again, by (2.19), it is enough to show that the second term in (3.25) without the absolute values in the exponential, tends to zero.By Feynman-Kac's inequality, where the supremum is taken over densities with respect to ν N α .Now, rewrite R N,ε (G s , η) using a telescopic sum to write the differences η i (x + j 1 e 1 ) − η i (x): Fix 0 ≤ ≤ j 1 ≤ εN .Performing the change of variable (ξ, ω) → (ξ e1+x,( +1)e1+x , ω e1+x,( +1)e1+x ) := (ξ, ω) ,x , and using (A.1), dν N α (ξ, ω).
(3.27) To deal with the first term in (3.27) using inequality (2.18), this term is bounded by where B > 0 will be chosen later and where we used that f is a density with respect to ν N α in the last line.Note that (ξ, ω) → f ((ξ, ω) ,x ) is not a density but we can deal with the last integral term as follows: where in the second line we used the change of variable formula (A.1) and, in the last line, the fact that f is a density with respect to ν N α and that R e1+x,( +1)e1+x i,j ( α) = O(N −1 ) so bounded by C/N where C is a constant.Therefore, for N large enough, for all 0 ≤ ≤ j 1 ≤ εN , We are left with (3.29) This, combined with (3.26) as well as Lemma 1 yields: taking N → ∞, then ε → 0, and then A → ∞, the result follows.

Energy estimates
In view of the proof of uniqueness of the limit of the sequence of probability measures (Q θ N ) N ≥1 , we state that any limiting measure Q θ is concentrated on a trajectory belonging to a specific functional space.This allows to define the hydrodynamic limit at the boundary.Proposition 5. Let θ ∈ (R + ) 2 and Q θ be a limit point of sequence of probability measures (Q θ N ) N ≥1 .Then, the probability measure Q θ is concentrated on paths ρ(t, u)du such that for every This follows from the Lemma below and the Riesz Representation Theorem.Lemma 2. For any θ ∈ (R + ) 2 , there is a constant K θ > 0 such that for every 1 ≤ i ≤ 3, where the supremum is carried over functions For the proof of Lemma 2, which we do not detail here, one can follow the arguments in [17, Lemma 7.2 Chapter 5].First prove (3.31) for a dense and countable set of elements of C 0,2 c ([0, T ] × B) thanks to Feynmann-Kac's inequality.Then, use an integration by parts to deal with the spatial derivatives in H, as well as the change of variable (A.1).To recover Proposition 5 from Lemma 2 and the Riesz Representation Theorem, we also refer to [17, Chapter 5, Theorem 7.1 ].

Characterization of the limit point in the (Dirichlet ; Robin) mixed regime
In order to show that the limit point of the sequence of probability measures (Q θ N ) N ≥1 lies on the trajectory with density profile the unique solution of the hydrodynamic equation associated to θ and γ, we give a characterization result (see Proposition 6).We will focus on the (Dirichlet ; Robin) mixed regime since the (Neumann ; Robin) mixed regime can be proved following the same strategy.Therefore, take θ ∈ [0, 1) and θ r = 1.
As mentioned in the introduction, in one dimension, the macroscopic trajectories are continuous in space and their values at the boundaries are defined in the classical sense.This is no longer valid in higher dimension.To deal with this difficulty we use the regularity of the trajectories proved in Proposition 5: the trajectories lie in L 2 ([0, T ], H 1 (B)) so their values at the boundary are defined via the trace operator (see Lemma 3).Proposition 6.If Q θ is a limit point of the sequence of probability measures (Q θ N ) N ≥1 , then where I G ( ρ) was defined in (2.21).
Proof.The fact that any limit point is concentrated on trajectories which are absolutely continuous with respect to the Lebesgue measure comes from Proposition 5. Let Q θ be a a limit point of the sequence of probability measures (Q θ N ) N ≥1 , To prove (3.32), it is enough to show that for any fixed δ > 0 and G ∈ C 1,2 0,− , Here, note that for s ∈ [0, T ] and r ∈ Γ, ρ i (s, r) stands for Tr(ρ)(s, r) which is well defined since ρ is in . By the triangular inequality, it suffices to prove that for any 1 ≤ i ≤ 3, As usual, we would like to approximate ρ by a convolution of its associated empirical measure with an approximation of the identity.Indeed, that convolution product can then be written in terms of the mean value of the configuration in a microscopic box.This is straightforward in the bulk, however, for the boundary terms, we need to justify that such an approximation works (see (3.44)).Without loss of generality, let us deal with i = 1.We turn to our martingales (3.1) M N 1,t ( G) and recall that we have proved that its quadratic variation vanishes as N ↑ ∞.For ε > 0, introduce the set We now use Proposition 2 to replace the local functions of η by functions of the particle density: From Proposition 3 and Proposition 4, the martingale M N 1,t ( G) can be rewritten as On the other hand, by (3.5) recall that lim sup Now, introduce the following approximations of the identity on B: (3.37) Note that for ε > 0, 1 where for any trajectory π and for any t ∈ [0, T ], where functions F i , i = 1, 2, 3 are defined in (2.22).By approximating Lebesgue integrals by Riemann sums, on the bulk and at the boundary, we obtain lim sup ε→0 lim sup where for any trajectory π and for any t ∈ [0, T ], with π , for each ε > 0, we get for any limit point Q θ of the sequence of probability measures To conclude the proof, it remains to prove that we may replace the convolutions appearing in the functional by the associated density of the trajectory.By Proposition 5, Q θ is concentrated on paths ( π(t, dr)) t∈[0,T ] = ( ρ(t, r)dr) t∈[0,T ] which are absolutely continuous with respect to the Lebesgue measure and such that for every 1 ≤ i ≤ 3, ρ i belongs to L 2 ([0, T ], H 1 (B)).For the replacement of the convolution with the density in the bulk, since u ε is an approximation of the identity in L 1 (B) and the functions F i are Lipschitz, the random variables For the other terms in F G,t 1, , by the dominated convergence Theorem, for almost every trajectory ( π(t, dr)

Uniqueness of the limit points
In order to finish the proof of the hydrodynamic limit specific to each regime we are left to show that each boundary valued problem (2.23) and (2.27) with fixed initial data admits a unique solution.For that, we use the standard method which consists in decomposing the difference of two solutions on the orthonormal basis of a well chosen eigenvectors of the Laplacian.The choice of the family of eigenvectors is not necessarily intuitive and depends on the boundary conditions of the mixed regime considered.We thus give details for both the (Neumann; Robin) and (Dirichlet; Robin) mixed regimes, for which the family of eigenvectors are different.As we are working in dimension d ≥ 1, we will need to control integral terms on the boundary.Therefore, we will make use of the following result regarding the continuity of the trace operator.We refer to [8, Part II Chapter 5] for a detailed survey of the trace operator.
There is a constant C tr > 0 depending only on Ω and p such that for any ϕ ∈ C ∞ (Ω), where .L p (∂Ω) denotes the L p norm on ∂Ω and .W 1,p the Sobolev norm on Ω given by where In particular, C tr = 1.
In the sequel we only make use of (3.46) but we stated Theorem 3 for the sake of completeness.Proof.By Liouville's Theorem stated for instance in [8], there is a countable system {V n , α n , n ≥ 1} of eingensolutions for the problem in H 1 (B) and containing all possible eigenvalues.The set {V n , n ≥ 1} forms a complete, orthonormal system in the Hilbert space L 2 (B) and the eigenvalues have finite multiplicity.Note that for any U, W ∈ H 1 (B), One can check that since we are working on (−1, 1) We have Note that by abuse of notations we indexed the family V k by N \ {0} instead of N × (N \ {0}) d−1 but this is not a problem because we can give an order to elements of N × (N \ {0}) d−1 .
Consider ρ 1 and ρ 2 two solutions of (2.27) associated to the same initial profile and for n ∈ N and t > 0, introduce For that, apply the weak formulation (2.27) with V k : for any 1 .dS(r)ds.
(3.54) Therefore ρ 1 i (t, .)− ρ 2 i (t, .),V k is time differentiable with derivative: and so is G n , with for any A > 0, where we used both the Cauchy-Schwarz and (2.18) inequalities in the last line.By (3.49), (3.50) and (3.52), the right-hand side of (3.56) converges to By the trace inequality (3.46), (3.58) Furthermore, using that ρ 1 and ρ 2 take their values in [0, 1] 3 , there is a constant C := C(λ 1 , λ 2 , r, d) > 0 such that for any ρ a , ρ b ∈ [0, 1] 3 and 1 ≤ i ≤ 3, Then, by Cauchy-Schwarz's inequality, there is a constant C > 0 such that for any 1 ≤ i ≤ 3, Grönwall's inequality and the fact that G(0) = 0 yields G(t) = 0 at any time.Proof.The proof follows the same lines as the previous one except that we consider another family of eigenfunctions of the Laplacian.Indeed, consider the following boundary-eigenvalue problem for the Laplacian: (3.61) Again, one can check that the countable system of eigensolutions {W n , γ n , n ≥ 1} given below (in (3.63)) for the problem (3.61) contains all possible eigenvalues and is a complete, orthonormal system in the Hilbert space L 2 (B), that the eigenvalues γ n have finite multiplicity and that with Again, by abuse of notations we have indexed the W k 's by N * instead of (N * ) d .
As before, take ρ 1 and ρ 2 two solutions of (2.23) with same initial data and introduce and (3.66) Using the weak formulation (2.25) with W k , we get that for any 1 where the Wk = Vk are defined in (3.51).Then, we conclude following exactly the same lines as the proof of Theorem 4.

Uniqueness of the solution in the other regimes
In order to prove uniqueness in the other regimes, one can follow the same classic method used above.The orthonormal basis used to decompose the difference of two solutions as in (3.53) or (3.65) then depends on the boundary conditions.For the (Dirichlet ; Dirichlet) regime, the decomposition is carried out on the eigenvectors of the following boundary-eigenvalue problem for the Laplacian: (3.68) for which the associated family of eigenvectors is with eigenvalues given by where the Ǔk = Vk are defined in (3.51).

Hydrostatic limit
In this section, we prove Theorem 2 which states that when the parameters r, λ 1 , λ 2 , d, D satisfy certain conditions, starting from an invariant measure, the system converges to the stationary profile of the corresponding hydrodynamic equation.Precisely, recall that in Section 2, for θ ∈ (R + ) 2 we defined µ ss N ( θ) as the sequence of unique invariant measures for the irreducible dynamics defined by (2.11).The hydrostatic principle states that this sequence is associated to the unique stationary solution of the hydrodynamic equation, if existence and uniqueness of such a solution hold.For the proof, we were inspired by [11] and the key argument relies on the convergence of all the trajectories satisfying the hydrodynamic equation to the unique stationary profile of these equations.In [11], the convergence of trajectories is established thanks to a comparison principle.The difficulty here is that we are dealing with a system of coupled equations and we need to define a specific order for which such a comparison principle holds.Now in [19,Theorem 4.1], it has been proved that at the microscopic level, the generalised contact process is attractive only for the following order: Note that in the corresponding state space Σ N , the order above translates into (0, 1) < (0, 0) < (1, 1) < (1, 0).
Attractiveness for the order (4.1) means that given two configurations η ≤ ∼ η, it is possible to build a coupling between (η t ) t≥0 and ( ∼ η t ) t≥0 where both these processes evolve according to the dynamics given by (2.11), such that η 0 ≤ ∼ η 0 and almost surely, for all t ≥ 0, η t ≤ ∼ η t pointwise in the sense of (4.1).Note that using [4, Theorem 2.4], one can show that the system remains attractive when adding an exchange and reservoir dynamics.It is then natural to think that attractiveness also holds at the macroscopic level through a comparison principle.A comparison principle means that if two profiles are such that at a certain time, one is smaller than the other almost everywhere, then the same is true at any later time.Considering the microscopic order (4.1) it is natural to consider that the largest state at the macroscopic level corresponds to (ρ 1 = 1, ρ 2 = 0, ρ 3 = 0) and the smallest state to (ρ 2 = 1, ρ 1 = ρ 3 = 0).We will work under the following change of coordinates: which is consistent with the fact that (1, 1, 1) corresponds to the largest profile (ρ 1 = 1, ρ 2 = 0, ρ 3 = 0) and (0, 0, 0) with the lowest one (ρ 2 = 1, ρ 1 = ρ 3 = 0).In the sequel, we will say that given two profiles ρ and φ, ρ ≤ φ if: almost everywhere.Note that this new order adapted at the microscopic level, i.e is not equivalent to the one given in (4.1).Indeed, consider the configuration η full of 3's and η full of 0's.Then η ≤ η for the order (4.1) but not for the order (4.4).However, one can check that if η ≤ η for the order (4.4), then η ≤ η for the order (4.1), so the macroscopic order is consistent with the microscopic one but it is weaker, so we can compare fewer profiles.However, the notable fact, which we will prove, is that we have monotonicity under this new order, i.e. a comparison principle holds under the change of coordinates (4.2).To prove that, as previously, since we are working in any dimension d ≥ 1 with mixed boundary conditions, some care must be taken to deal with the integral terms on Γ.For that, we strongly rely on analytical tools stated in [24].
Under the change of coordinates (4.2), the coupled equations in the bulk become, : We will see that the comparison principle stated and proved in Lemma 3 yields the following Theorem which is used to prove Theorem 2.
Theorem 6. Suppose that conditions (H 1 ) hold.Then, there exists a unique stationary solution ρ D,R , resp.ρ Ne,R of (2.23), resp.(2.27).Furthermore, for any solution ρ D,R , resp.ρ Ne,R to the boundary value problem Note that this result can be equivalently formulated in the change of coordinates (4.2) and we will prove it in that setting in the next subsection.Remark 4. One could ask if conditions on the parameters are necessary to establish existence and uniqueness of the stationary solution of the hydrodynamic equation.Could we not generalize the result to all parameters?
In order to answer that, we simulated the solutions to the equation in the (Neumann ; Neumann) regime for which the constant profile (ρ 1 = 0, ρ 2 = r r+1 , ρ 3 = 0) is stationary.Indeed, , 0 = 0 and it corresponds to the extinction regime, that is, there are no more wild insects.We observed (see below in the Appendix B) that in dimensions 1, for parameters λ 1 = 1, λ 2 = 0.75 and D = r = 1, for which conditions (H 1 ) are not satisfied, the solution of the hydrodynamic equation starting from ρ 1 = 1, ρ 2 = ρ 3 = 0 converges to a constant profile which is not (0, r r+1 , 0) so uniqueness does not hold.Simulations confirm that Theorem 6 does not hold in all generality and that conditions on the parameters are necessary, although conditions (H 1 ) might not be the optimal ones.

Proof of the hydrostatic limit
Let us prove Theorem 2. We prove the first point, the second one follows in the same way.Denote A T ⊂ D([0, T ], M + 3 ) the set of trajectories { ρ(t, u)du, 0 ≤ t ≤ T } whose density ρ = (ρ 1 , ρ 2 , ρ 3 ) satisfies conditions (2.24) and (2.26) of the definition of a weak solution of (2.23) for some initial profile ρ 0 .Consider Q * ss ( θ) a limit point of the sequence (Q N µ ss N ( θ) ) N ≥1 associated to the invariant measures.By Theorem 1, and Then, one concludes thanks to (4.7) in Theorem 6 and dominated convergence theorem.

Proof of Theorem 6
In order to prove Theorem 6 we first establish a comparison principle (Lemma 3).Then, we show that the difference between the largest solution and the smallest solution vanishes (Lemma 4).Using an integration by parts, it is useful to rewrite the weak formulations (2.25) and (2.25), in the following suitable forms: for any 0 Taking n → ∞ and using (3.50) and (3.52) using the W k s and W k s instead of the V k s and V k s we get L 2 and where we used the trace inequality (3.46) in the second inequality.Taking T → ∞, and using that D ≥ 1 we get that By Corollary 1, R 1 is almost surely decreasing and R 0 increasing therefore R 1 t −R 0 t is almost surely decreasing and the above inequality implies We are now left to show that lim We proceed following the same steps as for A n .To lighten notations we will not write the subscript t in the computations, when there is no confusion.Let us compute the second term. .
Integrating this between 0 and T and using the Cauchy-Scwharz inequality we are left with   Using Lemma 3 and the Cauchy-Schwarz inequality, we get: Integrating this between 0 and T and using the Cauchy-Scwharz inequality we are left with:     In the figures shown below (Figures 2 and 3), the x axis corresponds to the one dimensional space B = [0, 1] and the y axis is the space of values of the density profiles ρ 1 , ρ 1 + ρ 3 , 1 − (ρ 2 + ρ 3 ).

Figure 1 :
Figure 1: Mixed boundary conditions depending on the values of θ and θ r .The letters D, resp.R, resp.N e denote a Dirichlet, resp.Robin, resp.Neumann boundary condition.For instance (D ; N e ) denotes a left-hand side Dirichlet boundary condition and a right-hand side Neumann boundary condition.

Proposition 4 .
For any sequence of probability measures (µ N ) N ≥0 on Σ N , for any G ∈ C 1,2 ([0, T ] × B), for all i ∈ {1, 2, 3} and any t ∈ [0, T ], , x/N )(η εN i,s (x) − η i,s (x))ds = 0. (3.21) 24) where again, to get the last line, we used a Taylor expansion of G and with C G (N ) uniformly bounded in N by a constant C G depending only on G. Using (3.23) and (3.24), we get 3.30) with C 4 > 0, a constant that is uniform in N and ε.Taking B = N 4A G ∞ and putting together (3.26), (3.25) and (3.30) yields .39) Here we will only make use of(3.38) and the first relation in (3.39) since we need to replace elements in the bulk and the right-hand side boundary of the system to recover the weak formulation of the equation in the (Dirichlet; Robin) regime.For regimes where a replacement is needed on the left-hand side boundary, we use the second relation in(3.39)  in the same way.We may thus replace in (3.35) and (3.5), η εN i by π N i * u ε in the bulk and η εN i by π N i * u right ε at the right boundary.Therefore, for any δ > 0. lim sup ε→0 lim sup , 3,s (r) drds .(3.43)For the replacement of the convolution at the boundary we use the following result which follows from [9, Section 5.3]: for any H ∈ H 1 (B) lim ε→0 H * u right ε = Tr(H) a.s in Γ + .(3.44)

3. 5 . 1 Theorem 4 .
Uniqueness of the solution in the (Neumann ; Robin) mixed regime Here, we choose a basis of eigenvectors satisfying Neumann conditions on both boundaries (3.47).There exists a unique solution to the Neumann + Robin boundary problem (2.27).

3. 5 . 2 Theorem 5 .
Uniqueness of the solution in the (Dirichlet ; Robin) mixed regime Here, we choose a basis of eigenvectors satisfying a Dirichlet boundary condition on the left and Neumann boundary condition on the right (3.61).There exists a unique solution to the Dirichlet + Robin boundary problem (2.23).