# Switching interacting particle systems:
scaling limits, uphill diffusion and boundary layer

Simone Floreani1, Cristian Giardinà2, Frank den Hollander3,
Shubhamoy Nandan4, Frank Redig5

(August 27, 2024)

###### Abstract

In this paper we consider three classes of interacting particle systems on
${\mathbb{Z}}$: independent random walks, the exclusion process, and the
inclusion process. We allow particles to switch their jump rate (the rate
identifies the type of particle) between $1$ (fast particles) and
$\epsilon\in[0,1]$ (slow particles). The switch between the two jump rates
happens at rate $\gamma\in(0,\infty)$. In the exclusion process, the
interaction is such that each site can be occupied by at most one particle of
each type. In the inclusion process, the interaction takes places between
particles of the same type at different sites and between particles of
different type at the same site.

We derive the macroscopic limit equations for the three systems, obtained
after scaling space by $N^{-1}$, time by $N^{2}$, the switching rate by
$N^{-2}$, and letting $N\to\infty$. The limit equations for the macroscopic
densities associated to the fast and slow particles is the well-studied double
diffusivity model. This system of reaction-diffusion equations was introduced
to model polycrystal diffusion and dislocation pipe diffusion, with the goal
to overcome the limitations imposed by Fick’s law. In order to investigate the
microscopic out-of-equilibrium properties, we analyse the system on
$[N]=\\{1,\ldots,N\\}$, adding boundary reservoirs at sites $1$ and $N$ of
fast and slow particles, respectively. Inside $[N]$ particles move as before,
but now particles are injected and absorbed at sites $1$ and $N$ with
prescribed rates that depend on the particle type. We compute the steady-state
density profile and the steady-state current. It turns out that uphill
diffusion is possible, i.e., the total flow can be in the direction of
increasing total density. This phenomenon, which cannot occur in a single-type
particle system, is a violation of Fick’s law made possible by the switching
between types. We rescale the microscopic steady-state density profile and
steady-state current and obtain the steady-state solution of a boundary-value
problem for the double diffusivity model.

_Keywords:_ Switching random walks, fast and slow particles, duality, scaling
limits, uphill diffusion, Fick’s law.

_MSC2020:_ Primary: 60J70 60K35. Secondary: 82C26 92D25.

_Acknowledgement:_ The research in this paper was supported by the Netherlands
Organisation for Scientific Research (NWO) through grant TOP1.17.019. S.F.
thanks Antonio Agresti for several enlightening discussions and Mark Veraar
for useful suggestions.

1) Delft Institute of Applied Mathematics, TU Delft, Delft, The Netherlands,
<EMAIL_ADDRESS>
2) Modena and Reggio Emilia University, Modena, Italy,
<EMAIL_ADDRESS>
3) Mathematical Institute, Leiden University, Leiden, The Netherlands,
<EMAIL_ADDRESS>
4) Mathematical Institute, Leiden University, Leiden, The Netherlands,
<EMAIL_ADDRESS>
5) Delft Institute of Applied Mathematics, TU Delft, Delft, The Netherlands,
<EMAIL_ADDRESS>
##  1 Introduction

Section 1.1 provides the background and the motivation for the paper. Section
1.2 defines the model. Section 1.3 identifies the dual and the stationary
measures. Section 1.4 gives a brief outline of the remainder of the paper.

###  1.1 Background and motivation

Interacting particle systems are used to model and analyse properties of _non-
equilibrium systems_ , such as macroscopic profiles, long-range correlations
and macroscopic large deviations. Some models have additional structure, such
as duality or integrability properties, which allow for a study of the fine
details of non-equilibrium steady states, such as microscopic profiles and
correlations. Examples include zero-range processes, exclusion processes, and
models that fit into the algebraic approach to duality, such as inclusion
processes and related diffusion processes, or models of heat conduction, such
as the Kipnis-Marchioro-Presutti model [9, 21, 22, 30, 36]. Most of these
models have indistinguishable particles of which the total number is
conserved, and so the relevant macroscopic quantity is the _density_ of
particles.

Turning to more complex models of non-equilibrium, various exclusion processes
with _multi-type particles_ have been studied [24, 25, 39], as well as
reaction-diffusion processes [7, 8, 19, 17, 18], where non-linear reaction-
diffusion equations are obtained in the hydrodynamic limit, and large
deviations around such equations have been analysed. In the present paper, we
focus on a reaction-diffusion model that on the one hand is simple enough so
that via duality a complete microscopic analysis of the non-equilibrium
profiles can be carried out, but on the other hand exhibits interesting
phenomena, such as _uphill diffusion_ and _boundary-layer effects_. In our
model we have two types of particles, fast and slow, that jump at rate $1$ and
$\epsilon\in[0,1]$, respectively. Particles of identical type are allowed to
interact via exclusion or inclusion. There is no interaction between particles
of different type that are at different sites. Each particle can change type
at a rate that is adapted to the particle interaction (exclusion or
inclusion), and is therefore interacting with particles of different type at
the same site. An alternative and equivalent view is to consider two layers of
particles, where the layer determines the jump rate (rate $1$ for bottom
layer, rate $\epsilon$ for top layer) and where on each layer the particles
move according to exclusion or inclusion, and to let particles change layer at
a rate that is appropriately chosen in accordance with the interaction. In the
limit as $\epsilon\downarrow 0$, particles are immobile on the top layer.

We show that the _hydrodynamic limit_ of all three dynamics is a linear
reaction-diffusion system known under the name of double diffusivity model,
namely,

(1.1)
$\begin{cases}\partial_{t}\rho_{0}=\Delta\rho_{0}+\Upsilon(\rho_{1}-\rho_{0}),\\\
\partial_{t}\rho_{1}=\epsilon\Delta\rho_{1}+\Upsilon(\rho_{0}-\rho_{1}),\end{cases}$

where $\rho_{i}$, $i\in\\{0,1\\}$, are the macroscopic densities of the two
types of particles, and $\Upsilon\in(0,\infty)$ is the scaled switching rate.
The above system was introduced in [1] to model polycrystal diffusion (more
generally, diffusion in inhomogeneous porous media) and dislocation pipe
diffusion, with the goal to overcome the restrictions imposed by Fick’s law.
Non-Fick behaviour is immediate from the fact that the total density
$\rho=\rho_{0}+\rho_{1}$ does not satisfy the classical diffusion equation.

The double diffusivity model was studied extensively in the PDE literature [2,
33, 35], while its discrete counterpart was analysed in terms of a single
random walk switching between two layers [34]. The same macroscopic model was
studied independently in the mathematical finance literature in the context of
switching diffusion processes [46]. Thus, we have a family of interacting
particle systems whose macroscopic limit is relevant in several contexts.
Another context our three dynamics fit into are models of interacting active
random walks with an internal state that changes randomly (e.g. activity,
internal energy) and that determines their diffusion rate and or drift [15,
28, 32, 37, 41, 43, 3].

An additional motivation to study two-layer models comes from population
genetics. Individuals live in colonies, carry different genetics types, and
can be either active or dormant. While active, individuals resample by
adopting the type of a randomly sampled individual in the same colony, and
migrate between colonies by hopping around. Active individuals can become
dormant, after which they suspend resampling and migration, until they become
active again. Dormant individuals reside in what is called a _seed bank_. The
overall effect of dormancy is that extinction of types is slowed down, and so
genetic diversity is enhanced by the presence of the seed bank. A wealth of
phenomena can occur, depending on the parameters that control the rates of
resampling, migration, falling asleep and waking up [6, 31]. Dormancy not only
affects the long-term behaviour of the population quantitatively. It may also
lead to qualitatively different equilibria and time scales of convergence. For
a panoramic view on the role of dormancy in the life sciences, we refer the
reader to [40].

From the point of view of non-equilibrium systems driven by boundary
reservoirs, switching interacting particle systems have not been studied. On
the one hand, such systems have both reaction and diffusion and therefore
exhibit a richer non-equilibrium behaviour. On the other hand, the macroscopic
equations are linear and exactly solvable in one dimension, and so these
systems are simple enough to make a detailed microscopic analysis possible. As
explained above, the system can be viewed as an interacting particle system on
two layers. Therefore duality properties are available, which allows for a
detailed analysis of the system coupled to reservoirs, dual to an absorbing
system. In one dimension the analysis of the microscopic density profile
reduces to a computation of the absorption probabilities of a simple random
walk on a two-layer system absorbed at the left and right boundaries. From the
analytic solution, we can identify both the density profile and the current in
the system. This leads to two interesting phenomena. The first phenomenon is
_uphill diffusion_ (see e.g. [12, 16, 13, 20, 38]), i.e., in a well-defined
parameter regime the current can go against the particle density gradient:
when the total density of particles at the left end is higher than at the
right end, the current can still go from right to left. The second phenomenon
is _boundary-layer behaviour_ : in the limit as $\epsilon\downarrow 0$, in the
macroscopic stationary profile the densities in the top and bottom layer are
equal, which for unequal boundary conditions in the top and bottom layer
results in a _discontinuity_ in the stationary profile. Corresponding to this
jump in the macroscopic system, we identify a boundary layer of size
$\sqrt{\epsilon}\,\log(1/\epsilon)$ in the microscopic system where the
densities are unequal. The quantification of the _size_ of this boundary layer
is an interesting corollary of the exact macroscopic stationary profile that
we obtain from the microscopic system via duality.

###  1.2 Three models

For $\sigma\in\\{-1,0,1\\}$ we introduce an interacting particle system on
${\mathbb{Z}}$ where the particles randomly switch their jump rate between two
possible values, $1$ and $\epsilon$, with $\epsilon\in[0,1]$. For $\sigma=-1$
the particles are subject to the exclusion interaction, for $\sigma=0$ the
particles are independent, while for $\sigma=1$ the particles are subject to
the inclusion interaction. Let

$\displaystyle\eta_{0}(x)$ $\displaystyle:=\text{number of particles at site
}x\text{ jumping at rate }1,$ $\displaystyle\eta_{1}(x)$
$\displaystyle:=\text{number of particles at site }x\text{ jumping at rate
}\epsilon.$

The configuration of the system is

$\eta:=\\{\eta(x)\\}_{x\in{\mathbb{Z}}}\in\mathcal{X}=\begin{cases}\\{0,1\\}^{\mathbb{Z}}\times\\{0,1\\}^{\mathbb{Z}},&\text{if
}\sigma=-1,\\\
{\mathbb{N}}_{0}^{{\mathbb{Z}}}\times{\mathbb{N}}_{0}^{{\mathbb{Z}}},&\text{if
}\sigma=0,1,\end{cases}$

where

$\eta(x):=(\eta_{0}(x),\eta_{1}(x)),\qquad x\in{\mathbb{Z}}.$

We call $\eta_{0}=\\{\eta_{0}(x)\\}_{x\in{\mathbb{Z}}}$ and
$\eta_{1}=\\{\eta_{1}(x)\\}_{x\in{\mathbb{Z}}}$ the configurations of fast
particles, respectively, slow particles. When $\epsilon=0$ we speak of dormant
particles (see Fig. 1).

(a) Representation via _slow_ and _fast_ particles moving on the one-layer
graph ${\mathbb{Z}}$.

(b) Representation via particles moving on the two-layer graph
${\mathbb{Z}}\times I$.

Figure 1: Two equivalent representations of switching independent random walks
($\sigma=0$).

######  Definition 1.1.

[Switching interacting particle systems] For $\epsilon\in[0,1]$ and
$\gamma\in(0,\infty)$, let $L_{\epsilon,\gamma}$ be the generator

(1.2) $\displaystyle L_{\epsilon,\gamma}:=L_{0}+\epsilon L_{1}+\gamma
L_{0\updownarrow 1},$

acting on bounded cylindrical functions $f\colon\mathcal{X}\to{\mathbb{R}}$ as

