Phase transitions in the driven lattice gas (TASEP) with repulsive energies

We report on non-equilibrium phase transitions in the driven lattice gas model totally asymmetric exclusion processes caused by nearest-neighbor (NN) repulsions between particles, ɛ > 0, employing extensive Monte Carlo simulations and mean-field calculations. The diagra of phase regimes has five different regimes which are separated by distinct boundaries. We discuss the density of particles and the current which depend differently in the various regimes on the system’s entry and exit rates of particles, α and β, respectively. The steady state can be described by quasi-particles called ‘soft dimerons’, where each dimeron consists of a particle dragging an empty NN site on its left. The asymptotic state of dimerons at ε→∞ is equivalent to the state of hard dimers at ɛ = 0. Since the model is a uni-directional fluid, it has an internal hydrodynamic pressure which increases with increasing inter-particle repulsion ɛ. This leads to a collective effect which causes phase transitions between high and low densities of particles at certain critical points ɛ c .

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The totally asymmetric exclusion processes (TASEP) model belongs to a class of driven diffusive systems [1,2] which are of importance in various fields of physics and biology [3][4][5]. They often serve as simple models for biological transport phenomena [6,7], traffic flow [5,8] and fast ionic conductors [9]. In the simplest version of the model, the hopping transport of the particles through a finite one dimensional lattice is governed by the rates: (i) the entry-rate, α, for a particle to enter the first site of the lattice on the left, if it is empty, (ii) the hop-rate for a particle to move right to its empty nearest-neighbor (NN) site in the bulk and (iii) the exit-rate, β, for a particle to leave the system from the last site of the lattice. If α and β are expressed in the units of the hop-rate in the bulk, then the latter can be taken to be unity. The system is then characterized by only two parameters α and β. In the stationary state, the phase diagram in the (α, β)-plane [10] consists of three distinct phases namely, low-density (LD) phase corresponding to α < 1/2 and β > α, high-density (HD) phase corresponding to β < 1/2 and α > β and maximal current (MC) phase corresponding to α, β > 1/2. The phase diagram is sketched in figure 1 for convenience. The point, (α * , β * ) = (1/2.1/2), in the phase diagram at which these three phases meet may be called the 'triple point'. This phase diagram and the related behavior of the current and the density in the system have been firmly established by exact mathematical analysis [1,2].
In the present study we provide numerical and theoretical results for non-equilibrium phase transitions in the TASEP model which are brought about by tuning the strength of the NN repulsion (ε > 0) in an open one-dimensional driven system of particles. This is in contrast to one-dimensional systems with NN interactions in equilibrium [11] where phase transitions do normally not exist. A second motivation is to shed some lights to the mechanism of fast ion conduction which is of importance in science and technology [12][13][14][15][16]. Possible implications on fast ion transport is briefly addressed in our conclusions.
The outline of the paper is as follows. First we describe in section 2 the model and its corresponding simulation techniques. In the subsequent section 3 we describe and discuss a new diagram of phase regimes which consists of five different regimes with their corresponding boundaries between them. In the following two sections we describe and discuss the transitions between the different regimes. In particular in section 5 we address the phase transitions driven by the hydrodynamic pressure which is induced by the NN repulsion ε. Since we have phase transitions in the system, we discuss briefly in appendix B the dynamics of the density profiles and shock waves close to their critical points. Our main focus are on phase transitions induced by the hydrodynamic pressure which is discussed in section 5.

Model and simulation technique
We consider a finite one-dimensional lattice of length L occupied randomly by N particles with mutual hard-core repulsion. In contrast to equilibrium systems, the bulk behavior is sensitive to the boundaries [17] which induce phase transitions in one-dimensional system and which may result in complex phase behavior. In the present model, the 'source' of particles into the system is not provided by a reservoir at a certain density, but rather by assigning an entry-rate α, for a particle to enter the site x = 1. Similarly for the 'sink' of the system at x = L where the particles disappear and where the exit-rate is β for a particle to leave the system from its site x = L. Beside the commonly used mutual hard-core repulsion between particles we introduce a tunable NN repulsion ε > 0 which leads to considerable differences as compared to the classical case ε = 0. Our simulation method is based on the kinetic Monte Carlo (KMC) algorithm that was first proposed for simulating the kinetic evolution of a lattice-gas system [18] and also a variety of reaction-diffusion systems [19,20]. More recently, we have used the same generalized version of the algorithm for studying a TASEP model of hard rods (extended) particles [21,22]. In the present study the particle is a monomer which does not occupy more than a single site of the lattice.
We consider N particles with NN hard-core repulsion ε > 0. We associate with each site of the lattice an occupancy variable, n i = 0, 1 (i = 1, 2, · · · , L). The entry-rate, ω 0 , for a particle to occupy the first site, if empty (i.e., n 1 = 0), is equal to α or α e −ε depending on whether n 2 = 0 or 1. The hop-rate, ω i , for a particle to move in the bulk from site i to site i + 1 is (i) unity if n i−1 = n i+2 = 0, (ii) e −ε if n i−1 = 0, n i+2 = 1, (iii) e ε if n i−1 = 1, n i+2 = 0. Finally, the exitrate, ω L , is equal to β or β e ε depending on whether n L−1 = 0 or 1. Similar transition rates had been introduced recently [23]. The simulation process consists of a series of one of the three events namely, the entry, the hopping and the exit of a particle. The choice of an event and the corresponding time-update are based on the standard KMC protocol involving the cumulative rates, R j ≡ ∑ j i =0 w i , for j = 0, 1, 2, · · · N, so that the ratio,R j ≡ R j /R N , is the probability for the j th particle to move. Note that the probability for a particle to enter the system is given byR 0 . Let r 1 and r 2 be two uniform random numbers in the range (0, 1). If r 1 ⩽R 0 , then a particle enters the system; else, the value of j for which the inequalityR j−1 < r 1 R N ⩽R j holds points to the particle j that makes a jump. Under the implicit assumption that it is a Poisson process, the time associated with this event is given by t m = − ln(r 2 )/R N where m is the mth KMC cycle. The total time is t = ∑ m t m . This way of choosing a particle to make a jump is repeated a large number of times, say ∼10 7 -10 8 times (KMC cyles) in a single KMC run. Such a time-update protocol helps us study the kinetic evolution towards the stationary state of the system parameters-for example, the density and current profiles, the mean squared displacement of a tagged particle. The flow of particles is described by the global current J(ε). The global current is defined by J = (⟨n⟩ i + ⟨n⟩ o )/2⟨t⟩ where ⟨n⟩ i and ⟨n⟩ o are the average numbers of entering and exiting particles, respectively, during the average time ⟨t⟩.

