# Distributed fixed resources exchanging particles: Phases of an asymmetric
exclusion process connected to two reservoirs

Sourav Pal<EMAIL_ADDRESS>Theory Division, Saha Institute of Nuclear
Physics, A CI of Homi Bhabha National Institute, 1/AF Bidhannagar, Calcutta
700064, West Bengal, India  Parna Roy<EMAIL_ADDRESS>Shahid Matangini
Hazra Government General Degree College for Women, Purba Medinipore 721649,
West Bengal, India  Abhik Basu<EMAIL_ADDRESS><EMAIL_ADDRESS>Theory Division, Saha Institute of Nuclear Physics, A CI of Homi Bhabha
National Institute, 1/AF Bidhannagar, Calcutta 700064, West Bengal, India

###### Abstract

We propose and study a conceptual one-dimensional model to explore how the
combined interplay between fixed resources and particle exchanges between
different parts of an extended system can affect the stationary densities in a
current carrying channel connecting different parts of the system. To this
end, we consider a model composed of a totally asymmetric simple exclusion
process (TASEP) connecting two particle reservoirs without any internal
dynamics but which can directly exchange particles between each other,
ensuring nonvanishing currents in the steady states. The total particle number
in the system that defines the “resources” available, although is kept
constant by the model dynamics, can take any value independent of the model
parameters that define the dynamics of the model. We show how the resulting
phase diagrams of the model are controlled by the parameters, which define the
various dynamical update rules together with the total available resources.
These control parameters can be tuned to make the density on the TASEP lane
globally uniform or piecewise continuous with localized domain walls, and can
also control populations of the two reservoirs. In general, the phase diagrams
are quite different from a TASEP with open boundaries. In the limit of large
amount of resources, the phase diagrams in the plane of the control parameters
become topologically identical to that for an open TASEP together with
delocalization of domain walls.

##  I INTRODUCTION

Far-from-equilibrium phenomena are ubiquitous in nature, which, in contrast to
their equilibrium counterparts, are yet to be fully explored. Studies on
driven diffusive systems are fundamentally important to understand the rich
behaviour of low dimensional collective particle transport in nonequilibrium
systems. A particularly important model is the Totally Asymmetric Simple
Exclusion Process (TASEP), originally proposed to study protein synthesis in
biological cells macdonald . Subsequently, it has been used as a paradigmatic
one-dimensional (1D) open model that shows boundary-induced phase transitions
krug ; krug1 ; krug2 . A TASEP consists of a 1D lattice with $L$ sites having
open boundaries. The dynamics is stochastic in nature, and involves
unidirectional particle hopping, subject to exclusion at all the sites, i.e.,
any site can contain at most one particle at a time. In a TASEP, a particle
enters at the left boundary ($i=1$) of the 1D lattice at a given rate
$\alpha_{T}$, hops unidirectionally to the following sites at rate unity
subject to exclusion until it reaches the last site at the other end ($i=L$),
from which it leaves the system at another specified rate $\beta_{T}$.
Controlled by $\alpha_{T}$ and $\beta_{T}$, the TASEP exhibits three distinct
phases: the steady state densities in bulk can be either less than half, which
is the low-density (LD) phase, or more than half, forming the high-density
(HD) phase, or just half, which is the maximal current (MC) phase of TASEP.
Furthermore, in the $\alpha_{T}-\beta_{T}$ plane there is a coexistence line
separating the LD and HD phases, where nonuniform densities in the form of a
domain wall (DW) can be found. This DW is delocalized, a property that is
attributed to the particle nonconservering dynamics governing the open TASEP
tasep-rev .

A periodic TASEP on a uniform ring necessarily has a constant mean density,
owing to the translational invariance of the system. Macroscopically
nonuniform densities can be expected only when the translational invariance is
broken. This can be achieved by, for instance, the presence of bottlenecks,
point lebo ; niladri-tasep  or extended (either a slower segment or a segment
executing diffusive motion) mustansir ; hinsch ; tirtha-niladri ; parna-anjan
; tirtha-qkpz , or a combination of both sm-atri3 . A TASEP with open
boundaries does not have any particle number conservation, whereas a TASEP on
a ring obviously conserves the total particle number. A third possibility is
an open TASEP with a global constraint on the total particle number, or
equivalently, a closed system of a TASEP connected to a particle reservoir of
finite capacity or a particle storage; see, e.g., Refs. reser1 ; reser2 . More
than one TASEP connected to a reservoir has also been considered; see Ref.
reser3 . A closed system containing a TASEP and a reservoir manifestly breaks
the translation invariance, and hence is a potential candidate for
inhomogeneous steady state densities. These are also expected to be relevant
in related biological processes of protein synthesis in cells reser1  and also
in the context of traffic traffic1 ; see also Ref. limited  for a similar
study. Detailed studies, both numerical MCS and analytical MFT reveal rich
nonuniform steady state density profiles including domain walls in these
models reser1 ; reser2 ; reser3 . The ideas behind a TASEP connected to a
particle storage of finite capacity have been applied to various different
problems, e.g., limited resources in driven diffusive systems lim-driv ,
different biological contexts like mRNA translations and motor protein
dynamics in cells lim-bio1 ; lim-bio2  and traffic problems traffic1 . Notable
experimental studies relevant to these model systems include studies on
spindles in eukaryotic cells lim-exp1 ; lim-exp2 .

A fundamentally important question in nonequilibrium systems is the
consequence of coupling driven with equilibroum processes on the steady
states. In particular, we are interested in how asymmetric exclusion process,
a paradigmatic nonequilibrium process, couples with equilibrium diffusive
exchanges to produce nontrivial steady states. Diffusion is quite common in
cell biology diff-cell . Molecular motors in cells undergo driven motion along
the microtubules, which are quasi-1D channels, or diffusion in the cell
cytoplasm embedding the microtubules, providing a natural example to study the
mutual interplay between driven and equilibrium dynamics. Such interplay
between driven and diffusive dynamics has been explored by studying a variety
of models. These are characterised by specific descriptions of particle
diffusion, how it couples to the driven part (e.g., TASEP) of the dynamics and
the role of any global conservation laws. For instance, the consequences of
the effects of diffusion on the steady states of the filament (TASEP) has been
explored in Refs. klump-PRL ; klump-JSTAT ; klump-PRE ; ciandrini-2019 ;
dauloudet-2019 , by generally modeling the system as a filament (microtubule)
executing TASEP confined in a three-dimensional (3D) particle reservoir
representing the cell cytoplasm, where the motors diffuse around. In an
approach that is complementary to these 3D model studies, a two-lane system
has been considered with one lane being a TASEP and the other executing simple
exclusion process (SEP), a 1D representation of diffusion. This has been
considered in closed hinsch ; parna-anjan , half-open frey-graf-PRL ; bojer-
graf-frey-PRE  and open geometries sm , each revealing the nontrivial steady
states borne by the interplay between TASEP and SEP.

Although significant progress has been made in the research in TASEPs with
finite resources, there are still unresolved questions. In a recent study on a
similar model consisting of a single TASEP lane connected to two reservoirs at
both its ends sourav-1 , it has been shown that in general for a given total
particle number, a conserved quantity, not all the phases of an open TASEP may
be realised by the TASEP in that model. Further, it has been argued in Ref.
sourav-1  that when the particle exchanges between the two reservoirs compete
with the particle hopping in the TASEP, the effects of fluctuations remain
significant even in the thermodynamic limit, in contrast to when the particle
exchange processes dominate over the TASEP current. As a result, the phase
diagrams in the two cases can be quite different from each other. Furthermore,
not only there is a strict global conservation law, there is also a strict
upper limit on the number of particles that can be accommodated, given a set
of the model parameters that define the dynamical updates, in the model of
Ref. sourav-1 , for which the model is completely filled with the TASEP bulk
density reaching unity. A consequence of this is that this model never reduces
to an ordinary TASEP, even for a very large admissible number of total
particles (possible for specific choices of the model parameters). Since the
model in Ref. sourav-1  is effectively a two-constraint model, a pertinent
question is what happens when one of them is lifted. To that end, we construct
and study a conceptual model whose dynamics necessarily conserves the total
particle number, but has no upper limit on the admissible total particle
number for any value of the dynamical model parameters.

In this work, we focus on a TASEP connected to two particle reservoirs, which
can also directly exchange particles between them. The latter process ensures
a finite steady state current through the TASEP. The actual entry to and exit
rates from the TASEP, being dependent on the “supply side” and “receiving
side” instantaneous reservoir populations through some functions of these
populations parametrised by two constants $\alpha$ and $\beta$, are
dynamically derived quantities; see next Section for details. Our model is
constructed in such a way that it (i) strictly conserves the total particle
number (covering the TASEP lane and the two reservoirs), but (ii) can
accommodate any number of particles (including a diverging number of particles
approaching infinity) for any values of the model parameters which define the
dynamical update rules. This stands in stark contrast to a similar model with
a single TASEP channel connected to two reservoirs exchanging particles
sourav-1 , where unlimited accommodation of particles is generally not
possible. This has a strong consequence: in the present study, the phase
diagram of our model smoothly reduces to one that is topologically identical
to that for an open TASEP in the limit of large (formally diverging) number of
total particles for any choices of the model parameters that define the update
rules. This does not hold true for the model in Ref. sourav-1 . Lastly, with
increasing total available particle number or “resources”, a localized DW in
this model gradually delocalizes. We study our model in two different limits:
(i) when the particle exchanges between the reservoirs compete with the
hopping along the TASEP, called weak coupling limit, and (ii) when the former
overwhelms the latter, designated as the strong coupling limit of the model.

The remaining parts of this article are organised as follows. In Section II,
we define our model and give the dynamical update rules. It is then followed
by setting up a mean-field analysis of the steady state densities in the weak
coupling limit given in Section III. Next, in Section IV phase diagrams are
obtained using analytical mean-field theory corroborated by Monte Carlo
simulation. The steady state densities and phase boundaries are derived in
Section V; wherein in Section V.6 the nature of the phase transitions of the
model is discussed. We summarise the findings in Section VIII. Our results on
the strong coupling limit are given in the Appendix.

##  II THE MODEL

Our model consists of a 1D lattice $T$ having $L$ sites, connected to two
particle reservoirs $R_{1}$ and $R_{2}$ at both ends. The reservoirs are
considered to be point reservoirs without any internal structure or dynamics.
The hard-core exclusion principle enforces that each site of $T$ can hold a
maximum of one particle at any given moment. The lattice $T$ executes TASEP
dynamics sourav-1 . Provided the first site of $T$ is empty, a particle from
$R_{1}$ enters $T$ with a rate $\alpha_{\text{eff}}$ that depends on the
instantaneous population of $R_{1}$. The particle then hops clockwise to the
next available site. The hopping rate in the bulk of $T$ is taken to be unity
to set the time scale. Eventually, after reaching the last site of $T$, the
particle jumps into $R_{2}$ with a rate $\beta_{\text{eff}}$ that depends on
the instantaneous population of $R_{2}$. The closed geometry of the model
ensures the overall particle number (including both the TASEP lane and
reservoirs), $N_{0}$, to be conserved. This can be expressed as

$N_{0}=N_{1}+N_{2}+\sum_{i=1}^{L}n_{i},$ (1)

where $N_{1},N_{2}$ are the number of particles in $R_{1},R_{2}$ respectively,
and $n_{i}$’s are binary variables denoting the occupation number of sites
$i=1,2,...,L$ in the TASEP lane. We choose $n_{i}=0$ for vacant sites and
$n_{i}=1$ for occupied ones. The dynamical dependence of effective rates
$\alpha_{\text{eff}}$ and $\beta_{\text{eff}}$ on the instantaneous
populations $N_{1}$ and $N_{2}$ of the reservoirs $R_{1}$ and $R_{2}$ can be
expressed as following:

$\alpha_{\text{eff}}=\alpha f(N_{1});\hskip 14.22636pt\beta_{\text{eff}}=\beta
g(N_{2}),$ (2)

Figure 1: Schematic diagram of the model: Here, $R_{1}$ and $R_{2}$ are two
particle reservoirs, both connected to a TASEP lane $T$ with $L$ sites at its
both ends. Particles, represented by black solid circles, enter the first site
($i=1$) of $T$ from $R_{1}$ with a rate $\alpha_{\text{eff}}$, followed by
hopping to the sites on the right side with a rate of unity, subject to
exclusion. Eventually, particles exit from the last site ($i=L$) of $T$ to go
to $R_{2}$ with a rate of $\beta_{\text{eff}}$. Furthermore, particle
exchanges occur directly between $R_{1}$ and $R_{2}$ with rates $k_{1}$ and
$k_{2}$, modelling diffusion in the system (see text).

where $f$ and $g$ are functions of $N_{1}$ and $N_{2}$ respectively. Inflow of
particles from $R_{1}$ to $T$ increases with the rising population of $R_{1}$,
whereas the outflow of particles from $T$ to $R_{2}$ should be hindered by the
growing population in $R_{2}$. This suggests that $f$ and $g$ should be
monotonically increasing and decreasing functions of $N_{1}$ and $N_{2}$
respectively. We choose the following simple functions which captures these
criteria:

$f(N_{1})=\frac{N_{1}}{N_{0}},\hskip
14.22636ptg(N_{2})=1-\frac{N_{2}}{N_{0}}.$ (3)

Similar functions have been used elsewhere before; see, e.g., Refs. astik-
parna ; arvind1 ; arvind2 ; arvind3 . Notice that the choices (3) are
different from those used in Refs. reser1 ; reser2 ; reser3 ; sourav-1 . The
rate functions in the latter studies have the TASEP size $L$ appearing as a
normalisation constant, whereas in the present study as well as those in Refs.
astik-parna ; arvind1 ; arvind2 ; arvind3 , $N_{0}$ appears as the
normalisation constant. While this may appear as a minor technical matter, our
results show that depending upon these choices, the results can change
significantly.