$\displaystyle(L_{0}f)(\eta)=\sum_{|x-y|=1}\Big{\\{}\eta_{0}(x)(1+\sigma\eta_{0}(y))\,\big{[}f((\eta_{0}-\delta_{x}+\delta_{y},\eta_{1}))-f(\eta)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad+\eta_{0}(y)(1+\sigma\eta_{0}(x))\,\big{[}f((\eta_{0}+\delta_{x}-\delta_{y},\eta_{1}))-f(\eta)\big{]}\Big{\\}},$
$\displaystyle(L_{1}f)(\eta)=\sum_{|x-y|=1}\Big{\\{}\eta_{1}(x)(1+\sigma\eta_{1}(y))\,\big{[}f((\eta_{0},\eta_{1}-\delta_{x}+\delta_{y}))-f(\eta)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad+\eta_{1}(y)(1+\sigma\eta_{1}(x))\,\big{[}f((\eta_{0},\eta_{1}+\delta_{x}-\delta_{y}))-f(\eta)\big{]}\Big{\\}},$
$\displaystyle(L_{0\updownarrow
1}f)(\eta)=\gamma\sum_{x\in{\mathbb{Z}}^{d}}\Big{\\{}\eta_{0}(x)(1+\sigma\eta_{1}(x))\,\big{[}f((\eta_{0}-\delta_{x},\eta_{1}+\delta_{x}))-f(\eta)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad+\eta_{1}(x)(1+\sigma\eta_{0}(x))\,\big{[}f((\eta_{0}+\delta_{x},\eta_{1}-\delta_{x}))-f(\eta)\big{]}\Big{\\}}.$

The Markov process $\\{\eta(t)\colon\,t\geq 0\\}$ on state space $\mathcal{X}$
with

$\eta(t):=\\{\eta(x,t)\\}_{x\in{\mathbb{Z}}}=\big{\\{}(\eta_{0}(x,t),\eta_{1}(x,t))\big{\\}}_{x\in{\mathbb{Z}}},$

hopping rates $1,\epsilon$ and switching rate $\gamma$ is called switching
exclusion process for $\sigma=-1$, switching random walks for $\sigma=0$ (see
Fig. 1), and switching inclusion process for $\sigma=1$. $\spadesuit$

###  1.3 Duality and stationary measures

The systems defined in (1.2) can be equivalently formulated as jump processes
on the graph (see Fig. 1) with vertex set $\\{(x,i)\in{\mathbb{Z}}^{d}\times
I\\}$, with $I=\\{0,1\\}$ labelling the two layers, and edge set given by the
nearest-neighbour relation

$(x,i)\sim(y,j)\quad\text{ when }\quad\begin{cases}|x-y|=1\text{ and }i=j,\\\
x=y\text{ and }|i-j|=1.\end{cases}$

In this formulation the particle configuration is

$\eta=(\eta_{i}(x))_{(x,i)\in{\mathbb{Z}}\times I}$

and the generator $L$ is given by

(1.3) $\displaystyle(Lf)(\eta)=\sum_{i\in
I}\sum_{|x-y|=1}\epsilon^{i}\eta_{i}(x)(1+\sigma\eta_{i}(y))\,[f(\eta-\delta_{(x,i)}+\delta_{(y,i)})-f(\eta)]$
$\displaystyle\qquad\qquad\qquad\quad+\epsilon^{i}\eta_{i}(y)(1+\sigma\eta_{i}(x))\,[f(\eta-\delta_{(y,i)}+\delta_{(x,i)})-f(\eta)]$
$\displaystyle\qquad\qquad+\sum_{i\in
I}\gamma\sum_{x\in{\mathbb{Z}}}\eta_{i}(x)(1+\sigma\eta_{1-i})\,[f(\eta-\delta_{(x,i)}+\delta_{(x,1-i)})-f(\eta)].$

Thus, a single particle (when no other particles are present) is subject to
two movements:

  * i)

Horizontal movement: In layer $i=0$ and $i=1$ the particle performs a nearest-
neighbour random walk on ${\mathbb{Z}}$ at rate $1$, respectively, $\epsilon$.

  * ii)

Vertical movement: The particle switches layer at the same site at rate
$\gamma$.

It is well known (see e.g. [44]) that for these systems there exists a one-
parameter family of reversible product measures

$\big{\\{}\mu_{\theta}=\otimes_{(x,i)\in{\mathbb{Z}}\times
I}\nu_{(x,i),\theta}\colon\,\theta\in\Theta\big{\\}}$

with $\Theta=[0,1]$ if $\sigma=-1$ and $\Theta=[0,\infty)$ if
$\sigma\in\\{0,1\\}$, and with marginals given by

(1.4)
$\displaystyle\nu_{(x,i),\theta}=\begin{dcases}\text{Bernoulli}\,(\theta),&\sigma=-1,\\\\[4.26773pt]
\text{Poisson}\,(\theta),&\sigma=0,\\\\[4.26773pt] \text{Negative--
Binomial}\,(1,\tfrac{\theta}{1+\theta}),&\sigma=1.\end{dcases}$

Moreover, the classical self-duality relation holds, i.e., for all
configurations $\eta,\,\xi\in\mathcal{X}$ and for all times $t\geq 0$,

${\mathbb{E}}_{\eta}[D(\xi,\eta_{t})]={\mathbb{E}}_{\xi}[D(\xi_{t},\eta)],$

with $\\{\xi(t):\ t\geq 0\\}$ and $\\{\eta(t):\ t\geq 0\\}$ two copies of the
process with generator given in (1.2) and self-duality function
$D\colon\,\mathcal{X}\times\mathcal{X}\to{\mathbb{R}}$ given by

(1.5) $D(\xi,\eta):=\prod_{(x,i)\in{\mathbb{Z}}^{d}\times
I}d(\xi_{i}(x),\eta_{i}(x)),$

with

(1.6) $d(k,n):=\frac{n!}{(n-k)!}\frac{1}{w(k)}\,\mathbf{1}_{\\{k\leq n\\}}$

and

(1.7) $w(k):=\begin{cases}\frac{\Gamma(1+k)}{\Gamma(1)},&\sigma=1,\\\
1,&\sigma=-1,0.\end{cases}$

######  Remark 1.2.

[Possible extensions] Note that we could allow for more than two layers,
inhomogeneous rates and non-nearest neighbour jumps as well, and the same
duality relation would still hold (see e.g. [26] for an inhomogeneous version
of the exclusion process). More precisely, let
$\\{\omega_{i}(\\{x,y\\})\\}_{x,y\in{\mathbb{Z}}}$ and
$\\{\alpha_{i}(x)\\}_{x\in{\mathbb{Z}}}$ be collections of bounded weights for
$i\in I_{M}=\\{0,1,\ldots,M\\}$ with $M<\infty$. Then the interacting particle
systems with generator

(1.8) $\displaystyle(L_{D,\gamma}f)(\eta)$
$\displaystyle=\sum_{i=0}^{M}D_{i}\sum_{|x-y|=1}\omega_{i}(\\{x,y\\})\,\Big{\\{}\eta_{i}(x)\,(\alpha_{i}(y)+\sigma\eta_{i}(y))\,\big{[}f(\eta-\delta_{(x,i)}+\delta_{(y,i)})-f(\eta)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\eta_{i}(y)\,(\alpha_{i}(x)+\sigma\eta_{i}(x))\,\big{[}f(\eta-\delta_{(y,i)}+\delta_{(x,i)})-f(\eta)\big{]}\Big{\\}}$
$\displaystyle\quad+\sum_{i=0}^{M-1}\gamma_{\\{i,i+1\\}}\,\sum_{x\in{\mathbb{Z}}}\Big{\\{}\eta_{i}(x)\,\big{[}f(\eta-\delta_{(x,i)}+\delta_{(x,i+1)})-f(\eta)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad+\eta_{i+1}(x)\,\big{[}f(\eta-\delta_{(x,i+1)}+\delta_{(x,i)})-f(\eta)\big{]}\Big{\\}},$

with $\eta=(\eta_{i}(x))_{(x,i)\in{\mathbb{Z}}\times I_{M}}$,
$\\{D_{i}\\}_{i\in I_{M}}$ a bounded decreasing collection of weights in
$[0,1]$ and $\gamma_{\\{i,i+1\\}}\in(0,\infty)$, are still self-dual with
duality function as in (1.5), but with $I$ replaced by $I_{M}$ and single-site
duality functions given by
$d_{(x,i)}(k,n)=\frac{n!}{(n-k)!}\frac{1}{w_{(x,i)}(k)}\,\mathbf{1}_{\\{k\leq
n\\}}$ with

$\displaystyle w_{(x,i)}(k)\
:=\begin{dcases}\frac{\alpha_{i}(x)!}{(\alpha_{i}(x)-k)!}{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\\{k\leq\alpha_{i}(x)\\}},&\sigma=-1,\\\
\alpha_{i}(x)^{k},&\sigma=0,\\\
\frac{\Gamma(\alpha_{i}(x)+k)}{\Gamma(\alpha_{i}(x))},&\sigma=1.\end{dcases}$

In the present paper we prefer to stick to the two-layer homogeneous setting
in order not to introduce extra notations. However, it is straightforward to
extend many of our results to the inhomogeneous multi-layer model.
$\spadesuit$

###  1.4 Outline

Section 2 identifies and analyses the hydrodynamic limit of the system in
Definition 1.1 after scaling space, time and switching rate diffusively. We
thereby exhibit a class of interacting particle systems whose microscopic
dynamics scales to a macroscopic dynamics called the double diffusivity model.
Moreover, we provide a discussion on the solutions of this model, connecting
mathematical literature applied to material science and to financial
mathematics. Section 3 looks at what happens, both microscopically and
macroscopically, when boundary reservoirs are added, resulting in a non-
equilibrium flow. Here the possibility of uphill diffusion becomes manifest,
which is absent in single-layer systems, i.e., the two layers interact in a
way that allows for a violation of Fick’s law. We characterise the parameter
regime for uphill diffusion. Moreover, we show that, in the limit as
$\epsilon\downarrow 0$, the macroscopic stationary profile of the type-1
particles adapts to the microscopic stationary profile of the type-0
particles, resulting in a _discontinuity_ at the boundary in the case of
unequal boundary conditions on the top layer and the bottom layer. Appendix A
provides the inverse of a certain boundary-layer matrix.

##  2 The hydrodynamic limit

In this section we scale space, time and switching diffusively, so as to
obtain a hydrodynamic limit. In Section 2.1 we scale space by $1/N$, time by
$N^{2}$, the switching rate by $1/N^{2}$, introduce scaled microscopic
empirical distributions, and let $N\to\infty$ to obtain a system of
macroscopic equations. In Section 2.2 we recall some known results for this
system, namely, there exists a unique solution that can be represented in
terms of an underlying diffusion equation or, alternatively, via a Feynman-Kac
formula involving the switching diffusion process.

###  2.1 From microscopic to macroscopic

Let $N\in{\mathbb{N}}$, and consider the scaled generator
$L_{\epsilon,\gamma_{N}}$ (recall (1.2)) with $\gamma_{N}=\Upsilon/N^{2}$ for
some $\Upsilon\in(0,\infty)$, i.e., the reaction term is slowed down by a
factor $N^{2}$ in anticipation of the diffusive scaling we are going to
consider.

In order to study the collective behaviour of the particles after scaling of
space and time, we introduce the following empirical density fields, which are
Radon measure-valued càdlàg processes:

$\mathsf{X}^{N}_{0}(t):=\frac{1}{N}\sum_{x\in{\mathbb{Z}}}\eta_{0}(x,tN^{2})\,\delta_{x/N},\qquad\mathsf{X}^{N}_{1}(t):=\frac{1}{N}\sum_{x\in{\mathbb{Z}}}\eta_{1}(x,tN^{2})\,\delta_{x/N}.$

In order to derive the hydrodynamic limit for the switching interacting
particle systems, we need the following set of assumptions. In the following
we denote by $C^{\infty}_{c}({\mathbb{R}})$ the space of infinitely
differentiable functions with values in ${\mathbb{R}}$ and compact support, by
$C_{b}({\mathbb{R}};\sigma)$ the space of bounded and continuous functions
with values in ${\mathbb{R}}_{+}$ for $\sigma\in\\{0,1\\}$ and with values in
$[0,1]$ for $\sigma=-1$, by $C_{0}({\mathbb{R}})$ the space of continuous
functions vanishing at infinity, by $C^{2}_{0}({\mathbb{R}})$ the space of
twice differentiable functions vanishing at infinity and by $M$ the space of
Radon measure on ${\mathbb{R}}$.

######  Assumption 2.1.

[Compatible initial conditions] Let $\bar{\rho}_{i}\in
C_{b}({\mathbb{R}};\sigma)$ for $i\in\\{0,1\\}$ be two given functions, called
initial macroscopic profiles. We say that a sequence
$(\mu_{N})_{N\in{\mathbb{N}}}$ of measures on $\mathcal{X}$ is a sequence of
compatible initial conditions when:

  * (i)

For any $i\in\\{0,1\\}$, $g\in C^{\infty}_{c}({\mathbb{R}})$ and $\delta>0,$

$\lim_{N\to\infty}\mu_{N}\left(\left|\langle\mathsf{X}^{N}_{0}(t),g\rangle-\int_{\mathbb{R}}{\rm
d}x\ \bar{\rho}_{i}(x)g(x)\right|>\delta\right)=0.$

  * (ii)

There exists a constant $C<\infty$ such that

(2.1) $\sup_{(x,i)\in{\mathbb{Z}}\times
I}{\mathbb{E}}_{\mu_{N}}[\eta_{i}(x)^{2}]\leq C.$

$\spadesuit$

Note that Assumption 2.1(ii) is the same as employed in [11, Theorem 1,
Assumption (b)] and is trivial for the exclusion process.

######  Theorem 2.2.

[Hydrodynamic scaling] Let $\bar{\rho}_{0},\bar{\rho}_{1}\in
C_{b}({\mathbb{R}};\sigma)$ be two initial macroscopic profiles, and let
$(\mu_{N})_{N\in{\mathbb{N}}}$ be a sequence of compatible initial conditions.
Let ${\mathbb{P}}_{\mu_{N}}$ be the law of the measure-valued process

$\\{X^{N}(t)\colon\,t\geq 0\\},\qquad X^{N}(t):=(X_{0}^{N}(t),X_{1}^{N}(t)),$

induced by the initial measure $\mu_{N}$. Then, for any $T,\delta>0$ and $g\in
C^{\infty}_{c}({\mathbb{R}})$,

$\lim_{N\to\infty}{\mathbb{P}}_{\mu_{N}}\left(\sup_{t\in[0,T]}\left|\,\langle\mathsf{X}^{N}_{i}(t),g\rangle-\int_{\mathbb{R}}{\rm
d}x\,\rho_{i}(x,t)g(x)\,\right|>\delta\right)=0,\qquad i\in I,$

where $\rho_{0}$ and $\rho_{1}$ are the unique continuous and bounded strong
solutions of the system

(2.2)
$\begin{cases}\partial_{t}\rho_{0}=\Delta\rho_{0}+\Upsilon(\rho_{1}-\rho_{0}),\\\
\partial_{t}\rho_{1}=\epsilon\Delta\rho_{1}+\Upsilon(\rho_{0}-\rho_{1}),\end{cases}$

with initial conditions

(2.3) $\begin{cases}\rho_{0}(x,0)=\bar{\rho}_{0}(x),\\\
\rho_{1}(x,0)=\bar{\rho}_{1}(x).\end{cases}$

###### Proof.

The proof follows the standard route presented in [45, Section 8] (see also
[19, 11]). We still explain the main steps because the two-layer setup is not
standard. First of all, note that the macroscopic equation (2.2) can be
straightforwardly identified by computing the action of the rescaled generator
$L^{N}=L_{\epsilon,\Upsilon/N^{2}}$ on the cylindrical functions
$f_{i}(\eta):=\eta_{i}(x),\ i\in\\{0,1\\}$, namely ,

$(L^{N}f_{i})(\eta)=\epsilon^{i}\left[\eta_{i}(x+1)-2\eta_{i}(x)+\eta_{i}(x-1)\right]+\frac{\Upsilon}{N^{2}}\left[\eta_{1-i}(x)-\eta_{i}(x)\right]$

and hence, for any $g\in C^{\infty}_{c}({\mathbb{R}})$,

$\displaystyle\int_{0}^{tN^{2}}{\rm
d}s\,L^{N}(\langle\mathsf{X}^{N}_{i}(s),g\rangle)$
$\displaystyle=\int_{0}^{tN^{2}}{\rm
d}s\,\frac{\epsilon^{i}}{N}\sum_{x\in{\mathbb{Z}}}\eta_{i}(x,s)\,\tfrac{1}{2}\,\big{[}g((x+1)/N)-2g(x/N)+g((x-1)/N)\big{]}$
$\displaystyle\quad+\int_{0}^{tN^{2}}{\rm
d}s\,\frac{1}{N}\sum_{x\in{\mathbb{Z}}}g(x/N)\,\frac{\Upsilon}{N^{2}}\,\big{[}\eta_{1-i}(x,s)-\eta_{i}(x,s)\big{]},$

where we moved the generator of the simple random walk to the test function by
using reversibility w.r.t. the counting measure. By the regularity of $g$, we
thus have

$\displaystyle\int_{0}^{tN^{2}}{\rm
d}s\,L^{N}(\langle\mathsf{X}^{N}_{i}(s),g\rangle)=\int_{0}^{t}{\rm
d}s\,\langle\mathsf{X}^{N}_{i}(sN^{2}),\epsilon\Delta
g\rangle+\int_{0}^{tN^{2}}{\rm
d}s\,\frac{\Upsilon}{N^{2}}\,\big{[}\langle\mathsf{X}^{N}_{1-i}(s),g\rangle-\langle\mathsf{X}^{N}_{i}(s),g\rangle\big{]}+o(\tfrac{1}{N^{2}}),$

which is the discrete counterpart of the weak formulation of the right-hand
side of (2.2), i.e., $\int_{0}^{t}{\rm d}s\int_{\mathbb{R}}{\rm
d}x\,\rho_{i}(x,s)\Delta g(x)+\Upsilon\int_{0}^{t}{\rm
d}s\int_{\mathbb{R}}{\rm d}x\,[\rho_{1-i}(x,s)-\rho_{i}(x,s)]g(x)$. Thus, as a
first step, we show that

$\lim_{N\to\infty}{\mathbb{P}}_{\mu_{N}}\left(\sup_{t\in[0,T]}\left|\langle\mathsf{X}^{N}_{i}(t),g\rangle-\langle\mathsf{X}^{N}_{i}(0),g\rangle-\int_{0}^{t}{\rm
d}s\,\langle\mathsf{X}^{N}_{i}(sN^{2}),\epsilon^{i}\Delta
g\rangle\right.\right.\\\ \left.\left.-\int_{0}^{tN^{2}}{\rm
d}s\,\frac{\Upsilon}{N^{2}}\,\big{[}\langle\mathsf{X}^{N}_{1-i}(s)-\mathsf{X}^{N}_{i}(s),g\rangle\big{]}\right|>\delta\right)=0.$

In order to prove the above convergence we employ the Dynkin’s formula for
Markov processes, which gives that the process defined as

$\displaystyle
M^{N}_{i}(g,t):=\langle\mathsf{X}^{N}_{i}(t),g\rangle-\langle\mathsf{X}^{N}_{i}(0),g\rangle-\int_{0}^{tN^{2}}{\rm
d}s\,L^{N}(\langle\mathsf{X}^{N}_{i}(s),g\rangle)$

is a martingale w.r.t. the natural filtration generated by the process
$\\{\eta_{t}\\}_{t\geq 0}$ and with predictable quadratic variation expressed
in terms of the carré du champ, i.e.,

$\langle M^{N}_{i}(g,t),M^{N}_{i}(g,t)\rangle=\int_{0}^{t}{\rm
d}s\,{\mathbb{E}}_{\mu_{N}}\left[\Gamma^{N}_{i}(g,s)\right]$

with

$\Gamma^{N}_{i}(g,s)=L^{N}\left(\langle\mathsf{X}^{N}_{i}(s),g\rangle\
\right)^{2}-\langle\mathsf{X}^{N}_{i}(s),g\rangle
L^{N}\left(\langle\mathsf{X}^{N}_{i}(s),g\rangle\ \right).$

We then have, by Chebyshev’s inequality and Doob’s martingale inequality,

$\displaystyle{\mathbb{P}}_{\mu_{N}}\left(\sup_{t\in[0,T]}\left|\langle\mathsf{X}^{N}_{0}(s),g\rangle-\langle\mathsf{X}^{N}_{0}(s),g\rangle-\int_{0}^{t}{\rm
d}s\,\langle\mathsf{X}^{N}_{0}(sN^{2}),\epsilon\Delta g\rangle\right.\right.$
(2.4)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\quad\qquad\qquad\qquad\left.\left.-\int_{0}^{tN^{2}}{\rm
d}s\,\frac{\Upsilon}{N^{2}}\,\big{[}\langle\mathsf{X}^{N}_{1}(s),g\rangle-\langle\mathsf{X}^{N}_{0}(s),g\rangle\big{]}\right|>\delta\right)$
$\displaystyle\leq\frac{1}{\delta^{2}}{\mathbb{E}}_{\mu_{N}}\left[\sup_{t\in[0,T]}\left|M^{N}_{i}(g,s)\right|^{2}\right]\leq\frac{4}{\epsilon^{2}}{\mathbb{E}}_{\mu_{N}}\left[\left|M^{N}_{i}(g,T)\right|^{2}\right]=\frac{4}{\epsilon^{2}}{\mathbb{E}}_{\mu_{N}}\left[\langle
M^{N}_{i}(g,T),M^{N}_{i}(g,T)\rangle^{2}\right]$ (2.5)
$\displaystyle=\frac{4}{\delta^{2}N^{2}}{\mathbb{E}}_{\mu_{N}}\left[\int_{0}^{N^{2}T}{\rm
d}s\,\sum_{x\in{\mathbb{Z}}^{d}}\eta_{i}(x,s)(1+\sigma\eta_{i}(x\pm
1,s))\left(g\left(\frac{x\pm
1}{N}\right)-g\left(\frac{x}{N}\right)\right)\right]$
$\displaystyle+\frac{4\Upsilon}{\delta^{2}N^{4}}{\mathbb{E}}_{\mu_{N}}\left[\int_{0}^{N^{2}T}{\rm
d}s\,\sum_{x\in{\mathbb{Z}}^{d}}(\eta_{i}(x,s)+\eta_{1-i}(x,s)+2\sigma\eta_{i}(x,s)\eta_{1-i}(x,s))g^{2}\left(\frac{x}{N}\right)\right],$

where in the last equality we explicitly computed the carré du champ. Let
$k\in{\mathbb{N}}$ be such that the support of $g$ is in $[-k,k]$. Then, by
the regularity of $g$, (2.1) is bounded by

(2.6) $\frac{4}{\delta^{2}N^{2}}(N^{2}T)(2k+1)N\frac{\|\nabla
g\|_{\infty}}{N^{2}}\sup_{x,\in{\mathbb{Z}},s\in[0,N^{2}T]}{\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,s)(1+\sigma\eta_{i}(x+1,s))\right]\\\
+\frac{4\Upsilon}{\delta^{2}N^{4}}(N^{2}T)(2k+1)N\|g\|_{\infty}\sup_{x,\in{\mathbb{Z}},s\in[0,N^{2}T]}{\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,s)+\eta_{1-i}(x,s)+2\sigma\eta_{i}(x,s)\eta_{1-i}(x,s)\right].$

We now show that, as a consequence of (2.1), for any $(x,i),\
(y,j)\in{\mathbb{Z}}\times I,$

(2.7) ${\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,s)\right]\leq
C,\quad{\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,s)\eta_{j}(y,s)\right]\leq C,$

from which we obtain

$\displaystyle{\mathbb{P}}_{\mu_{N}}\left(\sup_{t\in[0,T]}\left|\langle\mathsf{X}^{N}_{0}(s),g\rangle-\langle\mathsf{X}^{N}_{0}(s),g\rangle-\int_{0}^{t}{\rm
d}s\,\langle\mathsf{X}^{N}_{0}(sN^{2}),\epsilon\Delta g\rangle\right.\right.$
(2.8)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\quad\qquad\qquad\qquad\left.\left.-\int_{0}^{tN^{2}}{\rm
d}s\,\frac{\Upsilon}{N^{2}}\,\big{[}\langle\mathsf{X}^{N}_{1}(s),g\rangle-\langle\mathsf{X}^{N}_{0}(s),g\rangle\big{]}\right|>\delta\right)$
$\displaystyle\leq\frac{8T}{\delta^{2}N}(2k+1)\|\nabla
g\|_{\infty}C+\Upsilon\frac{16T}{\delta^{2}N}(2k+1)\|g\|_{\infty}C,$

and the desired convergence follows. In order to prove (2.7), first of all
note that, by the Cauchy-Schwartz inequality, it follows from (2.1) that, for
any $(x,i),\ (y,j)\in{\mathbb{Z}}\times I$,

(2.9) ${\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x)\eta_{j}(y)\right]\leq C.$

Moreover, recalling the duality functions given in (1.5) and defining the
configuration $\xi=\delta_{(x,i)}+\delta_{(y,j)}$ for $(x,i)\neq(y,j)$, we
have that

$\displaystyle{\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,t)\eta_{j}(y,t)\right]={\mathbb{E}}_{\mu_{N}}[D(\xi,\eta_{t})]=\int_{\mathcal{X}}{\mathbb{E}}_{\eta}[D(\xi,\eta_{t})]\,{\rm
d}\mu_{N}(\eta)$
$\displaystyle=\int_{\mathcal{X}}{\mathbb{E}}_{\xi}[D(\xi_{t},\eta)]\,{\rm
d}\mu_{N}(\eta)={\mathbb{E}}_{\xi}\left[{\mathbb{E}}_{\mu_{N}}[D(\xi_{t},\eta)]\right]$

and, labeling the particles in the dual configuration as $(X_{t},i_{t})$ and
$(Y_{t},j_{t})$ with initial conditions $(X_{0},i_{0})=(x,i)$ and
$(Y_{0},j_{0})=(y,j)$, we obtain

$\displaystyle{\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,t)\eta_{j}(y,t)\right]={\mathbb{E}}_{\xi}\left[{\mathbb{E}}_{\mu_{N}}[\eta_{i_{t}}(X_{t})\eta_{j_{t}}(Y_{t}){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{(X_{t},i_{t})\neq(Y_{t},j_{t})}]+{\mathbb{E}}_{\mu_{N}}[\eta_{i_{t}}(X_{t})(\eta_{i_{t}}(X_{t})-1){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{(X_{t},i_{t})=(Y_{t},j_{t})}]\right]$ (2.10)
$\displaystyle\leq{\mathbb{E}}_{\xi}\left[{\mathbb{E}}_{\mu_{N}}[\eta_{i_{t}}(X_{t})\eta_{j_{t}}(Y_{t})]\right]\leq{\mathbb{E}}_{\xi}\left[\sup_{(x,i),(y,j)\in{\mathbb{Z}}\times\\{0,1\\}}{\mathbb{E}}_{\mu_{N}}[\eta_{i}(x)\eta_{j}(y)]\right]\leq
C,$

where we used (2.9) in the last inequality. Similarly, for
$\xi=\delta_{(x,i)}$ and $(X_{t},i_{t})$ the dual particle with initial
condition $(X_{0},i_{0})=(x,i)$, we have that
${\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,t)\right]\leq{\mathbb{E}}_{\mu_{N}}\left[D(\xi,\eta_{t})\right]={\mathbb{E}}_{\xi}[{\mathbb{E}}_{\mu_{N}}[\eta_{i_{t}}(X_{t})]]$.
Using that $\eta_{i}(x)\leq\eta_{i}(x)^{2}$ for any
$(x,i)\in{\mathbb{Z}}\times I$ and using (2.1), we obtain (2.7).

The proof is concluded after showing the following:

  * i)

Tightness holds for the sequence of distributions of the processes
$\\{\mathsf{X}^{N}_{i}\\}_{N\in{\mathbb{N}}}$, denoted by
$\\{Q_{N}\\}_{N\in{\mathbb{N}}}$.

  * ii)

All limit points coincide and are supported by the unique path
$\mathsf{X}_{i}(t,{\rm d}x)=\rho_{i}(x,t)\,{\rm d}x$, with $\rho_{i}$ the
unique weak (and in particular strong) bounded and continuous solution of
(2.2).

While for (i) we provide an explanation, we skip the proof of (ii) because it
is standard and based on PDE arguments, namely, the existence and the
uniqueness of the solutions in the class of continuous-time functions with
values in $C_{b}({\mathbb{R}},\sigma)$ (we refer to [45, Lemma 8.6 and 8.7]
for further details), and the fact that Assumption 2.1(i) ensures that the
initial condition of (2.2) is also matched.

Tightness of the sequence $\\{Q_{N}\\}_{N\in{\mathbb{N}}}$ follows from the
compact containment condition on the one hand, i.e., for any $\delta>0$ and
$t>0$ there exists a compact set $K\subset M$ such that
${\mathbb{P}}_{\mu_{N}}(X^{N}_{i}\in K)>1-\delta$, and the equi-continuity
condition on the other hand, i.e.,
$\limsup_{N\to\infty}{\mathbb{P}}_{\mu_{N}}(\omega(\mathsf{X}_{i}^{N},\delta,T))\geq\mathfrak{e})\leq\mathfrak{e}$
for
$\omega(\alpha,\delta,T):=\sup\\{d_{M}(\alpha(s),\alpha(t)):s,t\in[0,T],|s-t|\leq\delta\\}$
with $d_{M}$ the metric on Radon measures defined as

$d_{M}(\nu_{1},\nu_{2}):=\sum_{j\in{\mathbb{N}}}2^{-j}\left(1\wedge\left|\int_{\mathbb{R}}\phi_{j}{\rm
d}\nu_{1}-\int_{\mathbb{R}}\phi_{j}{\rm d}\nu_{2}\right|\right)$

for an appropriately chosen sequence of functions
$(\phi_{j})_{j\in{\mathbb{N}}}$ in $C^{\infty}_{c}({\mathbb{R}})$. We refer to
[45, Section A.10] for details on the above metric and to the proof of [45,
Lemma 8.5] for the equi-continuity condition. We conclude by proving the
compact containment condition. Define

$K:=\Big{\\{}\nu\in M\ s.t.\ \exists\,k\in{\mathbb{N}}\ s.t.\
\nu[-\ell,\ell]\leq A(2\ell+1)\ell^{2}\ \forall\
\ell\in[k,\infty]\cap{\mathbb{N}}\Big{\\}}$

with $A>0$ such that $\frac{C\pi}{6A}<\epsilon$. By [45, Proposition A.25], we
have that $K$ is a pre-compact subset of $M$. Moreover, by the Markov
inequality and Assumption 2.1(ii), it follows that

$\displaystyle
Q_{N}(\bar{K}^{c})\leq\sum_{\ell\in{\mathbb{N}}}{\mathbb{P}}_{\mu_{N}}\left(\mathsf{X}_{i}^{N}([-\ell,\ell])\geq
A(2\ell+1)\ell^{2}\right)\leq\sum_{\ell\in{\mathbb{N}}}\frac{1}{A(2\ell+1)\ell^{2}}{\mathbb{E}}_{\mu_{N}}\left[\mathsf{X}_{i}^{N}([-\ell,\ell])\right]$
$\displaystyle=\sum_{\ell\in{\mathbb{N}}}\frac{1}{A(2\ell+1)\ell^{2}}\sum_{x\in[-\ell,\ell]\cap\frac{{\mathbb{Z}}}{N}}{\mathbb{E}}_{\mu_{N}}\left[\eta_{i}(x,tN^{2})\right]\leq\sum_{\ell\in{\mathbb{N}}}\frac{1}{A(2\ell+1)\ell^{2}}\frac{2\ell
N+1}{N}C\leq\frac{C}{A}\sum_{\ell\in{\mathbb{N}}}\frac{1}{\ell^{2}}<\epsilon,$

from which it follows that $Q_{N}(\bar{K})>1-\epsilon$ for any $N$. ∎

######  Remark 2.3.

[Total density] (i) If $\rho_{0},\rho_{1}$ are smooth enough and satisfy
(2.2), then by taking extra derivatives we see that the total density
$\rho:=\rho_{0}+\rho_{1}$ satisfies the thermal telegrapher equation

(2.11)
$\partial_{t}\left(\partial_{t}\rho+2\gamma\rho\right)=-\epsilon\Delta(\Delta\rho)+(1+\epsilon)\Delta\left(\partial_{t}\rho+\rho\right),$

which is second order in $\partial_{t}$ and fourth order in $\partial_{x}$
(see [2, 33] for a derivation). Note that (2.11) shows that the total density
does not satisfy the usual diffusion equation. This fact will be investigated
in detail in the next section, where we will analyse the non-Fick property of
$\rho$.
(ii) If $\epsilon=1$, then (2.11) simplifies to the heat equation
$\partial_{t}\rho=\Delta\rho$.
(iii) If $\epsilon=0$, then (2.11) reads

$\partial_{t}\left(\partial_{t}\rho+2\lambda\rho\right)=\Delta\left(\partial_{t}\rho+\rho\right),$

which is known as the strongly damped wave equation. The term
$\partial_{t}(2\lambda\rho)$ is referred to as frictional damping, the term
$\Delta(\partial_{t}\rho)$ as Kelvin-Voigt damping (see [10]). $\spadesuit$

