# Asymmetric exclusion process with global hopping

Chikashi Arita, Jérémie Bouttier, P. L. Krapivsky, Kirone Mallick

► To cite this version:

Chikashi Arita, Jérémie Bouttier, P. L. Krapivsky, Kirone Mallick. Asymmetric exclusion process with global hopping. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2013, 88 (4), pp.042120. 10.1103/PhysRevE.88.042120. hal-00871967

HAL Id: hal-00871967

https://hal.science/hal-00871967

Submitted on 20 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Asymmetric exclusion process with global hopping

Chikashi Arita,¹,² Jérémie Bouttier,¹,³ P. L. Krapivsky,¹,⁴ and Kirone Mallick¹

¹Institut de Physique Théorique, IPhT, CEA Saclay and URA 2306, CNRS, 91191 Gif-sur-Yvette Cedex, France

²Theoretische Physik, Universität des Saarlandes, 66041 Saarbrücken, Germany

³Département de Mathématiques et Applications, École Normale Supérieure, 75230 Paris Cedex 05, France

⁴Department of Physics, Boston University, Boston, Massachusetts 02215, USA

(Received 16 July 2013; published 10 October 2013)

We study a one-dimensional totally asymmetric simple exclusion process with one special site from which particles fly to any empty site (not just to the neighboring site). The system attains a nontrivial stationary state with a density profile varying over the spatial extent of the system. The density profile undergoes a nonequilibrium phase transition when the average density passes through the critical value $1 - [4(1 - \ln 2)]^{-1} = 0.185\,277...$, viz., in addition to the discontinuity in the vicinity of the special site, a shock wave is formed in the bulk of the system when the density exceeds the critical density.

DOI: 10.1103/PhysRevE.88.042120

PACS number(s): 02.50.—r, 05.40.—a, 05.70.Ln

# I. INTRODUCTION

The asymmetric simple exclusion process (ASEP) is a well-studied model of low-dimensional transport of particles with hard-core interactions. This model has become a standard tool in the context of low-dimensional transport, and it is commonly used to represent the motion of molecular motors or more generally enzymes along cytoskeletal fibers. It is interesting to recall that the ASEP was originally proposed in 1968 as a model of the kinetics of biopolymerization in RNA templates [1]. Since then, the ASEP has been extensively studied and it has achieved the status of a paradigm in statistical physics (see [2–7] and references therein). It is fair to say that the ASEP is one of the simplest nontrivial models of a nonequilibrium process: it is an interacting N-body system without detailed balance, which reaches a non-Gibbsian stationary state with nonvanishing currents. In the past 40 years, exact solutions for the ASEP in various contexts have been obtained thanks to increasingly sophisticated and elegant mathematical methods [8–15]. Conversely, these theoretical studies together with experimental progress in micromanipulations have triggered a renewed interest in biophysical applications of the ASEP in recent years [16–24] (see [25] for a review).

The presence of a nonvanishing current in a stationary state results in the transport of information from one part of a system to another, leading to long-range correlations. Henceforth, the behavior of the ASEP is highly sensitive to boundary couplings or to the presence of local defects. In particular, the ASEP can undergo boundary-induced transitions even in one dimension [8], which cannot occur in one-dimensional equilibrium systems with short-range interactions. Hence, simple deformations of the ASEP are often very challenging. These models are important since they enable us to gain insight into macroscopic and microscopic behaviors. One such seemingly innocent generalization is to introduce blockages. This is easy to realize for the ASEP on a periodic one-dimensional lattice where particles undergo a biased nearest-neighbor hopping to empty sites; if the ring is homogeneous, a stationary state with all permissible configuration being equiprobable is reached (on the macroscopic level, the density is uniform in the stationary state). Suppose that the hopping rate between one pair of neighboring sites is suppressed by a certain factor. This

inhomogeneous model with two parameters (the suppression factor and the average density) is a well-known challenge that has not been exactly solved [26–28]; see [29] for a review and recent progress. The chief message is clear: A local change of the environment (one special blockage bond) can result in a global change, namely the system segregates into high- and low-density phases. However, details of the density profiles as well as the stationary-state correlations are still unknown [29].

In the present work, we consider a variant of the ASEP on a one-dimensional lattice, namely on a ring with L sites. The total number of particles N is conserved by dynamics. The system is assumed to be homogeneous apart from one special site from which the particle can fly to any other site. Hence, in contrast to Refs. [26–28], we investigate the effect of a special site rather than a special bond. Our basic motivation is to define a simple model with a *unique* parameter, the average density $\rho = N/L$ (or equivalently the fraction of the vacant sites $v = 1 - N/L$), which can exhibit different regimes with respect to the different values of $\rho$. We show that introducing one special site with global hopping generates a nontrivial phase diagram and leads to a subtle selection mechanism of the stationary density profile.