In our model, the upper limit on each of $N_{1},\,N_{2}$ can be $N_{0}$, the
total particle number, or the “total available resources” in the system.
Further, we do not restrict $N_{0}$ to be finite, i.e., $N_{0}$ can take any
value including approaching infinity. Thus $N_{1}$ and $N_{2}$ can also grow
indefinitely, independent of the other model parameters. This is in stark
contrast to the model studied in Ref. sourav-1 . The rate functions $f$ and
$g$ have the following properties: (i) $f(0)=0$, (ii) $f(\infty)=1$, (iii)
$g(0)=1$, and (iv) $g(\infty)=0$. In general, the rate functions $f$ and $g$
have values between 0 and 1. In addition to particle hopping motion in $T$,
the reservoirs can directly exchange particles between them with certain rates
to maintain a nonzero current in $T$: $R_{1}(R_{2})$ releases particles at
rate $k_{1}(k_{2})$ which $R_{2}(R_{1})$ receives instantaneously. We define a
parameter $\mu=N_{0}/L$ to denote the availability of particles in the system.
Since the total particle number, $N_{0}$, has no upper limit, $\mu$ can be any
positive real number: $0\leq\mu<\infty$. Taken together, our model has five
control parameters: $\alpha$, $\beta$, $\mu$, $k_{1}$, and $k_{2}$, same as
the model studied in Ref. sourav-1 .

##  III MEAN-FIELD THEORY

We now set up the mean-field theory (MFT) for our model by following the same
line of reasoning as in previous studies blythe . MFT involves neglecting the
correlation effects and replacing averages of products by the products of
averages. The ensuing MFT equations take the same form as those in Ref.
sourav-1 , which we revisit here for the convenience of the readers.
Considering the density $\rho_{i}$ as the time-average of occupation number
$n_{i}$ at site $i$, $\rho_{i}\equiv\langle n_{i}\rangle$, we have the
following equations of motion for the bulk sites $1<i<L$ in MFT:

$\displaystyle\frac{d\rho_{i}}{dt}=\rho_{i-1}(1-\rho_{i})-\rho_{i}(1-\rho_{i+1}).$
(4)

For the boundary sites ($i=1$ and $i=L$), we have

$\displaystyle\frac{d\rho_{1}}{dt}=\alpha_{\text{eff}}(1-\rho_{1})-\rho_{1}(1-\rho_{2}),$
(5)
$\displaystyle\frac{d\rho_{L}}{dt}=\rho_{L-1}(1-\rho_{L})-\beta_{\text{eff}}\rho_{L}.$
(6)

Similarly, the time-evolution of reservoir populations $N_{1}$ and $N_{2}$ are
represented by the following equations:

$\displaystyle\frac{dN_{1}}{dt}=k_{2}N_{2}-k_{1}N_{1}-J_{T-in},$ (7)
$\displaystyle\frac{dN_{2}}{dt}=k_{1}N_{1}-k_{2}N_{2}+J_{T-out},$ (8)

where $J_{T-in}=\alpha_{\text{eff}}(1-\rho_{1})$ and
$J_{T-out}=\beta_{\text{eff}}\rho_{L}$ are, respectively, the incoming and
outgoing currents in the TASEP.

Now, introduce a quasi-continuous variable $x=i/L\in[0,1]$ where
$\epsilon=1/L$ is the lattice constant sourav-1 . The density at a position
$x$ at time $t$ is denoted by $\rho(x,t)\equiv\rho_{i}(t)$. As shown in Ref.
sourav-1 , the discrete equation of motion [Eq. (4)] takes a conservation law
form

$\displaystyle\frac{\partial\rho(x,t)}{\partial
t}=-\epsilon\partial_{x}(\rho(x)[1-\rho(x)])+{\cal O}(\epsilon^{2}),$ (9)

supplemented by the boundary conditions
$\rho(0)=\alpha_{\text{eff}},\,\rho(1)=1-\beta_{\text{eff}}$. Eq. (9) has a
conservation law form, reflecting the conserving nature of the TASEP dynamics
in the bulk, and thus allows us to extract a mean-field particle current
$J_{T}$ in the bulk

$J_{T}=\rho(1-\rho),$ (10)

which becomes a constant in the steady state, along with
$J_{T-in}=J_{T-out}=J_{T}$ in the steady states.

In the nonequilibrium steady states, the reservoir populations are constants
over time in MFT. Hence,

$\frac{dN_{1}}{dt}=\frac{dN_{2}}{dt}=0.$ (11)

Substituing Eq. (11) in Eqs. (7) and (8), we get

$k_{2}N_{2}=k_{1}N_{1}+J_{T}.$ (12)

Eq. (12) is a flux balance equation relating the particle current $J_{T}$ with
the reservoir populations $N_{1}$ and $N_{2}$.

As discussed in Ref. sourav-1 , the solution of Eq. (10) gives the steady
state densities

$\rho=\frac{1}{2}\left(1\pm\sqrt{1-4J_{T}}\right):=\rho_{\pm}.$ (13)

Thus, $\rho$ can exhibit spatial uniformity with a low-density (LD) phase
characterised by $\rho=\rho_{-}\leq 1/2$ and a high-density (HD) phase
characterised by $\rho=\rho_{+}>1/2$. Additionally, a maximum current (MC)
phase can arise with $J_{T}=1/4$ when $\rho=1/2$. Another possibility is the
presence of a macroscopic nonuniformity in density in the form of a domain
wall connecting the LD and HD domains, which is the domain wall (DW) phase or
a shock phase.

When solved for $N_{1}$ and $N_{2}$, Eq. (12) along with PNC,
$N_{0}=N_{1}+N_{2}+N_{T}$, where $N_{T}\equiv L\int_{0}^{1}\rho(x)\,dx$ is the
number of particles in the TASEP lane, gives sourav-1

$N_{1}=\frac{k_{2}(N_{0}-N_{T})-J_{T}}{k_{1}+k_{2}},$ (14a)
$N_{2}=\frac{k_{1}(N_{0}-N_{T})+J_{T}}{k_{1}+k_{2}}.$ (14b)

Thus, $N_{1},\,N_{2}$ depend upon $J_{T}$ explicitly, and also implicitly
through the dependence of $N_{0}$ on $J_{T}$.

Similar to Ref. sourav-1 , we consider two distinct cases. We have $J_{T}<1$
($J_{\text{max}}=1/4$), whereas $N_{1}$ and $N_{2}$ rise with $L$, since
$N_{0},\,N_{T}$ scale with $L$ for a given $\mu$. Thus, if $k_{1,2}\sim{\cal
O}(1)$, we can neglect $J_{T}$ safely in Eq. (12) in the thermodynamic limit
(TL). Similarly, we note that in Eqs. (14a) and (14b), where $J_{T}<1$, but
$N_{0}$ and $N_{T}$ should rise with $L$ indefinitely. Thus, if
$k_{1,2}\sim{\cal O}(1)$, we can neglect $J_{T}$ in TL, giving
$k_{1}N_{1}=k_{2}N_{2}$ asymptotically exact, independent of $J_{T}$.
Therefore, in TL, the relative population of the two reservoirs is controlled
by just the ratio $k_{2}/k_{1}$, and independent of $J_{T}$, $\alpha$, and
$\beta$. This effectively eliminates one of the reservoirs, enabling us to
describe the steady state in terms of an effective single reservoir. This is
the diffusion dominated regime, named strong coupling limit below. Another
regime of interest, where diffusion process competes with hopping, is also
discussed. We call this weak coupling limit, which can be achieved by
introducing a mesoscopic scaling $k_{1(2)}=k_{10(20)}/L$ where
$k_{10,20}\sim\mathcal{O}(1)$. In this case, the general solutions for the
reservoir populations $N_{1}$ and $N_{2}$ in the MFT are given by

$N_{1}=\frac{k_{20}(N_{0}-N_{T})-LJ_{T}}{k_{10}+k_{20}},$ (15a)
$N_{2}=\frac{k_{10}(N_{0}-N_{T})+LJ_{T}}{k_{10}+k_{20}},$ (15b)

revealing their explicit $J_{T}$-dependence, which do not vanish in the limit
$L\rightarrow\infty$, since both $N_{0},\,N_{T}$ also scale linearly with $L$
for a fixed $\mu$. This ensures competition between $J_{T}$ and the other
contributions in (15a) and (15b). In the main text of the paper, phase
diagrams and density profiles in weak coupling limit are discussed
elaborately. We discuss the strong coupling limit case in the Appendix.
Although these classifications of weak and strong coupling cases of our model
are identical to those in Ref. sourav-1 , as in both the models these are done
in terms of the scaling of $k_{2},\,k_{2}$ with $L$, nonetheless, we show
below that the results from the present study are significantly different from
those obtained in Ref. sourav-1  due to the fundamentally different forms for
the rate functions $f$ and $g$.

##  IV MEAN-FIELD PHASE DIAGRAMS IN THE WEAK COUPLING CASE

Figure 2: Phase diagrams in the $\alpha-\beta$ plane with different values of
the filling factor $\mu$ in the weak coupling limit of the model, with
exchange rates $k_{10}=0$ and $k_{20}=0.95$. The LD, HD, MC, and DW phases are
denoted respectively by green, gray, yellow, and blue regions with solid black
lines as the phase boundaries, according to MFT; see Eqs. (34), (35), (56),
and (57). The phase boundaries obtained from MCS are represented by discrete
colored points: red (LD-MC), blue (HD-MC), orange (LD-DW), and magenta (HD-
DW). Noticeable deviations between the theoretical and simulated phase
diagrams are observed for intermediate values of $\mu$. Either two or four
phases can be simultaneously observed for small or intermediate values of
$\mu$. Interestingly, for large values of $\mu$ (e.g., $\mu=100$), the phase
diagram of our model exhibits a structure similar to an open TASEP without any
noticeable discrepancies between MFT and MCS results, with LD, HD, and MC
phases become visible and the DW phase region contracts to an inclined line.

Before presenting the phase diagrams of the present model, we briefly revisit
the usual open TASEP model blythe . The phase diagram of an open boundary
TASEP, with $\alpha_{T},\,\beta_{T}<1$ as the entry and exit rates which are
the control parameters, exhibits three distinct phases: an entry-limited low-
density (LD) phase with bulk density $\rho_{\text{LD}}<1/2$, an exit-limited
high-density (HD) phase with bulk density $\rho_{\text{HD}}>1/2$, and a bulk-
limited maximal current (MC) phase with bulk density $\rho_{\text{MC}}=1/2$.
Conditions for the occurrence of these phases are as follows: for the LD
phase, we have $\rho_{\text{LD}}=\alpha_{\text{T}}<\beta_{\text{T}}$ and
$\alpha_{\text{T}}<1/2$; for the HD phase, we have
$1-\rho_{\text{HD}}=\beta_{\text{T}}<\alpha_{\text{T}}$ and
$\beta_{\text{T}}<1/2$, and for the MC phase, we have
$\alpha_{\text{T}},\beta_{\text{T}}\geq 1/2=\rho_{\text{MC}}$. On the
coexistence line between (0,0) and (1/2,1/2), which corresponds to points with
$\alpha_{\text{T}}=\beta_{\text{T}}<1/2$, the bulk density in $T$ exhibits a
piecewise discontinuity, forming a domain wall (DW) that connects the LD and
HD domains together. This domain wall spans the entire lattice instead of
being confined to a specific region and hence is delocalized. This is
attributed to the particle number nonconservation in the model due to the
stochastic entry-exit events. The conditions for transitions between different
phases in an open TASEP can be obtained by equating the currents in the
respective phases. Specifically, the transition from LD to HD phases occurs
when $\alpha_{\text{T}}=\beta_{\text{T}}<1/2$. Likewise, the transitions from
LD to MC and from HD to MC occur when $\alpha_{\text{T}}=1/2<\beta_{\text{T}}$
and $\beta_{\text{T}}=1/2<\alpha_{\text{T}}$, respectively. In the
$\alpha_{T}-\beta_{T}$ plane, the three phases intersect at a common point
$(1/2,1/2)$. The LD-HD transition manifests as a sudden change in the bulk
density, indicating a first-order phase transition. In contrast, the
transitions from either LD or HD to the MC phase involve a continuous
variation in the bulk density, indicating a second-order phase transition.

Fig. 2 gives the phase diagrams of our model in the $\alpha-\beta$ plane in
the weak coupling limit for various $\mu$ together with exchange rates
$k_{10}=0$ and $k_{20}=0.95$, where the different phases are delineated by
comparing them with those found in an open TASEP. The structure of these phase
diagrams in Fig. 2 is quite sensitive to the values of $\mu$. The steady state
bulk density in the LD phase is $\rho_{\text{LD}}=\alpha_{\text{eff}}$, while
in the HD phase, $\rho_{\text{HD}}=1-\beta_{\text{eff}}$. In an open TASEP,
phase transitions between different phases can be determined by equating the
currents corresponding to each phase. Extending that principle to the present
case, the transition from the LD to the HD phase occurs when
$\alpha_{\text{eff}}=\beta_{\text{eff}}<1/2$. Similarly, the transitions
between the LD and MC phases, as well as between the HD and MC phases, are
characterised by the conditions $\alpha_{\text{eff}}=1/2<\beta_{\text{eff}}$
and $\beta_{\text{eff}}=1/2<\alpha_{\text{eff}}$, respectively.

The phase diagram with $\mu=0.5$ [see Fig. 2 (top left)], there are only LD
and DW phases. Indeed, with such a low value of $\mu$, there are fewer
particles and hence the existence of the HD and MC phases in the TASEP lane is
not possible. We find that this particular phase pattern persists up to
$\mu=0.76$ for $k_{10}=0$ and $k_{20}=0.95$. As $\mu$ exceeds this threshold,
the MC and HD phases start to emerge. Consequently, the phase diagrams with
$\mu>0.76$ consist of all four phases, namely LD, HD, MC, and DW phase, see
the phase diagrams corresponding to $\mu=1$ and $\mu=1.5$ in Fig. 2. The phase
boundary between the LD(HD) and MC phases appears as a straight line parallel
to the $\beta$-($\alpha$-)axis in both MFT and Monte Carlo simulation (MCS)
studies. Noticeable deviations, however, between MFT and MCS outcomes are
observed in obtaining these phase boundaries as $\mu$ grows. As $\mu$
increases further, the HD phase expands while the DW region shrinks
accordingly, see the phase diagram with $\mu=1.5$. Surprisingly, when $\mu$ is
raised to a significantly large value, such as $\mu=100$, only three distinct
phases, namely LD, HD, and MC, are visible in the phase space. The DW phase
region now gets contracted into an inclined line extending from (0,0) to
(1/2,1/2). The domain wall spans the whole of the TASEP lane without any
specific or unique position, and hence is a delocalized domain wall (DDW).