###  2.2 Existence, uniqueness and representation of the solution

The existence and uniqueness of a continuous-time solution
$(\rho_{0}(t),\,\rho_{1}(t))$ with values in $C_{b}({\mathbb{R}},\sigma)$ of
the system in (2.2) can be proved by standard Fourier analysis. Below we
recall some known results that have a more probabilistic interpretation.

##### Stochastic representation of the solution.

The system in (2.2) fits in the realm of switching diffusions (see e.g. [46]),
which are widely studied in the mathematical finance literature. Indeed, let
$\\{i_{t}\colon\,t\geq 0\\}$ be the pure jump process on state space
$I=\\{0,1\\}$ that switches at rate $\Upsilon$, whose generator acting on
bounded functions $g\colon\,I\to{\mathbb{R}}$ is

$(Ag)(i):=\Upsilon(g(1-i)-g(i)),\qquad i\in I.$

Let $\\{X_{t}\colon\,t\geq 0\\}$ be the stochastic process on ${\mathbb{R}}$
solving the stochastic differential equation

${\rm d}X_{t}=\psi(i_{t})\,{\rm d}W_{t},$

where $W_{t}=B_{2t}$ with $\\{B_{t}\colon\,t\geq 0\\}$ standard Brownian
motion, and $\psi\colon\,I\to\\{D_{0},D_{1}\\}$ is given by

$\psi:=D_{0}\,\boldsymbol{1}_{\\{0\\}}+D_{1}\,\boldsymbol{1}_{\\{1\\}},$

with $D_{0}=1$ and $D_{1}=\epsilon$ in our setting. Let
$\mathcal{L}=\mathcal{L}_{\epsilon,\Upsilon}$ be the generator defined by

$(\mathcal{L}f)(x,i):=\lim_{t\downarrow
0}\frac{1}{t}\,{\mathbb{E}}_{x,i}[f(X_{t},i_{t})-f(x,i)]$

for $f\colon\,{\mathbb{R}}\times I\to{\mathbb{R}}$ such that $f(\cdot,i)\in
C^{2}_{0}({\mathbb{R}})$. Then, via a standard computation (see e.g. [29,
Eq.(4.4)]), it follows that

$\displaystyle(\mathcal{L}f)(x,i)$ $\displaystyle=\psi(i)(\Delta
f)(x,i)+\Upsilon[f(x,1-i)-f(x,i)]$ $\displaystyle=\begin{cases}\Delta
f(x,0)+\Upsilon\,[f(x,1)-f(x,0)],&i=0,\\\ \epsilon\Delta
f(x,1)+\Upsilon\,[f(x,0)-f(x,1)],&i=1.\end{cases}$

We therefore have the following result that corresponds to [29, Chapter 5,
Section 4, Theorem 4.1](see also [46, Theorem 5.2]).

######  Theorem 2.4.

[Stochastic representation of the solution] Suppose that
$\bar{\rho}_{i}\colon\,{\mathbb{R}}\to{\mathbb{R}}$ for $i\in I$ are
continuous and bounded. Then (2.2) has a unique solution given by

$\rho_{i}(x,t)={\mathbb{E}}_{(x,i)}[\bar{\rho}_{i_{t}}(X_{t})],\qquad i\in I.$

Note that if there is only one particle in the system (1.2), then we are left
with a single random walk, say $\\{Y_{t}\colon\,t\geq 0\\}$, whose generator,
denoted by $\mathsf{A}$, acts on bounded functions
$f\colon\,{\mathbb{Z}}\times I\to{\mathbb{R}}$ as

$(\mathsf{A}f)(y,i)=\psi(i)\left[\sum_{z\sim
y}[f(z,i)-f(y,i)]\right]+\Upsilon\,[f(y,1-i)-f(y,i)].$

After we apply the generator to the function $f(y,i)=y$, we get

$\mathsf{(}Af)(y,i)=0,$

i.e., the position of the random walk is a martingale. Computing the quadratic
variation via the carré du champ, we find

$\mathsf{A}(Y_{t}^{2})=\psi(i_{t})[(Y_{t}+1)^{2}-Y_{t}^{2}]+\psi(i_{t})[(Y_{t}-1)^{2}-Y_{t}^{2}]=2\psi(i_{t}).$

Hence the predictable quadratic variation is given by

$\int_{0}^{t}{\rm d}s\,2\psi(i_{s}).$

Note that for $\epsilon=0$ the latter equals the total amount of time the
random walk is not dormant up to time $t$.

When we diffusively scale the system (scaling the reaction term was done at
the beginning of Section 2), the quadratic variation becomes

$\int_{0}^{tN^{2}}{\rm d}s\,\psi(i_{N,s})=\int_{0}^{t}{\rm d}r\,\psi(i_{r}).$

As a consequence, we have the following invariance principle:

  * Given the path of the process $\\{i_{t}\colon\,t\geq 0\\}$,

$\lim_{N\to\infty}\frac{Y_{N^{2}t}}{N}=W_{\int_{0}^{t}{\rm
d}r\,\sqrt{\psi(i_{r})}},$

where $W_{t}=B_{2t}$ with $\\{B_{t}\colon\,t\geq 0\\}$ is standard Brownian
motion.

Thus, if we knew the path of the process $\\{i_{r}\colon\,r\geq 0\\}$, then we
could express the solution of the system in (2.2) in terms of a time-changed
Brownian motion. However, even though $\\{i_{r}\colon\,r\geq 0\\}$ is a simple
flipping process, we cannot say much explicitly about the random time
$\int_{0}^{t}{\rm d}r\,\sqrt{\psi(i_{r})}$. We therefore look for a simpler
formula, where the relation to a Brownian motion with different velocities is
more explicit. We achieve this by looking at the resolvent of the generator
$\mathcal{L}$. In the following, we denote by $\\{S_{t},\,t\geq 0\\}$ the
semigroup on $C_{b}({\mathbb{R}})$ of $\\{W_{t}:\,t\geq 0\\}$.

######  Proposition 2.5.

[Resolvent] Let $f\colon\,{\mathbb{R}}\times I\to{\mathbb{R}}$ be a bounded
and smooth function. Then, for $\lambda>0$, $\epsilon\in(0,1]$ and $i\in I$,

(2.12) $\displaystyle(\lambda\boldsymbol{I}-\mathcal{L})^{-1}f(x,i)$
$\displaystyle=\int_{0}^{\infty}{\rm d}t\,\frac{1}{\epsilon^{i}}\,{\rm
e}^{-\tfrac{1+\epsilon}{\epsilon}\ell(\Upsilon,\lambda)t}\left(\cosh(tc_{\epsilon}(\Upsilon,\lambda))+\frac{1-\epsilon}{\epsilon}\ell_{\epsilon}(\Upsilon,\lambda)\frac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\lambda)}\right)\,(S_{t}f(\cdot,i))(x)$
$\displaystyle+\int_{0}^{\infty}{\rm d}t\,{\rm
e}^{-\tfrac{1+\epsilon}{\epsilon}\ell(\Upsilon,\lambda)t}\left(\frac{\Upsilon}{\epsilon}\sinh(tc_{\epsilon}(\Upsilon,\lambda))\right)\,(S_{t}f(\cdot,1-i))(x),$

where
$c_{\epsilon}(\Upsilon,\lambda)=\sqrt{\left(\tfrac{1-\epsilon}{\epsilon}\right)^{2}\ell(\Upsilon,\lambda)^{2}+\frac{\Upsilon^{2}}{\epsilon}}$
and $\ell(\Upsilon,\lambda)=\frac{\Upsilon+\lambda}{2}$, while for
$\epsilon=0$,

(2.13)
$\displaystyle(\lambda\boldsymbol{I}-\mathcal{L})^{-1}f(x,i)=\int_{0}^{\infty}{\rm
d}t\,{\rm
e}^{-\lambda\tfrac{2\Upsilon+\lambda}{\Upsilon+\lambda}t}\left(\left(\tfrac{\Upsilon}{\lambda+\Upsilon}\right)^{i}(S_{t}f(\cdot,0))(x)+\left(\tfrac{\Upsilon}{\Upsilon+\lambda}\right)^{i+1}(S_{t}f(\cdot,1))(x)\right).$

###### Proof.

The proof is split into two parts.

Case $\epsilon>0$. We can split the generator $\mathcal{L}$ as

$\mathcal{L}=\psi(i)\tilde{\mathcal{L}}=\psi(i)\left(\Delta+\frac{1}{\psi(i)}A\right)=\psi(i)(\Delta+\tilde{A}),$

i.e., we decouple $X_{t}$ and $i_{t}$ in the action of the generator. We can
now use the Feynman-Kac formula to express the resolvent of the operator
$\mathcal{L}$ in terms of the operator $\tilde{\mathcal{L}}$. Denoting by
$\tilde{\mathbb{E}}$ the expectation of the process with generator
$\tilde{\mathcal{L}}$, we have, for $\lambda\in{\mathbb{R}}$,

$(\lambda\boldsymbol{I}-\mathcal{L})^{-1}f(x,i)=\left(\frac{\lambda\boldsymbol{I}}{\psi}-\tilde{\mathcal{L}}\right)^{-1}\left(\frac{f(x,i)}{\psi(i)}\right)\\\
=\int_{0}^{\infty}{\rm d}t\,\,\tilde{\mathbb{E}}_{(x,i)}\left[{\rm
e}^{-\int_{0}^{t}{\rm
d}s\,\frac{\lambda}{\psi(i_{s})}}\,\,\frac{f(X_{t},i_{t})}{\psi(i_{t})}\right],$

and by the decoupling of $X_{t}$ and $i_{t}$ under $\tilde{\mathcal{L}}$, we
get

$\displaystyle(\lambda\boldsymbol{I}-\mathcal{L})^{-1}f(x,i)$
$\displaystyle=\int_{0}^{\infty}{\rm d}t\,\,\tilde{\mathbb{E}}_{i}\left[{\rm
e}^{-\lambda\int_{0}^{t}{\rm
d}s\,\frac{1}{\psi(i_{s})}}\,\,\frac{\boldsymbol{1}_{\\{0\\}}(i_{t})}{\psi(i_{t})}\right](S_{t}f(\cdot,0))(x)+\int_{0}^{\infty}{\rm
d}t\,\,\tilde{\mathbb{E}}_{i}\left[{\rm
e}^{-\lambda\int_{0}^{t}\frac{1}{\psi(i_{s})}}\,\,\frac{\boldsymbol{1}_{\\{1\\}}(i_{t})}{\psi(i_{t})}\right](S_{t}f(\cdot,1))(x)$
$\displaystyle=\int_{0}^{\infty}{\rm d}t\,\,\tilde{\mathbb{E}}_{i}\left[{\rm
e}^{-\lambda\int_{0}^{t}{\rm
d}s\,\frac{1}{\psi(i_{s})}}\,\,\boldsymbol{1}_{\\{0\\}}(i_{t})\right](S_{t}f(\cdot,0))(x)+\frac{1}{\epsilon}\int_{0}^{\infty}{\rm
d}t\,\,\tilde{\mathbb{E}}_{i}\left[{\rm e}^{-\lambda\int_{0}^{t}{\rm
d}s\,\frac{1}{\psi(i_{s})}}\,\,\boldsymbol{1}_{\\{1\\}}(i_{t})\right](S_{t}f(\cdot,1))(x).$

Defining

$A:=\left[\begin{array}[]{c c c c}-\Upsilon&\Upsilon\\\ \Upsilon&-\Upsilon\\\
\end{array}\right],\qquad\psi_{\epsilon}:=\left[\begin{array}[]{c c c c}1&0\\\
0&\epsilon\\\ \end{array}\right],$

and using again the Feynman-Kac formula, we have

$(\lambda\boldsymbol{I}-\mathcal{L})^{-1}\left[\begin{array}[]{c c c
c}f(x,0)\\\ f(x,1)\\\ \end{array}\right]=\int_{0}^{\infty}{\rm
d}t\,K_{\epsilon}(t,\lambda)\left[\begin{array}[]{c c c
c}(S_{t}f(\cdot,0))(x)\\\ (S_{t}f(\cdot,1))(x)\\\ \end{array}\right]$

with $K_{\epsilon}(t,\lambda)={\rm
e}^{t\psi_{\epsilon}^{-1}(-\lambda\boldsymbol{I}+A)}\psi_{\epsilon}^{-1}$.

Using the explicit formula for the exponential of a $2\times 2$ matrix (see
e.g. [4, Corollary 2.4]), we obtain

(2.14) ${\rm e}^{t\psi_{\epsilon}^{-1}(-\lambda\boldsymbol{I}+A)}={\rm
e}^{-\tfrac{1+\epsilon}{\epsilon}\ell(\Upsilon,\lambda)t}\left[\begin{array}[]{c
c c
c}\cosh(tc_{\epsilon}(\Upsilon,\lambda))+\tfrac{1-\epsilon}{\epsilon}\ell(\Upsilon,\lambda)\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}&\Upsilon\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}\\\
\tfrac{\Upsilon}{\epsilon}\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}&\cosh(tc_{\epsilon}(\Upsilon,\lambda))-\tfrac{1-\epsilon}{\epsilon}\ell(\Upsilon,\lambda)\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}\\\
\end{array}\right]$

with
$c_{\epsilon}(\Upsilon,\lambda)=\sqrt{\left(\tfrac{1-\epsilon}{\epsilon}\right)^{2}\ell(\Upsilon,\lambda)^{2}+\frac{\Upsilon^{2}}{\epsilon}}$
and $\ell(\Upsilon,\lambda)=\frac{\Upsilon+\lambda}{2}$, from which we obtain
(LABEL:eq:_resolvent_epsilon_>0).

Case $\epsilon=0$. We derive $K_{0}(t,\lambda)$ by taking the limit
$\epsilon\downarrow 0$ in the previous expression, i.e.,
$K_{0}(t,\lambda)=\lim_{\epsilon\downarrow 0}K_{\epsilon}(t,\lambda)$. We thus
have that $K_{0}(t,\lambda)$ is equal to

$\displaystyle\lim_{\epsilon\downarrow 0}{\rm
e}^{-\tfrac{1+\epsilon}{\epsilon}\ell(\Upsilon,\lambda)t}\left[\begin{array}[]{c
c c
c}\cosh(tc_{\epsilon}(\Upsilon,\lambda))+\tfrac{1-\epsilon}{\epsilon}\ell(\Upsilon,\lambda)\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}&\tfrac{\Upsilon}{\epsilon}\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}\\\
\tfrac{\Upsilon}{\epsilon}\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}&\tfrac{1}{\epsilon}\cosh(tc_{\epsilon}(\Upsilon,\lambda))-\tfrac{1-\epsilon}{\epsilon^{2}}\ell(\Upsilon,\lambda)\tfrac{\sinh(tc_{\epsilon}(\Upsilon,\lambda))}{c_{\epsilon}(\Upsilon,\lambda)}\\\
\end{array}\right]$ $\displaystyle={\rm
e}^{-\lambda\tfrac{2\Upsilon+\lambda}{\Upsilon+\lambda}t}\left[\begin{array}[]{c
c c c}1&\tfrac{\Upsilon}{\Upsilon+\lambda}\\\
\tfrac{\Upsilon}{\Upsilon+\lambda}&\left(\tfrac{\Upsilon}{\Upsilon+\lambda}\right)^{2}\\\
\end{array}\right],$

from which (2.13) follows. ∎

######  Remark 2.6.

[Symmetric layers] Note that for $\epsilon=1$ we have

$\displaystyle(\lambda\boldsymbol{I}-\mathcal{L})^{-1}f(x,i)=\int_{0}^{\infty}{\rm
d}t\,{\rm e}^{-\lambda t}\left(\tfrac{1+{\rm e}^{-2\Upsilon
t}}{2}\,(S_{t}f(\cdot,i))(x)+\tfrac{1-{\rm e}^{-2\Upsilon
t}}{2}\,(S_{t}f(\cdot,1-i))(x)\right).$

$\spadesuit$

We conclude this section by noting that the system in (2.2) was studied in
detail in [2, 33]. By taking Fourier and Laplace transforms and inverting
them, it is possible to deduce explicitly the solution, which is expressed in
terms of solutions to the classical heat equation. More precisely, using
formula [33, Eq.2.2], we have that

(2.15) $\rho_{0}(x,t)={\rm e}^{-\Upsilon
t}\,(S_{t}\,\bar{\rho}_{0})(x)+\frac{\Upsilon}{1-\epsilon}{\rm e}^{-\Upsilon
t}\int_{\epsilon t}^{t}{\rm d}s\,\left(\left(\frac{s-\epsilon
t}{t-s}\right)^{1/2}I_{1}(\upsilon(s))\,(S_{s}\,\bar{\rho}_{0})(x)+I_{0}(\upsilon(s))\,(S_{s}\,\bar{\rho}_{1})(x)\right)$

and

(2.16) $\allowdisplaybreaks\rho_{1}(x,t)={\rm e}^{-\Upsilon t}(S_{\epsilon
t}\,\bar{\rho}_{1})(x)\\\ +\frac{\Upsilon}{1-\epsilon}{\rm e}^{-\sigma
t}\int_{\epsilon t}^{t}{\rm d}s\,\left(\left(\frac{s-\epsilon
t}{t-s}\right)^{-1/2}I_{1}(\upsilon(s))\,(S_{s}\,\bar{\rho}_{1})(x)+I_{0}(\upsilon(s))\,(S_{s}\,\bar{\rho}_{0})(x)\right),$

where $\upsilon(s)=\frac{2\Upsilon}{1-\epsilon}((t-s)(s-\epsilon t))^{1/2}$,
and $I_{0}(\cdot)$ and $I_{1}(\cdot)$ are the modified Bessel functions.

##  3 The system with boundary reservoirs

In this section we consider a finite version of the switching interacting
particle systems introduced in Definition 1.1 to which boundary reservoirs are
added. Section 3.1 defines the model. Section 3.2 identifies the dual and the
stationary measures. Section 3.3 derives the non-equilibrium density profile,
both for the microscopic system and the macroscopic system, and offers various
simulations. In Section 3.4 we compute the stationary horizontal current of
slow and fast particles both for the microscopic system and the macroscopic
system. Section 3.5 shows that in the macroscopic system, for certain choices
of the rates, there can be a flow of particles uphill, i.e., against the
gradient imposed by the reservoirs. Thus, as a consequence of the competing
driving mechanisms of slow and fast particles, we can have a flow of particles
from the side with lower density to the side with higher density.

###  3.1 Model

We consider the same system as in Definition 1.1, but restricted to
$V:=\\{1,\ldots,N\\}\subset{\mathbb{Z}}$. In addition, we set
$\hat{V}:=V\cup\\{L,R\\}$ and attach a left-reservoir to $L$ and a right-
reservoir to $R$, both for fast and slow particles. To be more precise, there
are four reservoirs (see Fig. 2):

(a) Representation via _slow_ and _fast_ particles moving on $V$.

(b) Representation via particles moving on $V\times I$.

Figure 2: Case $\sigma=0$, $\epsilon>0$ with boundary reservoirs: two
equivalent representations.

  * i)

For the fast particles, a left-reservoir at $L$ injects fast particles at
$x=1$ at rate $\rho_{L,0}(1+\sigma\eta_{0}(1,t))$ and a right-reservoir at $R$
injects fast particles at $x=N$ at rate $\rho_{R,0}(1+\sigma\eta_{0}(N,t))$.
The left-reservoir absorbs fast particles at rate $1+\sigma\rho_{L,0}$, while
the right-reservoir does so at rate $1+\sigma\rho_{R,0}$.

  * ii)

For the slow particles, a left-reservoir at $L$ injects slow particles at
$x=1$ at rate $\rho_{L,1}(1+\sigma\eta_{1}(1,t))$ and a right-reservoir at $R$
injects slow particles at $x=N$ at rate $\rho_{R,1}(1+\sigma\eta_{1}(N,t))$.
The left-reservoir absorbs fast particles at rate $1+\sigma\rho_{L,1}$, while
the right-reservoir does so at rate $1+\sigma\rho_{R,1}$.

Inside $V$, the particles move as before.

For $i\in I$, $x\in V$ and $t\geq 0$, let $\eta_{i}(x,t)$ denote the number of
particles in layer $i$ at site $x$ at time $t$. For $\sigma\in\\{-1,0,1\\}$,
the Markov process $\\{\eta(t)\colon\,t\geq 0\\}$ with

$\eta(t)=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V}$

has state space

$\mathcal{X}=\begin{cases}\\{0,1\\}^{V}\times\\{0,1\\}^{V},&\sigma=-1,\\\
{\mathbb{N}}_{0}^{V}\times{\mathbb{N}}_{0}^{V},&\sigma=0,1,\end{cases}$

and generator

(3.1) $L:=L_{\epsilon,\gamma,N}=L^{\mathrm{bulk}}+L^{\mathrm{res}}$

with

(3.2) $\displaystyle L^{\mathrm{bulk}}:=L^{\mathrm{bulk}}_{0}+\epsilon
L^{\mathrm{bulk}}_{1}+\gamma L^{\mathrm{bulk}}_{0\updownarrow 1}$

acting on bounded cylindrical functions $f\colon\,\mathcal{X}\to{\mathbb{R}}$
as

$\displaystyle(L^{\mathrm{bulk}}_{0}f)(\eta)=$
$\displaystyle\sum_{x=1}^{N-1}\Big{\\{}\eta_{0}(x)(1+\sigma\eta_{0}(x+1))\,\big{[}f(\eta_{0}-\delta_{x}+\delta_{x+1},\eta_{1})-f(\eta_{0},\eta_{1})\big{]}$
$\displaystyle+\eta_{0}(x+1)(1+\sigma\eta_{0}(x))\,\big{[}f(\eta_{0}-\delta_{x+1}+\delta_{x},\eta)-f(\eta_{0},\eta_{1})\big{]}\Big{\\}},$
$\displaystyle(L^{\mathrm{bulk}}_{1}f)(\eta)=$
$\displaystyle\sum_{x=1}^{N-1}\Big{\\{}\eta_{1}(x)(1+\sigma\eta_{1}(x+1))\,\big{[}f(\eta_{0},\eta_{1}-\delta_{x}+\delta_{x+1})-f(\eta_{0},\eta_{1})\big{]}$
$\displaystyle+\eta_{1}(x+1)(1+\sigma\eta_{1}(x))\,\big{[}f(\eta_{0},\eta_{1}-\delta_{x+1}+\delta_{x})-f(\eta_{0},\eta_{1}))\big{]}\Big{\\}},$
$\displaystyle(L^{\mathrm{bulk}}_{0\updownarrow 1}f)(\eta)=$
$\displaystyle\sum_{x=1}^{N}\Big{\\{}\eta_{0}(x)(1+\sigma\eta_{1}(x))\,\big{[}f(\eta_{0}-\delta_{x},\eta_{1}+\delta_{x})-f(\eta_{0},\eta_{1})\big{]}$
$\displaystyle+\eta_{1}(x)(1+\sigma\eta_{0}(x))\,\big{[}f(\eta_{0}+\delta_{x},\eta_{1}-\delta_{x})-f(\eta_{0},\eta_{1}))\big{]}\Big{\\}},$