The outline of this work is as follows. In Sec. II, we explain the dynamical rules of our model. In Sec. III, we write down the governing equations characterizing large-scale stationary hydrodynamic behavior. In Secs. IV and V, we investigate two types of density profiles that can occur, and we show that naive reasoning does not provide a suitable selection mechanism for the mean-field equations. A more precise analysis (Sec. VI) allows us to establish the phase diagram of the system, which, despite the simplicity of the model, turns out to be rather rich. Concluding remarks are given in Sec. VII.

# II. THE MODEL

We consider the totally asymmetric simple exclusion process (TASEP) on a one-dimensional lattice with L sites, which are labeled by $j = 1, ..., L$. Each site j is either occupied by one particle ($\tau_j = 1$) or empty ($\tau_j = 0$). In an infinitesimal time interval $dt$, a particle on site $j$ ($1 \le j \le L - 1$) attempts to hop forward to site $j + 1$ with probability
<figure><img src="image_2.png" /><figcaption>FIG. 1. The TASEP with global hopping from one site.</figcaption></figure>

$dt$; the attempt is successful only when the target site is empty. A particle occupying the rightmost site $L$ can fly to any site $j = 1, \dots, L-1$; a flying attempt is made with probability $dt$ and is allowed if the target site (chosen uniformly at random) is empty; see Fig. 1. These rules define our basic model. The total number of particles $N$ is a conserved quantity. The only free parameter in this model is the density $\rho = N/L$, or equivalently the fraction $v = 1 - \rho$ of vacant sites.

Monte Carlo simulations show that the system attains an interesting stationary state which exhibits a qualitative change when the average density passes through the critical value; we will show below that the critical density equals $1 - [4(1 - \ln 2)]^{-1} = 0.185277...$. Our approach can be extended to a generalized system in which the special site is unique not only since the hopping from this site is global, but also the hopping rate is arbitrary. More precisely, we postulate that in an infinitesimal time interval $dt$, a particle at the special site $L$ attempts to fly with rate $\beta dt$, where $\beta$ is a positive real number, and the target site is chosen uniformly at random among the $vL$ empty sites. We shall derive the phase diagram of the full model in the $(v, \beta)$ plane. The original model corresponds to the case $\beta = v$. The goal of this article is to characterize analytically the phase diagram of the model. We note that, by regarding holes as particles and particles as holes, the process we study is a many-particle generalization of the stochastic-resetting problem [30].

### III. LARGE-SCALE HYDRODYNAMIC BEHAVIOR

A macroscopic (or hydrodynamic) description is valid when the average density slowly varies. In such situations, instead of the density per each site $\rho_j = \langle \tau_j \rangle$, one can study the quantity $F(x)$, which depends on the macroscopic variable $x = j/L$. The macroscopic density can be defined as

$$F(x) = \frac{1}{2\ell + 1} \sum_{k=j-\ell}^{j+\ell} \rho_k, \quad (1)$$

where $1 \ll \ell \ll L$. We study $F(x)$ by employing a hydrodynamic approach, mostly the Eulerian inviscid treatment. Such an approach cannot provide a fair description of the regions inside shock waves and boundary layers, but it can be trusted away from such narrow regions. In the stationary state, we have

$$\frac{d}{dx} F(1-F) = \frac{\beta}{v} \rho_L(1-F). \quad (2)$$

The average density $\rho_L$ at the special point $j = L$ requires a special treatment: A boundary layer could appear near this point, and the limiting hydrodynamic density (we denote it by $R$) may be different from the true microscopic density, i.e., we

have generically

$$\lim_{x \to 1} F(x) = R \neq \rho_L. \quad (3)$$

On the left boundary, we have

$$F(0) = 0. \quad (4)$$

This boundary condition can be understood as follows. A site located near the left boundary can be occupied only if the particle leaving the special site at $L$ lands on it. This is a rare event (of probability $\sim 1/L$) and therefore the microscopic density $\rho_i$ ($i = 1, 2, \dots$) vanishes in the limit $L \to \infty$. This observation is also supported by numerical simulations.

Integrating (2) over the interval $0 < x < 1$ and using (4), we obtain

$$\beta \rho_L = R(1-R). \quad (5)$$

This is the current of flying particles from the site $L$. Inserting (5) back to (2), we arrive at the equation

$$\frac{d}{dx} F(1-F) = \frac{R(1-R)}{v} (1-F), \quad (6)$$

which will be our chief governing differential equation. Note that $\beta$ does not appear in this equation.

Equation (6) is akin to the one obtained by Parmeggiani, Franosh, and Frey (PFF) in their study of the TASEP with open boundaries and Langmuir kinetics (adsorption onto empty sites and desorption from occupied sites) in the bulk [31,32] where the total number of particles is not conserved. We emphasize that in our system there is neither injection nor extraction of particles (see also [33–36] for related models). A more precise discussion of the relation between our system and the PFF model will be given in Sec. VI.

In the following two sections, we shall study the possible solutions to Eq. (6). The selection of the correct density profile $F(x)$ will be achieved by going beyond the simple hydrodynamic limit (6). We shall discuss this in Sec. VI, where the parameter $\beta$ appears again.

### IV. SMOOTH DENSITY PROFILE