In the following sections, we calculate the steady state density profiles in
the different phases and the boundaries separating them in MFT.

##  V STEADY STATE DENSITIES AND PHASE BOUNDARIES

###  V.1  Low-density phase

In the LD phase, steady state bulk density is given by

$\rho_{\text{LD}}=\alpha_{\text{eff}}=\alpha\frac{N_{1}}{N_{0}}<\frac{1}{2}.$
(16)

Substituting $N_{1}$ from Eq. (14a) in Eq. (16) and identifying the steady
state particle current in the LD phase as
$J_{\text{LD}}=\rho_{\text{LD}}(1-\rho_{\text{LD}})$, we obtain a quadratic
equation in $\rho_{\text{LD}}$:

$\rho_{\text{LD}}=\alpha\bigg{[}\frac{k_{2}}{k_{1}+k_{2}}\bigg{(}1-\frac{\rho_{\text{LD}}}{\mu}\bigg{)}-\frac{\rho_{\text{LD}}(1-\rho_{\text{LD}})}{N_{0}(k_{1}+k_{2})}\bigg{]}.$
(17)

Now in the weak coupling limit, particle exchange rates
$k_{1(2)}=k_{10(20)}/L$. Substituting that in Eq. (17) and recalling
$N_{0}=\mu L$, we get

$\rho_{\text{LD}}=\alpha\bigg{[}\frac{k_{20}}{k_{10}+k_{20}}\bigg{(}1-\frac{\rho_{\text{LD}}}{\mu}\bigg{)}-\frac{\rho_{\text{LD}}(1-\rho_{\text{LD}})}{\mu(k_{10}+k_{20})}\bigg{]}.$
(18)

Solving Eq. (18) for $\rho_{\text{LD}}$, we get two solutions:

$\displaystyle\begin{split}\rho_{\text{LD}}^{\pm}&=\bigg{(}\frac{1+k_{20}}{2}+\frac{\mu(k_{10}+k_{20})}{2\alpha}\bigg{)}\\\
&\quad\pm\bigg{[}\bigg{(}\frac{1+k_{20}}{2}+\frac{\mu(k_{10}+k_{20})}{2\alpha}\bigg{)}^{2}-\mu
k_{20}\bigg{]}^{\frac{1}{2}}.\end{split}$ (19)

The condition that $\rho_{\text{LD}}\rightarrow 0$ when $\mu\rightarrow 0$
establishes the following as the physically acceptable solution for the LD
phase density:

$\displaystyle\begin{split}\rho_{\text{LD}}&=\rho_{\text{LD}}^{-}\\\
&=\bigg{(}\frac{1+k_{20}}{2}+\frac{\mu(k_{10}+k_{20})}{2\alpha}\bigg{)}\\\
&\quad-\bigg{[}\bigg{(}\frac{1+k_{20}}{2}+\frac{\mu(k_{10}+k_{20})}{2\alpha}\bigg{)}^{2}-\mu
k_{20}\bigg{]}^{\frac{1}{2}}.\end{split}$ (20)

Clearly, $\rho_{\text{LD}}$ is independent of the exit rate parameter $\beta$.
Eq. (16) provides us with $N_{1}$ in terms of $\rho_{\text{LD}}$, which has
already been expressed in terms of the control parameters $\alpha$, $k_{10}$,
$k_{20}$, and $\mu$ of the model in Eq. (20). The PNC condition, which states
$N_{0}=N_{1}+N_{2}+L\rho_{\text{LD}}$ in the LD phase, can be used to find
$N_{2}$ once $N_{1}$ has been obtained. Following are the expressions for
$N_{1}$ and $N_{2}$:

$\displaystyle N_{1}=L\frac{\mu}{\alpha}\rho_{\text{LD}},$ (21) $\displaystyle
N_{2}=L\bigg{(}\mu-\frac{\mu+\alpha}{\alpha}\rho_{\text{LD}}\bigg{)}.$ (22)

From Eq. (20), it can be seen that

$\rho_{\text{LD}}\xrightarrow{\mu\rightarrow\infty}\frac{\alpha
k_{20}}{k_{10}+k_{20}}.$ (23)

We now investigate two special cases: one with equal particle exchange rates,
i.e., $k_{10}=k_{20}$, and another when $k_{10}=0,k_{20}\neq 0$. For the first
case, $\rho_{\text{LD}}$ goes to $\alpha/2$ as $\mu$ approaches $\infty$, see
Eq. (23). Substituting $\rho_{\text{LD}}=\alpha/2$ in the infinite-$\mu$ limit
in Eqs. (21) and (22), we find $N_{1}=N_{0}/2$ and $N_{2}\approx N_{0}/2$.
This implies, for infinitely large value of $\mu$ with equal particle exchange
rates, majority of particles are accumulated in the reservoirs with almost
equal proportion. For the second case when $k_{10}=0$, $\rho_{\text{LD}}$
approaches $\alpha$ as $\mu$ tends to $\infty$, see Eq. (23). When substituted
in Eqs. (21) and (22), $\rho_{\text{LD}}=\alpha$ gives $N_{1}=N_{0}$ and
$N_{2}\approx 0$ in the infinite-$\mu$ limit. Thus, particles are mostly in
reservoir $R_{1}$ with reservoir $R_{2}$ (and also the TASEP lane $T$)
containing far fewer particles. The phase diagram of the model with
$\mu\rightarrow\infty$ in the $\alpha-\beta$ plane is topologically identical
to that of an open TASEP for all nonzero values of $k_{10},k_{20}$.

We now determine the range of $\mu$ over which the LD phase can exist. For
sufficiently small values of $\mu$, there will certainly be an LD phase, the
density of which will gradually rise with an increase in $\mu$, as shown in
Eq. (20). Eventually, $\rho_{\text{LD}}$ approaches $\alpha
k_{20}/(k_{10}+k_{20})$ as $\mu\rightarrow\infty$, as indicated in (23). Thus,
the LD phase can be found for any value of $\mu$ in our model
$0\leq\mu<\infty$. This can be seen in the phase diagrams of Fig. 2, where LD
phase appears for small values of $\mu$ ($\mu=0.5$) as well as large values
($\mu=100$).

Furthermore, the condition $\rho_{\text{LD}}<1/2$ when imposed on Eq. (20)
gives

$\alpha<\bigg{(}\frac{k_{10}+k_{20}}{2k_{20}-\frac{\frac{1}{2}+k_{20}}{\mu}}\bigg{)}.$
(24)

This implies that, for a given set of values for $k_{10}$, $k_{20}$, and
$\mu$, the LD phase can only be obtained for $\alpha$ values satisfying the
condition mentioned above. This is illustrated in the phase diagrams of Fig.
2. Consider, for example, the case with $k_{10}=0$, $k_{20}=0.95$, and $\mu=1$
(upper right diagram), where the condition (24) indicates that the LD phase
can exist only for $\alpha<2.11$ according to MFT.

Fig. 2 gives the phase diagram of the present model for different values of
$\mu$ with particle exchange rates $k_{10}=0$ and $k_{20}=0.95$. The LD phase
is observed for all the values of $\mu$, including large values such as
$\mu=100$.

###  V.2 High-density phase

The steady state bulk density in HD phase is given by

$\rho_{\text{HD}}=1-\beta_{\text{eff}}=1-\beta\bigg{(}1-\frac{N_{2}}{N_{0}}\bigg{)}>\frac{1}{2}.$
(25)

Plugging $N_{2}$ obtained above in Eq. (14b) into Eq. (25) and using the
expression for steady state current in the HD phase as
$J_{\text{HD}}=\rho_{\text{HD}}(1-\rho_{\text{HD}})$, we obtain a quadratic
equation in $\rho_{\text{HD}}$:

$\rho_{\text{HD}}=1-\beta+\beta\bigg{[}\frac{k_{1}}{k_{1}+k_{2}}\bigg{(}1-\frac{\rho_{\text{HD}}}{\mu}\bigg{)}+\frac{\rho_{\text{HD}}(1-\rho_{\text{HD}})}{N_{0}(k_{1}+k_{2})}\bigg{]}.$
(26)

With $k_{1(2)}=k_{10(20)}/L$ in the weak coupling case and $N_{0}=\mu L$, Eq.
(26) translates into

$\rho_{\text{HD}}=1-\beta+\beta\bigg{[}\frac{k_{10}}{k_{10}+k_{20}}\bigg{(}1-\frac{\rho_{\text{HD}}}{\mu}\bigg{)}+\frac{\rho_{\text{HD}}(1-\rho_{\text{HD}})}{\mu(k_{10}+k_{20})}\bigg{]},$
(27)

solving which we get the following two solutions for $\rho_{\text{HD}}$:

$\displaystyle\begin{split}\rho_{\text{HD}}^{\pm}&=\bigg{(}\frac{1-k_{10}}{2}-\frac{\mu(k_{10}+k_{20})}{2\beta}\bigg{)}\pm\bigg{[}\bigg{(}\frac{1-k_{10}}{2}\\\
&\quad-\frac{\mu(k_{10}+k_{20})}{2\beta}\bigg{)}^{2}+\frac{\mu
k_{10}}{\beta}-\mu\bigg{(}1-\frac{1}{\beta}\bigg{)}k_{20}\bigg{]}^{\frac{1}{2}}.\end{split}$
(28)

This has two solutions. Physically, in the limit when $N_{2}\rightarrow 0$,
any reservoir crowding effect ceases to operate and as a result
$\beta_{\text{eff}}\rightarrow\beta$. Furthermore, in that case if there are
infinite resources, i.e., $\mu\rightarrow\infty$, most of the particles are to
be found in $R_{1}$, giving $N_{1}\approx N_{0}$ and
$\alpha_{\text{eff}}\approx\alpha$. Such a situation can be achieved when
$k_{10}=0$ and $k_{20}\neq 0$ along with $\mu\rightarrow\infty$. We use these
considerations to determine which of the two solutions in (28) is physically
acceptable. Indeed, we find that between the two solutions of HD phase density
in Eq. (28), the one with positive discriminant gives

$\rho_{\text{HD}}\xrightarrow{\mu\rightarrow\infty}\bigg{(}1-\frac{\beta
k_{20}}{k_{10}+k_{20}}\bigg{)}.$ (29)

Thus, when $k_{10}=0$, $\rho_{\text{HD}}$ approaches ($1-\beta$) in the limit
$\mu\rightarrow\infty$. This suggests the following as the HD phase density in
the MFT:

$\displaystyle\begin{split}\rho_{\text{HD}}&=\rho_{\text{HD}}^{+}\\\
&=\bigg{(}\frac{1-k_{10}}{2}-\frac{\mu(k_{10}+k_{20})}{2\beta}\bigg{)}+\bigg{[}\bigg{(}\frac{1-k_{10}}{2}\\\
&\quad-\frac{\mu(k_{10}+k_{20})}{2\beta}\bigg{)}^{2}+\frac{\mu
k_{10}}{\beta}-\mu\bigg{(}1-\frac{1}{\beta}\bigg{)}k_{20}\bigg{]}^{\frac{1}{2}}.\end{split}$
(30)

Unsurprisingly, $\rho_{\text{HD}}$ is independent of $\alpha$. To find the
reservoir populations in the HD phase, we first obtain $N_{2}$ from Eq. (25)
in terms of $\rho_{\text{HD}}$ (expression of which is obtained in Eq. (30) in
terms of model parameters), after which PNC relation
$N_{0}=N_{1}+N_{2}+L\rho_{\text{HD}}$ in the HD phase can be used to find
$N_{1}$. We thus find

$\displaystyle
N_{1}=L\bigg{[}\frac{\mu(1-\rho_{\text{HD}})}{\beta}-\rho_{\text{HD}}\bigg{]},$
(31) $\displaystyle
N_{2}=L\bigg{[}\mu-\frac{\mu(1-\rho_{\text{HD}})}{\beta}\bigg{]}.$ (32)

As considered in our MFT for the LD phase above, we again analyse with two
specific cases, $k_{10}=k_{20}$ and $k_{10}=0,k_{20}\neq 0$, in the HD phase.
For the first case, i.e., when $k_{10}=k_{20}$, Eq. (29) indicates that
$\rho_{\text{HD}}$ approaches ($1-\beta/2$) when $\mu\rightarrow\infty$.
Consequently, substituting $\rho_{\text{HD}}=(1-\beta/2)$ in Eqs. (31) and
(32) when $\mu\rightarrow\infty$ yields $N_{1}\approx N_{0}/2$ and
$N_{2}=N_{0}/2$. For the second case, i.e., when $k_{10}=0,k_{20}\neq 0$, Eq.
(29) gives $\rho_{\text{HD}}$ approaching ($1-\beta$) in the limit of
infinite-$\mu$. Thus, substituting $\rho_{\text{HD}}=(1-\beta)$ in Eqs. (31)
and (32) when $\mu\rightarrow\infty$ results into $N_{1}\approx N_{0}$ and
$N_{2}=0$.

When $\mu$ is reduced, HD phase is less likely to occur due to insufficient
supply of particles. Eventually, below a certain value of $\mu$, HD phase will
disappear. This lower threshold is determined later in Section V.3. As
$\mu\rightarrow\infty$, our model reduces to an open TASEP with
$\rho_{\text{HD}}=1-\beta$(for $k_{10}=0,k_{20}\neq 0$), see Eq. (29). Thus,
there will be no upper limit of $\mu$ for HD phase to occur. Using the fact
$\rho_{\text{HD}}>1/2$, one gets the following condition for HD phase
existence:

$\beta<\bigg{(}\frac{k_{10}+k_{20}}{2k_{20}+\frac{k_{10}-\frac{1}{2}}{\mu}}\bigg{)},$
(33)

which implies that, for specific values of $k_{10}$, $k_{20}$, and $\mu$, HD
phase can occur when the condition mentioned above is satisfied. This is shown
in the phase diagrams of Fig. 2. For instance, consider the case with
$k_{10}=0$, $k_{20}=0.95$, and $\mu=1$ (upper right phase diagram), where (33)
suggests that the HD phase will exist for $\beta<0.68$.

The phase diagrams in Fig. 2 illustrate the occurrence of the HD phase for
different values of the parameter $\mu$, specifically $\mu=1$, $\mu=1.5$, and
$\mu=100$. However, for $\mu=0.5$, the phase diagram includes only the LD and
DW phases.

###  V.3 Maximal current phase

The MC phase is characterised by $J_{\text{max}}=1/4$ which is the maximum
possible steady state particle current, and is obtained when the particle
density $\rho=1/2$ throughout the bulk of $T$ blythe . Thus, in this phase,
TASEP lane is half-filled and the PNC relation reads $N_{0}=N_{1}+N_{2}+L/2$.

Analogous to an open TASEP, conditions for the existence of the MC phase in
the present model are $\alpha_{\text{eff}}>1/2$ and $\beta_{\text{eff}}>1/2$.
We determine the boundaries between the LD and MC phases, and between the HD
and MC phases by setting $\rho_{\text{LD}}=1/2$ and $\rho_{\text{HD}}=1/2$ in
Eqs. (20) and (30) respectively, which yields the following as LD-MC and HD-MC
phase boundaries, respectively:

$\displaystyle\alpha=\frac{\mu(k_{10}+k_{20})}{(2\mu-1)k_{20}-\frac{1}{2}},$
(34) $\displaystyle\beta=\frac{\mu(k_{10}+k_{20})}{k_{10}+2\mu
k_{20}-\frac{1}{2}}.$ (35)

Considering that $\alpha$ and $\beta$ are non-negative real numbers in Eq.
(34) and Eq. (35), we obtain

$\displaystyle\mu>\mu_{1}=\bigg{(}\frac{1}{2}+\frac{1}{4k_{20}}\bigg{)},$ (36)
$\displaystyle\mu>\mu_{2}=\bigg{(}\frac{1}{4k_{20}}-\frac{k_{10}}{2k_{20}}\bigg{)}.$
(37)

Since $\mu_{1}>\mu_{2}$ for any positive value of $k_{10}$ and $k_{20}$,
meeting (36) automatically infers fulfillment of (37). Thus the lower
threshold of $\mu$ for MC phase existence is

$\mu_{l}^{\text{MC}}=\bigg{(}\frac{1}{2}+\frac{1}{4k_{20}}\bigg{)}.$ (38)

This is also the lower threshold for HD phase existence, since if
$\mu<\mu_{l}^{\text{MC}}$, the LD-MC phase boundary in Eq. (34) yields
$\alpha<0$, which is unphysical, thus ruling out the emergence of the HD phase
in the phase diagram. This is shown in Fig. 2 for $\mu=0.5$ (upper left phase
diagram), where the exchange rate $k_{20}=0.95$ yielding the lower threshold
of $\mu$ for the existence of the MC and HD phases to be $0.76$. As seen
before, the present model reduces to an open TASEP in the limit
$\mu\rightarrow\infty$. Thus, MC phase can be present for arbitrarily large
values of $\mu$ as well. Taken together, following is the range within which
MC and HD phases appear:

$\bigg{(}\frac{1}{2}+\frac{1}{4k_{20}}\bigg{)}<\mu<\infty.$ (39)

The phase diagrams of the model obtained by MFT and MCS with exchange rates
$k_{10}=0$ and $k_{20}=0.95$ are shown in Fig. 2. The lower threshold for the
existence of the MC phase is given by the condition (39) as $\mu=0.76$ for
these specific exchange rates, according to MFT. This prediction is confirmed
by MCS simulations, where the MC phase does not exist for values of $\mu$ less
than $0.76$, such as $\mu=0.5$.

###  V.4 Domain wall phase

In the DW phase, the density profile is piecewise continuous forming a DW that
connects the LD and HD regions of a stationary density profile. Analogous to
an open TASEP one has $\rho_{\text{LD}}=1-\rho_{\text{HD}}<1/2$ in the DW
phase. However, unlike in an open TASEP, the precise location of the DW in our
model can be obtained by using PNC. Therefore, we get a localized domain wall
(LDW), which is formed, say at a position $x_{w}$, in the bulk of the lane.
The nonuniform bulk density in the DW phase can be represented in general as

$\rho(x)=\rho_{\text{LD}}+\Theta(x-x_{w})(\rho_{\text{HD}}-\rho_{\text{LD}}),$
(40)

where $\rho_{\text{LD}}$ and $\rho_{\text{HD}}$ are the densities respectively
in LD and HD domain of the DW and $\Theta$ is the Heaviside step function,
defined as $\Theta(x)=1$ for $x>0$ and $\Theta(x)=0$ for $x<0$.

Particle number in $T$ is given by

$N_{T}=\sum_{i=1}^{L}\rho_{i}=L\int_{0}^{1}\rho(x)dx.$ (41)

Plugging $\rho(x)$ from Eq. (40) into Eq. (41), we get after simplifying

$N_{T}=L\bigg{[}\alpha\frac{N_{1}}{N_{0}}(2x_{w}-1)+1-x_{w}\bigg{]}.$ (42)

Substituting $N_{T}$ from Eq. (42) into the PNC relation, and the expression
for the steady state current
$J_{\text{DW}}=\rho_{\text{LD}}(1-\rho_{\text{LD}})$ or
$J_{\text{DW}}=\rho_{\text{HD}}(1-\rho_{\text{HD}})$ in the DW phase into Eq.
(14a), we obtain two equations coupled by $N_{1}/N_{0}$ and $x_{w}$:

$\displaystyle\frac{N_{1}}{N_{0}}\bigg{[}\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}+\alpha(2x_{w}-1)\bigg{]}-x_{w}+1=0,$
(43)
$\displaystyle\begin{split}\frac{N_{1}}{N_{0}}&=\frac{k_{2}}{k_{1}+k_{2}}\bigg{[}1-\frac{1}{\mu}\bigg{\\{}\alpha\frac{N_{1}}{N_{0}}(2x_{w}-1)+1-x_{w}\bigg{\\}}\bigg{]}\\\
&\quad-\frac{1}{N_{0}(k_{1}+k_{2})}\alpha\frac{N_{1}}{N_{0}}\bigg{(}1-\alpha\frac{N_{1}}{N_{0}}\bigg{)}.\end{split}$
(44)

These equations provide two solutions for $N_{1}/N_{0}$ as follows:

$\displaystyle\begin{split}\bigg{(}\frac{N_{1}}{N_{0}}\bigg{)}^{\pm}&=\bigg{(}\frac{1}{2\alpha}+\frac{\mu
k_{10}}{2\alpha^{2}}+\frac{\mu k_{20}}{2\alpha\beta}\bigg{)}\\\
&\quad\pm\bigg{[}\bigg{(}\frac{1}{2\alpha}+\frac{\mu
k_{10}}{2\alpha^{2}}+\frac{\mu k_{20}}{2\alpha\beta}\bigg{)}^{2}-\frac{\mu
k_{20}}{\alpha^{2}}\bigg{]}^{\frac{1}{2}},\end{split}$ (45)

With these two possible solutions in Eq. (45), the density in the LD region of
domain wall becomes

$\displaystyle\begin{split}\rho_{\text{LD}}&=\alpha\bigg{(}\frac{N_{1}}{N_{0}}\bigg{)}^{\pm}\\\
&=\bigg{(}\frac{1}{2}+\frac{\mu k_{10}}{2\alpha}+\frac{\mu
k_{20}}{2\beta}\bigg{)}\\\ &\quad\pm\bigg{[}\bigg{(}\frac{1}{2}+\frac{\mu
k_{10}}{2\alpha}+\frac{\mu k_{20}}{2\beta}\bigg{)}^{2}-\mu
k_{20}\bigg{]}^{\frac{1}{2}}.\end{split}$ (46)

At the LD-DW phase boundary, the density corresponding to the LD phase [Eq.
(20)] must be identical to the acceptable density in the LD part of domain
wall [Eq. (46)]. Considering the solution of Eq. (46) with negative
discriminant, we found that it is identical to the LD phase density given in
Eq. (20) under the following condition:

$\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}+\alpha=0.$ (47)

Eq. (47) is exactly the same equation that defines the boundary between LD and
DW phases, which we will obtain later [Eq. (56)] in another way. Thus the
acceptable solution for $N_{1}/N_{0}$ is

$\displaystyle\begin{split}\frac{N_{1}}{N_{0}}&=\bigg{(}\frac{1}{2\alpha}+\frac{\mu
k_{10}}{2\alpha^{2}}+\frac{\mu k_{20}}{2\alpha\beta}\bigg{)}\\\
&\quad-\bigg{[}\bigg{(}\frac{1}{2\alpha}+\frac{\mu
k_{10}}{2\alpha^{2}}+\frac{\mu k_{20}}{2\alpha\beta}\bigg{)}^{2}-\frac{\mu
k_{20}}{\alpha^{2}}\bigg{]}^{\frac{1}{2}}.\end{split}$ (48)

In the DW phase

$\displaystyle\rho_{\text{LD}}+\rho_{\text{HD}}=1,$ (49)
$\displaystyle\implies$ $\displaystyle\alpha_{\text{eff}}=\beta_{\text{eff}},$
(50) $\displaystyle\implies$
$\displaystyle\alpha\frac{N_{1}}{N_{0}}=\beta\bigg{(}1-\frac{N_{2}}{N_{0}}\bigg{)}.$
(51)

One then gets the expression of $N_{2}$ from that of $N_{1}$:

$\displaystyle\begin{split}\frac{N_{2}}{N_{0}}&=1-\frac{\alpha}{\beta}\frac{N_{1}}{N_{0}}\\\
&=1-\bigg{(}\frac{1}{2\beta}+\frac{\mu k_{10}}{2\alpha\beta}+\frac{\mu
k_{20}}{2\beta^{2}}\bigg{)}\\\
&\quad+\bigg{[}\bigg{(}\frac{1}{2\beta}+\frac{\mu
k_{10}}{2\alpha\beta}+\frac{\mu k_{20}}{2\beta^{2}}\bigg{)}^{2}-\frac{\mu
k_{20}}{\beta^{2}}\bigg{]}^{\frac{1}{2}}.\end{split}$ (52)

Consequently, the DW position ($x_{w}$) is obtained by solving Eq. (43). We
obtain

$\begin{split}x_{w}=\frac{1-\frac{N_{1}}{N_{0}}\bigg{[}\mu\bigg{(}\frac{\alpha}{\beta}-1\bigg{)}+\alpha\bigg{]}}{1-2\alpha\frac{N_{1}}{N_{0}}}.\end{split}$
(53)

Substituting $N_{1}/N_{0}$ from Eq. (48) into Eq. (53), one finds $x_{w}$ in
terms of the control parameters, $\alpha$, $\beta$, $\mu$, $k_{10}$, and
$k_{20}$, of the model. Furthermore, densities corresponding to the LD and HD
regions of the DW phase can be calculated as follows:

$\displaystyle\rho_{\text{LD}}$ $\displaystyle=\alpha\frac{N_{1}}{N_{0}}$
$\displaystyle=\bigg{(}\frac{1}{2}+\frac{\mu k_{10}}{2\alpha}+\frac{\mu
k_{20}}{2\beta}\bigg{)}$
$\displaystyle\quad-\bigg{[}\bigg{(}\frac{1}{2}+\frac{\mu
k_{10}}{2\alpha}+\frac{\mu k_{20}}{2\beta}\bigg{)}^{2}-\mu
k_{20}\bigg{]}^{\frac{1}{2}},$ (54)

with $\rho_{\text{HD}}=1-\rho_{\text{LD}}$. The DW height ($\Delta$) is then

$\displaystyle\Delta$ $\displaystyle=\rho_{\text{HD}}-\rho_{\text{LD}}$
$\displaystyle=\bigg{[}\bigg{(}1+\frac{\mu k_{10}}{\alpha}+\frac{\mu
k_{20}}{\beta}\bigg{)}^{2}-4\mu k_{20}\bigg{]}^{\frac{1}{2}}-\bigg{(}\frac{\mu
k_{10}}{\alpha}+\frac{\mu k_{20}}{\beta}\bigg{)}.$ (55)

Figure 3: Plot of an LDW in the weak coupling limit of the model. System size
is $L=1000$ and $2\times 10^{9}$ Monte Carlo steps are performed to get the
density profiles. Lines (solid and dotted) represent the MFT predictions,
whereas discrete points are the MCS results. For smaller values of $\mu$ we
get an LDW. Two different sets of parameters are chosen: $\alpha=0.2$,
$\beta=0.12$, $\mu=1.3$, $k_{10}=0$, and $k_{20}=0.95$ (red); $\alpha=0.2$,
$\beta=0.12$, $\mu=1.7$, $k_{10}=0$, and $k_{20}=0.95$ (blue). Clearly,
increasing the value of $\mu$ shifts the LDW position to the left [see Eq.
(53)].

To determine the boundaries between the LD as well as the HD and DW phases, we
make use of the fact that the domain wall is situated at the far right and
left ends of the TASEP lane at the point of LD-DW and HD-DW transitions,
respectively. By setting $x_{w}=1$ and $x_{w}=0$ in Eq. (53), we obtain the
following equations that define the LD-DW and HD-DW phase boundaries
respectively:

$\displaystyle\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}+\alpha=0,$ (56)
$\displaystyle\frac{N_{1}}{N_{0}}\bigg{[}\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}-\alpha\bigg{]}+1=0,$
(57)

with $N_{1}/N_{0}$ given in Eq. (48). The DW position $x_{w}$ obtained in Eq.
(53) spans from 0 to 1 with $x_{w}=0$ and $x_{w}=1$ corresponding to the HD-DW
and LD-DW phase boundaries, respectively. Analogous to the open TASEP, the DW
phase in the present model is characterised by
$\alpha_{\text{eff}}=\beta_{\text{eff}}<1/2$. Since
$\alpha_{\text{eff}}=\alpha N_{1}/N_{0}<1/2$, the denominator in Eq. (53) is
always positive. The condition $x_{w}>0$ then is met only if the numerator in
Eq. (53) is positive, which leads to

$\displaystyle\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}>\alpha-\frac{N_{0}}{N_{1}}.$
(58)

Also, the condition $x_{w}<1$ when applied in Eq. (53) yields the following:

$\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}<-\alpha.$ (59)

Within the DW phase, the right-hand side of (58), $\alpha-N_{0}/N_{1}$, is
negative. Moreover, the right-hand side of (59) is clearly negative ($\alpha$
being positive always) indicating the left-hand side, $\mu(1-\alpha/\beta)$,
to be negative. Given $\mu>0$, this implies $(1-\alpha/\beta)$ to be negative
or $\alpha>\beta$ over the DW region. We recast these inequalities (58) and
(59) as follows:

$\displaystyle\mu\bigg{(}\frac{\alpha}{\beta}-1\bigg{)}<\frac{N_{0}}{N_{1}}-\alpha,$
(60) $\displaystyle\mu\bigg{(}\frac{\alpha}{\beta}-1\bigg{)}>\alpha.$ (61)

Taken together, (60) and (61) delineate the range of $\mu$ within which the DW
phase exists:

$\bigg{(}\frac{\alpha}{\frac{\alpha}{\beta}-1}\bigg{)}<\mu<\bigg{(}\frac{\frac{N_{0}}{N_{1}}-\alpha}{\frac{\alpha}{\beta}-1}\bigg{)}.$
(62)

The thresholds of $\mu$ for DW phase existence obtained in (62) are not fixed
but can be varied with the control parameters involved. Only the region in the
$\alpha-\beta$ plane for a given $\mu$ that satisfies both (60) and (61)
simultaneously admits an LDW. Notice that the line $\alpha=\beta$ satisfies
(60) but not (61), meaning the line $\alpha=\beta$ falls outside the DW phase
region in the phase diagram. From (60) and (61), we find that as $\mu$
approaches infinity, both LD-DW and HD-DW boundaries converge to the
$\alpha=\beta$ line, which is below them for any finite $\mu$; see the phase
diagrams in Fig. 2.

See Fig. 3 for stationary densities in the form of LDWs. In MFT, PNC gives a
solution for the position $x_{w}$ of the DW, corresponding to an LDW. However,
our results from the MCS studies reveal that for smaller values of $\mu$, an
LDW is obtained, whereas for larger values of $\mu$, the DW gets delocalized.
We attribute this to fluctuations; see Section VII below for additional
results and discussions.

Figure 4: (Left) Plots of $x_{w}$ vs. $\mu$ in the weak coupling limit, with
parameter values: $\alpha=0.5$, $\beta=0.05$, $k_{10}=0$, and $k_{20}=0.95$
(red); $\alpha=0.5$, $\beta=0.2$, $k_{10}=0$, and $k_{20}=0.95$ (blue). Solid
lines depict MFT results, while discrete points show MCS results. MFT and MCS
agree better for lower parameter values, as seen in the first red set. (Right)
Plots of $\Delta$ vs. $\mu$ in the weak coupling limit using the same
parameter values as in the left figure. These plots demonstrate that the DW
position $x_{w}$ diminishes as $\mu$ increases [see Eq. (53)], whereas the DW
height $\Delta$ remains relatively constant with changes in $\mu$ [see Eq.
(55)].

Figure 5: (Left) Plots of $x_{w}$ vs. $k_{20}$ in the weak coupling limit.
Parameter values are: $\alpha=0.5$, $\beta=0.05$, $\mu=1$, and $k_{10}=0$ for
the main plot (red); $\alpha=0.5$, $\beta=0.1$, $\mu=1$, and $k_{10}=0$ for
the inset plot (blue). (Right) Plots of $\Delta$ vs. $k_{20}$ in the weak
coupling limit. Parameter values are same as the left figure. MFT results [see
Eq. (53) and (55)] are shown as solid lines, while MCS outcomes are depicted
as discrete points.

Our MCS studies reveal that increasing the particle number causes them being
accumulated near the entrance end of $T$, thus positioning the DW towards left
side. This is illustrated in Fig. 4 (left), where $x_{w}$ is observed to
decrease with increasing $\mu$. The DW height $\Delta$, however, remains
nearly unchanged with varying $\mu$, see Fig. 4 (right). In Fig. 5, the
variation of $x_{w}$ and $\Delta$ with $k_{20}$ for a fixed value of $\mu$ is
shown.

###  V.5 Phase boundaries meet at a common point

Within a specific range of $\mu$ where all four phases can be present in the
$\alpha-\beta$ plane, the phase boundaries converge at a common point known as
the (four-phase) multicritical point. This is the point where four phase
boundaries [see Eqs. (34), (35), (56), and (57)] intersect:

$(\alpha_{c},\beta_{c})=\bigg{(}\frac{\mu(k_{10}+k_{20})}{(2\mu-1)k_{20}-\frac{1}{2}},\frac{\mu(k_{10}+k_{20})}{k_{10}+2\mu
k_{20}-\frac{1}{2}}\bigg{)}.$ (63)

Thus $(\alpha_{c},\beta_{c})$ depends explicitly on $k_{10}$, $k_{20}$, and
$\mu$. In the limit $\mu\rightarrow\infty$, the multicritical point approaches

$(\alpha_{c},\beta_{c})_{\mu\rightarrow\infty}=\bigg{(}\frac{k_{10}+k_{20}}{2k_{20}},\frac{k_{10}+k_{20}}{2k_{20}}\bigg{)}.$
(64)

Thus, for $k_{10}=k_{20}$,
$(\alpha_{c},\beta_{c})_{\mu\rightarrow\infty}=(1,1)$ and for
$k_{10}=0,k_{20}\neq 0$ we have
$(\alpha_{c},\beta_{c})_{\mu\rightarrow\infty}=(1/2,1/2)$. Clearly from (63)
this multicritical point exists for

$\mu>\bigg{(}\frac{1}{2}+\frac{1}{4k_{20}}\bigg{)},$ (65)

which is the lower threshold of $\mu$ for MC phase existence; see the
condition (39).

The distance $d$ between the origin (0,0) and the multicritical point
$(\alpha_{c},\beta_{c})$ is calculated as

$d=\sqrt{\bigg{(}\frac{\mu(k_{10}+k_{20})}{(2\mu-1)k_{20}-\frac{1}{2}}\bigg{)}^{2}+\bigg{(}\frac{\mu(k_{10}+k_{20})}{k_{10}+2\mu
k_{20}-\frac{1}{2}}\bigg{)}^{2}}.$ (66)

This implies that the distance $d$ diverges as $\mu$ approaches the threshold
$(1/2+1/4k_{20})$ from above. We get

$d\xrightarrow{\mu\rightarrow\infty}\frac{k_{10}+k_{20}}{\sqrt{2}k_{20}}.$
(67)

Consequently for $k_{10}=k_{20}$ one gets $d\rightarrow\sqrt{2}$, while
$d\rightarrow 1/\sqrt{2}$ when $k_{10}=0,k_{20}\neq 0$.

###  V.6 NATURE OF THE PHASE TRANSITIONS

As shown in Fig. 2, the phase diagrams in weak coupling limit of the present
model are markedly different from that for an open TASEP. This is attributed
to the global particle number conservation and dynamically controlled
effective entry and exit rates for the TASEP channel. Depending on the values
of $\mu$, there are either two or four phases in the $\alpha-\beta$ plane; see
Fig. 2. We now briefly discuss the nature of the phase transitions. As in an
open TASEP, we consider the steady state bulk density ($\rho$) as the order
parameter. We recall that in open TASEP, the transition between LD and HD
phases is marked by a sudden jump in $\rho$, thus implying it a first-order
phase transition. On the other hand, the transitions between either LD or HD
phases and the MC phase involve continuous changes in $\rho$, indicative of
second-order phase transitions. In the present model, there are no phase
transitions between the LD and HD phases for any finite $\mu$. Instead, there
are transitions from the LD to DW and DW to HD phases, both of which are
second order transitions. Similarly, all the other phase transitions, such as
LD-MC and HD-MC transitions in the present model are also second-order in
nature characterised by a continuous change in $\rho$. Finally, in the
asymptotic limit of $\mu\rightarrow\infty$, the LD-DW and DW-HD phase
boundaries merge and overlap with the line $\alpha=\beta$, effectively giving
a transition between the LD and the HD phases. This transition is first order
as in an open TASEP.

Figure 6: Plot of the reservoir population ratio $N_{1}/N_{2}$ versus $\mu$ in
the weak coupling limit of the model for fixed (Left) $k_{10}=k_{20}=0.95$
with good agreement between MFT and MCS visible, (Middle)
$k_{10}=0,k_{20}=0.95$ with discrepancies between MFT and MCS results visible.
(Right) Plot of $N_{1}/N_{2}$ with $k_{10}/k_{20}$ for $\mu=1$ in the weak
coupling limit. In these plots, the entry and exit rates are carefully
selected to ensure that the TASEP lane remains only in one phase (which is
here the LD phase). The values of $\alpha=0.1$ and $\beta=1$ are used. Notice
the mismatches between MFT and MCS results increasing with decreasing
$k_{10}/k_{20}$. MFT results in these plots are obtained from Eqs. (20), (21),
and (22).

In Fig. 6, we show the variation of the reservoir population ratio
($N_{1}/N_{2}$) with respect to either $\mu$ with constant exchange rates or
the ratio of exchange rates ($k_{10}/k_{20}$) with a fixed $\mu$ in the weak
coupling limit in the LD phase, see Eqs. (20), (21), and (22). Discrepancies
between MFT and MCS results are clearly visible in some cases.

##  VI ROLE OF FLUCTUATIONS

Figure 7: Disagreements in the MFT and MCS results on the density profiles in
the weak coupling limit of the model. Solid line represents MFT predictions
(see Sections V.1, V.2, V.3, and V.4) while discrete points show MCS results.
Two sets of MCS results, labeled as MCS1 (red circles) and MCS2 (blue
squares), correspond to time-averaged densities with $2\times 10^{9}$ and
$2\times 10^{10}$ runs, respectively. Notably, fluctuations persist even for
longer runs. The system size is fixed at $L=1000$. (Left) For parameters
$\alpha=1.5$, $\beta=1.2$, $\mu=1$, $k_{10}=0$, and $k_{20}=0.95$, MFT
predicts an LD phase with density $\rho_{\text{LD}}=0.44$, while MCS suggests
an MC phase with density $\rho_{\text{MC}}=0.5$. (Middle) In the case of
$\alpha=2$, $\beta=0.3$, $\mu=1$, $k_{10}=0$, and $k_{20}=0.95$, MFT and MCS
predict DW and HD phases (with density $\rho_{\text{HD}}=0.69$), respectively.
(Right) With parameters $\alpha=1.5$, $\beta=0.25$, $\mu=0.5$, $k_{10}=0$, and
$k_{20}=0.95$, both MFT and MCS predict a DW phase. However, MFT suggests an
LDW with different low and high densities and positions compared to MCS. See
text in Section VI.

A significant deviation, particularly for small and intermediate values of
$\mu$, between our analytical mean-field and simulated Monte Carlo results has
been observed in the weak coupling limit of the model; see Fig. 2, Fig. 3,
Fig. 4, Fig. 5, Fig. 7, and Fig. 8. In contrast, we have observed excellent
agreement between MFT and MCS results in the strong coupling case, see Fig. 9,
Fig. 10, Fig. 11, Fig. 12, and Fig. 14 in the Appendix. This Section explains
the reason behind such discrepancies. The argument is similar to a recent
study on a similar model consisting of a single TASEP lane connected to two
reservoirs at both its ends sourav-1  except for large $\mu$-values where the
present model reduces to an open TASEP model. MFT often produces satisfactory
outcomes as long as the fluctuation effects can be neglected.

In Fig. 7, discrepancies between MFT predictions and MCS results in steady
state density profiles in the weak coupling limit are shown for a set of
control parameters mentioned in the caption. With parameters specified in the
left plot, MFT suggest an LD phase with density $\rho_{\text{LD}}=0.44$, while
MCS suggests an MC phase. In the center plot, MFT predicts an LDW, whereas MCS
yields HD phase with density $\rho_{\text{HD}}=0.69$. In the right plot, both
MFT and MCS suggest occurrence of an LDW but with different densities and
positions. The system size is $L=1000$ and time-average with $2\times
10^{9}$(red circles) and $2\times 10^{10}$(blue squares) number of runs are
performed. No indication of improvement between theoretical and simulative
results in the density profiles are observed. We now argue for these
discrepancies.

We now set up heuristic arguments to explain the fluctuation-induced
disagreements between the MFT and MCS results in the thermodynamic limit by
following and extending the arguments developed in Ref. sourav-1 . Fluctuation
induces deviations in the effective value of reservoir populations from their
respective MFT values. One can thus generalise the relations
$\alpha_{\text{eff}}=\alpha N_{1}/N_{0}$ and
$\beta_{\text{eff}}=\beta(1-N_{2}/N_{0})$ to

$\tilde{\alpha}_{\text{eff}}=\frac{\alpha}{N_{0}}N_{1}^{\text{eff}},\hskip
14.22636pt\tilde{\beta}_{\text{eff}}=\beta\bigg{(}1-\frac{N_{2}^{\text{eff}}}{N_{0}}\bigg{)}$
(68)