and

(3.3) $\displaystyle
L^{\mathrm{res}}:=L^{\mathrm{res}}_{0}+L^{\mathrm{res}}_{1}$

acting as

$\displaystyle(L^{\mathrm{res}}_{0}f)(\eta)=\eta_{0}(1)(1+\sigma\rho_{L,0})\,\big{[}f(\eta_{0}-\delta_{1},\eta_{1})-f(\eta_{0},\eta_{1})\big{]}$
$\displaystyle\quad+\rho_{L,0}(1+\sigma\eta_{0}(1))\,\big{[}f(\eta_{0}+\delta_{1},\eta_{1})-f(\eta_{0},\eta_{1})\big{]}$
$\displaystyle+\eta_{0}(N)(1+\sigma\rho_{R,0})\,\big{[}f(\eta_{0}-\delta_{N},\eta_{1})-f(\eta_{0},\eta_{1})\big{]}+\rho_{R,0}(1+\sigma\eta_{0}(N))\,\big{[}f(\eta_{0}+\delta_{N},\eta)-f(\eta_{0},\eta_{1})\big{]},$
$\displaystyle(L^{\mathrm{res}}_{1}f)(\eta)=\eta_{1}(1)(1+\sigma\rho_{L,1})\,\big{[}f(\eta_{0},\eta_{1}-\delta_{1})-f(\eta_{0},\eta_{1})\big{]}$
$\displaystyle\quad+\rho_{L,1}(1+\sigma\eta_{1}(1))\,\big{[}f(\eta_{0},\eta_{1}+\delta_{1})-f(\eta_{0},\eta_{1})\big{]}$
$\displaystyle+\eta_{1}(N)(1+\sigma\rho_{R,1})\,\big{[}f(\eta_{0},\eta_{1}-\delta_{N})-f(\eta_{0},\eta_{1})\big{]}+\rho_{R,1}(1+\sigma\rho_{R,N})\,\big{[}f(\eta_{0},\eta_{1}+\delta_{N})-f(\eta_{0},\eta_{1})\big{]}.$

###  3.2 Duality and stationary measures

In [9] it was shown that the partial exclusion process, a system of
independent random walks and the symmetric inclusion processes on a finite set
$V$, coupled with proper left and right reservoirs, are dual to the same
particle system but with the reservoirs replaced by absorbing sites. As
remarked in [27], the same result holds for more general geometries,
consisting of inhomogeneous rates (site and edge dependent), and for many
proper reservoirs. Our model is a particular instance of the case treated in
[27, Remark 2.2]), because we can think of the rate as conductances attached
to the edges.

More precisely, we consider the system where particles jump on two copies of

$\hat{V}:=V\cup\\{L,R\\}$

and follow the same dynamics as before in $V$, but with the reservoirs at $L$
and $R$ absorbing. We denote by $\xi$ the configuration

$\xi=(\xi_{0},\xi_{1}):=(\\{\xi_{0}(x)\\}_{x\in\hat{V}},\\{\xi_{1}(x)\\}_{x\in\hat{V}}),$

where $\xi_{i}(x)$ denotes the number of particles at site $x$ in layer $i$.
The state space is
$\hat{\mathcal{X}}={\mathbb{N}}_{0}^{\hat{V}}\times{\mathbb{N}}_{0}^{\hat{V}}$,
and the generator is

(3.4)
$\hat{L}:=\hat{L}_{\epsilon,\gamma,N}=\hat{L}^{\mathrm{bulk}}+\hat{L}^{L,R}$

with

$\displaystyle\hat{L}^{\mathrm{bulk}}:=\hat{L}_{0}^{\mathrm{bulk}}+\epsilon\hat{L}_{1}^{\mathrm{bulk}}+\gamma\hat{L}^{\mathrm{bulk}}_{0\updownarrow
1}$

acting on cylindrical functions $f\colon\,\mathcal{X}\to{\mathbb{R}}$ as

$\displaystyle(\hat{L}_{0}^{\mathrm{bulk}}f)(\xi)=\sum_{x=1}^{N-1}\Big{\\{}\xi_{0}(x)(1+\sigma\xi_{0}(x+1))\,\big{[}f(\xi_{0}-\delta_{x}+\delta_{x+1},\xi_{1})-f(\xi_{0},\xi_{1})\big{]}$
$\displaystyle\quad+\xi_{0}(x+1)(1+\sigma\xi_{0}(x))\,\big{[}f(\xi_{0}-\delta_{x+1}+\delta_{x},\xi_{1})-f(\xi_{0},\xi_{1})\big{]}\Big{\\}},$
$\displaystyle(\hat{L}_{1}^{\mathrm{bulk}}f)(\xi)=\sum_{x=1}^{N-1}\Big{\\{}\xi_{1}(x)(1+\sigma\xi_{1}(x+1))\,\big{[}f(\xi_{0},\xi_{1}-\delta_{x}+\delta_{x+1})-f(\xi_{0},\xi_{1})\big{]}$
$\displaystyle\quad+\xi_{1}(x+1)(1+\sigma\xi_{1}(x))\,\big{[}f(\xi_{0},\xi_{1}-\delta_{x+1}+\delta_{x})-f(\xi_{0},\xi_{1})\big{]}\Big{\\}},$
$\displaystyle(\hat{L}^{\mathrm{bulk}}_{0\updownarrow
1}f)(\eta)=\sum_{x=1}^{N}\Big{\\{}\xi_{0}(x)(1+\sigma\xi_{1}(x))\,\big{[}f(\xi_{0}-\delta_{x},\xi_{1}+\delta_{x})-f(\xi_{0},\xi_{1})\big{]}$
$\displaystyle\quad+\xi_{1}(x)(1+\sigma\xi_{0}(x))\,\big{[}f(\xi_{0}+\delta_{x},\xi_{1}-\delta_{x})-f(\xi_{0},\xi_{1})\big{]}\Big{\\}},$

and

$\displaystyle\hat{L}^{L,R}=\hat{L}_{0}^{L,R}+\hat{L}_{1}^{L,R}$

acting as

$\displaystyle(\hat{L}_{0}^{L,R}f)(\xi)=\xi_{0}(1)\,\big{[}f(\xi_{0}-\delta_{1},\xi_{1})-f(\xi_{0},\xi_{1})\big{]}+\xi_{0}(N)\,\big{[}f(\xi_{0}-\delta_{N},\xi_{1})-f(\xi_{0},\xi_{1})\big{]},$
$\displaystyle(\hat{L}_{1}^{L,R}f)(\xi)=\xi_{1}(1)\,\big{[}f(\xi_{0},\xi_{1}-\delta_{1})-f(\xi_{0},\xi_{1})\big{]}+\xi_{1}(N)\,\big{[}f(\xi_{0},\xi_{1}-\delta_{N})-f(\xi_{0},\xi_{1})\big{]}.$

######  Proposition 3.1.

[Duality] [9, Theorem 4.1] and [27, Proposition 2.3] The Markov processes

$\displaystyle\\{\eta(t)\colon\,t\geq
0\\},\qquad\eta(t)=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V},$
$\displaystyle\\{\xi(t)\colon\,t\geq
0\\},\qquad\xi(t)=\\{\xi_{0}(x,t),\xi_{1}(x,t)\\}_{x\in\hat{V}},$

with generators $L$ in (3.1) and $\hat{L}$ in (3.4) are dual. Namely, for all
configurations $\eta\in\mathcal{X}$, $\xi\in\hat{\mathcal{X}}$ and times
$t\geq 0$,

${\mathbb{E}}_{\eta}[D(\xi,\eta_{t})]={\mathbb{E}}_{\xi}[D(\xi_{t},\eta)],$

where the duality function is given by

$D(\xi,\eta):=\left(\prod_{i\in
I}d_{(L,i)}(\xi_{i}(L))\right)\times\left(\prod_{x\in
V}d(\xi_{i}(x),\eta_{i}(x))\right)\times\left(\prod_{i\in
I}d_{(R,i)}(\xi_{i}(R))\right),$

where, for $k,n\in{\mathbb{N}}$ and $i\in I$, $d(\cdot,\cdot)$ is given in
(1.6) and

$d_{(L,i)}(k)=\left(\rho_{L,i}\right)^{k},\qquad
d_{(R,i)}(k)=\left(\rho_{R,i}\right)^{k}.$

The proof boils down to checking that the relation

$\hat{L}D(\cdot,\eta)(\xi)=LD(\xi,\cdot)(\eta)$

holds for any $\xi\in\mathcal{X}$ and $\xi\in\hat{\mathcal{X}}$, as follows
from a rewriting of the proof of [9, Theorem 4.1].

###  3.3 Non-equilibrium stationary profile

Also the existence and uniqueness of the non-equilibrium steady state has been
established in [27, Theorem 3.3] for general geometries, and the argument in
that paper can be easily adapted to our setting.

######  Theorem 3.2.

[Stationary measure] [27, Theorem 3.3(a)] For $\sigma\in\\{-1,0,1\\}$ there
exists a unique stationary measure $\mu_{stat}\ $ for $\\{\eta(t)\colon\,t\geq
0\\}$. Moreover, for $\sigma=0$ and for any values of $\\{\rho_{L,0},\
\rho_{L,1},\ \rho_{R,0},\ \rho_{R,1}\\}$,

(3.5) $\mu_{stat}=\prod_{(x,i)\in V\times
I}\nu_{(x,i)},\qquad\nu_{(x,i)}=\mathrm{Poisson}(\theta_{(x,i)}),$

while, for $\sigma\in\\{-1,1\\}$, $\mu_{stat}$ is in general not in product
form, unless $\rho_{L,0}=\rho_{L,1}=\rho_{R,0}=\rho_{R,1}$, for which

(3.6) $\mu_{stat}=\prod_{(x,i)\in V\times I}\nu_{(x,i),\theta},$

where $\nu_{(x,i),\theta}$ is given in (1.4).

###### Proof.

For $\sigma=-1$, the existence and uniqueness of the stationary measure is
trivial by the irreducibility and the finiteness of the state space of the
process. For $\sigma\in\\{0,1\\}$, recall from [27, Appendix A] that a
probability measure $\mu$ on $\mathcal{X}$ is said to be tempered if it is
characterized by the integrals $\big{\\{}{\mathbb{E}}_{\mu}[D(\xi,\eta)]\,:\,\
\xi\in\hat{\mathcal{X}}\big{\\}}$ and that if there exists a
$\theta\in[0,\infty)$ such that
${\mathbb{E}}_{\mu}[D(\xi,\eta)]\leq\theta^{|\xi|}$ for any
$\xi\in\hat{\mathcal{X}}$. By means of duality we have that, for any
$\eta\in\mathcal{X}$ and $\xi\in\hat{\mathcal{X}}$,

$\displaystyle\lim_{t\to\infty}{\mathbb{E}}_{\eta}[D(\xi,\eta_{t})]=\lim_{t\to\infty}\hat{{\mathbb{E}}}_{\xi}[D(\xi_{t},\eta)]$
(3.7)
$\displaystyle=\sum_{i_{0}=0}^{|\xi|}\sum_{i_{0,L}=0}^{i_{0}}\sum_{j_{1,L}=0}^{|\xi|-i_{0}}\rho_{L,0}^{i_{0,L}}\
\rho_{R,0}^{i_{0}-i_{0,L}}\ \rho_{L,1}^{i_{1,L}}\
\rho_{R,1}^{|\xi|-i_{0}-i_{1,L}}$ (3.8)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\times\hat{{\mathbb{P}}}_{\xi}\left(\xi_{\infty}=i_{0,L}\delta_{(L,0)}+(i_{0}-i_{0,L})\delta_{(R,0)}+i_{1,L}\delta_{(L,1)}+(|\xi|-i_{0}-i_{1,L})\delta_{(R,1)}\right),$

from which we conclude that
$\lim_{t\to\infty}{\mathbb{E}}_{\eta}[D(\xi,\eta_{t})]\leq\max\\{\rho_{L,0},\rho_{R,0},\rho_{L,1},\rho_{R,1}\\}^{|\xi|}$.
Let $\mu_{s}$ be the unique tempered probability measure such that for any
$\xi\in\hat{\mathcal{X}},$ ${\mathbb{E}}_{\mu_{stat}}[D(\xi,\eta)]$ coincides
with (3.3). From the convergence of the marginal moments in (3.3) we conclude
that, for any $f:\mathcal{X}\to{\mathbb{R}}$ bounded and for any
$\eta\in\mathcal{X}$,

$\lim_{t\to\infty}{\mathbb{E}}_{\eta}[f(\eta_{t})]={\mathbb{E}}_{\mu_{stat}}[f(\eta)].$

Thus, a dominated convergence argument yields that for any probability measure
$\mu$ on $\mathcal{X}$,

$\lim_{t\to\infty}{\mathbb{E}}_{\mu}[f(\eta_{t})]={\mathbb{E}}_{\mu_{stat}}[f(\eta)],$

giving that $\mu_{stat}$ is the unique stationary measure. The explicit
expression in (3.5) and (3.6) follows from similar computations as in [9],
while, arguing by contradiction as in the proof of [27, Theorem 3.3], we can
show that the two-point truncated correlations are non-zero for
$\sigma\in\\{-1,1\\}$ whenever at least two reservoir parameters are
different. ∎

####  3.3.1 Stationary microscopic profile and absorption probability

In this section we provide an explicit expression for the stationary
microscopic density of each type of particle. To this end, let $\mu_{stat}$ be
the unique non-equilibrium stationary measure of the process

$\\{\eta(t)\colon\,t\geq
0\\},\qquad\eta(t):=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V},$

and let $\\{\theta_{0}({x}),\theta_{1}({x})\\}_{x\in V}$ be the stationary
microscopic profile, i.e., for $x\in V$ and $i\in I$,

(3.9) $\theta_{i}({x})={\mathbb{E}}_{\mu_{stat}}[\eta_{i}(x,t)].$

Write ${\mathbb{P}}_{\xi}$ (and ${\mathbb{E}}_{\xi}$) to denote the law (and
the expectation) of the dual Markov process

$\\{\xi(t)\colon\,t\geq
0\\},\qquad\xi(t):=\\{\xi_{0}(x,t),\xi_{1}(x,t)\\}_{x\in\hat{V}},$

starting from $\xi=\\{\xi_{0}(x),\xi_{1}(x)\\}_{x\in\hat{V}}$. For $x\in V$,
set

(3.10) $\displaystyle\vec{p}_{x}$ $\displaystyle:=\big{[}\begin{array}[]{c c c
c}\hat{p}(\delta_{(x,0)},\delta_{(L,0)})&\hat{p}(\delta_{(x,0)},\delta_{(L,1)})&\hat{p}(\delta_{(x,0)},\delta_{(R,0)})&\hat{p}(\delta_{(x,0)},\delta_{(R,1)})\end{array}\big{]}^{T},$
$\displaystyle\vec{q}_{x}$ $\displaystyle:=\big{[}\begin{array}[]{c c c
c}\hat{p}(\delta_{(x,1)},\delta_{(L,0)})&\hat{p}(\delta_{(x,1)},\delta_{(L,1)})&\hat{p}(\delta_{(x,1)},\delta_{(R,0)})&\hat{p}(\delta_{(x,1)},\delta_{(R,1)})\end{array}\big{]}^{T},$

where

(3.11)
$\hat{p}(\xi,\tilde{\xi})=\lim\limits_{t\to\infty}{\mathbb{P}}_{\xi}(\xi(t)=\tilde{\xi}),\quad\xi=\delta_{(x,i)}\text{
for some }(x,i)\in V\times I,\
\tilde{\xi}\in\\{\delta_{(L,0)},\delta_{(L,1)},\delta_{(R,0)},\delta_{(R,1)}\\},$

and let

(3.12) $\vec{\rho}:=\big{[}\begin{array}[]{c c c
c}\rho_{(L,0)}&\rho_{(L,1)}&\rho_{(R,0)}&\rho_{(R,1)}\end{array}\big{]}^{T}.$

Note that $\hat{p}(\delta_{(x,i)},\cdot)$ is the probability of the dual
process, starting from a single particle at site $x$ at layer $i\in I$, of
being absorbed at one of the four reservoirs. Using Proposition 3.1 and
Theorem 3.2, we obtain the following.

######  Corollary 3.3.

[Dual representation of stationary profile] For $x\in V$, the microscopic
stationary profile is given by

(3.13) $\begin{aligned} \theta_{0}(x)&=\vec{p}_{x}\,\cdot\,\vec{\rho},\\\
\theta_{1}(x)&=\vec{q}_{x}\,\cdot\,\vec{\rho},\end{aligned}\qquad
x\in\\{1,\ldots,N\\},$

where $\vec{p}_{x},\vec{q}_{x}$ and $\vec{\rho}$ are as in (3.10)–(3.12).

We next compute the absorption probabilities associated to the dual process in
order to obtain a more explicit expression for the stationary microscopic
profile $\\{\theta_{0}({x}),\theta_{1}({x})\\}_{x\in V}$. The absorption
probabilities $\hat{p}(\cdot\,,\cdot)$ of the dual process satisfy

$(\hat{L}\hat{p})(\cdot,\tilde{\xi})(\xi)=0\qquad\forall\,\xi\in\hat{\mathcal{X}},$

where $\hat{L}$ is the dual generator defined in (3.4), i.e., they are
harmonic functions for the generator $\hat{L}$.

In matrix form, the above translates into the following systems of equations:

(3.14) $\displaystyle\vec{p}_{1}$
$\displaystyle=\frac{1}{2+\gamma}\,(\vec{p}_{0}+\vec{p}_{2})+\frac{\gamma}{2+\gamma}\,\vec{q}_{1},$
$\displaystyle\vec{q}_{1}$
$\displaystyle=\frac{\epsilon}{(1+\epsilon)+\gamma}\,\vec{q}_{2}+\frac{1}{(1+\epsilon)+\gamma}\,\vec{q}_{0}+\frac{\gamma}{(1+\epsilon)+\gamma}\,\vec{p}_{1},$
$\displaystyle\vec{p}_{x}$
$\displaystyle=\frac{1}{2+\gamma}\,(\vec{p}_{x-1}+\vec{p}_{x+1})+\frac{\gamma}{2+\gamma}\,\vec{q}_{x},$
$\displaystyle x\in\\{2,\ldots,N-1\\},$ $\displaystyle\vec{q}_{x}$
$\displaystyle=\frac{\epsilon}{2\epsilon+\gamma}\,(\vec{q}_{x-1}+\vec{q}_{x+1})+\frac{\gamma}{2\epsilon+\gamma}\,\vec{p}_{x},$
$\displaystyle x\in\\{2,\ldots,N-1\\},$ $\displaystyle\vec{p}_{N}$
$\displaystyle=\frac{1}{2+\gamma}\,(\vec{p}_{N-1}+\vec{p}_{N+1})+\frac{\gamma}{2+\gamma}\,\vec{q}_{N},$
$\displaystyle\vec{q}_{N}$
$\displaystyle=\frac{\epsilon}{(1+\epsilon)+\gamma}\,\vec{q}_{N-1}+\frac{1}{(1+\epsilon)+\gamma}\,\vec{q}_{N+1}+\frac{\gamma}{(1+\epsilon)+\gamma}\,\vec{p}_{N},$

where

$\displaystyle\vec{p}_{0}:=\big{[}\begin{array}[]{c c c
c}1&0&0&0\end{array}\big{]}^{T},$
$\displaystyle\vec{q}_{0}:=\big{[}\begin{array}[]{c c c
c}0&1&0&0\end{array}\big{]}^{T},$
$\displaystyle\vec{p}_{N+1}:=\big{[}\begin{array}[]{c c c
c}0&0&1&0\end{array}\big{]}^{T},$
$\displaystyle\vec{q}_{N+1}:=\big{[}\begin{array}[]{c c c
c}0&0&0&1\end{array}\big{]}^{T}.$

We divide the analysis of the absorption probabilities into two cases:
$\epsilon=0$ and $\epsilon>0$.

##### Case $\epsilon=0$.

######  Proposition 3.4.

[Absorption probability for $\epsilon=0$] Consider the dual process

$\\{\xi(t)\colon\,t\geq
0\\},\qquad\xi(t)=\\{\xi_{0}(x,t),\xi_{1}(x,t)\\}_{x\in V},$

with generator $\hat{L}_{\epsilon,\gamma,N}$ (see (3.4)) with $\epsilon=0$.
Then for the dual process, starting from a single particle, the absorption
probabilities $\hat{p}(\cdot,\cdot)$ (see (3.11)) are given by

(3.15) $\displaystyle\hat{p}(\delta_{(x,0)},\delta_{(L,0)})$
$\displaystyle=\frac{1+\gamma}{1+2\gamma}\left(\frac{(1+N)+(1+2N)\,\gamma}{1+N+2N\gamma}-\frac{1+2\gamma}{1+N+2N\gamma}\,x\right),$
$\displaystyle\hat{p}(\delta_{(x,0)},\delta_{(L,1)})$
$\displaystyle=\frac{\gamma}{1+2\gamma}\left(\frac{(1+N)+(1+2N)\,\gamma}{1+N+2N\gamma}-\frac{1+2\gamma}{1+N+2N\gamma}\,x\right),$
$\displaystyle\hat{p}(\delta_{(x,0)},\delta_{(R,0)})$
$\displaystyle=\frac{1+\gamma}{1+2\gamma}\left(\frac{-\gamma}{1+N+2N\gamma}+\frac{1+2\gamma}{1+N+2N\gamma}\,x\right),$
$\displaystyle\hat{p}(\delta_{(x,0)},\delta_{(R,1)})$
$\displaystyle=\frac{\gamma}{1+2\gamma}\left(\frac{-\gamma}{1+N+2N\gamma}+\frac{1+2\gamma}{1+N+2N\gamma}\,x\right),$
(3.16) $\displaystyle\hat{p}(\delta_{(1,1)},\delta_{(L,0)})$
$\displaystyle=\frac{\gamma\,(N-\gamma+2N\gamma)}{(1+2\gamma)(1+N+2N\gamma)},\quad\hat{p}(\delta_{(1,1)},\delta_{(L,1)})=\frac{1+N+(1+3N)\gamma-(1-2N)\gamma^{2}}{(1+2\gamma)(1+N+2N\gamma)},$
$\displaystyle\hat{p}(\delta_{(1,1)},\delta_{(R,0)})$
$\displaystyle=\frac{\gamma(1+\gamma)}{(1+2\gamma)(1+N+2N\gamma)},\quad\hat{p}(\delta_{(1,1)},\delta_{(R,1)})=\frac{\gamma^{2}}{(1+2\gamma)(1+N+2N\gamma)},$