We first suppose that the density profile is smooth in the bulk. Solving (6) subject to (4), we obtain

$$2F + \ln(1-F) = \frac{R(1-R)}{v} x. \quad (7)$$

This is valid for all $0 < x < 1$. Combining (7) with (3), we get

$$2R + \ln(1-R) = \frac{R(1-R)}{v}. \quad (8)$$

It is straightforward to check that Eq. (8) may only be satisfied for

$$v \ge v_c = \frac{1}{4(1 - \ln 2)} = 0.8147228 \dots, \quad (9)$$

where $v = v_c$ corresponds to $R(v_c) = \frac{1}{2}$. Conversely, for any $v \ge v_c$, there exists a unique smooth density profile obtained as follows. We let $R = R(v)$ be the unique solution of (8) with $R \in (0, \frac{1}{2}]$. Then, since the function $2F + \ln(1-F)$ is increasing when $F \in (0, \frac{1}{2})$, Eq. (7) determines consistently $F(x)$ on the entire spatial range $0 < x < 1$. [Note that, since
<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color online) Smooth profiles. The top curve is the critical profile, i.e., $v = v_c$, and the lower curve is a subcritical profile with $v = 0.95$. They are realized by Monte Carlo simulations. The symbols □ and ◦ correspond to simulation results with parameters $(v, \beta) = (v_c, 0.6)$ and $(v, \beta) = (0.95, 0.6)$, respectively. We have set $L = 10^5$ sites and took an average over time steps $10^{11} \le T \le 2 \times 10^{11}$.</figcaption></figure>

$2F + \ln(1 - F)$ decreases when $F \in (\frac{1}{2}, 1)$, a smooth solution cannot exceed the density $\frac{1}{2}$. When $v = v_c$, we have

$$2F + \ln(1 - F) = x(1 - \ln 2). \qquad (10)$$

In Fig. 2, we plot the critical profile (10) and a subcritical profile (7) arising at a certain $v > v_c$. Note that the critical density profile (10) has a singularity near $x = 1$, viz., $\frac{1}{2} - F(x) \sim \sqrt{1 - x}$. The analytical expressions (7) and (10) are confirmed by Monte Carlo simulations (see further discussion in Sec. VI).

We performed simulations using the following algorithm. In each time step $T$, we randomly choose one site $j$ among $L$ sites. If site $j \le L - 1$ is chosen and $(\tau_j, \tau_{j+1}) = (1, 0)$, we move the particle at site $j$ as $(\tau_j, \tau_{j+1}) = (0, 1)$. If site $L$ is chosen and $\tau_L = 1$, we make the particle at site $L$ fly with probability $\beta$ to a randomly chosen empty site. [When $\beta > 1$, we have to modify the strategy. If site $j \le L - 1$ is chosen and $(\tau_j, \tau_{j+1}) = (1, 0)$, we move the particle at site $j$ with probability $p$. If site $L$ is chosen and $\tau_L = 1$, we make the particle at site $L$ fly with probability $p\beta$. We need to choose $p$ such that $p\beta \le 1$.]

## V. SHOCK WAVE FORMATION

When the fraction of vacant sites is below the critical value, $v < v_c$, the density profile is no longer described by Eq. (7) on the entire spatial range $0 < x < 1$. In principle, there may be several patches where the density increases but remains smaller than $\frac{1}{2}$, alternating with patches where the density decreases and remains higher than $\frac{1}{2}$. (Note that, in the argument of this section, we do not need to impose the restriction $v < v_c$.)

The only consistent stationary arrangement has precisely two such patches, a low-density patch on the left and a high-density patch on the right. To appreciate this assertion, we recall that in asymmetric exclusion processes, shock waves result from jumps from low to high density, while jumps

from high to low density yield rarefaction waves which are nonstationary (the width of the region covered by a rarefaction wave increases linearly in time). The densities on both sides of a stationary shock wave at position $s$ are related by current conservation (the Rankine-Hugoniot condition), which for the TASEP implies that the densities before and after the shock sum up to 1. We denote the densities before and after the shock by $r$ and $1 - r$, respectively; from the above discussion, we expect that $r \le \frac{1}{2}$. We confirm this bound by Monte Carlo simulations.

In the region before the shock, we have, by (6) and (4),

$$2F + \ln(1 - F) = \frac{R(1 - R)}{v} x \quad (0 < x < s). \qquad (11)$$

The location of the shock $s$ and the density $r$ just before the shock are related via

$$2r + \ln(1 - r) = \frac{R(1 - R)}{v} s. \qquad (12)$$

In the complementary spatial interval $(s, 1)$, the solution to Eq. (6) subject to the boundary condition $\lim_{x \downarrow s} F(x) = 1 - r$ is

$$2(F - 1 + r) + \ln \frac{1 - F}{r} = \frac{R(1 - R)}{v} (x - s). \qquad (13)$$

Combining this solution with the boundary condition (3), we obtain

$$2(R - 1 + r) + \ln \frac{1 - R}{r} = \frac{R(1 - R)}{v} (1 - s). \qquad (14)$$