Diagram of phase regimes and boundaries
We present in the following section the 'diagram of phase regimes' for ε > 0 which is the analog to the classical phase diagram for ε = 0 (figure 1). The five regimes are denoted by P α , P β , P m , P λ and P αβ (figure 2) where each regime encompasses the corresponding phases, α, β, m, αβ, respectively, which depend on ε. Phase diagrams as function of ε had been considered previously in various publications [23][24][25]. In the present paper, however, for the sake of clarity and completeness, we define a 'diagram of phase regimes' in figure 2 which consists of regimes rather than of 'phases' as defined previously in the literature. In general, a driven system of repulsively interacting monomers, in particular in the limit of very strong repulsion in its stationary states, has a rich generic (or, full) phase diagram. Such a full phase diagram in the (α, β, ε) space with proper phase boundaries [23][24][25] would present a comprehensive summary of the ε-dependent stationary phases as function of α and β. In the present work, however, we present a more simpler and different approach employing a diagram of phase regimes (figure 2). The diagram of phase regimes can be considered as a type of projection of the full phase diagram onto the (α, β) plane, but without phase boundaries.
The boundaries between the regimes are denoted by γ * and α m and are independent of ε. The transitions between the various regimes are characterised by the density of particles, ρ(ε), the density of dimeron, D(ε), and the current, J(ε). In order to understand the transitions of the densities it is instructive and useful to introduce the concept of a 'quasi-particle' which we call 'dimeron'. This quasi-particle consists of a particle dragging an empty NN site on its left. This dimeron has to be distinguished from hard dimers [21,26,27] which is a particle of length 2 consisting of two lattice sites. In the limit of very large repulsion ε → ∞ the empty left NN site is tagged with the monomer with probability 1 and hence the effective size of a dimeron is 2, i.e. a hard dimer. In the limit ε → 0, the left NN is tagged with probability 0. Hence, the asymtotic state at ε → ∞ is described as an ordered sequence of dimerons. Accordingly, we define D(ε) as the density of dimerons D = N D /L where N D is the number of dimerons N D per size of the system L, or equivalently, D(ε) = ⟨01⟩/L where ⟨01⟩ denotes the average occupation of such NN sites. With increasing repulsion ε the system exhibits a rearrangment of the configurations of the particles from random disordered configurations towards ordered configurations (01) thereby maintaining the fluid character of the particles exemplified by the global current J(ε). Even in the limit ε ≫ 1 the flux of particles is maintained as long as at the exit rate β > 0. The diagram for ε > 0 (figure 2) is different in various aspects from the classical phase diagram (figure 1). However, similar as in the classical case for ε = 0 there exists a coexistence line between α = β = 0 and the 'triple point' α = β = γ * (black thick line) separating the low density regime P λ from the P αβ regime. The transition between the P λ and the P αβ phases is similar as in the classical case ε = 0 (figure 1) where the HD 1 −LD 1 transition is of first order. The dynamics of the particle density ρ(ε) across phase transitions, such as the coexistence line, is dicussed briefly in appendix B where we present and discuss density profiles and shocks. In contrast to the classical case ε = 0 the triple point is smaller than 1/2, but γ * ≈ 0.4142 which is corroborated in all cases ε > 0 which are considered in the following sections. The range of the maximal regime P m is limted by β > γ * and by α > α m ≈ 2. The two other regimes, P α and P β , depend only on α and β, respectively. The regime P αβ depends on both α and β. The transitions at α m , which are denoted by the vertical red dotted line, are continuous. The transition P α − P m is signified by the maximum of the dimeron density D α (ε). The transition P β − P m is characterized by maximum of the current J. This is discussed in details in sections 4 and 5.
In order to locate and describe the boundaries it is useful to consider the densities and current in their asymptotic state at ε → ∞ which we denote by '(∞)'. In this limit many exact results can be employed from known results of 'hard dimers' [26]. Numerically the asymtotic state is reached for ε ≳ 4. Boundary Regimes P β −P m . The horizontal boundary at β = γ * and α > γ * separates the regimes of P α and P m from the regimes P αβ and P β . The transition between P β and the maximal regime P m as function of β is shown in figure 3. There the manifestation of dimerons in the limit of ε → ∞ is corroborated by results for 'hard dimers' [21,26,27] (equation (2)): in the maximal regime P m we cannot have a pair of NN particles and hence the density is maximal and becomes the same as the density of hard dimers [26], and hence the density of dimerons is D m (∞) ( figure 3). In addition Lakatos and Chou [26] have shown that for hard dimers the relation between the density ρ(β) and the current J(β) as function of β is which is in perfect agreement with our simulation data (figure 3). In the limit of very strong repulsion ε → ∞ and with exit rates β > γ * the current becomes maximal, J β (∞) = 0.1716 (equation (4)), and the density minimal, ρ β (∞) = 0.2929 (equation (2)). This is known from the Tonks gas in equilibrium for hard dimers [26] (equation (11) in their paper) that in the limit of very strong repulsion and with increasing exit rate β the current becomes maximal, ). Therefore, in the limit ε → ∞ the densities of particles and the density of dimerons are identical, From the intersection point at (1 − β)/2 = 1/(2 + √ 2) we obtain the exact relation Employing equations (1) and (2) provides the current It should be emphasised that the 'triple' point γ * is valid for all ε > 0 and is independent of ε which is corroborated in the following sections 4 and 5.
Boundary Regimes P λ −P α . In the following paragraph we provide evidence for the vertical boundary at α = γ * and β > γ * which separates the regimes P λ and P α (figure 4, left panel). The low density regime P λ corresponds to the parameter space α < γ * and β > γ * . In this regime the density ρ and the current J increases linearly until they reach at α > γ * and β > γ * the regime P α . In the limit ε → ∞ and α ≫ 1 the density and the current assume their values valid in the maximal regime P m . We have The currents in the P λ and the P α regimes are given by where the current in the maximal regime P m is the same as for hard dimers [21,26]. At α = γ * the current matches J λ (∞) = J α (∞) which corroborates equation (3). The data exhibit a short crossover in the range 0.4 < α < 0.5 indicating the transition around γ * . The P λ − P α transition for weak repulsion (0 < ε < 2) is shown in figure 4(b) in order to show that the transition is continuous: the density ρ(ε) and current J(ε) start for ε ⩾ 0 to deviate continuously from their values at ε = 0.