with

$N_{1}^{\text{eff}}=N_{1}+\delta N_{1},\hskip
14.22636ptN_{2}^{\text{eff}}=N_{2}+\delta N_{2},$ (69)

where $N_{1}$, $N_{2}$ are the solution for the populations of reservoirs
$R_{1}$, $R_{2}$ in MFT [i.e., solutions in Eqs. (15a) and (15b)], $\delta
N_{1}$, $\delta N_{2}$ are the corresponding deviations due to the
fluctuations from their respective mean-field values, and
$\tilde{\alpha}_{\text{eff}},\tilde{\beta}_{\text{eff}}$ are the effective
entry, exit rates in MCS studies. Since $\delta N_{1}$ and $\delta N_{2}$ are
fluctuations, they can be positive or negative, making
$\tilde{\alpha}_{\text{eff}}-\alpha_{\text{eff}}$ and
$\beta_{\text{eff}}-\tilde{\beta}_{\text{eff}}$ positive or negative
accordingly.

In the weak coupling limit, as shown in Eqs. (15a) and (15b), $N_{1}$ and
$N_{2}$ have dependence on TASEP current $J_{T}$ that scales as $LJ_{T}$.
$J_{T}$ being a bounded quantity ($J_{T_{\text{max}}}=1/4$), fluctuation in
$J_{T}$ is $\mathcal{O}(1)$. Hence the reservoir populations $N_{1},N_{2}$ as
well as their fluctuations $\delta N_{1},\delta N_{2}$ must be
$\mathcal{O}(L)$ quantities. Indeed, the ratio of $N_{1},N_{2}$ depends
explicitly on $J_{T}$. This immediately gives
$\tilde{\alpha}_{\text{eff}}-\alpha_{\text{eff}}\sim{\cal O}(1)$ and
$\beta_{\text{eff}}-\tilde{\beta}_{\text{eff}}\sim{\cal O}(1)$ even in the
thermodynamic limit. Furthermore, the sum of the two reservoir populations
$N_{R}=N_{0}-N_{T}$ should scale with $L$ (for a fixed $\mu$) and hence should
have typical fluctuations of the size ${\cal O}(\sqrt{N_{R}})\sim{\cal
O}(\sqrt{L})\ll$ fluctuations in $N_{1}$ or $N_{2}$ in the limit of large $L$.
Given that $N_{1}+N_{2}=N_{R}$, we must have $\delta N_{1}\sim-\delta N_{2}$
to the leading order in $L$ for sufficiently large $L$. This in turn means
when $N_{1}^{\text{eff}}$ is larger (smaller) than $N_{1}$,
$N_{2}^{\text{eff}}$ is smaller (larger) than $N_{2}$, giving
$\tilde{\alpha}_{\text{eff}}>(<)\alpha_{\text{eff}}$ and
$\tilde{\beta}_{\text{eff}}>(<)\beta_{\text{eff}}$. This suggests that in the
weak coupling case, fluctuations persist in the thermodynamic limit, causing
the transition between LD or HD phases to occur at values of $\alpha$ and
$\beta$ distinct from those predicted by MFT. Contrastingly, in the strong
coupling regime, $N_{1}$ and $N_{2}$ become independent of $J_{T}$ in the
thermodynamic limit and are found to maintain a fixed ratio,
$N_{1}/N_{2}=k_{2}/k_{1}$, see Appendix. The particle exchange dynamics
between reservoirs effectively resemble equilibrium dynamics, leading to
$\delta N_{1}$ and $\delta N_{2}$ scaling as ${\cal O}(\sqrt{L})$. As the MF
expressions for $N_{1}$ and $N_{2}$ scale with $L$, in the thermodynamic
limit, $\tilde{\alpha}_{\text{eff}}\rightarrow\alpha_{\text{eff}}$ and
$\tilde{\beta}_{\text{eff}}\rightarrow\beta_{\text{eff}}$. This convergence
ensures a strong quantitative agreement between MFT and MCS results in the
strong coupling case.

Interestingly, the phase diagrams in both weak and strong coupling limit
reveals that MFT and MCS results are in excellent agreement for large values
of $\mu$, see Fig. 2 ($\mu=100$) and Fig. 9 ($\mu=1000$). This is argued
below. With $N_{0}=\mu L$, Eq. (68) leads to

$\displaystyle\tilde{\alpha}_{\text{eff}}=\frac{\alpha(N_{1}+\delta
N_{1})}{N_{0}}=\alpha_{\text{eff}}+\frac{1}{\mu}\frac{\alpha\delta N_{1}}{L},$
(70)
$\displaystyle\tilde{\beta}_{\text{eff}}=\beta\bigg{(}1-\frac{N_{2}+\delta
N_{2}}{N_{0}}\bigg{)}=\beta_{\text{eff}}-\frac{1}{\mu}\frac{\beta\delta
N_{2}}{L},$ (71)

where entry and exit rates with and without tilde corresponds to MCS and MFT
respectively. As $\delta N_{1},\delta N_{2}$ are $\mathcal{O}(L)$ quantities
in the weak coupling limit, fluctuations in the MCS values of entry and exit
rates (second part of Eqs. (70) and (71)) are clearly $\mathcal{O}(1/\mu)$
quantities that vanishes in the limit $\mu\rightarrow\infty$, thereby giving
$\tilde{\alpha}_{\text{eff}}\rightarrow\alpha_{\text{eff}}$ and
$\tilde{\beta}_{\text{eff}}\rightarrow\beta_{\text{eff}}$ in the thermodynamic
limit.

##  VII Delocalization of the LDW

We have discussed above the conditions to form LDWs in the weak coupling limit
of the model in the mean-field theory. The MFT predicts unique LDW positions
[see Eq. (53)], which can be evaluated by using PNC. For low values of $\mu$,
the MFT agrees qualitatively with the MCS results, as can be seen in Fig. 3.
However, as mentioned above, for large values of $\mu$, the LDW clearly
delocalizes; see Fig. 8 (top); see also Fig. 13 (left). More intriguingly, for
$k_{10}>>k_{20}$, the tendency to delocalize vanishes even with very large
values of $\mu$; see Fig. 8 (bottom). We now explore these systematically
below by using our analysis in Section VI above.

Figure 8: Plots of domain walls in the weak coupling limit of the model with
system size $L=1000$. (Top) Increasing $\mu$ sufficiently (of the order of
$L$), we can get DDW for some choices of the particle exchange rates.
Parameter values are: $\alpha=0.25$, $\beta=0.2499$, $\mu=1000$, $k_{10}=0$,
and $k_{20}=0.95$. (Bottom) Parameter values: $\alpha=10$, $\beta=7$,
$\mu=100$, $k_{10}=0.95$, and $k_{20}=0.01$ (red); $\alpha=23$, $\beta=22$,
$\mu=1000$, $k_{10}=0.95$, and $k_{20}=0.01$ (blue). With these sets of
parameter values where $k_{10}>>k_{20}$, only LDWs are formed even with very
large values of $\mu$; see text in Section VII.

We have in general, $N_{1},N_{2}\sim\mathcal{O}(\mu L)$. This means typical
fluctuations in $N_{1},N_{2}$ $\sim{\cal O}(\sqrt{\mu L})$. Now, for a single
DDW-like domain wall to exist in TASEP lane, domain wall position $x_{w}$ must
fluctuate over a scale comparable to $L$, which gives a typical number
fluctuation in TASEP as ${\cal O}(L)$. Number conservation implies number
fluctuations in the TASEP must be the same as the number fluctuations in
$N_{1},\,N_{2}$. This means that fluctuations in $N_{1}$ and $N_{2}$
$\sim\mathcal{O}(L)$, for a DDW to exist, which immediately gives $\mu\sim L$
for DDW to occur parna-anjan . Our MCS results in Fig. 8 (top) satisfies this.
We however note that Fig. 8 (bottom) is not a DDW, in spite of having a large
$\mu$. This can be explained by considering the Fokker-Planck equation for the
probability $P$ to find the DW at a point $x_{w}$ in the bulk tirtha-qkpz ;
astik-parna ; erwin-tobias , which has the form

$\frac{\partial P}{\partial t}=D\frac{\partial^{2}P}{\partial x^{2}}.$ (72)

Here,
$D\equiv[\alpha_{\text{eff}}(1-\alpha_{\text{eff}})+\beta_{\text{eff}}(1-\beta_{\text{eff}})]$
is the diffusivity. For a DW, $\alpha_{\text{eff}}=\beta_{\text{eff}}$, which
for Fig. 8(bottom) become very small. This means the typical time-scale for
the DW to span the entire TASEP channel due to fluctuations becomes very
large, which effectively localizes the DW. This explains an LDW in Fig.
8(bottom).

##  VIII CONCLUSIONS AND OUTLOOK

In this work, we have studied the generic nonequilibrium steady state density
profiles and phase diagrams of a quasi-1D current-carrying channel $T$
connecting two particle reservoirs $R_{1}$ and $R_{2}$ with unlimited storage
capacity. Particles from the left reservoir $R_{1}$ enter the channel $T$,
followed by a unidirectional (from left to right) hopping throughout $T$
subject to exclusion and eventually reach the right reservoir $R_{2}$. To
model diffusion in our system, the reservoirs are allowed to exchange
particles directly between them. The later process ensures a finite current,
making it essentially a nonequilibrium system to study. In this model, entry
and exit of particles into and from the lattice are dynamically controlled by
the instantaneous population of reservoirs, unlike the open TASEP where these
two rates are constant and uncorrelated. In this model, there is no upper
bound for the total particle number that is held conserved by the dynamics.
This sets our model apart from the particle nonconserving dynamics of open
TASEP. Actual or effective entry and exit rates in our model are controlled by
the parameters $\alpha$ and $\beta$, together with the two functions $f$ and
$g$. We choose simple forms of $f$ and $g$, see Eq. (3), where $f$ and $g$
respectively are the monotonically increasing and decreasing functions of
$N_{1}$ and $N_{2}$ with values between 0 and 1, $N_{1}$ and $N_{2}$ being the
populations of reservoirs $R_{1}$ and $R_{2}$. Taken together, the dynamics of
the system is influenced by five control parameters – one entry and exit rate
($\alpha$ and $\beta$), two particle exchange rates ($k_{1}$ and $k_{2}$), and
a filling fraction $\mu$. The behaviour of the system is particularly
interesting when particle hopping through the lattice sites and direct
particle exchange between the reservoirs compete. This corresponds to the weak
coupling limit, captured by scaling $k_{1}$ and $k_{2}$ with
$\mathcal{O}(1/L)$, where $L$ represents the lattice size. We have also
studied the strong coupling limit, where direct particle exchange between
reservoirs dominates particle current through TASEP lane, i.e.,
$k_{1},k_{2}\sim\mathcal{O}(1)$.

We have analysed our model at the mean-field level and corroborated the MFT
analytical results with Monte Carlo simulations. Depending on the value of
$\mu$, generally either two or four phases appear simultaneously in the
$\alpha-\beta$ phase plane, with continuous transitions between them. Above a
lower threshold of $\mu$, all these phases meet at a single point.
Interestingly, in the limit of infinite $\mu$, the phase diagram of our model
mimics that of an open TASEP. Indeed, when $\mu\rightarrow\infty$, the phase
space of the present model includes three phases – LD, HD, and MC — with a DW
region in the form of an inclined line connecting LD phase to the HD phase.
Initially manifested as localized domain walls (LDWs) for smaller $\mu$
values, these are now delocalized domain walls (DDW) in this limit. While PNC
together with MFT still predicts an LDW, large fluctuations around the mean-
field position of the DW effectively delocalizes the domain wall, giving a
DDW.

It is useful to compare our results with a recent work on a similar model with
one TASEP lane connected to two reservoirs directly exchanging particles that
uses different forms for the functions $f$ and $g$, Ref. sourav-1 . First of
all, the model in Ref. sourav-1  has a particle-hole symmetry for equal
exchange rates between the reservoirs in both weak and strong coupling limits.
In contrast, in the present model, there is none. As mentioned in the
Introduction Section, in the model of Ref. sourav-1 , the maximum particle
content is limited by the exchange rates between the reservoirs: Consequently,
$\mu$ in the model of Ref. sourav-1  can have a maximum value
$\mu_{\text{max}}=(2+k_{20}/k_{10})$. In contrast, in the model of the present
study, it is independent or free of any other model parameters. Thus the
present model has less constraints. In addition, the model in Ref. sourav-1 ,
admits particle-hole symmetry in the special case of equal exchange rates,
i.e., $k_{1}=k_{2}$, unlike for the present model that has no particle-hole
symmetry for any choice of the model parameters. In contrast to the model in
Ref. sourav-1 , for any finite values of the exchange rates, the DW in the
present model can delocalize if sufficiently large number of particles are
available. This leads to the significant differences in the phase diagrams in
the two models for large “resources” or large number of particles, i.e., for
large values of $\mu$. Unsurprisingly, for fewer particles, i.e., for small
$\mu$, the differences in the phase diagrams are small. On the whole, our
studies here together with Ref. sourav-1  reveal how the presence or absence
of constraints on the dynamics of the models can strongly influence the
stationary densities and the phase diagrams.

We have used analytical MFT, supplemented by extensive studies to obtain the
stationary densities and the phase diagrams. In the weak coupling limit, we
find significant disagreements between the MFT and MCS results for
intermediate values of $\mu$. Such disagreements disappear in the strong
coupling limit of the model. We further find that the domain wall can
delocalize for some choices of $k_{10},k_{20}$, but does not delocalize for
other choices even for sufficiently larger values of $\mu$. We are able to
explain these results in terms of a set of heuristic arguments that show the
crucial role of fluctuations that survive in the thermodynamic limit in the
weak coupling limit, but not in the strong coupling case.