We thus have two equations, (12) and (14), for three variables $s, r, R$. We must also keep in mind the bounds on these variables (following from the above discussion):

$$0 \le s \le 1, \quad 0 \le r \le \frac{1}{2} \le R \le 1 - r. \qquad (15)$$

Adding (12) and (14), we obtain

$$4r + \ln \frac{1 - r}{r} = \Psi(R, v), \qquad (16)$$

where we have used the shorthand notation

$$\Psi(R, v) = \frac{R(1 - R)}{v} + 2(1 - R) - \ln(1 - R). \qquad (17)$$

The variable $s$ can then be recovered from either (12) or (14).

The left-hand side of Eq. (16) is a monotonically decreasing function of $r$, thus $r$ is determined according to a given value of $R$. Then $s$ is fixed by inserting the values of $R$ and $r$ into (12) or (14). On the other hand, the right-hand side $\Psi$ of Eq. (16) is not always monotonic. For $v \ge \frac{1}{2}$, $\Psi$ monotonically increases, where the maximal $r$ is achieved by $R = \frac{1}{2}$. For $v < \frac{1}{2}$, $\Psi$ decreases on the interval $\frac{1}{2} < R < 1 - v$ and increases on the interval $1 - v < R < 1$. Thus $r$ takes its maximal value at $R = 1 - v$. For example, the dotted lines for $v = 0.1$ in the top and middle graphs of Fig. 3 show that $\Psi$ takes the minimal value and $r$ takes the maximal value at $R = 1 - v = 0.9$. In both cases ($v \ge \frac{1}{2}$ and $v < \frac{1}{2}$), the infimum $r = 0$ is achieved by the limit $R \to 1$, where the shock position $s$ approaches $v$. In the middle and bottom graphs of Fig. 3, we notice that the dashed curves for $v = 0.9$ reach $r = 0.5$ and $s = 1$, respectively. There exists a minimum $R_m$ such that $r$
<figure><img src="image_4.png" /><figcaption>FIG. 3. The right-hand side of Eq. (16) (top), the density just before the shock r (middle), and the position of the shock s (bottom) are determined by R. The fraction of the vacant sites is chosen as the following four values: v = 0.1 (dotted), 0.5 (bold), v_c, and 0.9 (dashed).</figcaption></figure>

and s satisfy the condition (15): for v < 1/2, R_m is 1/2, and for v ≥ 1/2, R_m is the solution to

$$2(1 - R_m) + \ln R_m = \frac{R_m(1 - R_m)}{v}. \quad (18)$$

We have seen that the density r before the shock and the shock position s are specified by R. In other words, once one gives a value for R ∈ [R_m, 1], the shape of the density

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) The density profiles [Eqs. (11) and (13)] with a shock arising at v = 0.6 and R = 0.9, 0.7, 0.5. They are realized by simulations with β = 0.1 (□), 0.3 (○), 0.7 (Δ), respectively, where L = 10⁵ and the averages over time steps 10¹¹ ≤ T ≤ 2 × 10¹¹ were taken.</figcaption></figure>

profile is uniquely fixed. In this sense, shock profiles construct a one-parameter (R) family of solutions. Figure 4 provides shock profiles with v = 0.6 fixed and R varied. These profiles arise in Monte Carlo simulations, as we will demonstrate in Sec. VI.

The results of the previous section tell us that the shape of the smooth profile is uniquely determined when v ≥ v_c. However, in this region, one can also construct shock profiles by patching two smooth solutions. For v < v_c, we have shown that the smooth profile ceases to exist and the solution ought to have a shock. Shock profiles are characterized by the parameter R. A selection mechanism has to be found to explain when and why the system “prefers” a continuous or discontinuous profile, and how the parameter R is determined by (v, β).

We have argued that the correct solution either has no shock or a single shock. Multiple shocks may arise in different models: for instance, in a system with open boundaries and Langmuir kinetics in the bulk, solutions with two shocks have been detected for a certain exclusion process with interaction between neighboring particles [34].

## VI. THE SELECTION MECHANISM

We have obtained a smooth solution and a one-parameter family of shock solutions of the Euler equation (6). In a certain range of filling fractions, we also have a competition between smooth and discontinuous profiles.

First, we remark that the analysis carried out in the previous sections was based on the differential equation (6) with the boundary condition (4); in particular, the flying rate β has played no role in discussing the two types of solutions. Therefore, we must incorporate in our analysis the boundary condition in the vicinity of the site L by taking into account the sum rule (5).

To determine the phase diagram of the system, we shall keep a subleading viscosity term to the mean-field equation, construct the right boundary layer, and apply a stability analysis to it. The stationary hydrodynamic equation with
viscosity is given by (see the Appendix for a derivation)

$$ \frac{1}{2L} \frac{d^2}{dx^2} F - \frac{d}{dx} F(1-F) + \frac{\beta}{v} \rho_L(1-F) = 0. \quad (19) $$

We are interested in what happens in the neighborhood of $x = 1$. Defining the local variable $x = 1 - \frac{\xi}{2L}$ and writing $f(\xi) = F(1 - \frac{\xi}{2L})$, we get