and

(3.17)
$\hat{p}(\delta_{(x,1)},\delta_{(\beta,i)})=\hat{p}(\delta_{(x,0)},\delta_{(\beta,i)}),\qquad
x\in\\{2,\ldots,N-1\\},\,(\beta,i)\in\\{L,R\\}\times I,$

and

(3.18) $\displaystyle\hat{p}(\delta_{(N,1)},\delta_{(L,0)})$
$\displaystyle=\hat{p}(\delta_{(1,1)},\delta_{(R,0)}),\quad\hat{p}(\delta_{(N,1)},\delta_{(L,1)})=\hat{p}(\delta_{(1,1)},\delta_{(R,1)}),$
$\displaystyle\hat{p}(\delta_{(N,1)},\delta_{(R,0)})$
$\displaystyle=\hat{p}(\delta_{(1,1)},\delta_{(L,0)}),\quad\hat{p}(\delta_{(N,1)},\delta_{(R,1)})=\hat{p}(\delta_{(1,1)},\delta_{(L,1)}).$

###### Proof.

Note that, for $\epsilon=0$, from the linear system in (3.14) we get

(3.19) $\begin{aligned}
&\vec{p}_{x+1}-\vec{p}_{x}=\vec{p}_{x}-\vec{p}_{x-1},\\\
&\vec{q}_{x}=\vec{p}_{x},\end{aligned}\quad x\in\\{2,\ldots,N-1\\}.$

Thus, if we set $\vec{c}=\vec{p}_{2}-\vec{p}_{1}$, then it suffices to solve
the following 4 linear equations with 4 unknowns
$\vec{p}_{1},\vec{c},\vec{q}_{1}$, $\vec{q}_{N}$:

(3.20)
$\displaystyle\vec{p}_{1}=\frac{1}{2+\gamma}\,(\vec{p}_{0}+\vec{p}_{1}+\vec{c})+\frac{\gamma}{2+\gamma}\,\vec{q}_{1},$
$\displaystyle\vec{q}_{1}=\frac{1}{1+\gamma}\,\vec{q}_{0}+\frac{\gamma}{1+\gamma}\,\vec{p}_{1},$
$\displaystyle\vec{p}_{1}+(N-1)\vec{c}=\frac{1}{2+\gamma}\,(\vec{p}_{1}+(N-2)\vec{c}+\vec{p}_{N+1})+\frac{\gamma}{2+\gamma}\,\vec{q}_{N},$
$\displaystyle\vec{q}_{N}=\frac{1}{1+\gamma}\,\vec{q}_{N+1}+\frac{\gamma}{1+\gamma}\,(\vec{p}_{1}+(N-1)\vec{c}).$

Solving the above equations we get the desired result. ∎

As a consequence, we obtain the stationary microscopic profile for the
original process $\\{\eta(t)\colon\,t\geq 0\\},\
\eta(t)=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V}$ when $\epsilon=0$.

######  Theorem 3.5.

[Stationary microscopic profile for $\epsilon=0$]
The stationary microscopic profile $\\{\theta_{0}(x),\theta_{1}(x)\\}_{x\in
V}$ (see (3.9)) for the process $\\{\eta(t)\colon\,t\geq 0\\}$ with
$\eta(t)=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V}$ with generator
$L_{\epsilon,\gamma,N}$ (see (3.1)) and $\epsilon=0$ is given by

(3.21) $\displaystyle\theta_{0}(x)$
$\displaystyle=\frac{1+\gamma}{1+2\gamma}\left[\left(\tfrac{(1+N)+(1+2N)\,\gamma}{1+N+2N\gamma}-\tfrac{1+2\gamma}{1+N+2N\gamma}\,x\right)\rho_{L,0}+\left(\tfrac{-\gamma}{1+N+2N\gamma}+\tfrac{1+2\gamma}{1+N+2N\gamma}\,x\right)\rho_{R,0}\right]$
$\displaystyle\qquad+\frac{\gamma}{1+2\gamma}\left[\left(\tfrac{(1+N)+(1+2N)\,\gamma}{1+N+2N\gamma}-\tfrac{1+2\gamma}{1+N+2N\gamma}\,x\right)\rho_{(L,1)}+\left(\tfrac{-\gamma}{1+N+2N\gamma}+\tfrac{1+2\gamma}{1+N+2N\gamma}\,x\right)\rho_{(R,1)}\right]$

and

(3.22) $\displaystyle\theta_{1}(1)$
$\displaystyle=\frac{\gamma}{1+\gamma}\,\theta_{0}(1)+\frac{1}{1+\gamma}\,\rho_{(L,1)},$
$\displaystyle\theta_{1}(x)$ $\displaystyle=\theta_{0}(x),$ $\displaystyle
x\in\\{2,\ldots,N-1\\},$ $\displaystyle\theta_{1}(N)$
$\displaystyle=\frac{\gamma}{1+\gamma}\,\theta_{0}(N)+\frac{1}{1+\gamma}\,\rho_{(R,1)}.$

###### Proof.

The proof directly follows from Corollary 3.3 and Proposition 3.4. ∎

##### Case $\epsilon>0$.

We next compute the absorption probabilities for the dual process and the
stationary microscopic profile for the original process when $\epsilon>0$.

######  Proposition 3.6.

[Absorption probability for $\epsilon>0$] Consider the dual process

$\\{\xi(t)\colon\,t\geq
0\\},\qquad\xi(t)=\\{\xi_{0}(x,t),\xi_{1}(x,t)\\}_{x\in V},$

with generator $\hat{L}_{\epsilon,\gamma}$ (see (3.4)) with $\epsilon>0$. Let
$\hat{p}(\cdot,\cdot)$ (see (3.11)) be the absorption probabilities of the
dual process starting from a single particle, and let
$(\vec{p}_{x},\vec{q}_{x})_{x\in V}$ be as defined in (3.10). Then

(3.23) $\begin{aligned}
\vec{p}_{x}&=\vec{c}_{1}\,x+\vec{c}_{2}+\epsilon(\vec{c}_{3}\,\alpha_{1}^{x}+\vec{c}_{4}\,\alpha_{2}^{x}),\\\
\vec{q}_{x}&=\vec{c}_{1}\,x+\vec{c}_{2}-(\vec{c}_{3}\,\alpha_{1}^{x}+\vec{c}_{4}\,\alpha_{2}^{x}),\end{aligned}\quad
x\in V,$

where $\alpha_{1},\alpha_{2}$ are the two roots of the equation

(3.24) $\epsilon\alpha^{2}-(\gamma(1+\epsilon)+2\epsilon)\,\alpha+\epsilon=0,$

and $\vec{c}_{1},\vec{c}_{2},\vec{c}_{3},\vec{c}_{4}$ are vectors that depend
on the parameters $N,\epsilon,\alpha_{1},\alpha_{2}$ (see (A.4) for explicit
expressions).

###### Proof.

Applying the transformation

(3.25)
$\vec{\tau}_{x}:=\vec{p}_{x}+\epsilon\vec{q}_{x},\qquad\vec{s}_{x}:=\vec{p}_{x}-\vec{q}_{x},$

we see that the system in (3.14) decouples in the bulk (i.e., the interior of
$V$), and

(3.26)
$\vec{\tau}_{x}=\frac{1}{2}(\vec{\tau}_{x+1}+\vec{\tau}_{x-1}),\qquad\vec{s}_{x}=\frac{\epsilon}{\gamma(1+\epsilon)+2\epsilon}(\vec{s}_{x+1}+\vec{s}_{x-1}),\qquad
x\in\\{2,\ldots,N-1\\}.$

The solution of the above system of recursion equations takes the form

(3.27)
$\vec{\tau}_{x}=\vec{A}_{1}x+\vec{A}_{2},\quad\vec{s}_{x}=\vec{A}_{3}\alpha_{1}^{x}+\vec{A}_{4}\alpha_{2}^{x},$

where $\alpha_{1},\alpha_{2}$ are the two roots of the equation

(3.28) $\epsilon\alpha^{2}-(\gamma(1+\epsilon)+2\epsilon)\,\alpha+\epsilon=0.$

Rewriting the four boundary conditions in (3.14) in terms of the new
transformations, we get

(3.29) $\big{[}\begin{array}[]{c c c
c}\vec{A}_{1}&\vec{A}_{2}&\vec{A}_{3}&\vec{A}_{4}\end{array}\big{]}=(1+\epsilon)(M_{\epsilon}^{-1})^{T},$

where $M_{\epsilon}$ is given by

(3.30) $M_{\epsilon}:=\left[\begin{array}[]{c c c c}0&1&\epsilon&\epsilon\\\
1-\epsilon&1&(\epsilon-1)\alpha_{1}-\epsilon&(\epsilon-1)\alpha_{2}-\epsilon\\\
N+1&1&\epsilon\alpha_{1}^{N+1}&\epsilon\alpha_{2}^{N+1}\\\
N+\epsilon&1&-\alpha_{1}^{N}(\epsilon\alpha_{1}+(1-\epsilon))&-\alpha_{2}^{N}(\epsilon\alpha_{2}+(1-\epsilon))\end{array}\right].$

Since $\vec{p}_{x}=\frac{1}{1+\epsilon}(\vec{\tau}_{x}+\epsilon\vec{s}_{x})$
and $\vec{q}_{x}=\frac{1}{1+\epsilon}(\vec{\tau}_{x}-\vec{s}_{x})$, by setting

$\vec{c}_{i}=\frac{1}{1+\epsilon}\vec{A}_{i},\quad i\in\\{1,2,3,4\\},$

we get the desired identities. ∎

Without loss of generality, from here onwards, we fix the choices of the roots
$\alpha_{1}$ and $\alpha_{2}$ of the quadratic equation in (3.24) as

(3.31)
$\alpha_{1}=1+\frac{\gamma}{2}\left(1+\frac{1}{\epsilon}\right)-\sqrt{\left[1+\frac{\gamma}{2}\left(1+\frac{1}{\epsilon}\right)\right]^{2}-1},\quad\alpha_{2}=1+\frac{\gamma}{2}\left(1+\frac{1}{\epsilon}\right)+\sqrt{\left[1+\frac{\gamma}{2}\left(1+\frac{1}{\epsilon}\right)\right]^{2}-1}.$

Note that, for any $\epsilon,\gamma>0$, we have

(3.32) $\alpha_{1}\alpha_{2}=1.$

As a corollary, we get the expression for the stationary microscopic profile
of the original process.

######  Theorem 3.7.

[Stationary microscopic profile for $\epsilon>0$]
The stationary microscopic profile $\\{\theta_{0}(x),\theta_{1}(x)\\}_{x\in
V}$ (see (3.9)) for the process $\\{\eta(t)\colon\,t\geq 0\\}$ and
$\eta(t)=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V}$ with generator
$L_{\epsilon,\gamma,N}$ (see (3.1)) with $\epsilon>0$ is given by

(3.33) $\begin{aligned}
\theta_{0}(x)&=(\vec{c}_{1}\,.\,\vec{\rho})x+(\vec{c}_{2}\,.\,\vec{\rho})+\epsilon(\vec{c}_{3}\,.\,\vec{\rho})\alpha_{1}^{x}+\epsilon(\vec{c}_{4}\,.\,\vec{\rho})\alpha_{2}^{x},\\\
\theta_{1}(x)&=(\vec{c}_{1}\,.\,\vec{\rho})x+(\vec{c}_{2}\,.\,\vec{\rho})-(\vec{c}_{3}\,.\,\vec{\rho})\alpha_{1}^{x}-(\vec{c}_{4}\,.\,\vec{\rho})\alpha_{2}^{x},\end{aligned}\quad
x\in V,$

where $(\vec{c}_{i})_{1\leq i\leq 4}$ are as in (A.4), and

$\vec{\rho}:=\big{[}\begin{array}[]{c c c
c}\rho_{(L,0)}&\rho_{(L,1)}&\rho_{(R,0)}&\rho_{(R,1)}\end{array}\big{]}^{T}.$

###### Proof.

The proof follows directly from Corollary 3.3 and Proposition 3.6. ∎

######  Remark 3.8.

[Symmetric layers] For $\epsilon=1,$ the inverse of the matrix $M_{\epsilon}$
in the proof of Proposition 3.6 takes a simpler form. This is because for
$\epsilon=1$ the system is fully symmetric. In this case, the explicit
expression of the stationary microscopic profile is given by

(3.34) $\displaystyle\theta_{0}(x)=$
$\displaystyle\frac{1}{2}\left(\frac{N+1-x}{N+1}+\frac{\alpha_{2}^{N+1-x}-\alpha_{1}^{N+1-x}}{\alpha_{2}^{N+1}-\alpha_{1}^{N+1}}\right)\rho_{L,0}+\frac{1}{2}\left(\frac{x}{N+1}+\frac{\alpha_{2}^{x}-\alpha_{1}^{x}}{\alpha_{2}^{N+1}-\alpha_{1}^{N+1}}\right)\rho_{R,0}$
$\displaystyle+\frac{1}{2}\left(\frac{N+1-x}{N+1}-\frac{\alpha_{2}^{N+1-x}-\alpha_{1}^{N+1-x}}{\alpha_{2}^{N+1}\alpha_{1}^{N+1}}\right)\rho_{(L,1)}+\frac{1}{2}\left(\frac{x}{N+1}-\frac{\alpha_{2}^{x}-\alpha_{1}^{x}}{\alpha_{2}^{N+1}-\alpha_{1}^{N+1}}\right)\rho_{(R,1)}$

and

(3.35) $\displaystyle\theta_{1}(x)=$
$\displaystyle\frac{1}{2}\left(\frac{N+1-x}{N+1}-\frac{\alpha_{2}^{N+1-x}-\alpha_{1}^{N+1-x}}{\alpha_{2}^{N+1}-\alpha_{1}^{N+1}}\right)\rho_{L,0}+\frac{1}{2}\left(\frac{x}{N+1}-\frac{\alpha_{2}^{x}-\alpha_{1}^{x}}{\alpha_{2}^{N+1}-\alpha_{1}^{N+1}}\right)\rho_{R,0}$
$\displaystyle+\frac{1}{2}\left(\frac{N+1-x}{N+1}+\frac{\alpha_{2}^{N+1-x}-\alpha_{1}^{N+1-x}}{\alpha_{2}^{N+1}-\alpha_{1}^{N+1}}\right)\rho_{(L,1)}+\frac{1}{2}\left(\frac{x}{N+1}+\frac{\alpha_{2}^{x}-\alpha_{1}^{x}}{\alpha_{2}^{N+1}-\alpha_{1}^{N+1}}\right)\rho_{(R,1)}.$

However, note that

$\theta_{0}(x)+\theta_{1}(x)=2[(\vec{c}_{1}.\vec{\rho})x+(\vec{c}_{2}.\vec{\rho})]-(1-\epsilon)[(\vec{c}_{3}\,.\,\vec{\rho})\alpha_{1}^{x}-(\vec{c}_{4}\,.\,\vec{\rho})\alpha_{2}^{x}],$

which is linear in $x$ only when $\epsilon=1$, and

$\theta_{0}(x)-\theta_{1}(x)=(1+\epsilon)[(\vec{c}_{3}\,.\,\vec{\rho})\alpha_{1}^{x}+(\vec{c}_{4}\,.\,\vec{\rho})\alpha_{2}^{x}],$

which is purely exponential in $x$. $\spadesuit$

####  3.3.2 Stationary macroscopic profile and boundary-value problem

In this section we rescale the finite-volume system with boundary reservoirs,
in the same way as was done for the infinite-volume system in Section 2 when
we derived the hydrodynamic limit (i.e., space is scaled by $1/N$ and the
switching rate $\gamma_{N}$ is scaled such that
$\gamma_{N}N^{2}\to\Upsilon>0$), and study the validity of Fick’s law at
stationarity on macroscopic scale. Before we do that, we justify below that
the current scaling of the parameters is indeed the proper choice, in the
sense that we obtain non-trivial pointwise limits (macroscopic stationary
profiles) of the microscopic stationary profiles found in previous sections,
and that the resulting limits (when $\epsilon>0$) satisfy the stationary
boundary-value problem given in (2.2) with boundary conditions
$\rho^{\text{stat}}_{0}(0)=\rho_{L,0},\ \rho^{\text{stat}}_{0}(1)=\rho_{R,0},\
\rho^{\text{stat}}_{1}(0)=\rho_{L,1}$ and
$\rho^{\text{stat}}_{1}(1)=\rho_{R,1}$.

We say that the macroscopic stationary profiles are given by functions
$\rho^{\text{stat}}_{i}:(0,1)\to{\mathbb{R}}$ for $i\in I$ if, for any
$y\in(0,1)$,

(3.36) $\lim_{N\to\infty}\theta_{0}^{(N)}(\lceil
yN\rceil)=\rho^{\text{stat}}_{0}(y),\quad\lim_{N\to\infty}\theta_{1}^{(N)}(\lceil
yN\rceil)=\rho^{\text{stat}}_{1}(y).$

######  Theorem 3.9.

[Stationary macroscopic profile] Let
$(\theta_{0}^{(N)}(x),\theta_{1}^{(N)}(x))_{x\in V}$ be the stationary
microscopic profile (see (3.9)) for the process $\\{\eta(t)\colon\,t\geq
0\\},\ \eta(t)=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V}$ with generator
$L_{\epsilon,\gamma_{N},N}$ (see (3.1)), where $\gamma_{N}$ is such that
$\gamma_{N}N^{2}\to\Upsilon$ as $N\to\infty$ for some $\Upsilon>0$. Then, for
each $y\in(0,1)$, the pointwise limits (see Fig. 3)

(3.37) ${\rho_{0}^{\rm{stat}}}(y):=\lim_{N\to\infty}\theta_{0}^{(N)}(\lceil
yN\rceil),\quad{\rho_{1}^{\rm{stat}}}(y):=\lim_{N\to\infty}\theta_{1}^{(N)}(\lceil
yN\rceil),$

exist and are given by

(3.38) $\displaystyle{\rho_{0}^{\rm{stat}}}(y)$
$\displaystyle=\rho_{L,0}+(\rho_{R,0}-\rho_{L,0})y,$ $\displaystyle
y\in(0,1),$ $\displaystyle{\rho_{1}^{\rm{stat}}}(y)$
$\displaystyle={\rho_{0}^{\rm{stat}}}(y),$ $\displaystyle y\in(0,1),$

when $\epsilon=0$, while

(3.39) $\displaystyle{\rho_{0}^{\rm{stat}}}(y)$
$\displaystyle=\frac{\epsilon}{1+\epsilon}\left[\frac{\sinh\left[B_{\epsilon,\Upsilon}\,(1-y)\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(L,0)}-\rho_{(L,1)})+\frac{\sinh\left[B_{\epsilon,\Upsilon}\,y\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(R,0)}-\rho_{(R,1)})\right]$
$\displaystyle+\frac{1}{1+\epsilon}\left[\rho_{(R,0)}\,y+\rho_{(L,0)}\,(1-y)\right]+\frac{\epsilon}{1+\epsilon}\left[\rho_{(R,1)}\,y+\rho_{(L,1)}\,(1-y)\right],$
(3.40) $\displaystyle{\rho_{1}^{\rm{stat}}}(y)$
$\displaystyle=\frac{1}{1+\epsilon}\left[\frac{\sinh\left[B_{\epsilon,\Upsilon}\,(1-y)\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(L,1)}-\rho_{(L,0)})+\frac{\sinh\left[B_{\epsilon,\Upsilon}\,y\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(R,1)}-\rho_{(R,0)})\right]$
$\displaystyle+\frac{1}{1+\epsilon}\left[\rho_{(R,0)}\,y+\rho_{(L,0)}\,(1-y)\right]+\frac{\epsilon}{1+\epsilon}\left[\rho_{(R,1)}\,y+\rho_{(L,1)}\,(1-y)\right],$

when $\epsilon>0$, where
$B_{\epsilon,\Upsilon}:=\sqrt{\Upsilon(1+\tfrac{1}{\epsilon})}$. Moreover,
when $\epsilon>0$, the two limits in (3.37) are uniform in $(0,1)$.

###### Proof.

For $\epsilon=0$, it easily follows from (3.21) plus the fact that
$\gamma_{N}N^{2}\to\Upsilon>0$ and $\tfrac{\lceil yN\rceil}{N}\to y$ uniformly
in $(0,1)$ as $N\to\infty$, that

$\lim_{N\to\infty}\sup_{y\in(0,1)}\left|\theta_{0}^{(N)}(\lceil
yN\rceil)-[\rho_{(L,0)}+(\rho_{(R,0)}-\rho_{(L,0)})\,y]\right|=0,$

and since $\theta_{1}(x)=\theta_{0}(x)$ for all $x\in\\{2,\ldots,N-1\\}$, for
fixed $y\in(0,1)$, we have

$\lim_{N\to\infty}\theta_{1}^{(N)}(\lceil
yN\rceil)={\rho_{0}^{\rm{stat}}}(y).$

When $\epsilon>0,$ since $\gamma_{N}N^{2}\to\Upsilon>0$ as $N\to\infty,$ we
note the following:

(3.41) $\displaystyle\gamma_{N}\overset{N\to\infty}{\longrightarrow}0,$
$\displaystyle\lim\limits_{N\to\infty}\alpha_{1}=\lim\limits_{N\to\infty}\alpha_{2}=1,$
$\displaystyle\lim\limits_{N\to\infty}\alpha_{1}^{N}={\rm
e}^{-B_{\epsilon,\Upsilon}},\quad\lim\limits_{N\to\infty}\alpha_{2}^{N}={\rm
e}^{B_{\epsilon,\Upsilon}}.$

Consequently, from the expressions of $(\vec{c}_{i})_{1\leq i\leq 4}$ defined
in (A.4), we also have