The α-regime P α and maximal regime P m
The range α > γ * and β > γ * defines the regime of P α (α < α m ) which reaches the maximal regime P m for α > α m and ε > ε D . This is shown in figures 5(a) and (b): for the various αphases (given in different colors) the densities of particles ρ α (ε) and the density of dimerons D α (ε) are given. The density and the dimeron profiles are shown for various values of α in the range 0.45 ⩽ α ⩽ 7. Each curve contains data for different values of 0.5 ⩽ β ⩽ 1 demonstrating its β-independence. First we note that all the profiles converge for ε → ∞ to one single profile ρ m (∞) = D m (∞) = [2 + √ 2] −1 ∼ 0.2929. At α = α m ≈ 2.0 the system enters the 'maximal regime' P m which corresponds to the classical 'maximal current' phase MC for ε = 0. The shapes of the density profiles ρ α (ε) change from a convex shape for α < α m to a concave shape for α > α m . The range of convexity, ε 1 < ε < ε 2 , which is marked by two weak transition points, ε 1 and ε 2 , narrows with increasing α and the two transition points finally merge for α ≫ 1. As an example in figure 5, the arrows at ε 1 and ε 2 are the transition points for α = 0.8. At α m the density profile is linear with negative slope, as indicated by a line (guide-to-the eye) in figure 5(a) for the intermediate range 1 ⩽ ε ⩽ 1.7. This indicates the transition between the P α and P m regimes. For α ≫ 1 the shape of ρ α (ε) assumes a limiting shape with an inflection point at ε D which is shown in figure 5(a) by the black curve for α = 7.0. In contrast, figure 5(b) shows that for α ≫ 1 the shapes of the dimeron profiles D α (ε) converge to a limiting shape with a 'cusp' at its maximum at ε D which separates the regimes of D α (ε) and D m (ε). It is observed that the profiles of D α (ε) become more flat for α → γ * = 0.414 which indicates their approach to the low density regime P λ (α < γ * ). We assume, and it is conceivable, that the 'cusp' of D α (ε) indicates a phase transition between the regimes P α and P m , but without any type of criticality at ε D . For values ε < ε 1 the density of particles ρ α (ε) is very high and hence the density of dimerons is very small. Due to increasing repulsion ε the space between the particles increases and hence the number of dimerons increases as well. For ε 1 < ε < ε 2 there exits a mixture of dimerons with configurations (0,1,0) and (0,1,1), where each particle is dragging an empty NN site on its left, but may or may not have an occupied NN site on its right. The first and the latter configuration we may call hard dimerons and soft dimerons, respectively. This picture changes for ε > ε 2 where for strong repulsion we have an 'ordered fluid' of hard dimerons. There is a flow of hard dimerons, even at very strong repulsion. This is so because as soon as the right-part particle of a dimeron drops out at the exit gate the flow is still able to continue by a locomotive motion where particles hop subsequently.
The current profiles J α (ε) (figure 7) are different from the density profiles ρ α (ε) and D α (ε) as presented above in figure 5. Each profile of J α (ε) has its maximum at ε 1 for given α, and points of inflection at ε 2 . For α → α m the maxima becomes ε 1 → ε J and the inflection points becomes ε 2 → ε D . In contrast to D α (ε), which is characterized by its maximal value of D max ≈ 3.9 at ε D ≈ 1.7, the current is characterized by its maximal value of J max ≈ 0.3 at ε J ≈ 0.68 for α ⩾ α m . For decreasing values of α the peak positions of J α (ε) shift towards ε = 0 and their peak values decrease to J α (0) = 1/4. In the limit ε → ∞, J α (ε) saturates to a value [ √ 2 + 1] −2 ≈ 0.1716. We have calculated employing a cluster mean-field theory (appendix A) the  (7)). The transition points ε 1 and ε 2 for α = 0.8 (red arrows) are the same as in figure 5. For comparison the black dotted curve is the maximal profile of the dimeron density Dα(ε) for α > αm which indicates the difference between the profiles of current and dimeron. Each curve contains data from 4 ⩽ β ⩽ 6 which indicates their independence of β.
profile of the maximal current J m , which is the limiting profile for α ≫ 1 (black curve in figure 7), where ε 2 is the inflection point of J m .