$$ \frac{d^2}{d\xi^2} f + \frac{d}{d\xi} f(1-f) = 0 \quad (20) $$

as $L \to \infty$ (after neglecting terms of order $L^{-1}$). Integrating this equation with the condition $f(+\infty) = R$, we obtain

$$ \frac{df}{d\xi} = R(1-R) - f(1-f) = (R-f)(1-R-f). \quad (21) $$

This is the boundary layer equation that has to be solved subject to the initial condition $f(0) = \rho_L = \frac{R(1-R)}{\beta}$. We note that (21) has two fixed points $R$ and $1-R$. Setting $f = R+h$ (with $h \ll 1$), we have

$$ \frac{d}{d\xi}h = (2R - 1)h. \quad (22) $$

Therefore, if $R < 1/2$, the fixed point $f = R$ is stable and $f = 1 - R$ is unstable. If $R > 1/2$, the fixed point $f = R$ is unstable and $f = 1 - R$ is stable. We also recall (Secs. IV and V) that smooth profiles are characterized by the fact that $R \le 1/2$, whereas for shock solutions we always have $R \ge 1/2$.

To determine the phase diagram, consider first the situation in which a shock profile is selected ($R \ge 1/2$). If $R > 1/2$, then, according to Eq. (22), $f = R$ is unstable and no boundary layer matching $f(0) = \rho_L = \frac{R(1-R)}{\beta}$ with $R$ can exist. Only one possibility remains: $R = \rho_L = \frac{R(1-R)}{\beta}$. This implies

$$ R = 1 - \beta, \quad (23) $$

which is realized only when $\beta < 1/2$ because of the assumption $R > 1/2$.

If the bulk value is $R = 1/2$, it can be connected to the boundary value $\rho_L = \frac{1}{4\beta}$ through a boundary layer. Using Eq. (21), we observe that $f(\xi)$ is increasing with $\xi$, i.e., $\rho_L < R = 1/2$, which in turn implies $\beta > 1/2$; thus, for $\beta$ greater than $1/2$, the selected shock profile is the one with $R = 1/2$.

To summarize, when the selected profile has a shock, the parameter $R$ is given by

$$ R = \begin{cases} \frac{1}{2} & (\beta > \frac{1}{2}), \\ 1 - \beta & (\beta \le \frac{1}{2}), \end{cases} \quad (24) $$

where the true density of the rightmost site is given by

$$ \rho_L = \begin{cases} \frac{1}{4\beta} \neq R & \text{(boundary layer)} \quad (\beta > \frac{1}{2}), \\ 1 - \beta = R & \text{(no boundary layer)} \quad (\beta \le \frac{1}{2}). \end{cases} \quad (25) $$

To understand when the shock profile with parameter $R$ given by (24) can be selected, we set the shock position $s = 1$ and we observe that a shock can be found only when $v < v_c(\beta)$,

<figure><img src="image_6.png" /><figcaption>FIG. 5. The phase diagram of the model. The full lines represent the boundaries between the three phases of the model: a phase with a smooth density profile and two domains where a shock appears. The two shock phases differ at the right boundary: When $R \neq \rho_L$, a boundary layer is formed; when $R = \rho_L = 1 - \beta$, the solution of the inviscid hydrodynamic equation provides an excellent description of the density profile of the entire system.</figcaption></figure>

with

$$ v_c(\beta) = \begin{cases} v_c = \frac{1}{4(1-\ln 2)} & (\beta > \frac{1}{2}), \\ \frac{\beta(1-\beta)}{2\beta + \ln(1-\beta)} & (\beta \le \frac{1}{2}), \end{cases} \quad (26) $$

which is drawn in Fig. 5.

By contraposition, only the smooth profile is selected for $v \ge v_c(\beta)$. As we have shown in Sec. IV, a smooth profile cannot exist for $v < v_c$. Hence, at this stage, we have not yet determined which profile is selected in the region $v_c < v < v_c(\beta)$, and $\beta < 1/2$. We now show that a smooth profile is not allowed for $v < v_c(\beta)$, $\beta < 1/2$. Let us assume, by reductio ad absurdum, that a smooth profile exists. The parameter $R < 1/2$ is a stable fixed point of (21), and must be connected to $\rho_L = \frac{R(1-R)}{\beta}$ through a boundary layer satisfying Eq. (21). But Eq. (21) has two fixed points $R$ and $1 - R$. Because this is a first-order equation, no solution can cross a fixed point, i.e., the boundary layer that connects $R$ and $\rho_L$ cannot pass through the value $1 - R$. This implies $\rho_L \le 1 - R$, i.e., $R \le \beta$. This condition gives $v \ge v_c(\beta)$. In other words, a smooth profile cannot exist when $v < v_c(\beta)$.