(3.42) $\displaystyle\lim\limits_{N\to\infty}N\vec{c}_{1}$
$\displaystyle=\frac{1}{1+\epsilon}\big{[}\begin{array}[]{c c c
c}-1&-\epsilon&1&\epsilon\end{array}\big{]}^{T},\
\lim\limits_{N\to\infty}\vec{c}_{2}=\frac{1}{1+\epsilon}\big{[}\begin{array}[]{c
c c c}1&\epsilon&0&0\end{array}\big{]}^{T},$
$\displaystyle\lim\limits_{N\to\infty}\vec{c}_{3}$
$\displaystyle=\frac{1}{1+\epsilon}\big{[}\begin{array}[]{c c c c}\frac{{\rm
e}^{B_{\epsilon,\Upsilon}}}{{\rm e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}&-\frac{{\rm e}^{B_{\epsilon,\Upsilon}}}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm e}^{-B_{\epsilon,\Upsilon}}}&-\frac{1}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm e}^{-B_{\epsilon,\Upsilon}}}&\frac{1}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}\end{array}\big{]}^{T},$
$\displaystyle\lim\limits_{N\to\infty}\vec{c}_{4}$
$\displaystyle=\frac{1}{1+\epsilon}\big{[}\begin{array}[]{c c c c}-\frac{{\rm
e}^{-B_{\epsilon,\Upsilon}}}{{\rm e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}&\frac{{\rm e}^{-B_{\epsilon,\Upsilon}}}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm e}^{-B_{\epsilon,\Upsilon}}}&\frac{1}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm e}^{-B_{\epsilon,\Upsilon}}}&-\frac{1}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}\end{array}\big{]}^{T}.$

Combining the above equations with (3.33), and the fact that $\tfrac{\lceil
yN\rceil}{N}\to y$ uniformly in $(0,1)$ as $N\to\infty$, we get the desired
result. ∎

######  Remark 3.10.

[Non-uniform convergence] Note that for $\epsilon>0$ both stationary
macroscopic profiles, when extended continuously to the closed interval
$[0,1]$, match the prescribed boundary conditions. This is different from what
happens for $\epsilon=0$, where the continuous extension of
${\rho_{1}^{\rm{stat}}}$ to the closed interval $[0,1]$ equals
${\rho_{0}^{\rm{stat}}}(y)=\rho_{L,0}+(\rho_{R,0}-\rho_{L,0})y$, which does
not necessarily match the prescribed boundary conditions unless
$\rho_{(L,1)}=\rho_{(L,0)}$ and $\rho_{(R,1)}=\rho_{(R,0)}$. Moreover, as can
be seen from the proof above, for $\epsilon>0,$ the convergence of
$\theta_{i}$ to $\rho_{i}$ is uniform in $[0,1]$, i.e.,

$\lim_{N\to\infty}\sup_{y\in[0,1]}\left|\,{\rho_{0}^{\rm{stat}}}(y)-\theta_{0}^{(N)}(\lceil
yN\rceil)\,\right|=0,\quad\lim_{N\to\infty}\sup_{y\in[0,1]}\left|\,{\rho_{1}^{\rm{stat}}}(y)-\theta_{1}^{(N)}(\lceil
yN\rceil)\,\right|=0,$

while for $\epsilon=0$, the convergence of $\theta_{1}$ to $\rho_{1}$ is not
uniform in $[0,1]$ when either $\rho_{(L,0)}\neq\rho_{(L,1)}$ or
$\rho_{(R,0)}\neq\rho_{(R,1)}$.

Also, if $\rho_{i}^{\text{stat},\epsilon}(\cdot)$ denotes the macroscopic
profile defined in (3.39)$-$(3.40), then for $\epsilon>0$ and $i\in\\{0,1\\}$,
we have

(3.43) $\lim\limits_{\epsilon\to
0}\rho_{i}^{\text{stat},\epsilon}(y)\to\rho_{i}^{\text{stat},0}(y)$

for fixed $y\in(0,1)$ and $i\in\\{0,1\\}$, where
$\rho_{i}^{\text{stat},0}(\cdot)$ is the corresponding macroscopic profile in
(3.38) for $\epsilon=0$. However, this convergence is also not uniform for
$i=1$ when $\rho_{(L,0)}\neq\rho_{(L,1)}$ or $\rho_{(R,0)}\neq\rho_{(R,1)}$.
$\spadesuit$

In view of the considerations in Remark 3.10, we next concentrate on the case
$\epsilon>0.$ The following result tells us that for $\epsilon>0$ the
stationary macroscopic profiles satisfy a stationary PDE with fixed boundary
conditions and also admit a stochastic representation in terms of an absorbing
switching diffusion process.

######  Theorem 3.11.

[Stationary boundary value problem] Consider the boundary value problem

(3.44) $\begin{cases}0=\Delta u_{0}+\Upsilon(u_{1}-u_{0}),\\\ 0=\epsilon\Delta
u_{1}+\Upsilon(u_{0}-u_{1}),\end{cases}$

with boundary conditions

(3.45) $\begin{cases}u_{0}(0)=\rho_{L,0},\ u_{0}(1)=\rho_{R,0},\\\
u_{1}(0)=\rho_{L,1},\ u_{1}(1)=\rho_{R,1},\end{cases}$

where $\epsilon,\Upsilon>0$, and the four boundary parameters
$\rho_{(L,0)},\,\rho_{(L,1)},\,\rho_{(R,0)},\,\rho_{(R,1)}$ are also positive.
Then the PDE admits a unique strong solution given by

(3.46) $u_{i}(y)={\rho_{i}^{\rm{stat}}}(y),\qquad y\in[0,1],$

where $({\rho_{0}^{\rm{stat}}}(\cdot),{\rho_{1}^{\rm{stat}}}(\cdot))$ are as
defined in (3.37). Furthermore,
$({\rho_{0}^{\rm{stat}}}(\cdot),{\rho_{1}^{\rm{stat}}}(\cdot))$ has the
stochastic representation

(3.47)
${\rho_{i}^{\rm{stat}}}(y)={\mathbb{E}}_{(y,i)}[\Phi_{i_{\tau}}(X_{\tau})],$

where $\\{i_{t}\colon\,t\geq 0\\}$ is the pure jump process on state space
$I=\\{0,1\\}$ that switches at rate $\Upsilon$, the functions
$\Phi_{0},\Phi_{1}\colon\,I\to{\mathbb{R}}_{+}$ are defined as

$\Phi_{0}=\rho_{(L,0)}\,\boldsymbol{1}_{\\{0\\}}+\rho_{(R,0)}\,\boldsymbol{1}_{\\{1\\}},\qquad\Phi_{1}=\rho_{(L,1)}\,\boldsymbol{1}_{\\{0\\}}+\rho_{(R,1)}\,\boldsymbol{1}_{\\{1\\}},$

$\\{X_{t}\colon\,t\geq 0\\}$ is the stochastic process $[0,1]$ that satisfies
the SDE

${\rm d}X_{t}=\psi(i_{t})\,{\rm d}W_{t}$

with $W_{t}=B_{2t}$ and $\\{B_{t}\colon\,t\geq 0\\}$ standard Brownian motion,
the switching diffusion process $\\{(X_{t},i_{t})\colon\,t\geq 0\\}$ is killed
at the stopping time

$\tau:=\inf\\{t\geq 0\,\colon\,X_{t}\in I\\},$

and $\psi\colon\,I\to\\{1,\epsilon\\}$ is given by
$\psi:=\boldsymbol{1}_{\\{0\\}}+\epsilon\,\boldsymbol{1}_{\\{1\\}}$.

###### Proof.

It is straightforward to verify that for $\epsilon>0$ the macroscopic profiles
$\rho_{0},\rho_{1}$ defined in (3.39)$-$(3.40) are indeed uniformly continuous
in $(0,1)$ and thus can be uniquely extended continuously to $[0,1]$, namely,
by defining
${\rho_{i}^{\rm{stat}}}(0)=\rho_{(L,i)},\,{\rho_{i}^{\rm{stat}}}(1)=\rho_{(R,i)}$
for $i\in I$. Also ${\rho_{i}^{\rm{stat}}}\in C^{\infty}([0,1])$ for $i\in I$
and satisfy the stationary PDE (3.44), with the boundary conditions specified
in (3.45).

The stochastic representation of a solution of the system in (3.44) follows
from [29, p385, Eq.(4.7)]. For the sake of completeness, we give the proof of
uniqueness of the solution of (3.44). Let $u=(u_{0},u_{1})$ and
$v=(v_{0},v_{1})$ be two solutions of the stationary reaction diffusion
equation with the specified boundary conditions in (3.45). Then
$(w_{0},w_{1}):=(u_{0}-v_{0},u_{1}-v_{1})$ satisfies

(3.48) $\begin{cases}0=\Delta w_{0}+\Upsilon(w_{1}-w_{0}),\\\ 0=\epsilon\Delta
w_{1}+\Upsilon(w_{0}-w_{1}),\end{cases}$

with boundary conditions

(3.49) $w_{0}(0)=w_{0}(1)=w_{1}(0)=w_{1}(1)=0.$

Multiplying the two equations in (3.48) with $w_{0}$ and $w_{1}$,
respectively, and using the identity

$w_{i}\Delta w_{i}=\nabla\cdot(w_{i}\nabla w_{i})-|\nabla w_{i}|^{2},\quad
i\in I,$

we get

(3.50) $\begin{cases}0=\nabla\cdot(w_{0}\nabla w_{0})-|\nabla
w_{0}|^{2}+\Upsilon(w_{1}-w_{0})w_{0},\\\ 0=\epsilon\nabla\cdot(w_{1}\nabla
w_{1})-\epsilon|\nabla w_{1}|^{2}+\Upsilon(w_{0}-w_{1})w_{1}.\end{cases}$

Integrating both equations by parts over $[0,1],$ we get

(3.51) $\displaystyle 0$ $\displaystyle=-[w_{0}(1)\nabla
w_{0}(1)-w_{0}(0)\nabla w_{0}(0)]-\int_{0}^{1}dy\,|\nabla
w_{0}(y)|^{2}+\Upsilon\int_{0}^{1}dy\,(w_{1}(y)-w_{0}(y))w_{0}(y),$
$\displaystyle 0$ $\displaystyle=-\epsilon[w_{1}(1)\nabla
w_{1}(1)-w_{1}(0)\nabla w_{1}(0)]-\epsilon\int_{0}^{1}dy\,|\nabla
w_{1}(y)|^{2}+\Upsilon\int_{0}^{1}dy\,(w_{0}(y)-w_{1}(y))w_{1}(y).$

Adding the above two equations and using the zero boundary conditions in
(3.49), we have

(3.52) $\int_{0}^{1}dy\,|\nabla w_{0}(y)|^{2}+\epsilon\int_{0}^{1}dy\,|\nabla
w_{1}(y)|^{2}+\Upsilon\int_{0}^{1}dy\,[w_{1}(y)-w_{0}(y)]^{2}=0.$

Since both $w_{0}$ and $w_{1}$ are continuous and $\epsilon>0,\Upsilon>0$, it
follows that

(3.53) $w_{0}=w_{1},\quad\nabla w_{0}=\nabla w_{1}=0,$

and so $w_{0}=w_{1}\equiv 0$. ∎

Note that, as a result of Theorem 3.11, the four absorption probabilities of
the switching diffusion process $\\{(X_{t},i_{t})\,:\,t\geq 0\\}$ starting
from $(y,i)\in[0,1]\times I$ are indeed the respective coefficients of
$\rho_{(L,0)},\rho_{(L,1)},\rho_{(R,0)}$, $\rho_{(R,1)}$ appearing in the
expression of ${\rho_{i}^{\rm{stat}}}(y)$. Furthermore note that, as a
consequence of Theorem 3.11 and the results in [33, Section 3], the time-
dependent boundary-value problem

(3.54)
$\begin{cases}\partial_{t}\rho_{0}=\Delta\rho_{0}+\Upsilon(\rho_{1}-\rho_{0}),\\\
\partial_{t}\rho_{1}=\epsilon\Delta\rho_{1}+\Upsilon(\rho_{0}-\rho_{1}),\end{cases}$

with initial conditions

(3.55) $\begin{cases}\rho_{0}(x,0)=\bar{\rho}_{0}(x),\\\
\rho_{1}(x,0)=\bar{\rho}_{1}(x),\end{cases}$

and boundary conditions

(3.56) $\begin{cases}\rho_{0}(0,t)=\rho_{L,0},\ \rho_{0}(1,t)=\rho_{R,0},\\\
\rho_{1}(0,t)=\rho_{L,1},\ \rho_{1}(1,t)=\rho_{R,1},\end{cases}$

admits a unique solution given by

(3.57)
$\begin{cases}\rho_{0}(x,t)=\rho^{\text{hom}}_{0}(x,t)+{\rho_{0}^{\rm{stat}}}(x),\\\
\rho_{1}(x,t)=\rho^{\text{hom}}_{1}(x,t)+{\rho_{1}^{\rm{stat}}}(x),\end{cases}$

where

(3.58) $\rho_{0}^{\text{hom}}(x,t)={\rm e}^{-\Upsilon
t}h_{0}(x,t)+\frac{\Upsilon}{1-\epsilon}{\rm e}^{-\Upsilon t}\int_{\epsilon
t}^{t}{\rm d}s\,\left(\left(\frac{s-\epsilon
t}{t-s}\right)^{1/2}I_{1}(\upsilon(s))h_{0}(x,s)+I_{0}(\upsilon(s))h_{1}(x,s)\right),$
(3.59) $\rho_{1}^{\text{hom}}(x,t)={\rm e}^{-\Upsilon t}h_{1}(x,\epsilon
t)+\frac{\Upsilon}{1-\epsilon}{\rm e}^{-\Upsilon t}\int_{\epsilon t}^{t}{\rm
d}s\,\left(\left(\frac{s-\epsilon
t}{t-s}\right)^{-1/2}I_{1}(\upsilon(s))h_{1}(x,s)+I_{0}(\upsilon(s))h_{0}(x,s)\right),$

$\upsilon(s)=\frac{2\Upsilon}{1-\epsilon}((t-s)(s-\epsilon t))^{1/2}$,
$I_{0}(\cdot)$ and $I_{1}(\cdot)$ are the modified Bessel functions,
$h_{0}(x,t)$, $h_{1}(x,t)$ are the solutions of

(3.60) $\begin{cases}\partial_{t}h_{0}=\Delta h_{0},\\\
\partial_{t}h_{1}=\Delta h_{1},\\\
h_{0}(x,0)=\bar{\rho}_{0}(x)-{\rho_{0}^{\rm{stat}}}(x),\\\
h_{1}(x,1)=\bar{\rho}_{1}(x)-{\rho_{1}^{\rm{stat}}}(x),\\\
h_{0}(0,t)=h_{0}(1,t)=h_{1}(0,t)=h_{1}(1,t)=0,\end{cases}$

and ${\rho_{0}^{\rm{stat}}}(x)$, ${\rho_{1}^{\rm{stat}}}(x)$ are given in
(3.40).

We conclude this section by proving that the solution of the time-dependent
boundary-value problem in (3.54) converges to the stationary profile in
(3.40).

######  Proposition 3.12.

[Convergence to stationary profile] Let ${\rho_{0}^{\rm{hom}}}(x,t)$ and
${\rho_{1}^{\rm{hom}}}(x,t)$ be as in (3.58) and (3.59), respectively, i.e.,
the solutions of the boundary-value problem (3.54) with zero boundary
conditions and initial conditions given by
${\rho_{0}^{\rm{hom}}}(x,0)=\bar{\rho}_{0}(x)-{\rho_{0}^{\rm{stat}}}(x)$ and
${\rho_{1}^{\rm{hom}}}(x,0)=\bar{\rho}_{1}(x)-{\rho_{1}^{\rm{stat}}}(x)$.
Then, for any $k\in{\mathbb{N}}$,

$\lim_{t\to\infty}\Big{[}\|{\rho_{0}^{\rm{hom}}}(x,t)\|_{C^{k}(0,1)}+\|{\rho_{1}^{\rm{hom}}}(x,t)\|_{C^{k}(0,1)}\Big{]}=0.$

###### Proof.

We start by showing that

(3.61)
$\lim\limits_{t\to\infty}\,\Big{[}\|{\rho_{0}^{\rm{hom}}}(x,t)\|_{L^{2}(0,1)}+\|{\rho_{1}^{\rm{hom}}}(x,t)\|_{L^{2}(0,1)}\Big{]}=0.$

Multiply the first equation of (3.54) by $\rho_{0}$ and the second equation by
$\rho_{1}$. Integration by parts yields

(3.62) $\begin{cases}\partial_{t}\left(\int_{0}^{1}{\rm
d}x\,\rho_{0}^{2}\right)=-\int_{0}^{1}{\rm
d}x\,|\partial_{x}\rho_{0}|^{2}+\Upsilon\int_{0}^{1}{\rm
d}x\,(\rho_{1}\rho_{0}-\rho_{0}^{2}),\\\ \partial_{t}\left(\int_{0}^{1}{\rm
d}x\,\rho_{1}^{2}(x,t)\right)=-\epsilon\int_{0}^{1}{\rm
d}x\,|\partial_{x}\rho_{1}|^{2}+\Upsilon\int_{0}^{1}{\rm
d}x\,(\rho_{0}\rho_{1}-\rho_{1}^{2}).\end{cases}$

Summing the two equations and defining $E(t):=\int_{0}^{1}{\rm
d}x\,(\rho_{0}^{2}+\rho_{1}^{2})$, we obtain

(3.63) $\partial_{t}E(t)=-\left(\int_{0}^{1}{\rm
d}x\,|\partial_{x}\rho_{0}|^{2}+\epsilon\int_{0}^{1}{\rm
d}x\,|\partial_{x}\rho_{1}|^{2}\right)-\Upsilon\int_{0}^{1}{\rm
d}x\,(\rho_{0}-\rho_{1})^{2}.$

The Poincaré inequality implies that $\int_{0}^{1}{\rm
d}x\,|\partial_{x}\rho_{0}|^{2}+\epsilon\int_{0}^{1}{\rm
d}x\,|\partial_{x}\rho_{1}|^{2}\geq C_{p}E(t)$, in particular,
$\partial_{t}E(t)\leq-\epsilon C_{p}E(t)$, from which we obtain

$E(t)\leq{\rm e}^{-C_{p}t}E(0),$

and hence (3.61).

From [42, Theorem 2.1] it follows that

$A:=\left[\begin{array}[]{c c c c}\Delta-\Upsilon&\Upsilon\\\
\Upsilon&\epsilon\Delta-\Upsilon\\\ \end{array}\right],$

with domain $D(A)=H^{2}(0,1)\cap H_{0}^{1}(0,1)$, generates a semigroup
$\\{\mathcal{S}_{t}\colon\,t\geq 0\\}$. If we set
$\vec{\rho}(t)=\mathcal{S}_{t}(\vec{\bar{\rho}}^{\text{hom}})$, with
$\vec{\bar{\rho}}^{\text{hom}}=\vec{\bar{\rho}}-\vec{\rho}^{\text{stat}}$,
then by the semigroup property we have

$\vec{\rho}(t)=\mathcal{S}_{t-1}(\mathcal{S}_{1/k})^{k}(\vec{\bar{\rho}}^{\text{\text{hom}}}),\qquad
t\geq 1,$

and hence
$A^{k}\vec{\rho}(t)=\mathcal{S}_{t-1}(A\mathcal{S}_{1/k})^{k}(\vec{\bar{\rho}}^{\text{hom}})$.
If we set $\vec{p}:=(A\mathcal{S}_{1/k})^{k}(\vec{\bar{\rho}}^{\text{hom}})$,
then we obtain, by [42, Theorem 5.2(d)],

$\|A^{k}\vec{\rho}(t)\|_{L^{2}(0,1)}\leq\|\mathcal{S}_{t-1}\vec{p}\|_{L^{2}(0,1)},$

where $\lim_{t\to\infty}\|\mathcal{S}_{t-1}\vec{p}\|_{L^{2}(0,1)}=0$ by the
first part of the proof. The compact embedding

$D(A^{k})\hookrightarrow H^{2k}(0,1)\hookrightarrow C^{k}(0,1),\qquad
k\in{\mathbb{N}},$

concludes the proof. ∎

###  3.4 The stationary current

In this section we compute the expected current in the non-equilibrium steady
state that is induced by different densities at the boundaries. We consider
the microscopic and macroscopic systems, respectively.

##### Microscopic system.

We start by defining the notion of current. The microscopic currents are
associated with the edges of the underlying two-layer graph. In our setting,
we denote by ${\cal J}^{0}_{x,x+1}(t)$ and ${\cal J}^{1}_{x,x+1}(t)$ the
instantaneous current through the horizontal edge $(x,x+1)$, $x\in V$, of the
bottom layer, respectively, top layer at time $t$. Obviously,

${\cal J}^{0}_{x,x+1}(t)=\eta_{0}(x,t)-\eta_{0}(x+1,t),\qquad\qquad{\cal
J}^{1}_{x,x+1}(t)=\epsilon[\eta_{1}(x,t)-\eta_{1}(x+1,t)].$

We are interested in the stationary currents ${J}^{0}_{x,x+1}$, respectively,
${J}^{1}_{x,x+1}$, which are obtained as

(3.64)
$J^{0}_{x,x+1}={\mathbb{E}}_{stat}[\eta_{0}(x)-\eta_{0}(x+1)],\qquad\qquad
J^{1}_{x,x+1}=\epsilon{\mathbb{E}}_{stat}[\eta_{1}(x)-\eta_{1}(x+1)],$

where ${\mathbb{E}}_{stat}$ denotes expectation w.r.t. the unique invariant
probability measure of the microscopic system $\\{\eta(t)\colon\,t\geq 0\\}$
with $\eta(t)=\\{\eta_{0}(x,t),\eta_{1}(x,t)\\}_{x\in V}$. In other words,
$J^{0}_{x,x+1}$ and $J^{1}_{x,x+1}$ give the average flux of particles of type
$0$ and type $1$ across the bond $(x,x+1)$ due to diffusion.

Of course, the average number of particle at each site varies in time also as
a consequence of the reaction term:

$\displaystyle\frac{d}{dt}{\mathbb{E}}[\eta_{0}(x,t)]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}[{\cal J}^{0}_{x-1,x}(t)-{\cal
J}^{0}_{x,x+1}(t)]+\gamma({\mathbb{E}}[\eta_{1}(x,t)]-{\mathbb{E}}[\eta_{0}(x,t)]),$
$\displaystyle\frac{d}{dt}{\mathbb{E}}[\eta_{1}(x,t)]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}[{\cal J}^{1}_{x-1,x}(t)-{\cal
J}^{1}_{x,x+1}(t)]+\gamma({\mathbb{E}}[\eta_{0}(x,t)]-{\mathbb{E}}[\eta_{1}(x,t)]).$

Summing these equations, we see that there is no contribution of the reaction
part to the variation of the average number of particles at site $x$:

$\frac{d}{dt}{\mathbb{E}}[\eta_{0}(x,t)+\eta_{1}(x,t)]={\mathbb{E}}[{\cal
J}_{x-1,x}(t)-{\cal J}_{x,x+1}(t)].$

The sum

(3.65) $J_{x,x+1}=J^{0}_{x,x+1}+J^{1}_{x,x+1},$