Phase transitions of the density of particles
In the following section we describe and discuss phase transitions of the particle density ρ(ε) in the αβ-regime P αβ (β < γ * , γ * < α < α m ) and in the β-regime P β (β < γ * , α > α m ), which are shown in figure 8. The density ρ αβ (ε) is shown in figure 8(a) for β = 0.1 and 0.1 ⩽ α ⩽ 2.0. There are two transition points, ε c1 and ε c2 , which change their locations towards each other with increasing α and merge to ε c (β) at α ≈ α m where the density enters the P β regime. The transition points ε c (β) depend on β. They are the transition points in the P β regime ( figure 8(b)) between the high density phase ρ β (0) = 1 − β which decreases continuously to the maximal phase ρ m (ε) (black curve). It should be noted that simulation data published previously in Teimouri et al [23] agreed with our data in figure 8(a). For example, the transition in figure 8(a) at ε c1 ≈ 0.7 for α = 0.2 corresponds to the high-to-low density transition in their phase diagram figure 2 [23] for ε = 0.6 at the same value α = 0.2 and β = 0.1, which corroborates the validity of our data. Phase diagrams for other values of ε as function of (α, β) could be constructed but are out of scope in the present work. (a) Density of particles ρ αβ (ε) between the αβ-regime P αβ and the β-regime P β for various α > β at β = 0.1. The coexistence point between high and low density phases is denoted by ε c1 . The green dotted lines correspond to the particle densities ρ λ (ε) in the low density regime P λ (α < γ * ), and the black dotted lines correspond to ρα(ε) in the Pα regime (α > γ * . ρ 0 (ε) denotes the density in the ground state. Note that the change between the dotted green and the dotted black curves takes place at α = γ * ≈ 0.414. In practice the density assumes saturation ρ β (∞) for ε > 4. (b) Density ρ β (ε) in the β-regime P β as function of ε and for various β. Each curve contains data from 2 ⩽ α ⩽ 7 which shows the independence of α. ε D denotes the location of the cusp of the dimeron density Dα(ε) in figure 5(b), and ρm(ε) is the limiting shape as shown in figure 5(a). The pink arrow indicates, as an example, the critical point εc(β) for β = 0.3.
The high density phase, where ρ αβ (0) = 1 − β, exists in the range ε < ε c1 and changes to the low density phase at ε c1 . The low density phase is denoted by ρ 0 (ε) and exists in the range ∆ε ≡ ε c1 < ε < ε c2 . This intermediate phase, ρ 0 (ε), we call the 'ground state'. It is intriguing to propose a 'ground state' similar as in the classical theory of phase transitions [11]. Instead of phase transitions between 'high' and 'low temperature phases', here we encounter genuine phase transitions at ε c1 between 'high' and 'low density phases'. Figure 8(a) shows green and black dotted lines. Both green and black dotted lines match the ground state ρ 0 (ε) for certain α in certain intervals ∆ε. It turns out that the green dotted lines correspond to the densities of particles ρ λ (ε) in the low density regime P λ (α < γ * , β > α) and matches ρ αβ (ε) = ρ λ (ε) in the range ∆ε. The black dotted lines correspond to the particle densities ρ α (ε) in the P α regime (α > γ * , β > γ * , compare figure 5(a)) and match ρ αβ (ε) = ρ α (ε) in the range ∆ε. Hence, the set of green and black dotted lines match ρ αβ (ε) in the range ∆ε. Both the sets are separated by their initial values ρ α (0) = 0.5 and ρ λ (0) = α at ε = 0, and separated at α = γ * ≈ 0.414. It is observed that the critical point ε c1 depends both on α and β, whereas the ground state ρ 0 is independent of β and is an increasing function of α which is shown in figure 9(a) and can be summarized as where ⟨.....⟩ ε denotes the average over the range ∆ε. The transitions betweeen high and low densities at ε c1 is a discontinuous (first-order) phase transition exhibiting (within numerical accuracy) discontinuities at ε c1 . The critical points have been estimated and shown in figure 9(b). This is corroborated by jumps in the density profiles which is discussed in detail in appendix B. Such 'shock' profiles (or synonymously: phase boundaries, domain walls) separate regions of different densities [10,[28][29][30][31][32]. Such phase boundaries which separates high  density and low density phases can be very sharp, hence called a 'shock', i.e. a discontinuity in the density profile

Transitions of dimerons and current
As a corollary and for the sake of completeness we present in the following paragraph results for the density of dimerons D(ε) and the global current J(ε) in the P αβ regime as discussed above for the particle density ρ(ε) in section 5. We present in figure 10 the dimeron profiles for the P αβ and P β regimes. The profiles in figure 10(a) corroborate the phase transitions as shown above in figure 8(a) at the points ε c1 and ε c2 . The two critical points merge to one critical point ε c (β) as α → α m when the system enters the P β regime. This behavior is similar to the density ρ αβ (ε) as shown in figure 8(a). The asymptotic ('ordered') state of dimerons, D β (∞), is reached for ε > ε c2 in variance with the densities ρ αβ (ε) given in figure 8(a). The profiles of D β , as shown in figure 10(b), are equivalent to the density profiles in figure 8(b). They exhibit only weak signatures of phase transitions at ε c (β). As another corollary we present in figure 11 the current profiles for the P αβ and P β regimes. Figure 11(a) corroborates the phase transitions as shown above in figure 8(a) at the points ε c1 and ε c2 . Figure 11(b), however, do not exhibit signatures of phase transitions as compared to figure 8(b). In the maximal regime P m the current assumes a limiting shape, equation (7), which is depicted in the figure by the black curve. A notable difference between the regimes P α (figure 7) and P β is that, at ε = 0, J α (ε) is in the P α region and has the same value 1/4 whereas J β (ε) in the P β regime decrease to zero as β → 0. The latter situation is expected because, as β decreases, the system tends to get increasingly jammed, because the blockage 1/β increases, which leads thus to a vanishingly small values of J(0); in fact, at β = 0, J(0) = 0. Moreover, in the limit ε → ∞, J α (ϵ) saturates to a value [ √ 2 + 1] −2 ∼ 0.1716 whereas J β (ϵ) saturates to β-dependent values.

The hydrodynamic pressure
More insights into the physical reasons for the appearance of the intermediate phase ρ 0 (ε) can be obtained by comparing various characteristic quantities near the transition at ε c1 as ρ(ε), D(ε), and the density of energy E(ε). Defining ⟨11⟩ and ⟨01⟩ as the average local configurations of two neighboring sites, either ('bonded' particles) or one-empty-one-occupied neighboring sites ('unbonded' particles), respectively, then one finds that the point of phase transition at ε c1 is related to the interplay between the density of dimerons D(ε) = ⟨01⟩/L and the contact 'energy' E(ε) = ⟨11⟩/L. More precisely, if we define the number of bonds per particles by E/ρ = ⟨11⟩/N, then we found heuristically that the point of intersection of D(ε) and E(ε)/ρ(ε) determines the location of ε c1 which is shown in figure 12 for the case (without loss of generality) α = 0.8 and β = 0.2. This can be expressed by Since ρ = ⟨01⟩/L + ⟨11⟩/L = D + E, equation (9) can be rewritten as The physical reasons for the intersection of D(ε) and E(ε)/ρ(ε) can be interpreted as the competition between the populations of NN configurations of (0, 1) and (11). Since weak repulsion leads to an increased space between the particles and hence to an increased number (0, 1) configurations, the density of dimerons D(ε) (red line, figure 12) increases accordingly. On the other hand, the average energy density per particle ('bonded particles') E(ε)/ρ(ε) (blue line, figure 12) decreases with stronger repulsion.
The physical reason of the competing behavior of D and E/ρ is due to the increasing internal repulsion ε among the particles. This implies that there exists a concomitant increasing internal pressure. The TASEP model is a uni-directional fluid and has a hydrodynamic pressure which is determined entirely by the motion and density of the fluid, where ρ and v are the density and velocity fields, respectively. If we approximate the fields v(x, t) = J(ρ(x, t))/ρ(x, t) [8] by their steady state values, v = J/ρ, this yields, and by employing the general relation J = ρ(1 − ρ), The replacement of the space-time dependent density and current by their steady-state values is an approximation. Equation (12) addresses the stationary regimes. The pressure Π(ε) is shown in figure 13 for various β from above to below β = γ * . For all cases of β the pressure increases with increasing repulsion ε until the critical point ε c1 (β) and remains at its saturation for ε > ε c1 in the P α regime until a weak transition at ε c2 ( figure 13(a)). The internal pressure increases thereby reducing the density until the pressure reaches a certain threshold at ε c1 where the density matches a low density which is compatible with densities in the P α regime. Since by further increasing the repulsion ε ≫ ε c2 the system has to change from its ground state to reach its asymptotic state at ρ β (∞) or ρ α (∞) where dimerons with configurations (0, 1) are prevailing. This transition at ε c2 is probably continuous and disappears for β > γ * . There is a second transition at ε c2 , where the density of the ground state ρ 0 (ε) and its pressure Π(ε) has to adjust to the pressure and the density in the asymptotic state at ε → ∞, which can be higher or lower depending on α and β, and which is discussed above in section 3. Similar observations of the transitions between the P αβ regime and the low density regime P λ are obtained, but their presentation are omitted here. Similarly, there are transitions of Π β (ε) in the P β regime which are shown in figure 13(b). Similar as in the P αβ phase ( figure 13(a)) the pressure increases until it reaches the value of the pressure in the maximal phase P m (maximal pressure ≈ 0.0732) at ε c (β). This corresponds to the asymptotic case of 'hard dimerons' D m (∞) = ρ m (∞) ( figure 8(b)). Still there is a weak β-dependency of Π β (ε) comparable to observations for ρ β (ε) in figure 13(b). It should be noted that, according to figure 13(b), the pressure exhibits a weak, but significant, transition at ε D from Π β (ε) to Π m (ε).
In summary, since the system in the P αβ regime, which starts at ε = 0 at high density ρ = 1 − β, the hydrodynamic pressure causes a dilution of particles with increasing ε such that at a critical point ε c1 the particles exhibit collectively phase transitions to their ground states which corresponds to the low densities of either the P λ or the P α regimes depending on the combination of α and β. The entry and exit probability of particles, α and β, respectively, determine the critical points.

Conclusions
We have investigated the TASEP model with NN repulsive interactions of strength ε > 0. The diagram of phase regimes has five different regimes separated by boundaries. Among the five different phases there are four new regimes as compared to the classical case with ε = 0 which has three main phases. The present work contains a detailed analysis of energy-driven phase transitions which are genuine phase transitions between high and low densities phases. The phase transitions are caused by the collective inter-particle repulsions ε which yields an increasing hydrodynamic pressure with increasing ε. The steady state is described by quasiparticles called 'soft dimerons' which reach their asymptotic state at ε → ∞ and which are equivalent to the state of hard dimers at ε = 0.
One of the motivations for the present work is the analogy between equilibrium systems where the steady state during energy-driven phase transitions is given by the Boltzmann distribution, and nonequilibrium systems which are kept permanently far from equilibrium and for which no general macroscopic theory analogous to equilibrium thermodynamics exists [3,5,33].
Another motivation and aspect of the present work is related to the inter-particle repulsion and its physical realisation. If one assumes that the screening of ion-ion repulsion can be considered as an ε -dependent particle-particle repulsion, at least for small ε, then the present study may have some impact on our understanding of fast ion conduction in various systems. In fact, the present model shares some similarities with models used to describe the mechanism of fast ion conduction which is of importance in science and technology [12,13] as e.g. in nanofluidic devices such as nanodiodes [14], diffusion of colloids in narrow channels [16], and biological ion channels [15,34,35].
It is known that the single-file transport of ions through a narrow channel, for example through a biological ion channel [15], is quite fast despite the coulombic repulsion between them [36]. Motivated by the desire to know whether strong repulsion alone is responsible for speeding up the flow through a channel we found that the unidirectional steady-state flow of particles is considerably enhanced in a certain range of ε relative to the case ε = 0. The maxima of the current J(ε) at the critical point ε = ε c (β) (compare figure 11) have been estimated relative to their values J(0) = 1/4 and J(0) = β(1 − β) in the regimes P α , P m and P αβ , respectively, for all α We denote by A the 'relative amplification factor' which is shown in figure 14. A good fit to the data in the P α regime is provided by A(α) = 1.2 tanh[(α + 0.05)/0.4]. We see that A in the P α regime saturates to a value ∼4/3 as α → α m . In the same figure, we have also shown the ratio for the P β regime as a function of β, the divergence of which is due to J(0) = β(1 − β) in the denominator. It is of interest to note that the maximal amplification is observed for large entry rate α and small exit rate β (large blockage 1/β), which is encountered in the P β regime ( figure 14). For example, for α = 3.0 and β = 0.1 we obtained A ≈ 2.7.
In conclusion, our results indicate that during ion conduction in nanofluidic devices the flow of ions can be considerably enhanced for moderate repulsions ε > 0 as compared to neutral particles with ε = 0.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Acknowledgments
A B would like to acknowledge the hospitality at the Departamento de Fisica at the Unversidad de Extremadura, Badajoz (Spain).

Appendix A. The maximal current: mean field theory
The behavior of the current had been discussed previously employing a cluster mean-field theory and simulations at moderate repulsive and attractive interactions [23,37]. Here we present detailed calculations of the maximal current employing a different cluster mean-field theory MCAK [22]. In the MC-phase, which we denote by P m (α ⩾ α m ≡ 2; β > γ * ), the stationary current is independent of the entry and the exit rates and, therefore, depends only on the repulsive NN interaction, ε. It is of interest to derive a mean-field expression for the current profile, J(ϵ). Since it depends on the local environment of the hopping particle, we have We have omitted the site indices because, for a given value of ε, the stationary current and density in the system are the same (i) throughout a closed system (lattice-ring), (ii) in the interior (far from the boundary regions) of an open systen. Using the Marlov Chain approximation [25,[38][39][40], the current can be expressed in terms of the average density, ρ, and the NN correlation, C(ρ; ϵ) ≡ ⟨(11)⟩: In the limit of very strong repulsion (ε → ∞), J(ρ) = ρ(1 − 2ρ)/(1 − ρ) and so, the current is zero at the density value, ρ = 1/2. For ρ > 1/2, the negative value of J wrongly suggests a change in the direction of particle transport whereas the particles are still moving in the same direction. This can be resolved by invoking the particle-hole symmetry, which implies an equivalent picture of holes (i.e. empty sites) hopping in the opposite direction. So, replacing ρ by (1 − ρ) for ρ > 1/2, we have the current given by J(ρ > 1/2) = (1 − ρ)(2ρ − 1)/ρ which maintains the direction of particle transport and the J − ρ profile has a double-hump shape with a minimum (zero) at ρ = 1/2.
The maximum values of J are given by The critical value of ϵ, say ϵ c , at which the current changes its shape from unimodal to bimodal can be obtained by differentiating J(ρ, v), given by equation (15), twice with respect to ρ, setting its value equal to zero at ρ = 1/2 and then solving for ϵ-i.e. solving the equation J ′ ′ (ρ = 1/2, ϵ) = 0 for ϵ, we get ϵ c = ln(2 + √ 5) ≈ 1.445. In the case of an open system in the MC-phase (α > α * and β > β * ), the current is independent of the boundary rates, and changes its value from 1/4 to [1 + √ 2] −2 ∼ 0.1716 as ϵ increases from 0 to ∞. In order to obtain the v-dependence of the maximal current, we proceed as follows.
For ϵ > ϵ c , the above definition of R leads to 4x(1 − e −2ϵ ) = R(2 − R) which when substituted in equation (21) leads to Writing it this way helps us find the other two maxima of the (J, ρ) diagram besides the everpresent extremum (minimum in this case) at ρ = 1/2: so that ρ = 1/2 is always an extremum. Since which on simplifying leads to the following quadratic equation for R. .
Using the definition of R, we have Since the LHS is negative, the RHS has to be 1 − √ 2(1 − e −ϵ ). It will be equal to zero if ϵ = ln 2, the value at which J(ϵ) peaks. At this value, the LHS leads to the quadratic equation, 3ρ 2 − 3ρ + 1 = 0, that does not have a real root. Hence, the root is valid for ϵ > ln 2. Squarring both sides of this root and using the definition of R, we get Solving this quadratic equation for ρ, namely ρ 2 − ρ − f(ϵ) = 0, we have the roots, where ρ − and ρ + correspond to the first and the second maxima of the bimodal (J − ρ) diagram. This can be inferred by noting that In this limit, ρ − corresponds to the stationary density of non-interacting dimers in the MC-phase, see [26].
We focus attention on ρ − as this is what leads to the (J, ρ) diagram for an open system. In fact, particle-hole symmetry implies that ρ ± (1 − ρ ± ) = x = −f(ϵ). Using the definition of f(ϵ), we have ] which, in turn, leads to Substituting this expression for R in equation (22), we get At ϵ c = ln(2 + √ 5), J 2 (ϵ c ) = 1/4. So, J 1 (ϵ) and J 2 (ϵ) have the same value 1/4 at ϵ = ϵ c thus ensuring continuity for the J(ϵ) profile. Putting them together, we have the current profile, J m (ϵ) in the MC-phase, It can shown that there is an expression for the stationary current in the maximal phase P m , which relates the current J and the density of dimerons D(ε) for given values of ε and which is given by equation (15) and can be rewritten in the form is just the density of a NN pair, (01) [or equivalently (10)], which in the limit ε → ∞ is the density of dimerons (i.e. virtual hard dimers).

Appendix B. Density profiles and shocks
In this appendix we present some results of the dynamics of the densities of particles near the phase transition ε c1 as discussed above in section 5. It had been shown in many publications that density profiles can exhibit 'shock' profiles (or synonymously: phase boundaries, domain walls) separating regions of different densities [10,[28][29][30][31][32]. Such a phase boundary which separates high density and low density phases [30,32,41] can be very sharp, hence called a 'shock', i.e. a discontinuity in the density profile which can perform random walks along the x-axis. The phenomenon of phase separation in equilibrium systems has its non-equilibrium analog in models of the driven lattice gas [42][43][44][45].
In figure 15 we present data for the average density profiles ρ(x) around the high-to-low density transition near the coexistence point ε c1 between the high and the low density phases as discussed in section 5 on the coexistence in the P αβ regime ( figure 8(a)). The density profiles in figure 15 represent the situation in the steady state reached at very long times and is an average over both stochastic dynamics and initial distributions. It is observed that the density profiles ρ(x) for different ε change from a concave shape below ε c1 ≈ 1.30 to a convex shape above the critical point. Most interestingly very near the critical point ε c1 the shape exhibits a linear behavior. The profiles ρ(x) can be fitted by (28) which are shown in the figure by the dotted curves. The correlation length ξ(ε) describes the correlation between entering particles at x = 0 and exiting particles at x = L and can be estimated from the function equation (28). The length ξ(ε) has been estimated for the critical point ε c1 for different values of α and β. For all considered cases the correlation length diverges according to Figure 15. Transition of the average density profiles ρ(x) from concave to convex shape for various interaction parameters ε above and below the critical point ε c1 ≈ 1.30. The innermost curve in black corresponds to ε = 1.297 which is very close to the critical point. The left column correspond to the upper curves ε < ε c1 , the right column to ε > ε c1 .
which yields an estimate ν ≈ 1. This seems to be consistent with earlier theoretical predictions [17] for ε = 0. It should be noted that figure 15 provides a generic picture which seems to be valid for other cases of α and β in the P αβ regime close to ε c1 . It is generally believed and it has been shown that shocks can perform random walks along the x-axis [42,43,45]. At the microscopic level the structure, location and diffusion of a shock can be examined by simulations quite accurately and is discussed in the following. The shock position in a finite system is confined to a region compatible with the average density ρ(ε). Therefore, the fluctuations of the shock itself are limited to a finite region around the average location of the shock x s (t). Of course, at a very short time the precise location is blurred by local fluctuations of the density around x s (t). Therefore one has to average the shock over a short time interval ∆t, typically 700 < ∆t < 2000, in order to define the location x s (t) at time t. For the present simulations we used typically about 10 8 moves of the particles in the system. This single evolution of the system we call a 'trajectory'. Each average shock profile ρ s (x, t) comprises approximately 400 trajectories. The time scale is defined according to the continuous-time Monte Carlo algorithm [18] (section II). In figure 16 the shapes (a) and the location (b) of shock profile are presented. Figure 16(a) shows different shocks after a certain time t. A typical trajectory x s (t) of one the profiles shown in (a) is depicted in (b) which exemplifies the randomness of the trajectory. All shocks represent a jump between two densities, ρ + and ρ − . This is marked in figure 16 for ρ(ε) by isolated data points between 0.7 > ρ(ε) > 0.3 for the case α = 0.5 and β = 0.1 in the P αβ regime near the critical point ε c1 . One can define a shock accordingly by its lower and its upper 'bound' ρ − and ρ + , respectively, as Since according to figure 16(b) the shock moves randomly between the left and the right, i.e. the entrance and exit gates, respectively, this implies that the density fluctuates between ρ − and ρ + . The positions of a shock wave perform a diffusive motion around an average drift. It should be noted that drift in the present case of α = 0.5, β = 0.1 is very small and hence negligible. This is because a shock travels with a velocity (drift) [28] v s ∝ (1 − ρ + − ρ − ) as required by  mass conservation. Investigations of the localization and the stability of microscopic shocks had performed previously using second class particles and employing computer simulations [46,47]. The mean-squared displacement (MSD) of the position of the shock, x 2 s (t) ≡ [x s (0) − x s (t)] 2 , is shown in figure 17. At a characteristic time τ there is a crossover from diffusion x 2 s (t) ∼ D s t to a constant value x 2 s (∞) which depends on the distance ∆ε ≡ |ε − ε c1 | to the critical point and depends on the system size L. D s is the diffusion coefficient of the shock.
With increasing distance to the critical point the plateaus x 2 s (∞) decrease. This indicates a restriction of the random motions of the shocks to a range closer to either the entrance or the exit gates which is in variance with the average density profiles shown in figure 15. Since the fluctuations of the density are achieved by the random motions of the shocks, and since the diffusion of the shocks increases linerly in time, it is conceivable that near the critical point with L → ∞ the fluctuations of the density increase accordingly indefinitely. This one may consider as an analog to classical second order phase transitions. From this point of view it is seems to be in agreement with the divergent correlation length equation (29) of the density profile near the critical point ε c1 . The diffusion coefficient of the shock has been conjectured by Spohn [48] and proven by several authors [49,50] where q and p are the jump rates to the right and the left, respectively, and ρ + and ρ − are the average density at the right and the left side of the shock, respectively. The theoretical prediction of the diffusion coefficient is in reasonable agreement with the numerical estimate from our data D s ≈ 0.42 where we assumed q − p ≈ α − β and ρ + = 0.7 and ρ − = 0.3.