The full phase diagram of the system is represented in Fig. 5. There are three phases: a phase with a smooth density profile and two domains where a shock appears. The shock phases differ with respect to the boundary layer. When $\beta > 1/2$, a boundary layer is formed; for $\beta < 1/2$, the profile is smooth in the vicinity of $x = 1$ (see Fig. 6). The dotted line in the smooth phase is $v = \frac{\beta(1-\beta)}{2(1-\beta)+\ln\beta}$, which corresponds to the accidental case in which there is no boundary layer $\rho_L = R = 1 - \beta$. Hence, generically, in the smooth phase a boundary layer always exists.

From this full phase diagram, we can deduce the phase diagram of the original problem in which we were interested. There, the particle at site $L$ attempts to hop to a randomly chosen site with rate unity, and the jump is accepted if the target site is empty. This corresponds to the diagonal line $\beta = v$ in the phase diagram in Fig. 5. We observe that there
<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) A smooth profile and a shock profile are solutions to Eq. (6) with $v = 0.85$ (solid lines). For the shock profile, we have set the parameter $R = 0.75$. The symbols ◦ correspond to the simulation result with $(v, \beta) = (0.85, 0.25)$, where $L = 10^5$ and the average over time steps $10^{11} \le T \le 2 \times 10^{11}$ was taken. We observe from the simulation that the shock profile is selected with $R = 1 - \beta$.</figcaption></figure>

are also three phases for this special problem: for $v \ge v_c = [4(1 - \ln 2)]^{-1} = 0.8147228...$, the profile is smooth; for $1/2 < v < v_c$, the profile is of shock type with $R = 1/2$ and there exists a boundary layer; for $0 < v \le 1/2$, the profile is of shock type with $R = 1 - v$ and there is no boundary layer. We notice that $R$ is chosen so that $r$ is maximized.

Finally, we mention a relation between our system and the PFF model [31,32]. The PFF model is defined on a finite open lattice by the following rules:

$$ 
\begin{array}{l@{\quad}l@{\quad}l}
10 & \rightarrow & 01 \quad \text{at rate } 1, \\
0 & \rightarrow & 1 \quad \text{at rate } \Omega_a/L, \\
1 & \rightarrow & 0 \quad \text{at rate } \Omega_d/L
\end{array}
$$

in the bulk, and

$$ 0 \rightarrow 1 \quad \text{at rate } \alpha, \quad 1 \rightarrow 0 \quad \text{at rate } \beta $$

at the left and right boundaries, respectively. The rates $\alpha$ and $\beta$ are related to the densities of the left and right reservoir, respectively. We note that the hydrodynamic equation that describes the macroscopic density profiles of our model can be formally mapped into the one corresponding to the PFF model by setting

$$ \alpha = \Omega_d = 0, \quad \Omega_a = \frac{R(1-R)}{v}. \qquad (27) $$

This ensures that at the hydrodynamic level [Eq. (6)], the models coincide. However, we emphasize that the factor $\frac{R(1-R)}{v}$ is not determined by external conditions but rather by $R = F(1)$ and by Eq. (6). Besides, the system we consider is conservative and therefore the mapping to the PFF model can only be approximate at the mean-field level. This fact can be observed in Monte Carlo simulations, which show significant differences between the shock profiles of the two models for a given system size; see Fig. 7. Thus, it would be therefore of interest to explore further the connection between our model and the PFF model [31,32]. We emphasize that whereas the number of particles is strictly conserved in the process we study,

<figure><img src="image_8.png" /><figcaption>FIG. 7. (Color online) Comparison with a special case of the PFF model. The parameters of our model are chosen as $(\beta, v) = (0.6, 0.95)$ (◦, smooth profile) and $(\beta, v) = (0.6, 0.3)$ (□, shock profile), corresponding to the PFF model with $(\beta, \Omega_a) = (0.6, 0.09645...)(\times)$ and $(0.6, 0.8)(+)$, respectively. In both models we have set $L = 10^3$ (and we took an average over time steps $5 \times 10^8 \le T \le 5 \times 10^9$), but we observe a difference around the shock position.</figcaption></figure>

the PFF model exchanges particles with its surroundings. The correspondence cannot be perfect and it is also that for nonequilibrium systems, where the equivalence of ensembles (such as canonical versus grand-canonical) is a subtle issue, especially when shocks or phase separation are involved [37,38]. In particular, fluctuations can be significantly different.

## VII. CONCLUDING REMARKS

We studied an exclusion process with a special site from which nonlocal hoppings are allowed. We showed that a shock is formed when the average density exceeds the critical value. We established a phase diagram of the model by using a hydrodynamic approach, a boundary layer analysis, and an asymptotic matching. It would be interesting to devise a physical principle, perhaps reminiscent to the extremal current principle suggested by Krug [3] (see also [39–42]), that would lead to the same phase diagram. Overall our model provides another manifestation of a remarkable phenomenon, namely that a single defect can have a drastic effect on nonequilibrium steady states.

Although our analysis relies on a mean-field approximation, the results appear asymptotically exact and they are in excellent agreement with Monte Carlo simulations. We performed simulations for large sizes (up to $L = 10^5$) at over 100 points $(v, \beta)$ in the phase diagram. Our model demonstrates again that localized defects in far-from-equilibrium systems can produce global effects that drastically alter the phenomenology of a model. Systems far from equilibrium are therefore very sensitive to boundary conditions and to defects.