with $J^{0}_{x,x+1}$ and $J^{1}_{x,x+1}$ defined in (3.64), will be called the
stationary current between sites at $x,x+1$, $x\in V$, which is responsible
for the variation of the total average number of particles at each site,
regardless of their type.

######  Proposition 3.13.

[Stationary microscopic current] For $x\in\\{2,\ldots,N-1\\}$ the stationary
currents defined in (3.64) are given by

(3.66)
$J^{0}_{x,x+1}=-\tfrac{1+\gamma}{1+N+2N\gamma}[\rho_{(R,0)}-\rho_{(L,0)}]-\tfrac{\gamma}{1+N+2N\gamma}[\rho_{(R,1)}-\rho_{(L,1)}],\quad
J^{1}_{x,x+1}=0,$

when $\epsilon=0$ and by

(3.67) $\displaystyle J^{0}_{x,x+1}$
$\displaystyle=-\vec{c}_{1}\cdot\vec{\rho}-\epsilon[(\vec{c}_{3}\cdot\vec{\rho})\alpha_{1}^{x}(\alpha_{1}-1)+(\vec{c}_{4}\cdot\vec{\rho})\alpha_{2}^{x}(\alpha_{2}-1)],$
$\displaystyle J^{1}_{x,x+1}$
$\displaystyle=-\epsilon\vec{c}_{1}\cdot\vec{\rho}+\epsilon[(\vec{c}_{3}\cdot\vec{\rho})\alpha_{1}^{x}(\alpha_{1}-1)+(\vec{c}_{4}\cdot\vec{\rho})\alpha_{2}^{x}(\alpha_{2}-1)],$

when $\epsilon>0,$ where $\vec{c}_{1},\vec{c}_{3},\vec{c}_{4}$ are the vectors
defined in (A.4) of Appendix A, and $\alpha_{1},\alpha_{2}$ are defined in
(3.31). As a consequence, the current $J_{x,x+1}=J^{0}_{x,x+1}+J^{1}_{x,x+1}$
is independent of $x$ and is given by

(3.68)
$J_{x,x+1}=-\tfrac{1+\gamma}{1+N+2N\gamma}[\rho_{(R,0)}-\rho_{(L,0)}]-\tfrac{\gamma}{1+N+2N\gamma}[\rho_{(R,1)}-\rho_{(L,1)}]$

when $\epsilon=0$ and

(3.69)
$J_{x,x+1}=-(1+\epsilon)\left[C_{1}\,(\rho_{R,0}-\rho_{L,0})+\epsilon\,C_{2}\,(\rho_{R,1}-\rho_{L,1})\right]$

when $\epsilon>0$, where

(3.70) $\displaystyle C_{1}$
$\displaystyle=\frac{[\alpha_{1}(1-\epsilon)(\alpha_{1}^{N-1}-1)+\epsilon\,(\alpha_{1}^{N+1}-1)]}{\alpha_{1}(1-\epsilon)(\alpha_{1}^{N-1}-1)(N+1)+2\epsilon\,(\alpha_{1}^{N+1}-1)(N+\epsilon)},$
$\displaystyle C_{2}$
$\displaystyle=\frac{(\alpha_{1}^{N+1}-1)}{\alpha_{1}(1-\epsilon)(\alpha_{1}^{N-1}-1)(N+1)+2\epsilon\,(\alpha_{1}^{N+1}-1)(N+\epsilon)}.$

###### Proof.

From (3.64) we have

(3.71) $J^{0}_{x,x+1}=\theta_{0}(x)-\theta_{0}(x+1),\qquad\qquad
J^{1}_{x,x+1}=\epsilon[\theta_{1}(x)-\theta_{1}(x+1)],$

where $\theta_{0}(\cdot),\theta_{1}(\cdot)$ are the average microscopic
profiles. Thus, when $\epsilon=0$, the expressions of
$J^{0}_{x,x+1},J^{1}_{x,x+1}$ and consequently $J_{x,x+1}$ follow directly
from (3.21).

For $\epsilon>0$, using the expressions of
$\theta_{0}(\cdot),\theta_{1}(\cdot)$ in (3.33), we see that

(3.72) $\displaystyle J^{0}_{x,x+1}$
$\displaystyle=\theta_{0}(x)-\theta_{0}(x+1)=-\vec{c}_{1}\cdot\vec{\rho}-\epsilon[(\vec{c}_{3}\cdot\vec{\rho})\alpha_{1}^{x}(\alpha_{1}-1)+(\vec{c}_{4}\cdot\vec{\rho})\alpha_{2}^{x}(\alpha_{2}-1)],$
$\displaystyle J^{1}_{x,x+1}$
$\displaystyle=\epsilon[\theta_{1}(x)-\theta_{1}(x+1)]=-\epsilon\vec{c}_{1}\cdot\vec{\rho}+\epsilon[(\vec{c}_{3}\cdot\vec{\rho})\alpha_{1}^{x}(\alpha_{1}-1)+(\vec{c}_{4}\cdot\vec{\rho})\alpha_{2}^{x}(\alpha_{2}-1)],$

where $\vec{c}_{1},\vec{c}_{3},\vec{c}_{4}$ are the vectors defined in (A.4)
of Appendix A, and $\alpha_{1},\alpha_{2}$ are defined in (3.31). Adding the
two equations, we also have

(3.73)
$J_{x,x+1}=J^{0}_{x,x+1}+J^{1}_{x,x+1}=-(1+\epsilon)\,\vec{c}_{1}\cdot\vec{\rho}=(1+\epsilon)\left[C_{1}\,(\rho_{R,0}-\rho_{L,0})+\epsilon\,C_{2}\,(\rho_{R,1}-\rho_{L,1})\right],$

where $C_{1},C_{2}$ are as in (3.70). ∎

##### Macroscopic system.

The microscopic current scales like $1/N$. Indeed, the currents associated to
the two layers in the macroscopic system can be obtained from the microscopic
currents, respectively, by defining

(3.74) $J^{0}(y)=\lim_{N\to\infty}NJ^{0}_{\lfloor yN\rfloor,\lfloor
yN\rfloor+1},\qquad J^{1}(y)=\lim_{N\to\infty}NJ^{1}_{\lfloor
yN\rfloor,\lfloor yN\rfloor+1}.$

Below we justify the existence of the two limits and thereby provide explicit
expressions for the macroscopic currents.

######  Proposition 3.14.

[Stationary macroscopic current] For $y\in(0,1)$ the stationary currents
defined in (3.74) are given by

(3.75) $J^{0}(y)=-\left[(\rho_{R,0}-\rho_{L,0})\right],\quad J^{1}(y)=0,$

when $\epsilon=0$ and by

(3.76) $\displaystyle J^{0}(y)$ $\displaystyle=\frac{\epsilon
B_{\epsilon,\Upsilon}}{1+\epsilon}\left[\frac{\cosh\left[B_{\epsilon,\Upsilon}\,(1-y)\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(L,0)}-\rho_{(L,1)})-\frac{\cosh\left[B_{\epsilon,\Upsilon}\,y\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(R,0)}-\rho_{(R,1)})\right]$
$\displaystyle\quad-\frac{1}{1+\epsilon}\left[(\rho_{(R,0)}-\rho_{(L,0)})+\epsilon(\rho_{(R,1)}-\rho_{(L,1)})\right]$

and

(3.77) $\displaystyle J^{1}(y)$ $\displaystyle=-\frac{\epsilon
B_{\epsilon,\Upsilon}}{1+\epsilon}\left[\frac{\cosh\left[B_{\epsilon,\Upsilon}\,(1-y)\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(L,0)}-\rho_{(L,1)})-\frac{\cosh\left[B_{\epsilon,\Upsilon}\,y\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(R,0)}-\rho_{(R,1)})\right]$
$\displaystyle\quad-\frac{\epsilon}{1+\epsilon}\left[(\rho_{(R,0)}-\rho_{(L,0)})+\epsilon(\rho_{(R,1)}-\rho_{(L,1)})\right]$

when $\epsilon>0$. As a consequence, the total current
$J(y)=J^{0}(y)+J^{1}(y)$ is constant and is given by

(3.78)
$J(y)=-\left[(\rho_{R,0}-\rho_{L,0})+\epsilon\,(\rho_{R,1}-\rho_{L,1})\right].$

###### Proof.

For $\epsilon=0$ the claim easily follows from the expressions of
$J^{0}_{x,x+1},J^{1}_{x,x+1}$ given in (3.66) and the fact that $\gamma_{N}\to
0$ as $N\to\infty$.

When $\epsilon>0$, we first note the following:

(3.79)
$\displaystyle\gamma_{N}N^{2}\overset{N\to\infty}{\longrightarrow}\Upsilon>0,$
$\displaystyle\lim\limits_{N\to\infty}\alpha_{1}=\lim\limits_{N\to\infty}\alpha_{2}=1,$
$\displaystyle\lim\limits_{N\to\infty}N(\alpha_{1}-1)=-B_{\epsilon,\Upsilon},\
\lim\limits_{N\to\infty}N(\alpha_{2}-1)=B_{\epsilon,\Upsilon},$
$\displaystyle\lim\limits_{N\to\infty}\alpha_{1}^{N}={\rm
e}^{-B_{\epsilon,\Upsilon}},\quad\lim\limits_{N\to\infty}\alpha_{2}^{N}={\rm
e}^{B_{\epsilon,\Upsilon}}.$

Consequently, from the expressions for $(\vec{c}_{i})_{1\leq i\leq 4}$ defined
in (A.4), we also have