Our studies can be extended in various ways. One can consider several TASEP
lanes connecting the two reservoirs. More interestingly, one can consider a
“mixed model”, where one of the reservoirs has the structure of the reservoirs
in Ref. sourav-1 , and the other has a form similar to the present study. This
can lead to possible unexpected competitions between the reservoir structures
and the resulting phase diagrams. We hope our studies here will stimulate
further research along these lines.

##  IX Acknowledgement

A.B. thanks SERB (DST), India for partial financial support through the CRG
scheme [file: CRG/2021/001875].

##  Appendix A STRONG COUPLING LIMIT

In this Section, we focus on the strong coupling limit of the model, and study
the phase diagrams and density profiles. The strong coupling limit is
characterised by the particle exchange rates independent of system size $L$,
i.e., $k_{1},k_{2}\sim\mathcal{O}(1)$. This means the direct exchange of
particles between reservoirs dominates over the particle transport via hopping
through the TASEP. Therefore, in this limit, the reservoir populations $N_{1}$
and $N_{2}$ maintain a fixed ratio for a given set of particle exchange rates.
This considerably simplifies the corresponding MFT calculations, which agrees
quantitatively with the MCS results, contrasting the weak coupling limit.

###  A.1 Mean-field phase diagrams and steady state densities

Figure 9: Phase diagrams in the $\alpha-\beta$ plane in the strong coupling
limit of the model for different values of the filling factor $\mu$, with
$k_{1}=0$ and $k_{2}=0.95$. MFT predicts up to four distinct phases (LD, HD,
MC, and DW), represented by the green, gray, yellow, and blue regions,
respectively, with black solid lines separating them, see Eqs. (87), (88),
(97), and (98) for MFT phase boundaries. The phase boundaries obtained from
MCS are depicted by colored points: red (LD-MC), blue (HD-MC), orange (LD-DW),
and magenta (HD-DW). The MFT and MCS results exhibit excellent match. In the
limit of large $\mu$ (say for $\mu=1000$), the phase diagram of the present
model is identical to the familiar phase diagram of an open TASEP.

Phase diagrams in the strong coupling limit are shown in Fig. 9 for a set of
representative values of $\mu$ with $k_{1}=0$ and $k_{2}=0.95$. As explained
in Section III above, in the strong coupling limit

$\frac{N_{1}}{N_{2}}=\frac{k_{2}}{k_{1}}$ (73)

in the thermodynamic limit ($L\rightarrow\infty$). To calculate the steady
state densities in the different phases in the strong coupling limit, we
follow the logic and calculational scheme outlined above for the weak coupling
limit of the model.

####  A.1.1  Low-density phase

We start by using Eq. (17). Last part in the right-hand side of Eq. (17) is
$\mathcal{O}(1/L)$ and hence can be neglected in the thermodynamic limit
($L\rightarrow\infty$). In that limit, the quadratic equation in (17) becomes
linear in $\rho_{\text{LD}}$, solving which we obtain

$\rho_{\text{LD}}\approx\frac{\alpha
k_{2}}{k_{1}+\bigg{(}1+\frac{\alpha}{\mu}\bigg{)}k_{2}}.$ (74)

Clearly, $\rho_{\text{LD}}$ is independent of $\beta$. Considering Eq. (74) we
find the following:

$\displaystyle\rho_{\text{LD}}\xrightarrow{\mu\rightarrow 0}0,$ (75)
$\displaystyle\rho_{\text{LD}}\xrightarrow{\mu\rightarrow\infty}\frac{\alpha
k_{2}}{k_{1}+k_{2}}.$ (76)

Thus, for $\mu\rightarrow\infty$, we have
$\rho_{\text{LD}}\rightarrow\alpha/2$ when $k_{1}=k_{2}$ and
$\rho_{\text{LD}}\rightarrow\alpha$ when $k_{1}=0$.

The population of reservoirs, $N_{1}$ and $N_{2}$, can be obtained in the
strong coupling limit likewise we did in the weak coupling limit, i.e., using
$\rho_{\text{LD}}=\alpha N_{1}/N_{0}$ together with the PNC relation that
relates $N_{2}$ with $N_{1}$. We find:

$\displaystyle N_{1}=\frac{Lk_{2}\mu}{k_{1}+(1+\frac{\alpha}{\mu})k_{2}},$
(77) $\displaystyle
N_{2}=\frac{Lk_{1}\mu}{k_{1}+(1+\frac{\alpha}{\mu})k_{2}}.$ (78)

As previously stated, the relative population of the reservoirs is determined
solely by the ratio of particle exchange rates in the strong coupling limit:
$N_{1}/N_{2}=k_{2}/k_{1}$. This effectively eliminates one of the reservoirs
in the strong coupling limit.

Similar to the weak coupling limit, LD phase can be present for any (positive)
value of $\mu$ in the strong coupling limit. This is demonstrated in Fig. 9,
where the LD phase can be found for small $\mu$ values as well as large $\mu$
values. Since $\rho_{\text{LD}}<1/2$, one gets the following as a condition
for LD phase existence according to MFT:

$\alpha<\frac{k_{1}+k_{2}}{\bigg{(}2-\frac{1}{\mu}\bigg{)}k_{2}}.$ (79)

This is illustrated in Fig. 9. For instance, when $k_{1}=0$, $k_{2}=0.95$, and
$\mu=1$ (upper right phase diagram), the condition (79) infers an LD phase for
$\alpha<1$. Similarly when $\mu$ is very large, say $\mu=1000$ (lower right
phase diagram), LD phase is limited up to $\alpha=0.5$ along the
$\alpha$-axis.

####  A.1.2  High-density phase

We agin proceed with the reasoning presented in the weak coupling limit.
Neglecting the last term of order $\mathcal{O}(1/L)$ in Eq. (26) in the
thermodynamic limit, the HD phase density $\rho_{\text{HD}}$ in the strong
coupling limit can be obtained approximately as:

$\rho_{\text{HD}}\approx\frac{k_{1}+(1-\beta)k_{2}}{\bigg{(}1+\frac{\beta}{\mu}\bigg{)}k_{1}+k_{2}}.$
(80)

As expected, $\rho_{\text{HD}}$ is independent of $\alpha$. It can be shown
from Eq. (80) that

$\rho_{\text{HD}}\xrightarrow{\mu\rightarrow\infty}\bigg{(}1-\frac{\beta
k_{2}}{k_{1}+k_{2}}\bigg{)},$ (81)

which for the two cases $k_{1}=k_{2}$ and $k_{1}=0$, yields
$\rho_{\text{HD}}\rightarrow(1-\beta/2)$ and
$\rho_{\text{HD}}\rightarrow(1-\beta)$, respectively.

The reservoir populations in the HD phase can be obtained as

$\displaystyle
N_{1}=\frac{Lk_{2}(\mu-1+\beta)}{\bigg{(}1+\frac{\beta}{\mu}\bigg{)}k_{1}+k_{2}},$
(82) $\displaystyle
N_{2}=\frac{Lk_{1}(\mu-1+\beta)}{\bigg{(}1+\frac{\beta}{\mu}\bigg{)}k_{1}+k_{2}}.$
(83)

Again, the ratio of reservoir populations remains constant:
$N_{1}/N_{2}=k_{2}/k_{1}$.

The lower threshold of $\mu$ for HD phse existence is determined in Section
A.1.3. Similar to the weak coupling limit case, there will be no finite upper
threshold of $\mu$ for HD phase occurence. In the phase diagrams of Fig. 9,
the occurrence of HD phase is displayed. When the value of $\mu$ is
sufficiently small, such as $\mu=0.3$, HD phase is absent due to an inadequate
number of particles being available. Since $\rho_{\text{HD}}>1/2$, we get the
following condition for HD phase:

$\beta<\frac{k_{1}+k_{2}}{2k_{2}+\frac{k_{1}}{\mu}},$ (84)

which implies that, for a fixed value of $k_{1}$, $k_{2}$, and $\mu$ HD phase
exists if $\beta$ satisfies the above condition. This can be readily seen in
the phase diagrams Fig. 9 where we have taken $k_{1}=0,k_{2}=0.95$. In this
case, HD phase lies under the line $\beta=1/2$ for any value of $\mu$.

####  A.1.3  Maximal current phase

With a uniform hopping rate of 1, the steady state bulk density
$\rho_{\text{MC}}=1/2$ in the MC phase. The PNC relation reads
$N_{0}=N_{1}+N_{2}+L/2$ in this phase. The maximum current $J_{\text{max}}$
associated with the MC phase is given by $J_{\text{max}}=\rho(1-\rho)=1/4$.
Analogous to the open TASEP, following are the boundary conditions in MC
phase:

$\displaystyle\alpha_{\text{eff}}=\alpha\frac{N_{1}}{N_{0}}>\frac{1}{2},$ (85)
$\displaystyle\beta_{\text{eff}}=\beta\bigg{(}1-\frac{N_{2}}{N_{0}}\bigg{)}>\frac{1}{2}.$
(86)

The boundaries between the LD and MC phases, and between the HD and MC phases
can be determined by substituting the given values of $\rho_{\text{LD}}$ and
$\rho_{\text{HD}}=1/2$ into equations (74) and (80), respectively. We thus
obtain following equations as the LD-MC and HD-MC boundaries respectly:

$\displaystyle\alpha=\frac{\mu(k_{1}+k_{2})}{(2\mu-1)k_{2}},$ (87)
$\displaystyle\beta=\frac{\mu(k_{1}+k_{2})}{k_{1}+2\mu k_{2}}.$ (88)

Considering the non-negativity of $\alpha$ in Eq. (87), one finds the MC phase
to exist for any $\mu>1/2$ in the strong coupling limit. Nevertheless, $\beta$
in Eq. (88) is inherently positive. This infers no upper threshold of $\mu$
for MC phase to exist. Following is the range:

$\frac{1}{2}<\mu<\infty.$ (89)

Following the same line of argument as given in the weak coupling limit, (89)
is also the range of $\mu$ within which HD phase appears in the phase diagram,
see Fig. 9.

The existence of MC phase is illustrated in Fig. 9, where it does not appear
for $\mu<1/2$(, such as $\mu=0.3$) .

####  A.1.4  Domain wall phase

In addition to the above phases with uniform density, one can also expect the
emergence of a domain wall (DW) phase with non-uniform density, where the LD
and HD phases are connected by a domain wall in the bulk of the TASEP lane.
The precise location of this DW can be obtained by PNC, inferring it as an
LDW. At this point, we determine the position $x_{w}$ and height $\Delta$ of
the LDW in a similar manner as done for the weak coupling limit case. By
neglecting terms of order ${\cal O}(1/L)$ in the thermodynamic limit, we can
solve the coupled equations (43) and (44) to obtain expressions for
$N_{1}/N_{0}$ and $x_{w}$. The solutions in the strong coupling limit are as
follows:

$\displaystyle\frac{N_{1}}{N_{0}}=\frac{k_{2}}{k_{1}+\frac{\alpha}{\beta}k_{2}},$
(90) $\displaystyle
x_{w}=\frac{1-\bigg{(}\frac{k_{2}}{k_{1}+\frac{\alpha}{\beta}k_{2}}\bigg{)}\bigg{[}\mu\bigg{(}\frac{\alpha}{\beta}-1\bigg{)}+\alpha\bigg{]}}{1-2\alpha\frac{k_{2}}{k_{1}+\frac{\alpha}{\beta}k_{2}}}.$
(91)

Figure 10: Steady state density profile in the DW phase in the strong coupling
limit of the model. Model parameters are: $\mu=3/2$, $k_{1}=k_{2}=0.7$,
$\alpha=0.2$ and $\beta=0.1$. MCS results agree well with the MFT predictions
given by Eqs. (91), (93), and (94).

The reservoir populations, $N_{1}$ and $N_{2}$, are related by the equation
$\alpha_{\text{eff}}=\beta_{\text{eff}}$ or $\alpha
N_{1}/N_{0}=\beta(1-N_{2}/N_{0})$ in the DW phase, which yields the expression
for $N_{2}$:

$\frac{N_{2}}{N_{0}}=1-\frac{\alpha}{\beta}\frac{N_{1}}{N_{0}}=1-\frac{k_{2}}{\frac{\beta}{\alpha}k_{1}+k_{2}}.$
(92)

The steady state low and high densities in the DW phase are

$\displaystyle\rho_{\text{LD}}=\alpha\frac{N_{1}}{N_{0}}=\frac{\alpha
k_{2}}{k_{1}+\frac{\alpha}{\beta}k_{2}},$ (93)
$\displaystyle\rho_{\text{HD}}=1-\rho_{\text{LD}}=1-\frac{\alpha
k_{2}}{k_{1}+\frac{\alpha}{\beta}k_{2}}.$ (94)

Then the domain wall height $\Delta$ is given by

$\Delta=\rho_{\text{HD}}-\rho_{\text{LD}}=1-\frac{2\alpha\beta k_{2}}{\beta
k_{1}+\alpha k_{2}}.$ (95)

The expression of $\Delta$ in the strong coupling limit, as given in Eq. (95),
depends only on parameters $\alpha$, $\beta$, $k_{1}$, and $k_{2}$. Notably,
it is independent of the filling factor $\mu$.

In Fig. 11, the graph illustrates a linear decrease in the position of the
domain wall (DW) as $\mu$ increases, while the height of the DW remains
constant throughout the range of $\mu$ values. On the other hand, in Fig. 12,
both the position and height of the DW exhibit no significant variation as
$k_{2}$ changes.

Following the same arguments given in weak coupling limit of the model, the
range of $\mu$ over which DW phase exists can be identified in strong coupling
limit. Using the fact that DW position, $x_{w}$, should be between 0 and 1, we
find this range as below:

$\bigg{(}\frac{\alpha}{\frac{\alpha}{\beta}-1}\bigg{)}<\mu<\bigg{(}\frac{\frac{k_{1}}{k_{2}}+\frac{\alpha}{\beta}-\alpha}{\frac{\alpha}{\beta}-1}\bigg{)}.$
(96)