Models with periodic and open boundaries have been studied extensively, yet little is known about inhomogeneous systems [26–29]. We believe that the interplay of impurities and currents can lead to rich and unexpected behavior (note that even in thermal equilibrium, the presence of a magnetic impurity can result in the highly nontrivial Kondo effect).
In biological applications, inhomogeneities are unavoidable, and various semirealistic inhomogeneous extensions of the TASEP have been analyzed; see, e.g., [43–45]. It would be very interesting to find and to analyze exactly solvable extensions of the TASEP that are both far from equilibrium and nonhomogeneous.

ACKNOWLEDGMENTS

We thank Cécile Appert-Rolland, Arvind Ayyer, Thierry Bodineau, and Tridib Sadhu for helpful discussions and suggestions. C. Arita acknowledges support from the JSPS fellowship program for research abroad.

APPENDIX: DERIVATION OF THE HYDRODYNAMIC EQUATIONS

Here we explain how the hydrodynamic limit (2) and the viscosity correction can be obtained starting from the microscopic dynamics. This is a standard procedure (see, e.g., [3,6]) which we recall for the sake of completeness.

Exact dynamical equations can be written for the expectation value $\rho_j = \langle \tau_j \rangle$ of $\tau_j$ (i.e., the density of site j) at time t as

$$ \frac{d\rho_j}{dt} = \langle \tau_{j-1}(1-\tau_j) \rangle - \langle \tau_j(1-\tau_{j+1}) \rangle + \frac{\beta}{L v} \langle \tau_L(1-\tau_j) \rangle \quad (j=2,3,\dots,L-1), \quad (A1) $$

$$ \frac{d\rho_1}{dt} = \frac{\beta}{L v} \langle \tau_L(1-\tau_1) \rangle - \langle \tau_1(1-\tau_2) \rangle, \quad (A2) $$

$$ \frac{d\rho_L}{dt} = \langle \tau_{L-1}(1-\tau_L) \rangle - \beta \rho_L. \quad (A3) $$

As usual, Eqs. (A1)–(A3) are part of a hierarchy that couples correlations of a given order to higher-order correlations. Here, this hierarchy will be closed by assuming that in the hydrodynamic limit $L \to \infty$, the local density profile $F(x)$ defined in (1) satisfies a mean-field equation. This may sound like a very bold assumption, but it is known to lead to sound results for exclusion processes.

However, we shall go beyond the Eulerian inviscid treatment by keeping track of viscosity terms that scale as $L^{-1}$. These higher-order contributions will have an important physical role within the boundary layers.

In the bulk, we have

$$ \frac{\partial F}{\partial t} = \frac{1}{2L^2} \frac{\partial^2 F}{\partial x^2} - \frac{1}{L} \frac{\partial}{\partial x} F(1-F) + \frac{\beta}{L v} \rho_L(1-F), \quad (A4) $$

where $x = j/L$ and $\rho_L \equiv \langle \tau_L \rangle$. Rescaling the time as $t \to Lt$ and going to the stationary limit, we obtain

$$ \frac{1}{2L} \frac{d^2 F}{dx^2} - \frac{d}{dx} F(1-F) + \frac{\beta}{v} \rho_L(1-F) = 0. \quad (A5) $$

The viscosity term is a singular perturbation that can be neglected in the regions where $F$ varies smoothly (i.e., outside the shock regions and the boundary layers); Eq. (2) is then recovered. Near $x = 1$, a boundary layer may exist and this fact plays a crucial role in selecting the global density profile: this is precisely what is studied in Sec. VI. Near $x = 0$, we observe $F(0) = 0$ and there is no boundary layer. This is supported by the fact the stationary current near $x = 0$ vanishes in the limit $L \to \infty$. For example, from Eq. (A2) one indeed finds that the current from site 1 to site 2 vanishes as the system size diverges, viz., $\langle \tau_1(1-\tau_2) \rangle \sim L^{-1} \to 0$ when $L \to \infty$.

[1] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, Biopolymers 6, 1 (1968); C. T. MacDonald and J. H. Gibbs, ibid. 7, 707 (1969).

[2] H. Spohn, Large Scale Dynamics of Interacting Particles (Springer-Verlag, New York, 1991).

[3] J. Krug, Phys. Rev. Lett. 67, 1882 (1991).

[4] T. Halpin-Healy and Y.-C. Zhang, Phys. Rep. 254, 215 (1995).

[5] G. M. Schütz, Exactly Solvable Models for Many-Body Systems Far From Equilibrium, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic, London, 2000), Vol. 19.

[6] R. A. Blythe and M. R. Evans, J. Phys. A 40, R333 (2007).

[7] B. Derrida, J. Stat. Mech. (2007) P07023.

[8] B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier, J. Phys. A 26, 1493 (1993).

[9] G. M. Schütz and E. Domany, J. Stat. Phys. 72, 277 (1993).