(3.80) $\displaystyle\lim\limits_{N\to\infty}N\vec{c}_{1}$
$\displaystyle=\frac{1}{1+\epsilon}\,\big{[}\begin{array}[]{c c c
c}-1&-\epsilon&1&\epsilon\end{array}\big{]}^{T},$
$\displaystyle\lim\limits_{N\to\infty}\vec{c}_{3}$
$\displaystyle=\frac{1}{1+\epsilon}\,\big{[}\begin{array}[]{c c c c}\frac{{\rm
e}^{B_{\epsilon,\Upsilon}}}{{\rm e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}&-\frac{{\rm e}^{B_{\epsilon,\Upsilon}}}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm e}^{-B_{\epsilon,\Upsilon}}}&-\frac{1}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm e}^{-B_{\epsilon,\Upsilon}}}&\frac{1}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}\end{array}\big{]}^{T},$
$\displaystyle\lim\limits_{N\to\infty}\vec{c}_{4}$
$\displaystyle=\frac{1}{1+\epsilon}\,\big{[}\begin{array}[]{c c c
c}-\frac{{\rm e}^{-B_{\epsilon,\Upsilon}}}{{\rm
e}^{B_{\epsilon,\Upsilon}}-{\rm e}^{-B_{\epsilon,\Upsilon}}}&\frac{{\rm
e}^{-B_{\epsilon,\Upsilon}}}{{\rm e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}&\frac{1}{{\rm e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}&-\frac{1}{{\rm e}^{B_{\epsilon,\Upsilon}}-{\rm
e}^{-B_{\epsilon,\Upsilon}}}\end{array}\big{]}^{T}.$

Combining the above equations with (3.67), we have

(3.81) $\displaystyle J^{0}(y)$
$\displaystyle=\lim_{N\to\infty}NJ^{0}_{\lfloor yN\rfloor,\lfloor
yN\rfloor+1}$ $\displaystyle=-\epsilon
B_{\epsilon,\Upsilon}\Big{[}\Big{(}\lim\limits_{N\to\infty}\vec{c}_{4}\cdot\vec{\rho}\Big{)}\,{\rm
e}^{B_{\epsilon,\Upsilon}y}-\Big{(}\lim\limits_{N\to\infty}\vec{c}_{3}\cdot\vec{\rho}\Big{)}\,{\rm
e}^{-B_{\epsilon,\Upsilon}y}\Big{]}-\Big{(}\lim\limits_{N\to\infty}N\vec{c}_{1}\cdot\vec{\rho}\Big{)}$
$\displaystyle=\frac{\epsilon
B_{\epsilon,\Upsilon}}{1+\epsilon}\left[\frac{\cosh\left[B_{\epsilon,\Upsilon}\,(1-y)\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(L,0)}-\rho_{(L,1)})-\frac{\cosh\left[B_{\epsilon,\Upsilon}\,y\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(R,0)}-\rho_{(R,1)})\right]$
$\displaystyle\quad-\frac{1}{1+\epsilon}\left[(\rho_{(R,0)}-\rho_{(L,0)})+\epsilon(\rho_{(R,1)}-\rho_{(L,1)})\right]$

and, similarly,

(3.82) $\displaystyle J^{1}(y)$
$\displaystyle=\lim_{N\to\infty}NJ^{1}_{\lfloor yN\rfloor,\lfloor
yN\rfloor+1}$ $\displaystyle=\epsilon
B_{\epsilon,\Upsilon}\Big{[}\Big{(}\lim\limits_{N\to\infty}\vec{c}_{4}\cdot\vec{\rho}\Big{)}\,{\rm
e}^{B_{\epsilon,\Upsilon}y}-\Big{(}\lim\limits_{N\to\infty}\vec{c}_{3}\cdot\vec{\rho}\Big{)}\,{\rm
e}^{-B_{\epsilon,\Upsilon}y}\Big{]}-\epsilon\,\Big{(}\lim\limits_{N\to\infty}N\vec{c}_{1}\cdot\vec{\rho}\Big{)}$
$\displaystyle=-\frac{\epsilon
B_{\epsilon,\Upsilon}}{1+\epsilon}\left[\frac{\cosh\left[B_{\epsilon,\Upsilon}\,(1-y)\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(L,0)}-\rho_{(L,1)})-\frac{\cosh\left[B_{\epsilon,\Upsilon}\,y\right]}{\sinh\left[B_{\epsilon,\Upsilon}\right]}\,(\rho_{(R,0)}-\rho_{(R,1)})\right]$
$\displaystyle\quad-\frac{\epsilon}{1+\epsilon}\left[(\rho_{(R,0)}-\rho_{(L,0)})+\epsilon(\rho_{(R,1)}-\rho_{(L,1)})\right].$

Adding $J^{0}(y)$ and $J^{1}(y)$, we obtain the total current

(3.83)
$J(y)=J^{0}(y)+J^{1}(y)=-\left[(\rho_{R,0}-\rho_{L,0})+\epsilon\,(\rho_{R,1}-\rho_{L,1})\right],$

which is indeed independent of $y$. ∎

######  Remark 3.15.

[Currents] Combining the expressions for the density profiles and the current,
we see that

$J^{0}(y)=-\frac{d\rho_{0}}{dy}(y),\qquad
J^{1}(y)=-\epsilon\frac{d\rho_{1}}{dy}(y).$

$\spadesuit$

###  3.5 Discussion: Fick’s law and uphill diffusion

In this section we discuss the behaviour of the boundary-driven system as the
parameter $\epsilon$ is varied. For simplicity we restrict our discussion to
the macroscopic setting, although similar comments hold for the microscopic
system as well.

In view of the previous results, we can rewrite the equations for the
densities $\rho_{0}(y,t),\rho_{1}(y,t)$ as

$\begin{cases}\partial_{t}\rho_{0}=-\nabla
J^{0}+\Upsilon(\rho_{1}-\rho_{0}),\\\ \partial_{t}\rho_{1}=-\nabla
J^{1}+\Upsilon(\rho_{0}-\rho_{1}),\\\ J_{0}=-\nabla\rho_{0},\\\
J_{1}=-\epsilon\nabla\rho_{1},\end{cases}$

which are complemented with the boundary values (for $\epsilon>0$)

$\begin{cases}\rho_{0}(0,t)=\rho_{L,0},\ \rho_{0}(1,t)=\rho_{R,0},\\\
\rho_{1}(0,t)=\rho_{L,1},\ \rho_{1}(1,t)=\rho_{R,1}.\end{cases}$

We will be concerned with the total density $\rho=\rho_{0}+\rho_{1}$, whose
evolution equation does not contain the reaction part, and is given by

(3.84) $\begin{cases}\partial_{t}\rho=-\nabla J,\\\
J=-\nabla(\rho_{0}+\epsilon\rho_{1}),\end{cases}$

with boundary values

(3.85) $\begin{cases}\rho(0,t)=\rho_{L}=\rho_{L,0}+\rho_{R,0},\\\
\rho(1,t)=\rho_{R}=\rho_{R,0}+\rho_{R,1}.\end{cases}$

##### Non-validity of Fick’s law.

From (3.84) we immediately see that Fick’s law of mass transport is satisfied
if and only if $\epsilon=1$. When we allow diffusion and reaction of slow and
fast particles, i.e., $0\leq\epsilon<1$, Fick’s law breaks down, since the
current associated to the total mass is not proportional to the gradient of
the total mass. Rather, the current $J$ is the sum of a contribution $J^{0}$
due to the diffusion of fast particles of type $0$ (at rate 1) and a
contribution $J^{1}$ due to the diffusion of slow particles of type $1$ (at
rate $\epsilon$). Interestingly, the violation of Fick’s law opens up the
possibility of several interesting phenomena that we discuss in what follows.

##### Equal boundary densities with non-zero current.

In a system with diffusion and reaction of slow and fast particles we may
observe a non-zero current when the total density has the same value at the
two boundaries. This is different from what is observed in standard diffusive
systems driven by boundary reservoirs, where in order to have a stationary
current it is necessary that the reservoirs have different chemical
potentials, and therefore different densities, at the boundaries.

Let us, for instance, consider the specific case when
$\rho_{L,0}=\rho_{R,1}=2$ and $\rho_{L,1}=\rho_{R,0}=4$, which indeed implies
equal densities at the boundaries given by $\rho_{L}=\rho_{R}=6$. The density
profiles and currents are displayed in Fig. 3 for two values of $\epsilon$,
which shows the comparison between the Fick-regime $\epsilon=1$ (left panels)
and the non-Fick-regime with very slow particles $\epsilon=0.001$ (right
panels).

On the one hand, in the Fick-regime the profile of both types of particles
interpolates between the boundary values, with a slightly non-linear shape
that has been quantified precisely in (3.39)–(3.40). Furthermore, in the same
regime $\epsilon=1$, the total density profile is flat and the total current
$J$ vanishes because $J^{0}(y)=-J^{1}(y)$ for all $y\in[0,1]$.

On the other hand, in the non-Fick-regime with $\epsilon=0.001$, the
stationary macroscopic profile for the fast particles interpolates between the
boundary values almost linearly (see (3.43)), whereas the profile for the slow
particles is non-monotone: it has two bumps at the boundaries and in the bulk
closely follows the other profile. This non-monotonicity in the profile of the
slow particles is due to the non-uniform convergence in the limit
$\epsilon\downarrow 0$, as pointed out in the last part of Remark 3.10. As a
consequence, the total density profile is not flat and has two bumps at the
boundaries. Most strikingly, the total current is $J=-2$, since now the
current of the bottom layer $J^{0}$ is dominating, while the current of the
bottom layer $J^{1}$ is small (order $\epsilon$).

Figure 3: Macroscopic profiles of the densities for slow and fast particles
(top panels), macroscopic profile of the total density (central panels), and
the currents (bottom panels). Here,
$\rho_{(L,0)}=2,\,\rho_{(L,1)}=4,\,\rho_{(R,0)}=4$ and
$\rho_{(R,1)}=2,\Upsilon=1$. For the panels in the left column, $\epsilon=1$
and for the panels in the right column, $\epsilon=0.001$.

##### Unequal boundary densities with uphill diffusion.

As argued earlier, since the system does not always obey Fick’s law, by tuning
the parameters $\rho_{(L,0)},\rho_{(L,1)},\rho_{(R,0)},\rho_{(R,1)}$ and
$\epsilon$, we can push the system into a regime where the total current is
such that $J<0$ and the total densities are such that $\rho_{R}<\rho_{L}$,
where $\rho_{R}=\rho_{(R,0)}+\rho_{(R,1)}$ and
$\rho_{L}=\rho_{(L,0)}+\rho_{(L,1)}$. In this regime, the current goes uphill,
since the total density of particles at the right is lower than at the left,
yet the average current is negative.

For an illustration, consider the case when $\rho_{L,1}=6,\rho_{R,0}=4$ and
$\rho_{L,0}=\rho_{R,1}=2$, which implies $\rho_{L}=8$ and $\rho_{R}=6$ and
thus $\rho_{R}<\rho_{L}$. The density profiles and currents are shown in Fig.
4 for two values of $\epsilon$, in particular, a comparison between the Fick-
regime $\epsilon=1$ (left panels) and the non-Fick-regime with very slow
particles $\epsilon=0.001$ (right panels). As can be seen in the figure, when
$\epsilon=1,$ the system obeys Fick’s law: the total density linearly
interpolates between the two total boundary densities 8 and 6, respectively.
The average total stationary current is positive as predicted by Fick’s law.
However, in the _uphill_ regime, the total density is non-linear and the
gradient of the total density is not proportional to the total current,
violating Fick’s law. The total current is negative and is effectively
dominated by the current of the fast particles. It will be shown later that
the transition into the uphill regime happens at the critical value
$\epsilon=\tfrac{|\rho_{(R,0)}-\rho_{(L,0)}|}{|\rho_{(R,1)}-\rho_{(L,1)}|}=\tfrac{1}{2}$.
In the limit $\epsilon\downarrow 0$ the total density profile and the current
always get dominated in the bulk by the profile and the current of the fast
particles, respectively. When $\epsilon<\tfrac{1}{2}$, even though the density
of the slow particles makes the total density near the boundaries such that
$\rho_{R}<\rho_{L}$, it is not strong enough to help the system overcome the
domination of the fast particles in the bulk, and so the effective total
current goes in the same direction as the current of the fast particles,
producing an uphill current.

Figure 4: Macroscopic profiles of the densities for slow and fast particles
(top panels), macroscopic profile of the total density (central panels), and
the currents (bottom panels). Here,
$\rho_{(L,0)}=2,\,\rho_{(L,1)}=6,\,\rho_{(R,0)}=4$ and
$\rho_{(R,1)}=2,\Upsilon=1$. For the panels in the left column, $\epsilon=1$
and for the panels in the right column, $\epsilon=0.001$.

##### The transition between downhill and uphill.

We observe that for the choice of reservoir parameters
$\rho_{L,1}=6,\rho_{R,0}=4$ and $\rho_{L,0}=\rho_{R,1}=2$, the change from
downhill to uphill diffusion occurs at
$\epsilon=\tfrac{|\rho_{(R,0)}-\rho_{(L,0)}|}{|\rho_{(R,1)}-\rho_{(L,1)}|}=\tfrac{1}{2}$.
The density profiles and currents are shown in Fig. 5 for two additional
values of $\epsilon$, one in the “mild” downhill regime $J>0$ for
$\epsilon=0.75$ (left panels), the other in the “mild” uphill regime $J<0$ for
$\epsilon=0.25$ (right panels). In the uphill regime (right panel), i.e., when
$\epsilon=0.75$, the “mild” non-linearity of the total density profile is
already visible, indicating the violation of Fick’s law.

Figure 5: Macroscopic profiles of the densities for slow and fast particles
(top panels), macroscopic profile of the total density (central panels), and
the currents (bottom panels) in the “mild” downhill and the “mild” uphill
regime. Here, $\rho_{(L,0)}=2,\,\rho_{(L,1)}=6,\,\rho_{(R,0)}=4$ and
$\rho_{(R,1)}=2,\Upsilon=1$. For the panels in the left column,
$\epsilon=0.75$ and for the panels in the right column, $\epsilon=0.25$.

##### Identification of the uphill regime.

We define the notion of uphill current below and identify the parameter ranges
for which uphill diffusion occurs.

######  Definition 3.16.

[Uphill diffusion] For parameters
$\rho_{(L,0)},\rho_{(L,1)},\rho_{(R,0)},\rho_{(R,1)}$ and $\epsilon>0,$ we say
the system has an uphill current in stationarity if the total current $J$ and
the difference between the total density of particles in the right and the
left side of the system given by $\rho_{R}-\rho_{L}$ have the same sign, where
it is understood that $\rho_{R}=\rho_{(R,0)}+\rho_{(R,1)}$ and
$\rho_{L}=\rho_{(L,0)}+\rho_{(L,1)}$. $\spadesuit$

######  Proposition 3.17.

[Uphill regime] Let $a_{0}:=\rho_{(R,0)}-\rho_{(L,0)}$ and
$a_{1}:=\rho_{(R,1)}-\rho_{(L,1)}$. Then the macroscopic system admits an
uphill current in stationarity if and only if

(3.86) $a_{0}^{2}+(1+\epsilon)\,a_{0}a_{1}+\epsilon a_{1}^{2}<0.$

If, furthermore, $\epsilon\in[0,1]$, then

  * (i)

either

$a_{0}+a_{1}>0\text{ with }a_{0}<0,\,a_{1}>0$

or

$a_{0}+a_{1}<0\text{ with }a_{0}>0,\,a_{1}<0,$

  * (ii)

$\epsilon\in\big{[}0,-\tfrac{a_{0}}{a_{1}}\big{]}$.

###### Proof.

Note that, by (3.78), there is an uphill current if and only if $a_{0}+a_{1}$
and $a_{0}+\epsilon a_{1}$ have opposite signs. In other words, this happens
if and only if

$(a_{0}+a_{1})(a_{0}+\epsilon\,a_{1})=a_{0}^{2}+(1+\epsilon)\,a_{0}a_{1}+\epsilon
a_{1}^{2}<0.$

The above constraint forces $a_{0}a_{1}<0$. Further simplification reduces the
parameter regime to the following four cases:

  * •

$a_{0}+a_{1}>0$ with $a_{0}<0,\,a_{1}>0$ and $\epsilon<-\tfrac{a_{0}}{a_{1}}$,

  * •

$a_{0}+a_{1}<0$ with $a_{0}>0,\,a_{1}<0$ and $\epsilon<-\tfrac{a_{0}}{a_{1}}$,

  * •

$a_{0}+a_{1}>0$ with $a_{0}>0,\,a_{1}<0$ and $\epsilon>-\tfrac{a_{0}}{a_{1}}$,

  * •

$a_{0}+a_{1}<0$ with $a_{0}<0,\,a_{1}>0$ and $\epsilon>-\tfrac{a_{0}}{a_{1}}$.

Under the assumption $\epsilon\in[0,1]$, only the first two of the above four
cases survive. ∎

###  3.6 The width of the boundary layer

We have seen that for $\epsilon=0$ the microscopic density profile of the fast
particles $\theta_{0}(x)$ linearly interpolates between $\rho_{L,0}$ and
$\rho_{R,0}$, whereas the density profile of the slow particles satisfies
$\theta_{1}(x)=\theta_{0}(x)$ for all $x\in\\{2,\ldots,N-1\\}$. In the
macroscopic setting this produces a continuous macroscopic profile
${\rho_{0}^{\rm{stat}}}(y)=\rho_{L,0}+(\rho_{R,0}-\rho_{L,0})y$ for the
bottom-layer, while the top-layer profile develops two discontinuities at the
boundaries when either $\rho_{(L,0)}\neq\rho_{(L,1)}$ or
$\rho_{(R,0)}\neq\rho_{(R,1)}$. In particular,

${\rho_{1}^{\rm{stat}}}(y)\to\big{[}\rho_{L,0}+(\rho_{R,0}-\rho_{L,0})y\big{]}\,\mathbf{1}_{(0,1)}(y)+\rho_{L,1}\mathbf{1}_{\\{1\\}}(y)+\rho_{R,1}\mathbf{1}_{\\{0\\}}(y),\quad\epsilon\downarrow
0,$

for $y\in[0,1]$. For small but positive $\epsilon,$ the curve is smooth and
the discontinuity is turned into a boundary layer. In this section we
investigate the width of the left and the right boundary layers as
$\epsilon\downarrow 0$. To this end, let us define

(3.87) $W_{L}:=|\rho_{(L,0)}-\rho_{(L,1)}|,\quad
W_{R}:=|\rho_{(R,0)}-\rho_{(R,1)}|.$

Note that, the profile $\rho_{1}$ develops a left boundary layer if and only
if $W_{L}>0$ and, similarly, a right boundary layer if and only if $W_{R}>0$.

######  Definition 3.18.

We say that the _left boundary layer_ is of size $f_{L}(\epsilon)$ if there
exists $C>0$ such that, for any $c>0,$

$\lim_{\epsilon\downarrow 0}\frac{R_{L}(\epsilon,c)}{f_{L}(\epsilon)}=C,$

where
$R_{L}(\epsilon,c)=\sup\left\\{y\in\big{(}0,\tfrac{1}{2}\big{)}:\left|\frac{d^{2}}{dy^{2}}{\rho_{1}^{\rm{stat}}}(y)\right|\geq
c\right\\}$. Analogously, we say that the _right boundary layer_ is of size
$f_{R}(\epsilon)$ if there exists $C>0$ such that, for any $c>0,$

$\lim_{\epsilon\downarrow 0}\frac{1-R_{R}(\epsilon,c)}{f_{R}(\epsilon)}=C,$

where
$R_{R}(\epsilon,c)=\inf\left\\{y\in\big{(}\tfrac{1}{2},1\big{)}:\left|\frac{d^{2}}{dy^{2}}{\rho_{1}^{\rm{stat}}}(y)\right|\geq
c\right\\}$.

The widths of the two boundary layers essentially measure the deviation of the
top-layer density profile (and therefore also the total density profile) from
the bulk linear profile corresponding to the case $\epsilon=0$. In the
following proposition we estimate the sizes of the two boundary layers.

######  Proposition 3.19.

[Width of boundary layers] The widths of the two boundary layers are given by

(3.88) $f_{L}(\epsilon)=f_{R}(\epsilon)=\sqrt{\epsilon}\,\log(1/\epsilon),$

where $f_{L}(\epsilon),f_{R}(\epsilon)$ are defined as in Definition 3.18.

###### Proof.

Note that, to compute $f_{L}(\epsilon)$, it suffices to keep $W_{L}>0$ fixed
and put $W_{R}=0$, where $W_{L},W_{R}$ are as in (3.87). Let
$\overline{y}(\epsilon,c)\in(0,\tfrac{1}{2})$ be such that, for some constant
$c>0$,

(3.89) $\left|\frac{d^{2}}{dy^{2}}{\rho_{1}^{\rm{stat}}}(y)\right|\geq c,$

or equivalently, since $\epsilon\Delta\rho_{1}=\Upsilon(\rho_{1}-\rho_{0})$,

(3.90)
$|{\rho_{1}^{\rm{stat}}}(y)-{\rho_{0}^{\rm{stat}}}(y)|\geq\frac{c\epsilon}{\Upsilon}.$

Recalling the expressions of ${\rho_{0}^{\rm{stat}}}(\cdot)$ and
${\rho_{1}^{\rm{stat}}}(\cdot)$ for positive $\epsilon$ given in
(3.39)$-$(3.40), we get

(3.91)
$\left|\frac{\sinh\big{[}\sqrt{\Upsilon(1+\tfrac{1}{\epsilon})}(1-y)\big{]}}{\sinh\big{[}\sqrt{\Upsilon(1+\tfrac{1}{\epsilon})}\big{]}}\,(\rho_{(L,0)}-\rho_{(L,1)})+\frac{\sinh\big{[}\sqrt{\Upsilon(1+\tfrac{1}{\epsilon})}\,y\big{]}}{\sinh\big{[}\sqrt{\Upsilon(1+\tfrac{1}{\epsilon})}\big{]}}\,(\rho_{(R,0)}-\rho_{(R,1)})\right|\geq\frac{c\epsilon}{\Upsilon}.$

Using (3.87) plus the fact that $W_{R}=0,$ and setting
$B_{\epsilon,\Upsilon}:=\sqrt{\Upsilon\left(1+\tfrac{1}{\epsilon}\right)}$, we
see that

(3.92)
$\sinh\left[B_{\epsilon,\Upsilon}(1-y)\right]\geq\frac{c\epsilon}{\Upsilon
W_{L}}\sinh\left[B_{\epsilon,\Upsilon}\right].$

Because $\sinh(\cdot)$ is strictly increasing, (3.92) holds if and only if

(3.93) $\overline{y}(\epsilon,c)\leq
1-\frac{1}{B_{\epsilon,\Upsilon}}\sinh^{-1}\left[\frac{c\epsilon}{\Upsilon
W_{L}}\sinh\left(\tfrac{B_{\epsilon,\Upsilon}}{2}\right)\right].$

Thus, for small $\epsilon>0$ we have

(3.94)
$R_{L}(\epsilon,c)=1-\frac{1}{B_{\epsilon,\Upsilon}}\sinh^{-1}\left[\frac{c\epsilon}{\Upsilon
W_{L}}\sinh\left(\tfrac{B_{\epsilon,\Upsilon}}{2}\right)\right],$

where $R_{L}(\epsilon,c)$ is defined as in Definition 3.18. Since
$\sinh^{-1}x=\log(x+\sqrt{x^{2}+1})$ for $x\in{\mathbb{R}}$, we obtain

(3.95) $\displaystyle R_{L}(\epsilon,c)$
$\displaystyle=\frac{\sqrt{\epsilon}}{\sqrt{\Upsilon(1+\epsilon)}}\log\left[\frac{N_{\epsilon,\Upsilon}+\sqrt{N_{\epsilon,\Upsilon}^{2}+1}}{\epsilon
CN_{\epsilon,\Upsilon}+\sqrt{(\epsilon CN_{\epsilon,\Upsilon})^{2}+1}}\right]$
$\displaystyle=\frac{\sqrt{\epsilon}}{\sqrt{\Upsilon(1+\epsilon)}}\log(1/\epsilon)+\frac{\sqrt{\epsilon}}{\sqrt{\Upsilon(1+\epsilon)}}\log\left[\frac{1+\sqrt{1+(1/N_{\epsilon,\Upsilon})^{2}}}{C+\sqrt{C^{2}+(1/(\epsilon
N_{\epsilon,\Upsilon}))^{2}}}\right]$
$\displaystyle=\frac{\sqrt{\epsilon}}{\sqrt{\Upsilon(1+\epsilon)}}\log(1/\epsilon)+R_{\epsilon,\Upsilon,W_{L}},$

where
$N_{\epsilon,\Upsilon}:=\sinh\Big{(}\tfrac{B_{\epsilon,\Upsilon}}{2}\Big{)},C:=\tfrac{c}{\Upsilon
W_{L}}$, and the error term is

$R_{\epsilon,\Upsilon,W_{L}}:=\frac{\sqrt{\epsilon}}{\sqrt{\Upsilon(1+\epsilon)}}\log\left[\frac{1+\sqrt{1+(1/N_{\epsilon,\Upsilon})^{2}}}{C+\sqrt{C^{2}+(1/(\epsilon
N_{\epsilon,\Upsilon}))^{2}}}\right].$

Note that, since $\epsilon N_{\epsilon,\Upsilon}\to\infty$ as
$\epsilon\downarrow 0$, we have

(3.96) $\lim\limits_{\epsilon\downarrow
0}\frac{R_{\epsilon,\Upsilon,W_{L}}}{\sqrt{\epsilon}}=\frac{1}{\sqrt{\Upsilon}}\log(1/C)<\infty.$

Hence, combining (3.95)$-$(3.96), we get

(3.97) $\lim\limits_{\epsilon\downarrow
0}\frac{R_{L}(\epsilon,c)}{\sqrt{\epsilon}\log(1/\epsilon)}=\lim\limits_{\epsilon\downarrow
0}\frac{1}{\sqrt{\Upsilon}(1+\epsilon)}+\lim\limits_{\epsilon\downarrow
0}\frac{R_{\epsilon,\Upsilon,W_{L}}}{\sqrt{\epsilon}\log(1/\epsilon)}=\frac{1}{\sqrt{\Upsilon}}$

and so, by Definition 3.18, $f_{L}(\epsilon)=\sqrt{\epsilon}\log(1/\epsilon)$.

Similarly, to compute $f_{R}(\epsilon)$, we first fix $W_{L}=0,W_{R}>0$ and
note that, for some $c>0$, we have, by using (3.91),

(3.98) $|\partial^{2}{\rho_{1}^{\rm{stat}}}(y)|\geq c\quad\text{ if and only
if
}\quad\sinh\left[B_{\epsilon,\Upsilon}\,y\right]\geq\tfrac{c\epsilon}{\Upsilon
W_{R}}\sinh\left[B_{\epsilon,\Upsilon}\right].$

Hence, by appealing to the strict monotonicity of $\sinh(\cdot),$ we obtain

(3.99)
$R_{R}(\epsilon,c)=\inf\left\\{y\in\big{(}\tfrac{1}{2},1\big{)}:\left|\frac{d^{2}}{dy^{2}}{\rho_{1}^{\rm{stat}}}(y)\right|\geq
c\right\\}=\frac{1}{B_{\epsilon,\Upsilon}}\sinh^{-1}\left[\frac{c\epsilon}{\Upsilon
W_{R}}\sinh\left(\tfrac{B_{\epsilon,\Upsilon}}{2}\right)\right].$

Finally, by similar computations as in (3.95)–(3.97), we see that

(3.100) $\lim\limits_{\epsilon\downarrow
0}\frac{1-R_{R}(\epsilon,c)}{\sqrt{\epsilon\log(1/\epsilon)}}=\frac{1}{\sqrt{\Upsilon}}$

and hence $f_{R}(\epsilon)=\sqrt{\epsilon}\log(1/\epsilon)$. ∎

##  Appendix A Inverse of the boundary-layer matrix

The inverse of the matrix $M_{\epsilon}$ defined in (3.30) is given by
($\alpha_{1}$ and $\alpha_{2}$ are as in (3.31))

(A.1)
$M_{\epsilon}^{-1}:=\frac{1}{Z}\begin{bmatrix}-m_{13}&-m_{14}&m_{13}&m_{14}\\\
m_{21}&m_{22}&m_{23}&m_{24}\\\
m_{31}(\alpha_{2})&m_{32}(\alpha_{2})&m_{33}(\alpha_{2})&m_{34}(\alpha_{2})\\\
-m_{31}(\alpha_{1})&-m_{32}(\alpha_{1})&-m_{33}(\alpha_{1})&-m_{34}(\alpha_{1})\end{bmatrix},$

where

(A.2) $\displaystyle Z$
$\displaystyle:=\alpha_{1}^{N+1}[\alpha_{2}(1-\epsilon)(\alpha_{2}^{N-1}+1)+2\epsilon(\alpha_{2}^{N+1}+1)]\,[\alpha_{2}(1+N)(1-\epsilon)(\alpha_{2}^{N-1}-1)+2\epsilon(N+\epsilon)(\alpha_{2}^{1+N}-1)],$
$\displaystyle m_{13}$
$\displaystyle:=\alpha_{1}^{N+1}[\alpha_{2}(1-\epsilon)(\alpha_{2}^{N-1}+1)+2\epsilon(\alpha_{2}^{N+1}+1)]\,[\alpha_{2}(1-\epsilon)(\alpha_{2}^{N-1}-1)+\epsilon(\alpha_{2}^{N+1}-1)],$
$\displaystyle m_{14}$
$\displaystyle:=\epsilon\,\alpha_{1}^{N+1}[\alpha_{2}(1-\epsilon)(\alpha_{2}^{N-1}+1)+2\epsilon(\alpha_{2}^{N+1}+1)]\,(\alpha_{2}^{N+1}-1),$
$\displaystyle m_{21}$
$\displaystyle:=(1+N)(1-\epsilon)^{2}(\alpha_{2}^{N-1}-\alpha_{1}^{N-1})-\epsilon(1-\epsilon)^{2}(\alpha_{2}-\alpha_{1})$
$\displaystyle\qquad+\epsilon^{2}(1+2N+\epsilon)(\alpha_{2}^{N+1}-\alpha_{1}^{N+1})+\epsilon(1-\epsilon)(2+3N+\epsilon)(\alpha_{2}^{N}-\alpha_{1}^{N}),$
$\displaystyle m_{22}$
$\displaystyle:=\epsilon\,[(1-\epsilon)(1+N)(\alpha_{2}^{N}-\alpha_{1}^{N})+\epsilon(1+2N+\epsilon)(\alpha_{2}^{N+1}-\alpha_{1}^{N+1})],$
$\displaystyle m_{23}$
$\displaystyle:=\epsilon\,(1-\epsilon)[(N+\epsilon)(\alpha_{2}-\alpha_{1})-(1-\epsilon)(\alpha_{2}^{N}-\alpha_{1}^{N})-\epsilon(\alpha_{2}^{N+1}-\alpha_{1}^{N+1})],$
$\displaystyle m_{24}$
$\displaystyle:=-\epsilon(1-\epsilon)[(1+N)(\alpha_{2}-\alpha_{1})+\epsilon\,(\alpha_{2}^{N+1}-\alpha_{1}^{N+1})],$

and the polynomials $m_{31}(z),m_{32}(z),m_{33}(z),m_{34}(z)$ are defined as

(A.3) $\displaystyle m_{31}(z)$
$\displaystyle:=-(1-\epsilon)^{2}\,z-\epsilon\,(1-\epsilon)+(1-\epsilon)(N+\epsilon)\,z^{N}-\epsilon(1-2N-3\epsilon)\,z^{N+1},$
$\displaystyle m_{32}(z)$
$\displaystyle:=-(1-\epsilon)(1+N)z^{N}-\epsilon\,(1-\epsilon)-\epsilon(1+2N+\epsilon)\,z^{N+1},$
$\displaystyle m_{33}(z)$
$\displaystyle:=(1-\epsilon)^{2}\,z^{N}+\epsilon\,(1-\epsilon)\,z^{N+1}-(1-\epsilon)(N+\epsilon)\,z+\epsilon(1-2N-3\epsilon),$
$\displaystyle m_{34}(z)$
$\displaystyle:=(1+N)(1-\epsilon)\,z+\epsilon\,(1-\epsilon)\,z^{N+1}+\epsilon(1+2N+\epsilon).$

We remark that most of the terms appearing in the inverse simplify because of
(3.32). We define the four vectors
$\vec{c}_{1},\vec{c}_{2},\vec{c}_{3},\vec{c}_{4}$ as the respective rows of
$M_{\epsilon}^{-1}$, i.e.,

(A.4) $\displaystyle\vec{c}_{1}$
$\displaystyle:=(M_{\epsilon}^{-1})^{T}\vec{e}_{1},\quad\vec{c}_{2}:=(M_{\epsilon}^{-1})^{T}\vec{e}_{2},$
$\displaystyle\vec{c}_{3}$
$\displaystyle:=(M_{\epsilon}^{-1})^{T}\vec{e}_{3},\quad\vec{c}_{4}:=(M_{\epsilon}^{-1})^{T}\vec{e}_{4},$

where

$\displaystyle\vec{e}_{1}:=\big{[}\begin{array}[]{c c c
c}1&0&0&0\end{array}\big{]}^{T},$
$\displaystyle\vec{e}_{2}:=\big{[}\begin{array}[]{c c c
c}0&1&0&0\end{array}\big{]}^{T},$
$\displaystyle\vec{e}_{3}:=\big{[}\begin{array}[]{c c c
c}0&0&1&0\end{array}\big{]}^{T},$
$\displaystyle\vec{e}_{4}:=\big{[}\begin{array}[]{c c c
c}0&0&0&1\end{array}\big{]}^{T}.$

## References

  * [1] E.C. Aifantis, A new interpretation of diffusion in high-diffusivity paths – a continuum approach, Acta Metallurgica 27 (1979) 683–691.
  * [2] E.C. Aifantis, J.M. Hill, On the theory of diffusion in media with double diffusivity I. Basic mathematical results, Q. J. Mech. Appl. Math. 33 (1980) 1–21.
  * [3] G. Amir, C. Bahadoran, O. Busani, E. Saada, Invariant measures for multilane exclusion process, preprint available at https://arxiv.org/abs/2105.12974.
  * [4] D.S. Bernstein, W. So, Some explicit formulas for the matrix exponential, IEEE Trans. Autom. Control Vol.38, 8 (1993) 1228–1232.
  * [5] J. Blath, A. González Casanova, N. Kurt, M. Wilke-Berenguer, The seed bank coalescent with simultaneous switching, Electron. J. Probab. 25 (2020) 1–21.
  * [6] J. Blath, N. Kurt, Population genetic models of dormancy, in: Probabilistic Structures in Evolution (eds. E. Baake and A. Wakolbinger) EMS Series of Congress Reports, European Mathematical Society Publishing House (2021) 247–265.
  * [7] T. Bodineau, M. Lagouge, Large deviations of the empirical currents for a boundary-driven reaction diffusion model, Ann. Appl. Probab. 22 (2012) 2282–2319.
  * [8] C. Boldrighini, A. De Masi, A. Pellegrinotti, Nonequilibrium fluctuations in particle systems modelling reaction-diffusion equation, Stoch. Proc. Appl. 42 (1992) 1–30.
  * [9] G. Carinci, C. Giardinà, C. Giberti, F. Redig, Duality for stochastic models of transport, J. Stat. Phys. 152 (2013) 657–697.
  * [10] M. Cavalcanti, V.D. Cavalcanti, L. Tebou, Stabilization of the wave equation with localized compensating frictional and Kelvin-Voigt dissipating mechanism, Elect. J. Diff. Equ. 83 (2017) 1–18.