Furthermore, the boundaries between the LD and DW, and between the HD and DW
phases can be determined by substituting $x_{w}=1$ and $x_{w}=0$ respectively
into Eq. (91). The resulting equations for LD-DW and HD-DW phase boundaries
are respectively:

$\displaystyle\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}+\alpha=0,$ (97)
$\displaystyle\bigg{(}\frac{k_{2}}{k_{1}+\frac{\alpha}{\beta}k_{2}}\bigg{)}\bigg{[}\mu\bigg{(}1-\frac{\alpha}{\beta}\bigg{)}-\alpha\bigg{]}+1=0.$
(98)

Figure 11: (Left) Plots of DW position ($x_{w}$) vs. the filling factor
($\mu$) in the strong coupling limit of the model. The exchange rates are set
to $k_{1}=0$ and $k_{2}=0.95$, while two different sets of entry and exit rate
parameters are considered: $\alpha=0.5$ with $\beta=0.05$ (red) and
$\alpha=0.5$ with $\beta=0.1$ (blue). Solid lines represent the predictions
from MFT [see Eq. (91)], while the discrete points correspond to the results
obtained through MCS. Within the given range of $\mu$, the $x_{w}$ linearly
decreases with increasing $\mu$. (Right) Plots of DW height ($\Delta$) as a
function of the filling factor ($\mu$) in the strong coupling limit of the
model. The parameter values for this plot are the same as those used in the
left plot, with matching colors. Notably, $\Delta$ remains unchanged as $\mu$
varies within the examined range, see Eq. (95). Both MFT and MCS results are
in good agreement.

Figure 12: (Left) Plots of DW position ($x_{w}$) vs. $k_{2}$ in the strong
coupling limit of the model when $k_{1}=0$. The red plot corresponds to
$\alpha=0.5$ and $\beta=0.1$, while the blue plot represents $\alpha=0.5$ and
$\beta=0.2$. The position of DW does not vary with $k_{2}$, see Eq. (91).
(Right) Plots of DW height ($\Delta$) vs. $k_{2}$ in the strong coupling limit
with $k_{1}=0$. The red and blue plots correspond to the same parameter values
as in the left graph. MCS outcomes satisfy MFT predictions, see Eq. (95).

###  A.2 Delocalization of LDW

Similar to the weak coupling limit, DWs in the strong coupling limit also
gradually get delocalized with increasing $\mu$. This is clearly shown in Fig.
13 where one can see that the domain walls get completely delocalized with
$\mu\sim L$ — a fact already argued in Section VII of the main paper.

Figure 13: Gradual delocalization of DW with increase in $\mu$. (Left) This
corresponds to the weak coupling limit. Parameter values are: $\alpha=0.25$,
$\beta=0.24$, $\mu=10$, $k_{10}=0$, and $k_{20}=0.95$ (red circles);
$\alpha=0.25$, $\beta=0.2488$, $\mu=100$, $k_{10}=0$, and $k_{20}=0.95$ (blue
squares); and $\alpha=0.25$, $\beta=0.24985$, $\mu=1000$, $k_{10}=0$, and
$k_{20}=0.95$ (orange triangles). (Right) This corresponds to the strong
coupling limit. Parameter values are: $\alpha=0.25$, $\beta=0.24985$,
$\mu=10$, $k_{10}=0$, and $k_{20}=0.95$ (red circles); $\alpha=0.25$,
$\beta=0.24985$, $\mu=100$, $k_{10}=0$, and $k_{20}=0.95$ (blue squares); and
$\alpha=0.25$, $\beta=0.24985$, $\mu=1000$, $k_{10}=0$, and $k_{20}=0.95$
(orange triangles). For both weak and strong coupling limits, system size
$L=1000$ and averages over $2\times 10^{9}$ times steps are done to get the
MCS density profiles. Clearly, the LDW for smaller values of $\mu$ gets
delocalized as $\mu\sim L$.

####  A.2.1 Phase boundaries meet at a common point

For specific values of $\mu$ when all four phases (LD, HD, MC, and DW) exist,
they intersect in the $\alpha-\beta$ plane at a common single point. This
point, which is named (four-phase) multicritical point, can be obtaind as the
meeting point of all four phase boundaries in Eqs. (87), (88), (97), and (98):

$(\alpha_{c},\beta_{c})=\bigg{(}\frac{\mu(k_{1}+k_{2})}{(2\mu-1)k_{2}},\frac{\mu(k_{1}+k_{2})}{k_{1}+2\mu
k_{2}}\bigg{)}.$ (99)

Depending explicitly on $k_{1}$, $k_{2}$, and $\mu$ this unique point exists
for any value of $\mu$ greater than 1/2, above which all four phases appear in
the phase diagram, see Fig. 9. The distance between the origin (0,0) and the
multicritical point ($\alpha_{c},\beta_{c}$) can be calculated using the
formula:

$d=\sqrt{\bigg{(}\frac{\mu(k_{1}+k_{2})}{(2\mu-1)k_{2}}\bigg{)}^{2}+\bigg{(}\frac{\mu(k_{1}+k_{2})}{k_{1}+2\mu
k_{2}}\bigg{)}^{2}}.$ (100)

Figure 14: (Left) Plot showing dependence of the reservoir population ratio
$N_{1}/N_{2}$ on $\mu$ in the strong coupling limit of the model. (Right) Plot
presenting the variation of the reservoir population ratio $N_{1}/N_{2}$ with
$k_{1}/k_{2}$ for a fixed filling factor $\mu=1$ in the strong coupling limit.
In these plots, the entry and exit rates are carefully selected to ensure that
the TASEP lane remains in the LD phase. The values of $\alpha=0.1$ and
$\beta=1$ are used. The results obtained by MFT and MCS are in good agreement.

##  Appendix B Simulation algorithm

In this study, the mean-field predicted densities and phase diagrams of the
model were validated through Monte Carlo simulations using random updates. The
simulation rules are as follows: (i) If the first site ($i=1$) of the TASEP
lane ($T$) is empty, a particle from reservoir $R_{1}$ enters $T$ with a rate
$\alpha_{\text{eff}}$; (ii) Particles in $T$ can hop with rate 1 to the next
site in the bulk of $T$ provided that site is empty; (iii) Upon reaching the
last site ($i=L$) of $T$, a particle exits with a rate $\beta_{\text{eff}}$
into reservoir $R_{2}$; (iv) Concurrently, particles in the reservoirs
($R_{1}$ and $R_{2}$) can jump directly to the other reservoir with rates
$k_{1}$ from $R_{1}$ to $R_{2}$ and $k_{2}$ from $R_{2}$ to $R_{1}$,The above
dynamical update rules were implemented with random sequential updates. After
sufficient number of iterations (typically $10\%$ to $30\%$ of the total
number of time steps) were performed to reach a steady state, time-averaged
density profiles were obtained.

## References

  * (1) C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, “Kinetics of biopolymerization on nucleic acid templates”, Biopolymers 6, 1 (1968).
  * (2) J. Krug, “Boundary-induced phase transitions in driven diffusive systems”, Phys. Rev. Lett. 67, 1882 (1991).
  * (3) B. Derrida, S. A. Janowsky, J. L. Lebowitz, and E. R. Speer, “Exact solution of the totally asymmetric simple exclusion process: Shock profiles”, J Stat Phys , 813–842.
  * (4) B. Derrida and M. R Evans, “Exact correlation functions in an asymmetric exclusion model with open boundaries”, Journal de Physique I 3, 311–322 (1993).
  * (5) T. Chou, K. Mallick, and R. K. P. Zia, “Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport”, Rep. Prog. Phys. 74, 116601 (2011).
  * (6) S. A. Janowsky and J. L. Lebowitz, “Finite-size effects and shock fluctuations in the asymmetric simple-exclusion process”, Phys. Rev. A 45, 618 (1992).
  * (7) N. Sarkar and A. Basu, “Nonequilibrium steady states in asymmetric exclusion processes on a ring with bottlenecks”, Phys. Rev. E 90, 022109 (2014).
  * (8) G. Tripathy and M. Barma, “Driven lattice gases with quenched disorder: Exact results and different macroscopic regimes”, Phys. Rev. E 58, 1911 (1998).
  * (9) H. Hinsch and E. Frey, “Bulk-driven nonequilibrium phase transitions in a mesoscopic ring”, Phys. Rev. Lett. 97, 095701 (2006).
  * (10) T. Banerjee, N. Sarkar, and A. Basu, “Generic nonequilibrium steady states in an exclusion process on an inhomogeneous ring”, J. Stat. Mech. J. Stat. Mech. P01024 (2015).
  * (11) P. Roy, A. K. Chandra, and A. Basu, “Pinned or moving: States of a single shock in a ring”, Phys. Rev. E 102, 012105 (2020).
  * (12) T. Banerjee and A. Basu, “Smooth or shock: Universality in closed inhomogeneous driven single file motions”, Phys. Rev. Research 2, 013025 (2020).
  * (13) A. Goswami, R. Chatterjee, and S. Mukherjee, “Defect versus defect: Stationary states of single file marching in periodic landscapes with road blocks”, arXiv: 2402.08499 (2024).
  * (14) D A Adams, B Schmittmann, and R K P Zia, “Far from equilibrium transport with constrained resources”, J. Stat. Mech.: Theory Exp 2008, P06009 (2008).
  * (15) L. Jonathan Cook and R. K. P. Zia, “Feedback and fluctuations in a totally asymmetric simple exclusion process with finite resources”, J. Stat. Mech.: Theory Exp 2009, P02012 (2009).
  * (16) L. Jonathan Cook, R. K. P. Zia, and B. Schmittmann, “Competition between multiple totally asymmetric simple exclusion processes for a finite pool of resources”, Phys. Rev. E 80, 031142 (2009).
  * (17) M. Ha and M. den Nijs, “Macroscopic car condensation in a parking garage”, Phys. Rev. E 66, 036118 (2002).
  * (18) C. A. Brackley, M. C. Romano, C. Grebogi, and M. Thiel, “Limited resources in a driven diffusion process”, Phys. Rev. Lett. 105, 078102 (2010).
  * (19) C. A. Brackley, M. C. Romano, and M. Thiel, “Slow sites in an exclusion process with limited resources”, Phys. Rev. E 82, 051920 (2010).
  * (20) C. A Brackley, M. C. Romano, and M. Thiel, “The dynamics of supply and demand in mRNA translation”, PLoS computational biology 7, e1002203 (2011).
  * (21) L. Ciandrini, I. Neri, J. C. Walter, O. Dauloudet, and A. Parmeggiani, “Motor protein traffic regulation by supply–demand balance of resources,” Phys. biol. 11.
  * (22) M. C. Good, M. D. Vahey, A. Skandarajah, D. A. Fletcher, and R. Heald, “Cytoplasmic volume modulates spindle size during embryogenesis”, Science 342, 856–860 (2013).
  * (23) J. Hazel, K. Krutkramelis, P. Mooney, M. Tomschik, K. Gerow, J. Oakey, and JC Gatlin, “Changes in cytoplasmic volume are sufficient to drive spindle scaling”, Science 342, 853 (2013).
  * (24) A. Robson, K. Burrage, and M. C. Leake, “Inferring diffusion in single live cells at the single-molecule level”, Philos Trans R Soc (Lond) B Biol Sci. 368(1611), 20120029 (2013) [doi: 10.1098/rstb.2012.0029].
  * (25) R. Lipowsky, S. Klumpp, and T. M. Nieuwenhuizen, “Random walks of cytoskeletal motors in open and closed compartments”, Phys. Rev. Lett. 87, 108101 (2001).
  * (26) S. Klumpp and R. Lipowsky, “Traffic of molecular motors through tube-like compartments”, Journal of Statistical Physics 113, 233 (2003).
  * (27) S. Klumpp and R. Lipowsky, “Asymmetric simple exclusion processes with diffusive bottlenecks”, 70, 066104 (2004).
  * (28) L.D. Fernandes and L. Ciandrini, “Driven transport on a flexible polymer”, Phys. Rev. E 99, 052409 (2019).
  * (29) O. Dauloudet et al., “Modelling the eﬀect of ribosome mobility on the rate of protein synthesis”, Eur. Phys. J. E 44, 19 (2021).
  * (30) I. R. Graf and E. Frey, “Generic Transport Mechanisms for Molecular Traffic in Cellular Protrusions”, Phys. Rev. Lett. 118, 128101 (2017).
  * (31) M. Bojer, I. R. Graf, and E. Frey, “Self-organized system-size oscillation of a stochastic lattice-gas model”, Phys. Rev. E 98, 012410 (2018).
  * (32) A. Goswami, M. Chatterjee, and S. Mukherjee, “Steady states and phase transitions in heterogeneous asymmetric exclusion processes”, J. Stat. Mech. - Th. Exp. 123209 (2022).
  * (33) S. Pal, P. Roy and, A. Basu, “Availability, storage capacity, and diffusion: Stationary states of an asymmetric exclusion process connected to two reservoirs”, arXiv: 2308.08384 .
  * (34) A. Haldar, P. Roy, and A. Basu, “Asymmetric exclusion processes with fixed resources: Reservoir crowding and steady states”, Phys. Rev. E 104, 034106 (2021).
  * (35) B. Pal and A. K. Gupta, “Reservoir crowding in a resource-constrained exclusion process with a dynamic defect”, Phys. Rev. E 106, 044130 (2022).
  * (36) N. Bhatia and A. K. Gupta, “Role of site-wise dynamic defects in a resource-constrained exclusion process”, Chaos, Solitons and Fractals 167, 113109 (2023).
  * (37) A. Gupta, B. Pal, and A. K. Gupta, “Interplay of reservoirs in a bidirectional system”, Phys. Rev. E 107, 034103 (2023).
  * (38) R A Blythe and M R Evans, “Nonequilibrium steady states of matrix-product form: a solver’s guide”, J. Phys. A 40, R333 (2007).
  * (39) T. Reichenbach, T. Franosch, and E. Frey, “Domain wall delocalization, dynamics and fluctuations in an exclusion process with two internal states”, Eur. Phys. J. E 27, 47 (2008).