[10] B. Derrida, M. R. Evans, and D. Mukamel, J. Phys. A 26, 4911 (1993).

[11] B. Derrida, M. R. Evans, and K. Mallick, J. Stat. Phys. 79, 833 (1995).

[12] G. M. Schütz, J. Stat. Phys. 88, 427 (1997); R. J. Harris, A. Rákos, and G. Schütz, J. Stat. Mech. (2005) P08003.

[13] S. Prolhac and K. Mallick, J. Phys. A 41, 175002 (2008); 42, 175001 (2009).

[14] C. Tracy and H. Widom, J. Stat. Phys. 132, 291 (2008); Commun. Math. Phys. 290, 129 (2009); J. Stat. Phys. 140, 619 (2010).

[15] B. Derrida, J. Stat. Mech. (2011) P01030.

[16] A. Garai, D. Chowdhury, D. Chowdhury, and T. V. Ramakrishnan, Phys. Rev. E 80, 011908 (2009).

[17] D. Chowdhury, A. Schadschneider, and K. Nishinari, Phys. Life Rev. 2, 318 (2005).

[18] F. Turci, A. Parmeggiani, E. Pitard, M. C. Romano, and L. Ciandrini, Phys. Rev. E 87, 012705 (2013).

[19] L. Ciandrini, I. Stansfield, and M. C. Romano, Phys. Rev. E 81, 051904 (2010).

[20] J. Tailleur, M. R. Evans, and Y. Kafri, Phys. Rev. Lett. 102, 118109 (2009).

[21] O. Campàs, C. Leduc, P. Bassereau, J. Casademunt, J.-F. Joanny, and J. Prost, Biophys. J. 94, 5009 (2008).

[22] D. Johann, C. Erlenkämper, and K. Kruse, Phys. Rev. Lett. 108, 258103 (2012).

[23] L. Reese, A. Melbinger, and E. Frey, Biophys. J. 101, 2190 (2011); A. Melbinger, L. Reese, and E. Frey, Phys. Rev. Lett. 108, 258104 (2012).

[24] M. Ebbinghaus, C. Appert-Rolland, and L. Santen, Phys. Rev. E 82, 040901 (2010).

[25] T. Chou, K. Mallick, and R. K. P. Zia, Rep. Prog. Phys. 74, 116601 (2011).
[26] S. A. Janowsky and J. L. Lebowitz, Phys. Rev. A 45, 618 (1992).

[27] G. M. Schütz, J. Stat. Phys. 71, 471 (1993).

[28] S. A. Janowsky and J. L. Lebowitz, J. Stat. Phys. 77, 35 (1994).

[29] O. Costin, J. L. Lebowitz, E. R. Speer, and A. Troiani, Bull. Inst. Math. Acad. Sin. (N.S.) 8, 49 (2013).

[30] M. R. Evans and S. N. Majumdar, Phys. Rev. Lett. 106, 160601 (2011).

[31] A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev. Lett. 90, 086601 (2003).

[32] A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev. E 70, 046101 (2004).

[33] M. R. Evans, R. Juhász, and L. Santen, Phys. Rev. E 68, 026117 (2003).

[34] V. Popkov, A. Rákos, R. D. Willmann, A. B. Kolomeisky, and G. M. Schütz, Phys. Rev. E 67, 066117 (2003).

[35] N. Mirin and A. B. Kolomeisky, J. Stat. Phys. 110, 811 (2003).

[36] R. Jiang, Y.-Q. Wang, A. B. Kolomeisky, W. Huang, M.-B. Hu, and Q.-S. Wu, Phys. Rev. E 87, 012107 (2013).

[37] P. F. Arndt, T. Heinzel, and V. Rittenberg, J. Phys. A 31, L45 (1998); J. Stat. Phys. 97, 1 (1999).

[38] N. Rajewsky, T. Sasamoto, and E. R. Speer, Physica A 279, 123 (2000).

[39] A. B. Kolomeisky, G. M. Schütz, E. B. Kolomeisky, and J. P. Straley, J. Phys. A 31, 6911 (1998).

[40] V. Popkov and G. M. Schütz, Europhys. Lett. 48, 257 (1999).

[41] J. S. Hager, J. Krug, V. Popkov, and G. M. Schütz, Phys. Rev. E 63, 056110 (2001).

[42] M. Dierl, P. Maass, and M. Einax, Phys. Rev. Lett. 108, 060603 (2012); M. Dierl, M. Einax, and P. Maass, Phys. Rev. E 87, 062126 (2013).

[43] T. Chou and G. Lakatos, Phys. Rev. Lett. 93, 198101 (2004).

[44] J. Dong, B. Schmittmann, and R. Zia, J. Stat. Phys. 128, 21 (2007).

[45] C. A. Brackley, M. C. Romano, C. Grebogi, and M. Thiel, Phys. Rev. Lett. 105, 078102 (2010); C. A. Brackley, M. C. Romano, and M. Thiel, Phys. Rev. E 82, 051920 (2010).